Young Planets around Young Accreting Stars: I. Migration and Inner Stalling Orbits

Planet formation models suggest that many exoplanets do not form in situ, but rather originate at larger radii in protoplanetary disks and migrate inward due to gravitational interactions with the disk (e.g., Lin & Papaloizou, 1986; Kley & Nelson, 2012). When a low-mass planet (typically less than a few Earth masses) is embedded in a gaseous disk, it undergoes so-called Type I migration, driven by the differential Lindblad and corotation torques exerted by the surrounding disk material (e.g., Tanaka et al., 2002; Tanaka & Ward, 2004). As the planet mass grows, it can carve a partial or full gap in the disk, transitioning into the Type II regime, where migration generally slows and is coupled to the viscous evolution of the disk (Lin & Papaloizou, 1986).

ZSC24 combined Type I and Type II migration timescale into one equation. Specifically, for a planet of mass ratio $q=M_{p}/M_{\star}$ at orbital radius $r_{p}$ in a disk with surface density $\Sigma$, aspect ratio $h=H/r$, and dimensionless viscosity parameter $\alpha$, the migration timescale (see Equation 41 in ZSC24, Dempsey et al. 2020) is given by

where $K\equiv q^{2}/(\alpha h^{5})$, $\Sigma$ is the background disk surface density without the planet, and all quantities are evaluated at the planet’s ___location.

In this formulation, when $K\lesssim 1/h$, the migration rate reduces to the standard Type I regime. Massive planets enhance $K$ which puts them in the Type II regime and slows their migration compared to Type I. When a gap is opened, a higher $\alpha$ value enhances the disk’s radial accretion velocity, accelerating Type II migration (e.g., Crida et al., 2006).

In addition to Type I or II migration, magnetorotational instability (MRI) turbulence introduces stochastic torques on embedded planets (Nelson & Papaloizou, 2004; Baruteau & Lin, 2010). These random fluctuations in the disk density and velocity field impart a random-walk on the planet’s semi-major axis. Earlier approaches (e.g., Laughlin et al., 2004; Johnson et al., 2006) assume the correlation time $\tau_{c}$, the typical timescale of fluctuations, to be approximately half an orbital period. With a circular orbit assumption, the correlation time is thus:

where $J$ is the orbital angular momentum $J=M_{p}(GM_{\star}r_{p})^{1/2}$ of the planet.

For a given torque $\Gamma$, whose fluctuating component is $\delta\Gamma$, the variance of the torque’s fluctuations is taken to be

where $C_{D}$ is a dimensionless coefficient depending on the strength of the turbulence; $\Gamma_{0}=2\pi G\Sigma M_{p}r_{p}$ is the natural torque scaling calculated using the force that a planet would feel if suspended above the disk.

Combining the torque’s variance and the correlation time $\tau_{c}$, we arrive to a diffusion coefficient, which describes the efficiency and rate at which angular momentum is randomly exchanged within the disk due to turbulent motions:

The corresponding diffusion timescale $t_{\mathrm{diff}}$ measures the characteristic time over which stochastic processes significantly alter the planet’s orbital angular momentum. With $t_{\mathrm{diff}}=J^{2}/D$, we have

After plugging in Equation 2, we have

A notable property of the diffusion timescale is that, when compared to the Type I/II timescale (Equation 1) and tidal migration timescale, the turbulent diffusion timescale does not depend on the planet’s mass. However, the values of both $\tau_{c}$ and $C_{D}$ remain poorly constrained by theory alone. Different works have adopted simplified assumptions or order-of-magnitude estimates, with both Johnson et al. (2006) and Adams & Bloch (2009) taking $\tau_{c}=\pi/\Omega$ and assuming $C_{D}$ to be independent of radius and on the order of $C_{D}\sim 0.05$. One of our goals in this work is, through analyzing the torque time series in the simulation, we could measure its autocorrelation function and thereby determine $\tau_{c}$, and by comparing the torque amplitude to $\Gamma_{0}$, we determine $C_{D}$ at different radial positions.

To proceed, we assume an exponential decay for the torque autocorrelation function, consistent with some numerical studies of MRI turbulence (Fromang & Papaloizou, 2006; Oishi et al., 2007). Following Rein & Papaloizou (2009), if $F_{\phi}(t)$ is the azimuthal force per unit mass acting on the planet, we write:

where $\langle F_{\phi}^{2}\rangle$ is root-mean-square (RMS) value, and $g(\Delta t)$ is the autocorrelation function with

and

Once the RMS value and $\tau_{c}$ have been determined, the observed stochastic nature of the torque can be modeled by using a simple one-dimensional Markov chain (Kasdin, 1995; Rein & Papaloizou, 2009) (Section 4.2).

Close to the star, stellar magnetic fields are so strong that the magnetic pressure becomes larger than the accretion ram pressure inside the Alfvén radius $R_{A}$, defined as

where $B_{\star}$ is the stellar magnetic field at the equator, $R_{\star}$ is the stellar radius, and $\dot{M}$ is the accretion rate. In ZSC24, $R_{A}\sim R_{T}$, where $R_{T}$ is the magnetospheric truncation radius, which is defined as the radius where the disk’s azimuthal velocity drops to half the Keplerian velocity $v_{K}$. It also corresponds to where the disk density changes sharply. $R_{0}$ is the code length unit, which is quite close to $R_{T}$ at the end of the simulation. The stellar radius, denoted as $R_{\star}$, is chosen as $0.1~{}R_{0}$, so that the magnetospheric truncation radius is roughly 10 times larger than it.

The initial disk has a midplane density profile as

with $\rho_{0}\equiv\rho_{d}(r=R_{0},z=0)=1$. The time unit $1/\Omega(r=R_{0})=T_{0}/2\pi$, in which $T_{0}$ is the orbital period at $R_{0}$.

The temperature is set to be constant at each cylindrical radius:

where $c_{s}=\sqrt{p/\rho}$ is the isothermal sound speed. By choosing $p=-2.25$ and $q=-1/2$ we get a disk surface density $\Sigma\propto r^{-1}$. For the scale height $H=c_{s}/\Omega_{K}$, we set $(H/r)_{r=R_{0}}=0.1$.

We adopt a locally isothermal equation of state using a fast cooling. More details for this treatment can be found in Zhu et al. (2015).

For the central object, we introduce a stellar structure with a gravitational pull outside the radius $R_{\star}$, and a fixed density, velocity structure and pressure when at $r<R_{\star}$. The star is set to be non-rotating. The magnetic field is set to be a dipole. In Athena++ the vacuum permeability constant is set to 1 so that the magnetic pressure is simply $B^{2}/2$ in code units. Values are chosen so that the initial plasma $\beta=2P/B^{2}$ at $r=R_{0}$ is 10. A density floor that varies with position is introduced in the highly magnetized magnetosphere to avoid small timesteps. For further details see Section 3.2 of ZSC24.

Although our simulation adopts a non-rotating star, we will also discuss planetary migration around rotating stars in Section 5. For a rotating star, we can define a corotation radius in the disk, $R_{C}$, where the disk’s Keplerian frequency matches the star’s spin rate. For stars in the spin equilibrium state, simulations by Zhu (2025) find $R_{T}=0.79~{}R_{C}$ which places $R_{C}$ a bit outside the magnetospheric cavity.

For the rest of the paper, we drop the code unit (e.g., $\rho_{0}$, $R_{0}$, $T_{0}$) in expressions for ease of reading.

Planets are introduced after the disk has reached a quasi-steady accretion stage. We have chosen $t=65.0$ or the 65th orbit at $r=R_{0}$ in ZSC24 as our starting time. The state of the disk just before the planet is introduced is shown in Figure 1. As in ZSC24’s simulations, within the magnetosphere region of $r<1$, disk material becomes filamentary, with gas being clumped together into high density accretion columns or ”fingers” due to interchange instability. While some columns follow the magnetic field lines, being lifted up from the midplane and falling onto the star at about $30^{\circ}$ from the poles, some could penetrate deep into the magnetosphere at the midplane. Thus, while the region inside the magnetospheric boundary has very low density and is fully magnetically dominated $\beta<<1$, it is not completely empty. When the star rotates sufficiently fast (even faster than the equilibrium spin state), the interchange instability will be suppressed and the cavity will be clean (Zhu, 2025).

Planets are added into the disk as orbiting gravitational potentials. The mass ratio between the planet and the central star is $q$. For the planet’s gravitational potential, we apply a smoothing radius $R_{s}=0.012$. This smoothing length approximately corresponds to Jupiter’s radius considering that the stellar surface is at $R_{*}$=0.1. This value is small enough to accurately capture the dynamical interactions between the planet and the nearby disk material.

For our simulations, we consider two distinct cases:

• •

  Case A: Massive planet with $q=0.01$, equivalent to $10~{}M_{J}$, located at $r_{p}=1$.

• •

  Case B: Three smaller planets with $q=10^{-4}$, each roughly equivalent to $30~{}M_{\oplus}$, located at $r_{p}=0.45$, $0.9$ and $1.8$, respectively.

In all these cases, the planets follow fixed circular orbits at the local Keplerian velocity $v_{K}$.

In case A we have a very massive planet in the system which may be able to moderate the accretion from the disk towards the star; a process which could potentially be detectable via observations of variable near-UV excess from accretion shocks or variable $H_{\alpha}$ line as the disk material undergoes magnetospheric accretion (e.g., Hartmann et al., 2016). Meanwhile, the low-mass planets in Case B will not affect the disk structure, allowing us to study the interaction between the disk and multiple planets at different radii simultaneously in one simulation. We thus put three planets inside the magnetosphere, close to the truncation radius, and embedded in the disk, respectively to study the planetary torque and migration behavior for the same type of planet at different starting conditions.

To measure the planetary torque, we exclude the mass within the planet’s Hill sphere as it is gravitationally bound to the planet and can be considered as part of it (e.g., Benítez-Llambay et al., 2015). Without sufficient resolution, the circumplanetary region is not well resolved. Due to its proximity, any structure in this region leads to considerable fluctuations when calculating the torque. For the same reason we remove the innermost cells around the star, excising the volume within 2.5 $R_{\odot}$. During each time snapshot, we measure the contributions of each cell toward the net force vector felt by the planet. Then, we transform the net force’s Cartesian components to the components in spherical coordinates and get the azimuthal component $F_{\phi}$: the planetary torque which is a main driver for planetary migration (Kley & Nelson, 2012; Lesur, 2021).

The planetary torque represents the gravitational torque exerted on a planet by the surrounding protoplanetary disk. It dictates the angular momentum exchange between planets and disks, which leads to the changes in the planet’s orbital radius over time (Papaloizou & Terquem, 2006).¹¹1We have neglected the magnetic torque which will be studied in future. During our simulations, we measure the net torque acting on each planet with time and label it as $\Gamma$. Since our planets are on circular orbits with no inclination or eccentricity, only the azimuthal (tangential) component of the gravitational force, $F_{\phi}$, would contribute to changes in the planet’s angular momentum. The resulting torque is then normalized by the natural scale of the turbulent torque $\Gamma_{0}$ as in Equation 3 (e.g., Johnson et al., 2006; Rein & Papaloizou, 2009). $\Gamma_{0}$ uses the initial azimuthally averaged surface density at each ___location.

Figure 3 shows the normalized torques $\Gamma/\Gamma_{0}$, felt by the single large planet of case A and the three planets of case B, along with their corresponding root-mean-square (RMS) value, which quantifies the typical magnitude of fluctuations. For the $q=0.01$ planet, the RMS value is close to unity - about four times larger than the values measured for the low-mass $q=10^{-4}$ planets. If the turbulent disk is not affected by the planet, the normalized torque should be independent from the planet mass. Thus, the larger RMS value for the massive planet indicates that the massive planet has some feedback to the background turbulent disks. One possible feedback is gap opening by the massive planet (e.g., Oishi et al., 2007). However, with the strong turbulence at the inner disk, a gap has not been induced by the planet. We suspect that the feedback arises from the magnetic interaction between the circumplanetary region and the disk, which deserves further study in future. Nevertheless, the two orders of magnitude change in planet mass only leads to a factor of 4 difference in RMS.

In addition to the stochastic torque, the planet should also feel a net Type I/II migration torque due to its interaction with the launched spirals (see Zhu et al., 2015). Thus, we average the torque profiles to derive the mean torques felt by the planets (shown in Figure 3). If we normalize the mean torque with the typical Type I torque $\Sigma_{0}r_{p}^{4}\Omega_{p}^{2}q^{2}(c_{0}/v_{p})^{-2}$ where $c_{0}$ and $v_{p}$ are the sound speed and the planet’s orbital speed, the massive planet ($q=0.01$) at $r_{p}$=1 registers an average normalized torque of $-0.92$, indicating a gentle inward pull. For the $q=10^{-4}$ planets, the innermost one yields $+0.57$, the one at $r_{p}=0.9$ shows a very strong outward torque of $+17.12$, and the outermost planet registers $-116.3$. Although this implies that the zero torque radius is between 0.9 and 1 which is consistent with trapping by the truncation radius, these averaged torques are strongly affected by the statistical error, not reflecting the true mean Type I torque. The ratio between the stochastic torque ($\sim 2\pi G\Sigma M_{p}r_{p}$) and the Type I torque ($\sim\Sigma r_{p}^{4}\Omega_{p}^{2}q^{2}h^{-2}$) is on the order of $2\pi h^{2}/q$. With $h=0.1$, the Type I torque is smaller than the stochastic torque by a factor of $\sim$ 6 for a $q=0.01$ planet, and a factor of $\sim$ 600 for a $q=10^{-4}$ planet. To derive the small mean torque from the larger stochastic torque, we thus need 36 and 3.6$\times 10^{5}$ independent torque measurements. Since the stochastic torque is correlated within the local orbital time (discussed below), the required simulation duration, especially for the low mass planet or planets further from the star, is much longer than our simulation time.

To analyze the temporal behavior of the stochastic torque, we calculate the autocorrelation function (Equation 7) for the torque time series. This analysis is carried out by comparing the torque dataset with a lagged version of itself. For a given lag, the data is shifted, and the pointwise product of the original and shifted datasets is computed. The time average of these products is then normalized by the variance of the dataset, ensuring the autocorrelation values range between -1 and 1. This process is repeated for different lag values to determine the autocorrelation across various timescales, called the autocorrelation function (ACRF).

The ACRFs of each planet, along with their corresponding exponential fits, as defined in Equation 8, are presented in Figure 4. Each graph also includes an ACRF from our synthetic torque and the corresponding exponential fit, which will be addressed in the following subsection.

In each case, the correlation time, $\tau_{c}$ (written on each graph and also in Table 1), is found to be on the order of $\approx 0.1\mathrm{-}0.2$ times the local orbital period. This aligns with the rapid torque fluctuations experienced by planets embedded in a highly turbulent flow, driven by MRI (see Figure 3), and thus $\tau_{c}$ also represents the eddy turnover time of MRI turbulence.

One new feature we notice in the autocorrelation function is that it shows high positive peaks that appear approximately every local orbit, suggesting that the torque felt by each planet has a periodic component or ”tug” that corresponds to the orbital period. Such a correlation timescale is absent in local shearing box MRI simulations (e.g., Oishi et al., 2007). This implies that some relatively long-lived global disk structures persist over a few orbits, imprinting a characteristic periodicity on the torque felt by the planet each orbit.

To further explore where most torque comes from, we present the planetary torque and its autocorrelation like those shown on Figures 3 and 4 but for individual segments of the disk. The disk has been divided into three parts: an inner, a middle, and an outer part; they correspond to regions at $r\leq 0.8$, $0.8<r\leq 1.5$, and $r>1.5$, respectively. We take the measurements of the torque arising from each part and compare it to the full-disk measurement, while also calculating the ACRFs. In Figure 5 we show these torque contributions as felt by the second planet in case B ($r_{p}=0.9$, $q=10^{-4}$), and their respective ACRFs. This analysis demonstrates that, for this planet located close to the truncation boundary, most of the torque is from the middle and outer regions of the disk, while only a tiny portion from the inner, mass-depleted, magnetosphere cavity, despite the presence of the dense accretion columns. However, regardless of the magnitude of the torque from each part, the periodic behavior previously observed in the full-disk analysis clearly arises from the outer disk, with its ACRF’s peaks being of the highest magnitude, and appearing clearly on each orbit. The same analysis for case B’s other planets, $r_{p}=0.45$ and $r_{p}=1.8$, can be found in Appendix §A.

As we will see in the next subsection, the periodic component of the torque is small in magnitude compared to the stochastic component. However, its periodic character is reminiscent of the time-varying perturbations experienced by circumbinary planets (e.g., Doolin & Blundell, 2011) or stars in barred galaxies (e.g., Contopoulos & Papayannopoulos, 1980), where the non-axisymmetric gravitational potential induces periodic modulations in orbital parameters such as eccentricity and precession. Based on the midplane density contours in Figure 1, the periodic torque is likely from the eccentric cavity or asymmetric disk features due to magnetic bubbles (Zhu et al., 2024). An example of one such bubble’s formation and evolution is shown in Figure 6.

Based on the measured ACRFs (color curves in Figure 4), we assume that the planetary torque can be composed of both a stochastic and a periodic component. Concerning the former, we first build a discrete first order Markov process using $g(\Delta t)$ to generate a stochastic torque with the correlation time $\tau_{c}$ (see Kasdin, 1995; Rein & Papaloizou, 2009), in which the torque evolution at the next step, for each planet, only depends on the local RMS value, correlation time (in units of local orbital time, $T_{p}$), and the present torque. Torques at previous steps are not stored, which makes the method “memoryless”, and the mean value is zero.

The first order Markov chain is thus built from

with

where $\mathrm{d}t$ is the time interval between the $k$ and $k+1$ time step, $k$ is an integer, and $r_{k}$ is a random variable following a normal distribution

The first term on the right hand side of Equation 14 corresponds to an exponential decay of $X$ with a timescale of $\tau_{c}$, while the second term corresponds to the excitation of new white noise (Rein & Papaloizou, 2009).

Yet, this stochastic model lacks the periodic perturbation and its ACRFs thus behave differently to those in the simulations (shown in Figure 4). Therefore, we conjure a deterministic sine function that serves as the quasi-periodic term and add it to Equation 14. The new equation ends up looking like:

where $\xi$, $\psi$, and $\delta$ are amplitude, frequency (cycles per unit time), and phase-shift coefficients of the sine function.

We can thus generate the planetary torque using Equation 16. To find the best fit parameters, we produce a set of synthetic torque time series with the amplitude $\xi$ between $10^{-3}$ and $0.2$, the frequency $\psi$ between $0.5$ and $1.5$, and the phase-shift $\delta$ between 0 and $2\pi$. It’s important to note that having a large amplitude $\xi\gtrsim 1$ will produce ACRFs that stay at high values for a long time before decaying, and having no phase shift will produce peaks that won’t be able to align themselves with those from the simulation data.

From the produced set of possible synthetic torque time series, we calculate their respective ACRFs. We then employ the mean-square-error (MSE) method to systematically compare each synthetic autocorrelation function with its counterpart derived from the simulation; in other words, we compute the squared difference between both curves at each time lag, then select the parameter set $(\xi,\,\psi,\,\delta)$ for Equation 16 that minimizes the total MSE. This provides a quantitative “best-fit” criterion for matching the synthetic and simulated autocorrelations. The selected synthetic ACRFs are those shown in black curves along with the originals in Figure 4.

Additionally, in Figure 7 we show the synthetic torque time series that correspond to the best-fitted ACRFs. The sinusoidal component of the synthetic torque is also overdrawn and enlarged.

For the particular iterations shown in this work, the RMS value $\sqrt{\langle(\Gamma/\Gamma_{0})^{2}\rangle}$ and correlation time $\tau_{c}$ for each of the simulated planets’ torque time series, along with the parameter set $(\xi,\,\psi,\,\delta)$ needed to reproduce their corresponding ”best-fit” synthetic torque time series and respective ACRFs are shown in Table 1. Planet A.1 is the single planet in case A, while B.1, B.2 and B.3 are the three planets in case B (from closest to farthest to the star).

We now compare the measured planetary torque variance with the normalization scale $\Gamma_{0}$, and use Equation 3 to obtain average values of $C_{D}$ to be $0.670$ for case A and $0.551$, $0.258$ and $0.066$ for each of the planets (from closest to farthest) for case B. These values with respect to radius are plotted in Figure 8. From this result, we presume a radial and planet mass dependence for $C_{D}$. The mass dependence results from the feedback from the planet to the disk, and this dependence is relatively weak. Using the smaller planets’ values, two radial relations for $C_{D}$ are proposed, for both inside $R_{T}$ and outside it:

which are drawn in Figure 8.

Taking the simulated disk’s vertically integrated surface density, its aspect ratio radial profile $h(r=R_{0})=0.1$, and the simulation’s typical accretion rate of $\dot{M}=-0.005$ (see ZSC24’s section 5.4), we utilize $\dot{M}=3\pi\nu\Sigma$ and $\nu=\alpha c_{s}H$, to estimate the simulation’s $\alpha$ parameter. We find that, for the data past the truncation radius, the profile

represents the data well. The simulation’s $\Sigma$ and both the $\alpha$ (Figure 18) and complementary $\Sigma$ profiles ($\Sigma=\dot{M}/3\pi\nu$) are shown in Figure 9.

These relationships permit us to propose an $\alpha$ dependence for $C_{D}$ outside of the magnetosphere region. Combining Equations 17 and 18, we get

which further links the disk’s local turbulent parameter to turbulent torque in planetary migration.

Additionally, by plugging $\Sigma=\dot{M}/3\pi\nu$, $C_{D}$ (Equation 19), and $\tau_{c}=0.1~{}T_{p}$ into Equation 5, we derive the migration timescale due to turbulent diffusion

where $h_{p}$ is the disk aspect ratio ($H/r_{p}$ or $c_{s}/v_{\phi}$) at the planet ___location. Please note that the turbulent diffusion timescale has no explicit dependence on $\alpha$. If $h\propto r^{1/4}$ as in our simulation, the diffusion timescale decreases quickly with respect to the planet’s distance to the star ($\propto r_{p}^{-1/2}$ ) when the planet is near the magnetospheric truncation radius.

Finally, Figure 10 presents the radial profiles of the surface density fluctuations along the azimuthal direction in the simulation, $\delta\Sigma/\Sigma$, with scaled $C_{D}$ values, where $\delta\Sigma$ is the each annulus’ surface density’s standard deviation. This correlation between $\delta\Sigma/\Sigma$ and $C_{D}$, alongside the previous analysis, gives further evidence that regions with greater surface density perturbations, driven by disk turbulence, lead to enhanced stochastic torque variations (e.g., Nelson & Papaloizou, 2004; Johnson et al., 2006).

By scaling our simulation to a realistic protoplanetary disk that undergoes magnetospheric accretion, we are able to calculate the migration timescale of bodies embedded in the disk through various migration mechanisms. We consider a typical protoplanetary disk with an accretion rate of $\dot{M}=10^{-8}$ $M_{\odot}$ yr^-1 around a Solar analogue (1 $R_{\odot}$, 1 $M_{\odot}$), possessing a 1 kG dipole magnetic field.

With the code length unit $R_{0}$ being set to be the truncation radius $R_{T}=10~{}R_{\odot}$, $R_{0}$ is thus 0.047 au. The time unit $1/\Omega_{0}$ is thus 0.0016 yr. By equating the steady accretion rate of the disk $\dot{M}=-0.005$ in code units to $\dot{M}=10^{-8}$ $M_{\odot}$ yr^-1, the mass unit becomes $3.19\times 10^{-9}$ $M_{\odot}$. The surface density unit is then 13.1 g cm^-2.

Figure 9 thus suggests that $\Sigma\sim 0.2$ g/cm² at $r\sim 0.2$ au. This low surface density is due to the large $\alpha$ value from the strong turbulence. As concluded by ZSC24, even with more realistic aspect ratios making the surface density $\sim 10$ times larger and dust being a small percentage of the total mass, the resulting mass of solids is orders of magnitude less than what is required to explain the abundant exoplanets close to the star. The MRI-active inner disk appears to be highly unfavorable for planet formation, indicating that planets are more likely to form in the more massive deadzone.

To study the migration timescales, we follow Zhu (2025) to assume a realistic disk aspect ratio as $h=0.05(r/R_{T})^{1/4}$. Additionally, for a more realistic disk that is in the spin equilibrium state, we adopt

as shown in Zhu (2025). Thus

which is coincidentally independent from the disk radius.

On the other hand, for any planet that somehow manages to go through the inner MRI turbulent disk and enter into the magnetosphere, the planet cannot interact with the disk anymore (Equation 1 becomes invalid). Instead, the planet would be subjected to strong aerodynamic drag and continue its migration towards the central star if the planet is inside the corotation radius. From ZSC24, the planet’s migration rate through aerodynamic drag inside the inner cavity is governed by

where $R_{p}$ and $\rho_{p}$ are the planet’s radius and its material density (fixed to 3 g cm^-3), respectively. $v_{\mathrm{rel}}=v_{K}-\Omega_{*}r_{p}$ is the relative velocity between the planet and the background magnetosphere plasma with $\Omega_{*}$ being the star’s spin rate. Assuming that the star is in the spin-equilibrium state with the disk, we have $R_{T}=0.79R_{C}$ (Zhu, 2025), where $R_{C}$ is the corotation radius. With $R_{C}$ given, we can compute $v_{\mathrm{rel}}$. This is different from ZSC24, where the non-rotating star leads to $v_{\mathrm{rel}}\sim v_{K}$.

With this disk model, we calculate the migration timescale for planets at different distances from the star through Type I/II migration (Equation 1) or aerodynamic drag (Equation 23), as illustrated in Figure 11. Given that the lifetime of protoplanetary disks is estimated to be between $<1$ and 40 Myr (see e.g., Pfalzner & Dincer, 2024), we find that for most low-mass planets the disk-driven migration timescale exceeds the disk’s lifetime. Only the massive planets are able to reach the magnetospheric truncation radius, with planets of mass ratio $q\approx 0.003$ being the most favored. This implies that, within the inner MRI-active region, significant migration is unlikely for most planets except at the very early times when the accretion rate or surface density is much higher. Once a planet reaches the truncation radius, it enters the aerodynamic drag regime, where its migration timescale increases by several orders of magnitude. With no additional torques, the massive planets become effectively trapped at the truncation radius. For gas giants, this truncation radius may represent the innermost limit of their migration (see e.g., Mendigutía et al., 2024).

Since low mass planets cannot migrate in the inner MRI active disk, they are likely to stay at the deadzone inner boundary (DZIB). Before migrating planets even enter the MRI-active region (active when $T\gtrsim 1000$ K), they must first pass through the deadzone inner boundary. The DZIB marks the transition between the inner, sufficiently ionized, highly turbulent “MRI-active” zone and the outer, much more massive, low-ionization, less turbulent “MRI-dead” zone. This abrupt transition can create a jump in the surface density, which, in turn, generates a positive corotation torque. Such a torque may be strong enough to act as a pebble trap (see e.g., Chatterjee & Tan, 2013) and potentially even function as a planet trap, particularly affecting intermediate-mass planets in the range of $q\sim 10^{-5}\text{---}10^{-4}$ (e.g., Hu et al., 2016; Faure & Nelson, 2016). Thus, low mass planets could form or be trapped at the DZIB. For T Tauri stars, the DZIB is located at $\sim$0.1 au (e.g., D’Alessio et al., 1998; Flock et al., 2019), which is slightly larger than the magnetospheric truncation radius. For Herbig Ae/Be stars, the DZIB can extend to 1 au (e.g., Dullemond & Monnier, 2010), which is significantly larger than the truncation radius, especially considering Herbig Ae/Be stars have much weaker magnetic fields (see e.g., Järvinen et al., 2019). This implies that low mass planets around Herbig Ae/Be stars should be further out in the disk.

Our ACRF analysis in Section 4 suggests that the correlation timescale $\tau_{c}$ is shorter than estimates in other works (see e.g., Johnson et al., 2006; Adams & Bloch, 2009), where values of half or a full-orbit are suggested; instead, we measure values in order of $\sim 0.1\>T_{p}$. This means that Equation 6 needs to be increased by a factor of 5 to:

We compare the turbulent diffusion timescale of Equation 24 with that of gas-driven Type I/II migration (Equation 1) and aerodynamic drag (Equation 23) for two different planet-star mass ratios, $q=0.003$ and $q=10^{-5}$, in Figure 12. These two $q$ values are representative of relatively quick and slow Type I migration rates, respectively. As for $C_{D}$, instead of a fixed value we utilize the radial profile of Equation 19. We find that, regardless of the precise value of $C_{D}$ or $\tau_{c}$ utilized, the turbulent diffusion timescale ends up being several orders of magnitude larger than that of gas-driven migration, being larger than the Hubble time in the whole region of interest. This is due to the small disk scale height and the small disk surface density from large $\alpha$. The ratio between turbulent diffusion timescale (Equation 5) and the Type I migration timescale (Equation 1) is

With $M_{p}=M_{\oplus}$, $\Sigma\sim 2$ g/cm², $C_{D}\sim 0.02$, $h\sim 0.05$ at 0.1 au, and $\tau_{c}=0.1T_{p}$, this ratio is 5$\times 10^{7}$, indicating turbulent diffusion is negligible for planets at the inner MRI turbulent disk.

As we move outwards, from the MRI active region into the deadzone, the ionization fraction decreases, prompting a decoupling of the magnetic fields from the gas. This leads to low or negligible turbulence, and the $\alpha$ parameter could settle to a much smaller value (e.g. $\alpha\sim 10^{-4}$). Turbulence will then be constrained to the active layers above and below the midplane (e.g., Oishi et al., 2007). This would result on the turbulent diffusion timescale at the midplane staying large or even increasing. Therefore, we can generally presume that stochastic torques have a negligible effect on changing a planet’s semi-major axis within several au.

In addition to the Type I/II migration and turbulent diffusion driven by planet-disk interaction, a planet can also migrate through tidal interaction with the central star. Figure 12 also displays the corresponding timescales for tidal and magnetic torques. Following Wei & Lin (2024), we use a standard tidal time-lag formalism (see e.g., Hut, 1981; Eggleton et al., 1998; Ogilvie & Lin, 2004, 2007) to express the torques due to star-on-planet and planet-on-star interactions

where the tidal time lags $\tau_{\star}$ and $\tau_{p}$ represent the delay between the tidal forcing and the tidal response in the star and planet, respectively.

Owing to the substantial mass difference between the star and the planets, and their short orbital separations, planets would typically synchronize their spin ($\omega_{p}\approx\Omega_{p}$) on timescales of a few million years. Consequently, the star-on-planet tidal torque in Equation 27 quickly becomes negligible before significant migration takes place, leaving the planet-on-star torque of Equation 26 dominant.

The tidal time lag is directly related to the tidal quality factor ($Q$), defined as $Q\simeq(\omega\tau)^{-1}$, where $\omega$ is the tidal frequency. For close-in planets, a representative tidal frequency is $\omega\sim 10^{-5}\mathrm{s}^{-1}$. Higher values of $Q$ indicate less efficient tidal dissipation and slower orbital evolution, while lower values imply stronger dissipation and more rapid migration. T Tauri stars generally have $Q\sim 10^{4}$, but as they evolve onto the main sequence, becoming less convective and slowing their rotation, $Q$-values increase significantly, diminishing tidal effects (Gallet et al., 2017). We choose a stellar tidal time lag of $\tau_{\star}=1s$, which is equivalent to the stellar tidal quality factor of $Q=10^{5}$ for 7 days periods.