Low-Eccentricity Migration of Ultra-Short Period Planets in Multi-Planet Systems

Recent studies suggest that ultra-short period planets (USPs), Earth-sized planets with sub-day periods, constitute a statistically distinct sub-sample of {\it Kepler} planets: USPs have smaller radii ($1-1.4R_\oplus$) and larger mutual inclinations with neighboring planets than nominal {\it Kepler} planets, and their period distribution is steeper than longer-period planets. We study a"low-eccentricity"migration scenario for the formation of USPs, in which a low-mass planet with initial period of a few days maintains a small but finite eccentricity due to secular forcings from exterior companion planets, and experiences orbital decay due to tidal dissipation. USP formation in this scenario requires that the initial multi-planet system have modest eccentricities ($\gtrsim 0.1$) or angular momentum deficit. During the orbital decay of the inner-most planet, the system can encounter several apsidal and nodal precession resonances that significantly enhance eccentricity excitation and increase the mutual inclination between the inner planets. We develop an approximate method based on eccentricity and inclination eigenmodes to efficiently evolve a large number of multi-planet systems over Gyr timescales in the presence of rapid (as short as $\sim 100$~years) secular planet-planet interactions and other short-range forces. Through a population synthesis calculation, we demonstrate that the"low-$e$ migration"mechanism can naturally produce USPs from the large population of {\it Kepler} multis under a variety of conditions, with little fine tuning of parameters. This mechanism favors smaller inner planets with more massive and eccentric companion planets, and the resulting USPs have properties that are consistent with observations.


INTRODUCTION
The existence of ultra-short period planets (USPs), defined to be small planets (R ≤ 2R ⊕ ), with sub-day periods (i.e. P ≤ 1 day) is one of the major surprises in exoplanetary astrophysics. The first example of such planets, CoRoT-7b, was discovered in 2009 (Léger et al. 2009); since then, about a hundred USPs have been found by various transit surveys (Sanchis-Ojeda et al. 2014), and the overall occurence rate of USPs sits at ∼ 1%, a figure that is similar to the census of hot Jupiters, i.e. giant planets with P ≤ 10 days (Cumming et al. 2008;Wright et al. 2012). A few notable USPs have received special attention: 55 Cnc e (Dawson & Fabrycky 2010) with P = 0.74 days was the first discovered Super-Earth, Kepler-10b with P = 0.83 days (Batalha et al. 2011a) was the first terrestrial planet discovered by Kepler, and dN/d log P ∝ P α with α 1.5 − 1.7 (see also Petigura et al. 2018;Weiss et al. 2018), USPs followed a steeper trend with α ∼ 3.0. In addition, the normalization of the period distribution may also be different: the planet occurrence rate is discontinuous across the P = 1 day boundary, with ∼ 50% more planets with periods just below P = 1 days than just above.
In addition to their period distribution, USPs have other statistical properties that differentiate them from longer-period planets. USPs have smaller radii, with the vast majority having 1R ⊕ ≤ R ≤ 1.4R ⊕ ), a fact which may be attributed to photo-evaporation or 'boil-off' as the planets are intensely irradiated. Compared with the other Kepler planets, planet systems with USPs also appear to have higher mutual inclinations: Dai et al. (2018) found that transiting Kepler planets with a semi-major axis to stellar radius ratio a/R < 5 had an inclination dispersion of ∆θ ≈ 6.7 ± 0.7 degrees, while planets with 5 < a/R < 12 had only ∆θ ≈ 2.0 ± 0.1 (consistent with the overall figures for Kepler multis, see e.g. Tremaine & Dong 2012;Fang & Margot 2012;Fabrycky et al. 2014). This observation is further corroborated by the fact that for FGK host stars, USPs feature a factor of ∼ 8 fewer co-transiting external companions compared with their merely 'short-period planet' (SP) counterparts Weiss et al. 2018), and when USPs do have external transiting companions, the period ratios between the USP and their closest companion is P 2 /P 1 15, a value that is nearly an order of magnitude above the typical period ratios of 1.3 − 4.0 seen in Kepler multis (see also Steffen & Farr 2013).
The statistical distinctness of USPs and their unusual locations so close to their host stars defy conventional understandings of planet formation, and the origins of these planets remains a mystery. USPs may sit in the short-period tail of the distribution of close-in rocky planets that formed in-situ through core accretion (Chiang & Laughlin 2013), or they may have migrated to their current locations from initially more distant orbits (Ida & Lin 2004;Schlaufman et al. 2010;Terquem 2014). In the latter scenario, they (like most of their SP bretheren) likely would have begun their lives as Super-Earths/Mini-Neptunes with a gaseous H/He envelope that was subsequently lost to photo-evaporation (Valencia et al. 2010;Owen & Wu 2013). To shove the planets very close to their host stars, some form of disk migration and/or tidal dissipation is required. Lee & Chiang (2017) considered stellar tides raised by the planet, treating the stellar tidal quality factor Q as a free parameter; since the orbital decay rate is proportional to the planet mass, the required Q value to induce significant decay of small planets would make hot Jupiters at P ∼ 1 day "disappear" on a short timescale. In the case of planetary tides, the proto-USP must maintain a finite eccentricity in order to undergo orbital decay. Petrovich et al. (2018) examined a high-eccentricity migration scenario in which a proto-USP attains large eccentricity due to secular chaos in a hierarchical system with N > 3 planets, followed by orbital circularization due to planetary tides. They also briefly explored the possibility of forming USPs through secular interactions with eccentric giant planet companions, but dismissed the possibility as unlikely: they found that producing USPs usually required eccentric giant planet companions with P ≤ 10 days, a requirement at odds with the observation that the presence of USPs do not seem to be correlated with the stellar host metallicity, and therefore by proxy the occurrence of giant planet companions (Winn et al. 2017). In short, although these previous ideas indeed may produce ultra-short period planets under some conditions or assumptions, there is yet no firm evidence that USP formation can be completely accounted for by any one of the aforementioned scenarios.
Indeed, the formation mechanism of USPs remains unclear and this is the question we aim to address. The main thesis of this work is that the combination of secular interactions and tidal dissipation in multi-planet systems is likely be the most natural and efficient way to generate USPs. This mechanism requires small, rocky planets to be born at moderate eccentricities (i.e. e 0.1) in multi-planet (N ≥ 3) systems, but otherwise requires little fine tuning of planet parameters. Empirical studies suggest the orbital eccentricities of Kepler multis have a dispersion of σ e ∼ 0.05 (Xie et al. 2016;Van Eylen & Albrecht 2015;Van Eylen et al. 2018), so USP formation in this mechanism would occur at the tail end of the eccentricity distribution. However, note that the currently observed eccentricity distribution has been damped over Gyrs by tidal dissipation (Hansen & Murray 2015), and the primordial eccentricities may be much larger.
Certainly, the idea of secular forcings coupled with tidal dissipation is not a new one; it has already been applied to short-period exoplanet systems in various contexts (e.g. Wu & Goldreich 2002;Mardling 2007Mardling , 2010Batygin et al. 2009;Hansen & Zink 2015;Petrovich et al. 2018), although this work is the first to tackle the problem in the context of USP formation in multi-planet systems. The mechanism studied in this paper has some similarity to the secular chaos mechanism proposed by Petrovich et al. (2018), but with important differences: Whereas Petrovich et al. (2018) rely on secular chaos driving small planets to attain large eccentricities (e.g. 1 − e 1) and thereby small pericenter distances to achieve USP formation, our mechanism requires the inner planet (initially at P = 1 − 3 days) to achieve only mild initial eccentricities (e ∼ 0.1) through secular interactions; although the mechanism proposed by Petrovich et al. (2018) allows for a more diverse proto-USP period, it also requires the presence of several well-separated exterior planets, whereas in our mechanism, the initial proto-USP period is more constrained, but the external planet companions are allowed more lee-way in terms of their spacing. In light of this fact, we call our proposed formation mechanism the 'low-e migration' of USPs.
In this paper, we present an investigation on the prospects of generating USPs through 'low-e migration'. The structure of the paper is as follows. In section 2, we present the mathematical formalism for the eccentricity and orbital evolution of a multi-planet system undergoing secular interactions and tidal dissipation. As we discuss below, a bruteforce approach to this problem is impractical, and we derive the evolution equations in the framework of eigenmodes in section 2.1. In section 3, we apply our formalism to the case of 2-planet systems, deriving semi-analytical results for the eccentricity and mode evolution; in section 3.2 we discuss the criterion for USP formation to occur in such 2-planet systems. We then extend these results to the case of 3-planet systems in section 4, and show that such systems allow for successful low-e secular migration under reasonable conditions. In section 5 we consider the inclination evolution of the planets, taking into account planet-planet coupling as well as interactions with stellar spin. The results of the preceding sections are synthesized into a population model in section 6 -readers who are most interested in observational implications of our results may skip to this section. In section 7, we discuss the feasibility of low-e secular migration and compare it against other proposed mechanisms. Finally, a summary of our work is provided in section 8.

ECCENTRICITY EVOLUTION AND ORBITAL DECAY: FORMALISM
Consider a N-planet system with individual planet mass m i , semi-major axis a i , initial orbital eccentricity e i , longitude of periapsis i , inclination θ i and longitude of the ascending node Ω i , where i ∈ [1, N] is the planet index, orbiting a host star with mass M and radius R . We assume that the host-star is Sun-like, i.e. M = M and R = R . The planet's semi-major axis is related to the orbital period by a 0.0196(P/day) 2/3 (M /M ) 1/3 au. The dynamical evolution of the system is governed by the interplay of several effects: planet-planet secular perturbations, General Relativistic (GR) periastron advance, spin-orbit coupling due to stellar oblateness, planetary tides and stellar tides. We define E i ≡ e i exp (i i ) to be the complex eccentricity of the ith planet, and define the eccentricity vector of the N-planet system as In the linear (Laplace-Lagrange) theory, the time evolution of ì E is governed by the equation where the coefficients of the time-varying N × N matrix H(t) is given by Here the complex "frequencies"ω i (taking into account the eccentricity damping due to tidal dissipation in the i-th planet) is defined as The quantities ω i j and ν i j are the quadrupole and octupole precession frequencies of the i-th planet driven by the actions of the j-th planet, given by and where we have defined a < = min(a i , a j ), a > = max(a i , a j ), α = a < /a > , L i = m i √ GM * a i is the (circular) angular momentum of the i-th planet, and b (n) 3/2 (α) are the usual Laplace coefficients given by In the limit α 1, the first-order expansion of the Laplace coefficients are b (2) 3/2 (α) 3α and b (2) 3/2 (α) 15α 2 /4. The general relativistic apsidal precession frequency is given by (for e i 1): where n i is the angular orbital frequency. The rate of periastron advance on the i-th planet due to its tidal bulge is given by with k 2,i being the tidal Love number of the i-th planet; in this work, we adopt a value of k 2,i = 1. Generally, as the planet moves inwards, the GR and tidal forces become increasingly important, while farther out, planet-planet secular interactions tend to prevail. We use the weak friction theory of equilibrium tides to describe tidal dissipation in the planet (Darwin 1880;Alexander 1973;Hut 1981). The eccentricity damping rate of the i-th planet due to tidal dissipation is where ∆t L,i is the tidal lag time of the i-th planet, and is related to the tidal quality factor Q i by Note that in this formalism, the value of Q i is not constant and instead varies with the planet's semi-major axis. In the Solar System, values of Q range from 10 − 500 for terrestrial planets and satillites, but the gas giants have values of Q that are much larger (Goldreich & Soter 1966;Ogilvie 2014). We assume proto-USPs to be predominantly rocky and adopt values of Q 1 in the range between 70 and 700, while the exterior planets are assumed to have H/He envelopes comprising a few percent of the planet's mass (but a large fraction of the radius) and therefore have much larger values of tidal Q i . Because the the outer planets (i ≥ 2) have much larger values of a i and Q i , the effects of tidal dissipation are much weaker for these planets. Therefore, we simplify the problem by only considering tidal dissipation for the innermost planet (i.e. by setting Q i = ∞ for i ≥ 2) in sections 2 through 5; the effects of tidal dissipation in the outer planets are included in our population synthesis study in section 6.
To completely determine the time evolution of the system, Eq. (2) for ì E should be supplemented by the evolution of the inner planet's semi-major axis a 1 : a 1 a 1 tide = −2γ 1 e 2 1 = −1.9 × 10 −9 k 2,1 ∆t L,1 100s We also consider the orbital decay driven by dissipation of tides raised on the host star by the planet (Goldreich & Soter 1966): where Q = 3Q /(2k 2, ) is the reduced tidal quality factor of the star. Empirical measurements by Penev et al. (2018) suggest a value of Q = 10 7 at a tidal forcing frequency of 2 days −1 , decreasing to Q = 10 5 when the forcing frequency is 0.5 day −1 . Lee & Chiang (2017) treated Q as a free parameter, and considered values of Q in the range of 10 6 − 10 8 . Thus, in general, the orbital decay rate of the inner-most planet is given by In sections 2 -5 we will focus on planetary tides and neglect the effect of stellar tides (i.e. by setting Q = ∞), although we will include stellar tides in our population synthesis study in section 6. We do this for two reasons: Firstly, for typical values of Q the effect of stellar tides is small over the lifetime of the system and only attains significance when an USP has already been produced (i.e. a 1 ≤ 0.02 au.), therefore its role is orthogonal to the aims of this work. Secondly, the addition of stellar tidal dissipation destroys the conservation of orbital angular momentum, an otherwise desirable property of Eqs.

Eccentrity Evolution in the Framework of Eigenmodes
A brute-force integration of Eqs.
(2) and (14) encounters difficulty: the relevant frequencies vary over many orders of magnitude with the orbital decay timescale (|a 1 / a 1 | 1 Gyr) being much longer than the precession timescales (2π/ω i ∼ 10 3 yrs); the "stiffness" of the equations make it impractical to integrate a large number of systems. Our approach is eschew calculating the phase of the eccentricity vector. We do this by decomposing the planet eccentricities into eigenmodes. We define the eigenvalue λ α and eigenvector ì E α (with modes denoted using Roman Numerals α ∈ [I, II, III...]) of the system as where Note that since H(t) depends on a 1 (the time evolution of a 2 , a 3 ... are negligible), the eigenvalue λ α (t) and eigenvector ì E α (t) evolve in time as a 1 decreases. We now introduce the matrices G(t) and V(t) formed from the eigenvalues and eigenvectors of H(t): and By definition, the matrices G(t) and V(t) satisfy the identity In general, the time evolution of ì E can be written as a superposition of eigenmodes, where ì A(t) ∈ C N is the vector of eigenmode amplitudes: The initial condition ì A(0) is given by Substituting Eq. (20) into Eq. (2), and using the identity HV ì A = VG ì A (which follows from Eq. 19), we find The above equation is exact, but still involves highly oscillatory complex mode amplitudes. To make further progress, we note that Eq. (23) yields the evolution equation for | A α (t)|: In general, A * α and A β contain fast varying phases. We now adopt the ansatz that when averaging over timescales much longer than |λ α −λ β | −1 but shorter than the orbital evolution time |a 1 / a 1 |, With this ansatz, Eq. (24) reduces to and the magnitude of the planet eccentricity is given by If we define B α (t) ≡ | A α (t)| ∈ + and let ì B(t) be the vector with components B α (t), Eq. (26) can be written in matrix form as Eq. (28) is much easier to solve numerically than the exact Eq. (23), because both the vector ì B and the matrix W are explicit functions of ì a 1 (t), and only vary with t as a 1 (t) varies. In particular, we can evaluate W from with a 1 given by Eq. (14). We solve Eq. (28) combined with Eqs. (14) and (27) to obtain the time evolution of the RMS eccentricity and semimajor axis. Although our formalism above does not capture the short-term oscillations in eccentricity, it is nonetheless possible to know the extent of these oscillations by computing the "instantaneous" maximum and minimum eccentricity. The instantaneous maximum eccentricity is given by and the minimum eccentricity is given by In the above RMS-averaged formulation, the ansatz leading to Eq. (28) is equivalent to assuming that the mode amplitude evolves adiabatically as a 1 decreases, i.e. we assume that the cross-terms corresponding to the mixing between modes average out to zero due to their incoherent phases, and only diagonal terms remain. In reality, the assumption in Eq. (25) may not hold in the later stages of orbital decay, as certain pairs of modes may become locked in either alignment or anti-alignment depending on the configuration of eigenvectors. In practice, this turns out not to be an issue since all but one mode will have decayed away by this point, leaving the question of how to handle the cross-mode terms moot, and we have found excellent agreements across the board between the approximate RMS-averaged formulation and the exact treatment. Nonetheless, the approximation made in Eq. (25) is the main source of uncertainty in our approximate formulation and may lead to errors in edge cases when modes do not vary sufficiently rapidly relative to the orbital decay timescale. For the blue curve, we also show its limiting cases: the thin, solid blue curve corresponds to the approximation given by e 1,forced /e 2 ν 12 /ω 12 , the dash-dotted curve corresponds to e 1,forced /e 2 ν 12 /ω 1,gr and the thin dashed line corresponds to e 1,forced /e 2 ν 12 /ω 1,tide .  Figure 2. Sample time evolution for a two-planet system with initial conditions a 1,0 = 0.03 au, a 2 = 0.08 au and e 1,0 = 0, e 2,0 = 0.15. The left panels have planet masses m 1 = M ⊕ , m 2 = 3M ⊕ while the right panels have m 1 = M ⊕ and m 2 = 30M ⊕ . Note here we adopt a large value of ∆t L,1 = 10 5 s in order to speed up numerical calculations. In the middle and bottom panels, the dashed curves show our approximate solution using Eq. (28) while the solid curves are obtained using a direct integration of Eqs.
(2) and (14). On the left panels, the orbital decay of m 1 is limited by the amount of total angular momentum deficit, as the two modes decay away before substantial orbital decay can take place, whereas on the right panels, the planets have sufficient AMD to undergo substantial orbital decay, and is instead limited by the rate of orbital decay.

Mode Properties and General Evolution Behaviors
We demonstrate the application of the formalism presented in section 2.1 by considering 2-planet proto-USP systems. In this case, analytic expressions for the modes can be derived explicitly, providing useful insight into the more general multi-planet systems. The complex eigenfrequencies λ I , λ II (see Eq. 17) are given by λ I = 1 2 ω 1 +ω 2 + ∆ω 2 + 4ν 12 ν 21 (32) where ∆ω ≡ω 1 −ω 2 (withω 1 = ω 1 + i γ 1 andω 2 ω 2 ), and the eigenvectors are In general, in order for the inner planet (m 1 ) to become an USP, one requires that L 1 L 2 (where L i = m i √ GM a i is the circular angular momentum). In this limit the above expressions simplify considerably: Since ω 2 , ν 12 , ν 21 ω 1 (recall that ω 2 ω 21 = ω 12 L 1 /L 2 ), the eigenfrequencies become The eigenvectors in this limit are given by Since ν 21 /ω 1 = (L 1 /L 2 )(ν 12 /ω 1 ) 1, it is clear that the mode α = I (II) is associated with the free oscillation (apsidal precession) of the inner (outer) planet. The damping rate rate of the two modes are given by Clearly, the decay of eigenmode II is substantially supressed relative to mode I. It is therefore safe to assume that any initial oscillation along mode I is quickly damped out, and the system is locked into mode II. At this stage e 1 is given by the forced eccentricity: e 1 = e 1,forced ν 12 ω 1 e 2 = ν 12 e 2 ω 12 + ω 1,gr + ω 1,tide .
Thus, for t γ −1 1 , the orbital evolution of the inner planet is governed by (neglecting stellar tides) Comparing | a 1 /a 1 | with γ II (Eq. 39), we see that the system may exhibit two possible outcomes, depending on the system parameters and initial conditions: (i) for L 1,0 /L 2 2e 2 2,0 (where the subscript '0' referring to the initial value), mode II does not experience significant damping (i.e. e 2 e 2,0 ), and the inner planet keeps undergoing orbital decay until its forced eccentricity is suppresed by GR and tides, dramatically slowing any further tidal decay. (ii) For L 1,0 /L 2 2e 2 2,0 , both modes are eventually damped out, preventing further tidal decay and leaving behind two planets with circular orbits and fixed semi-major axes. Figure 2 shows examples of time evolution for the two cases. The solid curves are direct integrations of Eq. (2) while the dashed curves utilize the approximate formulation in Eq. (28); to speed up the integration of Eq.
(2) we have adopted an unphysical value of ∆t L,1 = 10 5 s, corresponding to Q 1 = 0.07 (at P 1 = 1 day). Note the excellent agreement between the exact treatment and our approximation. The left three panels of Fig. 2 shows an example of case (ii) -we see that in this case the mode amplitudes ì B(t) decay away to zero well before the nominal orbital decay timescale |a 1 / a 1 |. In fact, this is a general barrier to USP formation: the planets must have sufficient initial eccentricities (i.e. angular momentum deficit or AMD, see below) in order to have substantial orbital decay prior to having all the mode amplitudes dissipated away.
The right three panels of Fig. 2 show an example of case (i). In this case the inner planet is able to undergo substantial decay to become an USP due to the greater amount of AMD stored in the more massive exterior planet, and substantial oscillation amplitude remains in mode II even after the inner planet has decayed to sub-day periods. We see that in this case the tidal decay is still self-limiting: As the inner planets decay further, the effects of short-range forces become important, which forces the inner planet to attain much lower e 1,forced (see Eq. 40) that dramatically slows down the rate of tidal decay.

Criteria for Orbital Decay
The analysis and examples shown in Section 3.2 show that generally, two criteria must be met in order for the inner planet to undergo substantial tidal decay: (i) the total angular momentum deficit (AMD) 1 of the system AMD L 2 e 2 2,0 /2 must be sufficiently large to allow the inner system to undergo orbital decay before all the eccentricities are decayed away. (ii) The inner planet must have sufficiently large forced eccentricity such that the orbital decay occurs within the lifetime of the system. The first criterion arises from the conservation of the total angular momentum. Indeed, for e 2 i 1, the total angular momentum 1 The AMD of a system is given by the difference between its total angular momentum if all planet orbits were circular and its actual The critical eccentricity e 2,0 needed to meet the AMD (dashed curves) and tidal decay time constraints (solid curves) for USP formation in a two-planet system is plotted as a function of the companion semi-major axis a 2 . The inner planet has initial semi-major axis a 1,0 = 0.03, mass m 1 = M ⊕ , radius R 1 = R ⊕ , and tidal lag time ∆t L,1 = 100 s. The dashed curve are given by Eq. (46), corresponding to a 1,min /a 1,0 = 1 2 , while the solid curves are given by Eq. (47), corresponding to | a 1 /a 1 | 10 −10 yr. The red, blue, green, magenta and cyan curves correspond to m 2 = 7, 10, 15, 30 and 300 M ⊕ respectively. For a given m 2 , in order for efficient orbital decay to occur, the outer planet's initial eccentricity must be above both curves of the corresponding color.
where we have omitted the orbital decay of the other planets. Substituting Eqs.
(2) -(3), and noting that L 1 ν 12 = L 2 ν 21 , we find that Thus L is constant when γ is negligible (which is the case until a 1 is already reduced to a value well below 0.02 au by planetary tides). The semi-major axis of the inner planet decreases at the expense of the planet eccentricities, while keeping the total angular momentum constant. Assuming initially e 1,0 = 0, we have For a given a 1,0 , a 2 (= const.) and e 2,0 , the minimum semimajor axis the inner planet can reach (after indefinite time) is given by Thus, the critical initial eccentricity of m 2 required for significant semi-major axis decay (a 1,min /a 1,0 1 2 ) is given by The second criterion pertains to the orbtial decay timescale. In order for a 1 to decrease significantly within 10 10 yrs, the inner planet must have (see Eq. 12) e 1 4.6 × 10 −3 ∆t L,1 100s The above requirement is in fact overly conservative, since the rate of semi-major axis decay tends to accelerate as a 1 decreases until short-ranged forces become dominant (see Eq. 12). Using e 1 e 1,forced (Eq. 51), Eq. (47) translates into another constraint on e 2,0 as a function of a 2 /a 1,0 .
Thus, in the 2-planet case, the formation of USPs is limited by two constraints given by Eq. (46) (the AMD constraint) and Eq. (47) (with e 1 = e 1,forced , the decay time constraint). These constraints are shown in Fig. 3. The combination of these two constraints make the formation of USPs from 2-planet progenitor systems a challenging prospect; one natural way around the two barriers is to consider the effects of an additional external planet -we examine the 3-planet case in section 4.

Set-up
In this section we perform a systematic study of USP formation in 3-planet systems. The goal here is to gain physical insights on the dynamical evolution of such systems. In section 6 we perform population synthesis model to assess whether our model can reproduce the observed USP demographics.
We consider a 3-planet systems where the proto-USP has an initial semi-major axis a 1,0 in the range between 0.02 and 0.04 au, corresponding roughly to a period of P 1,0 = 1 − 3 days. The proto-USP is assumed to have an Earthlike composition with radius R 1 ∈ [1, 1.4] R ⊕ with mass given by m 1 = (R 1 /R ⊕ ) 4 M ⊕ , a scaling consistent with purely rocky compositions (Zeng et al. 2016). The tidal lag time ∆t L,1 is taken to be 100 or 1000 s, corresponding to Q 1 = 70 and 7 (at P 1 = 1 day). The outer planets have semi-major axis a 2 , a 3 ∈ [0.04, 0.2] au, masses m 2 , m 3 ∈ [3, 30] M ⊕ and initial eccentricities e 2,0 = e 3,0 ∈ [0.05, 0.3], and their tidal dissipation is neglected. Note that these system parameters are chosen expediently to showcase the evolution behavior when mode mixing occurs. Some of the illustrated systems may be dynamically unstable; we discuss this issue in Sec. 6.
As discussed in section 2.1, due to the "stiffness" of Eqs.
(2) and (14), we integrate Eq. (28) for all our systems. For comparison purposes, we also integrate the set of systems using Eq. (2), but with an enhanced value of ∆t L,1 = 10 5 s, corresponding to an initial value of Q 1 = 0.07 (at P = 1 day); these integrations are compared to our approximate method (based on Eq. 28) using the same value of ∆t L,1 . We find that our approximate formulation achieves excellent results across the entire parameter space we consider.  , eccentricities e i (middle) and inner planet semi-major axis a 1 (bottom) as a function of scaled time for three different three-planet systems. Here, the time evolution is obtained from Eq. (28). The thick dashed, thin solid and thin dashed lines correspond to values of ∆t L,1 = 10 2 , 10 3 and 10 5 s respectively. The time is scaled ast ≡ (∆t L,1 /100s) t. For the top three panels, the red, green and blue curves correspond to modes I, II and III respectively, while for the middle panels they correspond to e 1 , e 2 and e 3 . The three columns have different values of a 2 (0.043, 0.059 and 0.089 au respectively, from left to right), but otherwise identical parameters. The planets have masses m 1 = M ⊕ and m 2 = m 3 = 17M ⊕ and initial semi-major axes a 1,0 = 0.03 au and a 3 = 0.10 au. The inner planet's radius is R 1 = R ⊕ , and stellar dissipation is neglected (γ = 0). The initial planet eccentricities are given by e 1,0 = 0 and e 2,0 = e 3,0 = 0.15, with the initial longitude of pericenter for all three planets set to i,0 = 0.
Note that the time evolution of systems with artificially reduced values of Q 1 cannot simply be considered timescaled versions of systems with more realistic values of Q 1 : when Q 1 1, the presence of a large imaginary component to the matrix H(t) substantially modifies the structure of eigenmodes, causing the inner planet eccentricity to become quantifiably different. We demonstrate this in section 4.2.

Time Evolution & Mode Mixing
For the 3-planet case, the evolution in the framework of eigenmodes becomes considerably more complicated, and explicit analytic expressions are no longer possible. Instead, one must resort to numerical solution for the eigenfrequencies and eigenmodes. Nonetheless, the general features from the 2-planet case carry over. The mode associated with the free eccentricity oscillation of the inner-most planet tends to be damped away rapidly, while the other two modes damp on much longer timescales. However, during the evolution the eigenvalues of the three modes may cross one another, leading to substantial mode mixing. Such mixings correspond to secular resonances, causing an enhancement in the eccentricity of the inner planet and potentially speeding up its orbital decay by orders of magnitude (see also Hansen & Murray 2015). This effect is most prominent when L 1 is much less than L 2 or L 3 ).
In Fig. 4, we depict some examples of the time-evolution of the mode amplitudes (B I , B II , B III ), the eccentricities (e 1 , e 2 , e 3 ) and the inner planet semi-major axis a 1 (t) for three hypothetical proto-USP systems. The three systems have the same parameters and initials conditions except for different values of a 2 . For each system, we consider three values of ∆t L,1 ∈ [10 2 , 10 3 , 10 5 ] s, corresponding to Q 1 = [70, 7, 0.07] (at a period P 1 = 1 day). Notice that for the cases with ∆t L,1 = 10 2 , 10 3 s, the curves (with a scaled time axis) lie right on top of each other; this shows that as long as Q 1 1, inner planets with different values of Q 1 will undergo identical time evolutions if the time is scaled ast = t(∆t L,1 /100s). On the other hand, the case with ∆t L,1 = 10 5 s (the thin dashed lines in Fig. 4), corresponding to an unphysical value of Q 1 = 0.07 (at P = 1 day), shows qualitatively different eccentricity and time evolutions. This demonstrates that for our parameter space, a naive approach of simply integrating Eq. (2) directly with re-scaled values of Q 1 1 would give rise to incorrect results. We now focus on the two cases with physical values of ∆t L,1 (10 2 and 10 3 s), i.e. the thick dashed and thin solid lines of Fig. 4. In Fig. 5 we show the evolution of the eigenvalues and eigenmodes for the same systems. In the left column, the system displays no mode mixing, each mode decays independently and the proto-USP reaches a final value of a 1,f = 0.016 au. In the middle column (with a 2 = 0.059 au), a resonance occurs between the two faster modes at a 1 0.027 au (see Fig. 5), this causes e 1 to increase temporarily, followed by a rapid decrease in both e 1 and a 1 . After passing this resonance, the system continues to evolve, with each mode decaying independently until reaching a final value of a 1,f 0.014 au 2 . In the right column (with a 2 = 0.087 au), the modes initially decay smoothly and independently of each other. A resonant mode crossing occurs at a 1 0.021 au between the two slower modes. This causes e 1 to reverse course and increase sharply, followed by a rapid decrease in both e 1 and a 1 . The inner planet reaches a final a 1,f 0.021 au.
To show that our method accurately captures the evolution of a resonant system, in Fig. 6 we show a comparison of the right panels of Fig. 4 -5 with the result obtained by a brute force computation using the exact Eq. (2), a result which took 5 days to complete on a Ryzen 1700 processor. Overall, there is an excellent agreement between our approximate results and the results obtained by the brute-force integration of Eq. (2).
Although the presence of the secular resonance can help to speed up tidal evolution of proto-USP systems, it is not a necessary condition to form USPs. In the next section, we discuss the conditions under which USPs may form in 3-planet systems.

Criteria for Orbital Decay
As in the 2-planet case (section 3), the formation of USPs is constrained by two factors: (i) The amount of the total AMD to sustain the orbital decay, and (ii) the amount of forced eccentricity of the proto-USP in order to have tidal decay occur within the lifetime of the system. Angular momentum conservation implies that the minimum semi-major axis that can be attained by the inner planet is (cf. Eq. 45) Therefore, to achieve a 1,min /a 1,0 1 2 , one requires (assuming e 2,0 ∼ e 3,0 ) e 2,0 ∼ e 3,0 0.77 At the same time, analogous to the 2-planet case, to have efficient orbital decay within the lifetime of the system, Eq.
(47) must be satisfied. In the case of 3 planets, the inner planet forced eccentricity can no longer be expressed in a simple expression; one must solve numerically the eigenvalues and eigenvectors; the forced eccentricity can be obtained from the amplitudes of the two slower decaying modes Since A α is determined from the initial values of e 2,0 and e 3,0 , the constraint on the inner planet eccentricity corresponds to a constraint on the external planet eccentricities. In the limit that L 3 L 2 , L 1 , an approximate expression for the forced eccentricity is given by (see  e 1 e 1,forced = ν 12 ω 2 + ν 12 ν 23 ω 1 ω 2 − ν 12 ν 21 e 3 . The above equation is more accurate when the planets are spaced evenly and well-separated, and does not fully capture the resonant mode crossings. In general, e 1,forced tends to be greater than given by the expression above, due to the contribution of other modes and aforementioned resonances. In Fig. 7, we show the two constraints for USP formation in three-planet systems; this figure is analogous to Fig. 3, except with the addition of a third planet (with m 3 = m 2 and fixed a 3 ). We find two important differences between the constraints for two-planet systems (see Fig. 2) and three-planet systems: (i) ceteris paribus, the presence of an additional planet lowers the eccentricity values (e 2,0 , e 3,0 ) required to meet the AMD constraint; (ii) the decay time constraint can be met by a larger set of values of a 2 and e 3,0 , since the presence of two secular resonances makes it possible for e 1,forced to be large even for smaller values of e 2,0 and e 3,0 .
In general, for three-planet systems, the AMD constraint is more stringent than the decay time constraint. To illustrate this, we show the final value of a 1,f reached after 10 Gyr of evolution as a function of a 2 in Fig. 8 for threeplanet systems with varying initial values of a 1,0 , with the planet masses, a 3 and e 2,0 = e 3,0 fixed. The dashed curves in Fig. 8 correspond to the minimum possible value of a 1,f given by the AMD constraint, while the solid curves are their actual values at the end of the evolution. For systems with e 3,0 0.1, the solid curve comes very close to the dashed curve, indicating that the orbital decay of a 1 is being stalled by a lack of AMD. For systems with e 3,0 0.15, the orbital decay instead becomes time-limited. Mode I Mode II Mode III 10 -12 10 -11 10 -10 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 Im(λ α ) 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 E α1 /E   . Note that there is some disagreement between the maximum and minimum values of our approximate method and the brute force calculation due to the fact that the system has insufficient time to reach the theoretical long-term extrema due to the rapid orbital decay.
The fact that USP production is more constrained by AMD has certain observational implications. One would expect USPs to be systematically lower in mass, as lower-mass inner planets are more likely to meet the AMD constraint (see Eq. 49). At the same time, we expect the external companions of USPs to have systematically larger masses, although giant planet companions are not required. To generate USPs efficiently, we also require the primordial planet eccentricities to be e 2,0 ∼ e 3,0 0.1, although their final values can be much lower due to tidal dissipation. The observational implications are explored in more detail in section 6, where we develop a population model for USP generation.
In this section. we have explored USP formation from three-planet systems. At first glance, there is a tension between USP generation from multi-planet systems and the fact that observed USPs have a dearth of exterior transiting companions compared with their non-USP counterparts. This prima facie contradiction can be rectified when we consider the mutual inclination evolution of USP-forming systems, in section 5.

INCLINATION EVOLUTION
We are interested in the inclination evolution of the proto-USP system because the evolution of the mutual inclination of planets determine the extent to which USPs will transit simultaneously with their companions, a quantity that can be observationally constrained (see Sec. 6). As we shall demonstrate in this section, there exists a secular mutual inclination 'resonances' that roughly coincides with the secular eccentricity resonance; systems that result in large excitations in eccentricity (and therefore forming USPs) should also expect large excitations in mutual inclination.
The inclination evolution of the proto-USP system proceeds in an analogous fashion as the eccentricity evolution. We define the complex variable I j = θ j exp(iΩ j ) for each planet, where θ j is the orbital inclination (relative to the initial orbital plane) and Ω j is the longitude of the ascending node. The mutual inclination θ i j between planets i and j is given by For convenience we define the inclination of the N-planet system as   47) and (50), with the eigenvectors being solved numerically and assuming that e 2,0 = e 3,0 . The two dips in the solid curves correspond to the two resonant mode crossings discussed in section 4.2. For a given m 2 = m 3 , in order for efficient orbital decay to occur, the outer planet's initial eccentricities must be above both curves of the corresponding color. Note that some values of a 2 may result in dynamically unstable systems.
For small inclinations (i.e. θ i 1) the time evolution of ì I is governed by where I is the complex obliquity of stellar spin. In Eq. (54), the first term in the RHS is due to secular planet-planet interactions, while the second term accounts for the nodal precession driven by the stellar spin; the vector ì ω is given by The N × N matrix H (t) is given by where ω i j is given by Eq. (5), and We also need to account for the evolution of stellar spin, governed by The nodal precession rate of the i-th planet driven by the stellar spin-induced quadrupole is where M , R and Ω = 2π/P are the stellar mass, radius and angular rotation frequency respectively. The constants k and k q are defined through the star's moment of inertia and quadrupole moment: I 3 = k M R 2 and I 3 − I 1 = k q Ω 2 M R 2 whereΩ 2 = Ω (GM /R 3 ) −1/2 . Typical values for solar type stars are k 0.06 and k q 0.01 (e.g. Lai et al. 2018). The ratio of the stellar spin angular momentum S = I 3 Ω to the orbital angular momentum of the i-th planet L i is The precession rate of the stellar spin driven by the i-th planet is For small planets (m 1 35M ⊕ ) at P 1 ∼ 1 day periods, the stellar spin angular momentum is much greater than the orbital angular momentum. Thus we can assume the stellar spin axis constantly points towards theẑ-axis, i.e. I 0. The rotation rate of the star Ω decreases due to magnetic braking. According to Skumanich (1972), Ω ∝ −Ω 3 , so that the time evolution of the spin rate is given by where Ω ,0 is the initial spin rate, and α MB is a constant, calibrated such that the rotation period reaches ∼ 30 days at an age ∼ 5 Gyr. For this section, we adopt a constant value of Ω in lieu of the Skumanich law to better control for the effect of stellar spin; the effect of a time-dependent stellar spin period is left for section 6. The above equations, coupled with the time-evolution of the planet eccentricities e i and inner planet semi-major axis a 1 fully describes the inclination evolution of the system in the linear regime (e i , θ i 1). Analogous to the case of eccentricity evolution, Eq. (54) involves terms that oscillate rapidly compared with the timescale of orbital decay, leading to a 'stiff' set of equations that resists brute-force simulations. In section 5.1 we address this issue by recasting the problem in the framework of eigenmodes.

Inclination Evolution in the Framework of Eigenmodes
In the discussion below we assume that the stellar spin axis is always along theẑ-axis. In this case Eq. (54) simplifies to We define the eigenvalue λ α and eigenvector I α with modes denoted using Roman Numerals (α ∈ [I, II, III...]) of the system as  Figure 10. Same as Fig. 9, except with the stellar spin period fixed at P = 1 day.
where we have define the vector eccentricity vector ì I to be As in section 2.1, we introduce the matrices G (t) and V (t) constructed from the eigenvalues (λ α ) and eigenvectors (I α ) of H (t): and The time evolution of ì I can be written as a superposition of eigenmodes where ì C(t) is the vector of eigenmode amplitudes: The time evolution of ì C is governed by The above equation is exact. Similar to the case of the eccentricity evolution (section 2.1), we bypass the stiffness of the above equation by focusing only the evolution of D α ≡ |C α |, whose evolution is given by [deg] Figure 11. Similar to the right panels of Fig. 9 (corresponding to m 2 = 17M ⊕ ), except we fix the value of a 2 = 0.055 au, and instead vary the value of P , which is fixed in time during the evolution.
Note here that W (t) can also depend on the spin-down of the star, i.e.
The instantaneous RMS inclination is given by The "instantaneous" maximum inclination is given by while the minimum inclination is Using the mode solution, we can also obtain the mutual inclination between planets. The RMS mutual inclination between planets i and j is given by while the maximum and minimum mutual inclinations are respectively given by

Resonance Crossing and Mutual Inclination Excitation
Significant mutual inclinations between the inner two planets can be excited when the planet system crosses a secular inclination resonance. When spin-orbit coupling is negligible (i.e. ω i is small), and L 3 L 1 , L 2 , the resonance occurs when the dimensionless "coupling parameter" 12 , defined by is of order unity (see Lai & Pu 2017). As the innermost planet decays in semi-major axis, the system may transition from 12 1 (strong coupling between the inner two planets) to 12 1 (weak coupling), crossing the resonance and generating appreciable mutual inclination θ 12 . If the orbits of m 1 and m 2 are initially co-planar, and m 3 is initially inclined with the inner two planets at an angle θ 3 , then the typical mutual inclination excited is of order ) To illustrate the possibility of resonance, in Fig. 9 we show the final values of the mutual inclination between the inner two planets (θ 12,f ), and the angle between the inner planet's orbit and the spin axis of its host star (θ 1 ,f ) 3 for a USP-forming three-planet system. Note that in the examples shown in Fig. 9, the stellar spin period is P = 30 days, corresponding to a case where the spin-orbit coupling is negligible (ω i is small). We find that indeed, when the system crosses 12 1 during orbital decay, large mutual inclinations can be excited between the innermost planet and its companion. Moreover, this excitation is larger when the ratio m 1,0 /m 2,0 (and thereby L 1 /L 2 ) is smaller: we find that systems with m 1 = M ⊕ and m 2 = 17M ⊕ achieved a maximum value of θ 2 12,f 1/2 ≈ 16 deg., compared to θ 2 12,f 1/2 ≈ 11 deg.
for the case of m 2 = 7M ⊕ .
When there is more substantial spin-orbit coupling, the mutual inclination evolution is similar, except that the resonance occurs at higher values of 12 . To illustrate this, in Fig. 10 we show the same examples as Fig. 9, except with the stellar spin period fixed at P = 1 day, corresponding to strong spin-orbit coupling (large ω i ). In this case, there is still a resonant excitation in the mutual inclination, except that it happens at much larger values of 12 .
Another way to look at the role of spin-orbit coupling is to consider what happens when P changes while fixing the other parameters. In Fig. 11, we show the final value of θ 12,f and θ 1 ,f as a function of the value of P (fixed in time during the evolution) for a three-planet system under-going In general, for systems with 12,0 < 1 and 12, f ∼ 1, the "transition" in Fig. 11 occurs when P reaches a value such that ω 1 ,f ω 12,f at the end of orbital decay. In other words, if ω 1 ,f ω 12,f , then spin-orbit coupling will suppress any resonant mutual inclination excitation between the innermost planets.
In real systems, the stellar rotation period increases over time, thus the importance of the spin-orbit coupling depends on the timescale of the proto-USP orbital decay: if the orbital decay occurs well with-in a Gyr, then spin-orbit coupling can be important. Otherwise, the star would have already spun down by the time the final USP semi-major axis is reached, and the effect of spin-orbit coupling is small.

POPULATION SYNTHESIS MODEL
We synthesize the results of sections 2 -5 by performing a population synthesis calculation of USPs generated through the low-e migration mechanism. Given the inherent uncertainties in various population statistics (of both USPs and larger-period planets), the purpose of this study is not to accurately reproduce all the observed population of USPs. Instead, our goal is to illustrate the statistical trends that would be expected when USPs are generated by low-e migration.
The initial semi-major axis of the inner-most planet a 1 is drawn from a power-law distribution given by in the interval [P min , P max ], with P max = 8 days and P min = 0.5, 1.0, 2.0 and 3.0 days in four separate experiments. The planet's mass m 1 is drawn from a log-uniform distribution between 0.5 and 3.0 M ⊕ . The inner planet's composition is assumed to be Earth-like, with a radius of R 1 = (m 1 /M ⊕ ) 1/4 R ⊕ (Zeng et al. 2016); its tidal lag time ∆t L,1 is chosen to be 1000, 100 or 10s, corresponding to Q 1 = 7, 70 and 700 for P 1 = 1 day. The outer planet masses are drawn from a log-uniform distribution between 3 and 20M ⊕ . We assign these outer planets (i ≥ 2) a rocky compositions with a H/He envelope comprising a few percent of its mass, with radii given by R i = R i,core + R i,env = R ⊕ [(m i /M ⊕ ) 1/4 + 1.5], and tidal lag times ∆t L,i = 1 or 10 sec, corresponding to Q i = 7 × 10 3 , 7 × 10 4 for P i = 10 days. The initial semimajor axis of the outer planets are given by the ratios a 2 /a 1 and a 3 /a 2 , chosen independently on a log-uniform distribution between 1.41 and 3.0, corresponding to period ratios between 1.67 and 5.2. The initial eccentricities of all planets are equal toē, which is chosen from a Rayleigh distribution with scale parameter σ e = 0.10. The initial complex inclinations I j of each planet are chosen from a 2-D Gaussian distribution with mean µ = 0 and variance σ θ =ē/2; the resulting ratio e/θ 2 is consistent with equipartition of random velocities suggested by numerical simulations of accreting planetesimals (Kokubo & Ida 2002). This choice of inclinations is equivalent to a Rayleigh distribution for |I i | with scale parameter equal toē/2 and with the complex argument uniformly distributed between 0 and 2π. We include the effect of tidal decay due to stellar tides, as given by Eq.
(13). The value of Q is chosen to be 10 6 , 10 7 or 10 8 . We adopt an initial stellar spin of P ,0 = 8 days, subject to the Skumanich law (Eq. 62) with a MB = 3.2 × 10 −14 yr −1 such that the stellar spin period lengthens to P = 30 days at t = 5 Gyr. We account for the possibility of dynamically unstable systems. A system of N planets on initially circular orbits is stable up to τ ≡ t/P 1 orbits if the spacing satisfies (a i+1 −a i ) ≥ k c R H , where k c is a parameter that depends on N and log τ, and R H is the mutual Hill radius given by For mildly eccentric systems, the same criterion as above can be applied, but instead of the semi-major axis difference (a i+1 − a i ) one should use the pericenter-apocenter distance (Pu & Wu 2015). We adopt a value of k c ∼ 7 (Smith & Lissauer 2009), applicable for N = 3 and τ ∼ 10 8 (the typical eccentricity damping timescale); when systems fail to meet this stability criterion, they are regarded as potentially dynamically unstable. We find that systems can indeed become potentially dynamically unstable before forming USPs. 17.4% of systems that formed USPs and 6.3% of systems that did not form USPs became dynamically unstable at some point of their evolution; systems that form USPs are more likely to become unstable due to their larger initial eccentricities, so dynamical instability may be an impediment to USP formation, although the effect is minor.
We evolve our systems for 10 Gyr. In some cases, the inner planet's semi-major axis can shrink to a value less than R ; when this occurs, we assume the planet is tidally disrupted and/or engulfed by the star, and we remove it from the system and halt the simulation.
We found that our initial population of planet systems indeed formed USPs during its evolution, with statistical properties similar to the observed population. The USP population show substantial statistical differences with the longer period planets. We summarize their main properties below.
• Final period: Our initial population of 3-planet systems is capable of producing USPs, with the final periods attained being as low as P 1 ∼ 5 hrs. Most notably, we can reproduce both the sudden change in the slope of the period distribution at P 1 ∼ 1 day and the mild excess of planets around P 1 ∼ 1 day in our simulations (see Fig. 12). This trend persists over a moderate range of of planet masses m 1 , and inner planet tidal Q 1 , but depends sensitively on the initial planet eccentricities e 2,0 , e 3,0 (both equal toē; see above) and the stellar tidal Q . Note that in our simulations, the final distribution for P 1 is not the same as the distribution for P (of all the planets), since in some cases P 2 can also be in the range [1,8] days, although this mixing does not affect our conclusions.
• The value of Q : The final period distribution of USPs has a strong dependence on the value of Q . We show a histogram of the initial and final periods for various choices of Q in Fig. 12. For systems with Q = 10 6 , the period distribution of USPs is strongly carved by stellar tides over Gyr timescales, which results in fewer USP planets observed at smaller periods. On the other hand, systems with Q = 10 8 are not strongly affected by stellar tides, resulting in a much larger fraction of planets surviving at smaller periods. Our simulations suggest that for ∆t L,1 = 100s, a value of Q between 10 6 and 10 7 best matches the power-law period distribution of USPs given by Petigura et al. (2018) and Lee & Chiang (2017), and Q 10 8 is incompatible with observations in our scenario.
• The value of Q 1 : The final period distribution of USPs also depends moderately on the inner planet's tidal Q 1 (see Fig. 13). We find that as expected, a larger value of Q 1 leads to fewer USPs: systems with Q 1 = 70 and 7 feature 2 and 3.5 times more USPs respectively than systems with Q 1 = 700.
• Inner planet mass: Less massive inner planets are more likely to become USPs. Over our sample, the inner planet's mass is smaller for USPs, with m 1 = 1.25M ⊕ for USPs versus 1.5M ⊕ for the entire population (see Fig. 14). Systems with 0.5M ⊕ < m 1 < 0.75M ⊕ were 60% more likely to form USPs than systems with 1.75M ⊕ < m 1 < 2.25M ⊕ . This is because USP formation is limited by the amount of angular momentum deficit (section 3.2), and systems with less massive inner planets have an easier time reaching the required amount of AMD.
• Outer planet masses: Conversely, we find that USP production favors systems with more massive outer planets, although this is a weak effect. The average mass for the exterior planets across all samples is 11.5M ⊕ , and 12.1M ⊕ for the subset that ended up producing USPs.
• Initial eccentricities: We find that USP generation is strongly dependent on the initial eccentricities, with the fraction of systems producing USPs roughly doubling for every 0.05 increase inē: systems with 0.075 <ē < 0.125 and 0.125 <ē < 0.175 produce 1.7 and 3.8 times more USPs respectively than systems with 0.025 <ē < 0.075. We show the dependence of the final period distribution for various initial eccentricities in Fig. 15.
• Initial inner period cutoff P min : Our results can be used to constrain the minimum period P min for the initial planet population. We show the dependence of the USP period distribution on P min in Fig. 16. The systems with P min = 0.5 day show little difference compared to those with P min = 1 day, because virtually all planets with initial P 1 ≤ 1 day spiral into their host stars through a combination of planetary and stellar tidal dissipation. In other words, the low-e USP migrationmechanism is not sensitive to planets with initial P 1 day. On the other hand, the results of experiments with P min = 2 or 3 days show a substantial deviation from the P min = 1 day case, and disagrees with the observed period distribution. Thus, planets must be formed in the 1 < P < 3 days range to reproduce the currently observed period distribution of USPs, although we cannot rule out the possibility of planets forming in-situ at P 1 1 day.
• Mutual inclinations: Observationally, USPs show substantially larger mutual inclinations with their closest neighbors compared with typical Kepler multis (Dai et al. 2018). We find that our low-e formation mechanism for USPs naturally generates larger mutual inclinations between the 0.5 1.0 1.99 3.97 10 -2 10 -1 10 0 PDF Q = 10 6 Q = 10 7 Q = 10 8 Initial Figure 12. Histogram of the initial and final period of the innermost planet P 1 for systems with inner planet tidal lag time ∆t L,1 = 100 s. The blue bars shows the initial distribution of P 1 , while the blue, green and red lines show the final USP period distribution for values of Q = 10 6 , 10 7 and 10 8 respectively. The two solid black lines are given by the power-law distribution d N /d log P 1 ∝ P α 1 , where α = 3.0 for P 1 ≤ 1 day and α = 1.5 for 1 < P 1 < 8 days; we also adopt the discontinuous "bump" at P 1 = 1 day corresponding to an excess of 50% more planets just below P = 1 day as proposed by Lee & Chiang (2017). The normalization of the black lines is chosen so that the total probability density integrates to unity. inner planets. We show a histogram of the final RMS mutual inclinations between the inner planets θ 12,f after 10 Gyr of low-e migration in Fig. 17. The final value of θ 12,f for systems that produced USPs is θ 2 12,f 1/2 ≈ 18 deg., which is more than double the value of θ 2 12,f 1/2 ≈ 8 deg. for systems that did not end up producing USPs. A condition for the i-th planet to transit its host star is that the orbital plane be inclined relative to the line of sight by less than arcsin R i +R a i . Using this criterion, we find that 17.5% of USPs had transiting companions, compared with 63.5% of inner planets with P 1 > 1 day. This result is consistent with empirical studies, which found USPs to have a transiting companion fraction of 4 -12 %, compared to 43 − 59% for small planets with 1 ≤ P ≤ 3 days .
• Inner Pair Period Ratio: We find that systems which resulted in USPs have substantially larger period ratios P 2 /P 1 : USP systems have a mean period ratio of P 2 /P 1 = 14, whereas for non-USP systems the mean period ratio is only 3.5. Fig. 18 shows the PDF of the initial and final period ratios. We also find the period ratio P 2 /P 1 increases as P 1 decreases: the mean period ratio is 4.0, 5.2 and 7.0 for P 1 = 3, 2 and 1 day respectively. This is consistent with the observation that USPs and their companions have period ratios ≥ 15, while non-USPs have a broader period ratio range between 1.4 P 2 /P 1 5  Initial Figure 13. Same as Fig. 12 except that we fix the value of Q = 10 7 and instead vary ∆t L,1 = 10, 100, 1000 s, corresponding to tidal Q 1 = 700, 70, and 7 (at P 1 = 1 day), for the blue, green and red lines respectively. Initial Figure 14. Same as Fig. 12 except that we fix the value of Q = 10 7 , ∆t L,1 = 100 s and instead vary m 1 = 0.5 ± 0.25, 1.0 ± 0.25 and 2.0 ± 0.25 M ⊕ , for the blue, green and red lines respectively.

DISCUSSION
We have studied the formation of USPs through low-e tidal dissipation driven by secular forcings from exterior (super-Earth/mini-Neptune) companions of proto-USPs. In this section we evaluate this proposed formation mechanism in light of the observations of USPs and their population statistics. We then discuss some specific USP sources, potential uncertainties and future extensions to this work.

Low-e USP migration and observations
As discussed in section 1, USPs have a number of distinct properties compared to the bulk of longer-period Super-Earth systems (see Winn et al. 2018). Our study shows that our low-e migration scenario produce USPs with the ob- served properties under a variety of initial conditions (see section 6). USPs are preferentially formed from smaller terrestrial planets with more eccentric external companions. A fiducial set of systems, with inner planet masses M ⊕ < m 1 < 3M ⊕ , exterior planet masses 3M ⊕ < m < 20M ⊕ , inner planet ∆t L,1 = 10s (corresponding to Q 1 = 700 at P 1 = 1 day) and Q = 10 7 ended up producing a posterior USP period distribution that qualitatively matched the observed one, without any fine-tuning. Note that this combination of parameters is not the only one that can produce the observed P 1 distribution; there is a hyper-surface of possible initial system parameters that can fit the observations. For example, a set of systems with Q = 3 × 10 6 and ∆t L,1 = 100 s would fit the . PDF of the final initial and final period ratio of the inner planets P 2 /P 1 distribution for systems in our population synthesis. The blue curve is the initial period ratio, while the green and blue curves are the final period ratios for systems that resulted in USPs and no USPs respectively.
observations similarly well. Nevertheless, the fact that our set-up was able to reproduce observations without tuning of parameters, and that similar looking distributions can be obtained when varying the parameters Q 1 , Q , m 1 andē (see Figs. 12 -16) over factors of a few lends us confidence in the robustness of this mechanism. Even more encouragingly, this formation mechanism naturally produces higher mutual inclinations between USPs and their closest companions, a trend which has been observed by empirical studies. Petrovich et al. (2018) found that a mutual inclination of θ 2 12,f 1/2 20 deg is needed to account for the observed dearth of transiting companions to USPs. In our population synthesis model, we found that systems that produced USPs featured a value of θ 2 12,f 1/2 ≈ 18 deg, which more than doubles the amount for systems that did not end up with USPs and is comparable to the value required by empirical studies. USPs formed in our mechanism have a transiting companion fraction of 18% compared with 64% for planets with 1 ≤ P 1 ≤ 3 days, close to empirical values of 4-12% and 43-59% respectively . We also reproduce the observation that the period ratio between USPs and its companions tends to be large (P 2 /P 1 15) and increases with decreasing P 1 (see also Steffen & Farr 2013).
The feasibility of the low-e formation mechanism for USPs hinges mostly on the magnitude of the initial eccentricities of multi-planet systems: for systems with initial proto-USP periods between 1 -3 days, our population model suggests that an initial eccentricity ofē 0.1 is required. In contrast, present-day Kepler multis have typical eccentricities of σ e = 0.05 − 0.08 (Van Eylen & Albrecht 2015; Xie et al. 2016;Van Eylen et al. 2018), although with the caveat that the presently observed planet eccentricities may have suffered damping over Gyrs by tidal dissipation, and their initial values may be larger.
It is interesting to compare our mechanism with the secular chaos "high-eccentricity" mechanism proposed by Petrovich et al. (2018). These two mechanisms require different initial conditions and produce USPs with distinct final configurations. In the Petrovich et al. (2018) scenario, proto-USPs with a 1,0 between 0.05 -0.1 au attain large eccentricities (1 − e 1 1) through to secular interactions with exterior planets; as the proto-USP pericenter reaches ∼ 2R , the planet is tidally captured and eventually circularized, becoming an USP. In contrast, our low-e migration mechanism requires a proto-USP with a 1,0 between 0.02 − 0.04 au, driven to mild eccentricities (e 1 ∼ 0.05 − 0.2) through secular interactions with exterior planets, followed by tidal decay. USPs formed via secular chaos have a smaller ratio a 1,f /a 1,0 , and therefore a larger amount of AMD is needed (see Eq. 48): one typically requires N ≥ 3 exterior planets with m i 10M ⊕ , e i 0.1 and period ratios P i+1 /P i 3.0, and planet companions with a 2 0.2 au are ruled out due to dynamical instability. In contrast, our low-e migration mechanism requires less AMD to succeed: N ≥ 2 exterior planets with m i 3M ⊕ can satisfy the AMD constraint, and there is no need for large values of P i+1 /P i . Observations (e.g Steffen & Hwang 2015) suggest that Kepler multis have typical period ratios 1.4 ≤ P i+1 /P i ≤ 3.0, and do not support the existence of large numbers of sparsely-spaced (i.e. P i+1 /P i 3.0) multi-planet systems. In addition, while Petrovich et al. (2018) did not attempt to reproduce the final period distribution of USPs formed in their scenario, the low-e migration mechanism can robustly reproduce the observed period distribution over a range of initial parameters (see section 6).
Another difference between the high-e and low-e migration is the final distribution of the USP inclination (θ 1 ). For a system under-going secular chaos, whenever e 1 grows to a very large value so too will the value of θ 1 due to the equipartition principle (e.g. Lithwick & Wu 2014). In the absence of strong spin-orbit coupling, Petrovich et al. (2018) found that USPs can often reach very large values of inclination, with potentially a third of systems attaining θ 1 ≥ 30 deg. In contrast, the low-e migration scenario produces USPs with more mild inclinations (θ 1 ∼ 18 deg), although this value is still enhanced relative to non-USPs.

Specific sources
We comment below on the feasibility of USP low-e migration for two well-studied USP systems.

Are USPs photo-evaporated cores of mini-Neptunes?
The observed population of USPs have radii that are mostly within the range 1.0R ⊕ < R < 1.4R ⊕ , and there is a dearth of planets with intermediate radius 2R ⊕ < R < 4R ⊕ with subday periods, despite such planets being ubiquitous amongst Kepler's longer-period planet population. One common explanation for this observation is the scenario that USPs were initially mini-Neptunes that have had their envelopes stripped due to photo-evaporation (e.g. Winn et al. 2017). This picture may be incompatible with our model, and an alternative explanation might be preferred. Our results show that because low-e USP formation is generally AMDlimited (section 4.3), more massive planets are severely disfavored from becoming USPs. This would preclude highermass super-Earths or mini-Neptunes from becoming USPs.
As a result, USP formation is limited to smaller mass planets (m 1 M ⊕ ), which would have a hard time maintaining their atmospheres against various escape mechanisms.
Another factor that can potentially explain the lack of larger-radius USPs is the dichotomy in tidal Q 1 between rocky planets and those with more extended gaseous envelopes. Our results show a reduction in USP formation efficiency by a factor of ∼ 2, when Q 1 is increased by a factor of 10. In the Solar System, values of tidal Q 1 are in the range of 10 − 500 for terrestrial planets and satellites, but planets with substantial gaseous envelopes (such as Jupiter, Saturn, Uranus and Neptune) have values of Q 1 that are hundreds of times larger (Goldreich & Soter 1966;Lainey 2016). If this trend can be extrapolated to exoplanetary systems, then this dichotomy in tidal Q 1 between planets with and without gaseous envelopes can also explain the lack of USPs with R 1 2R ⊕ .

Uncertainties and Future Work
In carrying out this work, we made several simplifications, which may cause additional uncertainties; we discuss them below.
• Effects of Mean Motion Resonance: One source of uncertainty is the role of mean-motion resonance (MMR) in modulating the secular interactions between planets. In our population synthesis model, we considered planet systems with semi-major axes ratios 1.41 ≤ a 2 /a 1 ≤ 3. In many cases, as the inner planet migrates in-wards the system may encounter MMRs (see also Hansen & Murray 2015). A careful study of the effect of MMRs on the secular interactions of multi-planet systems is beyond the scope of this work. MMRs can excite the eccentricities of the planets, independent of secular interactions. One example is Kepler-80, a system containing an USP accompanied by 5 external planets. MacDonald et al. (2016) found that the outer 4 planets of Kepler-80 are interlocked in 4 sets of three-body meanmotion resonances, each with a libration of around a few degrees. The resulting librations may provide the entire system with an additional source of AMD that ameliorates the effect of tidal dissipation. Another possibility is that MMRs can result in 'resonant repulsion', which would cause the semi-major axes of the two planets in resonance to suddenly diverge outside of the MMR (Batygin & Morbidelli 2013;Lithwick & Wu 2012). MMRs could bring about unexpected and interesting interactions in proto-USP systems and deserves to be the subject of further study.
• Secular Chaos and Dynamical Instability: In this work, we adopted a linear theory in the planet eccentricities and inclinations (by assuming e i , θ i 1). In this linear regime, the eccentricity and mutual inclination evolution of the planet orbits are decoupled. In reality, planet systems that produce USPs will often have inner planets with moderately large values of e 1 0.3. Such values would make higher-order terms in e and θ important, and our linear theory would break down. A non-linear coupling between planet eccentricities and inclinations can bring about secular chaos (Lithwick & Wu 2014), which can enhance the inner planet eccentricities even further as AMD diffuses throughout the system.
Another issue is that as planet eccentricities increase, their orbits may become dynamically unstable leading to orbit crossings. In our population study, we found that a small proportion (∼ 17%, Sec. 6) of systems that became USPs may become dynamically unstable. For these inner systems, the dominant final outcome of dynamical instability is physical collision between the two unstable planets. Once two planet have crossing orbits, for large eccentricities and inclinations (i.e. e 1 , θ 12 [(m 1 + m 2 )/3M ] 1/3 ) the timescale to the first physical collision is given by (e.g. Ida & Nakazawa 1989): Since the collisional timescale is much shorter than the eccentricity damping and orbital decay timescale, once two planets cross orbits, they will quickly undergo a physical collision, which can potentially inhibit USP formation. The extent to which these dynamical instabilities occur requires investigations using numerical N-body simulations and is outside the scope of this work.
• Effect of Additional Planets: In this work we have limited our attention to USP formation in systems with 2 or 3 planets. What happens when additional planets are present? Our framework for 3-planet proto-USP systems can be easily generalized to systems with more than 3 planets. In general, the generation of USPs is constrained by the dual criteria that the system must have sufficient AMD (Eq. 46), and the forced eccentricity on the inner planet must be sufficiently large (Eq. 47). In section 4, we found that for planet systems withē 0.1, the AMD criterion is usually more stringent. The presence of additional exterior planets only help to overcome this constraint and bolster the chances of USP generation, since having more outer planets will increase the total reservoir of AMD to maintain the tidal decay of the inner planet. Moreover, the presence of additional planets (and thereby eigenmodes) increases the likelihood of hitting one of eccentricity secular resonances that can speed up the tidal orbital decay timescale. As a result, we expect USP formation in systems with N ≥ 4 planets to be similar to systems with N = 3 planets, albeit at an enhanced rate.

SUMMARY AND CONCLUSION
In this paper we have studied a "low-eccentricity" migration scenario for the formation of USPs. In this scenario, a lowmass (m 1 ∼ M ⊕ ) inner planet with initial period of a few days is accompanied by several external planets in configurations typical of Kepler multi-planetary systems; the companion planets excite and maintain the eccentricity of the innermost planet, which then experiences orbital decay due to tidal dissipation and eventually becomes a USP. Tidal dissipation in the host star further enhances this orbital decay when the inner planet reaches a sufficiently small period. We find that this low-e mechanism naturally produce USPs from the large population of Kepler multis, and can explain most of the observed population properties of USPs. The key findings of this paper are: • We study analytically the condition for orbital decay of the inner planet induced by secular forcing from the outer planetary companions for systems with N = 2 (section 3) or N = 3 (section 4) planets. USP formation is governed by two criteria (section 3.2): (i) the total system angular momentum deficit (AMD) must be sufficiently large, and (ii) the forced eccentricity on the inner planet must be sufficiently large so that decay occurs within the lifetime of the system. We find that it is difficult for 2-planet systems to simultaneously satisfy both criteria due to the suppression of forced eccentricity on the inner planet by short-range forces. On the other hand, 3-planet systems have a much easier time forming USPs (section 4), as the presence of the 3rd planet introduces secular resonances that can boost the inner planet eccentricity, in addition to enhancing the AMD reservoir.
• Although the basic equations (based on secular Laplace-Lagrange theory) for eccentricity excitation and orbital decay in multi-planet systems are standard, in practice they are computationally difficult to evolve for long periods of time due to the "stiffness" of the equations: whereas orbital decay occurs on Gyr timescales, secular interactions proceed on timescales as short as ∼ 100 yr. To resolve this, we develop an approximate method based on the evolution of eigenmodes (section 2.1). We find that eigenmode crossing during orbital decay can lead to secular resonances, which can excite large eccentricities in the inner planet.
• We extend our analysis to the mutual inclination evolution in section 5. We find that secular inclination resonances can also excite mutual inclinations between the innermost planet and its companions as it undergoes tidal decay. Moreover, the range of parameters for which the secular inclination resonance and secular eccentricity resonance occur usually coincides with one another, which results in large mutual inclinations being generated whenever a USP is formed.
• Using our approximate "eigenmode" method, we carry out a large population synthesis study to examine the statistical properties of USPs formed in the low-e migration scenario (section 6). We find that USPs can be robustly produced from typical Kepler multis under a range of initial conditions. This formation mechanism favors smaller inner planets, and requires the initial eccentricities of the companion planets to beē 0.1. We find that the final USP period distribution depends on the values of planet tidal Q 1 and stellar tidal Q ; in particular, a configuration with proto-USP mass M ⊕ < m 1 < 3M ⊕ and tidal lag time ∆t L,1 ∼ 10 − 1000s, outer planet masses 3M ⊕ < m i < 20M ⊕ (i ≥ 2) and Q ∼ 10 6 − 10 7 produces USPs with a final period distribution that matches closely with the observed one.
• Confronting with observations, we find that our lowe migration mechanism can reproduce the empirical population properties of USPs. The final period distribution of USPs matches with the empirical distribution, and the radius distribution of USPs are biased towards small, Earthlike planets, in agreement with observations. Moreover, we find that in our low-e formation mechanism, systems with USPs have more than twice as large mutual inclinations between the innermost planets as do systems without USPs, in agreement with other empirical studies. Our mechanism also reproduces the empirical fraction of USPs with transiting companions, as well as the period ratio distribution of such USPs, without fine tuning of initial parameters.
Overall, we conclude that the low-e migration mechanism can more robustly produce the observed USPs than some of the other proposed mechanism (see section 7). For some systems (e.g. Kepler-10), our scenario makes specific predictions for the existence of unseen planets which can be tested by future observations.