Inflated Eccentric Migration of Evolving Gas Giants I – Accelerated Formation and Destruction of Hot and Warm Jupiters

The formation of gas-giant planets is thought to occur either through the gradual bottom-up growth of planetary embryos, followed by rapid gas accretion (core accretion), or through the direct formation of gas giants, following gravitational instability (Mizuno 1980; Boss 1997; Armitage 2010). However, the observed proximity of hot and warm Jupiters (HJs and WJs, correspondingly) to their host star sets severe constraints on both these mechanisms. The high temperatures, high velocities, low disk masses, and solids that characterize these environments constrict the predicted formation rates from the in situ channel (Bodenheimer et al. 2000; Rafikov 2005). Though in situ formation may still potentially explain a non-negligible fraction of the population of WJs residing farther out (Huang et al. 2016; Anderson et al. 2020), the formation of HJs and eccentric WJs is not likely to occur in situ. These have therefore been suggested to form at larger separations from the star, before migrating inward to their current locations (see Dawson & Johnson 2018 for a review).

The currently suggested migration channels can generally be divided between disk migration channels and high-eccentricity tidal migration channels (Dawson & Johnson 2018). While such migration scenarios could potentially explain the origins of HJs and WJs, detailed studies of these migration scenarios have struggled to reproduce the rates and properties of the formed HJs and WJs (Dawson & Johnson 2018; Zhu & Dong 2021). However, such models usually did not self-consistently account for role of the thermal evolution of the gas giants, from their initially hot inflated state to their later contraction (and possible heating), during their eccentric migration.

In high-eccentricity migration mechanisms, the semimajor axis shrinks due to energy dissipation induced by tidal forces. This process requires initial high eccentricity, as tidal coupling strongly depends on the distance from the star, which becomes small at pericenter approaches of highly eccentric orbits. The planet can initially be excited to such high eccentricities by planet–planet scattering (Rasio & Ford 1996; Chatterjee et al. 2008; Jurić & Tremaine 2008), through the von-Ziepel–Lidov–Kozai (ZLK) mechanism (von Zeipel 1910; Kozai 1962; Lidov 1962) coupled to tidal evolution (e.g., Kiseleva et al. 1998; Wu & Murray 2003; Fabrycky & Tremaine 2007; Naoz et al. 2011; Petrovich 2015a; Grishin et al. 2018; Vick et al. 2019), or coupled multiplanet secular evolution (e.g., Wu & Lithwick 2011; Beaugé & Nesvorný 2012; Hamers et al. 2017), or their combination. In this paper, we focus on excitation by planet–planet scattering, and we will discuss the secular channel in a future study. The current tidal migration models show that the typical migration timescales could be long, and the production rates of HJs and WJs (or both) are too low in comparison with observations (particularly for WJs). However, as we discuss below, the initial hot inflated state of a gas giant at birth, and its later radiative and tidal heating, could change this picture.

Core accretion formation of gas giants proceeds through the runaway accretion of gas from the protoplanetary disk typically onto a few Earth-mass solid cores (Perri & Cameron 1974; Bodenheimer & Pollack 1986). At the end of the core accretion process, the recently formed gas giants reach large radii, possibly up to 10 R_J for dusty planets (Ginzburg & Chiang 2019). After a very rapid contraction phase, when the planet's radius contracts to 4 R_J (Guillot et al. 1996), the contraction slows down, and takes place at typical Kelvin–Helmholtz timescales (tens of Myrs), as the planet contracts asymptotically to reach its (effective) final size. The contraction depends on the mass and, as we discuss in depth, on the externally injected energies, due to tides and/or irradiation, with only weak dependence on the exact initial conditions at birth (leaving aside the uncertainties/degeneracy in the models of the early evolution of gas giants; e.g., Marley et al. 2007). It should be noted that dustier planets could lead to even slower contraction (Ginzburg & Chiang 2019).

Typically, eccentric migration scenarios assume that the radii of the migrating gas giants are constant throughout their evolution, and are taken to be their asymptotic radii at late times, ∼ 1R_J , neglecting the initial inflated radii at birth and the later contraction. Since tidal forces strongly depend on the radius of the affected object, initially inflated planets should migrate faster than planets that have already contracted to their final radii. The dynamical and thermal evolutions cannot be decoupled, due to their mutual strong dependence on each other. To date, these issues have only been partially studied by Wu et al. (2007), Miller et al. (2009), and Petrovich (2015a), and even these studies have only addressed limited aspects of coupled evolution models, using simplified approaches.

In this paper and Glanz et al. (2021; hereafter, paper II), we study the coupled dynamical–thermal evolutions of HJs and WJs, which evolve via inflated eccentric migration, after the end of the core accretion phase, and compare them with corresponding cases of initially noninflated planets and the inefficient heat conduction that is typically studied in the literature. We present a semi-analytical model, described in this paper, and compare it with a numerical model, discussed in paper II. Both approaches couple the dynamical and tidal evolutions of the planets with their thermal evolutions. The former makes use of a simple analytic approach for the thermal evolution and its effect on the radius evolution, and the latter follows the evolution through a numerical model, using the stellar and planetary evolution code MESA (Paxton et al. 2011, 2013), coupled with the AMUSE framework (Portegies Zwart et al. 2009). We find excellent agreement between the model results, in terms of the overall evolutions and properties of the systems, as we discuss below and in paper II. Consequently, and due to the simplicity of the semi-analytical approach relative to the computationally costly numerical one, we discuss the results from the numerical evolution only for specific cases, and show comparisons to the analytic model (see below and in paper II). The efficient semi-analytical approach allows us to study both the large parameter space of the possible initial conditions and the formed population of HJs/WJs, which are inaccessible to the numerical approach, since it is computationally expensive and has limitations in terms of long timescales/small radii (see paper II for further details).

In the following, we describe our semi-analytical approach, and present the various aspects of inflated eccentric migration and its outcomes. We then consider the overall distribution of the initial conditions, the resulting rates and properties of the formed HJs and WJs, and their dependence on these initial distributions.

In order to present our approach, we first discuss the conditions for inflated eccentric migration arising from planet–planet scattering. We then discuss tidal migration in general in Section 3, and elaborate the weak and dynamical tidal models and their use in our models. In Section 4, we describe the semi-analytical model for inflated migration and demonstrate the formation of HJs and WJs in specific examples. In Section 5, we demonstrate the use of the semi-analytical model in several examples, and compare to the numerical model. In Section 6, we discuss the population synthesis of the formation of HJs and WJs. In Section 7, we present the results of the population synthesis, the choice of the parameters, and the role played by them. In Section 8, we discuss our findings and the further implications. In Section 9, we summarize our results.

In this study, we focus on planetary systems in which the initial conditions for high-eccentricity migration are dictated by planet–planet scattering. There are other channels for eccentricity excitation, such as secular Lidov–Kozai (LK) evolution in triple systems (e.g., Wu & Murray 2003; Naoz et al. 2011; Petrovich 2015a, 2015b; Grishin et al. 2018; Vick & Lai 2018) and secular resonances (e.g. Wu & Lithwick 2011; Hamers et al. 2017). Such processes typically operate on longer timescales, lead to intermittent/quasiperiodic high eccentricities, and are also sensitive to precession induced by tidal interactions, which can partially quench the level of eccentricity excitation. In this work, we consider only planet–planet scattering, but inflated eccentric migration is likely to be important for those other channels, such as the ones involve secular evolution (see Wu et al. 2007; Petrovich 2015a). It should be noted that in general the timescales of these channels could be longer, suggesting a priori that the effect of the initial inflation may be smaller. However, initially inflated planets will still leave a signature on the dynamical evolution. First, a larger fraction of the initial planets will be prone to tidal disruption. Moreover, the final distribution of the semimajor axes is expected to change accordingly (Petrovich 2015a). Another potential effect of initial inflation is a tidal quenching of the LK mechanism, which passes with the contraction of the planet. Although several studies have been done on the secular channel, to our knowledge the coupling there between the thermal and dynamical evolution is not self-consistent; however, self-consistent modeling could easily be implied in our semi-analytical model. This channel is out of the scope of this paper, and is left for follow-up studies.

In multiplanet systems, mutual gravitational interactions between the planets perturb their orbits and may destabilize the system, leading to the ejection of planets, mutual collisions, collisions with the star, and more general excitation of the planets' eccentricities and inclinations (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Chatterjee et al. 2008). The gravitational encounters could also involve tidal circularization during the scattering process, which might enhance the fractions of shorter-period and high-eccentricity gas giants (Nagasawa et al. 2008). More generally, it has been found that strong planet–planet scatterings eventually give rise to a Rayleigh distribution of the eccentricities of the surviving planets (e.g., Jurić & Tremaine 2008). It should be noted that although good fits exist, the exact final eccentricity distribution is still unknown in analytic terms, and usually a preliminary N-body simulation is needed to determine the distribution at the end of the planet–planet scattering. In our calculations, we assume that destabilized systems evolve through planet–planet scattering, leading to the eccentricity distribution given by

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{de}}\propto e\exp \left[-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{e}{2{\sigma }_{e}^{2}}\right)}^{2}\right],\end{eqnarray} \tag{ 1 }$$

where we adopt a value of σ_e = 0.5 in our fiducial model (see Table 1 and Section 6). We note that some studies have considered somewhat different distributions, with even higher fractions of highly eccentric orbits, which might improve our results (e.g., Nagasawa et al. 2008; Carrera et al. 2019).

The Rayleigh distribution of the eccentricities arising from the planet–planet scattering determines the fractions of the planets that would evolve through inflated eccentric migration to become HJs/WJs, and would not be tidally disrupted by the host star.

The tidal disruption radius is given by ${r}_{\mathrm{dis}}=\eta {R}_{p}{\left({M}_{\star }/{M}_{p}\right)}^{1/3}$, where η = 2.7 (Guillochon et al. 2011), and planets with pericenters approaching these values are assumed to be disrupted and, naturally, are not considered to be HJs/WJs. Planets with too high pericenter approaches, which are not affected by tides, or those that do migrate through inflated eccentric migration, but do not become HJs or WJs in the relevant time considered, i.e., the age considered for the HJs/WJs, are similarly not considered to be HJs and WJs.

The timescale for the planet–planet scattering sets the initial radius of the planet in the migration stage. Gas giants usually form with a typical radius that can exceed even 4 R_J , from which there will be a rapid cooling phase up to 4 R_J (Guillot et al. 1996). The following contraction phase is slower, along the Hayashi track. This timescale could take less than a Myr, to reach a typical radius of 1.5 R_J . The decoupling timescale from the planet–planet scattering could stray from less then a Myr to a few Myrs, and even more (e.g., Dawson & Johnson 2018, Figure 3), such that the initial radius for the migration phase could be large (corresponding to a short decoupling time) or small (corresponding to a long decoupling time), depending on the planet–planet scattering conditions. A more detailed calculation, which sets out more accurately the initial distribution for our semi-analytical-based population synthesis, is out of the scope of this paper, and is left for a future study. However, we do take this into account by considering several possible initial radii distributions.

Tidal migration is a dissipative process, where the tides raised on the planet by the host star extract energy from the planet's orbit, typically leading to its inward migration into shorter-period orbits. Tides raised on the star by the planet can also contribute to tidal migration, but these are typically negligible compared to the effects of the tides raised on the planet, though they might become important under some circumstances (Ginzburg & Sari 2017). In the following, we only consider the tides raised on the planet, and postpone the consideration of the tides on the star to later studies. These tides on the star will lead to corrections in Equation (3), which may be included (e.g., Miller et al. 2009).

In order for the tidal dissipation to be effective, and lead to a significant migration, a small pericenter approach should be considered. This dictates the basic initial conditions for eccentric migration. High-eccentricity tidal migration could therefore be roughly divided into two separate steps: reducing the planet's angular momentum and reducing the planet's energy. In the first step, the HJ/WJ progenitor is excited into an eccentric orbit via planet–planet scattering (Rasio & Ford 1996; Chatterjee et al. 2008; Jurić & Tremaine 2008), as we discuss here, or through other channels for eccentricity excitation (e.g., via the ZLK mechanism; von Zeipel 1910; Kozai 1962; Lidov 1962; Wu & Murray 2003; Fabrycky & Tremaine 2007; Naoz et al. 2011; Petrovich 2015a). In the second step, energy extraction via tides leads to the migration and circularization of the planet's orbit. The energy extracted from the orbit (calculated over an orbital period) per period is dissipated in the planet, and therefore the injected energy heats the planet. The injected energy per unit time is given by

$$\begin{eqnarray}&&{L}_{\mathrm{tide}}=-\displaystyle \frac{E}{a}\displaystyle \frac{{da}}{{dt}},\end{eqnarray} \tag{ 2 }$$

where E is the orbital energy and a is the semimajor axis. The angular momentum is approximately conserved. For HJs, the final orbit is usually circular, and given the angular momentum conservation one can estimate the final semimajor axis of the HJ, finding ${a}_{\mathrm{final}}={a}_{0}(1-{e}_{0}^{2})$, where a₀ and e₀ are the initial semimajor axis and the eccentricity, correspondingly. The same consideration could be applied to any other given final eccentricity (or any eccentricity during the evolution).

Modeling tides is not trivial; in particular, their strong dependence on the internal structure of the planet and other physical aspects of the problem raise many complications. Here, we consider two tide models: one is the widely used tidal model of weak/equilibrium tides (Darwin 1879; Goldreich & Soter 1966; Alexander 1973; Hut 1981), the second one is a model for dynamical tides (Zahn 1977; Mardling 1995a, 1995b). The latter could be especially important, and more efficient during the early migration phases when the planet's orbit is still highly eccentric; and, in that sense, considering only the weak tides model is potentially conservative in terms of the efficiency of the eccentric migration and the long timescales that result (Lai 1997). Our approach is general, and any other tide model could potentially be incorporated, such as the chaotic dynamical tides model (Ivanov & Papaloizou 2004, 2007; Vick & Lai 2018; Wu 2018; Vick et al. 2019), which can potentially further shorten the migration timescales, due to its more efficient extraction of energy. The tide models are discussed below. Here, we consider only nonchaotic dynamical tides. It should be noted that the very same equations are used in the detailed numerical model in paper II. To keep this paper as succinct as possible, and to minimize the overlap with paper II, we refer the reader to paper II for further discussion of tidal evolution, which is also relevant to the semi-analytical model.

3.1. Equilibrium Tide Model

The equilibrium tide model is used in many physical scenarios and has been widely discussed (e.g., Darwin 1879; Goldreich & Soter 1966; Alexander 1973; Hut 1981).

Let us consider an orbit-averaged time evolution of the eccentricity and semimajor axis. If we assume that pseudosynchronization of the planetary spin and the orbit occurs on a short timescale, and that the angular momentum is conserved during the migration, one finds that (e.g., Hut 1981; Hamers & Tremaine 2017)

$$\begin{eqnarray}&&\displaystyle \frac{{da}}{{dt}}=-21{k}_{\mathrm{AM}}{n}^{2}{\tau }_{p}\displaystyle \frac{{M}_{\star }}{{M}_{p}}{\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{5}{{ae}}^{2}\displaystyle \frac{f(e)}{{\left(1-{e}^{2}\right)}^{15/2}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{de}}{{dt}}=-\displaystyle \frac{21}{2}{k}_{\mathrm{AM}}{n}^{2}{\tau }_{p}\displaystyle \frac{{M}_{\star }}{{M}_{p}}{\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{5}e\displaystyle \frac{f(e)}{{\left(1-{e}^{2}\right)}^{13/2}},\end{eqnarray} \tag{ 4 }$$

where M_⋆ is the mass of the host star, M_p , R_p , e, a, and n are the mass, radius, orbital eccentricity, orbital semimajor axis, and mean motion of the gas giant, correspondingly, ${\tau }_{p}=0.66\ \sec$ is the planetary tidal lag time, and k_AM = 0.25 is the planetary apsidal motion constant (k_AM and τ_p are taken from Hamers & Tremaine 2017). f(e) is defined by

$$\begin{eqnarray}&&f(e)=\displaystyle \frac{1+\tfrac{45}{14}{e}^{2}+8{e}^{4}+\tfrac{685}{224}{e}^{6}+\tfrac{255}{448}{e}^{8}+\tfrac{25}{1792}{e}^{10}}{1+3{e}^{2}+\tfrac{3}{8}{e}^{4}}.\end{eqnarray} \tag{ 5 }$$

The energy extracted per period, and hence the migration and circularization rate, scales as R_p ⁵. Consequently, the migration timescales of initially inflated gas giants, where the planetary radius R_p > R_J , are shortened relative to the migration timescales of noninflated gas giants, i.e., with a constant radius of R_J . The contraction timescales are long enough to maintain inflated gas giants along a significant part of their dynamical evolution, such that the initial radius of an HJ/WJ will leave a signature on its expected final parameters, which could be also observed.

3.2. Dynamical Tides

Tidal forcing from the star might excite the internal energy modes of the planet (mainly the fundamental f-mode), which might induce an enhanced response (Mardling 1995a, 1995b; Lai 1997; Ogilvie 2014). The energy is mostly extracted during the pericenter approach, and the extraction is more efficient compared with the equilibrium tide model, potentially leading to the even more rapid circularization and migration of the planet. The eccentricity decay is accompanied by pseudosynchronization with the angular frequency of the star, and the excitation of the oscillations in the planet becomes less pronounced as the orbital eccentricity decreases. The energy dissipation by the various modes is gradually suppressed, until the transition to the regime in which the equilibrium tides become more dominant. The quadrupole order of the energy dissipation can be written as follows (Press & Teukolsky 1977; Moe & Kratter 2018):

$$\begin{eqnarray}&&{\rm{\Delta }}E={f}_{\mathrm{dyn}}\displaystyle \frac{{M}_{\star }+{M}_{p}}{{M}_{p}}\displaystyle \frac{{{GM}}_{\star }^{2}}{{R}_{p}}{\left(\displaystyle \frac{a(1-e)}{{R}_{p}}\right)}^{-9},\end{eqnarray} \tag{ 6 }$$

with f_dyn = 0.1, unless stated otherwise (Moe & Kratter 2018), being taken for our case of the tidal response of a gas giant. Combining this prescription with the equations of orbital energy and angular momentum, and assuming a constant pericenter, leads to the following equations of the orbital semimajor axis and eccentricity along the migration (Moe & Kratter 2018):

$$\begin{eqnarray}&&\displaystyle \frac{{da}}{{dt}}=\displaystyle \frac{a}{P}\displaystyle \frac{{\rm{\Delta }}E}{E},\displaystyle \frac{{de}}{{dt}}=\displaystyle \frac{1-e}{a}\displaystyle \frac{{da}}{{dt}}.\end{eqnarray} \tag{ 7 }$$

While dynamical tides dominate for large eccentricities, weak tides constitute a more physical description for low ones. The ratio of the migration rate due to dynamical tides to the migration rate due to weak tides is given by:

$$\begin{eqnarray}&&\beta ({R}_{p},a,e):= \displaystyle \frac{{da}/{dt}{| }_{\mathrm{dyn}}}{{da}/{dt}{| }_{\mathrm{weak}}}=\displaystyle \frac{2{f}_{\mathrm{dyn}}{R}_{p}^{3}{\left(1-{e}^{2}\right)}^{15/2}}{21{{GM}}_{p}{k}_{\mathrm{AM}}{\tau }_{p}{\left(1-e\right)}^{9}{{Pe}}^{2}f(e)},\end{eqnarray} \tag{ 8 }$$

where P is the period of the planet. The transition between the dynamical and weak tides occurs roughly at β ∼ 1, and we set a lower artificial cutoff at e = 0.2, such that the transition occurs at $\max \{0.2,e{| }_{\beta =1}\}$, to avoid the divergence of dynamical tides at e = 0.

We model inflated eccentric migration by coupling the orbital equations, governed by the tidal migration model, with the thermal evolution of the planet, dictated by heating due to tides, irradiation, and thermal cooling. The thermal evolution leads to a planetary contraction due to cooling, which can be slowed by external heating sources and, in extreme cases, even be stopped/partially reversed. The tidal evolution is strongly affected by the radius of the planet, such that the thermal evolution changes the migration rate and timescale.

The virial theorem states a relation between the total and potential energies E = − (3γ − 4)U/(3γ − 3), such that $U\,\propto 3{{GM}}_{p}^{2}/R(5-\tilde{n})$, where $\tilde{n}=1/(\gamma -1)$ is the polytropic index, γ is the heat capacity ratio, taken as γ = 5/3, and the proportion constant is determined by the numerical gauge. The change in the planet's luminosity is mostly determined by the change in the thermal energy of the ions, which are not degenerate, and their equation of state is given by the ideal gas equation E = (M_p /μ)k_B T_c , where μ is the molecular weight, taken as the proton mass, and T_c is the central temperature of the planet. Following this relation, we derive the following equation for the radius change due to heating and cooling:

$$\begin{eqnarray}&&\displaystyle \frac{{{dR}}_{p}}{{dt}}\propto \displaystyle \frac{5\gamma -6}{3\gamma -4}\displaystyle \frac{{R}_{p}^{2}}{{{GM}}_{p}^{2}}\left({L}_{\mathrm{ext}}-{L}_{\mathrm{cool}}\right).\end{eqnarray} \tag{ 9 }$$

It should be noted that the energy deposition and loss are generally not constant, and could vary with time and the changes of the orbital parameters—the semi-analytical approach takes this into account. In a similar manner, one can calculate the evolution of the central temperature of the planet:

$$\begin{eqnarray}&&\displaystyle \frac{{{dT}}_{c}}{{dt}}=\displaystyle \frac{{m}_{p}}{{M}_{p}{k}_{B}}\left({L}_{\mathrm{ext}}-{L}_{\mathrm{cool}}\right),\end{eqnarray} \tag{ 10 }$$

where m_p is the proton mass, L_cool describes the cooling of the planet, and L_extra describes the further external sources of energies injected, such as tidal heating and irradiation, or any other general source of heating. The external luminosity could generally be a function of the optical depth τ_dep in which it is deposited, or, equivalently, of the pressure P_dep. The effect of the deposition is strongest for τ_dep = τ_c , i.e., at the center, and it decreases as τ_dep decreases. Although deposition over layers has different physical consequences than deposition at the center (Youdin & Mitchell 2010; Spiegel & Burrows 2013; Komacek & Youdin 2017), for our purposes, the energy injected at outer layers could be translated to a reduced deposition at the center, as shown by Ginzburg & Sari (2016) for the case of a power-law distribution.

To reinflate the gas giant, it can be seen directly from Equation (9) that the minimal external luminosity added should be comparable to the radiating luminosity. For a tidal disruption, a larger energy is required, of the order of magnitude of the binding energy of the planet.

The reduction can be considered by an overall multiplication factor, depending on the depth of the deposition (Ginzburg & Sari 2015, 2016), and the effective temperature could be inferred from the effect of the irradiation of the host star. The exact multiplication factor depends on the not-yet-well-understood heat transfer processes in the planet (which we discuss in more detail below). The planet cools through a blackbody emission, which is given by

$$\begin{eqnarray}&&{L}_{\mathrm{cool}}\propto 4\pi {R}_{p}^{2}{\sigma }_{\mathrm{SB}}{T}_{\mathrm{eff}}^{4},\end{eqnarray} \tag{ 11 }$$

where a correction for blackbody radiation arises from the radiative–convective boundary (RCB) and the transition from isolation to insolation (Ginzburg & Sari 2015). While deposition at the photosphere does not significantly change the cooling rate of the planet, central deposition could lead to more prominent effects. A weak external energy source, such as deposition at the photosphere, could be thought of as "standard" Kelvin–Helmholtz cooling, up to a small correction arising from the external source, which becomes important when L_dep/L_cool ≳ 1. The effective temperature at late times is solely determined by the stellar irradiation, and it reaches a steady state, corresponding with the typical observationally inferred values for the effective temperatures of HJs and WJs. At early times, the temperature of the gas giant could exceed this final temperature, before cooling to a quasi-steady temperature, with its temperature then rising again to the one expected from stellar irradiation (see, for example, Figures 1(a) and 2(a)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While the thermal evolution, and hence also the evolution of the planetary radius, are governed by the central temperature of the planet, the effective temperature, which also appears in the cooling equation, could be changed significantly, without changing the central temperature.

The relation between the effective and central temperatures depends on the pressure gradient inside the planet:

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{\mathrm{eff}}}{{T}_{c}}\propto {\left(\displaystyle \frac{{P}_{\mathrm{RCB}}}{{P}_{c}}\right)}^{\gamma /(\gamma -1)};\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{P}_{c}\sim \displaystyle \frac{{{GM}}_{p}^{2}}{{R}_{p}^{4}},\ {P}_{\mathrm{RCB}}\sim \displaystyle \frac{{{GM}}_{p}}{\bar{\kappa }{R}_{p}^{2}},\end{eqnarray} \tag{ 13 }$$

where P_c is the central pressure of the planet, P_RCB is the pressure at the boundary layer between the radiative and convective regions, and $\bar{\kappa }$ is the mean opacity, averaged over the planet's atmosphere. We consider $\bar{\kappa }=5\times {10}^{-2}\ {\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, following the typical values in the literature (Ginzburg & Chiang 2019). The proportion factors could be larger than the order of unity, and they could be found more precisely from numerical simulations. The photosphere luminosity is determined by the effective boundary condition for L_cool. In general, the RCB layer is dictated by the internal adiabat of pressures. It is not located at the photosphere, and the photosphere is in fact an underestimation of it (see the further discussions in Thorngren et al. 2019 and Sarkis et al. 2021). It should be noted that the pressure of the RCB is dictated by both the internal adiabat and the external irradiation, as shown in Equation (12).

The opacity might vary strongly with temperature changes. For high-temperature planets, i.e., T ≳ 10⁴ K, the opacity is given by the Kramers bound-free opacity law κ ∝ ρ T^−3.5 (H⁻ opacity; Kippenhahn et al. 2012; Ginzburg & Sari 2015). Moreover, dustier planets could have higher initial opacity, and hence finish their core accretion with larger radii, which remain inflated for longer timescales, enhancing the effect discussed in this paper.

Coupling Equations (9) and (10) with orbit-averaged evolution equations defined by the tidal migration model (i.e., for weak tides, this is given by Equation (3), and for dynamical tides, by Equation (7)), provides a complete, consistent semi-analytical description of the coupled dynamical–thermal model.

In the following, we discuss the addition of external heating sources, which are discussed in more detail in paper II.

4.1. External Heating

Our semi-analytical model is general, and could effectively take account of any coupling between the thermal and dynamical evolution, as well as include external heating sources, such as tidal heating and irradiation. The efficiency of the heat deposition due to such sources depends both on the heating rate and on the depth of the deposition. Heat deposition in the central parts of a planet has a more significant effect compared with the limited effect of deposition at the photosphere, where radiative cooling effectively disposes of much of the heating.

External heating plays a dual role, keeping the planet inflated for longer time (and, in extreme cases, inflating it) and modifying the effective temperature. At late stages, the observed effective temperature of a gas giant is determined by the the irradiation/other external energy sources applied on the planet, such that the external heating, even if it does not play a role in inflation, should be taken into consideration in the effective temperature calculation. The relative role of the external heating can change during the planet's thermal–dynamical evolution. When the planet is sufficiently far away from the star, the heating due to irradiation is less than the internal heat, and hence it negligibly affects the thermal evolution. When the planet migrates closer to the star, the radiative heating (and tidal heating, when applicable) increases, and one then needs to consider irradiation when determining the planet's temperature. This can be determined by the planetary equilibrium temperature:

$$\begin{eqnarray}&&{T}_{\mathrm{eff},\mathrm{irr}}(t)\propto {\left(\displaystyle \frac{{L}_{\star }}{16\pi {\sigma }_{\mathrm{SB}}{r}^{2}(t)}\right)}^{1/4},\end{eqnarray} \tag{ 14 }$$

where r(t) is the instantaneous distance of the planet from the star. In order to include it in our model self-consistently, and in an efficient way, we consider the averaged effective temperature over an orbit (given that the heating rate is sufficiently small during dynamical times). Similar considerations would apply to the effects of other heating sources on the effective temperature.

In the following, we present the results of the semi-analytical model, and compare them with the results of the numerical model—we refer the reader to paper II for further examples of this comparison and an extensive description of the numerical approach, its limitations, and its applications.

The semi-analytical model could potentially take into account any general distribution of external energies, by the integration of the contributions from different optical depths (see Ginzburg & Sari 2016 for the case of a constant energy source with a power-law distribution in the optical depth). Currently, for simplicity, the injection of heat is simply given as a point source in the center, while in the numerical modeling one needs to distribute the energy injection over some finite region, so as not to lead to an unstable or divergent result from the nonphysical injection at a singular point. The current chosen distribution in the numerical prescription is a semi-Gaussian (with the dispersion changing with the location of the deposition), which injects heat with a deviation of 15%–50% from the intended value (see paper II for the technical details regarding the distribution). It should be noted that the actual physical distribution is not yet understood, since the external heating sources are likely to lead to a distribution with a preference toward the outer areas of the planet, and its extension depends on the exact nature of the heat transfer. Nevertheless, since the effect of a generally distributed external energy on the dynamical evolution could be translated into some appropriate amount of energy injected at the center, we are practically able to treat the effect of any external energy distribution, through an effective heat injection in the center, with an order-unity constant prefactor to calibrate between the distributed injection in the numerical model and in the semi-analytical one.

In Figure 1, we present an example of the coupled dynamical–thermal evolutions of eccentric 1 M_J gas giants with different initial radii, an initial semimajor axis of 1 au, and an initial eccentricity of 0.98. We compare the evolution of a constant 1 R_J with that of initially inflated planets. All the planets in our model experience photospheric heating induced by stellar irradiation.

We find that the semi-analytical and numerical approaches are in excellent agreement, and yield similar results (see the additional examples in paper II).

The typical migration and circularization timescales are correspondingly given by (for the weak equilibrium tide model):

$$\begin{eqnarray}&&{\tau }_{\mathrm{mig},\mathrm{weak}}=\displaystyle \frac{a}{| {\dot{a}}_{\mathrm{weak}}| }=\displaystyle \frac{{\left(1-{e}^{2}\right)}^{15/2}{M}_{p}}{21{k}_{\mathrm{AM}}{n}^{2}{\tau }_{p}{M}_{\star }{e}^{2}f(e)}{\left(\displaystyle \frac{a}{{R}_{p}}\right)}^{5},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\tau }_{\mathrm{circ},\mathrm{weak}}=\displaystyle \frac{e}{| {\dot{e}}_{\mathrm{weak}}| }=\displaystyle \frac{2}{21}\displaystyle \frac{{M}_{p}{\left(1-{e}^{2}\right)}^{13/2}}{{k}_{\mathrm{AM}}{M}_{\star }{n}^{2}{\tau }_{p}f(e)}{\left(\displaystyle \frac{a}{{R}_{p}}\right)}^{5}.\end{eqnarray} \tag{ 16 }$$

These timescales are somewhat optimistic, since the initial radius shrinks with time. However, it could be seen that initial larger radii lead to shorter migration/circularization timescales.

A larger initial radius shortens the migration and circularization timescales, which become 10 times shorter in this case, as expected from the strong dependence of the tidal forces on the radius. The radius contracts to a final radius of ≳1 R_J within a Hubble time. A significant part of the migration takes place with a radius larger than R_J , since the migration timescale is shorter than the contraction timescale, manifesting the key role played by the initial inflated radius on the evolution. It should be noted that the radius of the initially noninflated planet changes as well, but negligibly. In addition, the effective temperature is roughly constant at the early stages of the evolution, but increases as the planet gets closer to the star, since the irradiation from the star determines the effective equilibrium temperature of the planet.

In Figure 2, we present the coupled thermal–dynamical evolution of a formed WJ and HJ. Similar to Figure 1, the migration timescales are significantly shortened for a planet with the same initial separation, eccentricity, and mass, but a different initial radius.

We can see a good agreement between the semi-analytical model and the numerical one, throughout the evolution. Note that the numerical model terminates due to the limitations of the (older) MESA version used in this regime; for more details regarding the termination criteria of the numerical model, see paper II. In such cases, the numerical model could be extrapolated using the semi-analytical model.

We find that all the WJs produced through inflated eccentric migration are all still eccentric, though somewhat circularized, and that their final states are dictated by the initial angular momentum and the energy dissipation rate. Inflated eccentric migration enhances the migration rate, such that planets that would never otherwise have become HJs or WJs, when considering an initial and constant 1 R_J radii, migrate more efficiently, enabling WJs to become HJs and little-/nonmigrating planets to become WJs. Furthermore, inflated WJs, given the same initial conditions, would be less eccentric, since they proceed faster in their migration; some of the expected WJs from the 1 R_J case become HJs once inflated migration is accounted for. These transitions between regimes, the effective accelerated changes in the migration flows of gas giants, are discussed in more detail in Section 6.

The semi-analytical approach provides an efficient, simple, and computationally inexpensive approach to modeling the evolution, allowing us not only to consider a large phase space of initial conditions, but also to consider a detailed population synthesis study that describes the evolution and formation of Jupiters. Nevertheless, a priori, this simple approach cannot describe the detailed internal structure of a planet, and it might give rise to inaccurate modeling of the macroscopic properties of a planet and its evolution. However, we generally find an excellent agreement between the numerical approach and the semi-analytical one, which therefore allows us to use the simple semi-analytical approach robustly.

In order to calculate the occurrence rate of HJs and WJs and their properties, we use a population synthesis study. Given the excellent agreement that we found between the semi-analytical model and the numerical model (see Figures 1, 2, and the further discussion in paper II), we rely on the semi-analytical model to provide a fast and efficient evolution model, enabling us to study and analyze a wide range of initial conditions, and to explore the formed population of HJs and WJs through a population synthesis study, without the need of a costly (computational) numerical simulation.

We consider various plausible choices for the initial conditions of potential HJ/WJ progenitors (see Table 1), and randomly sample the parameter space of the initial conditions.

In order to characterize the population of HJs and WJs formed through regular eccentric migration and through inflated eccentric migration, we make several assumptions regarding the initial conditions. We also reiterate that in this study we are only considering eccentric migration following planet–planet scattering, while other secular evolution models for eccentricity excitation will be explored elsewhere.

We assume a continuous star formation rate in the Galactic disk, since the planet formation rate is proportional to the star formation rate (e.g., Behroozi & Peeples 2015), and we therefore sample the age (i.e., the evolution time) for each planet from a uniform distribution in the range 1–12 Gyr. We then evolve each planet in our sample using our semi-analytical approach, and examine its properties after a given time. We define an HJ as a gas giant with a final period shorter than 10 days, and we define a (migrating) WJ as a gas giant with a final period of 10–200 days, only considering cases where the semimajor axis has shrunk by at least a factor of 2 relative to the initial one, so as not to consider possible contributions from WJs that have potentially been formed in situ. A gas giant is assumed to be disrupted if its pericenter is smaller than the Roche radius. We define f_HJ/WJ as the fraction of formed HJs/WJs, as described above, and use it to derive the total expected frequency of HJ and WJ systems in the galaxy. The fraction of formed HJs and WJs for any given time t is given by

$$\begin{eqnarray}&&\begin{array}{l}{f}_{\mathrm{HJ}/\mathrm{WJ}}(t)=\displaystyle \int \displaystyle \frac{{dN}}{{{dm}}_{p}}\displaystyle \frac{{dN}}{{{de}}_{0}}\displaystyle \frac{{dN}}{{{da}}_{0}}\displaystyle \frac{{dN}}{{{dR}}_{p,0}}\\ \quad \times \chi ({m}_{p},{e}_{0},{R}_{p,0},{a}_{0},t){{de}}_{0}{{da}}_{0}{{dR}}_{p,0}{{dm}}_{p},\end{array}\end{eqnarray} \tag{ 17 }$$

where ${dN}/{{de}}_{0},{dN}/{{da}}_{0}$, and ${dN}/{{dR}}_{p,0}$ are the differential distributions of the initial eccentricities, the semimajor axes, and the initial planetary radii, respectively. χ is an indicator function determining whether a given system has evolved to become an HJ/WJ after a time t, based on the semi-analytical coupled orbital–thermal evolution. Calculating the fraction of formed HJs/WJs, we also take into consideration the possibility of tidal disruption, and exclude the disrupted population from the population of formed HJs/WJs. The tidal disruption radius is given by ${r}_{\mathrm{dis}}=\eta {R}_{p}{\left({M}_{\star }/{M}_{p}\right)}^{1/3}$, where η = 2.7 (Guillochon et al. 2011), and planets with a pericenter approaching these values are assumed to be disrupted, so, naturally, are not considered to be HJ/WJs. Planets with a sufficiently large pericenter approach, which are negligibly affected by tides, and do not migrate, are also not considered as migration-formed HJs/WJs, nor are those that do migrate due to inflated eccentric migration, but do not migrate close enough to the star to be considered as HJs or WJs in the relevant time considered.

In order to calculate the total occurrence rate, the frequency of HJs and WJs among the stellar systems is given by

$$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{HJ}/\mathrm{WJ}}={f}_{J}\times {f}_{2J}\times {f}_{\mathrm{unstable}}\times {f}_{{ecc},J}\times {f}_{\mathrm{HJ}/\mathrm{WJ}},\end{eqnarray} \tag{ 18 }$$

where f_J is the fraction of stellar systems hosting gas giants; f_2J is the probability of finding at least two gas giants (in most cases, in order to scatter a gas giant into a highly eccentric orbit that leads to eccentric migration, another planet as massive as Jupiter is needed); f_unstable is the fraction of unstable systems, in which planet–planet scattering is likely to occur; and f_ecc,J is the fraction of sufficiently eccentric gas giants that would evolve through eccentric migration to become HJs and WJs.

We study each choice of parameters (see the parameter sampling description below) for the initial distribution of ≳10⁴ semi-analytical simulations to determine the fractions of cases that successfully evolve to become an HJ/WJ. Calculating this at at a given time since birth (up to a specific chosen time) provides us with the delay-time distribution, i.e., the fractions of HJ or WJ as a function of time since planet formation. As time goes by, planets that were observed to be WJs could further migrate to become HJs, HJs could be disrupted, and planets that were too distant to be WJs could become them. The rate at which the different areas in the parameter space are filled is determined by the initial conditions of the planet (a₀, e₀, R₀, and m_p ), the tide model (weak, dynamical, etc.), and the external energy sources that could slow the contraction.

We consider a specific star formation history, and integrate the delay-time distribution (our Green function), weighted by the star formation rate, to obtain a realistic estimate of the current fractions of HJs/WJs in the galaxy at the current time, containing both young and old planetary systems. For disk stars, most relevant for the currently observed exoplanet hosts, we consider a continuous, uniform rate of star formation for the Galactic disk, as mentioned above. This is generally consistent with the inferred local star formation history of the Galactic disk stellar populations where most of the exoplanets have been observed to date. Note, however, that the current exoplanet samples are dominated by those identified by the Kepler mission. The age distribution for Kepler stars peaks at around 2.5 Gyr, and gradually falls off to larger ages (Berger et al. 2020).

In the following, we describe the parameters characterizing the initial conditions, the specific ranges of these parameters, and their motivations. In order to evolve a gas giant, we need to ascribe a planet with both physical and orbital properties, including mass, initial radius, initial separation, and initial eccentricity.

6.1. Parameter Space of the Initial Conditions

6.2. Normalization Factors

Our main results for the HJ/WJ occurrence rates from the various models considered are summarized in Table 1, in which we explore a large parameter space. We will note briefly the different parameters taken into consideration there: σ₁ = 0.5; σ₂ = 0.4; R₁ = 1 R_J ; R₂ = 2 R_J ; R₃ ∼ U[2, 4] R_J ; R₄ = 4 R_J ; a₁ ∼ LU[0.4, 5]au; d relates to the dynamical tides model; w relates to the weak tides model; c_y relates to the deposition of y% of the irradiation at the center and tidal heating (corresponding to the tides model), such that the absence of c_x stands for 0%; and f_x relates to a variation in f_dyn, i.e., f_x stands for f_dyn = x. U stands for a uniform distribution and LU for a log-uniform distribution.

The distributions of the orbital properties are only shown for our fiducial model, i.e., σ₁ R₂ a₁ df_0.1.

We considered a variety of models, which differ in terms of the initial radii distribution of the planets, the initial eccentricity distribution, and the initial semimajor axes distribution. In addition, we also consider different types of tidal evolution, either weak or dynamical, and different efficiencies of external heating (the fraction of heat injected into the center of the planet, which depends on the not-yet-understood heat transfer processes).

7.1. Occurrence Rates

7.2. Parameter-space Evolution

Figure 3 shows the resulting distribution of the eccentricities and orbital periods of our population synthesis models initialized with σ₁ R₃ a₁ df_0.1, where we find the majority of HJs formed via eccentric migration to have been circularized, whereas the WJs are eccentric, suggesting that their migration process has not yet terminated. The empty region on the left reflects the tidal disruption of planets with pericenters below the Roche radius. The empty region on the right reflects the conservation of angular momentum along the migration, such that an initial high eccentricity sets a lower bound on the semimajor axis; planets with too high pericenters are not affected by tides and do not migrate.

Zoom In Zoom Out Reset image size

In Figure 4, we present a histogram describing the final fractions of the possible outcomes of eccentric inflated migration: cold Jupiters, WJs, HJs, or tidal disruption. The population is dominated by cold Jupiters, i.e., gas giants with periods larger than 200 days or gas giants that did not migrate at least half of their initial semimajor axis. Inflated eccentric migration reduces the percentage of cold Jupiters, since the migration is more efficient in this model, leading to a more significant flow in the parameter space from the cold Jupiter regime to the WJs. The fraction of WJs increases, due to inflated eccentric migration, from similar considerations. However, the fraction of HJs also reduces, due to tidal disruption.

Zoom In Zoom Out Reset image size

In Figure 6, we present the eccentricity distribution of HJs and WJs. It can be seen that within this time, HJs tend to circularize, while WJs tend toward higher eccentricities, although they could obtain lower eccentricities, with the lowest at ≲0.1.

7.3. HJ and WJ Populations

HJs and WJs migrate faster via the inflated eccentric migration channel, as also manifested for a specific case in Figure 1. The HJ population formed via inflated eccentric migration can be distinguished from noninflated migrating giants. The distribution of the formed inflated HJs is shifted toward larger periods, as can be seen in Figure 5, in agreement with the findings of Petrovich (2015b), who conducted a population study of the cooling of initially inflated giants (where no external heating was considered in the study) for a secular evolution channel. The larger sizes of the inflated HJs make them susceptible to tidal disruption at larger pericenter approaches, as can be seen directly from the expression for the Roche limit ${r}_{\mathrm{dis}}=\eta {R}_{p}{\left({M}_{\star }/{M}_{p}\right)}^{1/3}$, due to the increased radius and small pericenter.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Indeed, the fraction of tidally disrupted planets increases significantly for models with inflated giants (e.g., see Figure 4). These disruptions arise from the initial conditions, and not the later evolution. Planets are initiated at their largest sizes following their formation, and they can therefore be disrupted at higher pericenter approaches at these times, but then the Roche radius increases as the planets contract.

Only in cases where significant external tidal/radiative heating is able to reinflate planets do tidal disruptions occur during the migration. Indeed, models with significant heating of the central parts (c1, c10) show even higher tidal disruption rates, and, moreover, such disruptions occur during the migration of the planets, following their tidal and radiative inflation, and not immediately after their initial scattering to high eccentricities.

Star formation modifies both the populations of HJs and WJs. As a convolution of single-time star formation events, it gives rise to the further formation of WJs, together with a smaller fraction of HJs, and it indicates that WJs might be younger than HJs, since gas giants that are currently observed as WJs could migrate in and become HJs. We therefore predict that, on average, WJs should reside in younger systems than HJs, if eccentric migration plays a significant role in their production.

7.4. Effects of External Heating on the Population

7.5. Dependence on Parameters

The final population and its properties depend strongly on the choice of the initial distributions and their parameters. While some of the distributions are well constrained from observations, others suffer from large uncertainties, which we account for by considering several choices of parameters. In the following, we discuss several possible choices of parameters and their effects on the final distributions.

In Figure 5, we present the dependence of the final periods on the initial radii distribution. The population formed via the inflated eccentric migration channel peaks in a larger period and gives rise to enhanced filling of the available parameter space of WJs.

As shown in Figure 7, the initial population of more massive gas giants gives rise to a postmigration population residing in smaller periods. This might be expected, given that migration is more efficient for massive planets, as can be seen directly from the tidal evolution equations (e.g., Equations (3) and (7), for weak and dynamical tides, correspondingly), although this is not trivial, given the transitions and flows between cold Jupiters, WJs, HJs, and tidally disrupted planets, which change in each model. It should be noted that the migration timescale of lower-mass planets is shorter than that for more massive ones (see also paper II). However, they are more vulnerable to tidal disruption, such that the overall effect is that lower masses are more efficient in the production of WJs, rather than HJs.

Zoom In Zoom Out Reset image size

The observational mass distribution sets more weight on the less massive gas giants, such that the total period distribution flattens to include a larger range of periods, from HJs to WJs.

In Figures 8 and 9, we show the time evolutions of the gas-giant population. As time goes by, more and more cold Jupiters migrate inward to become WJs, some of them migrate to become HJs, and some will migrate further and be disrupted. This can be seen by following the peak of the distribution. At early times, there is a quick rise in HJs, which form more rapidly; but at later times, WJs form, and the peak gradually moves toward larger periods. The vast majority of the HJs migrate via timescales shorter than 1 Gyr. Considering star formation leads to a continuous flow of formed HJs, and could basically be understood as a time convolution of the distribution derived from a single star formation event. In models without central heating, planets do not reinflate, and all tidal disruptions occur promptly, following the scattering of the planets, while tidal and radiative heating do not inflate the planets, which attain their maximal radii following their formation. When central heating is efficient, planets can reinflate and become larger than their original radii, and be more prone to tidal disruptions at later times. The exact amount of heat needed for reinflation can be roughly estimated by setting the total luminosity to be larger than 0, as can be seen in Equation (9). The effect of reinflation can also be seen from the decreased fraction obtained for gas giants with large amounts of energy injected in their centers (see Table 1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figures 10, 11, and 12, we present the dependence of the migration timescale on the initial semimajor axis and eccentricity, correspondingly. As expected, larger initial semimajor axes lead to larger migration timescales of HJs, and higher initial eccentricities lead to more efficient migration that finally leads to smaller migration timescales. Larger initial radii lead to extremely short migration timescales, while long timescales are cut out of the histogram, due to the elevated disruption rates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 13, we present the dependence of the eccentricity distribution on the dispersion. Lower-eccentricity distributions give rise to larger rates of WJs, but smaller rates of HJs, as expected. Eccentric tidal migration is more efficient, i.e., it extracts more energy from the orbit, when larger eccentricities are included. Since the WJs produced via inflated migration could be thought of as transient HJs that did not manage to end their migration after a given time, their fraction increases when a lower-eccentricity dispersion is taken into consideration.

Zoom In Zoom Out Reset image size

In Figure 14, we present the dependence of the dynamical tides model on f_dyn. Larger f_dyn corresponds to the more efficient extraction of energy via tidal force, which yields faster tidal migration. In terms of the population, the period distribution peak moves toward larger periods as the efficiency of the tides rises, and there is enhanced tidal disruption, together with enhanced production of HJs and WJs.

Zoom In Zoom Out Reset image size

In Figure 15, we present the dependence on the choice of the boundaries in the semimajor axis distribution, and consider the contributions from different semimajor axes. The overall distribution suggests a preference for HJ production from large initial semimajor axes, but the results seem robust under the choices of distributions.

Zoom In Zoom Out Reset image size

Key findings. As discussed and shown above, inflated eccentric migration following planet–planet scattering in planetary systems could potentially explain the whole population of eccentric WJs and a significant fraction of HJs (or even all of them, given the lowest inferred estimates). It also leads to a high disruption rate of systems, as it accelerates the parameter-space flow from HJs to tidally disrupted gas giants. Inflated eccentric migration could also play an important role in other systems, where high eccentricities are excited through secular processes in general and ZLK oscillations in particular—these will be discussed in future studies (see also Petrovich 2015b for a study of contracting planets in this context). Our results suggest that any modeling of planetary systems, and in particular of young planetary systems, should self-consistently account for thermal evolution. It should consider the evolving size of the planets and the coupling between the thermal and dynamical evolutions. These could play key roles in the system's dynamics, and in the final sculpting of its architecture.

The longer overall timescales and shorter time spent at eccentric orbits would potentially allow for more significant contraction of the gas giants before significant migration occurs, giving rise to a weaker, though still important, effect of inflated eccentric migration. It should be noted that it is important to model the ZLK–Octupole order, in this context, given the conditions and timescales involved (Naoz et al. 2011).

Inflated eccentric migration is most pronounced at the earliest times, when planets are still in their infancy/at a young age and are still far more inflated than at later times, after contraction. Moreover, when eccentricity excitation occurs through secular processes, tidal effects can induce tidal precession, which can quench the eccentricity excitation, giving rise to lower eccentricities that would be expected without the effects of tides, thereby leading to larger pericenter approaches and less effective tidal dissipation and migration.

Implications for giant-planet formation. Our results could have a wide range of implications with respect to various aspects of giant-planet formation and evolution. They shed light on and point to the critical role played by the physical evolution of the planets, and its coupling with the dynamical evolution, which significantly changes the behavior of eccentric migration processes. It gives rise to the efficient formation of both HJs and eccentric WJs, where the latter, in particular, are more difficult to form efficiently through previously studied eccentric migration. We provide their relative fractions overall, and as a function of the age of the systems, as well as detailed predictions of their physical properties. Furthermore, the shorter migration timescales due to inflated eccentric migration could give rise to very young HJs, Myrs old, which are typically suggested to form through disk migration (see the detailed review of disk migration in Baruteau et al. 2014) or via other channels that focus on the stages after gas dissipation (e.g., Wu et al. 2007). The HJs and WJs formed in the proposed channel could have a range of inclinations, even retrograde ones, but would generally have a preference for prograde orbits, given that their initial inclinations were excited by planet–planet scattering (Beaugé & Nesvorný 2012).

Potential caveats and challenges. In our models, all the WJs formed via eccentric tidal migration are effectively transient HJs that did not complete their migrations within a given time, i.e., eccentric WJs rather than circular ones (although they could reach relatively small eccentricities, and even ≲ 0.1), as can be seen in the eccentricity distributions in Figures 3 and 13. Even with the inflated radii, allowing for the formation of lower-eccentricity WJs, low-eccentricity (<0.6) WJs can hardly be formed through inflated migration, and were likely formed in situ and/or through disk migration, generally consistent with the analysis by Anderson et al. (2020). We note that the paucity of >0.9 eccentricity suggested by Socrates et al. (2012), and ruled out by Dawson et al. (2015), does not constrain our models, which indeed show >0.9 WJs to be very rare. That being said, the apparent low frequency of 0.6–0.9 eccentricity WJs is a challenge to inflated eccentric migration, and any other eccentric migration model. In fact, this is a potentially more general difficulty—any successful model of HJ/WJ production should also be able to suppress HJ/WJ formation through the various types of eccentric migration.

On the theoretical front, neither the tidal interactions of gas giants nor the heat transfer to the core are understood, giving rise to large uncertainties in the evolutionary models. Here, we tried to bracket these potential caveats, but, naturally, better understanding of these processes is critical for the assessment and modeling of any eccentric migration model.

Future work on eccentric migration. Though we have focused on the role of inflated migration for specific types of eccentric migration (initially by planet–planet scattering, and considering some specific models for weak and dynamical tides), the same coupled evolution is important for any suggested eccentric migration model. Follow-up papers may consider other such models and their variants. In addition, accounting for the inflated sizes of gas giants in their early phases is also important for the increased likelihood of their physical collisions with other planets (due to their larger cross sections), which, in turn, can give rise to the heating and further inflation of the collision product (Lin & Ida 1997), which could then affect the migration and further collisions. Eccentric inflated migration could also be coupled to other processes that play a role in planet formation and dynamics, such as photoevaporation (e.g., Tripathi et al. 2015). Since the geometric cross sections of initially inflated planets are larger than the cross sections of the noninflated ones, the role of photoevaporation might change accordingly. Furthermore, although we have discussed the eccentric tidal migration channel in this paper, initially inflated gas giants could also be discussed in the context of disk migration, where the radius plays a role in the evolution, too.

Central heating. We can use our model to set constraints on the effective amount of energy penetrating to the center, since it will affect the contraction timescale and hence the migration timescale. The exact amount of the energy deposit and its distribution are still unknown, but by using our population synthesis—assuming a given distribution, e.g., for simplicity, that all the energy is deposited at the center—we can estimate the energy amount needed to explain the observational results. It should be noted that the parameter space is large and includes many degeneracies. For example, the choice of the initial eccentricity distribution might shift the distribution in a similar direction as external heating would do.

In this paper, we have proposed the formation channel of hot and warm Jupiters via inflated eccentric migration, and discussed the implications on the population. Here, we have focused on the semi-analytical approach, and discussed specific examples and a detailed population synthesis. We have compared the semi-analytical approach with the numerical planet evolution models described in paper II, and found that the semi-analytical model is in good agreement with the numerical model in all the regimes available for the numerical model. This allows us to make use of the more efficient semi-analytical model to explore the dynamics of large populations of planets as well as the resulting populations of HJs/WJs for a wide range of initial conditions.

Our models are general and are able to include in principle any kind of external heating/dynamical evolution. For brevity, we have demonstrated the use of the semi-analytical model for several examples only. We have presented specific examples for inflated eccentric migration and provided a detailed population synthesis study, based on the semi-analytical model. We have studied the dependence of the resulting HJ/WJ populations on the assumptions made regarding the tidal models and heat transfer processes, as well as the initial properties of the progenitor planets and their initial orbits.

Using the detailed semi-analytical model and population synthesis results, we have showed that inflated eccentric migration could significantly shorten the migration timescales of HJs and WJs, generally form WJs more efficiently, and give rise to enhanced rates of tidal disruptions of gas giants. We have also considered the effect of the external energy injected into the migrating gas giants on their final formed population, and found that it leads to enhanced tidal disruption, together with a smaller fraction of HJs and a larger fraction of WJs compared with models without efficient central heating.

Inflated eccentric migration leads to significant differences in the final distribution of the parameters as compared with noninflated models, and suggests that inflated migration plays an important role in migration and should generally be accounted for in any eccentric migration models (and possibly also in disk migration models).

We gratefully acknowledge fruitful discussions with Sivan Ginzburg, Thaddeus D. Komacek, Michelle Vick, Nicholas C. Stone, and Eden Saig. M.R. acknowledges the generous support of the Azrieli fellowship.