© The Authors 2025

1 Introduction

Exoplanet research has rapidly expanded our understanding of planetary systems beyond our own (Winn & Fabrycky 2015). Astronomers have discovered over 5800 exoplanets so far, with nearly half occupying extremely close orbits around their host stars, completing their orbital revolutions in less than ten days (NASA Exoplanet Archive, consulted in April 2025). This population of close-in planets demonstrates remarkable diversity, ranging from small rocky worlds to giant gas planets. However, a notable and intriguing feature emerges within this population: a scarcity of Neptune-sized planets with orbital periods shorter than approximately three days. This phenomenon is known as the “hot Neptune desert” (e.g., Lecavelier des Etangs 2007; Szabó & Kiss 2011; Beaugé & Nesvorny 2012; Lundkvist et al. 2016).

Despite extensive research over the past decade, the fundamental origins and precise characteristics defining the boundaries of this desert are not entirely understood. Two primary hypotheses have emerged to explain this planetary population distribution. First, atmospheric evaporation driven by intense stellar radiation could play a crucial role (e.g., Owen & Lai 2018), potentially transforming Neptune-like planets into smaller, rocky bodies (Lecavelier des Etangs et al. 2004; Owen & Wu 2013; Lopez & Fortney 2013). Stellar X-ray and extreme ultraviolet (XUV) radiation can induce hydrodynamic atmospheric escape, particularly affecting planets in close proximity to their host stars (e.g., Lammer et al. 2003; Murray-Clay et al. 2009; Tripathi et al. 2015; Owen 2019; Vissapragada et al. 2022). Second, alternative theories propose that orbital migration mechanisms, including disk-driven migration (e.g., Goldreich & Tremaine 1979; Lin et al. 1996; Baruteau et al. 2016), secular chaos (e.g., Wu & Lithwick 2011; Hamers et al. 2017), planet-planet scattering (e.g., Ford & Rasio 2008; Nagasawa et al. 2008), and von Zeipel-Lidov-Kozai (ZLK, von Zeipel 1910; Lidov 1962; Kozai 1962) migration, might significantly influence planetary system architectures and composition.

To better understand the interaction between orbital migration and atmospheric escape for close-in exoplanets, there is a need for dedicated models to trace these processes over time, grounded in today’s observational constraints. While N-body integrators are highly accurate in simulating dynamical evolution (e.g., Rein & Liu 2012; Bolmont et al. 2020; Trani & Spera 2023), they often neglect planetary structure and fail to monitor long-term changes, limiting their ability to fully capture the history of older systems. On the other hand, atmospheric models provide detailed snapshots of planetary envelopes (e.g., Heng et al. 2012; Komacek & Youdin 2017; Johnstone et al. 2018) but struggle to account for temporal evolution and migration. The Joining Atmosphere and Dynamics for Exoplanets (JADE) code (Attia et al. 2021) was developed to address these limitations by integrating both orbital and atmospheric evolution, particularly focusing on late-stage migration scenarios in hot and warm Neptunes influenced by ZLK cycles. JADE is designed to simulate the intricate interplay between atmospheric erosion, orbital changes, and planetary structure in real time, offering a more comprehensive understanding of these phenomena in giant exoplanets.

In the first article (Attia et al. 2021), we benchmarked the JADE code on the intriguing case of GJ 436 b (Butler et al. 2004), a Neptune-sized planet still located at the fringes of the desert despite its current evaporation at tremendous rates (Bourrier et al. 2016). This planet further stands out because of its eccentric orbit (Trifonov et al. 2018), which should have been circularized a long time ago by tidal forces raised by the close star (as later discussed in Sect. 3.1). A possible scenario adduces a hidden distant companion generating a ZLK resonance resulting in a late migration (Bourrier et al. 2018; Attia et al. 2021). By forming farther out than its current position, GJ 436 b would have been trapped in the resonance for billions of years and migrated to a close-in orbit only recently; hence, it would have escaped the intense irradiation of the star when it was young and energetic (e.g., Jackson et al. 2012), thereby explaining its survival inside the desert. This narrative also naturally explains the residual eccentricity, as well as its polar spin-orbit angle (Bourrier et al. 2022), given the pervasive tendency of high-eccentricity migration mechanisms (HEM, which encompass ZLK) to elongate orbits and excite inclinations (e.g., Fabrycky & Winn 2009; Nagasawa & Ida 2011).

In the process, we revealed, for the first time in the context of ZLK cycles, a strong mutual feedback between the orbital and atmospheric history. Specifically, during the resonance, nearunity eccentricities periodically bring the planet in proximity to the star, heating and inflating its atmosphere to a significant extent, thereby inducing radius pulsations that are in tune with the orbital oscillations. While previous studies have explored coupled dynamics-atmosphere evolution involving radius inflation (e.g., Gu et al. 2003; Miller et al. 2009; Lopez & Fortney 2016; Glanz et al. 2022), the specific phenomenon of radius pulsations driven by the extreme eccentricity variations during ZLK cycles arising from the secularization of irradiation power that increases dramatically as e → 1 (Attia et al. 2021) has not been previously investigated. This is essential to account for in any secular investigation because of the substantial long-term effects of a varying radius in controlling the atmospheric structure as well as the intensity of tidal effects.

In light of these promising results, in this second article, we have resumed our the investigation of GJ 436 b, whose history might be prototypical of the evolution of intermediate-size planets at the periphery of the desert. Since the seminal definition of its borders by Mazeh et al. (2016), its inner edges have become less and less desertic, which led Castro-González et al. (2024) to update them, unveiling a surprising overdensity of warm Neptunes at the desert boundary (the so-called Neptunian “ridge”). As in the case of GJ 436 b, it is noteworthy that most planets populating this ridge possess eccentric orbits (Correia et al. 2020), coupled with substantial spin-orbit tilts (Bourrier et al. 2023; Attia et al. 2023). If all exo-Neptunes migrated early within their protoplanetary disk, we would expect the ridge to be dominated by circular and aligned systems. This is because tidal circularization and realignment are more effective at shorter periods, extinguishing any moderate primordial ellipticity and misalignment caused by disk-driven migration (e.g., Marzari & Nelson 2009). The ridge could consequently be the stable endpoint of the tumultuous HEM transporting Neptunian worlds from their supposed birthplace beyond the ice line, thereby solving the puzzle of the survival of their eroding atmospheres.

Assessing this scenario requires deploying advanced fully joint models of dynamical-atmospheric evolution for a large selection of Neptune-sized planets, which is made possible by the JADE code. While the first article (Attia et al. 2021) demonstrated the feasibility of the HEM-protracted erosion scenario for GJ 436 b, this study goes beyond by quantitatively constraining the system’s evolutionary pathway through a comprehensive parameter space exploration. Specifically, we seek to determine: (1) the planet’s formation location and initial mass; (2) the properties of the hypothetical companion GJ 436 c, providing observationally testable predictions for its mass and orbital period; and (3) the range of initial mutual inclinations compatible with the observed polar orbit. If a solution that is compatible with today’s observations would emerge, our narrative would be validated as viable, perhaps even hinting at similar histories for other ridge tenants. In Sect. 2, we recapitulate the building blocks of the JADE code and present its new features. Section 3 explains the preprocessing steps facilitating the application of JADE to GJ 436 b for its parameter space exploration, with the latter described in Sect. 4. Our results are then used to pinpoint the possible location of GJ 436 c within our ZLK framework in Sect. 5, while accounting for the existing observational constraints. We present our conclusions in Sect. 6.

2 The JADE code

2.1 Overview

The JADE code (Attia et al. 2021) simulates the evolution of a planet over secular timescales of several billion years (Gyrs) under the combined influence of realistic dynamical and atmospheric mechanisms. The dynamical evolution is governed by the action of an outer, distant third body in the system, capable of driving a ZLK orbital resonance of the close-in planet (e.g., Naoz 2016), counterbalanced by short-range forces (SRFs, encompassing tidal effects, spin distortion, and general relativity) binding the inner binary and quenching the resonance (e.g., Liu et al. 2015). The planet is modeled as composed of an iron core, a silicate mantle, and a H/He gaseous envelope, whose varying properties are self-consistently derived at each time step as the orbit, stellar irradiation, and inner planetary heating all evolve. The atmosphere easily erodes due to XUV-induced photoevaporation (e.g., Owen 2019) and internal released luminosity (e.g., Ginzburg et al. 2016), and JADE calculates this mass-loss rate using analytical formulae adjusted to detailed simulations of upper atmospheric structures (Salz et al. 2016b). This provides a better agreement with observed mass-loss rates than the commonly used energy-limited approximation (e.g., Watson et al. 1981; Lammer et al. 2003; Erkaev et al. 2007). Appendix A illustrates the overall functioning of the JADE code.

Fig. 1 Rosseland mean opacity, κR, as a function of the temperature, T, interpolated from the tables of Ferguson et al. (2005). Left: varying the metallicity, Z, at a fixed hydrogen mass ratio X = 0.8. Right: varying the hydrogen mass ratio, X, at a fixed metallicity Z = 0.0001.

2.2 Updates and new features

Since the initial release of the JADE code (Attia et al. 2021), we have implemented several new features for a more comprehensive and realistic treatment of close-in gas giants. While some of these improvements have been utilized in intermediate studies, we present here for the first time a complete description of all major updates, establishing this work as the authoritative reference for the current version of JADE. Accounting for the stellar Roche limit was the first of them, and succeeding JADE analyses insisted on its importance especially given the extreme adjacency of our targets to their hosts (e.g., Pezzotti et al. 2021; Almenara et al. 2022). Formally speaking, any simulation halts if the planet overflows inside the stellar Roche limit; in other words, if the following condition is reached $$a(1-e)<2\times R_{\mathrm{p}}{\left(\frac{M_{\star }}{M_{\mathrm{p}}}\right)}^{1/3}.$$(1)

The left-hand side corresponds to the distance at periastron, enforced by the (1 - e) factor, as tidal forces are maximal at that point of the orbit. The factor of 2 on the right-hand side is a conservative value to account for the fluidity of the tidally disruptable material (e.g., Chandrasekhar 1963; Teyssandier et al. 2019).

Moreover, a more rigorous calculation of the outer boundary pressure (i.e., the pressure at the photosphere or transit radius, where the integration starts) replaces the arbitrary former value of P_tr = 1 mbar. The procedure is the following: choosing a certain P_tr will set the gas density at the photosphere via the equations of state (EoS), which, in turn, will set an outer boundary opacity, κ_tr, through the collection of tabulated opacities JADE employs (Ferguson et al. 2005). Plus, the radiative conditions at the transit radius justify that P_tr is none other than the photospheric pressure¹ $$P_{\mathrm{t}\mathrm{r}}=\frac{\tau g}{{\kappa }_{\mathrm{t}\mathrm{r}}},$$(2)

where g is the surface gravity of the planet and τ is the optical depth (equal to 2/3 at the photosphere, following the Eddington approximation). Hence, P_tr is hence iterated over until convergence of Eq. (2), which is carried out every time the structure integrator is called since the outer boundary temperature dynamically evolves.

The layered prescription for the planet structure has also been improved. While the atmosphere remains unchanged, divided into an upper region collecting the high-energy incident flux and a lower region redistributing energy via radiation and convection, the solid material is now more faithful to the reality of planet interiors. It can be broken down into an iron core (α-Fe, ferrite) topped by a silicate mantle (MgSiO₃, perovskite), which is well-described by a temperature-independent modified polytropic EoS (Seager et al. 2007). Such a composition is reminiscent of Earth and provides a good trade-off between a simple-fast model and an honest agreement with complex interior models (e.g., Dorn & Lichtenberg 2021).

The final major update concerns the opacity treatment in the atmosphere. The JADE code employs the preferred formalism for such one-dimensional (1D) models: the Rosseland mean opacity (Rosseland 1924), which represents a weighted average across many monochromatic opacities. It simplifies the process by enabling models to perform a single calculation for all frequencies instead of adjusting the calculation for each one. Considering the complex impact of a gas on a spectrum, this approach makes the models far more tractable. Atmospheres simulated by the JADE code are composed of H/He, and the novelty now is that a trace amount of metals can also be introduced in the gaseous material, which leads to a massive change in the associated opacities. This is best visible in Fig. 1, where the Rosseland mean opacity κ_R is represented as a function of the gas temperature T for a variety of compositions, as interpolated from the tabular entries used by JADE (Ferguson et al. 2005). The left panel shows the influence of the atmospheric metal mass fraction Z, especially explicit in the planetary regime (lowest temperatures). Adding even a hint of heavy elements pumps the opacity from a negligible value to a physical one. If the presence of metals does not seem to be important for the hottest atmospheres, low temperatures on the other hand inhibit collision-induced absorption, making heavy dust grains absorbing light the main source of opacity. This is corroborated by the linear increase in κ_R with Z, seen for T < 1000 K. On the contrary, the right panel tells us the minimal sensitivity of κ_R to the hydrogen fraction X, all other things being equal, except for the highest temperatures, even though marginally. Indeed, atomic absorption becomes the biggest contributor to the opacity for the hottest gases, while a lower X means more helium Y, which translates into helium lines quickly dominating over hydrogen lines, leading to a lower κ_R for T > 15 000 K (e.g., Rogers & Nayfonov 2002; Mendoza et al. 2007).

Our atmospheric integrator as initially designed did not include metals, which resulted quite systematically in a nonconvergence for low-irradiated planets: the lack of dust grains made the atmosphere virtually transparent, and an intense, nonphysical surface gravity had to be set to counterbalance, via hydrostatic equilibrium, the very efficient internal radiation. Introducing trace quantities of metals (e.g., Z = 0.0001) has the benefit of restoring the atmospheric opacity without having to implement additional EoS, since the mean molecular weight stays almost unaltered. It should additionally be noted that the exact mixture of elements composing Z is largely unknown in real-life cases, so in the JADE code it is set to match solar abundances (according to Asplund et al. 2021) by default. Even though sometimes host stars’ compositions might be used as a proxy for the contents of their affiliated planets, such procedures remain hazardous because the correlation is not robust enough, especially for the smallest stars (Hinkel et al. 2024, and references therein).

3 The flagship case of GJ 436: Setting the scene

3.1 Motivation

The aim of our previous work (Attia et al. 2021) was to exhibit the potential of a HEM-delayed erosion scenario, enabled by our novel fully coupled approach. We sought to explore as many aspects as possible and explain the peculiarities of this system altogether. Here, the motivation is more quantitative: an exploration of the parameter space meant to predict the possible properties of the concealed companion, which would be responsible for a ZLK resonance that leads to a configuration compatible with the observed reality. In doing so, most of the use cases of Fig. A.2 are exploited, ultimately culminating in RE1 and RE2. Before diving into such an analysis, it is important to take a step back to examine our target and identify the questions to be answered. We refer to Attia et al. (2021) for a more comprehensive description of the system and we recapitulate here the main elements that make this tenant of the hot Neptune desert of particular interest in this paper.

An eccentric orbit. GJ 436 b has been the subject of a growing interest in the past couple of decades, making it the target of several high-precision radial-velocity (RV) campaigns. One of the undisputed results of such Doppler measurements is the elliptic nature of its orbit (eccentricity of e = 0.152 ± 0.009, Trifonov et al. 2018), which has been progressively confirmed with higher and higher confidence through the years (Butler et al. 2004, 2006; Maness et al. 2007; Knutson et al. 2014; Lanotte et al. 2014; Turner et al. 2016; Rosenthal et al. 2021). The close proximity to the star (semimajor axis of a = 0.02858 ± 0.00043 AU, Maxted et al. 2022) should have led tidal friction to circularize the orbit a long time ago, given the advanced age of the system (4-8 Gyr, Bourrier et al. 2018). A straightforward solution could be the possible weakness of tidal forces (Trilling 2000; Mardling 2008). However, this argument is solely grounded in our poor knowledge of the tidal dissipation factor Q_p, which is required to be abnormally high for this hypothesis to work. Alternatively, Beust et al. (2012) formulated a more convincing explanation based on ZLK cycles generated by a distant, unseen companion and ending in a late inward migration (inspired by Tong & Zhou 2009; Batygin et al. 2009), which would also have the merit of solving the other unanswered interrogations presented below.

A polar architecture. Motivated by the appealing HEM solution to the eccentricity conundrum, Bourrier et al. (2018) measured the spin-orbit angle of GJ 436 b to investigate if it validates this scenario. The result is crystal clear: the orbit is not only misaligned, but joins the striking sample of polar ones (Albrecht et al. 2021; Attia et al. 2023; Siegel et al. 2023), lending even more credence to HEM. Such a finding was later confirmed with more precise RV data and the advent of the advanced Rossiter-McLaughlin revolutions (Bourrier et al. 2021) analysis technique (3D spin-orbit angle of $\psi ={103}_{-12}^{+13}{}^{\circ }$, Bourrier et al. 2022).

An evaporating atmosphere. Thanks to spectroscopic transit observations in Lyman-α, giant clouds of escaping hydrogen have been revealed enshrouding GJ 436 b and trailing the warm Neptune along its orbit (Kulow et al. 2014; Ehrenreich et al. 2015; Lavie et al. 2017; dos Santos et al. 2019). Subsequent models took profit of these results and unraveled many important characteristics of the eroding atmosphere, such as its comet-like structure (e.g., Shaikhislamov et al. 2018; Khodachenko et al. 2019; Villarreal D’Angelo et al. 2021), the space weather around it (Vidotto & Bourrier 2017; Bellotti et al. 2023; Vidotto et al. 2023), and, crucially its mass-loss rate (~2.5 × 10⁸ g∕s, Salz et al. 2016a; Bourrier et al. 2016). While the latter value would require ~200 Gyr to erode just 1% of the planet’s total mass, the situation was vastly different during the star’s youth. During the stellar saturation phase in the first few hundred Myr, XUV luminosities are typically 100-1000× higher than present values (Jackson et al. 2012; Pezzotti et al. 2021). Since atmospheric escape rates scale approximately linearly with XUV flux (as seen in the energy-limited regime, e.g., Lammer et al. 2003), this implies mass-loss timescales on the order of 200 Myr for each percent of total mass during the saturation phase. Over the system’s lifetime (>4 Gyr), this would result in several percent of cumulative mass loss: sufficient to substantially alter the planet’s radius given that H/He atmospheres contribute disproportionately to planetary radii (e.g., Zeng et al. 2019). Such erosion could have pushed GJ 436 b well outside the hot Neptune desert boundaries. The HEM scenario thus reconciles the atmospheric escape we are witnessing today with the planet’s survival at the inner fringes of the desert, explaining how it retains a substantial volatile envelope despite its advanced age.

3.2 Bulk composition

A few key steps are required prior to any exploration of the parameter space (Appendix A). Determining GJ 436 b’s bulk structure and composition (RS1, RS2, Fig. A.2) needs to be carried out first, since it conditions the outcome of any other subsequent analysis. The preferred approach would be a Bayesian inference of all the relevant parameters, constrained by today’s observed properties (e.g., Acuña et al. 2021). Nevertheless, with only the planet mass and radius as observational inputs, the result is highly degenerate and no clear solution emerges. Instead, we chose to employ the following “standard” simplifying assumptions, which are not expected to heavily influence the secular evolution:

• The mantle-to-core mass ratio is set to 2:1, similar to that of Earth;

• The atmospheric helium mass fraction is set to Y = 0.2, much like Neptune (Hubbard et al. 1995; Helled et al. 2020);

• The atmospheric metal mass fraction is set to Z = 10^-5, a trace amount that still provides a dominant source of opacity at low temperatures (Fig. 1) without affecting the mean molecular weight (see discussion in Sect. 2.2).

In this way, only one degree of freedom remains when it comes to the planet structure: M_core, the core mass². This parameter is derived by running a large set of internal structures under the above assumptions, with a fixed total mass of M_p = 21.68 M_⊕ (Maxted et al. 2022) and the present-day orbital properties for all of them. Each simulation differs from the others by its core mass, and the model that is ultimately retained is the one that yields the observational radius, R_p = 3.943 ± 0.076 R_⊕ (Maxted et al. 2022). All the planet structures were simulated at t = 6 Gyr, median of the system age range (Bourrier et al. 2018), as the time is critical in terms of controlling the stellar luminosity (computed using GENEC, a stellar model adjusted to our system, Eggenberger et al. 2008) and planetary internal heating; and, thus, in influencing the global radiative budget. Accordingly, we found that 89% of today’s planet mass is enclosed in the solid material, which corresponds to M_core = 19.30 M_⊕. The core mass remains unchanged in any evolutionary simulation: we assume for simplicity that the iron and silicates do not dissolve in the volatile envelope; although it could be a possibility for such giant planets (Guillot et al. 2004; Wilson & Militzer 2012). Consequently, this invariable value of M_core is the one employed in all simulations of the exploration. The aforementioned assumptions are also imposed in the exploration, for a clear consistency.

3.3 Grid of internal structures

We note that a full exploration including hundreds of thousands of evolutionary models would be made unfeasible by the prohibitive computational cost of the structure integrator. Even for a single simulation, as the orbit evolves over billions of years, as the stellar XUV luminosity gradually declines, as the internal heating progressively fades out, and as the atmosphere evaporates, the planet’s internal structure (and, ultimately, the planet radius) need to be regularly adjusted to follow these changes, which intensively calls the structure integrator. To overcome this limitation, a grid of internal structure models was precomputed, covering a maximum number of configurations that could arise during the exploration (FA2, Fig. A.2). Once the “static” bulk parameters (RS1, RS2, Fig. A.2) have been fixed, only three physical properties influencing the planet radius remain variable: the planet mass, M_p, the equilibrium temperature, T_eq (related to stellar irradiation), and the intrinsic temperature, T_int (related to planetary internal luminosity, with both temperatures being rigorously defined in Attia et al. 2021). The constructed grid spans 30 M_p values, 42 T_eq values, and 28 T_int values, yielding a total of 35 280 distinct points. The chosen values for each of the three parameters result from a uniform sampling of their conservatively large intervals of definition (M_p ∈ [M_core; 80 M_⊕], T_eq ∈ [150K; 1500K], and T_int ∈ [10K; 450K]) later refined where convergence is difficult (i.e., M_p > M_core, and close to the definition boundaries for T_eq and T_int).

The key parameter derived from every point of the grid is the planet radius, which is essential for the proper calculation of the dynamical forces (tides in particular) and the mass-loss rate. We could simply stop at the creation of this grid to complete FA2, by exporting the computed radii and interpolating them as a function of M_p, T_eq, and T_int within the evolutionary simulations that are part of the exploration, which would still drastically decrease the run time in theory. However, a few points in the grid actually fail to converge, especially in the regions with a refined resolution, breaking down the grid-like structure. If interpolation is relatively straightforward on a regular mesh (via, e.g., trilinear interpolation), schemes defined for scattered points are much more complex, in particular, for a dimensionality higher or equal than three. In our case, numerical methods are not guaranteed to converge and would sometimes negate the expected run time gain in practice.

3.4 Mass-radius relationships

Instead, we derive analytical relationships for R_p as a function of the other parameters. Evaluating a scalar function (instead of interpolating) is the optimal solution from a numerical standpoint, but it requires “guessing” the shape the formula would take. We could, for example, try to derive it from first principles (as ingenuously executed by Turbet et al. 2020) and then compare its results to the grid radii. A simpler and more flexible solution is to utilize a generic parametric form, subsequently adjusted to fit the data (much like the vast majority of mass-radius relationships in the literature, e.g., Grasset et al. 2009; Swift et al. 2012; Zeng et al. 2016; Otegi et al. 2020; Parc et al. 2024), which is the approach chosen here. Attempts at constructing a radius function with an explicit dependency to all three parameters M_p, T_eq, and T_int unfortunately gave inconclusive results. Instead, we resorted to classical mass-radius relationships, whose coefficients depend on the two temperatures. The simplest form we might conceive of would be a power law $R_{\mathrm{p}}\propto M_{\mathrm{p}}^a$, where a is the power-law exponent. An observational fit to the telluric planets of the Solar System yields a = 0.33 (e.g., Wolfgang & Laughlin 2012) like for constant density spheres, whereas including gas giants gives higher values (e.g., Lissauer et al. 2011), but the most massive super-Jupiters are expected to have a negative exponent (Chabrier et al. 2009). To account for these different regimes, the best parametric form modulating the power law, which we found by trial and error, being initially inspired by Mordasini et al. (2012); Aguichine et al. (2021), is the following $$\log R_{\mathrm{p}}=a{(\log \frac{M_{\mathrm{p}}}{b})}^c+d,$$(3)

where a, b, c, and d are coefficients to be fitted. In this regard, they are determined for each isotherm (T_eq, T_int) in the grid by a least-squares minimization, so that Eq. (3) delivers the closest radius possible to the one computed by the JADE code³. Importantly, scrutinizing the results in the form of a pair-by-pair corner plot (Fig. 2) reveals that some coefficients are in fact correlated. This is a considerable advantage since it allows us to reduce the number of free parameters: c and d can be expressed as a function of a, a power law for the former and a simple linear relation for the latter, as explicitly seen in Fig. 2. Therefore, we can write $$\log c=c_0\log a+c_1,$$(4) $$d=d_{0}a+d_1,$$(5)

where c₀, c₁, d₀, and d₁ are, this time, unique coefficients (constant for any combination of T_eq, T_int ). Their exact values, also found with a least-squares minimization, are shown in Table 1, but strongly suggest that c ∝ 1/a and d ≃ 2 - a. Their extreme proximity to integer values might be indicative of a physical reason, but we did not find any convincing lead. As a matter of fact, fixing them to their integer counterparts degrades the goodness of the fit, so we decided to use the values of Table 1 as they are.

To summarize, the procedure to calculate the planet radius at any moment of any evolutionary simulation can be sketched this way:

1. Computing T_eq and T_int (Attia et al. 2021);

2. Bilinearly interpolating a and b from T_eq and T_int;

3. Computing c and d from a using Eqs. (4) and (5);

4. Injecting M_p, a, b, c, and d in Eq. (3).

Compared to one run of the structure integrator, determining the radius this way represents an immense run time gain, as illustrated in Appendix B. Figure 3 depicts some massradius relationships for a couple of combinations of T_eq and T_int. On each isotherm are overplotted a few radii directly computed by JADE’s structure integrator, showcasing the agreement with the analytical formulae, as also further explored in Appendix B.

Fig. 2 Corner plot for the fitted coefficients a, b, c, and d from Eq. (3). Each black point is the result of fitting one internal structure model from the precomputed grid (Sect. 3.3).

4 The flagship case of GJ 436: Exploration

4.1 Framework

The aim of this investigation is to constrain at the same time the initial conditions, after the disk dissipation, of planet b (RE1, Fig. A.2), as well as the properties of the hypothetical planet c (RE2, Fig. A.2). To this effect, we employed a novel framework that combines grid-based forward modeling with Bayesian inference. We termed this approach “semi-Bayesian” because even though we did use Bayesian methods (priors, likelihoods, and MCMC sampling) to explore the parameter space, the likelihood computation relies on precomputed grids of simulations rather than on-the-fly calculations. Specifically, we first generated a comprehensive grid of fully coupled simulations carried out over the upper age limit of the system, 8 Gyr, spanning a wide domain of the parameters to be constrained. The likelihood at any point in parameter space was then obtained by interpolating between these precomputed grid points.

When it comes to RE1, the two main properties that we explore, characterizing the infancy of GJ 436 b, are the initial semimajor axis, a₀, and the initial planet mass, M_p,0. The initial eccentricity, argument of periastron, and spin-orbit angle are left unexplored because, provided that the ZLK mechanism occurs, they are bound to oscillate. Their initial values, unless they are extreme, will be regularly reached during the long-term secular evolution of the system, because the ZLK characteristic timescale, τ_ZLK, imposing the period of the cycles within a factor on the order of unity (Beust & Dutrey 2006), is much smaller than the global timescale of the evolution. In other words, two systems with different initial eccentricities (or arguments of periastron and spin-orbit angles), all else being equal, will differ by a phase shift. While the previous statement is technically not entirely true, as very specific and well-chosen initial conditions can for instance restrict the amplitude of the cycles, we performed a few such simulations to corroborate its conclusion. In this way, all inner orbits were initialized as circular and aligned with the spin of the star. The planet axial tilt, monitored by the JADE code, also starts as null.

As for RE2, the perturber would be fully defined, ZLK-wise, via its mass, M_pert, semimajor axis, a_pert, eccentricity, e_pert, and initial mutual inclination, i_mut,0, with respect to the inner orbit. We assume for simplicity to truncate the Hamiltonian up to the quadrupolar order, which enables very fast computations without an excessive accuracy loss. Beust et al. (2012) carried out dynamical studies at various orders, showing that the two-fold ZLK evolution reported below (Fig. 4) still holds whatever the truncation order. Only the transition time, t_trans, appears sensitive to the truncation order (and defined in Sect. 4.2), with an overestimation of ~30% at the quadrupolar level (Fig. 4 of Beust et al. 2012). Within this hypothesis, we deduce from the equations of motion (Appendix A.3 of Attia et al. 2021) that the mass, semimajor axis, and eccentricity of the companion can in fact be combined into a single parameter controlling the intensity of the resonance, $$\mathrm{\Lambda }\equiv \frac{a_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}}^3{(1-e_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}}^2)}^{3/2}}{\mathcal{G}{M}_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}}},$$(6)

which is the parameter we explore in the end, allowing for a very welcome reduction in dimensionality. In all occurrences where numerical values of Λ are given, a_pert is assumed to be in AU, M_pert in M_⊕, and the gravitational constant, G = 1. Of course, we numerically checked that two perturbers with different properties leading to the same value of Λ actually generated an identical ZLK resonance. It should be noted that Λ is proportional to a time squared, and is highly reminiscent of the ZLK characteristic timescale. This is no coincidence, and the two quantities can in fact be linked through ${\tau }_{\mathrm{Z}\mathrm{L}\mathrm{K}}\simeq \mathrm{\Lambda }/P$, with P as the inner orbital period. All in all, three parameters needed to be constrained: a₀, M_p,0, and Λ. For this specific section and the following, the initial mutual inclination is fixed to i_mut,0 = 75° as a typical value so as to isolate its influence (Sect. 4.2), but it was also explored in a second time (Sect. 4.3).

A fundamental challenge arises when comparing our evolutionary models to observations regarding how we can evaluate the compatibility between simulations that predict time-varying planetary properties over 8 Gyr and observational data that represent a snapshot of the system at a single, uncertain epoch. The system’s age is only constrained to lie within 4-8 Gyr and properties such as eccentricity, mass, and spin-orbit angle all evolve throughout this window. Simply checking whether a simulation matches observations at a fixed time would be inadequate, as it would miss potentially valid evolutionary pathways that reproduce the observed state at different epochs within the age uncertainty. Therefore, we would need a statistical framework that accounts for both the temporal evolution of the simulated properties and the uncertainty involved when these properties would be match up with the observations.

Then, for each 8-Gyr fully joint simulation corresponding to a distinct point (a₀, M_p,0, Λ) in the explored 3D parameter space, a score function needs to be evaluated to quantify the compatibility with the currently observed parameters. To this effect, each simulation can be broken down into a temporal collection of evolving system properties (time, semimajor axis, eccentricity, mass, etc.). For every point i of the time series, a log-likelihood of the constraining data (today’s observed parameters) with respect to the model log p_i is computed using the distance between i and the posterior distribution underlying the constraining data. To make things clearer, we can take a concrete example, with only one constraining parameter: today’s eccentricity, with its error bar, e_• ± σ (e_•). Observational measurements are routinely assumed to take the shape of a Gaussian distribution, the log-distance between i and the constraining eccentricity would be - ((e_i - e_•)/σ (e_• ))²/2, where e_i is the eccentricity of point i of the simulation. The latter quantity needs to be complemented with a time information, since obtaining the correct eccentricity needs to happen within the correct time span as well. Here, the constraint we have on the system age is uniform between two bounds, t_•↓ and t_•↑, which gives the total log-likelihood for point i of the simulation constrained by the current eccentricity, $$\log p_i=-\frac{1}{2}{\left(\frac{e_i-e_{•}}{\sigma (e_{•})}\right)}^2+\log (1_{[t_{•\downarrow },t_{•\uparrow }]}(t_i)).$$(7)

This value can be obtained up to an additive constant, with t_i the time value of point i and 1 the indicator function that is equal to 1 when t_i falls within the allowed age range and 0 otherwise. This ensures that only configurations matching the observations within the system’s age constraints contribute to the total likelihood, with no preference for any particular age within that range. We can see from Eq. (7) that the i points maximizing their log p_i are the ones yielding the closest e_i to e, during the time interval [t_•↓, t_•↑], which is exactly the desired behavior. In practice, there are several constraining observational measurements; each contribution is added to log p_i, assuming it is independent from the others. Finally, the total log-likelihood of a (a₀, M_p,0, Λ) simulation is computed by numerically time-integrating the various p_i pertaining to the individual points of said simulation, $$p=\int p_i(t)\mathrm{d}t\simeq \sum _i{p}_i\mathrm{\Delta }{t}_i.$$(8)

From this, we can then calculate $$\log p=\log p_{\mathrm{m}\mathrm{a}\mathrm{x}}+\log \sum _i\exp (\log p_i-\log p_{\mathrm{m}\mathrm{a}\mathrm{x}})\mathrm{\Delta }{t}_i,$$(9)

with log p_max as the maximal value of the different log p_i, this formulation of Eq. (9) avoiding null values at machine precision in the evaluation of the exp function. It is important to note that this temporal integration approach inherently addresses the likelihood of observing the system in its current state. Simulations that match the observational constraints for longer durations within the allowed age range [t_•↓, t_•↑] naturally accumulate higher total log-likelihoods through the summation in Eq. (9). Conversely, configurations that only briefly pass through the observed parameter space, such as those undergoing rapid migration, receive proportionally lower scores. This weighting scheme thus automatically penalizes fine-tuned solutions that require precise timings to match current observations. We describe more rigorously the construction of the probability function in Appendix C.

After setting a uninformative prior for all three jump parameters, we were equipped with a grid of JADE simulations, paving the (α₀, M_p,0, Λ) parameter space, each one of them associated with a log p. Intermediate values, uncovered by JADE models, have their log p interpolated from the grid. Ultimately, the regions of high likelihood are sought for by employing an MCMC implementation using emcee (Foreman-Mackey et al. 2013). Because such regions can be highly multimodal, we found that using a custom “move” parametrization was best to efficiently explore the entirety of the parameter space (a move is the algorithm for proposing an update of the coordinates of the sampled posterior at each step). In our case, a well-behaved setup is a weighted combination of a differential evolution proposal (Nelson et al. 2013) with a probability of 90% and a snooker move (ter Braak & Vrugt 2008) with a probability of 10%, which has the crucial benefit of avoiding being stranded in a local maximum by occasionally allowing for large moves. In any case, the length of the MCMC chains (typically 10000), number of walkers (typically 64), and burn-in phase (typically 5000) are adjusted to ensure convergence.

Fig. 3 Analytical mass-radius relationships for GJ 436 b. The crosses represent some radii directly computed by the structure integrator to emphasize the goodness of the fit. Left: varying equilibrium temperature, Teq, with a fixed intrinsic temperature Tint = 100 K. Right: varying intrinsic temperature Tint, with a fixed equilibrium temperature Teq = 500 K.

Fig. 4 Secular evolution of the system corresponding to the maximum log-likelihood of the fiducial exploration (imut,0 = 75°, Sect. 4.2), as simulated by JADE. The orange shaded regions represent the observational constraints of Table 3. In all plots, ttrans is shown as a dashed vertical red line.

4.2 Results: Fiducial inclination

To maintain the readability of the results, we fixed the initial mutual inclination to a typical value i_mut,0 = 75° as a first step. In total, the base grid of simulations we carried out for this fiducial exploration comprises 113 986 points, each corresponding to a different (α₀, M_p,0, Λ) configuration. Table 2 shows the set of fixed parameters common to all simulations. When it comes to the companion properties, we imposed a circular orbit at 20 AU. This way, the explored values of Λ are unambiguously controlled by M_pert via Eq. (6), preserving in all cases the hierarchical triple assumption (a_pert ≫ a). This does not mean that the companion lies in fact at 20 AU, as a Λ solution is degenerate between an infinity of compatible combinations of a_pert, e_pert, and M_pert. We recall that the perturber remains unchanged during the entire simulation, acting like an angular momentum reservoir. Moreover, the inner planet’s spin rate is always initialized to be synchronized with the orbital motion, as its value is unknown. As a matter of fact, it is a safe assumption because the planet spin has a minimal impact on the secular evolution in our situation, although it could have a greater influence on energy dissipation in a more general case, as part of advanced tidal models (Makarov & Efroimsky 2013). Plus, Beust et al. (2012) showed that during the resonance phase, the planet’s spin tends to accelerate to synchronize with the orbital angular velocity at periastron in the high eccentricity phases of ZLK cycles, that is to say when SRFs are strongest. The choice of planet b’s initial argument of periastron is motivated by the fact that ω = 0° represents the locus of minimum eccentricity within a quadrupole approximation (Lithwick & Naoz 2011), consistent with the initial circular orbit.

Table 3 lists the various observational measurements constraining the evolutionary simulations and contributing to the calculation of log p (Eq. (9)). As seen in Fig. D.1, a clear solution emerges, embodied by the peaks of maximum likelihood in the 1D histograms. Configurations falling within the shaded regions, which can be interpreted as the uncertainty intervals for the jump parameters [a₀, M_p,0, Λ), yield outcomes satisfactorily close to the currently observed reality. This is best visible in Fig. 4, simulating the evolution of a system corresponding to the medians of the posterior distributions (their numerical values are reported on top of the figure, as well as above the 1D histograms in Fig. D.1). All the observational constraints of Table 3 are matched for this maximum-likelihood simulation, except for ψ. as discussed later. The typical two-fold ZLK evolution is explicit. In a first phase, dominated by the resonance, large eccentricity oscillations are seen, modulated by a shrinking envelope due to SRFs. The spin-orbit angle wildly oscillates between 0° and 150°, while the semimajor axis remains roughly constant. When the amplitude of the eccentricity cycles becomes negligible, GJ 436 b decouples from the companion: the resonance is over. This point in time t_trans is crucial as it marks the transition to the second phase, which is dominated by tidal damping, featuring an abrupt inward migration and locking the planet on a misaligned orbit.

The total impact of mass loss is extremely mild, to say the least, but the radius variations are worth commenting. The sharp decrease at t ~ 200 Myr coincides with the end of the stellar saturation phase, when the XUV flux suddenly drops by orders of magnitude, leading to rapid atmospheric contraction (similar behavior is shown in Figures 9-11 of Attia et al. 2021). The subsequent radius increase appears linked to T_int dropping below 100 K due to the sharp decrease in planetary luminosity at this epoch. While the exact mechanism causing this inflation at this specific temperature threshold remains unclear, it likely reflects a transition in the atmospheric opacity regime or convective properties that our mass-radius relationships capture empirically but whose physical origin merits further investigation. Then, the radius globally mimics the steady increase in the equilibrium temperature due to the slowly increasing mean eccentricity, before inflating at t_trans because of migration and the resulting additional heating. Finally, the atmosphere slowly contracts because of evaporation, despite it being minimal. The same radius pulsations in tune with the ZLK cycles, identified in Attia et al. (2021), are reported all along the first phase.

Globally speaking, the following conclusions can be drawn from a meticulous inspection of the corner plot (Fig. D.1) and the underlying simulations. First, configurations from inside the uncertainty margins of the posterior distributions are generally compatible with Table 3 within the error bars of the constraints (except perhaps for ψ_•, see discussion below), which is exactly the coveted behavior, lending credence to our semi-Bayesian model.

The observational constraint on the semimajor axis is the most restrictive especially for the inferred a₀, which can be intuited from a theoretical standpoint. The value of the final separation a_f of GJ 436 b is controlled by tidal damping, main dynamical engine after the decoupling from the companion. During this phase, a (1 - e²) is a roughly constant⁴ quantity thanks to the conservation of angular momentum (Socrates et al. 2012). Since the fate of the orbit is ultimately to circularize, the final separation is dictated by a_f = a(t = t_trans) (1 - e(t = t_trans)²). Nonetheless, a(t = t_trans) ≃ a₀ because tides, responsible for the orbit shrinkage, are largely surpassed by the resonance counteracting them during the first phase (as seen in Fig. 4). Furthermore, e(t = t_trans) can be seen as a proxy for the maximum eccentricity reached during the ZLK cycles e_max, which only depends on the initial inner eccentricity/argument of periastron and the initial mutual inclination to first order within the quadrupole approximation (Katz et al. 2011). All these parameters being fixed (Table 2), e_max should be roughly the same for all explored configurations (limited variations might subsist due to SRFs). Consequently, a_• (the desired value for a_f) constrains a₀ very strongly, in particular given the low uncertainty σ(a_•) due to the minute precision of CHEOPS (Benz et al. 2021). This argument is validated by the very peaked distribution of maximum eccentricities in the simulations where the resonance was triggered e_max = 0.941 ± 0.005, and the ensuing theoretical expectation $a_0=a_{•}/(1-e_{\mathrm{m}\mathrm{a}\mathrm{x}}^2)=0.249\pm 0.009\mathrm{A}\mathrm{U}$, which is satisfyingly consistent with the inferred value from the exploration (Fig. D.1).

The observational constraint on the eccentricity is actually not as restrictive as anticipated. It is understandable, a posteriori, because e_• represents but a certain nonzero value that is eventually browsed after the decoupling, if t_trans falls reasonably within [t_•↓; t_•↑]. This interval being quite large, the constraint brought by e_• on t_trans translates into a moderate constraint on Λ, which explains the variation of the latter of more than half a dex within 1 σ of its posterior.

Similarly to the case of the maximum-likelihood simulation (Fig. 4), evaporation has a minimal impact on compatible configurations in general. By forming farther out than its current orbit, GJ 436 b is spared from the infancy of the star, when it is most energetic (e.g., Jackson et al. 2012). When it migrates to its final position, it is already too late for evaporation to make a difference: the stellar XUV luminosity has already dropped by several orders of magnitude. In fact, the inferred M_{p, 0} along with its uncertainty are encompassed in the error bar of M_p•, consistent with a null net effect of mass loss. Nevertheless, accurately modeling the atmospheric structure is still pivotal because of the radius pulsations and their meaningful influence on tides ($\propto R_{\mathrm{p}}^5$, e.g., Jackson et al. 2008), to which t_trans is extremely sensitive (Beust et al. 2012; Attia et al. 2021).

The observational constraint on the spin-orbit angle is the most difficult to meet. Only a handful of simulated planets emerge from the resonance with a compatible ψ. We can check these results by examining the distribution of final spin-orbit angles, ψ_f, within simulations that have the expected two-fold behavior (Fig. D.2). The bulk of the orbits finish their runs on moderately tilted orbits between 15° and 60°, and a tiny fraction ends up on a retrograde orbit. The most striking feature is the total absence of intermediate values of ψ_f between these two modes, which is related to our choice of initial mutual inclination. This final distribution is similar to the one produced by Fabrycky & Tremaine (2007), with the major difference that their two modes spread enough to overlap, leaving no void between them. This is due to the fact that they also investigated a large number of i_mut,0. In any case, the role of ψ_• will become central when i_mut,0 is explored (Sect. 4.3).

The radius variations observed in our simulations warrant comparison with previous studies of coupled planetary evolution. While several works have investigated the effects of a dynamically evolving atmospheric structure, the specific radius pulsations we observe during ZLK cycles arise from a distinct mechanism. Gu et al. (2003) and Miller et al. (2009) focus on tidal heating’s impact on planetary radii, a process not implemented in JADE. However, since tidal heating scales with |ȧ/a| (neglecting spin variation, Eggleton et al. 1998; Mardling & Lin 2002) and the semimajor axis remains roughly constant during the ZLK phase (Fig. 4), this omission should minimally affect our radius pulsations, which instead originate from fluctuations in stellar irradiation. Lopez & Fortney (2016) include all relevant contributions to the planet energy budget except tides (similar to JADE, even though we do not implement them in a self-consistent way) but do not secularize the irradiation, focusing instead on radius evolution during the stellar postmain sequence. Our pulsations are fundamentally an eccentricity-driven phenomenon emerging from the secularization of irradiation power, which increases substantially as e approaches unity during ZLK cycles (Attia et al. 2021). Glanz et al. (2022) implement a comprehensive treatment including secularized irradiation and could in principle capture these pulsations, though their focus is on the tidal circularization phase. We acknowledge that during the rapid migration after t_trans, tidal heating becomes significant and our model surely underestimates the radius inflation during this phase.

We also note that the tidal dissipation constant Q_p is fixed to a fiducial value of 10⁵ for all simulations (Table 2). While a higher Q_p would assuredly delay the onset of tidal migration (shifting t_trans to later times), it would not extend the duration over which the simulated eccentricity matches the observed value. The eccentricity evolution would simply cross the observationally compatible region at a later time, maintaining a similar crossing duration. Crucially, this means that the total log-likelihood log p (Eq. (9)) would not be enhanced by adopting a higher Q_p; in other words, a slower migration has no intrinsic advantage in our framework as long as the crossing occurs within the allowed age range. Given the poor observational constraints on Q_p for Neptune-mass planets, we chose not to vary this parameter, focusing instead on the initial conditions that more directly control the evolutionary pathway.

4.3 Impact of the initial mutual inclination

Showing the results for a fiducial initial mutual inclination in this way (Sect. 4.2) allows us to scrutinize the influence of the other parameters in isolation and deepen our understanding of the ZLK mechanism. We now expand our analysis beyond this case to assess how the initial mutual inclination impacts our exploration. For each explored value of i_mut,0, we conducted the same comprehensive semi-Bayesian retrieval process described previously (Sect. 4.1), with each refined grid comprising approximately 100 000 simulations. All other simulation parameters (Table 2) and observational constraints (Table 3) remain unchanged.

The results of our (a₀, M_p,0, Λ) retrieval as a function of i_mut,0 are shown in Table 4. Globally, we observe the same behaviors and reach similar conclusions as in Sect. 4.2, particularly regarding the minimal impact of atmospheric evaporation, corroborated by the consistent derived values of M_p,0 across all explored initial mutual inclinations. Notable differences emerged when comparing higher initial mutual inclinations (80° and 85°) to the fiducial case. For these higher values, the maximum eccentricity reached during ZLK cycles e_max increases (Katz et al. 2011), necessitating larger initial separations a₀ due to angular momentum conservation as discussed in Sect. 4.2. This in turn leads to higher derived values of Λ to maintain a roughly constant ZLK timescale ${\tau }_{\mathrm{Z}\mathrm{L}\mathrm{K}}\propto \mathrm{\Lambda }/a^{3/2}$, as illustrated in Table 4. Such a behavior does not hold anymore for the two lowest inclinations (60° and 70°), which is due to the much more chaotic nature of their evolutions, where SRFs become comparable in magnitude to the ZLK perturbation. The clear two-fold behavior showcased by the three highest i_mut,0 (e.g., Fig. 4) is not seen anymore in the low-inclination regime. In any case, the retrieval seems to favor configurations with much shorter ZLK timescales for them, yet with quite large uncertainty intervals.

The impact of initial mutual inclination is particularly striking when examining the distribution of final spin-orbit angles ψ_f after decoupling (Fig. 5). As i_mut,0 increases and approaches 90°, the associated ψ_f distribution becomes wider, with a growing proportion of retrograde orbits, akin to the results of Fabrycky & Tremaine (2007). For the lowest investigated initial mutual inclinations (60° and 70°), no retrograde orbits emerged after decoupling from the perturber, which is why these distributions have been excluded from Fig. 5. Since these configurations cannot produce final spin-orbit angles compatible with the observed ψ•, they can be discarded as viable solutions, though we include their retrieval results in Table 4 for completeness.

To complement our analysis, we investigated whether mirror configurations with respect to i_{mut, 0} = 90° would yield similar results. The ZLK equations at the quadrupolar level exhibit symmetry with respect to i_mut,0 = 90°, suggesting that a system with initial mutual inclination i_mut,0 and another with 180° - i_mut,0 should follow identical evolutionary paths for most orbital elements during the resonance phase (Sidorenko 2018; Lei & Gong 2022). We tested this theoretical expectation by selecting 5 000 random configurations for each prograde i_mut,0 value (60°, 70°, 75°, 80°, and 85°) that successfully triggered a ZLK resonance and subsequently decoupled within the simulation timeframe. We then simulated their mirror configurations with retrograde i_mut,0 values (120°, 110°, 105°, 100°, and 95° respectively), totaling 25 000 additional simulations. Our numerical verification confirms the expected symmetry for most initial mutual inclinations during the resonance phase. The exception is at i_mut,0 = 60° (and its mirror at 120°), where we occasionally observe substantial differences - some systems would decouple while their mirrors remain locked in ZLK cycles throughout the simulation, for instance. This discrepancy results again from the more chaotic nature of lower-inclination configurations where SRFs become prevalent even during the coupled phase.

Furthermore, while the symmetry holds during the resonance, it breaks after decoupling as the subsequent evolution becomes dominated by tidal effects, whose governing equations lack this symmetry. Despite this, our analysis reveals that the distributions of final spin-orbit angles for mirror configurations are remarkably similar to their prograde counterparts, as shown in Fig. 5. While individual mirror systems do not necessarily produce identical final spin-orbit angles (i.e., we observed some jitter between paired systems and even some cases where a system with ψ_f ≃ 60° had a mirror counterpart with ψ_f ≃ 120° and vice versa), these differences average out when examining the global distribution. Kolmogorov-Smirnov tests confirm this similarity with p-values of >0.5, indicating no statistically significant difference between ψ_f distributions of mirror systems, for all pairs except i_mut,0 = 60° and 120°. Notably, configurations with i_mut,0 = 120° produced no retrograde final orbits, consistent with their prograde counterparts at 60°. Thus, they were also excluded from Fig. 5.

Based on our analysis of final spin-orbit angles (Fig. 5), we can now confidently identify the range of initial mutual inclinations that best satisfy the observational constraint on ψ_•. Configurations with i_mut,0 = 80° and 85° (and by extension, their retrograde counterparts at 100° and 95° ) produce substantial numbers of systems with final spin-orbit angles compatible with the observed value. We therefore conclude that initial mutual inclinations in the range 80°-100° most successfully reproduce GJ 436 b’s observed polar orbit while maintaining consistency with all other observational constraints. We did not directly explore initial mutual inclinations beyond 85° (or below 95° for retrograde configurations) because we observed increasing discrepancies between JADE simulations and N-body codes (Trani & Spera 2023) near i_mut,0 ≃ 90°. Nevertheless, the clear trend of increasing retrograde final orbits as i_mut,0 approaches 90° from either direction strongly suggests that this central value would also be compatible with observations.

Fig. 5 Distribution of final spin-orbit angles as a function of different initial mutual inclinations, for simulations where a resonance was triggered and then exited. The orange shaded region is the observational constraint on ψ.

5 Plausible location of GJ 436 c

5.1 Exploring the detection limits

It is finally worth putting our results in an observational context to pinpoint the possible location of the concealed companion GJ 436 c. Given the decades of RV data acquired for the system, a certain fraction of the parameter space can already be excluded. To perform such an analysis, we retrieved all publicly available HARPS (Mayor et al. 2003) data collected between 2003 and 2020 from the ESO archive⁵, encompassing three distinct observational programs: HARPS GTO (Mayor, 072.C-0488), “Transit of telluric exoplanets orbiting M dwarfs” (Bonfils, 082.C-0718), and “The HARPS search for planets around M dwarfs” (Bonfils, 1102.C-0339). The RVs were extracted using a templatematching code (NAIRA), specifically optimized for M dwarfs, previously applied in multiple HARPS-based exoplanet searches (e.g., Astudillo-Defru et al. 2015). The RVs were then binned per night and corrected from secular acceleration. Our final set consists of 171 nights from HARPS03, spanning 1532 days, and 18 nights from HARPS15, spread over/covering only 41 nights. This results in a total observational baseline of 5157 days, with a 3584-day gap. The offset between the two datasets is determined through the modeling of GJ 436 b. Since both datasets are extracted from different templates, the offset is statistically consistent with zero.

After removing the RV contribution of GJ 436 b (using the planetary parameters from Bourrier et al. 2018), we compute detection limits based on the Generalized Lomb-Scargle method proposed by Zechmeister & Kürster (2009). This involves an injection-recovery test of a sinusoidal signal in the RVs at 12 equispaced phases, for each tested period. For this study, we adopt the most conservative detection limit (see Mignon et al. 2024, for details), defined as the largest injected amplitude required for the periodogram signal to reach the 99% significance level at the period of the injected signal.

We also mention the collection of supplementary observational detection limits stemming from the direct images of NACO (Brandner et al. 2002) from the ESO archive (Apaï, 081.C-0430). However, they were discarded because less stringent than the RV constraints.

5.2 Constraining the parameter space

Figure 6 shows the possible values for the properties of planet c from our exploration. The red hatched region corresponds to the detection limits excluded from archival RV measurements (Sect. 5.1). The bright streak includes the retrieved values of Λ for i_mut,0 values of 80° and 85° (Table 4), as they are the only explored initial mutual inclinations satisfying the observational constraint on the spin-orbit angle (along with their mirror counterparts 100° and 95°, which would yield the same retrieved values of Λ thanks to the symmetry, see discussion in Sect. 4.3). The two compatible values of Λ were merged, assuming their dependence to i_mut,0 is continuous for this high-inclination regime, which seems justified in light of Table 4. Left of the bright streak and for a given a₀, the decoupling occurs too quickly: the orbit is already circularized after t_•↓. To the right of it, we can see the resonance is either too slow (no decoupling after t_•↑) or does not even occur at all, in line with the findings of Beust et al. (2012).

Our compatible streak is consistent with the initial proposal of Bourrier et al. (2018) for long periods. The shape of their solution departs from ours for the lowest periods, where our compatible streak diverges from their results. This divergence occurs because these configurations correspond to low ratios of outer to inner orbital angular momentum (of order a few), a regime where the secular approximation begins to break down (Mangipudi et al. 2022) and higher-order effects become significant. While our quadrupole-only treatment cannot capture such contributions, Beust et al. (2012) showed that the corrected compatible parameter space in this regime would follow $M_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}}\propto P_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}}^{-1/3}$ (locus of constant angular momentum). Crucially, this correction would push the viable solutions at P_pert < 10 yr into the region already excluded by our RV detection limits (as seen in Fig. 6). Therefore, while our quadrupolar approximation does introduce differences compared to Bourrier et al. (2018) at short periods, these differences minimally affect our conclusions about the viable parameter space for GJ 436 c. It is also worth mentioning that our solution is wider compared to Bourrier et al. (2018), which is due to our semi-Bayesian exploration framework yielding uncertainty margins (while they find compatible configurations by trial and error on a coarse grid). Plus, we virtually include the range of i_mut,0 between 80° and 85° (while they fix i_mut,0 = 85°). In reality, our compatible streak should be even larger so as to encompass systems with i_mut,0 up to 90°, which are out of our numerical reach.

Importantly, our compatible parameter space does intersect with the refined RV detection limits at P_pert > 200-300 yr, from which we effectively rule out stellar companions, as they would fall entirely within the excluded region. Even brown dwarf companions appear highly unlikely, as they would require fine-tuned combinations of mass and orbital parameters to remain consistent with both the ZLK constraints and the nondetection in RV data. The viable parameter space thus strongly favors (cold) planetary-mass companions. In any case, there is still a window of opportunity to eventually detect this rogue companion, possibly by future missions, or JWST (Gardner et al. 2006) via direct imaging. Figure 6 could serve for observers as a reliable basis for subsequent planet hunting campaigns.

Fig. 6 Possible parameter space where GJ 436 c could hide in the massperiod plane, assuming a circular orbit. The orange shaded streak is the compatible configuration from the derived values of Λ for imut, 0 values of 80° and 85°. The blue line is the result from Bourrier et al. (2018). The red hatched region, delimited by the red dashed line, is excluded by RV measurements.

6 Discussion and conclusion

The analysis described in this work comes with a certain number of caveats. The main limitation would be the first-order formalism we employed for tidal forces, which could be further refined to include frequency-dependent inertial waves flowing from the rheological properties of the planet (André et al. 2017, 2019; Dhouib et al. 2024). In relation to this, while our calculation of the intrinsic temperature, T_int, includes contributions from atmospheric and core cooling/contraction as well as radiogenic heating from heavy element decay (following tabulated models from Mordasini 2020), it currently neglects the heat generated by tidal friction. The planet’s internal energy budget in JADE thus accounts for three of the main contributions (cooling, radioactive decay, and irradiation), but missing the tidal compo-nent⁶. This omission means that T_int does not spike during the rapid migration phase (as might be expected in Fig. 4), potentially underestimating planet inflation at t_trans. Additionally, the planet structure handling could have been better than imposing default values for most of the parameters (Sect. 3.2), but the exploration tells us that erosion is notably ineffective. A more elaborate prescription would have just resulted in a systematic shift in the radius for all configurations, perhaps slightly changing the inferred values because of the feedback of tides. Finally, our treatment of the stellar bulk parameters is quite crude. Stars experience mass loss, inflation, and contraction along their main sequences, which should be taken into account given the sensitivity of tidal dissipation to such effects. In particular, a more realistic evolution of the spin-orbit angle is anticipated, should these processes be implemented.

In this work, we applied the novel semi-Bayesian exploration framework of JADE to GJ 436 b. It allowed us to constrain its formation location around 0.3 AU and pinpoint the properties of a potential distant companion, identified to be planet-mass, which could explain this enigmatic system’s observed characteristics on the basis of a late-stage HEM scenario. Our exploration reveals that GJ 436 b likely formed far enough from its host star to avoid substantial atmospheric erosion during the star’s highly active youth. The planet’s current eccentric and polar orbit strongly supports a scenario where a hidden companion induced ZLK oscillations with initial mutual inclinations of 80°-100°, eventually driving GJ 436 b to its present-day location. The minimal impact of atmospheric mass loss throughout this history demonstrates that atmospheric evolution and orbital dynamics must be treated simultaneously to accurately reconstruct the evolutionary pathways of close-in planets.

We draw attention to the fact that our framework accounts for the statistical likelihood of observing GJ 436 b in its current state. We acknowledge that our proposed late migration scenario results in the planet matching its observed eccentricity for less than 1% of the system’s lifetime (as seen in Fig. 4). However, this apparent fine-tuning is properly considered through our temporal integration of individual probabilities, expressed via Eq. (9), and also elaborated on in Appendix C. It naturally favors evolutionary scenarios that maintain compatibility with observations over extended periods, while penalizing those that only briefly match the current system parameters. Importantly, we find this late migration scenario substantially more conceivable than the alternative “no-migration” hypothesis. The latter would require an unrealistically high tidal dissipation factor (Q_p ≳ 10⁶, Mardling 2008) to maintain the observed eccentricity over Gyrs, offers no explanation for the polar spin-orbit angle, and would necessitate an implausibly massive initial planet to survive the intense evaporation during the star’s energetic infancy. Forming such massive planets at ~0.03 AU around M dwarfs appears to be disfavored by both theory and observations (Laughlin et al. 2004; Endl et al. 2006; Burn et al. 2021; Schlecker et al. 2022). In contrast, our late migration scenario explains all observed properties, albeit requiring an unseen companion, and set in a configuration that is consistent with current detection limits (Sect. 5).

Nonetheless, the inferred initial configuration of GJ 436 b at ~0.3 AU with a companion on a near-perpendicular orbit raises important questions about the system’s primordial architecture. Such a hierarchical two-planet system with high mutual inclination is certainly not a typical outcome of standard planet formation scenarios, possibly representing the weakest aspect of our hypothesis from a physical standpoint. While we do not directly address the origin of this primordial mutual inclination in this work, several mechanisms in the literature could potentially explain such a configuration. One possibility involves formation in a warped protoplanetary disk with misaligned inner and outer components. In this scenario, a massive planet carving a gap could break the disk into two misaligned sections (Nealon et al. 2018, 2019; Zhu 2019). Nevertheless, this mechanism preferentially tilts the inner disk, conflicting with our assumption of an initially aligned inner orbit, and would require an additional planet not considered in our analysis. Alternatively, planet-planet scattering could produce the required architecture if multiple cold giant planets underwent disruptive dynamical encounters, leaving a survivor on a highly inclined orbit (Rasio & Ford 1996; Chatterjee et al. 2008; Beaugé & Nesvorny 2012; Petrovich et al. 2014). However, the ability to generate near-perpendicular configurations through such events remains uncertain.

Stellar flybys during the early cluster phase could also explain such a configuration by violently tilting the outer orbit (Malmberg et al. 2011), although this would need to occur shortly after disk dispersal to match our inferred initial conditions, as the disk’s presence strongly dampens external perturbations (Picogna & Marzari 2014). Flybys could directly deposit planets on wide orbits as well, following episodes of enhanced scattering (Bailey & Fabrycky 2019), as demonstrated comprehensively by Izidoro et al. (2025) via instabilities that end up leveraging the dense birth environment. Such far-orbiting companions would be dynamically decoupled from inner planets, naturally allowing for substantial mutual inclinations. However, most studies of these cluster-based mechanisms have focused on Sun-like stars and infer typical orbital separations >100 AU. The relevant distances for cluster dynamics scale proportionally as $M_{\star }^{1/3}$ (Hill sphere), potentially bringing the effective range down to >70 AU and closer to our inferred companion orbital distances. Whether these mechanisms can operate efficiently at the closer end of our viable parameter space (a few tens of AU) remains an open question. Encouragingly, a growing number of exoplanet systems exhibit similar architectures, with close-in small planets being highly misaligned in relation to distant giants, including HAT-P11 (Winn et al. 2010; Yee et al. 2018), π Men (Jones et al. 2002; Xuan & Wyatt 2020; De Rosa et al. 2020), and HD 73344 (Sulis et al. 2024; Zhang et al. 2025). These observations confirm that our inferred configuration, while appearing exotic at first glance, could indeed be possible in nature. The challenge now is understanding the potentially multistage processes that produce such systems.

This study has important implications for understanding the broader population of exoplanets at the edge of the hot Neptune desert. The eccentric, misaligned systems forming the “Neptunian ridge” described by Castro-González et al. (2024) might share similar evolutionary histories involving late migration, explaining both their survival and their unusual orbital configurations. GJ 436 b could therefore be a prototype of a distinct population of planets whose histories are shaped by the complex interplay between atmospheric processes and dynamical perturbations. Future works will extend this analysis to a larger sample of Neptune-sized planets to test whether HEM is a common pathway for planets populating the boundaries of the desert (Bourrier et al., under review). While our current constraints on GJ 436 c provide a roadmap for future observational campaigns, direct detection would provide the ultimate validation of this scenario.

Acknowledgements

We are grateful to the anonymous referee for their review, which improved the quality of our manuscript. We offer our thanks to Elsa Bersier (ESBDI, Geneva) for designing the JADE logotype. We thank Aline Vidotto for insightful conversations. The authors made use of the Claude AI assistant (Anthropic, 2024) for code optimization. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project SPICE DUNE, grant agreement No 947634; grant agreement No 730890). This material reflects only the authors’ views and the Commission is not liable for any use that may be made of the information contained therein.

References

Appendix A JADE 101: Use cases

Fig. A.1 Schematic representation of the JADE code’s key components and processes. The system shows a central star (red) with two orbiting bodies: an inner planet (blue) with a gaseous envelope and a solid core, and a distant perturber (gray). The sketch illustrates the evolution of the system across two time steps, demonstrating secular changes in the inner planet’s orbital parameters and spin-axis orientation, while the outer orbit remains fixed. The zoomed insets depict the inner planet’s atmospheric structure, showing core heating (orange) on one hand, and atmospheric escape (purple) processes under stellar irradiation (yellow) on the other hand. The three-dimensional (3D) spin-orbit angle, ψ, is highlighted in light blue. See Attia et al. (2021) for more details.

The typical configuration simulated by JADE is depicted in Fig. A.1. One of the major strengths of the JADE code is its versatility, capable of tackling a wide spectrum of questions. Figure A.2 illustrates the variety of use cases (FD1, FD2, FA1...) that are achievable with it, covered by pedagogical tutorials within the first distribution of the software linked in the first page of this article.

The JADE code was originally designed to identify families of secular behaviors and pinpoint the influential mechanisms at play in the framework of a coupled dynamics-atmosphere treatment. The FD1, FA1, FA3 use cases, and essentially their combination FC (Fig. A.2, left panel) embody this qualitative philosophy. They represent the foundational blocks on which the other use cases are built. Starting from a specified set of orbital and bulk parameters, JADE can forwardly simulate the system’s dynamical evolution including all of the relevant contributions for ZLK-induced HEM (FD1), as well as the close-in planet’s evaporative evolution (FA3), characterize its atmospheric structure at a given time (FA1) by computing various profiles (e.g., temperature-pressure, mass distribution), or amalgamate the above for a fully coupled, realistic evolutionary history (FC).

Running sequentially several simulations pertaining to one of these use cases allows us to consider additional, more quantitative application scenarios. For example, by iterating FD1 over some system parameters in a grid-like structure, we can locate the regions of the parameter space where the ZLK resonance can surpass the action of SRFs and be effectively triggered, which forms the FD2 use case (an implementation is described in Hagelberg et al. 2023). In practice, such “ZLK viability” maps are generated using analytical expressions and later validated with JADE simulations coarsely sampling the parameter space, for efficiency purposes. On the atmospheric side, we can perform multiple simulations where evaporation is the exclusive operating process (FA3) launched on a range of different initial planet masses, assuming a fixed orbit for all of them; namely, the one corresponding to the present-day configuration. In doing so, the initial planet mass yielding a mass compatible with the current one at the age of the system can be interpreted as the maximum initial mass (FA4). Indeed, an unchanged close-in orbit during the entire planet’s lifetime (in situ formation or early-on disk migration) represents the locus of maximum possible erosion, and any initial mass that is too high to evaporate enough to match that of the present-day can definitely be excluded⁷ (an application will be soon available in Livingston et al. in prep. ). Finally, individual snapshots of what the planet structure looks like (FA1) for different sets of bulk and irradiation parameters can be used to construct a grid of atmospheric models (FA2, Sect. 3.3), which can prove helpful in deriving accurate mass-radius relationships (Sect. 3.4) and capitalizing on all the benefits they bring.

Fig. A.2 Overview of the different use cases for the JADE code, structured into forward and retrieval modeling. For every use case, an identification code is specified on its very left, and an icon shows whether it requires a single simulation or multiple ones. The shaded rectangles behind the use cases denote which integrators are involved.

Recently, JADE has been also employed in retrieval mode, meaning that today’s system characteristics allow us access to all sorts of hitherto unknown properties (Fig. A.2, right panel). Indeed, the JADE code is able to constrain a planet’s internal structure (iron core, silicate mantle, and volatile atmosphere mass fractions, RS1) as well as the composition of its H/He-dominated atmosphere (hydrogen X, helium Y, and metal Z mass fractions, RS2) with the observational mass and radius in a Markov chain Monte Carlo (MCMC) approach. This “static” (present-day) bulk characterization often serves as a first step prior to the “evolutionary” parameter space exploration, the latter ultimately enabling us to constrain the close-in planet’s birthplace (RE1) and outer companion (RE2) to reveal credible pathways compatible with the currently measured parameters. The semi-Bayesian infrastructure in which RE1 and RE2 are embedded is detailed in Sect. 4.1. Such an exploration is made possible among others through the upstream reduction of the parameter space to inspect via FD2 and FA4 by ruling out, from the outset, all configurations that quench ZLK cycles on the one hand, and that are too massive to sufficiently erode on the other hand. Critically, even the most limited of the parameter spaces would require a considerable number of simulations for a meaningful exploration due to the chaotic nature of three-body dynamics and, what is more, the nonlinear intervention of atmospheric feedback. This is precisely where FA2 intervenes so as to dramatically decrease the computation time: the atmospheric properties, in particular the self-consistently derived planet radius, are interpolated from a precomputed grid, bypassing the usage of the structure integrator, which is by far the biggest computational bottleneck of JADE. This way, a typical fully coupled run can be made in around ten minutes, authorizing large explorations on computer clusters for instance. The quantitative performance tests of such an approach are presented in Appendix B, and its results on our flagship system, GJ 436, are shown in Sect. 4.

Appendix B Performance of the mass-radius relationships

Fig. B.1 Performance assessment of the derived mass-radius relationship for GJ 436 b. Top: histogram of residuals. Are also shown the median, the 95%, and 99% quantiles. Bottom: histogram of run time gains. Are also shown the median, the 5%, and 1% quantiles.

Even though the analytical radii of Sect. 3.4 seem to be consistent with the ones directly generated by JADE upon visual inspection (Fig. 3), it is necessary to quantify their deviation to assess the performance of our mass-radius relationships. To this effect, a control sample of 10000 (M_p, T_eq, T_int) 3D-points are randomly drawn from the same domains of definition set for the internal structure grid (Sect. 3.3). For each point, the analytical radius R_{p, ana} is calculated from Eq. (3) and then compared to the “true” radius R_{p, mod}, modeled by the structure integrator of the JADE code, from which a residual metric is derived using the unsigned relative error between them $$\mathrm{R}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}=|\frac{R_{\mathrm{p},\mathrm{m}\mathrm{o}\mathrm{d}}-R_{\mathrm{p},\mathrm{a}\mathrm{n}\mathrm{a}}}{R_{\mathrm{p},\mathrm{m}\mathrm{o}\mathrm{d}}}|.$$(B.1)

Figure B.1 (top panel) shows the histogram of residuals for the 10 000 random points, manifesting an excellent accuracy. The median residual is 0.6%, which is very satisfying compared to the typical residuals such mass-radius relationships yield (usually considered satisfactory if lower than 5%). 95% (respectively 99%) of the control sample has a better residual than 1.9% (respectively 3.9%), meaning that the error resulting in using Eq. (3) instead of the integrator is always safely negligible. Further enhancement, the run time gain when computing R_p,ana versus R_p,mod is utterly dramatic (Fig. B.1, bottom panel). In 99% of cases, the former is faster than the latter by a factor of 10^3.1 ≃ 1 300, and on average, by a factor of 10^3.9 ≃ 8 000.

Appendix C Posterior distribution of the constrained parameters

We use ${\theta }_t^{\mathrm{I}}$ to denote the orbital parameters of the inner parameters at the time of the observations and y represents the data. The posterior distribution is $$p({\theta }_t^{\mathrm{I}}\mid y)=\frac{p(y\mid {\theta }_t^{\mathrm{I}})p({\theta }_t^{\mathrm{I}})}{p(y)},$$(C.1)

where $p({\theta }_t^{\mathrm{I}})$ is the prior distribution on the orbital elements. Eq. (C.1) gives the classical posterior distribution of the orbital elements fitted onto the data. In principle, the same could be done to obtain the posterior distribution of the initial position of the inner and outer planets (θ^I, θ^O), as well as the time from the initial starting point t

$$p({\theta }^{\mathrm{I}},{\theta }^{\mathrm{O}},t\mid y)=\frac{p(y\mid {\theta }^{\mathrm{I}},{\theta }^{\mathrm{O}},t)p({\theta }^{\mathrm{I}},{\theta }^{\mathrm{O}},t)}{p(y)}.$$(C.2)

Indeed, for each θ^I, θ^O, we could propagate the equation in time to get the estimated current properties of the inner planet Φ(θ^I, θ^O, t), setting ${\theta }_t^{\mathrm{I}}$ in Eq. (C.1) as this propagated value.

This approach is approximated here by considering ${\theta }_t^{\mathrm{I}}$ as the data, which follows a Gaussian distribution. Denoting by n the number of components of ${\theta }_t^{\mathrm{I}}$, denoting the uncertainty on each component k by σ_k and assuming that there is no correlation between the parameter components, the likelihood then becomes $$p({\theta }_t^{\mathrm{I}}\mid {\theta }^{\mathrm{I}},{\theta }^{\mathrm{O}},t)=\frac{1}{\sqrt{2\pi }\prod _{k=1}^n{\sigma }_k}\exp (-\frac{1}{2}\sum _{k=1}^n\frac{{({\theta }_{t,k}^{\mathrm{I}}-\mathrm{\Phi }({\theta }^{\mathrm{I}},{\theta }^{\mathrm{O}},t{)}_k)}^2}{{\sigma }_k^2}).$$(C.3)

The likelihood function employed in this work (Eq. (8)) is none other than Eq. (C.3) marginalized over time.

We further add a prior on θ^I, θ^O, t (uninformative in this study), so that their posterior distribution is approximated by $$p({\theta }^{\mathrm{I}},{\theta }^{\mathrm{O}},t\mid {\theta }_I^t)=\frac{p({\theta }_t^{\mathrm{I}}\mid {\theta }^{\mathrm{I}},{\theta }^{\mathrm{O}},t)p({\theta }^{\mathrm{I}},{\theta }^{\mathrm{O}},t)}{p(y)}.$$(C.4)

Integrated over time, Eq. (C.4) is used to constrain the different quantities described in Sect. 4.

Appendix D Additional figures

Fig. D.1 Corner plot for the jump parameters (a0, Mp,0, Λ) of the fiducial exploration (imut,0 = 75°, Sect. 4.2). Orange lines indicate the medians of the posterior distributions, and the orange shaded regions their 68% highest-density intervals (HDI). Their numerical values are reported on top of the 1D histograms. Here, log denotes the natural logarithm.

Fig. D.2 Distribution of final spin-orbit angles as a function of the jump parameters, for the simulations where a resonance was triggered and then exited in the fiducial exploration (imut,0 = 75°, Sect. 4.2). The orange shaded region is the observational constraint on ψ. An aggregated histogram of the final spin-orbit angles is also shown on the far right.

All Tables

All Figures

	Fig. 1 Rosseland mean opacity, κR, as a function of the temperature, T, interpolated from the tables of Ferguson et al. (2005). Left: varying the metallicity, Z, at a fixed hydrogen mass ratio X = 0.8. Right: varying the hydrogen mass ratio, X, at a fixed metallicity Z = 0.0001.
In the text

	Fig. 2 Corner plot for the fitted coefficients a, b, c, and d from Eq. (3). Each black point is the result of fitting one internal structure model from the precomputed grid (Sect. 3.3).
In the text

	Fig. 3 Analytical mass-radius relationships for GJ 436 b. The crosses represent some radii directly computed by the structure integrator to emphasize the goodness of the fit. Left: varying equilibrium temperature, Teq, with a fixed intrinsic temperature Tint = 100 K. Right: varying intrinsic temperature Tint, with a fixed equilibrium temperature Teq = 500 K.
In the text

	Fig. 4 Secular evolution of the system corresponding to the maximum log-likelihood of the fiducial exploration (imut,0 = 75°, Sect. 4.2), as simulated by JADE. The orange shaded regions represent the observational constraints of Table 3. In all plots, ttrans is shown as a dashed vertical red line.
In the text

	Fig. 5 Distribution of final spin-orbit angles as a function of different initial mutual inclinations, for simulations where a resonance was triggered and then exited. The orange shaded region is the observational constraint on ψ.
In the text

	Fig. 6 Possible parameter space where GJ 436 c could hide in the massperiod plane, assuming a circular orbit. The orange shaded streak is the compatible configuration from the derived values of Λ for imut, 0 values of 80° and 85°. The blue line is the result from Bourrier et al. (2018). The red hatched region, delimited by the red dashed line, is excluded by RV measurements.
In the text

Fig. A.1 Schematic representation of the JADE code’s key components and processes. The system shows a central star (red) with two orbiting bodies: an inner planet (blue) with a gaseous envelope and a solid core, and a distant perturber (gray). The sketch illustrates the evolution of the system across two time steps, demonstrating secular changes in the inner planet’s orbital parameters and spin-axis orientation, while the outer orbit remains fixed. The zoomed insets depict the inner planet’s atmospheric structure, showing core heating (orange) on one hand, and atmospheric escape (purple) processes under stellar irradiation (yellow) on the other hand. The three-dimensional (3D) spin-orbit angle, ψ, is highlighted in light blue. See Attia et al. (2021) for more details.

In the text

	Fig. A.2 Overview of the different use cases for the JADE code, structured into forward and retrieval modeling. For every use case, an identification code is specified on its very left, and an icon shows whether it requires a single simulation or multiple ones. The shaded rectangles behind the use cases denote which integrators are involved.
In the text

	Fig. B.1 Performance assessment of the derived mass-radius relationship for GJ 436 b. Top: histogram of residuals. Are also shown the median, the 95%, and 99% quantiles. Bottom: histogram of run time gains. Are also shown the median, the 5%, and 1% quantiles.
In the text

	Fig. D.1 Corner plot for the jump parameters (a0, Mp,0, Λ) of the fiducial exploration (imut,0 = 75°, Sect. 4.2). Orange lines indicate the medians of the posterior distributions, and the orange shaded regions their 68% highest-density intervals (HDI). Their numerical values are reported on top of the 1D histograms. Here, log denotes the natural logarithm.
In the text

	Fig. D.2 Distribution of final spin-orbit angles as a function of the jump parameters, for the simulations where a resonance was triggered and then exited in the fiducial exploration (imut,0 = 75°, Sect. 4.2). The orange shaded region is the observational constraint on ψ. An aggregated histogram of the final spin-orbit angles is also shown on the far right.
In the text