Halting Migration in Magnetospherically Sculpted Protoplanetary Disks

Our model for how the stellar magnetic field interacts with the protoplanetary disk is based on the work of Ghosh & Lamb (1979), previously applied to accreting compact objects. In this model, the disk has 3 distinct regions: the inner truncation, the transition region and the outer disk. A diagram of the overall disk structure can be found in fig 1. The three regions are determined based on the amount of stellar magnetic influence, with strongest to weakest.

The planetary population observed by the Kepler satellite orbits primarily F, G and K stars. So, we choose the host star to be a solar mass star in its pre-main sequence stage. The fiducial parameters we considered are shown in Table 1. The stellar surface magnetic field is assumed to be dominated by the dipole component. We vary the stellar surface magnetic field while holding the mass, radius and surface temperature at the value listed.

We acknowledge that pre-main sequence stars are known to have drastically contracting radius over a disk life time of a few million years (e.g. Chabrier & Baraffe (1997)). However, we did not consider the effect of a contracting host star in this project as the contraction should not be substantial while the accretion flow from the protoplanetary disk is still present (Siess et al. (1997)). Any subsequent contraction after the protoplanetary disk has been dispersed would not affect the planet in the scope of this project.

Planetary migration is determined by the gravitational interaction between the planet and the natal gas disk, and so planetary migration is also affected by the time evolution of the background protoplanetary disk. To account for this effect, we compute the accretion rate with time evolution following Alexander & Armitage (2007). This model describes a protoplanetary disk, evolving under the combined influence of viscosity and a photoevaporative wind off the disk, driven by the high energy radiation from the central star. The wind eventually opens up a gap at larger radii, cutting off the gas supply to the inner disk, resulting in a rapid drop in the accretion rate as the decoupled inner disk grains onto the star. The calculated accretion rate used for this work is shown in Figure 2. We have shaded the slow dissipating phase and the rapid dispersing phase respectively. The accretion rate initially decreases slowly on a timescale of several Myr, driven by the standard viscous evolution of the disk. However, at some point ($\sim$2 Myr in the case shown in the figure), the the mass loss rate due to photoevaporative winds becomes comparable to the disk accretion rate, and the winds disconnect the inner disk from the outer disk, shutting off the replenishment of material from larger radii. After this, the accretion rate onto the star drops precipitously and will have an important influence on when and where the migration will freeze out.

To calculate the temperature and surface density profiles, we adopt the steady-state disk model (e.g. Pringle (1981)):

where $\Sigma$ is the surface density of the disk, $\nu$ is the kinematic viscosity and is assumed to be equal to $\alpha H^{2}\Omega$. No sink other than the central star is considered and the same accretion rate is applicable throughout the disk.

Under the above assumption, we consider the heating of the disk to be caused by both viscous dissipation and passive heating from the host star. The disk is cooled radiatively. The passive heating is computed following Ueda et al. (2017), which consists of 4 regions. The regions, described in order of their radial proximity to the host star, are a dust free region, a dust halo region, a dust condensation front, and an opaque thick disk region. Our temperature profile and surface density profile are then described by the solution to the following equations:

where $T$ is the temperature of the disk at the mid-plane. $R$ is the radial location in the disk. $\kappa$ is given in Bell et al. (1997), which formulates their opacity in the form: $\kappa=\kappa_{0}\rho^{a}T^{b}$. $\sigma$ is the Stefan–Boltzmann constant. The local passive heating term, $\Gamma(R)$, is defined according to appendix A.2.1.

The final element of our model are torques that drive planetary migration. As we are mainly interested in Super Earths and Sub-Neptunes, we chose to focus on Type I migration; the mass range for Super Earths and Sub-Neptunes should not substantially modify the surface density of the disk for the majority of the disk’s lifetime. We have checked the applicability of Type I torque for our planets and found that it is valid until less than 0.1 Myrs before the complete dispersion of the disk. At that point, the disk surface density is too low to drive any substantial migration.

We combine all the components discussed in Sections 2, 3.1, 3.2, and 3.3 to develop a time-dependent, 1D disk model that contains the transition zone. The model is then used to compute the torque as described in Section 3.4. The computed torque is coupled with Mercury N-body simulator (Chambers, 1999) via the user defined force routine, practically. We use the Bulirsch-Stoer algorithm with in Mercury to perform the integration.

To lower computation cost, we pre-generate the disk models before loading them into Mercury. The pre-generated disk models all have radial resolution of $5\times 10^{-4}$AU below 1AU and $5\times 10^{-3}$AU beyond. When values in between radial resolution are required for calculation, they are evaluated via spline fitting. The same is performed for the evaluation of the power-law indexes of both temperature and surface density profiles. We update the model every 500 years before the simulation time reaches 1.5Myrs and switches to every 50 years to keep up with the faster disk evolution. The migration torque is, on the other hand, updated at each time step that Mercury evaluates orbital elements.

For each set of simulations, we applied the combination of the stellar magnetic field strength and the power-law index listed in table 2 to the disk model and performed 100 unique simulations. Each individual simulation contains a single planet and an unique initial location. Planet masses were chosen at integer values between 1 and 10 Earth masses ($M_{\oplus}$), 10 per mass value, and initial semi-major axes were chosen between 0.3 and 1.2 AU in intervals of 0.1AU, 10 per location. Initial eccentricities (e) and inclinations (i) are both set to 0.1. All other orbital angles are chosen randomly.

By employing the torque prescription as laid out in section 3.4.1, we illustrate the nature of the torques for different disk locations and planet masses, using the disk profile of our fiducial model, in Figure 4. The colour map indicates the strength and sign of the calculated torque assuming circular planetary orbits. The black region indicates that the torque is essentially zero (because the disk surface density is zero inside the magnetospheric cavity The colors from green to blue indicate negative torques, which move the planet inwards, towards the star. Colours from yellow to red indicate positive torques, which will cause the planet to move outwards. We see that there is a sharp transition from inwardly directed to outwardly directed torques at orbital periods $\sim 30$ days. This is the location of the peak of the surface density in Figure 3, and the boundary between the outer transition zone and the generic viscous disk zone. This torque reversal implies the existence of a ‘planet trap’, where the migration will stall.

The construction of the torque map is generally governed by the balance between the Lindblad torque and corotation torques. The blue and green area is dominated by the Lindblad torque while the red and yellow area is dominated by the corotation torques. The cause of these two distinct regions is related to the value and sign of the corotation torques. Through out the entire radial extent of the disk, Lindblad torque remains negative. Even when the power-law index turns from positive to negative moving into the transition zone, both the value and sign of the Lindblad torque remain largely the same. On the other hand, the corotation torques jump from negative at beyond the transition zone to positive when in the transition zone. In the transition zone, summation of the corotation torque and the Lindblad torque turns out to be positive and leads to a torque reversal. It is important to note that the above description applies mainly to low eccentricity cases of $e<0.01$. As eccentricity increases, the corotation torques weaken and Lindblad torque would once again dominate in the transition region and would effectively erase the torque reversal. However, in this paper, we only considered single planet cases which, in combination with eccentricity damping from disk interactions, limited the eccentricity to well below 0.01 for the majority of their life time in the disk.

In the previous section, we consider characteristic disk features at a certain accretion rate. Here we explore how these features evolve with time, following disk evolution.

When disk accretion rates decreases with time, two major effects related to the disk profile are seen. First, features of both the surface density and the temperature profile move outward. Second, the surface density decreases as the accretion rate decreases. As disk features are crucial in both interpreting the result and comparing with similar works, we have included in Figure 5 the time evolution of both the inner edge (black) and the outer edge of the transition zone (orange) throughout the disk’s lifetime. As the accretion rate drops, features in our model tend to move outward. For the inner edge, the outward movement is due to the magnetic pressure balance requirement. For the outer edge of the transition zone, the outward movement is due to the effective viscosity ($\alpha$) being inversely proportional to the accretion rate. So as the accretion rate drops, both of those features moves outward but at different rates and this leads to the overall widening of the transition zone. On the other hand, the steady decrease in surface density is a consequence of the steady-state assumption (equation 5). We will discuss the effect of these outward-moving features in section 4.3

Our disk model represents a refinement of the widespread model for disk-mediated planetary migration by adopting a stopping criterion motivated by a detailed physical model for the manner in which a stellar magnetic field sculpts the inner edge of a protoplanetary disk. More specifically, we refined the stopping locations of planetary migration in magnetically truncated disks; when the transition zone is taken into account, migrating planets can be halted further away from the central star.

To further illustrate the consequences of this effect, we over plotted in Figure 5 the typical time-evolutions of the locations of migrating planets for a 1 and 10 $M_{\oplus}$ planet. We denote them as the generic tracks.

The initial inward migration features planets moving from their starting position towards the planet trap located at the outer edge of the transition zone. Once planets reach the trap, migration halts and they largely follow the evolution of the disk and remain coupled to the trap. As the accretion is constantly decreasing in our model, the trap moves outward and brings along any planets that are coupled to it. At a later time, usually 0.5 to 1 Myrs before the gas disk is completely dispersed, the disk evolution time scale becomes shorter than the outward migration timescales of the planets. This marks the beginning of the decoupling process, where the planet migrating outward cannot catch up with the outward movement of the outer edge of the transition zone. The planets now lie interior to the location of the torque reversal and within the transition zone. So they migrate outward, but at a rate slower than that at which the trapping location and, overall, the transition zone, are moving outwards. Planets are thus exposed to areas of the disk with ever lower surface density and weaker outward torques. So the planets reach their final location asymptotically.

The difference in decoupling time between the 1 and 10 $M_{\oplus}$ case is caused by the planet mass dependence of the torque; the 10 $M_{\oplus}$ planet experiences a stronger Type I torque. Due to the stronger torque, higher mass planets are coupled to the planet trap for a longer time and decouple from the disk at more distant locations. Further discussion of the effect of planet masses on the final locations of planets can be found in section 4.5. The generic track is seen in all $n=3$ simulations, where the planet traps remained similar to the description in section 4.1 through out the lifetime of the disk.

Following the overall picture of the whole migration process, we will discuss in more detail the individual phases and mechanisms of the decoupling process.

Although the generic case describes over 90% of our simulated planets, some planets do exhibit variations in their orbital evolution. One such variation, which we denote as the late-arrival case, occurs when a planet arrives at the planet trap location at a late time (for instance, if it began migrating from further out). As the disk is already largely depleted, the strength of the torque reversal at the planet trap location is too weak to hold the planet. So the planet migrates inwards asymptotically to it’s final location, as shown in the grey points in Figure 5. This late-arrival track is seen in about 5% of all the simulations. This type of migration track is typically seen in our simulations when the spawning location of the planet is larger than 0.8AU and with planet masses less than or equal to 2 $M_{\oplus}$.

To illustrate the importance of factors such as planet mass and starting location on the final location, we focus on the simulations set B1000n3, which shares the same magnetic field as our fiducial scenario: $B=1000$G and $n=3$. Figure 6 summarizes the results which show mass and initial location dependency for the lower mass end.

First, we can clearly see that the spread in final location is correlated with planet masses — with lower mass planets located closer to the host star and the heavier planets located further out. The only exception to the general trend comes from the stragglers resulted from the late-arrival case and are represented by the 1 and 2 $M_{\oplus}$ planets as seen on the right of the rest of the distribution.

Second, the initial locations of planets have little effect on the final location of the planets except for the least massive planets. As mentioned in section 3.5, a range of initial locations are considered for each set of simulations. However, despite the difference in initial location, all planets with mass $>3$ $M_{\oplus}$, and the majority of the 2 $M_{\oplus}$ planets, share virtually identical final locations with their peers with the same mass. The exceptions to this trend are, again, the results of the late-arrival track which is a consequence of the slower migration speeds for lower mass planets.

A high value of magnetic diffusivity in the disk would allow the magnetic field to penetrate further outwards, making the field profile more shallow and lowering the value of $n$. Alternatively, a low value of diffusivity would make it difficult for the field to penetrate the disk, in the face of the inward advection of material, and therefore this would make the profile steeper, increasing the value of $n$. Using our normalisation to a fixed global magnetic flux, lowering $n$ would amount to a weaker magnetic field at locations within the disk, while larger $n$ imply stronger local magnetic fields. A discussion on the exact effect on the magnetic field within the disk can be found in appendix A.1. To demonstrate the effect of changing power law index ($n$) on both the surface density and temperature profiles, we have included the $n=2.75$ case as gray dash lines on fig 3. The main difference between the fiducial case and the $n=2.75$ case is that the peak in surface density moved outward as the power-law index decreased and the drop off in surface density towards the host star is slower in the $n=2.75$ case. On the other hand, the temperature profile is largely identical. These behaviors are consistently seen in all non-dipole ($n<3$) simulations. The outward shift in peak location and slower drop off rate in surface density have subtle but important effects on the migration process which we will discuss in the next two sections.

There are two effects of decreasing n value, namely an increase in radial extension of the transition zone and an overall decrease in torque strength in the vicinity of the outer edge of the transition zone. For example, the transition zone outer edge extended from 0.13AU to 0.14AU when we decrease from $n=3$ to $n=2.75$ (Figure 3). On the other hand, the torque value experiences a drop of about $5\%$ near the immediate vicinity of the trapping location as a result of a lower surface density and a shallower surface density profile. The overall consequences of the two competing effects above is a general earlier decoupling time and, hence, smaller final semi-major axis for the planets in the systems. We will discuss the mechanism that leads to above outcome and the difference in decoupling time between the two example cases in the following section.

For the majority of parameters we considered, the behavior of planets are consistent with the fiducial case — Most planets follow the generic migration track with a small number of them following the late-arrival track. This behavior in migration track is true for all magnetic field values we considered. While the same types of tracks are observed, the decoupling time is observed to be consistently earlier for a given planet mass and given magnetic field strength when n is lowered. We attribute this earlier decoupling time to the lowering in torque value in the immediate vicinity of the trapping location as described in section 5.2. For example, the 1 $M_{\oplus}$ planets in the $n=3$ case decouple at around 1.07Myrs while the same type of planets in the $n=2.75$ case decouple at around 0.225Myrs. This earlier decoupling time is consistently seen in all simulations with non-dipole ($n<3$) magnetic profile.

Furthermore, we see a deviation from the normal-and-late-arrival track scenario when the power-law index of the magnetic profile (n) drops below 2.65. The full scope of the behavior seen in that situation can be found in appendix C

With decoupling time tightly tied to the final location of planets, it is natural to expect an overall smaller final semi-major axis as n decreases. In the next section, we will discuss the outcome of our simulations with all effect included.

Of the 11 simulation sets, 9 of them, which span $100G\leqslant B\leqslant 1000G$ and $2.5<n\leqslant 3$, behaved similar to our fiducial case (section 4). The median final locations from individual simulation sets range from $7.1\times 10^{-2}$AU (7 days) to $3.1\times 10^{-1}$AU (64 days), indicating that the model places planets at a wide range of locations, depending on the strength of the magnetic field. These simulations contain roughly the same mix between generic type and late-arrival type migration tracks. This similarity leads them to all have the same mass distribution as described in section 4.5. The individual simulation sets statistic, including the $1~{}\sigma$ equivalences and the median location, can be found in Table 2.

Two sets of simulations yield a behaviour different from the others. These have power-law index of $n\leqslant 2.5$ and yield median final locations of $<0.05AU$ (4 days). The final locations distributions of these two simulations has little to no dependency in mass and are tightly packed at virtually the same location, because the planet trap is either partially deactivated or too weak to halt the migration to the inner truncation of the disk, respectively. These cases are discussed in more detail in § C.1 and C.2.

Grouping together all simulations we have, we found that the overall result can be described more clearly by designating each of the simulations with their Normalized B-field ($B_{nor}$, see equation 9) instead of their stellar surface B-field ($B_{0}$) or the local magnetic field at the final location of planets. Disks in simulations with the same $B_{nor}$ are threaded with the same total magnetic flux, even if they have different power laws $n$. This correlation may appear counter intuitive, as the normalized B-field is not the one that set the location of the edge of the transition zone – where the planet trap should reside. However, if one looks at section 4.3, it is common for planets to have further interactions with the dispersing disk after decoupling begins. Ultimately, it is perhaps not surprising that the final location correlates more with a global measure of the magnetic influence on this disk than with a localised one. In Figure 7, we plotted the final median location against the normalized B-field($B_{nor}$) of each individual simulation set. The accretion rate used in the $B_{nor}$ calculation is assumed to be $10^{-8}M_{\odot}/\textit{year}$. We can see that, when we plot all of our simulation in the normalized B-field and semi-major axis space, the high magnetic field cases, between 1000G and 100G, form a tight relationship in log-log space. The solid black fitted function is calculated based on the the $n=3$ cases only. Hence, the fit is independent of the normalization method used. Yet, the fit also appears applicable to most of our non-dipole simulations. The fit only fails at the lowest field values, when the magnetic effects are not strong enough to halt the migration and the planets migrate to the inner edge.

The fit function is of the following form:

Where $R_{final}$ is the median final location for the given normalized B-field ($B_{nor}$).

Based on the above result, we are confident that the range for which our model can produce a transition zone based planet trap exterior to the inner truncation includes the parameter space $100\leq B_{*}\leq 1000$ with $2.75\leq n\leq 3$.

We compare the final locations of our simulated planets to observation by directly comparing the result from B100n3 with the Super Earths and Sub Neptune occurrence rate outlined in the California-Kepler Survey IV. (CKSIV ) (Petigura et al., 2018) in Figure 8. The B100n3 simulation is chosen for this comparison as it does not have the complication of normalization as required for the non-dipole profiles while its median final location falls close to the observed drop-off in super Earth occurrence rate. One caveat of this comparison is that the CKSIV occurrence rates are calibrated in terms of planet radius, while our simulations are in terms of planet mass, with the relationship to size muddied by unknown envelope fractions. As the mass of Super Earths can be up to $10M_{\oplus}$ and that of Sub-Neptunes can be between $3$ and $20M_{\oplus}$ (Zeng et al., 2019), the 1 to $10M_{\oplus}$ mass range in our simulations does not represent a single type of planet but a mix between the two. Hence, it is not possible to compare to either one of the population quantitatively. However, a qualitative comparison between the observed occurrence rate and our result can still provide insight into the strengths and shortcomings of our model. Namely, the distribution of the final locations is of most interest.

In Figure 8, the CKSIV occurrence rates for Super Earth and Sub Neptune are shown as the blue and green curves along with the errors shaded, respectively. The kernel density estimate of the result of B100n3 is shown as the black curve. The Gaussian used in the estimate is chosen to have standard deviation of $0.001AU$. The fit using convention laid out in Howard et al. (2012) (HOW12) is shown in orange. It is clear that the final location for our simulation is positioned between the drop-offs of Super Earth and Sub Neptune and, hence, can qualitatively match with the drop-off location as seen in the CKSIV fits. We did not expect an exact match between our result and a specific planet type as we cannot distinguish between Super Earth and Sub Neptune in our simulations. However, given that our model can place planets in a distribution that mimics the observed population with only single planets in each systems, we consider that our model as promising.

Our single planet simulations match the characteristic feature of the turnover in the observed period distribution, but does not match the drop-off below 10 days and the plateau beyond 10 days. We will address these features in a future publication containing multiple planet simulations. We postulate that, once more planets are added to our simulation, the plateau area should be filled as interior planets are halted by the planet trap and interact gravitationally with exterior planets. The discrepancy of the drop-off at short periods could also potentially be explained by multiple planets as the gravitational influence of later arriving planets may force the innermost closer to the star. Furthermore, the tidal damping of any eccentricity resulting from dynamical interactions will also drive the innermost planet closer to the star (Hansen & Murray, 2013). Furthermore, using equation 8, we can then predict the parameter space that can place a planets around the location of the population drop-off. By restricting the median final location to be between 8 and 12 days orbit (i.e. $7.8\times 10^{-2}$ and $1.03\times 10^{-1}$AU), we get the applicable normalized B-field values between 67 and 113G. Mapping these normalized B-field range back into the stellar B-field and stellar power law index space, we get the result in Figure 9. The colored region represents the area of the parameter space that planets can be placed in the vicinity where a slow drop in occurrence is observed. The color code represents the orbital period at which the median of the planets would be located at. The minimum and maximum stellar magnetic field in this inferred parameter space is 67 and 180G respectively. As we don’t have any dipole based simulation with magnetic field below 100G, we conclude that we are only confident that a stellar magnetic field of 100-180G is able to reproduce observation.

Both the full functional range where migrations can be halted by the transition zone (100-1000G) and the observation reproducing range (100-180G) of the magnetic field for our model are supported by observations on T Tauri stars. As outlined in Johnstone et al. (2014), their observed stellar dipole field strength for T Tauri stars range from 80 to 1720G with a substantial fraction of them being below 400G. This range matches with the functional range of 100 to 1000G where our model can produce a planet trap associated with the outer edge of the transition zone. More importantly, a large fraction of the T Tauri star exhibiting dipole field below 400G indicates that our applicable range for matching with the observed super Earth and sub Neptune population is fairly common.

Our disk model shares sectional similarity with other works in terms of the accretion model considered. On the other hand, we employed a drastically different inner disk from other disk models and that resulted in different migration behavior from other disk models.

Our accretion model, which is based on Alexander & Armitage (2007), can be seen as a bridge between works that consider the long term evolution (e.g. Hueso & Guillot (2005), as HG05 below) and those that consider the final dispersion phase of the disk (e.g. Liu et al. (2017), as LOL17 below).

For most of the viscously dominated disk phase, our model behaves comparably with HG05, wherein the accretion rate is described as a simple power law between the accretion rate and the disk life time, motivated by observations. The full model in example 1 of HG05 span $10^{-7}$ to $10^{-9}M_{\odot}/\textit{yrs}$ and has form similar to our slow dissipating phase. This similarity is expected as their model aims at modeling the long term disk behavior which should be dominated by viscous evolution. However, the migration depends sensitively on how the accretion rate drops towards the end of the disk lifetime, and the drop in $\dot{M}$ in our model is based on a direct calculation of the truncation of accretion due to photoionization of the disk material (see section 3.2). This drastic decrease in accretion rate is more abrupt that that used in LOL17, who used an exponential decay model that transitions from a timescale of a few Myr to $10^{5}$yrs during dispersal. This is potentially an important distinction, especially for planets that remain coupled with the planet trap to the rapid dispersal phase (i.e. planets with masses higher than 10 $M_{\oplus}$), as it means that the complete freeze-out of migration is even faster in our model. However, given the upper limit of our mass range, this distinction is not observed in our current result. On the contrary, we found that the rapid dispersal phase might have minimal effect on planets of mass up to 10 $M_{\oplus}$. We will discuss the reasoning below based on the decoupling time seen in our simulations.

As most of our planets follow the generic track, planets typically decouple from the disk around 1 to 0.5 Myrs before the complete dispersion of the gas disk. During this time range, the accretion rate in our model corresponding to $\approx 10^{-9}M_{\odot}/\textit{yrs}$. Similar accretion rate occurs in the HG05 model at similar time — around 0.8 Myrs before the gas disk is assumed to be completely dispersed. So planets are expected to decouple around that time under the accretion model of HG05. From the above analysis, we speculate that the final stage of the dispersion process in common models, in combination with our torque treatment, has little to no effect on single planets similar to those in our simulations.

Another important distinction between our model and previous works is that our model for the magnetic field interaction leads to a shallower slope in the surface density in the inner parts of the disk. Comparing to HG05, which is derived with similar assumption on the heating sources, our temperature and surface density profiles agree largely at the 1AU and beyond scale. Interior to this, the magnetically enhanced viscosity in our model leads to a surface density is considerably lower with a positive power-law index while the surface density seen in HG05 remains higher than the rest of the disk with a close to 0 or negative power-law index. This difference in inner disk structures implies that migrating planets are expected to be trapped further out than those in a HG05 disk. Similar argument also applies to other works such as Emsenhuber et al. (2021), Izidoro et al. (2017) and LOL17 where simple cut-off functions with short scale drop-offs are employed. Hence, our disk model, under the same stellar magnetic field prescription, predicts a further out trapping location compared to other works.

While our simulations indicate both a strong mass–final location trend and normalized-magnetic-field trend, our result requires a few word of caution because of the simplicity of our model on a few aspects, namely, the mass evolution, the lack of multiplicity and viscous $\alpha$ prescription we chose.

First, the mass in our simulation would represent the primordial mass as the planet was still in the disk or right after the disk has been dispersed. Further mass modification process such as wind driven atmospheric lose or photoevaporation will likely further modify the mass and radius of the planets and lead to possibly a different mass trend that this work does not take into account. Further work focusing on atmospheric evolution both during and after the disk phase will likely lead to a more realistic mass and size distribution.

Second, we only considered single planets in this work while it is well-known that super Earths frequently come in multiples. So while our current result points to a rather weak stellar magnetic field, we speculate that the suitable magnetic field range will increase with the introduction of a second or third planet. This speculation is based on the understanding that planet-planet interaction will likely increase the eccentricity of both planets leading to a weakening of the corotation torques. So planets would be able to venture deeper into the transition zone. So, to maintain the median final location at 10-days orbit while still keeping the drop off below 10 days, the transition should occur further out to compensate for this effect. So the suitable stellar magnetic field value should increase with higher multiplicity. Future works are planned to explore the effect of multiplicity on the required stellar magnetic field, radial distribution, mass distribution and orbital architecture.

Third, our viscous $\alpha$ description does not accurately represent the effect from the stellar magnetic field at low magnetic field strength. As pointed out in section 2.1.2, the expression we used is based on a fit (Salvesen et al., 2016) up to $\beta\lessapprox 10^{5}$. However, in order to reach the required $\alpha=0.01$ for the edge of the transition zone, it requires $\beta\approx 5.5\times 10^{5}$. While Salvesen et al. (2016) does not explicitly state that the fit will break down beyond the suggested range, Suzuki et al. (2010)suggest a much shallower relationship should be in place for $10^{6}\lessapprox\beta\lessapprox 10^{5}$. However, since the effect from the viscous $\alpha$ profile is localized, we expect the modification to our median final locations to be modest. As the overall effect would require a detailed analysis into the exact balance between the transition zone outward extension and the corotation torque weakening, we leave this topic for future works.