Chaotic Type I Migration in Turbulent Discs

By performing global hydrodynamical simulations of accretion discs with driven turbulence models, we demonstrate that elevated levels of turbulence induce highly stochastic migration torques on low-mass companions embedded in these discs. This scenario applies to planets migrating within gravito-turbulent regions of protoplanetary discs as well as stars and black holes embedded in the outskirts of active galactic nuclei (AGN) accretion discs. When the turbulence level is low, linear Lindblad torques persists in the background of stochastic forces and its accumulative effect can still dominate over relatively long timescales. However, in the presence of very stronger turbulence, classical flow patterns around the companion embedded in the disc are disrupted, leading to significant deviations from the expectations of classical Type I migration theory over arbitrarily long timescales. Our findings suggest that the stochastic nature of turbulent migration can prevent low-mass companions from monotonically settling into universal migration traps within the traditional laminar disc framework, thus reducing the frequency of three-body interactions and hierarchical mergers compared to previously expected. We propose a scaling for the transition mass ratio from classical to chaotic migration $q\propto \alpha_R$, where $\alpha_R$ is the Reynolds viscosity stress parameter, which can be further tested and refined by conducting extensive simulations over the relevant parameter space.


INTRODUCTION
Planetary migration occurs within the early stages of planet formation, when planets tidally interact with their protoplanetary discs (PPDs).Type I migration primarily affects low-mass planets (Goldreich & Tremaine 1980;Ward 1997), typically those up to a few times the mass of Earth, when the companion lacks the mass required to create a gap in the disc density profile (Lin & Papaloizou 1986).It plays an important role in the formation of super Earths and cores of proto-gas-giant planets (Ida & Lin 2008;Benz et al. 2014;Liu et al. 2015).Early analytic analysis (Ward 1997) and numerical simulations (Paardekooper et al. 2010) of this process primarily focus on the determination of tidal torque between the disc and planets with fixed orbits.The pace and direction of this migration are determined by the sum of the differential (usually negative) Lindblad and unsaturated (often positive) corotation torque.The details of torque imbalance is set by the disc structure (Paardekooper et al. 2010;Kley & Nelson 2012), particularly the temperature and surface density gradients.Type I migration usually leads to inward motion and occurs much faster than both the growth of giant planet cores and the evaporation of the disc (Ward 1997;Tanaka et al. 2002).However, depending on the specific heating mechanisms and accretion-disc structures, the total torque can exhibit local positivity in certain regions of the PPD.This in turn can cause stalling of Type ★ yw505@leicester.ac.uk (YW), yc9993@princeton.edu(YC) I migration among a range of planetary mass (Kretke & Lin 2012), forming effective migration traps at the boundary of outward migration regions (Liu et al. 2015).The potentially rapid migration speed warrants self consistent simulations in which the embedded planet's orbit evolve together with the disc structure.Moreover, recent studies (Wafflard-Fernandez & Baruteau 2020;Wu et al. 2023) show that observational signatures generated by migrating planets may differ from those produced by stationary planets.
Stellar-mass Black Holes (sBHs) also migrate and evolve in active galactic nuclei (AGN) accretion discs around Supermassive Black Holes (SMBHs).The hierarchical mergers of sBHs could possibly contribute to LIGO-Virgo gravitational wave (GW) events (McKernan et al. 2014(McKernan et al. , 2018)).Early population studies of these sBHs and their GW features invoke classical Type I migration formulae to account for their orbital decay due to tidal interaction with the disc, driving them to converge radially towards certain migration traps existent in the AGN disc structure (Bellovary et al. 2016) and undergo frequent mergers and interactions (Tagawa et al. 2020).
Turbulence ubiquitously exist in both kinds of discs, albeit at different levels.In mid-planes of PPDs, Magneto-rotational Instability (MRI) turbulence level is proven to be low, with an amplitude  ≲ 10 −3 in the  presciption (Shakura & Sunyaev 1973) in numerical simulations (Bai & Stone 2013), due to low ionization level and non-ideal MHD effects.Such upper limits are also confirmed by observations of molecular line broadening (Flaherty et al. 2017(Flaherty et al. , 2020;;Rosotti 2023).This level of turbulence have been shown to drive stochastic Type I migration of low-mass planets in short-term MHD simulations (Nelson 2005), but it has been pointed out by Baruteau & Lin (2010) that such levels of turbulence is not able to disrupt the classical planetary wake structure, and the residue differential Lindblad torque can still dominate on the long term to drive inward migration.This justifies the use of an effective viscosity  in laminar disc simulations to mimick the effect of turbulent accretion in planet-disc interaction and conclude Type I differential Lindblad torque dependencies, as in the widely applied Paardekooper et al. (2010) formalism.
However, planetesimal or planet core formation might occur very early during PPDs' evolution (Xu 2022;Xu & Armitage 2023;Yamato et al. 2023;Han et al. 2023).During their infancy, PPDs are massive and dense.Marginal gravitation instability leads to amplification of local disturbance and gravito-turbulence induced effective viscosity corresponds to an average effective Reynolds ⟨  ⟩ as large as 0.05 (Gammie 2001;Johnson & Gammie 2003) or 0.1 (Lin & Pringle 1987;Deng et al. 2020).Outskirts of AGN discs are also subject to marginal gravitational instability (Paczynski 1978), gravitoturbulence (Lin et al. 1988), and intense star formation (Goodman 2003;Thompson et al. 2005;Jiang & Goodman 2011;Chen et al. 2023).With a swarm of embedded massive stars co-evolving inside the disc (Cantiello et al. 2021;Ali-Dib & Lin 2023;Huang et al. 2023), they generate density waves which interfere and dissipate, producing an effective turbulent viscosity whose magnitude positively correlates with the mass and number density of embedded objects (Goodman & Rafikov 2001).If strong turbulence is able to significantly alter the flow pattern around the companion, applying linear migration torques in population synthesis of these objects or adding simple stochastic torques on top of a linear component, e.g. in the prediction sBH evolution in AGN discs and GW statistics from their mergers (Secunda et al. 2019;Tagawa et al. 2020), may no longer be accurate.
In this letter, we revisit the chaotic migration problem as posed in Baruteau & Lin (2010), but focus on what occurs in the highly turbulent regime applicable to the fore-mentioned contexts.We emphasise that for ⟨  ⟩ ≳ 0.1, migration of low-mass companions can indeed become chaotic over long timescales as classical circum-companion flow pattern becomes disrupted.The Letter is organised as follows: in §2, we introduce our numerical setup for turbulent companiondisc interaction; We will present our results in §3 and give a brief conclusion in §4.

NUMERICAL SETUP
To explore the effect of strong turbulence on Type I migration, we apply a modified version of the hydrodynamic grid code FARGO (Masset 2000) with a phenomenological turbulence prescription that follows Laughlin et al. (2004) and Baruteau & Lin (2010).
The initial and boundary conditions are similar to that adopted by Baruteau & Lin (2010).The axisymmetric isothermal disc has aspect ratio ℎ = ℎ 0 (/ 0 ) 1/2 , with ℎ 0 = 0.03.The initial surface density is Σ = Σ 0 (/ 0 ) −1/2 , with Σ 0 a normalization constant irrelevant to the dynamics.The embedded companion with a mass of  p , or a mass ratio of  =  p / * compared to its host mass  * , is initially at the orbital radius of  0 = 1, with other code unit being  =  * = 1, the orbital frequency Ω 0 = Ω( 0 ) = 1.Self-gravity is neglected because turbulence is prescribed and not generated, and the instantaneous migration torque that the disc exerts on the planet's orbital angular momentum is normalized in units of Σ 0  4 0 Ω 2 0 .The planet's orbit selfconsistently evolves under the influence of the torque, but since the extent of migration is relatively small within a few thousand orbits, the magnitude of torque isn't affected even if we normalize with Σ 4 Ω 2 at real-time locations.
The simulation domain covers (0.4 − 1.8) 0 in the radial direction, and 0−2 in the azimuthal direction.We resolve the disc by   = 512 radial zones with logarithmic spacing, and   = 1536 azimuthal zones.wave-killing zones are adopted close to the boundaries to minimize unphysical wave reflections (de Val-Borro et al. 2006).To prevent gas velocity from diverging infinitely close to the companion, the planet gravity is softened by a smoothing (Plummer) length of  = 0.6ℎ 0  0 .
To mimic turbulence, a fluctuating potential Φ turb ∝  is applied to the disc, corresponding to the superposition of 50 wave-like modes (Laughlin et al. 2004) such that where  is the dimensionless characteristic amplitude of turbulence.
Each stochastic factor for the -th mode, Λ  , expressed as is associated with a randomly drawn wavenumber   from a logarithmically uniform distribution ranging from  = 1 to the maximum wave number  max .The initial radial   and azimuthal   location of the mode are drawn from uniform distribution.The modes'radial extent is   =   /4.Modes start at time  0, , and their lifetime is Δ  = 0.2  /  , with   being the local sound speed.Ω  is the Keplerian frequency at   , t =  −  0, ,and   is a dimensionless constant drawn from a Gaussian distribution of unit width.
As discussed in Baruteau & Lin (2010), the choice of parameters for this turbulence driver emulates the power spectrum of a typical Kolmogorov cascade, following the  −5/3 scaling law up to  max .An effective Reynold stress parameter ⟨  ⟩ around the companion location can be calibrated from the velocity fluctuations to be approximately i.e. proportional to the square of the turbulence amplitude (see Baruteau & Lin (2010) for details).

Dependence on Turbulence Strength
In Figure 1 we plot the running-time average of migration torques, measured in units of Σ 0  4 0 Ω 2 0 , for a fiducial set of parameters with companion mass ratio  = 5 × 10 −6 and ℎ 0 = 0.03.Different colors represent either runs with different turbulence amplitude  or with different selection of mode azimuthal wavenumbers.
The classical formula (Tanaka et al. 2002;Paardekooper et al. 2010) for type I migration torque scaling is Calibrating with the steady-state  = 0 result, we obtain  L ≈ 1.8 (Figure 1 red line).This quantity is dominated by the Lindblad torque because in the inviscid ( = 0) disc, corotation torque saturates (Goldreich & Tremaine 1980).Although  = 10 −4 brings in extra turbulence and effective viscosity, the running average of the torque nevertheless converges to the value corresponding to  L ≈ 1.8 after a few hundred orbits.By default we include all  up to grid scale, but we tested that truncating higher order modes makes very little difference (Figure 1, green and yellow lines), as the large-scale (low-) modes are most influential to migration.
Such behaviour is expected from the analysis of Baruteau & Lin (2010).When turbulence level is moderate, structure of density waves and flow patterns of classical planet-disc interaction can still manifest after averaging the density distribution, therefore the migration torque can be seen as the superposition of a random fluctuating component with magnitude of and a continuous component close to the classical type-I torque Eqn 4. In a quasi-steady state, the fluctuating amplitude is typically much (by nearly 2 order of magnitude) larger than the continuous component as plotted in Figure 2, while its accumulative effect will decay with time.
With the auto-correlation/eddy-turnover timescale of the turbulent component being ∼ one planetary orbit, the time for the accumulative effect of the linear torque to be significant is  conv ∼ (⟨ turb ⟩/ L ) 2 orbits.From the measured torque data, we estimate the standard deviation of Gaussian-distributed torques to be ⟨ turb ⟩ = 8 and 63 for  = 10 −4 and 10 −3 respectively, consistent with the histogram of total torques (Figure 3) that show a roughly Gaussian dispersion dominated by the fluctuating component with typical dispersion being ∼ ⟨ turb ⟩(/ℎ 0 ) 2 .In our numerical simulations, the running average torque for  = 10 −4 cases quickly converge with the  = 0 scenario within ≪ 100 orbits (Figure 1), corresponding to the expectation of  conv ∼ 20 orbits.While for the  = 10 −3 cases, the average torque quickly decays to values two orders of magnitude below the linear expectation, instead of showing any trend of settling towards the linear prediction.
The analysis that the total torque can be decomposed into a stationary component and a fast-varying component, with no significant coupling between both, assumes that differential Lindblad torque is not significantly altered by turbulence.If we apply the same tactics to the  = 10 −3 scenario, we would expect the timescale for the accumulative effect of the linear torque to be significant to be around 1200 orbits.However, for such strong turbulence corresponding to ⟨  ⟩ ∼ 0.04, structure of density waves will be disrupted even after averaging the density distribution across orbits, similar to the ⟨  ⟩ ∼ 0.1 experiments in Chen & Lin (2023), which led us to believe that the running-average showing no sign of convergence towards the  = 0 value is genuine evidence that random walk dominated this chaotic migration process.In support of this claim, we highlight the torque density distribution of the  = 10 −3 in Figure 4.The top panel represents the torque density of the  = 0 model at 2000 orbits in which the planet has migrated to  ≃ 0.9 0 .In this classical flow pattern picture, the main contribution of Lindblad torque is associated with the competition between the leading waves extending towards the upper left (positive torque) and the trailing waves extending towards the lower right (negative torque) (Tanaka et al. 2002).The middle panel is the snapshot of the  = 10 −3 model in which the planet's orbit has not evolved significantly due to the turbulent interruption of the planet's tidal wake.The bottom panel is a time average (over 200 orbits) torque density (for the  = 10 −3 model) which reduce the fluctuation but still show no sign of the density-wave structure emerging and is generally smaller than magnitude than the inviscid case.

Dependence on Planet Mass
In this section, we further compare the torque evolution of the fiducial case (with  = 5 × 10 −6 ) as well as companion mass ratios  = 10 −6 and  = 3 × 10 −5 .The scale height ℎ 0 = 0.03 is unchanged and all azimuthal wavenumbers are included.In Figure 5 we plot the standard deviation of the torque (normalized by Σ 0  4 0 Ω 2 0 )in these three cases over 500-2000 orbits, fitted by ⟨Γ⟩ ≈ 3.6 × 10 2 Σ 0  4 0 Ω 2 0 , or ⟨ turb ⟩ = 3.6 × 10 2  −1 ℎ 2 0  (6) conforming with Eqn 10 of Baruteau & Lin (2010) with a factor of %50 discrepancy.Figure 6 shows the evolution of the running-average torque for three different mass ratios under strong turbulence, also normalized by Σ 0  4 0 Ω 2 0 .For the low mass ratio  = 10 −6 case, the profile settles towards a fluctuating torque with an amplitude smaller than the linear torque, which is predicted to be ∼ 10 −9 in corresponding units (Σ 0  4 0 Ω 2 0 ), suggestive for chaotic migration similar to the fiducial model.Nevertheless, since |⟨ turb ⟩/ L | is larger than the fiducial case, we expect a robust convergence timescale of ∼ 2 × 10 4 orbits, only after which we can confidently tell if the linear background torque is completely suppressed.We did not run the simulation long enough due to computational constraints, but we comment that for very low masses, both the linear and fluctuating components are expected to be minimal, resulting in an overall small migration effect for all ⟨ R ⟩.
In the high mass case with marginal gap-opening  ∼ ℎ 3 0 = 2.7 × 10 −5 , the fluctuation (∝ ) is larger than the linear torque (∝  2 ) by only one order of magnitude now, suggesting it's much harder to drive chaotic migration of high-mass companions through turbulence (Baruteau & Lin 2010).Naturally we do expect gap-opening planets' migration to be less perturbed by turbulent effects, as shown in MHD simulations (Aoyama & Bai 2023).In Figure 6 we see the average of  = 3 × 10 −5 converging towards ∼ −10 −7 , which suggest that longterm evolution of this gap opening companion will be linear torque dominated.From the specific torque distribution averaged over the last 200 orbits for  = 3 × 10 −5 (Figure 7), we also observe the emergence of density wave structure (trails extending towards the upper left and lower right) in contrast to the fiducial case (Figure 7).Interestingly, the value of converged Γ is still about one order-ofmagnitude smaller than the expected linear torque from Eqn 4. This might be due to partial gap opening or non-linear enhancement of positive corotation torque, which is not the main focus of our study.
Extrapolating to different values of  and  or ⟨  ⟩, we generally expect a critical transition  for a given level of turbulence   above which the planet shifts from chaotic migration towards linear migration.Meanwhile, the mean value of ⟨Γ⟩ begins to become significant.If the marginal gap opening case is close to linear migration.We can determine for ⟨  ⟩ = 0.04, ℎ 0 = 0.03 the transition 5 × 10 −6 ≲  crit ≲ 3 × 10 −5 .The same paradigm can be applied to a subsequent parameter survey to conclude a comprehensive scaling for  crit .
Chen et al. ( 2022) investigated the critical planet orbital eccentricity above which circum-companion flow patterns are disrupted.They concluded that when the characteristic epicyclic impact veloc-    Since Lindblad resonances are effective at radial distances of a few pressure scale height ∼  = ℎ, we might expect strong chaotic fluctuation to significantly impact the type I torque when  ,  is smaller than some factor  of , i.e. when  ≲ ⟨ R ⟩ℎ 3 .
Here 4 <  < 25 is loosely constrained from our limiting test cases.More accurate values need to be obtained by extensive simulations.Such a determination is technically straightforward since the transition from classical to chaotic migration can be distinctively identified by plotting the running time average of torque over long enough timescales.

CONCLUSIONS
In this letter, we use hydrodynamical simulations with a driven turbulence prescription to show that low mass planets, or generally disc-embedded companions migrating in strongly turbulent discs are subject to highly stochastic driving forces, and the time-average residue migration torque becomes negligible after a few thousand orbits.For a  = 5 × 10 −6 companion, when the ⟨  ⟩ ∼ 4 × 10 −4 , the average Lindblad torque and corotation torque will dominate after the decay of the running-average of turbulent torque component, as in Baruteau & Lin (2010), and become subject to linear expectations.However, for stronger turbulence corresponding to ⟨ R ⟩ ∼ 4 × 10 −2 , the classical resonances are disrupted and the residue for runningtime average does not tend to saturate to the classical predictions of Type I torque.This randomizing effect is stronger for lower companion mass, and becomes less significant for marginally gap-opening companion mass.
This stochastic migration paradigm applies to planet migrating in outer PPDs dominated by strong gravito-turbulence, as well as sBHs in AGN discs under influence of gravito-turbulence and the density waves generated by the swarm of embedded objects themselves (Goodman & Rafikov 2001), including intermediate mass sBHs and massive stars (Cantiello et al. 2021).The ⟨  ⟩ ∼ 4 × 10 −2 requirement for our fiducial example appears to be reaching the upper limit of quasi-steady gravito-turbulence (Gammie 2001;Johnson & Gammie 2003;Deng et al. 2017).However, even if intense cooling triggers disc fragmentation, the gas would not be completely deposited into fragments before the disc region in between the fragments re-stablish quasi-steady state at the maximum turbulence level, provided that the fragments formed do not play a major role in the energy and angular momentum budget of the disc.This self-regulation can either be due to a prolonged cooling timescale from lower surface density (more relevant to gas pressure dominated PPDs), or auxiliary heating mechanisms such as large scale spiral shocks or star formation heating (more relevant to AGN discs) (Goodman 2003;Sirko & Goodman 2003;Thompson et al. 2005).Moreover, lower mass companions, whose flow patterns are more susceptible to disruption, would have lower transition turbulence parameters which widens the viable viscosity parameter space for chaotic migration.In light of these considerations, we propose that chaotic migration can indeed operate in a wide range of parameter space relevant to AGN discs, as well as a significant range of parameter space relevant to PPDs.Regarding the dependencies of the magnitude of the random torque ⟨ turb ⟩, our measurement is consistent with the scaling summarised by Baruteau & Lin (2010), which can be recast in terms of ⟨  ⟩ as: which serves as a quantitative formula for future population studies that examine the effect of random driving torque on AGN channel GW statistics.In terms of a reference gravitational stability parameter  0 ≡ ℎ 0 Ω 2 0 /Σ 0 , the net effect of turbulent migration will move the companion radially by a fraction of in either the inward or outward direction over each orbital timescale and its magnitude would increase with time by an amount ∼ (Ω 0 /2) 1/2 .This estimate (Eq.8) is consistent with the numerical simulations (Figures 1 & 4) which show much slower random-walk diffusion (with  = 10 −3 ) rather monotonic migration (with  = 0).
To first order, we expect chaotic migration to prevent sBHs from monotonically migrating towards a migration trap where they become dynamically crowded, with many migrating outwards as well as inwards.Consequently, three-body interactions and hierarchical mergers in migration traps will be less frequent as previously predicted (Tagawa et al. 2020).Although the self gravity of the disc's is neglected in these simulation, the scaling relation in Eqn. 8 indicates that the extend of Δ ≲  0 during the lifespan of massive stars (Ali-Dib & Lin 2023) in typical AGN discs with marginal gravitational stability (Paczynski 1978;Goodman 2003), modest gravitoturbulence (Deng et al. 2020), and small scale height (Starkey et al. 2023).This slow diffusion enables the heavy element pollution of AGN disc (Huang et al. 2023) by massive stars before they migrate into and consumed by the SMBH without the need of universal migration traps (Bellovary et al. 2016).It is worth noting that strong three-body interaction in migration traps has been invoked to explain the low mean effective spin seen in LIGO-Virgo observations (The LIGO Scientific Collaboration et al. 2021), assuming a significant portion of which is indeed from the AGN channel.On the other hand, Li et al. (2022);Chen et al. (2022); Chen & Lin (2023) has proposed other natural ways from the local gas accretion to reduce the magnitude as well as the mean of effective sBH spin from gas accretion, which does not rely on frequent interactions between sBH and sBH binaries.
We also propose a scaling for the transition companion mass ratio from classical to chaotic migration  crit ∝ ⟨ R ⟩ℎ 3 0 , which can be further tested and refined by conducting extensive parameter survey using numerical simulations.The distinct nature of the transition allows for straightforward testing.The parameter survey can be performed using either a comparable code setup or by incorporating more realistic turbulence treatments in subsequent works.

Figure 1 .
Figure 1.Time evolution of the running-time averaged specific torque (normalized by Σ 0  4 0 Ω 2 0 ) for a fixed-orbit companion with a fiducial mass ratio of  = 5 × 10 −6 (∼ 1.6 ⊕ ) compared to a solar mass host starting from 500 orbits.Different colors represented either runs with different turbulence amplitude , or with different selection of mode azimuthal wavenumbers.Discarding small scale turbulence (truncating  > 6 modes) does not introduce significant difference in our results.

Figure 2 .
Figure 2. The total migration torques in units of Σ 0  4 0 Ω 2 0 at around 1000 orbits (upper panel) and 2000 orbits (lower panel) respectively for the fiducial case.20 data points are collected per orbit.

Figure 5 .
Figure 5. Black symbols indicate the standard deviation of the migration torque ⟨Γ⟩ measured over 500-2000 orbits for different values of mass ratio .The values are closely fitted by Eqn 6 (dashed line).The magnitude of the negative linear torque (Eqn 4) with   = 1 is plotted with a solid line.

Figure 6 .
Figure 6.Same as Figure 1 but used to compare different planet-to-star mass ratio  with a fixed  = 1 × 10 −3 .No truncation of turbulent modes are performed.