Next Article in Journal
Exploring Neutrino Mass Orderings through Supernova Neutrino Detection
Previous Article in Journal
Lorentz Symmetry Violation Effects Caused by the Coupling between the Field fμγ5 and the Derivative of the Fermionic Field on One-Dimensional Potentials
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Secular Orbital Dynamics of the Possibly Habitable Planet K2-18 b with and without the Proposed Inner Companion

1
United States Naval Observatory, 3450 Massachusetts Ave. NW, Washington, DC 20392-5420, USA
2
Teza Technology, 150 N Michigan Ave, Chicago, IL 60601, USA
*
Author to whom correspondence should be addressed.
Universe 2023, 9(11), 463; https://doi.org/10.3390/universe9110463
Submission received: 6 October 2023 / Revised: 20 October 2023 / Accepted: 23 October 2023 / Published: 28 October 2023
(This article belongs to the Section Planetary Sciences)

Abstract

:
The transiting planet K2-18 b is one of the best candidates for a relatively nearby world harboring biological life. The long-term orbital evolution of this planet is investigated using theoretical and purely numerical techniques for two possible configurations: A single planet orbiting the host star, and a two-planet system including the proposed inner planet close to the 4:1 mean motion rationalization. The emphasis is made on the secular changes of eccentricity and orbital inclination, which are important for the climate stability of the planet. It is demonstrated that the secular orbital dynamics of planet K2-18 b with an internal companion are accurately represented by the periodic eccentricity and inclination exchange on the time scales of a few Kyr. A single planet is not expected to experience fast orbital changes, with the much weaker tidal and rotation-driven perturbations mostly reflecting in a slow periastron and nodal precession. The tidal decay of the orbit is too insignificant on the time scale of the stellar age. However, the conditions for the habitability of a single K2-18 b planet are much improved if, like the Earth, it rotates faster than the mean motion and its rotation angle is tilted by a hypothetical moon. Milanković’s cycles of the habitable planet’s climate are discussed for both configurations.

1. Introduction

K2-18 is an M2.8 main sequence dwarf at a distance of 38.1 pc from the Sun. It harbors at least one planet (K2-18 b) detected in the photometric light curves of the extended Kepler mission [1,2]. The orbital period of this planet is 32.94 d [3], and its estimated bulk density of 3300 kg m 3 suggests a terrestrial central body covered with a massive ocean [4]. Observing the planetary transits with the Hubble Space Telescope, ref. Benneke et al. [5], the presence of water vapor and clouds in the atmosphere of planet b were detected. Coupled with the estimated insolation flux, which is only slightly higher than the insolation of Earth [6], these properties make K2-18 b one of the best candidate habitable worlds with thriving biological life. Precision radial velocity observations confirmed the presence of planet b, but also indicated a possible inner companion with a period of 8.962 ± 0.008 d [4]. The mass and radius of this proposed planet c are only slightly smaller than the values for planet b. A downward revision of the projected mass M c sin i followed a more detailed analysis of the available radial velocity data from two separate instruments [7]. A novel technique to process combined spectroscopic data further strengthened the case for an inner planet and updated the estimated mass to 6.99 M and the period to 9.2 d [8]. Still, doubts linger in the reality of the inner planet, and it is listed as a “controversial” entry in the NASA Exoplanet Archive available online (https://exoplanetarchive.ipac.caltech.edu, accessed on 27 August 2023). The tentative inner planet is not transiting; thus, its orbital inclination should differ from that of planet b by at least 2 . The period ratio is approximately 3.7, which is close to the potential 4:1 mean motion resonance (MMR). Are such systems stable long-term?
Another consideration to be taken into account is the rate of rotation of the host star. The projected surface velocity from the broadening of the spectral lines is v sin i = 2.0 km s 1 [9], which yields an upper bound rotation period of 8.8 d for the estimated radius. This is close to the orbital period of the tentative planet c. Photospheric spots and plages can cause measurable modulation of the radial velocity curve, which, in most cases, manifests itself as an enhanced jitter because of the short life-times of these formations and their small effective area [10]. Magnetically active stars, on the other hand, are known to have persistent structures (such as spot groups) that can last for multiple rotations and occupy a sizable portion of the visible disk. They cause quasi-sinusoidal variations of the Doppler shift, which can be confused with an exoplanet signal. K2-18, however, does not show significant signs of magnetic activity or binarity in the available data. Persistent photometric structures are not expected for this star. Precision radial velocity measurements over extended timescales should resolve this open issue by confirming or rejecting the persistent phase of the radial velocity modulation.
The goal of this study is to investigate the consequences of both options, the single planet and two-planet system, for the secular orbital evolution of the possibly habitable planet b. Variations of the semimajor axis, eccentricity, and inclination are very important factors for the considerations of habitability and climate modeling. The plan of this paper is as follows. In Section 2, we investigate the basic secular exchange of eccentricity and inclination for the two-planet option. Our approach combines high-accuracy and long-term numerical simulations with the theoretical models derived from the classical Laplace–Lagrange analysis.

2. Eccentricity and Inclination Exchange in a Two-Planet K2-18 System

In this section, we investigate the fundamental process of the orbital momentum exchange between two planets and its manifestations in the secular orbital evolution of planet K2-18 b if the proposed inner planet c is present. A series of numerical simulations was performed using the symplectic integrator WHFast [11], which is part of the REBOUND package (https://rebound.readthedocs.io/en/latest/, accessed on 19 October 2023). This code provides high-speed computation in conservative dynamical systems. Most of the computations were made for 3 × 10 5 yr with 50 steps per period of the inner planet. The basic characteristics of the planet system adopted in our integration are listed in Table 1. The orbital inclination and eccentricity of the two orbits are practically unknown. Specifically, only a very wide upper limit on eccentricity can be derived from the available transit and radial velocity data, which is not helpful for the accurate characterization of orbital dynamics. It is logical to assume that both planets’ eccentricities are “small", while the mutual inclination is larger than a few degrees (from the fact that the inner planet is not transiting). We, therefore, investigate a grid of possible initial configurations. The initial eccentricity of the inner planet is always e 1 ( 0 ) = 0.01 , and the accepted values for the outer planet are e 2 ( 0 ) = 0.01 , 0.05 , 0.1 , 0.2 . Each of these configurations are also integrated on a grid of initial inclinations i 2 ( 0 ) = 5 , 10 , 20 , 40 , thus making a set of 16 different configurations. The initial orientation of the inner orbit is set to i 1 ( 0 ) = 0 , ω 1 ( 0 ) = 0 , and Ω ( 0 ) = 0 by the choice of the fixed non-rotating coordinate frame. We also set ω 2 ( 0 ) = 0 and Ω 2 ( 0 ) = 0 . These Euler angles are found to be rapidly varied on the time scale of 100 yr, and no loss of generality is expected from this choice.
Secular gravitational interactions between the planets and the star cause the orbital elements to vary with time. In the asymptotic limit of small eccentricity and mutual inclination, this behavior is predicted by the classical Laplace–Lagrange linear perturbation theory, which was initially developed for the Solar system dynamics, and specifically, for the Jupiter–Saturn pair [12,13]. Under these conditions, the corresponding model is
e j ( t ) = A 1 + A 2 cos ( 2 π t / p e + ϕ e + π j ) ,
i j ( t ) = B 1 + B 2 cos ( 2 π t / p i + ϕ i + π j ) ,
where j = 1 designates the inner planet, and j = 2 is for the outer planet. We note that in the special case e 1 ( 0 ) = 0 , the equations for e j ( t ) become equivalent to the mirrored sine model with A 1 = A 2 = A e 2 / 2 and ϕ e = 0 , which was used in ref. (Makarov et al., 2023, submitted). The phase shift π j assures that the oscillations of both eccentricities (and both inclinations) are in the opposite direction, so that when one planet reaches its maximum eccentricity, the other planet reaches the minimum value. When both initial eccentricities are nonzero, A 1 A 2 , and the curve looks like a sinusoid with sharper minima. The periods of oscillation p e and p i are equal for both planets. The common phase for the tested configurations is either 0 or π .
Equation (2) provides excellent fits to the osculating orbital elements for all trial initial eccentricity values and for initial inclinations i 2 ( 0 ) up to 20 . An example of the output simulated eccentricities and inclinations for both planets integrated over 30 Kyr is shown in Figure 1. The initial conditions are: e 1 ( 0 ) = 0.01 , e 2 ( 0 ) = 0.05 , i 1 ( 0 ) = 0 , i 2 ( 0 ) = 5 , ω 1 ( 0 ) = ω 2 ( 0 ) = 0 , and Ω 1 ( 0 ) = Ω 2 ( 0 ) = 0 , and only the first 10 Kyr of the data are shown with black dots sampled every 10 yr. The fits (shown with red curves) follow the output data so closely that the black dots are barely seen. The fit for i 1 ( t ) , with its minima close to the initial i 1 ( 0 ) = 0 , degenerates to i 1 ( t ) = 2 B 1 | sin ( 4 π t / p i ) | , as explained above. The full amplitudes of eccentricity and inclination variations are significantly smaller for planet 2 than for planet 1. This is expected from the total angular momentum preservation, and is true for all configurations. The “observed minus calculated” residuals are smaller than the oscillation amplitudes by two orders of magnitude in this example, but they reveal the presence of additional eigenfrequencies that are not captured by the model. These additional eigenfrequencies, however, are sensitive to the initial mutual inclination of the orbits, and become obvious at i 2 ( 0 ) = 20 .
The emerging picture is consistent with previous finds from long-term integrations of other exoplanet systems with period ratios closer to the 2:1 MMR. The main-mode periods p e and p i are weakly dependent on the initial e and i, as long as they are close to zero. The general dependence of p i on i 2 ( 0 ) , however, is exponential. For the configuration in Figure 1, for example, the fitted periods are p e = 4960 yr and p i = 2737 yr. With the same initial conditions, but with i 2 ( 0 ) set to 20 instead of 5 , the periods become 9050 yr and 3206 yr, respectively. Furthermore, if the initial inclination is set to 40 , model (2) is no longer valid for eccentricity, as shown in Figure 2. Both planets’ eccentricities vary non-periodically in a wide range. For the inclination of the outer planet, model (2) still provides an adequate, but imperfect, fit (right panel). The estimated period p i is 4722 yr.

3. Apsidal Precession and Nodal Recession in a Two-Planet K2-18 System

We are not aware of an accurate and tested analytical theory describing the secular and long-period evolution of the node and periapse position in the three-body non-restricted problem with nonzero eccentricities and inclinations. In some papers, e.g., [14], the impact of the inner planet is approximated with an additional quadrupole moment J 2 assigned to the star, and the well-developed Lagrange theory for multipole perturbations is then employed. Essentially, the inner planet is replaced with a stationary gravitating ring of a uniform density with a radius of a 1 . Our numerical experiment for K2-18 with two planets reveals that this approximation gives an incomplete and distorted picture of the complex behavior of the two orbits. Specifically, it predicts that the line of apsides precesses for the outer planet (if mutual inclination is less than arcsin ( 2 / 5 ) ), and the line of nodes recesses. The likely reason for this inadequacy is that the quadrupole approximation does not capture the exchange of angular momentum between the two planets.
From our benchmark simulation with the initial parameters e 1 ( 0 ) = 0.01 , e 2 ( 0 ) = 0.05 , i 1 ( 0 ) = 0 , i 2 ( 0 ) = 5 , ω 1 ( 0 ) = ω 2 ( 0 ) = 0 , and Ω 1 ( 0 ) = Ω 2 ( 0 ) = 0 , we find that the inner planet’s ω 1 ( t ) shows both types of secular perturbation: A linear prograde drift (precession) and a long-period variation. The rate of precession can be expressed as
ω 1 ( t ) = 2 π t / P ω 1 ,
where P ω 1 is the period of periapse circulation. For the tentative inner planet in this configuration, we estimate P ω 1 = 3790 yr, which is significantly shorter than the eccentricity exchange period p e = 4960 yr. Overlaid on this drift, there is a long-period (hereafter, LP) oscillation with a period of about 5000 yr and a full amplitude of ∼1.5 rad. The oscillation is markedly non-sinusoidal. Remarkably, the Ω 1 ( t ) curve does not include any LP variations but shows a steady linear recession with a circulation period P Ω 1 = 5500 yr. Perhaps, the results for the outer planet 2 are more important. The ω 2 ( t ) curve includes both the linear prograde drift and LP components. The estimated periapse circulation period is P ω 1 = 13,000 yr, which is much slower than the inner planet circulation. Figure 3, on the left panel, shows the LP component of ω 2 ( t ) after the subtraction of the linear precession term. The variations are vaguely periodic but non-sinusoidal, with uneven extrema and faster declining slopes than the inclining slopes. Remarkably, the line of nodes’ temporal behavior depicted in the right panel without any manipulation with the output data displays only a non-sinusoidal long-period variation. It looks much more regular and of a stable character than ω 2 ( t ) , and has the same period as the accurately estimated p i = 2737 yr. Thus, the librational parts of secular perturbation for the outer planet have the same period for inclination i, periastron argument ω , and nodal longitude Ω . In stark contradiction to the general predictions of the quadrupole approximation, there is no discernible nodal recession of the outer planet in this configuration.

4. Transit Time Variations in a Two-Planet K2-18 System

The chronology of planetary transits in front of the stellar disk gives us a rare possibility to observe the effects of orbital dynamics. K2-18 b is not known to show any variations in transit times, which are generally more difficult to detect in the extended mission data than in the main Kepler mission data. We approach this problem from the theoretical side, trying to predict the magnitude of possible variations and their time scale. In a two-planet system with the proposed inner companion (planet c, or planet 1 in this paper), the gravitational interaction between the planets results in the variability of all orbital elements (Section 2). Our high-accuracy numerical simulations of the orbital motion confirm that the semimajor axis (a) and, therefore, the mean motion (n), change in a very narrow interval with time. The exchange of orbital momentum between the planets is mostly seen in the periodic variation of eccentricity and inclination, and the more complex behavior of the geometric orientation angles Ω and ω . The former element is the angle between a fixed (inertial) direction of the x Cartesian axis and the ascending node of the orbit. The periapse argument ω is the angle in the instantaneous plane of the orbit between the ascending node and the periapse of the planet. In the classical Lagrange perturbation theory, the variations of each orbital element are assumed to be separable into three independent categories: short-period variations, long-period variations, and secular drift. To approximate the latter two, the complex perturbation equations are averaged over one orbital period. The secular drift is assumed to be linear in time. Our numerical results for a grid of initial configurations confirm that there is no secular drift in n, e, or i. The osculating angles Ω and ω , on the contrary, show both long-period variations, secular drifts, and sudden jumps. To elucidate the nature of these changes, we consider the behavior of the angular orbital momentum vector for planet 2.
For small eccentricity, the magnitude of orbital momentum is
h = G ( M s + M 2 ) a 2 ( 1 e 2 2 ) G ( M s + M 2 ) a 2 ( 1 e 2 2 / 2 ) .
Therefore, in the first-order approximation, h 2 changes periodically with time in a similar manner to e 2 but with the opposite phase. The amplitude of this variation is relatively small, whereas the tilt variation of the vector h is much more substantial. Our fixed inertial reference frame is defined by the orbital plane of planet 1 at time t = 0 . The axes x and y of this frame are in this plane, and the third axis z completes the right-hand triad. Thus, axis z is close to the average direction of the planets’ orbital momenta. The direction cosines of the vector, i.e., the projections of the normalized vector h = h / | | h | are computed as
h x = sin i 2 sin Ω 2 h y = sin i 2 cos Ω 2
from the osculating elements in the output file. The result for the initial configuration with e 1 ( 0 ) = 0.01 , e 2 ( 0 ) = 0.05 , i 1 ( 0 ) = 0 , i 2 ( 0 ) = 5 , ω 1 ( 0 ) = ω 2 ( 0 ) = 0 , and Ω 1 ( 0 ) = Ω 2 ( 0 ) = 0 is shown as a parametric curve for each planet in Figure 4. Both orbits describe complicated folded loops in this projection. The inner planet 1 (red curve) acquires significantly greater orbital inclinations than the outer planet. The picture is strikingly different from a regular nodal recession, which has been assumed in many analyses. A secular recession of the nodes with a constant inclination would show as a concentric circle on this plot. We note that the trajectory of the inner planet crosses the origin several times during the simulated interval of 30 Kyr. In these instances, the inclination of planet 1 nullifies, and the osculating ω and Ω jump by π . This could alternatively be represented by continuous functions allowing the inclination to acquire negative values (which formally violates the definition of the ascending node). Keeping this technical feature in mind, we can now compute the observable transit time variations (TTV) of planet 2.
The commonly used formula for consecutive transit time t k from [15] in the presence of periastron precession is
t k t 0 + P orb k e P orb π cos ω 0 + d ω d k k + O ( e 2 ) ,
where ω 0 is the argument of the periastron at time t 0 , and k is the scaled time (orbit counter) equal to t / P orb . This formula is not valid in our case. First, we note that it is derived as a first-order approximation specifically for the case i = 0 , and the ω appearing there is not the periastron argument but a compound variable u = Ω + ω , i.e., the true longitude of the periastron. It also implicitly assumes that the periastron precesses at a constant rate, while the nodal variation is nil. For a gravitationally perturbed planet, all the orbital elements are functions of time, and computing transit times are not straightforward. A transit time T can be defined as the instance when the planet’s trajectory projected on the fixed { x y } plane, which includes both the vernal equinox and the line of sight, crosses the line of sight. Without a loss of generality, the x axis can be chosen as the line of sight, and the transit equation is then
x ( T ) R 3 ( Ω ) R 1 ( i ) R 3 ( ω ) [ 1 , 0 , 0 ] T [ [ 1 ] ] = 0 ,
where R j ( α ) denotes the right-handed Euler rotation matrix around axis j through angle α , and the time arguments ( T ) of all the osculating elements are omitted for brevity. In the traditional trigonometric form, this equation is equivalent to
1 e 2 sin E ( cos Ω sin ω + cos i cos ω sin Ω ) + ( cos E e ) cos ω cos Ω cos i sin ω sin Ω ) = 0 ,
where E is the eccentric anomaly. Using the output states of the simulations files, which include the values Ω , ω , and i on a regular cadence, the corresponding values E ( T ) can be computed by finding the root of this equation by numerical methods. These values can be converted to the corresponding mean anomaly values by using Kepler’s equation:
M = E e sin E .
Finally, the transit time T is obtained from M ( T ) and the corresponding times of the periastron passage T 0 (also given in the numerical integration output) by
T = T 0 + M / n .
Equation (8) has two roots for each recorded state, because the planet is aligned with the line of sight twice during each revolution, when it transits in front of the stellar disk, and when it is eclipsed behind the star. If the ascending node is defined as the point on the celestial plane where the planet moves toward the observer, the transit time corresponds to the root where y ( T ) is positive.
The resulting curve T ( t ) is a wavy inclined line. It would be a perfectly straight line for an unperturbed single planet. To see the effects of secular gravitational perturbations in the observable intervals between transits, which are called transit time variations (TTV), the linear trend (which is not observable) can be removed by taking the modulo function of the T values:
Δ T = Mod [ T , P orb ] ,
where P orb is the empirically estimated average orbital period in the inertial frame (not to be confused with the often invoked anomalistic period). Thus, the computed periodic part of the TTV for the initial configuration e 1 ( 0 ) = 0.01 , e 2 ( 0 ) = 0.05 , i 1 ( 0 ) = 0 , i 2 ( 0 ) = 5 , ω 1 ( 0 ) = ω 2 ( 0 ) = 0 , and Ω 1 ( 0 ) = Ω 2 ( 0 ) = 0 , which has also been used in Figure 1 and Figure 4, is shown in Figure 5 for the planet of interest (b). The empirically estimated P orb = 32.940562 d is slightly longer than the nominal (unperturbed) orbital period 32.94 d listed in Table 1. At first glance, the result is surprising because we see quasi-sinusoidal variations with a period of ∼4960 yr and a full amplitude of several days. The secular TTV period equals the period of eccentricity libration P e , which is also the period of orbital momentum oscillation. This huge TTV variation, however, cannot be observed. The time span of currently available TTV data, mostly from the Kepler missions, is only a few to several years. The character of secular TTV oscillations is such that the planet spends most of the time on the ascending or descending segments of the curve, which are almost linear. A linear trend in transit times is not detectable, because it is automatically nullified by fitting a slightly adjusted P orb . Only the epochs of the highest curvature are of interest, because they may be detectable. Still, for K2-18 c, the TTV measurements should span > 300 yr for the curvature around the extrema of the transit time curve to become detectable at the current level of measurement precision (approximately, 20 min).
In terms of measurable effects, which might be useful to confirm the presence of a perturbing inner planet, the short-period variations of transit times (and the orbital momentum magnitude) are more interesting. These are smaller by a few orders of magnitude than the LP oscillations, but their time scale is down to ∼10 yr. To elucidate the issue of TTV detectability, we performed a special integration for 30 Kyr with state dumps every year and the same initial parameters. The output files allow us to look closer at the character of short-term variations in transit times. We processed the simulated data in the same way as described above. The result for a small section of the curve spanning 30 yr is shown in Figure 6. The linear trend is not measurable, because it would be absorbed in the fitted average transit period. The remaining variations look non-periodic and have an amplitude up to 10 min. They should be present in the TTV measurement data, but it is not easy to recognize them because of the absence of a specific pattern. Most likely, such variations would be interpreted as excessive observational noise that is not captured by the formal uncertainties.

5. Tidal Orbital Damping in a Single-Planet K2-18 System

If the inner planet c (planet 1 in this paper) is a false positive, and K2-18 is a single-planet system, the possible effects of the orbital perturbation scale down in magnitude and variety. A Keplerian orbit with fixed orbital elements is then a very good approximation, and the secular effects are caused by the finite size of the gravitating bodies and their rotation. The prevalence of low eccentricity orbits for detected exoplanets is in stark disparity with the statistics for binary stars, where a flat distribution of eccentricity is determined. It is widely assumed that the tidal dissipation mechanisms are responsible [16]. As the planet revolves around the star, the local gradient of the gravitational potential within the planet’s body creates a transient distortion (asphericity) of its shape, which is called the tidal bulge. A stationary tidal bulge is realized when the perturber (the star in this case) is motionless with respect to the planet’s surface. In a more realistic dynamical model, the bulge moves across the surface because of the combined effect of the perturber’s orbital motion and the planet’s rotation. The internal tidal friction is a nonlinear function of the relative angular velocity, and the resulting polar torque on the planet can be formally represented with an infinite series of Fourier harmonics in the more advanced Maxwell and Andrade rheological models [17,18,19]. As a mathematical abstraction, this can be viewed as a number of “waves" running around the circumference with diminishing amplitudes and harmonics of the main tidal frequency. The physical reality is a single tidal deformation, which is mostly an ellipsoidal shape, with overlaid periodically-changing higher-order spherical harmonics.
For simplicity, the perturber is considered to be a point mass, and the tidal torque acting on the planet is reciprocated by an orbital torque acting on the perturber. To complete the model, the tides raised by the planet on the star should be considered, in which case the planet is assumed to be a point mass. The total orbital action is the sum of the torques generated by the planetary and stellar tidal bulges. The differential equations of orbital evolution are derived in the framework of the Lagrange perturbation theory by averaging out the disturbing function over one orbital period. Thus, only the secular and long-period can be obtained in principle, but an additional averaging over one circulation period of the periastron leaves only the secular drifts. With these caveats and limitations, the resulting equations for the tides raised by the star on the planet are [20]:
d a d t n M s M 2 R 2 5 a 4 [ 3 ( 1 5 e 2 + 63 8 e 4 ) K 2 ( 2 n 2 f 2 ) + 3 8 e 2 ( 1 1 4 e 2 ) K 2 ( n 2 f 2 ) + 9 4 e 2 ( 1 + 9 4 e 2 ) K 2 ( n ) + 81 8 e 4 K 2 ( 2 n ) + 441 8 e 2 ( 1 123 28 e 2 ) K 2 ( 3 n 2 f 2 ) + 867 2 e 4 K 2 ( 4 n 2 f 2 ) ] ,
d e d t n e M s M 2 R 2 5 a 5 [ 3 16 ( 1 1 4 e 2 ) K 2 ( n 2 f 2 ) 3 4 ( 1 21 4 e 2 ) K 2 ( 2 n 2 f 2 ) + 9 8 ( 1 + 5 4 e 2 ) K 2 ( n ) + 81 16 e 2 K 2 ( 2 n ) + 147 16 ( 1 179 28 e 2 ) K 2 ( 3 n 2 f 2 ) + 867 8 e 2 K 2 ( 4 n 2 f 2 ) ] ,
where f 2 is the rotational frequency of planet 2, R 2 is its radius, K 2 is its frequency-dependent tidal quality function (kvalitet), n is the orbital mean motion, M s is the mass of the star, and the other parameters are as previously defined. The contributions from the tides raised by the planet on the star can be obtained by swapping the subscripts 2 and s everywhere. These are approximate equations that ignore all terms O ( ϵ 2 ) , O ( e 6 ) , and higher orders of obliquity ϵ and eccentricity, also taking into account only the main (quadrupole) term of the perturbing potential. Even with these simplifying approximations, the resulting secular trend is quite uncertain in general, because the parameters e, f, and K ( ν ) are practically unknown. It is safe, however, to neglect the contribution from the star, because the kvalitet values for old M-dwarfs are expected to be very low. For gas giant planets, these values can be as high as 10 5 10 4 , and even larger for semi-molten terrestrial compositions [21].
The existing models for K 2 ( ν ) suggest complex functions with anti-symmetric kinks at the main spin–orbit resonances 1:1, 3:2, etc. In the better theoretically developed Maxwell rheological model, the amplitude of the kinks and their position depends on the so-called Maxwell time, which is the ratio of the effective viscosity and the unrelaxed shear modulus [22,23]. The driving parameter is the viscosity of the energy-dissipating layer of the planet. While this parameter is estimated to be in the range of 10 4 10 2 Pa s for warm, pressurized water in the interiors of icy giants [24], pressurized ice has values of the order 10 15 Pa s. The steep dependence of the effective viscosity on temperature for silicates also defines a great dynamical range of Maxwell times, and, therefore, of tidal qualities, for planets of a terrestrial composition. Planet K2-18 b is presumed to have an Earth-like core covered by a massive water ocean. The ocean is unlikely to be of significance in dissipating the tidal energy, and most of the action should be allocated to the terrestrial core. Given these prerequisites, we can consider three distinct options for the tidal spin–orbit interaction.
The first option is for a rapidly rotating planet’s core of terrestrial rheology. If the solid body of K2-18 b rotates as fast as the Earth, f 2 > > n , we can omit the multiples of n in the tidal mode arguments of K 2 , which obtains, also omitting the relatively small term proportional to K 2 ( 2 n ) :
d a d t ( fast ) n M s M 2 R 2 5 a 4 [ 3 8 ( 8 + 108 e 2 + 573 e 4 ) K 2 ( 2 f 2 ) 9 4 e 2 ( 1 + 9 4 e 2 ) K 2 ( n ) ] ,
d e d t ( fast ) n e M s M 2 R 2 5 a 5 [ 3 4 ( 11 + 505 e 2 ) K 2 ( 2 f 2 ) 9 8 ( 1 + 5 4 e 2 ) K 2 ( n ) ] .
Even the sign of these derivatives is uncertain, because for sufficiently cold and inviscid planets, K 2 ( n ) > K 2 ( 2 f 2 ) . The closest analogy of such systems is the Earth–Moon system, which is expanding both in a and e. This is not, however, always the case, even for such cold and inviscid planets as Earth. When both frequencies n and 2 f 2 are far from the main resonances and the tidal quality peaks, K 2 ( ν ) is approximately inversely proportional to frequency. In this case,
K 2 ( n ) / K 2 ( 2 f 2 ) 2 f 2 / n .
A simple calculation reveals that if the rotational period of K2-18 b equals 1 d, for example, d a / d t becomes negative for e > 0.165 . Orbital decay for such rapidly rotating planets is precluded for any eccentricity only if f 2 15 n . The relative contribution of the two terms in Equation (14) is of the same order of maginitude, because they both are O ( e 0 ) , and we find that the planet’s orbit can circularize ( d e / d t < 0 ) at small eccentricities if f 2 4 n .
The Maxwell time and kvalitet values are quite uncertain for K2-18 b. We are using the K 2 ( n ) = 10 4 value as a benchmark, which is the upper envelope for astrometrically estimated values for Saturn and Jupiter. The computed derivatives are then, with the other parameters listed in Table 1, ( d a / d t ) ( fast ) 3 × 10 7 m d 1 and ( d e / d t ) ( fast ) 0 for e = 0 , and ( d a / d t ) ( fast ) 3 × 10 7 m d 1 and ( d e / d t ) ( fast ) 3 × 10 18 m d 1 for e = 0.1 . Thus, a secular expansion of the orbit and a slow circularization are possible at low eccentricity values, but the rates are marginally low, and it is therefore unlikely that the tidal torques play any role in the habitability conditions of the single K2-18 planet.
The second option considered in this paper is a 3:2 spin–orbit resonance. For planets of terrestrial compositions that are close to their host stars and have a finite orbital eccentricity, this is the more probable equilibrium state than the synchronous rotation (1:1 resonance). The closest example of a planet in a 3:2 spin–orbit resonance is Mercury [25], although temporary captures into other resonances are possible over the course of Mercury’s tumultuous dynamical history [26]. Higher chances of capture into the higher-order resonances are associated with larger orbital eccentricities, but the triaxial elongation of the permanent shape is also an important factor [27,28]. The differential equations of orbital evolution for the 3:2 resonance, without any additional omissions, become
d a d t ( 3 : 2 ) n M s M 2 R 2 5 a 4 [ ( 3 69 4 e 2 6639 16 e 4 ) K 2 ( n ) + 3 32 e 2 ( 4 109 e 2 ) K 2 ( 2 n ) ] ,
d e d t ( 3 : 2 ) n e M s M 2 R 2 5 a 5 [ 3 32 ( 20 + 1129 e 2 ) K 2 ( n ) + 3 64 ( 4 + 107 e 2 ) K 2 ( 2 n ) ] .
The orbit decays ( d a / d t is negative) for e 0.26 and expands for a smaller eccentricity. The orbital eccentricity always decreases, irrespective of its current value. The estimated rates for a cold terrestrial planet with K 2 ( n ) = 10 4 are too low to be of consequence. We estimate ( d a / d t ) ( 3 : 2 ) 4.5 × 10 6 m d 1 and ( d e / d t ) ( 3 : 2 ) 2.3 × 10 17 d 1 for e = 0.1 , and ( d a / d t ) ( 3 : 2 ) 2.7 × 10 6 m d 1 and ( d e / d t ) ( 3 : 2 ) 9.6 × 10 17 m d 1 for e = 0.2 .
Finally, the third option is a synchronous rotation, i.e., the 1:1 spin–orbit resonance. This equilibrium state is often assumed in the literature for habitable exoplanets. However, this equilibrium state corresponding to the global minimum of potential energy should not be taken for granted for planets of terrestrial composition with significant triaxiality of the inertia tensors [27]. In particular, it has been shown that higher-order spin–orbit resonances are more probable for specific two-planet systems with habitable planets [29]. Still, for the sake of consistency, we provide the approximate equations for this case:
d a d t ( 1 : 1 ) n M s M 2 R 2 5 a 4 [ 3 e 2 ( 19 + 79 e 2 ) K 2 ( n ) 3549 8 e 4 K 2 ( 2 n ) ] ,
d e d t ( 1 : 1 ) n e M s M 2 R 2 5 a 5 [ 3 8 ( 28 + 153 e 2 ) K 2 ( n ) 1815 16 e 2 K 2 ( 2 n ) ] .
Note that the largest terms for d a ( 1 : 1 ) / d t are O ( e 2 ) , which almost nullifies the orbit size evolution at a small or moderate eccentricity. The sign of d e ( 1 : 1 ) / d t changes from minus to plus at e 0.4 , but our formulae become inaccurate at such eccentricity values, and the 1:1 synchronization is hardly possible for terrestrial planets anyway. The estimated rates of eccentricity decay are still vanishingly low, amounting to ∼10 16 d 1 .
More complicated scenarios of tidal evolution are possible if planet K2-18 b is covered with a deep global water ocean. The liquid outer layer does not have a permanent figure, and it is likely to be rotating pseudosynchronously [23], if it is mechanically decoupled from the rocky bottom. The rocky body, on the other hand, is either synchronized or rotates faster than the ocean. Depending on the coefficient of friction between the ocean and the rocky core, a significant additional mechanism of energy dissipation can emerge. Similar models, but with a reversed structure, including a molten core and a solid outer shell, have been investigated for Mercury [30].
Our conclusion is that, irrespective of the rotational state of the planet, tidal dissipation can hardly be significant in shaping the orbit in the long term. If the eccentricity is indeed very small, there should be other reasons for that. The dynamically cool population orbiting M-dwarfs with small eccentricities [31] can be primordial, i.e., inherited from the early stages of migration in the protoplanetary disk. This result may change if the inner planet in K2-18 is real. Because of the explicit n ( R / a ) 5 R 5 / a 13 / 2 dependence of secular rates in Equation (12), the rates can be higher by 2–3 orders of magnitude for the inner companion. This makes the inner planets relatively efficient sinks of the orbital energy and angular momentum. The eccentricity exchange mechanism (Section 2) prevents circularization of the inner planet as long as the outer planet remains eccentric. Therefore, the tidal damping of the inner planet is effectively transferred to the outer planet as well. This makes the characteristic e-folding times somewhat longer, but they may be comparable with the lifetimes of M-dwarfs.

6. Secular Perturbations Caused by Stellar Rotation

The equilibrium shape of a rotating star is an axially symmetric oblate figure. The gravitational potential outside the stellar surface is formally represented by an expansion in spherical harmonics, where the lowest-order zonal (quadrupole) term is by far the greatest. The disturbing function generated by this deviation from spherical symmetry is proportional to the corresponding coefficient J 2 , which is the quadrupole moment. Ignoring the small difference between the total mass M s + M 2 and M s , the equation is
J 2 = k s 3 f s n 2 R s a 3 ,
where k s is the corresponding Love number and f s is the frequency of the stellar rotation. Using the nominal values in Table 1 and assuming P rot = 8.8 d, J 2 1.6 × 10 5 k s is estimated for the star K2-18. The Love number is quite uncertain, but we can assume it to be of the order of 0.1, which would bring J 2 to approximately 1.6 × 10 6 , which is almost eight times the time-average J 2 value of the Sun [32].
The expansions for the disturbing functions and Lagrange planetary equations have been worked out up to O ( e 4 ) and the J 6 zonal harmonics for the significantly more perturbed geopotential [33,34]. As usual, the orbital effects are supposed to be separable into long-period and secular terms. For the latter, limiting our consideration only to the leading terms, the equations are:
d Ω d t ( sec ) 3 n 2 p 2 cos ϵ J 2 d ω d t ( sec ) 3 n p 2 1 5 4 sin 2 ϵ J 2 d M d t ( sec ) n + 3 n 2 p 2 1 e 2 1 3 2 sin 2 ϵ J 2 ,
where ϵ is the obliquity of the planet’s orbit on the equator of the star, and p = ( 1 e 2 ) a / R s . Note that at obliquity values below 63.43 , the sign of ( d ω / d t ) ( sec ) is positive, which signifies a constant-rate precession of the line of apsides. The sign of ( d Ω / d t ) ( sec ) , however, is always negative. Thus, stellar rotation causes a nodal recession. The perturbation of the mean anomaly rate is a constant offset reflecting the additional gravitational acceleration of the planet in the radial direction, which is positive (faster mean motion) for ϵ < 54.74 . As we will see in the following, the equation for ( d ω / d t ) ( sec ) is of special importance, which, after substituting the formula for J 2 , becomes
d ω d t ( sec ) k s n ( 1 e 2 ) 2 f s n 2 R s a 5 ( 1 5 4 sin 2 ϵ ) .
This equation differs by a factor of two from the analogous formulae given in [35,36]. In these papers, the marginal case of ϵ = 0 is considered, and the quoted equations are given not for the periastron argument ω but for the compound variable u = Ω + ω . This subtle detail is crucially important for the general consideration with ϵ 0 . There are no secular drifts of other orbital elements.
Long-period variations emerge for all the orbital elements. The approximate equations are obtained by the integration of periodic terms of the corresponding derivatives over time. All these relations are proportional to e J 2 except the two most important elements in the context of habitability conditions, eccentricity, and the periastron argument. Limiting the expansion to the dominant (for small e) main harmonics in ω , the appropriate equations are:
Δ ω ( LP ) J 2 8 e p 2 ( 12 25 e 2 ) + 8 ( 1 + 4 e 2 ) I sin ω Δ e ( LP ) J 2 8 p 2 3 ( 4 + e 2 ) 8 ( 1 e 2 ) I cos ω ,
where I = sin 2 ϵ . In the right-hand parts of these equations, ω includes only the secular motion term. Note that the amplitude of ω oscillations can be much larger than the constant-rate precession over one precession period for e < 0.5 . At a very small e and ϵ , the LP variation becomes the dominant term. Its amplitude can be arbitrarily large for e 0 . This, however, does not make the LP perturbation more detectable, because the observable TTV effect is proportional to e.
Assuming J 2 = 1.6 × 10 6 for the star K2-18, we estimate the period of periapse precession driven by the rotational deformation to be approximately equal to 84 Myr. Such small effects cannot be measured today. The amplitude of LP oscillation of ω is a meager 5.4 × 10 7 rad for e 2 = 0.001 , and 1.1 × 10 7 rad for e 2 = 0.005 . The amplitude of the eccentricity oscillations is also negligibly small. We conclude that the orbit of planet K2-18 b should be extremely stable and nearly constant, if it is the only planet in the system.

7. Milanković Cycles of Habitable Exoplanets

Gravitational and tidal interactions of Earth with the Sun, the Moon, and other planets produce cyclic and secular perturbations of the orbit and the rotation axis, which modulate the climate (https://climate.nasa.gov/news/2948/milankovitch-orbital-cycles-and-their-role-in-earths-climate/, accessed on 6 October 2023). The fastest cycle with a period of ∼26 Kyr is related to the nodal precession caused by the gravitational pull of the Sun and the Moon on the oblate shape of Earth with a significant J 2 coefficient. The climate impact is limited to the time of the perihelion passage in the year. It combines with the much slower apsidal precession d ω / d t ( sec ) to produce a constant drift of the seasons within the year with a period of ∼21 Kyr, as well as the modulation of the relative duration of the seasons. The cold and warm seasons are equal in duration when ω = 0 or π . Perhaps of greater importance are the secular changes in the obliquity of the Earth’s equator on the ecliptic in the range 2.4 . A larger obliquity causes greater seasonal variations of insolation and temperature. The period of this oscillation is ∼41 Kyr. Finally, as the Earth’s orbit interacts with the other planets’ orbits, the true inclination with respect to the invariable plane oscillates with a period of approximately 100 Kyr. The climate impact is similar to the long-period variation of the obliquity. Although this effect is believed to be relatively small, the timescale is conspicuously close to the main cycle of glaciation periods [37]. The original Milanković model includes the oscillations of orbital eccentricity (of similar periods), and indeed, both numerical integrations and the Laplace–Lagrange theory show that Earth’s eccentricity varies between 0 and 0.06 and its inclination varies between 0 and 2.5 in a quasi-cyclic manner with characteristic periods close to 100 Kyr [38]. The shape of each cycle is remarkably similar to the curves in our Figure 1 for the inner planet K2-18 c, indicating that the same basic dynamics are at work. Empirically, these relatively small variations of orbital inclination and eccentricity emerge as the most important factors driving the ice ages on Earth.
Similar processes take place in systems with two or more exoplanets [39]. As far as the system K2-18 is concerned, the options with Milanković-like cycles are quite different for the single-planet and two-planet configurations. If planet b is a lone planet, both the tidal perturbation of the planet and the secular perturbations caused by the rotation of the star are too slow and too small to be of significance (Section 5 and Section 6). The remaining important parameters, which are also the most uncertain, are the rate of the rotation of the planet and its obliquity on the orbital plane. The planet may be covered by a deep global water ocean, which is most likely to rotate pseudosynchronously, i.e., slightly faster than the mean orbital motion. The rocky core may be dynamically decoupled from the ocean, and its rotation may initially be different, including a finite obliquity. The force of friction at the bottom of the ocean is probably sufficient to equalize this differential rotation over the age of the system. However, the core has a permanent figure with triaxial elongation, which precludes the pseudosynchronous equilibrium. It can remain captured in a 3:2 or 1:1 spin–orbit synchronicity. The emerging picture is a constant circulation of the global ocean superimposed with periodic flows caused by the forced librations of the solid core. In the context of climate changes, a water-world planet is unlikely to have seasons, and it would probably have permanent surface ice covering the polar caps, but the insolation flux and temperature is more uniformly distributed around the equatorial zone by the slow residual rotation of the ocean with respect to the direction to the star.
Possible scenarios of Milanković-like cycles in a two-planet K2-18 systems are rich and diverse. Apart from the secular nodal recession and apsidal precession, relatively fast periodic variations of the habitable planet b are predicted in eccentricity and inclination. Their characteristic periods are a few thousand years. The tidal damping of eccentricity is somewhat more efficient because of the constant eccentricity exchange between the planets, and the much higher dissipation rate for the inner planet 1. The constantly changing inclination precludes an equilibrium state with zero obliquity, which results in rapidly changing seasons, smaller ice caps around the poles, and a climate that is more temperate and amenable to life. If the amplitude of the eccentricity cycles is significant, the outer planet’s ocean will be warmed up and cooled down every 5000 yr. Furthermore, since the rate of pseudosynchronous rotation is defined by eccentricity, the long-term oscillation of the global ocean rotation and the current patterns should be expected.

8. Summary

We have investigated the dynamical evolution of the K2-18 exoplanet system using numerical and analytical methods in two possible configurations: A single planet b (detected via transits in front of the stellar disk), and a two-planet system including the suggested RV-detected planet c. Our main conclusions are:
  • In the two-planet option, a relatively rapid exchange of angular momentum, eccentricity, and inclination between the planet is found with a characteristic cycle period of a few Kyr. The integrated trajectories are in excellent agreement with the first-order Lagrange perturbation theory developed for the Solar system planets.
  • In the two-planet option, the lines of apsides of both planets are subject to both precession and long-period variations. The periods of the periodic components are non-commensurate. The line of nodes of the outer (habitable) planet undergoes long-period oscillations, but no significant secular drift.
  • In the two-planet option, the vectors of orbital momenta describe complicated trajectories in the inertial frame because of the simultaneously varying inclination and nodes. The resulting secular TTV changes are large in amplitude (several days peak-to-peak) but so slow that they cannot be detected. The short-term TTV effects on the time scale of 1 year are marginally detectable, but they lack a specific pattern.
  • In the single planet option, the tidal damping of the eccentricity and orbit size is found to be negligibly slow for three possible rotational states: viz., a synchronous rotation—a 3:2 spin–orbit resonance—and a fast or retrograde rotation. The estimated characteristic e-folding times are longer than the stellar age.
  • The main J 2 moment of the star K2-18 due to its sidereal rotation is estimated at 1.6 × 10 6 . It turns out to be too small to excite a secular nodal recession or periastron precession that could be measurable. The amplitude of long-period variations is very small too.
  • Milancović-like variations of planet b’s climate are not expected in the single planet mode, but habitability conditions can be harsh, unless a significant obliquity is maintained. A rich variety of climatic changes emerges in the two-planet option with a rapid succession of seasons, constant circulation of the global ocean, a more temperate equatorial zone, and a long-period partial glaciation.

Author Contributions

Conceptualization, V.V.M. and A.G.; methodology, V.V.M.; software, A.G.; validation, V.V.M.; formal analysis, V.V.M. and A.G.; investigation, V.V.M. and A.G.; resources, A.G.; writing—original draft preparation, V.V.M.; writing—review and editing, V.V.M. and A.G.; visualization, V.V.M.; supervision, V.V.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The simulation inputs and results discussed in this paper are available upon reasonable request to the corresponding author.

Conflicts of Interest

Author Alexey Goldin was employed by the company Teza Technology. The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

  1. Foreman-Mackey, D.; Montet, B.T.; Hogg, D.W.; Morton, T.D.; Wang, D.; Schölkopf, B. A Systematic Search for Transiting Planets in the K2 Data. Astrophys. J. 2015, 806, 215. [Google Scholar] [CrossRef]
  2. Montet, B.T.; Morton, T.D.; Foreman-Mackey, D.; Johnson, J.A.; Hogg, D.W.; Bowler, B.P.; Latham, D.W.; Bieryla, A.; Mann, A.W. Stellar and Planetary Properties of K2 Campaign 1 Candidates and Validation of 17 Planets, Including a Planet Receiving Earth-like Insolation. Astrophys. J. 2015, 809, 25. [Google Scholar] [CrossRef]
  3. Sarkis, P.; Henning, T.; Kürster, M.; Trifonov, T.; Zechmeister, M.; Tal-Or, L.; Anglada-Escudé, G.; Hatzes, A.P.; Lafarga, M.; Dreizler, S.; et al. The CARMENES Search for Exoplanets around M Dwarfs: A Low-mass Planet in the Temperate Zone of the Nearby K2-18. Astrophys. J. 2018, 155, 257. [Google Scholar] [CrossRef]
  4. Cloutier, R.; Astudillo-Defru, N.; Doyon, R.; Bonfils, X.; Almenara, J.M.; Benneke, B.; Bouchy, F.; Delfosse, X.; Ehrenreich, D.; Forveille, T.; et al. Characterization of the K2-18 multi-planetary system with HARPS. A habitable zone super-Earth and discovery of a second, warm super-Earth on a non-coplanar orbit. Astron. Astrophys. 2017, 608, A35. [Google Scholar] [CrossRef]
  5. Benneke, B.; Wong, I.; Piaulet, C.; Knutson, H.A.; Lothringer, J.; Morley, C.V.; Crossfield, I.J.M.; Gao, P.; Greene, T.P.; Dressing, C.; et al. Water Vapor and Clouds on the Habitable-zone Sub-Neptune Exoplanet K2-18b. Astrophys. J. 2019, 887, L14. [Google Scholar] [CrossRef]
  6. Benneke, B.; Werner, M.; Petigura, E.; Knutson, H.; Dressing, C.; Crossfield, I.J.M.; Schlieder, J.E.; Livingston, J.; Beichman, C.; Christiansen, J.; et al. Spitzer Observations Confirm and Rescue the Habitable-zone Super-Earth K2-18b for Future Characterization. Astrophys. J. 2017, 834, 187. [Google Scholar] [CrossRef]
  7. Cloutier, R.; Astudillo-Defru, N.; Doyon, R.; Bonfils, X.; Almenara, J.M.; Bouchy, F.; Delfosse, X.; Forveille, T.; Lovis, C.; Mayor, M.; et al. Confirmation of the radial velocity super-Earth K2-18c with HARPS and CARMENES. Astron. Astrophys. 2019, 621, A49. [Google Scholar] [CrossRef]
  8. Radica, M.; Artigau, É.; Lafreniére, D.; Cadieux, C.; Cook, N.J.; Doyon, R.; Amado, P.J.; Caballero, J.A.; Henning, T.; Quirrenbach, A.; et al. Revisiting radial velocity measurements of the K2-18 system with the line-by-line framework. Mon. Not. R. Astron. Soc. 2022, 517, 5050–5062. [Google Scholar] [CrossRef]
  9. Reiners, A.; Zechmeister, M.; Caballero, J.A.; Ribas, I.; Morales, J.C.; Jeffers, S.V.; Schöfer, P.; Tal-Or, L.; Quirrenbach, A.; Amado, P.J.; et al. The CARMENES search for exoplanets around M dwarfs. High-resolution optical and near-infrared spectroscopy of 324 survey stars. Astron. Astrophys. 2018, 612, A49. [Google Scholar] [CrossRef]
  10. Makarov, V.V.; Beichman, C.A.; Catanzarite, J.H.; Fischer, D.A.; Lebreton, J.; Malbet, F.; Shao, M. Starspot Jitter in Photometry, Astrometry, and Radial Velocity Measurements. Astrophys. J. 2009, 707, L73–L76. [Google Scholar] [CrossRef]
  11. Rein, H.; Tamayo, D. WHFAST: A fast and unbiased implementation of a symplectic Wisdom-Holman integrator for long-term gravitational simulations. Mon. Not. R. Astron. Soc. 2015, 452, 376–388. [Google Scholar] [CrossRef]
  12. Brouwer, D.; van Woerkom, A.J.J. The Secular Variations of the Orbital Elements of the Principal Planets; Astronomical Papers Prepared for the Use of the American Ephemeris and Nautical Almanac; US Government Printing Office: Washington, DC, USA, 1950; Volume 13, pp. 81–107.
  13. Murray, C.D.; Dermott, S.F. Solar System Dynamics; Cambridge University Press: Cambridge, UK, 1999. [Google Scholar]
  14. Miralda-Escudé, J. Orbital Perturbations of Transiting Planets: A Possible Method to Measure Stellar Quadrupoles and to Detect Earth-Mass Planets. Astrophys. J. 2002, 564, 1019–1023. [Google Scholar] [CrossRef]
  15. Giménez, A.; Bastero, M. A Revision of the Ephemeris-Curve Equations for Eclipsing Binaries with Apsidal Motion. Astrophys. Space Sci. 1995, 226, 99–107. [Google Scholar] [CrossRef]
  16. Pont, F.; Husnoo, N.; Mazeh, T.; Fabrycky, D. Determining eccentricities of transiting planets: A divide in the mass-period plane. Mon. Not. R. Astron. Soc. 2011, 414, 1278–1284. [Google Scholar] [CrossRef]
  17. Kaula, W.M. Tidal Dissipation by Solid Friction and the Resulting Orbital Evolution. Rev. Geophys. Space Phys. 1964, 2, 661–685. [Google Scholar] [CrossRef]
  18. Ferraz-Mello, S.; Rodríguez, A.; Hussmann, H. Tidal friction in close-in satellites and exoplanets: The Darwin theory re-visited. Celest. Mech. Dyn. Astron. 2008, 101, 171–201. [Google Scholar] [CrossRef]
  19. Efroimsky, M.; Williams, J.G. Tidal torques: A critical review of some techniques. Celest. Mech. Dyn. Astron. 2009, 104, 257–289. [Google Scholar] [CrossRef]
  20. Boué, G.; Efroimsky, M. Tidal evolution of the Keplerian elements. Celest. Mech. Dyn. Astron. 2019, 131, 30. [Google Scholar] [CrossRef]
  21. Makarov, V.V.; Berghea, C.T.; Efroimsky, M. Spin-orbital Tidal Dynamics and Tidal Heating in the TRAPPIST-1 Multiplanet System. Astrophys. J. 2018, 857, 142. [Google Scholar] [CrossRef]
  22. Efroimsky, M. Bodily tides near spin-orbit resonances. Celest. Mech. Dyn. Astron. 2012, 112, 283–330. [Google Scholar] [CrossRef]
  23. Makarov, V.V. Equilibrium Rotation of Semiliquid Exoplanets and Satellites. Astrophys. J. 2015, 810, 12. [Google Scholar] [CrossRef]
  24. French, M.; Nettelmann, N. Viscosity and Prandtl Number of Warm Dense Water as in Ice Giant Planets. Astrophys. J. 2019, 881, 81. [Google Scholar] [CrossRef]
  25. Pettengill, G.H.; Dyce, R.B. A Radar Determination of the Rotation of the Planet Mercury. Nature 1965, 206, 1240. [Google Scholar] [CrossRef]
  26. Noyelles, B.; Frouard, J.; Makarov, V.V.; Efroimsky, M. Spin-orbit evolution of Mercury revisited. Icarus 2014, 241, 26–44. [Google Scholar] [CrossRef]
  27. Goldreich, P.; Peale, S. Spin-orbit coupling in the solar system. Astron. J. 1966, 71, 425. [Google Scholar] [CrossRef]
  28. Makarov, V.V. Conditions of Passage and Entrapment of Terrestrial Planets in Spin-orbit Resonances. Astrophys. J. 2012, 752, 73. [Google Scholar] [CrossRef]
  29. Makarov, V.V.; Berghea, C. Dynamical Evolution and Spin-orbit Resonances of Potentially Habitable Exoplanets. The Case of GJ 667C. Astrophys. J. 2014, 780, 124. [Google Scholar] [CrossRef]
  30. Margot, J.L.; Peale, S.J.; Jurgens, R.F.; Slade, M.A.; Holin, I.V. Large Longitude Libration of Mercury Reveals a Molten Core. Science 2007, 316, 710. [Google Scholar] [CrossRef]
  31. Sagear, S.; Ballard, S. The orbital eccentricity distribution of planets orbiting M dwarfs. Proc. Natl. Acad. Sci. USA 2023, 120, e2217398120. [Google Scholar] [CrossRef]
  32. Eren, S.; Rozelot, J.P. Exploring the Temporal Variation of the Solar Quadrupole Moment J2. Astrophys. J. 2023, 942, 90. [Google Scholar] [CrossRef]
  33. Brouwer, D.; Clemence, G.M. Methods of Celestial Mechanics; Elsevier: Amsterdam, The Netherlands, 1961. [Google Scholar]
  34. Merson, R.H. The Motion of a Satellite in an Axi-symmetric Gravitational Field. Geophys. J. 1961, 4, 17–52. [Google Scholar] [CrossRef]
  35. Sterne, T.E. Apsidal motion in binary stars. Mon. Not. R. Astron. Soc. 1939, 99, 451–462. [Google Scholar] [CrossRef]
  36. Ragozzine, D.; Wolf, A.S. Probing the Interiors of very Hot Jupiters Using Transit Light Curves. Astrophys. J. 2009, 698, 1778–1794. [Google Scholar] [CrossRef]
  37. Hays, J.D.; Imbrie, J.; Shackleton, N.J. Variations in the Earth’s Orbit: Pacemaker of the Ice Ages. Science 1976, 194, 1121–1132. [Google Scholar] [CrossRef]
  38. Mogavero, F.; Laskar, J. Long-term dynamics of the inner planets in the Solar System. Astron. Astrophys. 2021, 655, A1. [Google Scholar] [CrossRef]
  39. Ershkov, S.; Levchenko, D.; Prosviryakov, E.Y. A novel type of ER3BP introducing Milankovitch cycles or seasonal irradiation processes influencing onto orbit of planet. Arch. Appl. Mech. 2023, 93, 822. [Google Scholar] [CrossRef]
Figure 1. Secular evolution of K2-18 planets in a two-planet configuration. Numerical integration of orbital motion was performed for 30 Kyr with the initial conditions: e 1 ( 0 ) = 0.01 , e 2 ( 0 ) = 0.05 , i 1 ( 0 ) = 0 , i 2 ( 0 ) = 5 , ω 1 ( 0 ) = ω 2 ( 0 ) = 0 , and Ω 1 ( 0 ) = Ω 2 ( 0 ) = 0 . The black curves show the actual output of the numerical simulations. The red curves show the best fits with the theoretical model (2). Upper left panel: inner planet K2-18 c eccentricity. Upper right panel: inner planet K2-18 c inclination in degrees. Lower left panel: outer planet K2-18 b eccentricity. Lower right panel: outer planet K2-18 b inclination in degrees.
Figure 1. Secular evolution of K2-18 planets in a two-planet configuration. Numerical integration of orbital motion was performed for 30 Kyr with the initial conditions: e 1 ( 0 ) = 0.01 , e 2 ( 0 ) = 0.05 , i 1 ( 0 ) = 0 , i 2 ( 0 ) = 5 , ω 1 ( 0 ) = ω 2 ( 0 ) = 0 , and Ω 1 ( 0 ) = Ω 2 ( 0 ) = 0 . The black curves show the actual output of the numerical simulations. The red curves show the best fits with the theoretical model (2). Upper left panel: inner planet K2-18 c eccentricity. Upper right panel: inner planet K2-18 c inclination in degrees. Lower left panel: outer planet K2-18 b eccentricity. Lower right panel: outer planet K2-18 b inclination in degrees.
Universe 09 00463 g001
Figure 2. Secular evolution of K2-18 planets in a two-planet configuration. The numerical integration of orbital motion was performed for 30 Kyr with the initial conditions: e 1 ( 0 ) = 0.01 , e 2 ( 0 ) = 0.05 , i 1 ( 0 ) = 0 , i 2 ( 0 ) = 40 , ω 1 ( 0 ) = ω 2 ( 0 ) = 0 , and Ω 1 ( 0 ) = Ω 2 ( 0 ) = 0 . The black curves show the actual output of the numerical simulations. The red curve show the best fits with the theoretical model (2). Left panel: outer planet K2-18 b eccentricity. Right panel: outer planet K2-18 b inclination in degrees.
Figure 2. Secular evolution of K2-18 planets in a two-planet configuration. The numerical integration of orbital motion was performed for 30 Kyr with the initial conditions: e 1 ( 0 ) = 0.01 , e 2 ( 0 ) = 0.05 , i 1 ( 0 ) = 0 , i 2 ( 0 ) = 40 , ω 1 ( 0 ) = ω 2 ( 0 ) = 0 , and Ω 1 ( 0 ) = Ω 2 ( 0 ) = 0 . The black curves show the actual output of the numerical simulations. The red curve show the best fits with the theoretical model (2). Left panel: outer planet K2-18 b eccentricity. Right panel: outer planet K2-18 b inclination in degrees.
Universe 09 00463 g002
Figure 3. Temporal behavior of K2-18 planets’ nodes and apsides in a two-planet configuration. The numerical integration of orbital motion was performed for 30 Kyr with the initial conditions: e 1 ( 0 ) = 0.01 , e 2 ( 0 ) = 0.05 , i 1 ( 0 ) = 0 , i 2 ( 0 ) = 5 , ω 1 ( 0 ) = ω 2 ( 0 ) = 0 , and Ω 1 ( 0 ) = Ω 2 ( 0 ) = 0 . Left panel: outer planet K2-18 b periastron argument after subtraction of a constant-rate precession of 0.0277 /yr. Right panel: outer planet K2-18 b longitude of the ascending node. Both angles are in radians.
Figure 3. Temporal behavior of K2-18 planets’ nodes and apsides in a two-planet configuration. The numerical integration of orbital motion was performed for 30 Kyr with the initial conditions: e 1 ( 0 ) = 0.01 , e 2 ( 0 ) = 0.05 , i 1 ( 0 ) = 0 , i 2 ( 0 ) = 5 , ω 1 ( 0 ) = ω 2 ( 0 ) = 0 , and Ω 1 ( 0 ) = Ω 2 ( 0 ) = 0 . Left panel: outer planet K2-18 b periastron argument after subtraction of a constant-rate precession of 0.0277 /yr. Right panel: outer planet K2-18 b longitude of the ascending node. Both angles are in radians.
Universe 09 00463 g003
Figure 4. The paths of the orbital axes of K2-18 planets over 30 Kyr in the inertial { x , y } plane. The red curve is for the inner planet 1 (K2-18 c), and the black curve is for the outer planet 2 (K2-18 b). The same initial conditions are used for numerical integration as in Figure 1.
Figure 4. The paths of the orbital axes of K2-18 planets over 30 Kyr in the inertial { x , y } plane. The red curve is for the inner planet 1 (K2-18 c), and the black curve is for the outer planet 2 (K2-18 b). The same initial conditions are used for numerical integration as in Figure 1.
Universe 09 00463 g004
Figure 5. Transit time variations in days for planet K2-18 b computed from the same simulation of a two-planet system as shown in Figure 1.
Figure 5. Transit time variations in days for planet K2-18 b computed from the same simulation of a two-planet system as shown in Figure 1.
Universe 09 00463 g005
Figure 6. Transit time variations in minutes for planet K2-18 b computed from a special simulation of a two-planet system with the same initial parameters as in Figure 1 for 30 Kyr with a state dump step of 1 yr. A small section of the data spanning 30 yr close to a minimum of the secular variation curve is shown.
Figure 6. Transit time variations in minutes for planet K2-18 b computed from a special simulation of a two-planet system with the same initial parameters as in Figure 1 for 30 Kyr with a state dump step of 1 yr. A small section of the data spanning 30 yr close to a minimum of the secular variation curve is shown.
Universe 09 00463 g006
Table 1. Assumed and adopted physical parameters of a two-planet K2-18 system.
Table 1. Assumed and adopted physical parameters of a two-planet K2-18 system.
ObjectMassPeriod, dRadius
K2-180.36 M s u n 0.46 R s u n
K2-18 b0.0272 M j u p 32.940.233 R j u p
K2-18 c0.0236 M j u p 8.9620.22 R j u p
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Makarov, V.V.; Goldin, A. Secular Orbital Dynamics of the Possibly Habitable Planet K2-18 b with and without the Proposed Inner Companion. Universe 2023, 9, 463. https://doi.org/10.3390/universe9110463

AMA Style

Makarov VV, Goldin A. Secular Orbital Dynamics of the Possibly Habitable Planet K2-18 b with and without the Proposed Inner Companion. Universe. 2023; 9(11):463. https://doi.org/10.3390/universe9110463

Chicago/Turabian Style

Makarov, Valeri V., and Alexey Goldin. 2023. "Secular Orbital Dynamics of the Possibly Habitable Planet K2-18 b with and without the Proposed Inner Companion" Universe 9, no. 11: 463. https://doi.org/10.3390/universe9110463

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop