Formation of hot Jupiters through disk migration and evolving stellar tides

Here we address the hot Jupiter (hJ) pile-up at 0.05 AU around young solar-type stars observed in stellar radial velocity surveys, the hJ longterm orbital stability in the presence of stellar tides, and the hJ occurrence rate of 1.2 (+-0.38)% in one framework. We calculate the combined torques on the planet from the stellar dynamical tide and from the protoplanetary disk in the type II migration regime. We model a 2D non-isothermal viscous disk parameterized to reproduce the minimum-mass solar nebula and simulate an inner disk cavity at various radial positions near the star. We choose stellar rotation periods according to observations of young star clusters. The planet is on a circular orbit in the disk midplane and in the star's equatorial plane. We show that the torques can add up to zero beyond the corotation radius around young stars and stop inward migration. Monte Carlo simulations predict hot Jupiter survival rates between ~3% (alpha disk viscosity of 1e-1) and 15% (alpha = 1e-3). Once the protoplanetary disk has been fully accreted, the surviving hJs are pushed outward from their tidal migration barrier and pile up near 0.05 AU, as we demonstrate using a numerical implementation of a stellar dynamical tide model. Orbital decay is negligible on a billion year time scale due to the contraction of the highly dissipative convective envelopes in young Sun-like stars. We find that the higher pile-up efficiency around metal-rich stars can at least partly explain the observed positive correlation between stellar metallicity and hJ occurrence. Combined with the observed hJ occurrence rate, our results for the survival rate imply that<8 % (alpha = 1e-3) to<43 % (alpha = 1e-1) of sun-like stars initially encounter an inward migrating hJ. This reconciles models and observations of young spinning stars with the observed hJ pile up and hJ occurrence rates.


Introduction
Soon after the surprising detection of Jupiter-mass planets in very close orbits around Sun-like stars (Mayor & Queloz 1995), it was proposed that these hot Jupiters cannot have formed in situ but that they must have migrated from the cold, icy regions of the protoplanetary disk at several AU from the star (Lin et al. 1996). Competing theories have been put forward as to what stops the inward migration of planets (for a recent review see Dawson & Johnson 2018): tidal halting (Trilling et al. 1998), magnetorotational instabilities that evacuate the close-in protoplanetary disk (Kuchner & Lecar 2002;Romanova & Lovelace 2006), planet-disk magnetic interactions (Terquem 2003), the Kozai mechanism of a distant perturber (Nagasawa et al. 2008), planet traps (Hasegawa & Pudritz 2010), planet-planet scattering (Naoz et al. 2011;Wu & Lithwick 2011), or high-eccentricity migration (Wang et al. 2017).
Although the tidal stopping mechanism offers the best agreement with observations (Plavchan & Bilinski 2013), none of the previous theories for hot Jupiter formation could explain the following observations at the same time: (1) the sharp pile-up of hot Jupiters at 0.05 AU around Sun-like stars observed in radial velocity (RV) surveys; (2) the increase of the hot Jupiter occurrence rate with stellar metallicity in the Kepler planet sample (Petigura et al. 2018); and (3) the longterm orbital stability of hot Jupiters under the effect of tidal dissipation in the star. These are the points that we address in this study.
At 0.1 AU, tidal dissipation in the star is sufficiently large to affect the planetary orbit. Planets on circular orbits with an orbital plane near the stellar equatorial plane and with orbital semi-major axes (a) larger than the stellar corotation radius (r co ) are repelled, whereas planets interior to r co are driven into an ever faster orbital decay until they are either tidally disrupted or they fall into the star. Most previous studies used equilibrium tide models with an assumed fixed tidal dissipation constant (Q , typically chosen between 10 5 and 10 6 ) (Lin et al. 1996;Trilling et al. 1998;Dobbs-Dixon et al. 2004;Rice et al. 2012) to parameterize tidal dissipation in the star and to evaluate the tidally driven orbital circularization and migration of close-in planets. The resulting tidal torque is too weak to stop a migrating Jupitermass planet.
Moreover, constant-Q (or rather constant angle or constant phase lag; Greenberg 2009) models with a host star that has a Sun-like rotation period predict a gradual infall of hot Jupiters into their stars on a billion year time scale. Although the equilibrium tide model is compatible with the low frequency of planets within 0.03 AU around Sun-like stars, it requires a delicate finetuning of the constant stellar dissipation factor or of the initial conditions in the protoplanetary disk (Rice et al. 2012)  the existence and even pile-up hundreds of known hot Jupiters at about 0.05 AU around Sun-like stars.
The efficiency of tidal dissipation in the star is determined by the presence and extent of the convective envelope of the star (Zahn 1977;Ogilvie & Lin 2007). While Sun-like stars on the main-sequence have their core-envelope boundary at about 0.7 solar radii (R ), pre-main-sequence stars are much larger than our Sun today and they can have much more extended envelopes. As a consequence, while solar-type stars on the main sequence respond to the tidal perturbation by close-in massive planets with a tidal dissipation function of 10 5 Q 10 8 , young stars are much more dissipative (Bolmont & Mathis 2016) with Q as low as ≈ 10 3 , depending on the exact structure of a given star. The source of this dissipation is in the so-called dynamical tide within the star's convective regions, which is due to inertial waves that are caused by the Coriolis acceleration (Ogilvie & Lin 2007). Inertial waves are driven as long as the modulus of the planet's orbital mean motion |n| < 2|Ω |, where Ω is the stellar spin rate. As the tidal dissipation, and therefore the tidally induced orbital decay of the planet, depends on both n and Ω , a consistent picture for the tidal migration of hot Jupiters requires a model of the stellar spin evolution and its effects on the transfer of rotation to orbital angular momentum. Toward the end of a star's lifetime, stellar mass loss and the engulfment of the planet within the extended gaseous envelope of the star can further affect a planet's orbital evolution under the effect of the dynamical tide (Rao et al. 2018).
With n and Ω affecting the tidal dissipation regime, the star's initial rotational spin-up and subsequent magnetic braking are key ingredients to hot Jupiter formation and evolution. Stassun et al. (1999) found that solar mass stars in the ∼1 Myr old Orion Nebula have typical rotation periods (P rot ) between about 0.5 d and 8 d (the upper threshold being uncertain due to observational biases), implying corotation radii between about 0.01 AU and 0.08 AU. The rotation periods of stars with ages between 1 Myr and 5 Myr cluster between breakup periods (∼0.5 d) and about 8 d with a long tail in the period distribution up to ∼30 d (Irwin et al. 2008). We use this information to parameterize our disk model and the resulting torques on the planet (Sect. 2.1).
Figure 1(a) shows the well-known pile-up in the semimajor axis distribution of all known exoplanets with masses > 0.1 Jupiter masses (M J ) around sun-like stars. The open and closed circles refer to hot Jupiters with known stellar spin-orbit obliquities (ψ), which illustrates that the pile up is present in both lowand high-obliquity hot Jupiters. 1 Appendix A shows the same data plotted along a linear abscissa. In this paper, we develop a theory of hot Jupiter formation that can, at least partly, explain the existence of a pile-up in the sample of Kepler planets with RV measurements.

Methods
In the early phase of planet formation, giant planets supposedly form beyond the circumstellar ice line at a few AU around Sunlike stars (Hayashi 1981). They then migrate to close-in orbits at about 0.1 AU or less within the proto-planetary disk. Before a protoplanet has accreted sufficient mass to open up a gap in the disk, its radial drift is referred to as type I migration and it is driven by the Lindblad torque (Γ LB ) and, as the case may be, the corotation torque (Goldreich & Tremaine 1979;Lin & Papaloizou 1986). Planets with masses similar to or larger than that of Jupiter, however, open up a gap in the disk, which then leads to type II migration on the viscous time scale of the disk (Ward 1997;Nelson et al. 2000) 2 . As we are interested in hot Jupiters in this study, we consider disk torques on the planet in the type II migration regime (Γ II ).
If the tidal dissipation in the star is strong enough, the inward migration of a planet may halt at a stellar distance where the torque on the planet exerted by stellar tide (Γ t ) compensates for Γ II , i.e. where Γ t + Γ II = 0. We refer to this distance as the tidal migration barrier.

Disk model
The disk torque on the planet depends on the local disk properties. We assume that the planet orbits the star in the disk midplane, which has a temperature T m . We model a rotationally symmetric, two-dimensional optically gray disk with a vertical temperature gradient determined by the disk's viscous heating and by the stellar irradiation. The disk effective temperature is given by (Hubeny 1990) where T i is the temperature due to stellar illumination, τ = κΣ p /2 with κ ext as the Rosseland mean extinction opacity, κ abs as the Rosseland mean absorption opacity, and Σ p as the disk gas surface density at the position of the planet. The use of mean Rosseland mean extinction opacities is justified because the Planck mean and Rosseland mean extinction opacities are comparable in fully mixed dusty disks (Pollack et al. 1994). Furthermore, following Pollack et al. (1994), we set κ abs = κ = κ ext because it has been shown that this is adequate for disk temperatures and optical depths around accreting stars (Menou & Goodman 2004). Hence, τ abs = τ = τ ext . We consider a two-faced disk in thermodynamic equilibrium so that its cooling rate Q − = 2σ SB T 4 eff,d is equal to its heating rate Q + = 9 4 νΣ p Ω 2 (Menou & Goodman 2004). We assume that viscous heating is by far dominant within 0.2 AU to the star (T i T m ), the regime we are interested in. We make use of the Shakura & Sunyaev (1973) relation that describes the disk viscosity as ν = αc 2 s /n with α as the disk's kinematic viscous efficiency parameter and n = (G[M p + M ]/a 3 ) 1/2 as the Keplerian orbital frequency at the orbital radius a. With the speed of sound given as c s = k B T/µ, where k B is the Boltzmann constant and µ is the mean molecular weight of the gas in units of the proton mass (m p + ), we transform Eq. (1) into The mean molecular mass of the gas is determined by the degree of ionization, which can be derived from the Saha equations. For the range of disk temperatures we are interested in (1 000 K T 6 000 K), the Saha equations predict 1.3 ≤ µ ≤ 2.4 (see Eq. 19 and Fig. 2 of D'Angelo & Bodenheimer 2013) in a disk with a composition similar to the protosolar nebula, that is, with hydrogen and helium mass fractions of X = 0.7 and Y = 0.28, respectively. We use µ = 1.85 in our nominal disk model.
The mass accretion rates of T Tauri stars, the variations of FU Orionis outbursts, dwarf nova, and X-ray transients suggest 10 −3 α 10 −1 (King et al. 2007). The lower end of this scale likely represents partially ionized disks while the upper end of this range can be reached in fully ionized disks (Martin et al. 2019). The accretion rates and total masses of protoplanetary disks suggest 10 −3 < α < 10 −2 beyond 10 AU from solar type stars. In these regions, the disks are only partly ionized (e.g. the outermost layers). The inner parts of a protoplanetary disk (< 0.1 AU), however, can be highly thermally ionized, and so α > 10 −2 is fairly reasonable. There is also strong evidence that α must be < 1 (Martin et al. 2019) and the substantial ionization of the disk in the hot Jupiter regime strongly suggests α > 10 −3 . All things combined, we chose to test three reference values of α ∈ {10 −3 , 10 −2 , 10 −1 } representing a range of possible disk viscous efficiencies, irrespective of whether they be dominated by eddies in the turbulent gaseous disk or by its magnetic buoyancy (Shakura & Sunyaev 1973). In our nominal disk model we use α = 10 −3 .
(2) for our nominal disk model. Three relations are shown for different disk opacities, κ ∈ {10 −5 , 10 −6 , 10 −7 } m 2 kg −1 . In the following, we use κ = 10 −7 m 2 kg −1 because (i) this curve reproduces the 2 000 K at 0.05 AU predicted by Lin et al. (1996); and (ii) it is in good agreement with the midplane temperatures predicted in a numerical 2D rotationally symmetric model by Bell et al. (1997) of a viscous disk for α ≈ 0.001 and a stellar accretion rate of 10 −7 M yr −1 , where M is the solar mass. We need to keep in mind, however, that in a more realistic scenario the dust opacity could be modelled as a function of the temperature itself, κ = κ(T ) (Henning & Stognienko 1996).
The stellar rotation determines the corotation radius and it can, in some scenarios, influence the location of the magnetic cavity. In Fig. 3 we show the spin period distribution of stars with masses between 0.75 M and 1.25 M observed in the 5 Myr young open cluster NGC 2362 by Irwin et al. (2008). The gray bars show the original data and the solid black line shows our fit of a squared normal distribution (N(µ, σ)) 2 to the data with µ = 8.3 d as the mean and σ = 5 d as the standard deviation of N. This rotational period distribution defines a probability density function that we use to randomize the stellar spin period.
Finally, we add the effects of a magnetic inner cavity in the disk, which is reflected in our model as a truncation of the disk torque at the inner 2:1 orbital mean motion resonance with the radial position of the magnetic cavity. The orbital radius at which the disk torque is truncated is located at 0.63 r mc . 3 We explore two scenarios for the location of the magnetic cavity. In the first scenario, the magnetic cavity is drawn from a normal distribution parameterized as r mc = 0.05 (±0.015) AU. This scenario is similar to the one proposed by Kuchner & Lecar (2002), in which the gas disk is truncated at a temperature of 1500 K by the onset of a magnetorotational instability. Kuchner & Lecar (2002) estimate that this temperature is reached at a distance of about 0.067 AU. Similarly, Romanova & Lovelace (2006) argued that the critical distance for a cavity to form is the Alfvén radius (r A ), which is at about 0.05 AU for solar-mass T Tauri stars. In a second scenario, we consider that the magnetic cavity is at the same radial distance as the corotation radius.  Irwin et al. (2008), the solid line is our analytical fit to the data using a squared normal distribution with µ = 8.3 d and σ = 5 d.

The tidal torque
The star's tidal torque on a nearby planet in the star's equatorial plane can be estimated as (Efroimsky & Makarov 2013 where k 2 is the star's 2nd degree tidal Love number and g is the instantaneous angular separation between the line connecting the stellar center with the planet and the line connecting the stellar center with the center of the star's tidal bulges. If we assume that g is frequency independent, that the orbital eccentricity is small, and that the planetary orbit is aligned with the stellar equatorial plane, then the tidal torque follows from the quadrupolar modes of the tidal potential and we can introduce a stellar dissipation factor Q as per sin(2 g ) ≈ 1/Q (Murray & Dermott 1999). We can then derive a frequency-averaged, dimensionless quantity for the stellar tidal dissipation as per k 2, /Q = D ω (Bolmont & Mathis 2016). Thus, We use the estimates for the dissipation by the stellar dynamical tide (Bolmont & Mathis 2016;Amard et al. 2016;Bolmont et al. 2017) to calculate Γ t (a) on a close-in planet. These stellar models consider frequency-averaged tidal dissipation in the star's convective envelope (Ogilvie 2013), which is dominated by the dynamical tide for planets with |n| < 2|Ω | and which is dominated by the stellar equilibrium tide otherwise. Variations of the stellar tidal dissipation over short frequency intervals are thus mitigated into the averaged, or effective, tidal dissipation factor Q (for details see Sect. 2.3). In our nominal parameterization of a star-planet-disk system, the stellar frequency-averaged tidal dissipation is D ω = 10 −3.25 , a typical value for a solar-type star during the first ≈ 10 Myr of its lifetime (Bolmont et al. 2017), which corresponds to a frequency-averaged tidal dissipation factor of about 10 3.4 . We consider a nominal stellar corotation radius at 0.02 AU (corresponding to a stellar rotation period of 1 d), a stellar radius of 2 R , and a Jupiter-mass planet. Beyond this nominal parameterization of the system, which is primarily meant to have an illustrative purpose, we also explore other corotation radii corresponding to the empirically determined rotation period distribution of young stars as per (Irwin et al. 2008).
The algebraic sign of Γ t is positive beyond the stellar corotation radius (a > r co ), where the star transfers its rotational angular momentum to the planet's orbital angular momentum. In turn, planets closer to the star than the corotation radius transfer angular orbital momentum to spin up the star and, hence, Γ t (a < r co ) < 0.

The disk torque
With L = M p Ga(M p + M ) as the planet's orbital angular momentum and without the effects of mass accretion onto the planet (dM p /dt ≡Ṁ p = 0; t being time), the disk torque in the type II migration regime is given as where (Ida & Lin 2004;Alibert et al. 2005) Note that in Eq. (6), the disk viscosity depends on the distance to the star, in particular through its coupling with the sound velocity and the midplane temperature in our disk model (Sect. 2.1).

Stellar evolution models
We also consider the longterm orbital evolution during the tidally dominated period, i.e. once the proto-planetary gaseous disk has gone. We use precomputed stellar evolution tracks (Lagarde et al. 2012;Bolmont et al. 2017) that were generated with the STAREVOL code (Amard et al. 2016) 4 . These models consider a 1D rotating star with the radiative core rotating at a different speed than the convective envelope, and they include centrifugal accelerations as well as the resulting chemical stratification. The initial spin period was set to 1.4 d and the incremental stellar angular momentum loss during each numerical integration time step was calculated according to a differential equation of the torque exerted by magnetic braking with the stellar wind (Bouvier et al. 1997). The effect of a nearby planet is not taken into account in these models, but according to Bolmont et al. (2017) the variation due to a hot Jupiter would be limited to about 1 day in the stellar rotation period after 5 Gyr. For the critical phase of the tidally-driven hot Jupiter pile-up, which happens on a time scale of 10 Myr, the tidal effects on the stellar spin can thus safely be neglected. The situation is somewhat different for stars that merge with a Jupiter-mass planet after its tidally driven infall (Bolmont et al. 2012), but we ignore this effect of star-planet mergers and rather focus on the hot Jupiter pile-up. The frequency-averaged tidal dissipation D ω is calculated during the stellar evolution assuming a simplified two-layer model of the star, i.e. a radiative core and a convective envelope. The analytical description developed by Ogilvie (2013) takes into account the dominating tidal frequencies as a function of the (evolving) stellar properties and is computationally very efficient (Mathis 2015). With D ω = ∞ −∞ dω Im[k 2 2 (ω)]/ω, this procedure is equivalent to calculating the imaginary part of the star's second degree tidal Love number along the stellar evolution track.
In the stellar evolution models, the star contracts and spins up for the first about 100 Myr until it reaches a minimum rotation of 0.25 d. In this early phase, the stellar corotation radius moves inward (Bolmont et al. 2012). Most important for our purpose, the quality factor is calculated consistently from the stellar interior evolution, i.e. from the extent of its radiative core and of its convective envelope. This is an important improvement to earlier attempts, which used fixed, nominal Q values for particular stages of stellar evolution (see, e.g. Trilling et al. 1998).

Orbital evolution due to tides
For our numerical simulations of the tidally driven orbital evolution, we use the precomputed, time-dependent tidal dissipation functions from the stellar evolution tracks and compute the incremental tidal evolution of the planet's orbital semimajor axis (a) and orbital eccentricity (e) with an adaptive time step dt = 10 −4 t as (Eggleton et al. 1998) where (Hut 1981;Hansen 2010) and and with G as the gravitational constant. Equations (7) and (8) are valid for arbitrary eccentricities (Leconte et al. 2010) and we tested systems with small initial eccentricities that did not result in qualitatively different behavior of the orbital evolution. Hence, we focus this report on e = 0 and assume that the stellar and planetary spin axes are aligned.
We read D ω from the precomputed stellar evolution models based on the frequency-averaged analytical expressions (Ogilvie 2013;Mathis 2015). The frequency-averaged tidal dissipation constant of the star is calculated asQ = 3/(2 D ω ). Equations (11) and (12) ensure that the planet excites tidal inertial waves in the stellar convective layer as long as n < |2Ω |, whereas tidal friction is more adequately modelled by the equilibrium tide for more close-in planets with n ≥ |2Ω | (Bolmont & Mathis 2016). Converting orbital and spin frequencies into orbital radii, Kepler's third law of planetary motion predicts the transition radius (r e↔d ) between the equilibrium tide and the dynamical tide regime as The calibrated tidal dissipation constants in equations (14) and (15) are taken from Hansen (2010Hansen ( , 2012. In order to compare this model to a pure equilibrium tide model with fixed Q and constant stellar rotation, we set up another suite of simulations, which assumes a sun-like stellar radius, mass, and rotation and estimates σ in Eq. (10) as an approximation that is only valid in the limit of small eccentricities (Bolmont & Mathis 2016). Figure 4 displays the total torque Γ tot ≡ Γ t (a) + Γ II (a) (black solid line) acting on a close-in Jupiter-mass planet in our nominal disk scenario and for a Sun-like star with its corotation radius at 0.02 AU. The purpose of Fig. 4 is to illustrate the general behavior of the torques for different possible radial arrangements of the disk's inner magnetic cavity (here fixed at 0.05 AU) and the stellar corotation radius. The 2:1 MMR with the cavity determines the truncation of the disk torque acting on the planet. The corotation radius determines the transition from the dynamical tide to the equilibrium tide regime. All things combined determine the existence and, as the case may be, the radial position of the tidal migration barrier. In Fig. 4(a), P rot = 1 d and the  Fig. 4. Torques exerted on a Jupiter-mass planet that is embedded in a proto-planetary disk around a young, solar-type star. (a) Nominal parameterization of the star-planet-disk model, with R = 2 R , P rot = 1 d, D ω = 10 −3.25 , Σ p,0 = 1 000 g cm −2 , α = 10 −3 , µ = 1.85, κ = 10 −7 m 2 kg −1 , and r mc = 0.05 AU. The blue dashed curve refers to the tidal torque (Γ t ), the red dotted line to the disk's type II migration torque (Γ II ), and the thick solid line to the total torque. The location of zero total torque at 0.031 AU is indicated with a black circle. Γ t , and therefore also Γ t + Γ II , swap signs at the stellar corotation radius, interior to which any planet would be rapidly pulled into the star. (b) Similar to (a) but with P rot = 8 d. corotation radius is interior to the magnetic cavity. In Fig. 4(b), P rot = 8 d and the inner cavity happens to be at 0.63 times the corotation radius, that is, in a 2:1 MMR with r co . In Fig. 4(a), the negative type II migration torque (red dotted line) dominates the positive torque from the stellar dynamical tide (blue dashed line) beyond the 2:1 MMR with the inner cavity, so that Γ tot < 0 and the planet migrates inward. Interior to the 2:1 MMR with the cavity, the disk torque is zero and so the total torque (black line) switches its algebraic sign according to the positive (i.e. repulsive) torque of the stellar dynamical tide. As a consequence, the total torque must go through a point of zero torque. This location of zero torque, which we refer to as the tidal migration barrier, is indicated with a black dot and it is the distance at which planet migration stops.

Tidal migration barrier at zero torque
Interior to the tidal migration barrier, where the planet cannot arrive via migration alone, the strong distance dependence of Γ t (a) ∝ a −6 leads to a very effective tidal repulsion. At the corota-tion radius then, the algebraic sign of the tidal torque (and therefore of the total torque as well) switches negative and the planet would fall towards the star relatively fast. At the 2:1 MMR with the corotation radius, however, the stellar tides transition from the dynamical tide to the equilibrium tide, which is much weaker than the dynamical tide. As a consequence, the magnitude of the negative tidal torque is reduced as well. This transition is located at at 0.02 AU × 0.63 = 0.0126 AU (see Eq. 16) in Fig. 4(a). It is well possible that hot Jupiters, which do not encounter a tidal migration barrier beyond r co may survive in this extremely close stellar vicinity on a time scale of 10 Myr -100 Myr if they can also withstand tidal disruption.
In Fig. 4(b), where P rot = 8 d, the corotation radius is beyond the inner magnetic cavity of the disk, which is at 0.05 AU as in Fig. 4(a). The stellar rotation period was intentionally chosen to be 8 d so that the transition between the dynamical tide and the equilibrium tide regime coincides with the radial position of the inner cavity, although it is entirely unclear at this point if the magnetic cavity can be affected directly by the stellar rotation at all.
Beyond our nominal parameterization of the star-planetdisk system in Fig. 4, we explore a range of possible realizations in Figs. 5-7. We keep the stellar radius, and the planetary mass fixed, while we draw D ω from a log-normal distribution as per log 10 ( D ω ) = −3.25 ± 0.5. Similarly, we vary log 10 (Σ p,0 /[g cm −2 ]) = 3 ± 1 and draw the mean molecular weights from a normal distribution as per µ = 1.85 ± 0.55. We also randomize P rot as per the probability density distribution shown in Fig. 3 and also randomly draw the location of the magnetic cavity from a normal distribution parameterized as r mc = 0.05 (±0.015) AU (see Sect. 2.1).
We generate 10,000 Monte Carlo simulations for α = 10 −1 (Fig. 5), α = 10 −2 (Fig. 6), and α = 10 −3 (Fig. 7), respectively. These randomized parameterizations are carried out to derive a plausible distribution of the tidal migration barrier for hot Jupiters around sun-like stars. The resulting rates of a tidal migration barrier outside of the corotation radius are 2.9 % for α = 10 −1 , 8.1 % for α = 10 −2 , and 15.6 % for α = 10 −3 . Panels (a) in Figs. 5-7 show the distribution of the corresponding tidal, disk, and total torques in analogy to Fig. 4. Panels (b) in Figs. 5-7 illustrate the distributions of the magnetic cavity (solid lines), tidal migration barrier (dashed lines), and corotation radius (dotted line). In comparing Figs. 5-7, note that the histograms of r mc and r co do not depend on the respective choice of α, while the distribution of r tmb does depend on α. The only reason for a magnetic cavity not to exist in our simulations is for it to be located within the stellar radius.
As an alternative scenario for the location of the magnetic cavity, we consider that the magnetic cavity coincides with the corotation radius. In Fig. 8 we plot the histograms of r mc , r tmb , and r co from the resulting Monte Carlo simulations. We find that the rates for the tidal migration barrier being located beyond the corotation radius are very similar in all three cases of the α viscosity parameter that we studied as summarized in the upper left corners of the panels in Fig. 8.
With f occ = 1.2 ± 0.38 % as the bias-corrected observed occurrence rate (Wright et al. 2012), f sur ∈ {2.9 %, 8.1 %, 15.6 %} as the survival rates from our simulations, f alt as the rate of all alternative mechanisms to stop a migrating Jupiter from falling into its star, and f for as the hot Jupiter formation frequency around sun-like stars we have f occ = ( f sur + f alt ) · f for . This is equivalent to f for = f occ /( f sur + f alt ) < f occ / f sur , for which we can estimate f sur from our simulations.
In the first scenario, where the radial position of the magnetic cavity is chosen from a normal distribution, we have f for < 41.4 ± 13.1 % for α = 10 −1 , f for = 14.8 ± 4.7 % for α = 10 −2 , and f for = 7.7 ± 2.4 % for α = 10 −3 . In the second scenario, where the magnetic cavity coincides with the corotation radius, we find f for < 42.9 ± 13.6 % for α = 10 −1 , f for = 14.8 ± 4.7 % for α = 10 −2 , and f for = 8.1 ± 2.6 % for α = 10 −3 . Hence, in both scenarios the resulting hot Jupiter survival rates and the formation rates are very similar. In other words, when combined with observations, our simulations suggest that less than about 8 % (α = 10 −3 ) to 43 % (α = 10 −1 ) of sun-like stars initially encounter an inward migrating hot Jupiter. These values are almost irrespective of whether the magnetic cavity is located at about 0.05 AU and independent from the stellar rotation or if the magnetic cavity is coupled to the corotation radius. All these results are summarized in Table 1.

Hot Jupiter pile up from tidally driven migration
Once the protoplanetary disk has been accreted onto the central star after about 10 Myr into the star's lifetime (Haisch et al. 2001), the disk torque vanishes and the orbit of a hot Jupiter evolves under the effect of stellar tides only, neglecting the possibility of interaction with other planets or nearby stars. Figure 9 shows the outcome of our numerical simulations for a solar metallicity star, which are based on differential equations for da/dt (assuming e = 0) as derived from the orbit-averaged torque for a tidally evolving two-body system. Panels (a)-(c) on the left illustrate the first 10 Myr of evolution, and panels (d)-(f) on the right side show the evolution over 4.5 billion years, that is, over the age of the solar system. Figures 9(a) and (d) demonstrate the variation ofQ over several orders of magnitude. The ini-tialQ ≈ 10 3.25 imply highly effective dissipation over some 10 Myr, whereas values ofQ > 10 8 after a billion years mean negligible dissipation. The initial spin-up due to contraction and the subsequent spin-down owing to magnetic braking are displayed in panels (b) and (e). The orbital evolution of 100 starplanet two-body systems is shown in panels (c) and (f).
In Fig. 9(c), colored lines refer to the planet's orbital evolution, with color encoding the initial orbital semi-major axis. Three mechanisms are readily visible: (1) the pile-up of hot Jupiters at about 0.05 AU after just about 10 Myr; (2) the rapid infall of planets interior to the stellar corotation radius (initially at 0.025 AU); and (3) the switch from the dynamical tide to the equilibrium tide regime at orbital frequencies n = |2Ω | (initially at about 0.015 AU). In panel (f), we compare a subset of these orbital evolution tracks to another set of tracks that we calculated using the conventional constant-Q model and assuming a constant stellar rotation period of 27 d. In this model, the corotation radius is fixed at about 0.18 AU and any planet interior to this relatively wide orbit will permanently fall into the star. Within 4.5 billion years, any hot Jupiter that started at 0.045 AU to the star or closer is destructed. A pile-up, however, is not reproduced.
This discrepancy becomes even more apparent in the histograms shown in Figs. 10(a) and (b). Here we show snapshots of the simulated hot Jupiter populations in the pure equilibrium tide model (a) and in the dynamical tide model with stellar evolution (b) at 10 Myr (empty bars), 100 Myr (striped bars), and 1 000 Myr (gray bars), respectively. While the equilibrium tide model suggests a steady removal of close-in planets over a billion years, the dynamical tide and stellar evolution model predicts that the hot Jupiter population is essentially formed after between 10 Myr and 100 Myr. Moreover, while the equilibrium tide model does not produce a pile-up, such a peak in the planet distribution occurs very naturally in the dynamical tide model with stellar evolution.
In Fig. 10(c)    tidal torque on their close-in massive planets. More explicitly, we find that the hot Jupiter pile up is much more pronounced about metal-rich stars, increasing by about a factor of two as metallicity is increased from Z = 0.004 (or [Fe/H] = −0.53, gray shaded histogram) to Z = 0.02 (or [Fe/H] = +0.28, empty histogram). Even within the early years of exoplanet observations it became ever more apparent that stellar metallicity affects the likelihood of a star harbouring a planet (Gonzalez 1997). It has then been debated whether the magnitude of this correlation may be different for different planetary masses, with larger planets possibly following a stronger correlation (Buchhave et al. 2014;Schlaufman 2015). Spectroscopy of Kepler exoplanet host stars delivered further evidence of a positive relation between stellar metallicity and hot Jupiter occurrence rate (Petigura et al. 2018). Ida & Lin (2004) argued that such a trend could originate in the protoplanetary disks since more metal-rich disks forming metal-rich stars should also have had more solids available to form planets. Alternatively, recent simulations showed that stellar metallicity can also affect the tidally driven migration of close-in planets (Bolmont et al. 2017). Here we propose that the stronger tidal migration barrier around metal-rich stars shown in Fig. 10(c) can also explain the positive correlation between stellar metallicity and hot Jupiter occurrence rate. The root of the Article number, page 9 of 14 A&A proofs: manuscript no. ms correlation would then not be in the initial formation rate of hot Jupiters but in their survival rate against destruction via orbital migration toward to the star.

Discussion
The very existence of a zero total torque location (the tidal migration barrier) is an important feature of our model and rooted in the tidal dissipation implied by the stellar evolution tracks (Amard et al. 2016;Bolmont et al. 2017), which is much higher and the resulting tidal dissipation factorQ much smaller than previously assumed. For comparison, using canonical tidal dissipation factors of 10 5 ≤Q ≤ 10 6 (Trilling et al. 1998;Trilling 2000;Pätzold & Rauer 2002;Trilling et al. 2002;Dobbs-Dixon et al. 2004;Fabrycky & Tremaine 2007;Zhou & Lin 2008;Jackson et al. 2008Jackson et al. , 2009Miller et al. 2009;Benítez-Llambay et al. 2011;Rice et al. 2012;Beaugé & Nesvorný 2012;Lanza & Shkolnik 2014), we find that the location of zero torque would be between 0.008 AU and 0.012 AU (at orbital periods between 0.40 d and 0.66 d), which is well inside the stellar corotation radius of the stellar evolution tracks and extremely close to the stellar surface. As a consequence, the lower stellar tidal dissipation assumed in previous studies cannot produce a tidal migration barrier in the type II disk migration regime.
Our nominal disk properties come with significant uncertainties. For example, literature values of Σ p,0 span orders of magnitude, ranging from a few times 10 g cm −2 (Menou & Goodman 2004) over some 100 g cm −2 (Klahr & Kley 2006;Gressel et al. 2013;Flock et al. 2017) to 1 000 g cm −2 and more at about 1 AU (Bell et al. 1997;Ida & Lin 2004;Kretke & Lin 2012;Kley & Nelson 2012). Details depend on the dimensionality and assumptions of the respective models, and we chose Σ p,0 = 1 000 g cm −2 to reproduce the minimum-mass solar nebula for a viscous, flaring disk with α = 10 −2 . Different models would yield somewhat different estimates of the tidal migration barrier and, as a consequence, of the hot Jupiter survival and formation rates.
In particular, our specific estimates of the tidal migration barrier rate around sun-like stars depend on our assumptions of the distribution of disk properties (log 10 (Σ p,0 /[g cm −2 ]) = 3 ± 1, log 10 (α) = −2 ± 1, and µ = 1.85 ± 0.55) and of the efficiency of stellar dissipation (log 10 ( D ω ) = −3.25 ± 0.5), which are free parameters in our model. We do have estimates for these quantities from simulations and observations, but detailed 3D magnetohydrodynamical simulations of the disk properties inside 0.1 AU around accreting sun-like stars and of the evolution of tidal dissipation are required to validate or improve them.
In Eqs. (5) and (6), we have adopted a conventional description of the disk torque on the planet in the type II migration regime, that is, we assume that the planet has separated the disk into an outer and an inner part with negligible flow between the two. Magnetic fields, however, could open up a magnetic cavity around the star, which would affect the disk torque and actually halt planet migration altogether (Lin et al. 1996;Kuchner & Lecar 2002). The critical distance for a cavity to form, referred to as the Alfvén radius (r A ), is determined by the equilibrium between the dynamic pressure of the disk matter and the magnetic pressure from the star's dipole field. For T Tauri stars of about 0.8 M and 2.5 R , the resulting value of r A ≈ 0.05 AU (Romanova & Lovelace 2006) coincides equally well with the observed hot Jupiter pile up as our theory of a tidal migration barrier. That said, 3D magnetohydrodynamic simulations showed that the mass inflow rate (and consequently Σ) interior to r A can be significant. As a consequence, planet migration might in fact not stop near r A (Romanova & Lovelace 2006) and an alternative mechanism would be required, which could be the star's tidal torque as we show.
The switch from positive to negative tidal torques at the stellar corotation radius, which is key to the tidal migration barrier, is only valid if both the eccentricity and the planet's spinorbit misalignment (its obliquity, ψ p ) are small. Indeed, this switch does not apply for large obliquities, and in particular for ψ p > 90 • (Barker & Ogilvie 2009;Damiani & Mathis 2018). In our calculations, the planet is assumed in the disk midplane and in the stellar equatorial plane. Many hot Jupiters, however, have actually been found in substantially misaligned orbits (Albrecht et al. 2012) and tides actually might have played a key role in their formation (Fabrycky & Tremaine 2007;Dawson 2014;Anderson & Lai 2018). Several ways to put misaligned hot Jupiters in the context of this study could involve planet-planet gravitational interaction (Naoz et al. 2011;Wu & Lithwick 2011), the Kozai-mechanism (Nagasawa et al. 2008), or a combination of stellar tides, the disk torque, and high-eccentricity migration (Wang et al. 2017) after an initial migration stranding at the tidal migration barrier. In fact, these would be compelling mechanisms to investigate in follow-up studies. Given our restriction to circular, non-oblique orbits, our model has limited applicability and in light of the various alternative mechanisms that could form close-in massive planets it might just be part of the solution of the formation of hot Jupiters.
In this paper we focus on the disk migration and stellar tidal interaction of hot Jupiters in the equatorial plane of their stars. A significant fraction of the hot Jupiters, however, exhibits substantial spin-orbit misalignments (Winn & Fabrycky 2015) and the corresponding tidal torques from their host stars differ significantly (Lai 2012) from the zero-obliquity framework of this paper. Either these misalignments are created after the co-planar formation of hot Jupiters, e.g. through planet-planet scattering or the Lidov-Kozai mechanism, or the planets already formed in a warped or tilted protoplanetary disk. Since these processes act on different time scales and since the resulting stellar torque in oblique orbits is quite different from the formulation used in this paper, it might be possible that they can be observationally tested although this treatment is beyond the scope of this study. For example, (Morton & Johnson 2011) studied the outcomes of the spin-orbit evolution in the Lidov-Kozai migration scenario (Fab-rycky & Tremaine 2007) and in the planet scattering scenario (Nagasawa et al. 2008) and found evidence for the scattering model being the better explanation of the observed distribution of projected obliquities.
We have ignored in our model the effects of possible mass loss from the planet due to the Roche lobe overflow (RLO). RLO initiates when the Roche lobe around the planet R R ≈ a(M p /M p ) 1/3 becomes as small as the planetary radius. Even in the case of a highly inflated, young gas giant with a radius of 2 R J , RLO at R R ≈ 2 R J will only be reached at a semimajor axis a ≈ 0.01 AU, that is, very close to the stellar surface, which is at 0.0093 AU for a 1 Myr young solar type star. Hence, there will be a very sensible fine tuning involved between extreme planetary inflation (R p > 2 R J ), a sufficiently small star (R < 2 R ), and sufficiently weak tidally-driven orbital decay to prevent an inward spiralling hot Jupiter from falling into the star by RLO. Some of the known planets might indeed be or have been affected by RLO (Jackson et al. 2016) and that the inclusion of RLO could be relevant interior to the inner cavity of our model. The modelling of RLO, however, requires some parameterization (or even modeling of) planetary structure, which is beyond the scope of this paper.
As for the stellar spin-up, in our scenario of planet migration under the effect of stellar tides, hot Jupiters would either fall into their star within the first ≈ 10 Myr of their lifetime if the tidal torque is too weak to stop migration. Or they would stop beyond the stellar corotation radius at the tidal migration barrier, where they would act to slow down the stellar rotation. In our picture hot Jupiters spin up their stars if they fall through the corotation radius and get swallowed by the star or they survive and act as a rotational brake during the first ≈ 100 Myr of their star's lifetime.
Tidally excited waves in the stellar radiation zones (Goodman & Dickson 1998;Ogilvie & Lin 2007;Chernov et al. 2017) or the wave breaking mechanism at the stellar center (Barker & Ogilvie 2010) could trigger additional tidal dissipation beyond the processes in the convective envelope considered in this study. Further refinements of this theory could be achieved through the consistent modeling of the radial profile of the stellar density, the latter of which is assumed to be constant (though different) in both the stellar core and the envelope in our model (Ogilvie 2013).
We have verified the findings of Bolmont et al. (2017) that more metal-rich stars exert a stronger tidal torque on their closein hot Jupiter companions. This result could, at least partly, explain the finding of a lower abundance of massive planets (M p > 4 M Jup ) in short-period orbits (< 10 d) around low-metallicity stars by Narang et al. (2018). In other words, this observed low abundance could be due to the lower tidal dissipation rates in young solar mass low-metallicity stars and the resulting lower migration survival rates of the planets.

Conclusions
We present a new model for the formation of hot Jupiters under the combined effects of the dynamical tide in the star's convective envelope and type II planet migration. First, we calculate the nominal tidal torques of a young, highly dissipative solar-type star and of a 2D viscous disk in thermal equilibrium with radial temperature, gas surface density, and viscosity dependences acting on a Jupiter-sized migrating planet. We study a system of a 2 R star with a range of possible spin periods according to observations of young star clusters with a disk gas surface density similar to the minimum-mass solar nebula scale The histogram of the observed hot Jupiter pile-up in Fig. 1 is shown on a logarithmically scaled abscissa with constant bin width of 0.01 AU. In this representation, bins appear wider in close-in orbits and thinner in more distant orbits. As a consequence, the proposed hot Jupiter pile-up could simply be a binning artefact. Figure A.1 shows the same data as Fig. 1 and again on a logarithmically scaled abscissa, but now using a logarithmic bin width as well. In this representation, the bins appear to have constant width in the plot, although the effective bin width really depends on the semimajor axis. Near the pile-up at 0.05 AU, for example, the bin width is about 0.002 AU, whereas near 1 AU the bin width is roughly 0.05 AU.
Technical details aside, the most important conclusion to be drawn from Fig. A.1 is that the hot Jupiter pile-up is not a binning artefact. In this representation using logarithmic bin width, the maximum of the observed hot Jupiter distribution near 0.05 AU is about 5 times as high as it is at around 0.03 AU or 0.1 AU.

Appendix B: Lognormal randomization in gnuplot
The Monte Carlo simulations of Sect. 3 and the illustration of Fig. 4(b) were generated from a single gnuplot script (pile-up.gp; Heller 2018) and using gnuplot version 5.2. Both the PDF file of Fig. 4(b) and the estimated hot Jupiter survival rate of 31 % are direct outputs from this script.
The implementation of the randomized drawings in gnuplot code deserves some explanation because, to the best of the author's knowledge, there exists no simple way in gnuplot to sample a random variable from a given probability distribution. The aim was to sample Σ p,0 , D ω , α, and µ based on the probability distributions of each of these random variables. As an example, consider µ = 1.85 ± 0.55, where 1.85 is the mean value and 0.55 is the standard deviation (1 σ) of a normal distribution around the mean. The symmetric interval of ± 1 σ around the mean contains about 68.27 % of all realizations for a large number of samples.
Although gnuplot does not have a built-in function to directly sample a probability distribution, it does have a built-in function rand(0) to generate a random real number within [0,1] with constant probability density throughout the interval. We can combine rand(0) with the built-in function invnorm(), which is the inverse function of the cumulative normal distribution e −x 2 /2 , (B.1) for our purpose. invnorm() is defined within [0,1] and in particular we have invnorm( norm(x) ) = x. The operation y = invnorm( rand (0) ) is then equivalent to a random sampling of values of y according to a probability density that is given by a normal distribution. In Fig. B.1, we show invnorm(x). Note that invnorm(0.5 ± 1σ/2) = ±1, lim x → 0 (invnorm(x)) = −∞, and lim x → 1 (invnorm(x)) = +∞. The 1σ confidence interval extends from x = 0.5 − 1σ ≈ 0.159 to x = 0.5 − 1σ ≈ 0.841 on the abscissa and from y = −1 to y = +1 on the ordinate. As a consequence, a large number of randomized realizations y = invnorm( rand(0) ) will produce a normal distribution of y with a mean of 0 and a standard deviation of 1. We can scale the width of the standard deviation by multiplication of invnorm( rand(0) ) with the desired 1σ value. As an example, the gnuplot implementation of our randomized drawings of µ = 1.85 ± 0.55 reads W = 0.55 * invnorm( rand(0) ) mu_RAND = (1.85+W) where the temporary variable W is one particular realization from the normal distribution.