Polar Neptunes are Stable to Tides
Emma Louden Sarah Millholland
 



0000-0003-3179-5320]Emma M. Louden
Department of Astronomy, Yale University, New Haven, CT 06511, USA
emma.louden@yale.edu

0000-0003-3130-2282]Sarah C. Millholland
Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
MIT Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

§ ABSTRACT 

There is an intriguing and growing population of Neptune-sized planets with stellar obliquities near ∼90^∘. One previously proposed formation pathway is a disk-driven resonance, which can take place at the end stages of planet formation in a system containing an inner Neptune, outer cold Jupiter, and protoplanetary disk. This mechanism occurs within the first ∼10 Myr, but most of the polar Neptunes we see today are ∼Gyrs old. Up until now, there has not been an extensive analysis of whether the polar orbits are stable over ∼Gyr timescales. Tidal realignment mechanisms are known to operate in other systems, and if they are active here, this would cause theoretical tension with a primordial misalignment story. In this paper, we explore the effects of tidal evolution on the disk-driven resonance theory. We use both N-body and secular simulations to study tidal effects on both the initial resonant encounter and long-term evolution. We find that the polar orbits are remarkably stable on ∼Gyr timescales. Inclination damping does not occur for the polar cases, although we do identify sub-polar cases where it is important. We consider two case study polar Neptunes, WASP-107 b and HAT-P-11 b, and study them in the context of this theory, finding consistency with present-day properties if their tidal quality factors are Q ≳ 10^4 and Q ≳ 10^5, respectively.






§ INTRODUCTION


One of the many unexpected discoveries within exoplanetary science is that a planet's orbital angular momentum vector need not be aligned with its host star's spin vector. The angle between these two vectors is called the “stellar obliquity” or stellar spin-orbit misalignment. The solar obliquity is 7^∘ when measured using the mean plane of the planets <cit.>. This relatively small value is consistent with expectations from planet formation. However, for exoplanets, a broad range of stellar obliquities have been observed, and they exhibit various trends with stellar and planetary properties. (See  for a review.) For instance, hot Jupiters orbiting hot stars have high obliquities <cit.>, as do super-Earths and sub-Neptunes orbiting hot stars <cit.>.

The so-called “polar planets” are a class of systems with ∼90^∘ stellar obliquities. <cit.> computed the 3-dimensional stellar obliquity (Ψ) for 57 systems by combining constraints on stellar inclination angles with Rossiter-McLaughlin measurements of sky-projected stellar obliquities (λ). The findings revealed a significant number of systems with nearly perpendicular orbits (Ψ≈ 80^∘ - 120^∘), indicating a statistical preference for polar orbits over a full range of obliquities when using frequentist tests. Recent works using Bayesian methods have conflicting evidence for an overabundance of perpendicular orbits. <cit.> do not find strong evidence for the polar planets whereas <cit.> finds evidence that such a peak does exist when restricting attention to the sample of sub-Saturn planets and hot Jupiters orbiting F stars.  
Firm conclusions on the statistical nature of the polar planet population await a larger sample size. Regardless, the known polar planet systems are still dynamically puzzling. 

About half of the polar planets are  warm Neptune or super-Neptune-sized
planets. These polar Neptunes share a variety of distinguishing characteristics: moderately eccentric (∼0.05 < e ≲ 0.3) and polar orbits, puffy atmospheres that show evidence for mass loss in several cases, and evidence for exterior massive companion planets in
some cases (see Figure <ref>). Examples include WASP-107 b <cit.>, HAT-P-11 b <cit.>, HD 3167 c <cit.>, and K2-290 b <cit.>. Moreover, most polar Neptunes reside in or near the hot Neptune desert <cit.>, a sparsely-populated region of period-radius space that is thought to be sculpted by photoevaporation <cit.> and potentially high-eccentricity migration <cit.>. 



The mechanism(s) responsible for the unusual orbits of the polar Neptunes should also explain these other features. However, there is not yet consensus as to the origins of the polar Neptunes. Several theoretical scenarios have been proposed to explain 3D obliquities near Ψ≈ 90^∘. 


<cit.> showed that tidal dissipation within the host star can strand a system at a 90^∘ obliquity. This mechanism works through inertial wave dissipation in stars with convective zones, but it probably doesn't explain all perpendicular systems. 
In addition, primordial magnetic interactions between the protostar and the protoplanetary disk can torque the disk to high inclinations, with the planets then forming in these misaligned states <cit.>.

Similarly, <cit.> recently proposed that the central star in a polar planet system may have started as tight binary system that subsequently merged, and the polar planet formed in a protoplanetary disk that was perpendicular to the original binary orbital plane. As for another possibility, von Zeipel-Kozai-Lidov cycles and planet-planet scattering have long been invoked to explain spin-orbit misalignments <cit.>, but they have gained extra attention recently in the context of the polar planets <cit.>.



Separate from these ideas, another compelling theory for generating ∼90^∘ stellar obliquities is disk-driven resonance, which was proposed by <cit.> (hereafter ). This mechanism is envisioned to operate within the first ∼ 10 Myr in a system containing an inner Neptune-sized planet, outer giant planet, and a transition disk. As the disk dissipates, a secular inclination resonance is encountered between the inner Neptune and distant giant planet, which tilts the Neptune's orbit up to 90^∘. This theory accounts for the high obliquities, the non-zero eccentricities, and the Jovian companions with large mutual inclinations sometimes found in the polar Neptune systems. Its ability to explain all of these features at once makes disk-driven resonance particularly compelling. This theory is our primary focus in this paper. 



A key question is whether disk-driven resonance creates polar orbits that are long-term stable. Disk-driven resonance occurs at the very early epochs of formation, but the polar Neptunes we see today are mostly billions of years old. Within those ∼Gyrs of time, the orbits can be modulated by effects like tidal and stellar evolution, and these effects can potentially be significant. For example, tidally-driven orbital realignment is thought to operate for hot Jupiters orbiting cool stars, explaining the differences in the hot star and cool star populations <cit.>. There can also be more subtle effects, like inclination damping driven by interactions between the planet's spin and orbit <cit.>. There has not yet been a thorough investigation of tidal effects in the context of the disk-driven resonance mechanism. It is crucial to determine whether there is any impact on the mechanism's viability and the polar planet's long-term stability.





In this paper, we explore the degree to which the disk-driven resonance theory for polar Neptunes is influenced by tidal evolution. We investigate both short-term and long-term effects on the orbital and spin evolution. The organization of the paper is as follows. Section <ref> describes the disk-driven resonance mechanism and the properties of a system that affect its outcome. Section <ref> focuses on the system's early evolution while the disk is still present. We explore how tides impact the initial excitation into a polar orbit. Next, in Section <ref>, we study the long-term evolution of the system. We specifically explore whether tidal effects may erase the initial polar orbit. Finally, we discuss the results in Section <ref> and explore two case studies of observed near-polar Neptunes, HAT-P-11 b and WASP-107 b.



§ FORMATION OF SLIGHTLY ECCENTRIC POLAR NEPTUNES THROUGH DISK-DRIVEN RESONANCE



 introduced the disk-driven resonance theory to explain polar Neptunes. Here a close-in (≲ 0.1 AU) Neptune-sized planet interacts with both an exterior (∼1-5 AU) giant planet and an outer, slowly-depleting protoplanetary disk of gas and dust with the inner ∼ 5 AU evacuated. This is motivated by the existence of observed transition disks with inner regions depleted of gas and dust <cit.>. Initially, perturbations to the outer planet's orbit are dominated by the disk, while the inner planet's orbit is perturbed primarily by the outer planet and the star's quadrupolar gravitational field (oblateness). The inner planet is assumed to have zero eccentricity at the start (see Figure <ref>).  

The dynamical transition to a polar orbit begins with the establishment of secular inclination resonance. The disk induces nodal precession of the outer planet's orbit, Ω̇_out, at a rate proportional to the disk mass, which we assume decreases smoothly and slowly. The outer planet, in turn, induces nodal precession of the inner planet's orbit, Ω̇_in.  Initially, |Ω̇_out| > |Ω̇_in|. As the transition disk dissipates, |Ω̇_out| decreases until the precession rates of the two planets match and
a scanning secular resonance <cit.> is crossed. The inclination of the
inner planet is forced to increase because it stays in resonance by maintaining the same nodal precession frequency as the outer planet, and the precession frequency slows as the inclination rises. 









The inner planet's inclination growth stops at a critical inclination I_crit, which, under certain conditions, is polar. This is a generalization of the critical inclination in the von Zeipel-Kozai-Lidov mechanism <cit.>. That is, as soon as the critical inclination is reached, the inner planet undergoes eccentricity excitation that detunes the secular resonance and halts the inclination growth. If I_crit is near ∼90^∘ or undefined, then the planet can reach a polar orbit before the resonance is detuned. 

The critical inclination is dictated by the relative perturbations on the inner planet from the outer planet, general relativity (GR), and stellar oblateness. In the absence of GR and stellar oblateness, the critical inclination would be equal to 39.2^∘. However, the short-range forces suppress eccentricity excitation due to the periapse precession they induce <cit.>, and thus a larger inclination is required to induce eccentricity excitation. The critical inclination can be calculated using the dimensionless quantities η_GR and η_⋆. Here η_GR measures the strength of GR perturbations relative to the planet-planet interactions, and η_⋆ similarly measures the strength of the stellar quadrupole perturbations relative to the planet-planet interactions. They are defined as 

    η_⋆ = 2 J_2 M_⋆/m_2R_⋆^2 a_2^3/a_2^5(1-e_2^2)^3/2


    η_GR = 8 G M_⋆/c^2a_2^3/a_1^4M_⋆/m_2(1-e_2^2)^3/2,

where M_⋆ and R_⋆ are the stellar mass and radius, c is the speed of light, a_1 is the semi-major axis of the inner Neptune, a_2 and e_2 are the semi-major axis and eccentricity of the outer cold Jupiter, and m_2 is its mass. In the equation for η_⋆, J_2 is the coefficient of the quadrupole moment of the star's gravitational field, and it depends on the stellar rotation period P_⋆, the Love number k_2⋆ and R_⋆ as <cit.>

    J_2 ≈k_2⋆/34π^2/P_⋆^2R_⋆^3/G M_⋆.





In terms of  η_⋆ and η_GR, the critical inclination is given by

    I_crit = sin^-1(4 + 4  + / 10 + 5 )^1/2.

Figure <ref> shows how the combined values of η_⋆ and η_GR determine the critical inclination. For ≥ 6 +, the argument of the arcsine is greater than 1 and thus all inclinations are stable, meaning the planet would reach a polar orbit. We also note that the behavior in the parameter space in Figure <ref> changes at I_crit = 63.4^∘, since apsidal precession from stellar oblateness is prograde for inclinations < 63.4^∘ and retrograde for > 63.4^∘. 


To summarize, the inner planet can reach a polar and possibly eccentric orbit, but the outcome is dependent upon a “tug-of-war” between the effects of the stellar oblateness, the GR effects, and the cold Jupiter. This competition is summarized succinctly in the quantities η_⋆ and η_GR, which are most sensitively dependent on six key parameters of the system: mass and radius of the star (M_⋆, R_⋆), rotation period of the star (P_⋆), mass and semi-major axis of the cold Jupiter (m_2, a_2), and semi-major axis of the inner Neptune (a_1).



§ EARLY EVOLUTION WITH TIDES



The disk-driven resonance theory could be affected by tidal evolution at multiple epochs in time. First, there is the early ≲10 Myr period during the disk phase when the resonance is encountered and the inclination (and sometimes eccentricity) excitation occurs. Second, after the polar orbit is established, there is the long ∼ Gyrs period during which orbital realignment might occur. We will investigate these phases separately, starting with the early evolution. Here we explore how the initial resonant encounter is impacted by tides. 

The initial resonance excitation involves a complicated set of interactions between the two planets, disk, and oblate star. We opt to perform direct N-body simulations to fully capture these interactions and the corresponding tidal and spin dynamics. Since such simulations are expensive, we only investigate a representative case study rather than assess the full parameter space. We use a numerical integrator that evolves the system using instantaneous accelerations in the framework of <cit.>. This code self-consistently accounts for the accelerations from (1) standard Newtonian gravity, (2) the oblateness of the host star, (3) constant-time lag tides, and (4) the outer protoplanetary disk. It evolves the planetary orbits and spins simultaneously. The code was originally developed and described in <cit.>, and further details can be found therein. The disk effects were added in <cit.>, but the disk was taken to be infinite. Here we model the disk as a finite transition disk between R_in and R_out with surface density Σ = Σ_0 (r/R_out)^-α. The gravitational potential is given by <cit.> as

    Φ   = --α + 2/1 - η^-α + 2G M_disk/R_out[ 1 - η^1 - α/1 - α. 
       . + -1 + η^-1 - α/1+αr_p^2/2 R_out^2(-1 + 3/2sin^2θ_p) ],

where η = R_in/R_out, M_disk is the mass of the disk, r_p is the radial position of the planet, and θ_p is the polar angle of the planet. The disk mass decreases over time according to M_disk = M_disk,0/[1 + t/t_disk)]^3/2.

We run two simulations of a representative synthetic system with all being equal except one simulation includes tides and the other does not. The parameters of the system are the same as the simulation in the rightmost panel of Figure 2 of , so as to aid comparison to their results. This case corresponds to a resonant encounter that yields a nearly polar and eccentric planet. The system consists of a Solar-mass star, a Neptune-mass and Neptune-radius planet at 0.07 AU with an initial eccentricity of 0.01 and inclination of 1^∘, a 4  M_Jup-mass planet at 2 AU in a circular orbit with an inclination of 5^∘, and a transition disk with edges R_in = 3 AU, R_out = 30 AU and an initial mass of 50  M_Jup. The power law exponent of the disk surface density is set to α = 1, and the disk decay timescale is set to t_disk=1 Myr. The star has a radius of 1.3  R_⊙, rotation period of 7 days, and Love number of k_2⋆ = 0.2. For the simulation of the planet that includes tides, we set the Love number of the planet to k_2p = 0.4 and tidal quality factor equal to Q = 10^5.



Figure <ref> shows the results of the simulations. We show the inclination and eccentricity evolution of the inner Neptune-mass planet. 

Focusing first on the case without tides, we see that this set of parameters corresponds to a case in which the inner Neptune encounters the inclination resonance and tilts to I_crit∼ 70^∘, at which point the eccentricity grows rapidly and the resonance is detuned. The inclination growth is then stalled just above 70^∘. The behavior is the very similar to that seen in , which simultaneously validates our code and stands as an encouraging sign since we are using a different numerical method (N-body, whereas  used secular simulations).




As for the simulation with tides turned on, the results are pretty similar overall. The inner planet's inclination evolution is nearly identical to the case without tides; it increases to the same near-polar critical inclination with a small lag. The eccentricity also increases when the critical inclination is reached, but the eccentricity growth is slightly delayed and muted by the tidal damping. We expect that further damping would be observed if the simulation was run for longer. For completeness, we note that we explored additional simulations with different system parameters (particularly the tidal quality factor of the inner planet) and found more complicated eccentricity behavior in some simulations. In these cases, the eccentricity decreased rapidly after its initial excitation and then stabilized at a lower value around e∼0.1.

Exploring these dynamics across a broader parameter space is beyond the scope of this work. However, we can at least take away from this case study that tidal evolution does not significantly affect the inclination evolution but it can impact the eccentricity excitation following the resonant encounter. A more exhaustive parameter exploration could be considered in future work.

To round out our analysis of the early evolution, we also explore the behavior of the inner Neptune's spin vector during the disk-driven resonance encounter. Tidal dissipation in a planet is driven not only by its eccentricity but also its “planetary obliquity” (the angle between the planet's spin vector and the orbit normal vector) <cit.>. Obliquity tides are theorized to drive orbital evolution and atmospheric inflation in a variety of exoplanetary systems <cit.>. It is thus relevant to investigate whether we expect an obliquity enhancement in the inner Neptune, which could affect the subsequent orbital evolution. Planetary obliquities can be excited by secular spin-orbit resonances, which occur when the planet's spin-axis precession frequency becomes commensurable with one of its orbit nodal precession frequencies. In Appendix <ref>, we investigate this scenario in detail. We find that the parameter space in which the inner Neptune would obtain a high obliquity through a secular spin-orbit resonance does not overlap with the parameter space for it to obtain a polar orbit through disk-driven resonance. Therefore, we do not anticipate that the inner Neptune would experience significant obliquity tides, at least not due to obliquity excitation from a simultaneous resonance crossing.



§ LONG-TERM SECULAR EVOLUTION WITH TIDES
 




With a clear picture of the impact of tides on the short-term initial phase of resonance capture, we now turn to the subsequent long-term ∼Gyrs period. The observed polar Neptunes are billions of years old and have been subject to tidal forces throughout their lifetime. 
Here we explore a wide range of parameters and model the long-term evolution to test if the tides cause a meaningful change like destabilization of the high inclination. Such an outcome would cast doubt on the theory that the polar orbits were established early on.



 §.§ Simulation set-up


The N-body integrator we used in the last section is too slow to run over the long timescales of interest here. Instead, we perform secular (orbit-averaged) simulations using [<https://github.com/djmunoz/kozaipy>].  is a  package for integrating the secular equations of motion of few-body, hierarchical systems including tidal friction. The equations are the same as in <cit.> and originally <cit.>. The code evolves the orbits and spins simultaneously and accounts for mutual perturbations due to the Newtonian gravitational interactions, quadrupolar distortion of the star due to rotation and tides, tidal friction in the constant-time lag approximation, and general relativity. 

We use a Monte Carlo approach to create a suite of simulation outcomes for different system configurations. As discussed in Section <ref>, the orbital configurations resulting from the resonance interactions are dictated by η_⋆ and η_GR. We thus choose key parameters to vary based on η_⋆ and η_GR.  All relevant variables are defined in Table <ref>, but not all of these are varied in each simulation, as we will describe in more detail in the next section. 

The viscous timescale t_v is related to the tidal quality factor Q of the planet by[For the stellar tidal quality factor, k_2,1 would be replaced by k_2⋆, m_1 and R_1 would be replaced by M_⋆ and R_⋆, and t_v,1 would be replaced by t_v,⋆. Our standard t_v,⋆ values correspond to Q_⋆≈ 10^6-10^7. ]




    Q_1 = 4/3k_2,1/(1+2k_2,1)^2Gm_1/R_1^3(a_1^3/GM_⋆)^1/2 t_v,1.

 

The simulations are designed to model the system's evolution just after the disk-driven resonance encounter has taken place. Given a set of parameters, we calculate the critical inclination (equation <ref>) and set this as the inner Neptune's initial inclination. To account for the fact that the initial inclination resonance occurs during the star's pre-main sequence phase while the subsequent long-term evolution is during the main sequence phase, we determine the critical inclination using R_⋆ = 1.3  R_⊙ and P_⋆ =7 days to reflect fiducial pre-main sequence values. Then, for the integration, we set the stellar parameters using the sampled variables and keep them constant. We do not evolve the stellar properties over time because we expect that tides and planet-planet interactions dominate the orbital evolution. The integrations run for 6 Gyr.




 §.§ “Tug-of-war” simulation suites





We define three simulation suites that span the “tug-of-war” between the outer Jupiter, stellar oblateness, and general relativity. In Suite 1, we sample all parameters according to the schemes outlined in Table <ref>. This suite is designed to represent the full range of outcomes given any plausible parameter choices. Suites 2 and 3 are more targeted. In Suite 2, we vary R_⋆ and a_1, and in Suite 3, we vary M_⋆, a_1, and m_2. For all three suites, we run ∼5,000 different simulations. 



The results of Suites 1, 2, and 3 are shown in Figure <ref>. 



To understand the degree of variation in the inclination evolution, we compute the standard deviation of the inner planet's inclination for each run. Recall that the planet starts at the critical inclination, so a small standard deviation would indicate that the misaligned orbit is long-term stable. It is clear from Figure <ref> that very small standard deviations are the usual outcome. Specifically, for Suites 1, 2, and 3, there are respectively 375, 292, and 133 runs where the standard deviation of the inner planet's inclination is greater than 1^∘.

Suite 2 shows a distinct curve of initial mutual inclination vs. a_1. This is a consequence of initializing the inner planet at the critical inclination. As described earlier, Suite 2 varies R_⋆ and a_1, but R_⋆ is set to 1.3  R_⊙ for the critical inclination calculation. Thus, there is only one degree of freedom in setting I_crit, and that is reflected in the curve of the inclination versus a_1. For η_GR > 6 + η_⋆, I_crit = 90^∘ and fully polar orbits are created. The curve is also seen as an approximate lower envelope in the Suite 1 and 3 simulations where more parameters are varied.



Overall, the simulation results offer a few notable takeaways. First, and most importantly for this work, the polar planets are extremely stable to tides in all three suites. While there is some deviation in the inner planet's inclination in some of the runs in each suite, the planets in polar orbits are all stable over time. This is further illustrated in Figure <ref>. The top (bottom) panel shows the subset of runs from all three suites with a standard deviation greater (lower) than one. The top panel includes no runs with planets in polar orbits. There are some near-polar runs that exhibit changes in inclination, but the truly polar orbits are stable. This conclusion is further cemented by the bottom panel. Even within the low standard deviation subset, the runs with the highest standard deviation are predominately concentrated around 50^∘. 



The second key result pertains to the location of the instability hotspot for all three suites. In Figure <ref>, the color bar shows that the cases with the highest standard deviation of the inner planet's inclination have an initial inclination in the range ∼45^∘-80^∘. Looking more closely at these higher standard deviation runs in Figure <ref>, we see an evolution reminiscent of high-eccentricity migration. The eccentricity slowly increases over time, and when it gets very large, the orbit undergoes rapid tidal migration and circularization. The planet ends up in a stable circular orbit closer to the star with a lower inclination. Such outcomes are the only mechanism in which we observed large-scale inclination variations.

One caveat of these simulations is that we do not know the exact eccentricities emerging after the initial resonant encounters. As shown in Figure <ref>, the eccentricity can be excited during the resonant encounter but then influenced by tidal evolution to stabilize at smaller values. We set the eccentricities somewhat arbitrarily but informed by the earlier simulation and the observed population of Neptune-sized planets. To test the sensitivity of this assumption, we ran a suite of simulations where the eccentricity was allowed to vary more widely and found that the stability of the polar planets as well as the location of the hotspot was unchanged. There was the same order of magnitude of runs with standard deviation of the inner planet's inclination greater than 1^∘ and the hotspot appears in the same range of initial mutual inclination. 

It is useful to contrast our results with  <cit.>, who also explored the secular dynamics of similar hierarchical two-planet systems but in a slightly different parameter regime. They identified a planetary spin-facilitated mechanism that induced unexpected and significant damping of the mutual inclination in long-term orbital dynamics. To validate our own simulations and compare with their results, we run another suite of long-term simulations similar to those presented in <cit.>. The set-up and results are described in Appendix <ref>. In short, we reproduce their findings and show that a particular regime of parameter space (namely, where the inner planet is more massive than Neptune) results in inclination damping by tens of degrees. The effect is significant, but it doesn't overlap with the parameter regime of interest in this work.







§ CASE STUDIES

Thus far, we have used generalized synthetic planetary systems with randomized parameters to explore the long-term dynamics following disk-driven resonance. Here we turn to observed polar Neptune systems to gain further insights. We analyze two case study planets in polar and slightly eccentric orbits: HAT-P-11 b and WASP-107 b. Both of these planets have distant giant planet companions, making them strong candidates for a disk-driven resonance origin scenario . Moreover, they are actively undergoing atmospheric mass loss <cit.>, which could be a side effect of enhanced internal heating driven by tides from their eccentric orbits <cit.>. We aim to examine the long-term evolution of the inner planets following the initial establishment of their polar orbits via disk-driven resonance. We will use dynamical simulations to track the stability of the orbits under the influence of tides. Moreover, we will explore whether the long-term orbital evolution allows us to place bounds on the planets' tidal quality factors. 

Before proceeding, we note that HAT-P-11 b and WASP-107 b were explored in two recent papers by <cit.> and <cit.>, who showed that these planets could be explained by a high-eccentricity migration pathway induced by von Zeipel-Kozai-Lidov (ZKL) cycles. Disk-driven resonance is a potential alternative to this. We do not necessarily favor a certain formation pathway, as both theories have advantages. One advantage of disk-driven resonance is that the inner and outer planets form coplanar. Conversely, for ZKL cycles to start, the planets must first obtain a large mutual inclination, which could be produced by planet-planet scattering or a stellar flyby. However, we also note that both processes could work in concert, with disk-driven resonance producing an initial tilt and then scattering or secular interactions leading to further eccentricity and inclination changes. The outer planet in HAT-P-11 has a substantial eccentricity and possibly also a misaligned stellar obliquity, which is indicative of a history of scattering <cit.>.






 §.§ HAT-P-11 b

The HAT-P-11 planetary system, comprising two known exoplanets, offers a compelling glimpse into the dynamical processes producing planets in slightly eccentric polar orbits. Closest to the host star is HAT-P-11 b, a Neptune-sized planet discovered utilizing the transit method <cit.>. With a mass approximately 26 times that of Earth and a radius four times larger, HAT-P-11 b orbits its host star every 4.88 days at a close distance of about 0.053 AU with a moderate eccentricity of 0.218^+0.034_-0031 <cit.>. The projected obliquity is λ = 103^+26_-10 degrees, and the 3D obliquity is ψ = 97^+8_-4 degrees <cit.>. <cit.> discovered the outer Jupiter-sized planet and suggested a large mutual inclination between the inner sub-Neptune and outer Jupiter, which was recently confirmed by <cit.>. 



Here, we explore the idea that HAT-P-11-b's polar orbit was established via a disk-driven resonance. We consider the post-disk evolution of the system and perform simulations that probe the long-term stability of HAT-P-11-b's orbit. The initial eccentricity following resonant excitation is unknown, as is the tidal Q of the planet. The planet may also have undergone tidal migration, ending up in a smaller orbit than it started at. Thus, we construct a suite of 5,000 simulations with different initial eccentricities, semi-major axes, and tidal Q values. Specifically, we allow the eccentricity to vary from  0.02 to 0.95, the semi-major axis from 0.05 AU to 0.06 AU, and Q from 10^3 to 10^7, as summarized in Table <ref>. For all simulations, we start the inner planet at the critical inclination calculated using using the fiducial pre-main sequence stellar parameters (R_⋆ = 1.3  R_⊙,  P_⋆ = 7  days). The system is evolved using  for 6 Gyr.


Of the 5,000 simulations, 2,089 were excluded because the planet was engulfed by the star before reaching the end of the simulation. 



The orbital evolutions of all simulations are shown in Figure <ref>. The inclination stays remarkably stable over a 6 Gyr timescale (0.003% average change in mutual inclination). The planet's eccentricity and semi-major axis evolves due to the tidal dissipation. In Figure <ref>, the blue lines show the 120 runs where the final semi-major axis and eccentricity match the observed values with tolerances of 0.005 AU and 0.05, respectively. These cases show some decrease in the semi-major axis of the inner planet. Only a certain range of Q values results in an orbit that agrees with the present-day configuration. We will discuss this further in Section <ref>.







 §.§ WASP-107 b




The WASP-107 system hosts the well-studied super-Neptune WASP-107 b <cit.>. 

This planet has become a focal point for understanding the atmospheric properties and orbital dynamics of puffy planets in close proximity to their host stars. With a mass only about 30 times that of Earth yet a nearly Jupiter-sized radius, WASP-107 b completes an orbit around its host star every 5.7 days with a semi-major axis of 0.055 AU and eccentricity of 0.06±0.04 <cit.>. The planet's orbit is polar, with a sky-projected obliquity of λ = 118^+38_-19 degrees <cit.> and a 3D obliquity of ψ = 92.6^+30.7_-1.8 degrees <cit.>. <cit.> detected a 0.36  M_Jup companion planet on a wide eccentric orbit. Alongside dynamical studies, WASP-107 b has provided insights into the nature of inflated super-Neptunes through studies of its escaping atmosphere <cit.> and high internal heat flux <cit.>, both suggestive of strong tidal heating <cit.>. The planet's unusual orbit and interior/atmosphere make it a prime subject for detailed characterization. 






Just like our case study of HAT-P-11, we run 5,000 simulations of the WASP-107 system that explore the long-term evolution of planet b after the hypothesized disk-driven resonance encounter established its polar orbit. We allow the eccentricity to vary from  0.02 to 0.95, the semi-major axis from 0.05 AU to 0.06 AU, and Q from 10^3 to 10^7. 
Of the 5,000 simulations, 379 were excluded because the planet fell into the star before reaching the end of the 6 Gyr.  

The system evolution is shown in Figure <ref>. Like HAT-P-11 b, the inclination stays extremely stable over a 6 Gyr timescale (0.005% average change in mutual inclination). The planet's eccentricity and semi-major axis evolve due to tidal dissipation. In Figure <ref>, the blue lines show the 28 runs where the final semi-major axis and eccentricity match the observed values with tolerances of 0.005 AU and 0.05, respectively.  



 §.§ Constraining Tidal Q







Our simulations used a range of values for the planet's tidal quality factor Q, but only some of the resulting systems matched the observed system properties. Here we explore this in more detail in an effort to obtain constraints on the Q values. The tidal quality factor is a metric for the efficiency of tidal energy dissipation, defined by 

    Q^-1 = -1/2π E_0∮dE/dt dt,

where E_0 encapsulates the peak energy stored within the tidal distortion of the planet, and the integral represents the energy dissipated in a complete cycle <cit.>. Smaller values of Q correspond to more efficient dissipation and vice versa. Within the Solar System, the terrestrial planets and satellites have Q ∼ 10^2 <cit.>, while Uranus and Neptune have Q ∼ 10^4 <cit.>. 


For Neptune-like exoplanets, we expect ∼10^4-10^5, although the detailed values would depend on the body's specific rheological properties.

For both case study systems, we do not know the inner planet's initial eccentricity nor its Q, but we do have measurements for the current eccentricity. Both planets have non-circular orbits, indicating that Q cannot be too low, otherwise the orbits would have circularized. We isolate the simulations that produce final eccentricities matching the observed values (blue lines in Figures <ref> and <ref>). By identifying their initial eccentricities and Q values, we can obtain constraints on the Q values for HAT-P-11 b and WASP-107 b within the framework of this theory. These simulations place the first Q constraints for either of these planets. As shown in Figure <ref>, we find that Q≳ 10^5 for HAT-P-11 b, and Q≳10^4 for WASP-107 b. The constraints depend on the initial eccentricity as well. Overall, these results suggest a level of dissipation that is relatively inefficient but consistent with expectations from the Solar System planets. Moreover, we expect that the constraints on Q would further increase if we included the evolution of the planet's radius in our simulations <cit.>, since more significant past radius inflation would need to be balanced by a larger Q to result in the same timescales of tidal evolution.



§ CONCLUSIONS











The polar Neptunes are a mysterious new class of exoplanets. To investigate their origins, we extended the study of the disk-driven resonance theory proposed by <cit.>. 

This scenario suggests that the orbit tilting happens early, within the first ∼10 Myr of the system's lifetime. Before this work, there had not yet been investigations of the planet's subsequent ∼Gyrs of orbital evolution, during which tides might feasibly realign the orbits. We performed a comprehensive analysis of the impact of tidal evolution on both the initial excitation of the polar orbits and their long-term stability. Our main takeaways are as follows:


  * Short-term N-body integrations of the initial disk-driven resonance encounter reveal that tidal effects do not change the inclination evolution resulting from the resonance, but they can impact the eccentricity, leading it to stabilize at a smaller value than seen in a tides-free simulation. A broader parameter space exploration would be valuable in future work. 

  * Long-term secular simulations show that polar Neptune orbits are extremely stable. The orbits do not realign as a result of tidal dissipation in the star or in the planet. This result implies that there is no constraint on when the polar misalignment must be initially established, and thus a primordial excitation theory such as disk-driven resonance is in no tension with the observations. 

  * On the other hand, systems emerging from the initial disk-driven resonance with a more moderate, non-polar obliquity in the range of ∼45^∘-80^∘ can often undergo substantial realignment due to tidal interactions. In some cases, these interactions are mediated by effects of spin-orbit coupling. 

  * HAT-P-11 and WASP-107 are exemplary case study systems, as they contain both polar Neptunes and outer giant planets. Long-term secular simulations show stable orbits for the polar Neptunes. Moreover, working within the framework of disk-driven resonance, we constrain the tidal quality factors of the planets to Q ≳ 10^5 for HAT-P-11 b and Q ≳ 10^4 for WASP-107 b, similar to Uranus and Neptune. Disk-driven resonance is not the only possible formation pathway for these planets, but it is consistent with available constraints.


Our theoretical exploration of disk-driven resonance can be complemented with further observational tests. Disk-driven resonance predicts that the polar Neptunes are accompanied by distant giant planets with small stellar obliquities and ∼90^∘ mutual inclinations with respect to the inner Neptunes . This can be tested with future observations that constrain the full 3D architectures of systems containing polar Neptunes, including the stellar obliquities of the outer giant planets and the mutual orbital inclinations. Such constraints are already possible through the usage of Hipparcos and Gaia astrometry <cit.>, but they will become significantly more prevalent with Gaia DR4 <cit.>. Moreover, in systems where outer giant planets have not yet been discovered, many more long-period perturbers may soon be identified using long-term radial velocity monitoring and Gaia astrometry <cit.>. There is thus a promising future ahead for a more complete understanding of polar Neptune systems and their enigmatic origins.
















§ ACKNOWLEDGMENTS

This material is based upon work supported by the National Science Foundation under Grant No. 2306391. We thank the NSF for their support. We thank Cristobal Petrovich, Guðmundur Stefánsson, Alexandre Correia, and James Owen for helpful conversations. We are grateful to the Yale Center for Research Computing for the use of the research computing infrastructure, and we especially thank Tom Langford for guidance and support using the Grace cluster. Additionally, we gratefully acknowledge access to computational resources through the MIT Engaging cluster at the Massachusetts Green High Performance Computing Center (MGHPCC) facility.







§ SECULAR SPIN-ORBIT RESONANCE



The early orbital evolution responsible for the secular inclination resonance encounter could potentially result in other resonances being crossed. In particular, we are interested in the secular spin-orbit resonance, which involves a commensurability between the orbit nodal precession frequency, g = Ω̇, and the planetary spin-axis precession frequency, α. This resonance excites the planetary obliquity, and a non-zero planetary obliquity results in another source of tidal dissipation (in addition to eccentricity). Secular spin-orbit resonances put planets into so-called “Cassini states” <cit.>, equilibrium configurations of a body's spin pole in a uniformly precessing orbit frame. These have been well-studied in both Solar System and exoplanetary system contexts <cit.>. 

To excite the planetary obliquity, the ratio of frequencies, |g|/α, must cross unity from above. That is, a necessary condition for resonance crossing is that |g|/α > 1 initially. Cassini states are strictly defined only for uniform orbital precession when g = Ω̇ and the inclination I are constant. In cases where there are multiple sources of perturbations, including the scenario considered here, the precession is non-uniform. The planet's inclination and node evolution comprises a superposition of several modes with multiple frequencies {g_i}. However, for well-separated frequencies, the resulting spin vector equilibria are “quasi-Cassini states,” which behave approximately as Cassini states with g equal to one of the g_i modes. 

In the system considered here, a secular spin-orbit resonance for the inner planet would involve the fastest g_i frequency since α is generally fast due to the inner planet's close-in orbit. The full set of g_i frequencies can only be calculated through a complete secular analysis of the system, but they can be approximated by pairwise interactions. The fastest frequency for the inner planet is g_⋆,1, the nodal precession frequency of the inner planet induced by stellar oblateness <cit.>.

This frequency decreases over time due to stellar contraction and rotational spin down, meaning that |g_⋆,1|/α_1 evolves in the right direction for resonance crossing. We can thus assess the conditions in which the inner planet encounters the resonance by considering the initial value of |g_⋆,1|/α_1 

across the parameter space. 

Figure <ref> shows the initial value of |g_⋆,1|/α_1 in a parameter space of a_1 on the x-axis and a combination of stellar parameters on the y-axis.

Obliquity excitation is expected when |g_⋆,1|/α_1>1 initially. The figure can be compared with  Figure 4, where the outcomes of the disk-driven resonance are shown in the same parameter space. We find that the region of parameter space where polar orbits are created correspond to |g_⋆,1|/α_1 ≪ 1, far outside of the region where secular spin-orbit resonance is expected. This indicates that planets that obtain polar orbits through disk-driven resonance probably do not also obtain high obliquities, at least not through a simultaneous secular spin-orbit resonance encounter. 













§ COMPARISON WITH <CIT.>
 

<cit.>, hereafter ,  studied the dynamics of hierarchical two-planet systems and found that there can be a significant and surprising source of inclination damping. This does not actually result from a direct consequence of tidal effects on the orbit but rather from a complex interaction between tidal effects on the spin and secular eccentricity oscillations that lead to a gradual pumping of eccentricity and damping of inclination. Since we don't see widespread inclination damping in our simulations, it is useful to compare to . They considered cases where the inner planet was more massive and on a wider orbit than the close-in Neptunes we are studying here. This suggests that we do not see inclination damping simply because our systems are in a different parameter regime. We can confirm this interpretation and simultaneously validate our simulations by exploring cases with parameters similar to those used in . 

We first reproduce a simulation performed by  of the system HD 74156, which contains two giant planets on eccentric orbits <cit.>. We use the same set-up as  Figure 6, and the results are shown in Figure <ref>. The orbital evolution is nearly identical to that seen in . The inner planet experiences significant semi-major axis, inclination, eccentricity damping as a result of the coupled spin-orbit and tidal evolution. The close agreement of our simulation results with  offers a strong validation of our numerical method. 

We further explore a comparison to 's findings using a larger set of simulations. Similar to Section <ref>, we create a simulation suite of the long-term orbital evolution of synthetic systems where the inner planet is initialized at the critical inclination. However, we use parameters that are more similar to the systems explored in . The parameters that differ include the inner planet's semi-major axis and the planet masses. The outcomes are consistent with  in that they show significant inclination damping. About 89% of the runs resulted in a standard deviation of the inner planet's inclination greater than 1^∘, compared to only 6% of our earlier runs in the disk-driven resonance regime that had this high level of inclination variation (e.g. Figure <ref>). Clearly, these systems are in a different regime than the polar Neptunes we explored in Section <ref>. 





aasjournal



