# Combined Effects of Disk Winds & Turbulence-Driven Accretion on Planet
Populations

Matthew Alessi1 and Ralph E. Pudritz1,211footnotemark: 1
1Department of Physics and Astronomy, McMaster University, Hamilton, ON L8S
4M1, Canada
2Origins Institute, McMaster University, Hamilton, ON L8S 4M1, Canada
E-mail<EMAIL_ADDRESS>(MA<EMAIL_ADDRESS>(REP)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Recent surveys show that protoplanetary disks have lower levels of turbulence
than expected based on their observed accretion rates. A viable solution to
this is that magnetized disk winds dominate angular momentum transport. This
has several important implications for planet formation processes. We compute
the physical and chemical evolution of disks and the formation and migration
of planets under the combined effects of angular momentum transport by
turbulent viscosity and disk winds. We take into account the critical role of
planet traps to limit Type I migration in all of these models and compute
thousands of planet evolution tracks for single planets drawn from a
distribution of initial disk properties and turbulence strengths. We do not
consider multi-planet models nor include N-body planet-planet interactions.
Within this physical framework we find that populations with a constant value
disk turbulence and winds strength produce mass-semimajor axis distributions
in the M-a diagram with insufficient scatter to compare reasonably with
observations. However, populations produced as a consequence of sampling disks
with a distribution of the relative strengths of disk turbulence and winds fit
much better. Such models give rise to a substantial super Earth population at
orbital radii 0.03-2 AU, as well as a clear separation between the produced
hot Jupiter and warm Jupiter populations. Additionally, this model results in
a good comparison with the exoplanetary mass-radius distribution in the M-R
diagram after post-disk atmospheric photoevaporation is accounted for.

###### keywords:

accretion, accretion discs – planets and satellites: composition – planets and
satellites: formation – protoplanetary discs – planet-disc interactions

††pubyear: 2020††pagerange: Combined Effects of Disk Winds & Turbulence-Driven
Accretion on Planet Populations–A

##  1 Introduction

How are the observed properties of exoplanet populations connected to physical
processes and planet formation in disks? Signatures of planet formation on
extant populations could be encoded in several ways such as the distributions
of their orbital characteristics, planetary radii, or in the compositions of
their atmospheres. As an example, a classic paper argues that the C/O ratios
of materials in planetary atmospheres could constrain where a planet may form
in a disk in relation to water, $CO$, or $CO_{2}$ ice lines (Öberg et al.,
2011). A difficulty in making a quantitative connection is that planets likely
migrate while they form so that atmospheric composition also reflects the
collection of materials that they acquire from different parts of the disk
(Cridland et al., 2016; Madhusudhan, 2018; Cridland et al., 2019a, b).

Surveys have clearly shown that exoplanetary populations span a large range of
properties as is apparent in the mass-semimajor axis (hereafter M-a) and mass-
radius (hereafter M-R) diagrams (Batalha et al. 2013; Rowe et al. 2014; Morton
et al. 2016; see also review by Winn & Fabrycky 2015). The statistical
properties of populations have been significantly improved over the last 5
years: the most complete catalogue of transients in the Kepler data contains
4034 planet candidates (Thompson et al., 2018); RV methods are in development
to reach 10 cm per second in order to detect Earth mass planets (Fischer et
al., 2016); and new surveys are measuring planet masses around other kinds of
stars such as M dwarfs (Reiners et al., 2018).

Planets form in protoplanetary disks and the surge of spatially-resolved disk
(sub)millimeter observations has revealed a striking amount of substructure in
disks’ solid distributions (i.e. ALMA Partnership et al. 2015; van Boekel et
al. 2017; Andrews et al. 2018). These structures may result as a consequence
of gap opening in the dust due to the presence of forming planets (e.g. Dong
et al. 2015; Fung & Chiang 2016). Young planets forming in such disks have now
been directly observed in the PDS70 system (Keppler et al., 2018; Müller et
al., 2018; Haffert et al., 2019). Recent disk surveys show that these disks
have an inherent scatter in their global quantities, such as masses and radii
(Ansdell et al., 2016; Pascucci et al., 2016; Tazzari et al., 2017; Ansdell et
al., 2018; Andrews, 2020; van der Marel & Mulders, 2021; Long et al., 2022).

The scatter in the exoplanet M-a and M-R diagrams embodies a wide range of
processes. These include the range of initial conditions for newly formed
disks (disk masses, radii, lifetimes, metallicities) that play a huge role in
planet formation; the dominant processes that dictate disk evolution and hence
planet formation within them (Armitage, 2018); the dynamical evolution of
planetary systems due to N-body interactions between planets during both the
formation phase in the disk (Emsenhuber et al., 2021) as well as after disk
dispersal (Chatterjee et al., 2008; Ford & Rasio, 2008) ; and the billions of
years of planetary evolution that shaped planetary atmospheres.

The central physical process that dictates the evolution of accretion disks
and planet formation is how disk angular momentum is transported. Planet
formation theory has for decades assumed that turbulent stress is the sole
agent of disk angular momentum transport in the core accretion picture. The
strength of the disks’ turbulent intensity is parameterized by the
$\alpha_{\rm{turb}}$ parameter. In the standard theory of viscous accretion
(Shakura & Sunyaev, 1973): the viscosity arising from such turbulence is
$\nu=\alpha_{\rm{turb}}c_{s}H$, where $c_{s}$ is the sound speed and $H$ is
the pressure scale height. The ratio of the rms turbulent to thermal pressure
in the disk, measured by the value of $\alpha_{\rm{turb}}$, has a significant
effect on disks’ evolution timescales in viscous disk theory. Its value
determines the disk accretion rate as well as the outward spread of material
in the outer regions where angular momentum is ultimately deposited. Its
magnitude also plays a central role in planet formation controlling planet
migration through co-rotation torques as well as gap opening processes ( for a
recent review, see Nelson (2018)). The physical mechanism that generates
turbulence was thought to be the magneto-rotational instability (MRI) Balbus &
Hawley (1991) although weaker turbulence can arise from a number of
hydrodynamic and thermal instabilities can arise in vertically stratified
disks (Klahr et al., 2018).

In order to test whether turbulent viscosity by itself is sufficient to
explain the data, it is necessary to measure $\alpha$. This, until recently,
has proven to be frustratingly difficult to do. Previous indirect attempts
include matching the related evolution timescale to observationally inferred
disk lifetimes (Hernández et al., 2007), or matching the related accretion
rates to those observationally inferred through H$\alpha$ emission (Hartmann
et al., 1998). These approaches result in estimates for
$\alpha_{\rm{turb}}\sim 10^{-4}-10^{-2}$. It is important to note that these
estimations assume that disk evolution is _solely_ driven by turbulence. The
first direct measurements of disk turbulence have now been acquired by ALMA
observations. The results show that disks exhibit a lower levels than expected
from the earlier inferences. Observed levels of disk line-broadening (Flaherty
et al., 2018) and studies of dust properties within pressure traps (Dullemond
et al., 2018; Rosotti et al., 2020) have constrained the strength of
turbulence in protoplanetary disks to an upper limit of
$\alpha_{\rm{turb}}\lesssim 0.007$. The observations of disk accretion rates
(Manara et al., 2019) have orders of magnitude greater variation than can be
accounted for by a fixed fiducial value of $\alpha=2\times 10^{-3}$ (Mordasini
et al., 2012) used in population synthesis modeling. Mulders et al. (2017)
argued that a dispersion in the value of $\alpha$ by 2 dex is needed to
explain the variations in the observations in comparison to the use of a
single value of this parameter. Limited amounts of outward disk spreading
observed for disks in the Lupus star forming region suggest low disk viscosity
(Trapman et al., 2019). Detailed study of the HL Tau disk has shown that
whereas the accretion rate onto its central star by viscous torques would
require $\alpha_{\rm{turb}}\simeq 10^{-2}$, the actual dust scale height is a
factor of 5 smaller than the gas scale height which puts an observational
upper limit of $\alpha_{\rm{turb}}\simeq 10^{-4}$ on the disk turbulence
(Pinte et al., 2016). Theoretical studies show that when the disk turbulence
is high, dust trapping in pressure maxima is far less efficient than in the
case of low turbulent amplitudes (Pinilla et al., 2020). The most recent
compilation of a wide range of data and methods finds that typically
$\alpha_{\rm{turb}}\sim 10^{-3}-10^{-4}$ (Pinte et al., 2022). This extensive
set of observations suggests that a range of turbulent viscosities is
intrinsic to protoplanetary disk populations and we incorporate this into the
evolution of our planetary populations.

Disks with low values of disk viscosity $\alpha_{\rm{turb}}\sim 10^{-4}$ have
particularly interesting implications for planetary migration. Low viscosity
allows outward directed co-rotation torques to push low mass forming planets
out to 10s of AU via Type I migration before they eventually turn inward as
they get caught up in the heat transition dynamical trap (Speedie et al.,
2022). Low viscosities are expected in more massive disks which are better
shielded from external sources of radiation which ionize disks and promote MRI
turbulence. These theoretical results may provide an explanation as to why
only a minority of disks are large ring/gap systems (van der Marel & Mulders,
2021).

Hydromagnetic disk winds have long been known to be effective in transporting
disk angular momentum. Crucially, they provide a robust explanation for the
ubiquitous appearance of jets and outflows during star formation. They are
connected with protostellar disks associated with the full range of stellar
masses (see recent reviews (Pudritz & Ray, 2019; Pascucci et al., 2022). Early
theoretical work predicted that centrifugally driven MHD disk-winds could
efficiently transport disk angular momentum through an outflow of disk
material along magnetic field lines threading the disk (Blandford & Payne,
1982; Pudritz & Norman, 1986; Pelletier & Pudritz, 1992; Ferreira & Pelletier,
1995; Ouyed & Pudritz, 1997; Spruit, 2010; Bai, 2016; Tabone et al., 2020).
These can also be very efficient in that a small wind mass loss rate can drive
much higher accretion rates. Recent ALMA observations have resolved and
measured the rotation of outflows and that they originate from extended disk
scales (Chen et al., 2016; Bjerkeli et al., 2016; Louvet et al., 2018; Zhang
et al., 2018; de Valon et al., 2020). These measurements show that up to
60$\%$ or more of the disk’s angular momentum is being carried out in some of
the observed rotating outflows.

Numerical MHD simulations are now sophisticated enough to follow the coupling
and evolution of magnetic fields in dense poorly ionized protoplanetary disks.
This involves careful simulations of non-ideal MHD processes; Ohmic
dissipation, ambipolar diffusion, and the Hall effect. Such works have
explored how disk evolution is tied to the detailed ionization structure that
depends on non-equilibrium chemistry (e.g. Bai & Stone (2013); Lesur et al.
(2014); Gressel et al. (2015); Bai et al. (2016); Gressel et al. (2020);
Rodenkirch et al. (2020)) and that disk winds carry off the bulk of the
angular momentum. As an example, MRI is suppressed by a combinateion of
ambipolar diffusion and Hall effect down to $\alpha_{\rm{turb}}\sim 10^{-4}$
at the disk midplane in the outer regions (5-60 AU) of the disk (Bai, 2015).
Several papers have investigated disk wind effects on disk structure including
Suzuki et al. (2016); Hasegawa (2016); Bai (2016). The lack of magnetically
driven turbulence in wide reaches of the disk implies that a variety of
hydrodynamic instabilities such as the vertical shear, convective
overstability, and Zombie Vortex instabilities can provide weak to modest
levels of turbulence in disks (see reviews Klahr et al. (2018); Lesur et al.
(2022)). These could be important to help drive planetesimal formation.

Turbulence and disk winds act quite differently in controlling disk evolution.
In the former angular momentum is transported radially outward resulting in
the gradual expansion of the disk while the latter moves it vertically and
away from the disk leading to disk contraction, known as advective disks
(Nelson, 2018). One of the most significant consequences of the advective
nature of wind driven disks is that they can in principle have higher surface
densities in their inner regions than turbulence models. This may be
significant in promoting planet formation there (Suzuki et al., 2016;
Chambers, 2019). A caveat to this is that the central surface densities will
also depend on the mass loss rate carried by the disk wind, a point discussed
in the next section. Migration in purely wind-driven disks has been analyzed
by McNally et al. (2017, 2018); Kimmig et al. (2020). Winds that are strong
enough can drive higher disk inflow speeds which in turn power co-rotation
torques that drive outward migration (Kimmig et al., 2020; Speedie et al.,
2022). For low mass planetary migration, the absence of turbulent diffusion in
wind driven disks means that local vortensity gradients, important for the
corotation torques on the planets, are not dissipated so that the direction
and magnitude of the wind driven torque depends upon the migration history of
the planet (McNally et al., 2017, 2018; Nelson, 2018; Paardekooper et al.,
2022).

Our earlier papers focused on MRI driven turbulence, and incorporated the
results into a planet population synthesis method. We assumed that angular
momentum transport from the inner dead zone of such disks would be carried off
by a disk wind (Alessi et al., 2017) which was not treated in any great
physical detail. We found correspondence with many features of the planetary
M-a and M-R relations by considering observationally-constrained ranges of
disk properties as inputs to the core accretion planet formation model. In
Alessi & Pudritz (2018), we inferred the value of forming planets’ envelope
opacities of $\sim$0.001 cm2 g-1 by comparing our models’ gas giant orbital
radius distributions with that of observations. Next, in Alessi, Pudritz &
Cridland (2020a), we incorporated a full treatment of dust evolution and
radial drift into our approach, and discovered having large super Earth
populations requires initial disk radii of the order $\sim$50 AU. This is in
accord with recent ALMA observations of disk radii, the majority of which are
rather compact in structure (van der Marel & Mulders, 2021). Finally, in
Alessi, Inglis & Pudritz (2020b), we incorporated a full solid disk chemistry
treatment to track planets’ compositions and also considered post-disk phase
atmospheric mass loss via photoevaporation to determine the fraction of
accreted atmospheres that escape. These both shaped our populations’ M-R
distributions, with photoevaporaton reducing the masses of our close-in Hot
Jupiters to give better agreement with the M-a diagram.

In this paper, we extend our previous investigations to address planet
formation and migration in chemically evolving disks whose dynamical evolution
is driven by _both_ turbulence and magnetohydrodynamic (hereafter MHD) disk
winds in much greater physical detail. The remainder of this paper is
structured as follows: In section 2; we will summarize the Chambers (2019)
self-similar analytical disk model that handles both viscous and disk wind
driven evolution. This is a new addition to our treatment and we detail the
additions and constraints we have made to it, summarize our planet formation
model, and define individual models’ parameter settings that will be used in
this paper. In section 3, we first examine the effects of different relative
disk winds and turbulence strengths on individual planet formation tracks
using fiducial disk parameters. This analysis is extended to disks with low
turbulence levels in Appendix A, where the disk outflow strength is examined.
We then extend our investigation to full planet population synthesis models in
section 4, where synthetic M-a and M-R distributions are shown that
incorporate a distribution of turbulent $\alpha$ settings as well as post-disk
atmospheric mass loss through photoevaporation. Lastly, in section 5, we
summarize our main conclusions and discuss extensions of this treatment that
will be considered in future work.

##  2 Protoplanetary Disk Model: Combined Evolution Through Turbulence and
Disk Winds

It is currently not feasible to compute the complete evolution of planet
formation and migration for the lifetime (3 - 10 Myr) of protoplanetary disks
by using the full 3D MHD treatments discussed in the Introduction. One useful
approach that captures the essential physics is to build on the simplicity of
the $\alpha$ formalism and to adopt two types of $\alpha$ for the disk - one
for the turbulence and another for a disk wind. These two-$\alpha$ models have
proven to be successful in creating better matches to observed populations in
the M-a diagram (Alessi et al., 2017; Ida et al., 2018; Bitsch et al., 2019;
Matsumura et al., 2021). The addition of disk winds appears to produce a pile
up of planets at around 1AU without invoking photoevaporation of the disk
(Matsumura et al., 2021) .

Here we focus entirely on the disk phase of planetary evolution. Planet
population synthesis models utilize the best and most basic physical models of
the various aspects of planet formation; the evolution of disk structure by
both turbulence and disk winds; disk astrochemistry; and the details of Type I
and II planetary migration. By computing thousands of planetary evolution
tracks which sample the initial distributions of mass, radius, and metallicity
in evolving disks, we build synthetic populations using a statistical approach
(Ida & Lin, 2004a, 2008; Alibert et al., 2011; Benz et al., 2014; Mordasini et
al., 2015; Alessi et al., 2017, 2020a, 2020b; Pudritz et al., 2018) which are
compared with the observations. As far as possible, we use the distributions
of disk initial conditions (mass, radius, metallicity) inferred via
observations.

Computing entire populations of thousands of planets forming in such evolving
disks is a daunting task when population synthesis is coupled with
astrochemical evolution of the disks. This is where the power of accurate
semi-analytical treatments for disk evolution comes into its own as they are
an effective and computationally inexpensive means to deduce the basic effects
of winds on disk evolution and planetary populations.

###  2.1 General Model Description

We adopt the formalism of Chambers (2019) to calculate the evolution of an
evolving disk under the combined action of turbulence and disk winds. This
analytic model is an extension of disk evolution via disk turbulence
(Chambers, 2009) and provides an elegant analytical treatment of the numerical
approach to disk winds explored by Suzuki et al. (2016). Our own approach
builds upon our previous successful integration of Chambers (2009) into our
existing planet population synthesis framework (Alessi & Pudritz, 2018; Alessi
et al., 2020a, b). Another advantage of our approach is that it allows for the
strengths of MRI and MHD-winds to be individually specified, so that their
relative effects can be ascertained.

There are three basic components that are linked together in our model. The
first is to include the physics of turbulent and disk wind torques on planets
to compute their migration and accretion of materials as a function of time
and disk radius. This analysis makes critical use of our disk astrochemistry
treatment. The second component is to perform population synthesis
computations involving thousands of models to confront our results of planets
formation under these combined torques with observed exoplanet populations in
the M-a diagram. Our final step is to to compute the related M-R diagrams by
solving planet structure and atmosphere equations using the cumulative
accreted materials and gases as well as photoevaporation effects on
atmospheres of close in planets . This is an ambitious program and its success
depends critically upon striking an optimal balance between the use of
evolving disk models that provide an accurate treatment of the physics of
planet formation and migration and the ability to compute thousands of model
planet formation histories over 10 Myr of disk evolution and planet formation.

Mass loss rates in disk winds are an important input into these models. It has
been pointed out that if disk winds are strong enough, they can carve out a
significant portion of the interior of the disk so that planet migration can
be slowed and a population of close in SuperEarths established (Ogihara et
al., 2018). This requires that wind mass loss rates are comparable to or
exceed the disk accretion rate; $\dot{M}_{\rm{wind}}/\dot{M}_{\rm{acc}}\simeq
1.0$. This limit of heavy disk wind mass loss conflicts with observations of
87 sources ranging from Class 0 to Class II protostars for which a large
number of which have a ratio of wind mass loss to accretion rates
$\dot{M}_{\rm{wind}}/\dot{M}_{\rm{acc}}\simeq 0.1$ (Watson et al., 2016;
Pascucci et al., 2022). Observations jet outflows with detectable [OI] line
emission show that $30\%$ of T Tauri stars with ages between 1-3 Myr have high
velocity component ( $v>40kms^{-1}$) emission wherein the average value of
$\dot{M}_{\rm{jet}}/\dot{M}_{\rm{acc}}$ is 0.07, with a spread of 0.01 - 0.5
(Nisini et al., 2018). While there is a class of disk wind solutions known as
tower flows that feature slow accretion flows and heavy wind mass loss rates
(Lesur et al., 2022) it is not clear that these are representative of typical
situations in protoplanetary disks. Thus, in this work we focus on examining
the limit of light fast winds that extract angular momentum very efficiently
with low values of jet mass loss to disk accretion rates (Pelletier & Pudritz,
1992; Zhu & Stone, 2018).

We now summarize in a non-technical way, the basic components of our models
that we used in the previous papers. The reader may refer to (Alessi et al.,
2020a) for all of the mathematical details. For disks undergoing purely
visccous stress we have always employed the Chambers (2009) self-similar disk
model for turbulent disk evolution. Most importantly, it enables a semi-
analytic treatment of disk and planet evolution that allows a numerical
approach capable of tracking planetary populations over disk life times (up to
10 Myrs). At any moment in time, the solution prescribes the spatial variation
of the disk surface density and temperature ranging from the inner, viscously
heated region that gradually transitions to the outer radiatively heated disk
region. The time evolution of the disk is dictated by the self similar
evolution and accretion due to viscous stress. We showed that these models
also correspond very well with observationally tested disk models by D’Alessio
et al. (1999). We also incorporated equilibrium chemistry in the evolving
disks using an effective Gibbs free energy solver that allows the composition
of materials accreted onto the forming planets Alessi & Pudritz (2018); Alessi
et al. (2020b) as well as the makeup of their atmospheres. We used this
technology in Alessi et al. (2020b) to compute the composition of planetary
populations as well as their atmospheres and employ it in this paper.

Throughout this and our earlier papers we have adopted a core accretion
picture of planet formation. Recent advances show that planetesimals are most
likely constructed by the rapid accretion of pebbles as a consequence of
streaming instability (Youdin & Goodman, 2005; Johansen et al., 2007; Bitsch
et al., 2015). In earlier papers we have employed dust models that have both a
constant dust to gas ratio (Alessi & Pudritz, 2018) as well as including the
radial drift of dust (Alessi et al., 2020a). Here we use a constant dust to
gas density which effectively leads to less compact planetary populations in
the M-a diagram. Given that planetesimal formation by the streaming
instability is so rapid (within $10^{5}$ yr) in comparison with the migration
and planetary accretion time scales, we assume the local surface density of
planetesimals follows that of the dust. This allows us to better focus on the
interplay of turbulent and disk wind torques.

Planetary accretion from the surrounding disk is treated by a series of
approximations adopted from standard planetesimal accretion scenarios. The
details of this are given in Alessi & Pudritz (2018) which we briefly
summarize. We start each computation of a planet formation track in a given
disk model from a so-called oligarch or planetary embryo whose mass is a
hundredth of an Earth mass, $10^{-2}M_{E}$. The first stage is the growth of
oligarchs that takes place via planetesimal accretion, which we compute
following the standard model of Kokubo & Ida (2002). The transition from the
oligarchic growth phase to gas accretion phases occurs when the planetesimal
accretion decreases to the point that the envelope pressure is insufficient to
support of the surrounding gas, which then accretes on to the planet,
following Ikoma et al. (2000); Ida & Lin (2008); Hasegawa & Pudritz (2012).

Accretion onto forming planets with masses exceeding this critical mass takes
place on the Kelvin-Helmholtz time-scale. The ability of a growing core to
accrete gas strongly depends on how well that gas can cool, so the opacity of
the atmosphere plays a major role. As already indicated in the Introduction,
we choose an envelop opacity for accreting planets that gives the best results
for matching the position of gas giant planets in the M-a diagram,
$\kappa_{o,env}\simeq 0.001$ cm2 g-1. For this, we use the fits. provided by
Mordasini et al. (2014) to numerical models of gas accretion to Kelvin-
Helmholtz scaling parameterization.

We always apply one set of K-H parameters per disk model regardless of where
the planets were accreting from. Coleman et al. (2017) computed the accretion
of gas onto planets of various masses, opacities, and locations in the disk
and found runaway gas accretion could be limited by the inability of gas to
cool sufficiently quickly in the hot inner disk regions. This work used
interstellar opacities from Bell & Lin (1994) which were artificially reduced
by up to 100 but still higher than our best fit. We surmise as we had in our
earlier work that an excellent physical analysis of atmospheric opacity, and
in particular the behaviour and contributions of dust grains within them, is
critical. Given these uncertainties, we expect that the masses of Hot Jupiters
formed in our dead zone may be upper limits with the consequence that one
might have more SuperEarths and fewer Hot Jupiters being produced than our
models would suggest (see Results).

It is well established that planetary migration is an inescapable consequence
of planet- disk interaction. At sufficiently low planetary mass, forming
planets do not open gaps and undergo well known Type I migration. Many works
have shown that Type I migration rates depend on the interplay of co-rotation
and Lindblad torques very near the planet and hence depend critically upon
local conditions there (Lyra et al., 2010; Hellary & Nelson, 2012; Dittkrist
et al., 2014; Baillié et al., 2016; Coleman & Nelson, 2016). Growing planets
rapidly migrate through disks until caught in regions of zero net torque where
outward co-rotation torques balance inward directed Linblad torques. In a
recently published paper that is in a sense connected with the present work,
we provide detailed simulations of planet evolution in evolving disks in
Speedie et al. (2022) that are based on detailed disk astrochemistry
calculations and using the detailed torques computed in (Paardekooper et al.,
2011). We address this further in the migration discussion below.

###  2.2 Inclusion of a disk wind into disk structure and planetary migration
theory

Disk winds have a number of effects on disk evolution. The magnetic torque
exerted by the rotating magnetized wind removes disk angular momentum driving
accretion through the disk (Blandford & Payne, 1982; Pudritz & Norman, 1986).
Disk winds also carry away mass from the disk although this is generally a
small fraction of the accreted mass. The extension of the self-similar
treatment of viscous disk evolution to include the possibility of a disk wind
was carried out by Chambers (2019), whose basic approach is an application of
detailed analysis presented in Suzuki et al. (2016) and Bai (2016).

Wind driven advection and disk mass loss contribute two new terms to the
standard equation that describes the evolution of the disk surface density
$\Sigma(r,t)$ :

$\frac{\partial\Sigma}{\partial t}=\frac{3}{r}\frac{\partial}{\partial
r}\left[r^{1/2}\frac{\partial}{\partial
r}\left(r^{1/2}\nu\Sigma\right)\right]+\frac{1}{r}\frac{\partial}{\partial
r}\left(rv_{w}\Sigma\right)-\dot{\Sigma}_{w}\;,$ (1)

Here, the first term is familiar in the standard treatment of disk surface
density evolution driven via MRI-turbulent viscosity $\nu$. Its self-similar
solution is used in the earlier purely viscous models (i.e. the Chambers 2009
disk model that we have used previously in Alessi & Pudritz 2018; Alessi et
al. 2020a, and Alessi et al. 2020b). The second and third terms correspond to
the additional effects of the disk wind that contribute both a stress driving
disk accretion (second term), and disk mass loss by outflow (third term). The
former is characterized by $v_{w}$ \- the inward radial velocity caused by the
disk wind, and the latter is characterized by $\dot{\Sigma}_{w}$ \- the rate
of surface density loss due to the wind outflow.

Chambers (2019) introduces three parameters that correspond with each of these
three physical effects:

  * •

$v_{0}$; the inward velocity of material at $r_{0}$ = 1 AU, and $T_{0}$=150 K
which is the disk temperature at $r_{0}$ caused _solely_ from radiation (see
equation 10 and following description). This parameter in large part sets the
initial disk accretion rate $\dot{M}_{0}$ which would be measured as the
accretion rate onto the star.

  * •

$f_{w}$; the fraction of $v_{0}$ that is caused by disk winds. Setting
$f_{w}$=1 corresponds to a pure disk winds scenario, while setting $f_{w}$=0
corresponds to a pure viscous evolution.

  * •

$K$; which characterizes the strength of the winds-driven outflow.

Solving the protoplanetary disk structure in this framework now introduces two
additional parameters beyond the one $\alpha$ viscosity parameter (Shakura &
Sunyaev, 1973) that is needed in the pure viscous scenario (i.e. Chambers
2009),

$\nu=\alpha c_{s}H\;,$ (2)

where $c_{s}$ is the sound speed and $H$ is the disk scale height. In pure
viscous models, $\alpha$ in equation 2 sets the strength of turbulent
viscosity. However, despite its definitions in the pure viscous scenario, we
have regarded this $\alpha$ parameter more generally as a description of all
contributions to disk evolution (an “effective” $\alpha$, see Alessi & Pudritz
2018).

The three Chambers (2019) model parameters are related to the three individual
disk evolution mechanisms in equation 1 as follows. The turbulent viscosity
$\nu$ is related to $f_{w}$ and $v_{0}$ as,

$\nu=\frac{2}{3}(1-f_{w})r_{0}v_{0}\left(\frac{r}{r_{0}}\right)^{3/2}\left(\frac{T}{T_{0}}\right)\;,$
(3)

where $T$ is the midplane temperature at radius $r$. The winds-driven velocity
through the disk is,

$v_{w}=f_{w}v_{0}\left(\frac{T}{T_{0}}\right)^{1/2}\;.$ (4)

Lastly, the surface density outflow rate depends on the parameter $K$ as,

$\dot{\Sigma}_{w}=\frac{Kf_{w}v_{0}\Sigma}{r_{0}}\left(\frac{r}{r_{0}}\right)^{-3/2}\;.$
(5)

This equation can be integrated to determine the total mass outflow rate,

$\dot{M}_{\rm{wind}}=\int_{R_{in}}^{R_{out}}2\pi r\dot{\Sigma}_{w}\rm{d}r\;,$
(6)

where $R_{in}$ and $R_{out}$ are the inner and outer disk radii, respectively.

It is important to note that the wind parameter is not freely variable, since
it is the wind that carries off a portion of the disk’s angular momentum. It
is safe to assume that over times much shorter than the disk’s viscous
evolution time scale, that the steady state version of the disk angular
momentum equation can be used. These involve the magnetic torque term of the
wind upon the disk, whose solution shows that there is a link between the mass
loss rate of the wind and the disk accretion rate (Pudritz & Norman, 1986;
Pelletier & Pudritz, 1992). We return to this point below, but for now, we
merely note that within the present formalism, the angular momentum equation
written in these variables yields (Chambers, 2019),

$K={1\over 2}\left(\frac{r}{r_{0}}\right)^{-1/4}\left(\frac{\Omega
r}{v_{esc}-\Omega r}\right)$ (7)

where $v_{esc}$ is the tangential velocity of the escaping disk wind.

A similar form of the disk evolution equation 1 was numerically solved in
Suzuki et al. (2016). While three parameters are still used in their formalism
(setting the strength of each of the three evolution mechanisms), the
individual strengths of turbulence and disk winds are instead set using the
standard $\alpha$ parameters; $\alpha_{\rm{turb}}$ and $\alpha_{\rm{wind}}$.
The following equations can be used to convert between these parameters and
the Chambers (2019) $f_{w}$ and $v_{0}$ parameters,

$\alpha_{\rm{turb}}=\frac{(1-f_{w})r_{0}v_{0}\Omega_{0}}{c_{s0}^{2}}\;,$ (8)

and,

$\alpha_{\rm{wind}}=\frac{f_{w}v_{0}}{c_{s0}}\;.$ (9)

Here, $\Omega_{0}$ is the Keplerian angular frequency at the reference radius
$r_{0}$ = 1 AU, and $c_{s0}$ is the sound speed at reference temperature
$T_{0}$ = 150 K. In viscous models using $\alpha_{\rm{turb}}=10^{-3}$ outside
the dead zone, we obtain an $\alpha_{\rm{wind}}\simeq 2\times 10^{-4}$
throughout the disk. However, even with the difference in their magnitudes,
roughly 80% of the disk angular momentum is carried by the wind in such a
model.

The midplane temperature $T$ is solved using,

$T^{4}=T_{0}^{4}\left(\frac{r}{r_{0}}\right)^{-2}+\left(\frac{3GM_{*}F}{8\pi\sigma_{\rm{sb}}r^{3}}\right)\left(\frac{3\kappa\Sigma}{8}\right)\;,$
(10)

where $G$ is the gravitational constant, $M_{*}$ is the host-star mass,
$F\simeq 3\pi\nu\Sigma$ is the mass flux due to turbulent viscosity-driven
accretion, $\sigma_{\rm{sb}}$ is the Stefan-Boltzmann constant, and $\kappa$
is the disk opacity. Here, the first term corresponds to heating from host-
star radiation, and the second to heating via viscous dissipation (Ruden &
Pollack, 1991). In the absence of viscous heating, a pure radiative
equilibrium $T_{\rm{req}}$ profile is obtained for a disk with a constant
aspect ratio $H/r$ (i.e. Chiang & Goldreich (1997)),

$T_{\rm{req}}=T_{0}\left(\frac{r}{r_{0}}\right)^{-1/2}\;,$ (11)

where $T_{0}$=150 K is the temperature at reference radius $r_{0}$ = 1 AU.

We note that, in equation 10, the viscous heating is only generated through
turbulent viscosity, and not through a general dissipation of gravitational
potential energy. In a pure winds scenario, then, there will be no viscous
heating contribution, and the disk midplane temperature will be the radiative
equilibrium profile. In this circumstance, the gravitational potential energy
lost by accreting material will be carried away by the wind. This process has
been recently investigated in Mori et al. (2019) using MHD simulations, who
find that this general heating from gravitational dissipation is small
compared to radiative or viscous heating ($\lesssim$ a 10% change to the
midplane temperature), which confirms this result of the treatment of equation
10.

The disk opacity $\kappa$ has an effect on the strength of viscous heating. As
was the case in the previously-considered pure viscous Chambers (2009) model,
the Chambers (2019) model also takes the disk opacity to be constant
$\kappa=\kappa_{0}$ throughout the majority of the disk’s radial extent. The
exception is in the innermost region above the evaporation temperature
$T_{\rm{evap}}$ = 1500 K where dust grains sublimate, and the opacity becomes
a steeply decreasing function of temperature (Stepinski, 1998). While Chambers
(2019) uses a small value of $\kappa_{0}\simeq$ 0.1 cm2 g-1 that may arise
following grain growth, we continue to use a disk opacity of $\kappa_{0}$ = 3
cm2 g-1 in accordance with our previous disk model’s treatment.

The disk accretion rate $\dot{M}_{\rm{acc}}$ is determined by calculating the
mass flux across an inner radius $R_{in}$, for which Chambers (2019) uses 0.05
AU. The disk accretion rate is therefore,

$\dot{M}_{\rm{acc}}=2\pi R_{in}\Sigma(R_{in})v(R_{in})\;,$ (12)

where the velocity of disk material at the inner radius $v(R_{in})$ scales
with the velocity at 1 AU $v_{0}$ (a model input parameter) following,

$v(R_{in})=v_{0}\left(\frac{T(R_{in})}{T_{0}}\right)^{1/2}\;,$ (13)

where we recall that $T_{0}$ = 150 K is the temperature at reference radius 1
AU solely due to radiation. However, $T(R_{in})$ is the total midplane
temperature at $R_{in}$ = 0.05 AU from combined viscous and radiative heating.

In addition to the parameters listed, one also needs to specify the initial
disk mass and radius for the model, which combine to set the characteristic
surface density. The initial disk mass can be directly input, using for
example $f_{M}=0.1$ M⊙, which is the average initial disk mass used in our
population synthesis models. The initial disk radius is handled in the
Chambers (2019) formalism through an exponential cutoff radius $r_{\rm{exp}}$
which scales with the outer disk radius. Chambers (2019) uses a setting of
$r_{\rm{exp}}$ = 15 AU that we also adopt. As the name suggests, the value of
$r_{\rm{exp}}$ indicates the radius where the surface density profile begins
to decrease sharply with further increase in $r$. As we will see, however,
significant disk surface densities can exist well outside $r_{\rm{exp}}$, so
this parameter does not immediately indicate the outer disk radius. Our
assumed value of $r_{\rm{exp}}$ = 15 AU results in an outer disk radius that
is roughly consistent with the disk radius evolution found in our previous
viscous models, using an initial disk radius of 50 AU111Since an outer disk
radius is not directly calculated in the Chambers (2019) framework, we
compared surface densities between this disk model and that of the viscous
Chambers (2009) model at its outer radius to compare the models’ radii
evolution.. Lastly, we consider all of our disk and planet formation models to
take place around a Solar mass star.

We refer the reader to section 3 of Chambers (2019) for a detailed listing of
the analytic equations that are solved (equations 34-43). While this model
provides the main framework for the disk models investigated throughout this
paper, we make the following additions.

First we can reduce the number of disk model parameters by relating the
strength of the outflow as parameterized by $K$ (i.e. equations 5 and 6) to
the disk accretion rate (equation 12) following a basic result of disk wind
theory. The mass loss rate of a disk wind torque depends upon its lever arm,
which we first briefly discuss. Consider a field line that has its foot point
at a radius $r_{o}$ on the disk. Follow the accelerating outflow along that
field line until you reach a point on it where the outflow speed equals the
Alfvén speed on that field line. This is the Alfven critical point in the
outflow and it marks the point where further acceleration basically ceases:
the field cannot enforce corotation with the disk because the outflow speed is
greater than the speed at which an Alfén wave can propagate back to the disk.
The radial distance to that point on the field line is the Alfvén radius for
that field line, $r_{A}(r_{o})$. The ratio $r_{A}/r_{o}$ is the lever arm of
the resulting torque that is exerted back on the disk. The set of all such
points marks the Alfvén surface of the disk wind - one of the most fundamental
aspects of any MHD wind. In a self-similar disk wind model such as Blandford &
Payne (1982), $r_{A}/r_{o}=const$ whose Alfvén surface is a plane parallel to
the disk midplane. Analysis of the disk angular momentum equation quickly
reveals a most useful scaling; the ratio of the wind mass loss rate to the
disk accretion rate is inversely proportional to the square of the lever arm:
$(r_{A}/r_{o})^{-2}$ (Pudritz & Norman, 1986; Pelletier & Pudritz, 1992). For
the efficient winds discussed in the Introduction, most theoretical and
numerical studies gives values for this lever arm of $r_{A}/r_{o}\simeq 3.$
(Ouyed & Pudritz, 1997; Zhu & Stone, 2018). Following the discussion of the
observations in the Introduction, we adopt this limit as the typical situation
for at least T-Tauri stage disk winds (less than 1 Myr old) and adopt
following constraint on the wind outflow rate,

$\frac{\dot{M}_{\rm{wind}}}{\dot{M}_{\rm{acc}}}=(r_{A}/r_{o})^{-2}\simeq
0.1\;.$ (14)

By solving equation 14 at time $t=0$ for a particular disk model’s
specification of $\alpha_{\rm{turb}}$ and $\alpha_{\rm{wind}}$, the constant
$K$ can be solved for as opposed to being an input parameter. This method
reduces our list of disk input parameters by one. We highlight that only low
settings of $K\lesssim 0.1$ are needed to solve equation 14. The large values
of $K=1$ that are used in example disk models in Chambers (2019) and Ogihara
et al. (2018) are in a different regime of what we call heavy disk winds than
provided by equation 14, since a setting of $K=1$ results in
$\dot{M}_{\rm{wind}}\gtrsim\dot{M}_{\rm{acc}}$. We discuss this regime in the
Appendix.

We also determine the location of the dead zone throughout disk evolution, as
resulting from Ohmic dissipation. Within the dead zone, the disk ionization
fraction is insufficient for the MRI instability to operate. To determine its
location, we have followed an approach first developed in Matsumura & Pudritz
(2003) and applied to planet populations in Alessi & Pudritz (2018). There we
showed that the best results for planetary populations arise when disk
ionization is driven by host-star X-rays. This makes physical sense because
cosmic rays will be largely scattered by shocks in the magnetized disk wind
before reaching the disk (Cleeves et al., 2013; Cleeves et al., 2015) . The
X-ray luminosity we use here is also the same in our previous work, namely
$L_{X}=10^{30}$ erg s-1 with typical X-ray energies of $E_{X}=4$ keV. This
approach results in the following criteria for an MRI-active disk, written in
terms of the Ohmic Elsasser number (Simon et al., 2013),

$\Lambda=\frac{v_{A}^{2}}{\eta_{O}\Omega}\lesssim 1\;,$ (15)

where $v_{A}$ is the Alvén speed, $\eta_{O}$ is the Ohmic diffusivity, and
$\Omega$ is the local Keplerian orbital frequency. As we will show in later
sections, this method results in the outer edge of the dead zone
$\lesssim$20-30 AU and evolving inwards with time, similar to its evolution
within the previously-considered Chambers (2009) framework (Alessi & Pudritz,
2018).

An improvement we make in our treatment of the dead zone is that, in the
Chambers (2019) framework, we can set different turbulence strengths with the
$\alpha_{\rm{turb}}$ parameter inside and outside the dead zone. For example,
following Hasegawa & Pudritz (2010), we reduce $\alpha_{\rm{turb}}$ by two
orders of magnitude inside the Ohmic dead zone, while maintaining a constant
disk accretion rate as set by the parameter $v_{0}$. What this means
physically, is that within the Ohmic dead zone, the stress related to disk
winds increases so as to maintain a radially-constant disk accretion rate
(i.e. following the result of Bai & Stone (2013)). We specify our choice of
$\alpha_{\rm{turb}}$ settings for our various disk models in section 2.4. This
treatment is an improvement of our previous handling of the dead zone, for
which the dead zone was “passive” in the sense that we determined its outer
edge’s location, but it had no physical effect on the disk structure so as to
not break the assumed self-similarity of the Chambers (2009) model. Here, the
dead zone has a more self-consistent effect on the disk structure, as its
outer edge separates two distinct regions with different strengths of
$\alpha_{\rm{turb}}$ and $\alpha_{\rm{wind}}$.

We have also investigated ambipolar diffusion (hereafter AD; another non-ideal
MHD effect) in its ability to affect the dead zone’s structure along the disk
midplane. Following Bai & Stone (2011), this investigation was done by solving
for the AD parameter Am throughout the disk,

$\rm{Am}=\frac{v_{A}^{2}}{\eta_{A}\Omega}\;,$ (16)

where $\eta_{A}$ is the ambipolar diffusivity. The parameter Am is AD’s
counterpart to the Ohmic Elsasser number, in that it quantifies how effective
AD and its diffusivity will be in suppressing MRI growth. This calculation
resulted in Am values of 100-1000 throughout the disk’s extent along the
midplane for fiducial disk settings and $\alpha_{\rm{turb}}=10^{-3}$. At this
setting of $\alpha_{\rm{turb}}$, Bai & Stone (2011) show that MRI would only
be suppressed at values of the plasma $\beta\equiv P_{\rm{gas}}/P_{\rm{B}}$ (a
ratio of gas to magnetic pressure) near 0.1-1. Therefore, AD will only
suppress the MRI and create an “AD dead zone” in the most tenuous regions of
the disk, such as in the disk’s outer extent, or well above the disk midplane.
This result is in accordance with the commonly-found conclusion that AD
affects MRI turbulence only in the lowest density regions of the disk (e.g.
Armitage 2011; Simon et al. 2013). We therefore do not include AD in our disk
models, and the dead zone location we determine is only due to the effect of
Ohmic dissipation.

###  2.3 Planet Formation & Migration

####  2.3.1 Detailed numerical treatment of planet migration under viscous
torques

The torque experienced by the planet depends on the planet’s mass as well as
important physical properties of its host disk such as its turbulence and
scale height. Our handling of migration is an extrapolation of recent detailed
torque calculations carried out in Speedie et al. (2022). This closely
followed the method developed for computing Type I torques by Paardekooper et
al. (2011) and visualized in what we call torque maps, introduced by Coleman &
Nelson (2014).

A key issue in any theoretical treatment of planet migration is how type I
migration is handled. Many models assume a smooth power-law disk model such as
the well known MMSN model. Such models cannot therefore deal with the problem
of rapid Type I migration (Goldreich & Tremaine, 1980) without resorting to
some ad hoc assumptions. The central point of physics here is that disks are
not smoothly varying power laws in column density or other quantities. They
have inhomogeneities due to transitions in turbulence levels (at dead zone
boundaries), opacities (eg ice lines where density changes occur), and in disk
heating rates (transition from viscous to radiative heating). It is well known
that co-rotation torques can push planets outwards in disks, and that
therefore there exist planet traps where the inward Lindblad torque balances
this outward co-rotation torque. These zero net torque solutions typically
occur at the inhomogeneities previously mentioned. Bitsch et al. (2015)
demonstrated that for smooth MMSN models the lack of such traps made the
formation of massive planets difficult in such models. Subsequent work by
Bitsch et al. (2019) showed that in MMSN disk models, very low viscosity
($\alpha\simeq 10^{-4}$ ) outward directed co-rotation torques saturate and
are overcome by inner directed Lindblad torque.

The fact that disks are observed to have a range of turbulence levels
motivated Speedie et al. (2022) to compute how planets grow and migrate within
disk models whose detailed evolving astrochemistry and the resulting detailed
Llindblad and co-rotation torques are carefully followed. They compute two
different cases: the conventional $\alpha=10^{-3}$, and a lower
$\alpha=10^{-4}$. In all of these simulations, a heat transition trap is
present. This trap has a far more extended radial structure and is not
localized to a specific radius. This is because there is a gradual transition
in the total thermal energy from the viscously heated inner region to the
outer radiation dominated zone. Nevertheless we can identify a fiducial radius
for the heat transition where the heating rates are formally equal. We refer
the readers to Figures 4, 7, and 8 in that paper which clearly show the torque
maps, traps, and planet evolution tracks in the M-a diagram. The important and
at first glance somewhat surprising result is that planet migration histories
bifurcate depending on the value of $\alpha$: outward directed co-rotation
torques can carry planets outwards in low viscosity disks($\alpha\simeq
10^{-4}$) where they are ultimately trapped in extended heat transition traps,
whereas forming planets in higher viscosity disks tend to be pushed inwards.
In regions beyond the extended heat transition zone, Speedie et al. (2022)
showed that the co-rotation torque does indeed saturate so that inward
directed Linblad torque force planets pushes the heat transition gradually
inwards as the disk reduces in mass and the inner viscously heating region
shrinks. Their work shows that all three traps are active but that in the
lowest viscosity case the trapping masses is quite low, dropping to several
Earth masses.

Statistical analysis of hundreds of protoplanetary disks shows that the dust
component of disks falls into two populations: (i) radially compact with no
resolved structure or (ii) radially extended with ample resolved gap-ring
structure (van der Marel & Mulders, 2021) Given that co-rotation torques in
low viscosity disks push planets outward to large disk radii before they
return to the inner regions and that planets in high viscosity disks on the
other hand, only migrate inwards,Speedie et al. (2022) noted that the level of
disk turbulence is playing a major role. Extended structured protoplanetary
disk systems with dust gaps and rings observed at 10s of au are likely disks
that have these low viscosities while the majority of disks with more compact
dust components have higher viscosity. A plausible physical reason for this is
that massive disks are better screened from their external radiation sources,
which impacts the level of MRI turbulence within them. Thus, the statistics of
disk structural properties quite naturally leads to the prediction that there
is an underlying distribution of turbulence levels in disks, directly linked
to the distribution of disk masses most probably deriving from their formation
conditions. These results and physical ideas motivates our approach in the
present paper to consider a population synthesis of planets whose host disks
have a distribution of turbulent $\alpha_{turb}$ values, that we take to be a
lognormal distribution.

####  2.3.2  A computationally advantageous approach

We follow the same approach of computing core accretion models and trapped
type-I migration as has been used in our previous works (Alessi & Pudritz,
2018; Alessi et al., 2020a, b). For a complete description, the reader may
consult Appendix B of Alessi et al. (2020a). The planet traps we include the
water ice line, heat transition, and outer edge of the dead zone. In the cases
of the ice line and heat transition, the method by which their locations are
calculated has been altered slightly in this new disk framework. In Speedie et
al. (2022) we have demonstrated that these traps are found in both low and
high viscosity disk models in detailed numerical treatments of planetary
migration that includes the evolution of disk astrochemistry. One difference
that arises in our simplified approach and these more detailed torque
simulations is that it is possible for low mass planets in some models to
leave their Type I planet trap before they have sufficient mass to open a gap
in the disk (see also Dittkrist et al. (2014) ). This is an important issue
because such a planet could undergo significant Type I migration before
reaching a mass in which a gap opens and slower Type II migration occurs.
Hasegawa (2016) has analyzed the mass scale at which the co-rotation torque
saturates - equivalent to the co-rotation mass scale derived in Speedie et al.
(2022) \- and shown that there are conditions when growing planet achieve the
saturation mass they arrive at the condition to start opening a gap and start
to enter the Type II state. We discuss this further in the Discussion section
5.

We define the ice line simply as the location along the disk midplane where
the temperature is 170 K. Our models previously computed the full disks’
equilibrium chemical structures to determine the phase-transition point of
water. We have found that in all cases the resulting ice line location
(defined where the abundance of water vapour and ice are equal) has a midplane
temperature of 170 K. While in the Chambers (2009) model the heat transition
separating viscous and radiative heating was directly calculated, in the new
Chambers (2019) model it is not. We define the heat transition at the point
where the temperature due to viscous heating $T_{\rm{vis}}=0.5\,T_{\rm{req}}$.
Since the total midplane temperature is
$T^{4}=T_{\rm{vis}}^{4}+T_{\rm{req}}^{4}$, this definition corresponds to
viscous heating contributing a $\sim$ 6% increase to the radiative equilibrium
temperature. The method of computing the dead zone’s outer edge remains the
same as our previous works, following Matsumura & Pudritz (2003) for an X-ray
ionized disk.

Details of our planet formation model remain unchanged and we continue to use
the best fit values related to forming planets’ envelope opacities found in
(Alessi & Pudritz, 2018). We consider a fiducial value of the parameter
$f_{\rm{max}}\equiv M_{\rm{max}}/M_{\rm{gap}}$ = 50 that determines the mass
at which gas accretion onto massive planets is terminated. The fiducial
setting of $f_{\rm{max}}$ applies to models where individual planet formation
tracks are computed. Following our approach in previous works, a full
$f_{\rm{max}}$ distribution, which we model as log-uniform between
$f_{\rm{max}}$ values of 1 and 500 is stochastically sampled when full planet
populations are being computed.

We simplify the dust physics by considering only a constant dust-to-gas ratio
of $f_{\rm{dtg},0}=0.01$ at Solar disk metallicity, and do not include dust
evolution effects (i.e. radial drift) in this work, whose effect on planet
populations was studied extensively in Alessi et al. (2020a). The setting of
the constant dust-to-gas ratio scales with disk metallicity (a varied
parameter in our population runs) as
$f_{\rm{dtg}}=f_{\rm{dtg},0}10^{[\rm{Fe}/\rm{H}]}$.The assumption of a
constant disk dust-to-gas ratio removes the computational expense of solving
the Birnstiel et al. (2012) dust evolution model. While radial drift will have
an effect on planet formation models and certainly adds a layer of complexity
to this problem, our goal is to focus first on understanding the basic effects
of disk evolution via winds on planet formation. In this regard, we are
following the approach we took for the turbulent disk model (Chambers, 2009)
where we first considered a constant dust-to-gas ratio model (Alessi &
Pudritz, 2018) before the more complex scenario of incorporating dust
evolution (Alessi et al., 2020a).

###  2.4 Disk Parameter Settings

We now define the various disk models and their parameter settings that we
investigate in the following results sections.

First, in section 3, we will compare disk models at two different relative
strengths of $\alpha_{\rm{turb}}$ and $\alpha_{\rm{wind}}$. In the first case,
which we refer to as the “turbulent dominated model” (or “turbulent model” for
short), we set $\alpha_{\rm{turb}}$ = 10-3 which is the same setting we have
used in our previous works, considering a pure viscous disk model. The second,
“combined turbulence & winds model” (which we refer to as the combined model,
for short) considers $\alpha_{\rm{turb}}$ = 10-4, an order of magnitude lower.
In both cases, the fiducial models are normalized using the setting of
parameter $v_{0}$ such that they have the same initial accretion rate
$\dot{M}_{0,\rm{fid}}\simeq 6\times 10^{-8}$ M⊙ yr-1. This value of
$\dot{M}_{0,\rm{fid}}$ compares reasonably with initial accretion rates
calculated using the previous disk model (Chambers, 2009), and also with
observationally-inferred accretion rates from Watson et al. (2016). Our
normalization of $\dot{M}$ between the two disk models results in them having
the same _effective_ $\alpha$,

$\alpha_{\rm{eff}}\simeq\frac{\dot{M}_{0,\rm{fid}}}{3\pi
c_{s0}H_{0}\Sigma_{0}}\;,$ (17)

where $H_{0}$ is the disk scale height, and subscripts of 0 follow the
convention of indicating temperature-dependent quantities whose values are
calculated as caused solely from radiation. The disk surface density
$\Sigma_{0}$ is set by the initial disk mass $f_{M}$ and exponential cutoff
radius $r_{\rm{exp}}$. As we will further detail later in this section, we set
these parameters such that the turbulence-dominated and combined models have
the same $\Sigma_{0}$ values, which again results in the two disk models
having the same setting of $\alpha_{\rm{eff}}$. The normalizations of
$\dot{M}$ and $\Sigma$ between the two models results in the turbulent model
having an $\alpha_{\rm{wind}}\simeq 2.5\times 10^{-4}$, and the combined model
$\alpha_{\rm{wind}}\simeq 2.7\times 10^{-4}$.

In figure 1, we summarize the turbulent-dominated and combined models by
plotting snapshots of their $\alpha$ parameters across the disks’ extents. In
both models, within the MRI dead zone, $\alpha_{\rm{turb}}$ is decreased by
two orders of magnitude. This reduction is shown clearly on the plots of
figure 1 as a transition at the outer edge of the dead zone $r_{dz}$, whose
location is shown for this example at time $t=0$ in the disk’s evolution. This
transition point, being the dead zone trap, evolves inwards with time as the
disk evolves. We see that only a small increase in $\alpha_{\rm{wind}}$ is
needed at $r_{dz}$ to maintain a radially-constant disk accretion rate despite
$\alpha_{\rm{turb}}$ decreasing by two orders of magnitude at this location.
As described above, these models have the same $\alpha_{\rm{eff}}$ parameters
due to the normalization of $\dot{M}_{0,\rm{fid}}$ and $\Sigma_{0}$.

Figure 1: We show snapshots at $t$=0 of the turbulent dominated (left) and
combined turbulence & winds (right) models by plotting the disk $\alpha$
parameters across their disks’ extents, where the decrease in
$\alpha_{\rm{turb}}$ by two orders of magnitude at the outer edge of the dead
zone ($r_{dz}$) is shown. This transition point (the dead zone trap) evolves
inwards with time throughout the disk’s evolution. We normalize both models by
using the same initial disk accretion rate $\dot{M}_{0}\simeq 6\times 10^{-8}$
M⊙ yr-1, resulting in both models having the same $\alpha_{\rm{eff}}$.

We note that in the Chambers (2019) model, quite a large disk winds fraction
$f_{w}\simeq 0.8$ is required to produce $\alpha_{\rm{turb}}=10^{-3}$. Using a
“pure” turbulence setting of $f_{w}=0$ leads to $\alpha_{\rm{turb}}\simeq
0.02$ at this value of $\dot{M}_{0,\rm{fid}}$, over an order of magnitude
larger than the setting we considered in previous chapters. Therefore, even in
the turbulent model, a substantial fraction ($\sim$ 80%) of disk accretion is
generated from disk winds and its related stress.

Other factors that affect the disks’ initial accretion rates are the settings
of the initial disk mass $f_{M}$ and radius, which scales with the exponential
cutoff radius $r_{\rm{exp}}$. Following Chambers (2019), we set $r_{\rm{exp}}$
= 15 AU for all models. We then set the initial disk mass $f_{M}$ such that
the initial surface density at reference radius 1 AU is,

$\Sigma_{0}\equiv\Sigma(r_{0}=1\,\rm{AU},t=0)=1500\,\rm{g}\,\rm{cm}^{-2}\equiv\Sigma_{0,\rm{fid}}\;.$
(18)

This value of $\Sigma_{0,\rm{fid}}$ corresponds to initial disk masses
$\simeq$ 0.05 M⊙ in both the turbulent and combined models. A surface density
of 1500 g cm-2 is similar to $\Sigma_{0}$ values in fiducial models
investigated in the previously considered Chambers (2009) framework.

In section 3, we also determine how the setting of $\Sigma_{0}$ affects
results planet formation, as a means of foreshadowing the outcomes of a full
population synthesis calculation. For both the turbulent and combined models,
in addition to the fiducial setting, we consider a high setting
$\Sigma_{0}=3\,\Sigma_{0,\rm{fid}}$ and a low setting
$\Sigma_{0}=(1/3)\,\Sigma_{0,\rm{fid}}$. While ultimately these changes are
achieved by altering the initial disk mass, we label these models in terms of
their $\Sigma_{0}$ value since the differences in surface densities can arise
from a combination of changes to the initial disk mass and radius.
Furthermore, the disk surface density is the physical parameter responsible
for setting planet formation timescales. We note that we would arrive at
similar planet formation results if we were to instead keep $f_{M}$ constant
and alter the disk radius $r_{\rm{exp}}$, provided the same values of
$\Sigma_{0}$ and $\dot{M}_{0}$ were considered. The wind outflow is
constrained according to equation 14 in all models presented in section 3.

In the Appendix A, we continue to investigate individual disk models and
planet formation tracks while shifting our focus to a winds-dominated disk
model to examine the effect of the outflow strength. In these models, we set
$\alpha_{\rm{turb}}=10^{-6}$, which is the turbulent strength within the dead
zone of the combined model in section 3. We continue to use an initial disk
mass of $f_{M}$ = 0.05 M⊙, the fiducial setting in both the turbulent and
combined models. Obtaining an initial disk accretion rate $\dot{M}_{0}=6\times
10^{-8}$ M⊙ yr-1 requires $\alpha_{\rm{wind}}\simeq 2.7\times 10^{-4}$. This
value is similar to settings of $\alpha_{\rm{wind}}$ in the previous models.
However, given the low setting of $\alpha_{\rm{turb}}=10^{-6}$, the relative
strength of disk winds is much higher.

In the first “constrained outflow” model, we follow our approach of
constraining the wind outflow parameter $K$ according to equation 3. This
constraint results in a small value of $K\simeq 0.05$. In the second model, we
do not account for the constraint provided by equation 14, and instead adopt a
high setting of $K=1$ as was used in Chambers (2019). We refer to this as the
“unconstrained outflow” model. While this model leads to very high wind
outflow rates that are comparable to, or higher than the disk accretion rate
(in contention with disk winds theory and observations), we will see that a
very efficient wind outflow has interesting effects on disk evolution, giving
rise to a surface density maximum within the disk.

Table 1: A summary of model parameters, pertaining to individual disks where planet formation tracks are computed. The turbulent-dominated and combined turbulence & winds models are investigated in results section 3, while the winds-dominated models are investigated in appendix A. Model | Outside $r_{dz}$ | Inside $r_{dz}$
---|---|---
Turbulent-Dominated | $\alpha_{\rm{turb}}=10^{-3}$ | $\alpha_{\rm{turb}}=10^{-5}$
| $\alpha_{\rm{wind}}=2.5\times 10^{-4}$ | $\alpha_{\rm{wind}}=2.8\times 10^{-4}$
Fiducial $\Sigma_{0}$ |  $\Sigma_{0}=\Sigma_{0,\rm{fid}}$ = 1500 g cm-2
High $\Sigma_{0}$ | $\Sigma_{0}=3\,\Sigma_{0,\rm{fid}}$
Low $\Sigma_{0}$ | $\Sigma_{0}=(1/3)\,\Sigma_{0,\rm{fid}}$
Combined Turbulence & Winds | $\alpha_{\rm{turb}}=10^{-4}$ | $\alpha_{\rm{turb}}=10^{-6}$
| $\alpha_{\rm{wind}}=2.7\times 10^{-4}$ | $\alpha_{\rm{wind}}=2.8\times 10^{-4}$
Fiducial $\Sigma_{0}$ |  $\Sigma_{0}=\Sigma_{0,\rm{fid}}$ = 1500 g cm-2
High $\Sigma_{0}$ | $\Sigma_{0}=3\,\Sigma_{0,\rm{fid}}$
Low $\Sigma_{0}$ | $\Sigma_{0}=(1/3)\,\Sigma_{0,\rm{fid}}$
Winds-Dominated | $\alpha_{\rm{turb}}=10^{-6}$ | No dead zone
| $\alpha_{\rm{wind}}=2.8\times 10^{-4}$ |
Fiducial $\Sigma_{0}$ |  $\Sigma_{0}=\Sigma_{0,\rm{fid}}$ = 1500 g cm-2
Constrained outflow |  $K$ from equation 14
Unconstrained outflow | $K=1$

Table 1 summarizes each of the models we investigate in section 3 and appendix
A, listing their parameters’ settings.

###  2.5 Population Synthesis Models

We first compute planet populations with constant values of the disk $\alpha$
parameters within each population. Our population synthesis models are
constructed by computing 3000 planet formation tracks (1000 within each of the
3 traps) where we stochastically sample from distributions of 4 input
parameters. Three of these input parameters pertain to the protoplanetary disk
and are set by disk formation (disk lifetime, surface density, and
metallicity) and 1 of which is a true model parameter ($f_{\rm{max}}$ which
controls the termination of gas accretion onto massive gas giants; see section
2.3). We use the same approach and distributions of our previous works (i.e.
Alessi & Pudritz 2018; Alessi et al. 2020a) with the exception of one minor
modification.

Following our convention when considering individual disk models outlined in
the previous section 2.4, we vary the disk surface density via changes in the
disk mass $f_{M}$ parameter. This could have also been accomplished by
changing the disk radius (through varying $r_{\rm{exp}}$), or a combination of
the two. However, the disk surface density is the true driver of physical
changes in disk evolution timescales and resulting planet formation. To vary
the disk surface density, we sample from a log-normal distribution. As
previously stated, a lower mean initial disk mass of 0.05 M⊙ is needed in the
Chambers (2019) disk model with our parameterization conventions to obtain an
initial $\Sigma_{0}$ value that is comparable to the fiducial setting from our
previous models using the Chambers (2009) disk (which used a mean initial disk
mass of 0.1 M⊙). Aside from this change in the mean initial disk mass, all
other distributions’ details remain the same, and can be found in section 2 of
Alessi & Pudritz (2018).

All other details of our population models in constant-$\alpha$ disks remain
the same as previously outlined in section 2.4. Specifically, our setting of
the initial disk accretion rate remains as $\dot{M}_{0,\rm{fid}}\simeq 6\times
10^{-8}$ M⊙ yr-1, with the setting of the $f_{w}$ parameter varied to obtain
the desired $\alpha_{\rm{turb}}$ setting. Finally, in all models the outflow
parameter $K$ is constrained to follow the outflow strength criterion
presented in equation 14.

We construct the distribution of $\alpha_{\rm{turb}}$ as a log-normal
distribution using a mean of $10^{-3}$. Our previous works that have used this
setting of $\alpha_{\rm{turb}}$ have found reasonable correspondence with the
observed M-a and M-R distributions, motivating our choice of the
distribution’s average. We use $\alpha_{\rm{turb}}$ values of 0.007 and
$10^{-4}$ as $\pm$1$\sigma$ of the normal distribution, respectively (that is,
the normal distribution is defined to be asymmetric). Additionally, we
restrict the range of this parameters’ setting to be between these two values.
The upper limit of 0.007 was chosen to correspond with the
$\alpha_{\rm{turb}}$ upper limit found in the TW-Hya disk in Flaherty et al.
(2018). A log normal form for the distribution of $\alpha_{\rm{turb}}$ is also
physically reasonable in that more massive disks are expected to be more
poorly ionized. Therefore we expect there to be a correspondence between the
initial disk mass distribution (which is roughly lognormal) and that of the
turbulence levels they contain (Speedie et al., 2022).

Figure 2: A population’s sampled $\alpha_{\rm{turb}}$ distribution is plotted
as a histogram along with our constructed distribution’s population density
that is sampled in our population synthesis models (dashed curve).

We show our constructed $\alpha_{\rm{turb}}$ distribution in figure 2,
plotting a histogram corresponding to a population’s sample of
$\alpha_{\rm{turb}}$ values along with the modelled distribution’s population
density. While our choice to use a log-normal within the range of [10-4,
0.007] as opposed to a log-uniform distribution may be somewhat arbitrary
given the lack of observational or theoretical support, we highlight that our
choice to make the upper and lower bounds of this range $\pm 1\sigma$,
respectively, results in $\alpha_{\rm{turb}}$ values near the average of 10-3
being sampled only slightly more frequently than values near each end of the
range. We therefore do not expect that changing to a log-uniform distribution
across this range to greatly affect our results.

Additionally, we note that in all populations, regardless of the setting of
$\alpha_{\rm{turb}}$, we maintain a $v_{0}$ parameter value of $\simeq$ 20 cm
s-1, which results in an accretion rate $\dot{M}=\dot{M}_{0,\rm{fid}}\simeq
6\times 10^{-8}$ M⊙ yr-1 at the fiducial value of surface density
$\Sigma_{0}$. However, the accretion rate will vary at different values of
disk surface density that are sampled in a population synthesis calculation.
By defining both $\alpha_{\rm{turb}}$ and $v_{0}$ for each disk sampled in a
population, $\alpha_{\rm{wind}}$ is a parameter that may be solved for but not
pre-defined in a given population. This choice is somewhat arbitrary, and an
alternate approach could be to set both $\alpha_{\rm{turb}}$ and
$\alpha_{\rm{wind}}$. However, our chosen approach allows for a transparent
means to compare with previous versions of our disk model. We note that both
$\alpha$ parameters are likely to show variance among different systems, and
the ratio of the turbulence to wind strength is an important quantity that
affects disk evolution, and therefore planet formation results. Since we are
directly investigating changes in $\alpha_{\rm{turb}}$, we will effectively be
changing the relative strength of the two drivers of disk evolution, and
therefore will be probing the effect of both parameters.

###  2.6 Planet & Atmospheric Modelling

We use the same core and atmospheric structure models as Alessi et al. (2020b)
to compute planet radii, and refer the reader to that work for a more detailed
description. A planet’s core radius will be dependent upon its solid
composition (i.e. Valencia et al. 2007), and we performed detailed solid disk
chemistry models to investigate this relationship in our previous paper
(Alessi et al., 2020b). However, when using the standard approach of binning
the solid refractories and minerals throughout the disk into bulk components
of irons, silicates, and water ice, the biggest cause of change in these
components’ compositions is their location with respect to the ice line. On
this basis, rather than computing a detailed chemistry model in this work, we
simplify this process by calculating the time-dependent location of the water
ice line (tracking the midplane temperature of 170 K), and fitting the radial
profiles of the bulk abundances of irons, silicates, and ices with a
reasonably high degree of accuracy to the profiles we obtained in Alessi et
al. (2020b) where a full equilibrium chemistry model was performed using Solar
abundances from Asplund et al. (2009). We fit to equilibrium chemistry models
corresponding to Solar C/O and Mg/Si ratios, which were investigated
separately in our previous work. The disk compositions are used to track
accreted materials onto planets throughout their formation, which is a
necessary input to the core radius calculation.

We also explored the importance of post-disk atmospheric photoevaporation on
our synthetic populations’ M-R distributions in Alessi et al. (2020b). We
continue to use the same atmospheric mass-loss model, which combines
approaches of Murray-Clay et al. (2009) and Jackson et al. (2012). We model
atmospheric mass loss as caused by photoevaporation, driven by host-stellar
X-rays and UV for up to 2 Gyr (or until convergence) following the dissipation
of the protoplanetary disk at its disk lifetime $t_{\rm{lt}}$. We did this by
combining the UV and X-ray driven models of these two papers, respectively. We
use the power law fits to measured integrated fluxes from Ribas et al. (2005)
for young solar-type stars in the X-ray (1-20 Å) and extreme ultra-violet
(EUV) (100-360 Å) wavelengths. In the early evolution of the planet X-ray
driven photoevaporation dominates due to the high x-ray fluxes from a young
star (Ribas et al., 2005; Jin et al., 2014). Mass loss is modelled by assuming
that the energy from incident photons is converted into work to remove gas
from the gravitational potential of the planet (Jackson et al., 2012). For the
extreme EUV driven regime we adopt the two regimes highlighted in the model of
Murray-Clay et al. (2009). Finally, between X-ray and UV-driven mass flows
there is an important transition in the flows. Above a certain UV flux, X-rays
are no longer able to penetrate the UV ionization front, resulting in a UV-
dominated flow (Owen & Jackson, 2012). We start the mass loss evolution of our
planets immediately after the protoplanetary disk evaporates, a parameter that
is stochastically-varied throughout our population of planets according to the
observed range of disk lifetimes. We evolve each planet forward until it is 1
Gyr old beyond which mass loss is negligible. We again refer the reader to
Appendix C in Alessi et al. (2020b) for a complete description of our
approach.

An atmospheric mass loss model is an important inclusion when considering out
populations’ M-R distributions, as the amount of gas accreted onto a planet
greatly affects its radius since it is the lightest constituent material out
of which planets form. Therefore, it is important to not only track the amount
of accreted gas during the disk phase (as we do with our planet formation
model), but also the fraction that is retained when exposed to high-energy
radiation from the host star.

##  3 Results I: Effect of relative strength of turbulence and disk winds

We begin by investigating individual disk models to compare the effects of
different relative turbulent and disk winds strengths on disk evolution and
planet formation. Specifically, we compare models at
$\alpha_{\rm{turb}}=10^{-3}$ and $\alpha_{\rm{turb}}=10^{-4}$, using otherwise
fiducial disk and model parameters that we have outlined in detail in section
2.4.

In figure 3, we compare the evolution of the turbulent-dominated and combined
disk models to discern the effect of the relative strength of turbulence and
disk winds. We plot the disks’ accretion rates, as well as radial profiles of
surface densities and midplane temperatures throughout 3 Myr of evolution at
the fiducial $\Sigma_{0}$ setting. For both models the disk accretion rate
decreases by roughly 1.5 orders of magnitude below the initial accretion rate
of 6$\times 10^{-8}$ M⊙ yr-1 after 3 Myr, decreasing slightly more in the
combined model. This range of $\dot{M}$ is quite comparable to those resulting
from the previous Chambers (2009) models with fiducial parameters.

The inhomogeneities present within the radial surface density and temperature
profiles correspond to the outer edge of the dead zone, where there is a local
change in the strength of $\alpha_{\rm{turb}}$. Comparing the two models’
surface density evolutions in the outer disk (at radii near 70-100 AU), we see
a key difference between disk evolution via turbulence and winds. In the
viscous evolution case, a small amount of spreading in the outer disk occurs -
a necessary consequence of angular momentum conservation within the viscous
evolution mechanism. Disk spreading is absent in the combined model where
winds-driven evolution is more prominent, and in this circumstance the disk,
in fact, slightly contracts. We recall that the former ‘turbulent’ model
required quite a high disk winds fraction $f_{w}\simeq 0.8$ to set
$\alpha_{\rm{turb}}=10^{-3}$. Since a small fraction ($\simeq$ 20%) of
accretion is a consequence of turbulent viscosity, there is only a small
amount of spreading in the outer disk. We note that a more extreme turbulent
model of $\alpha_{\rm{turb}}=10^{-2}$, as considered in Chambers (2019), shows
more significant amount of viscous spreading.

Figure 3: We compare protoplanetary disk evolution in the ‘turbulent’
($\alpha_{\rm{turb}}$=10-3; left column) and ‘combined’ turbulence and winds
($\alpha_{\rm{turb}}$=10-4; right column) cases. Both disks are normalized to
the fiducial setting of $\Sigma_{0}$. Time evolution of the disk accretion
rate is plotted in the top row, and profiles of surface density and midplane
temperature are plotted throughout disk evolution in the middle and bottom
rows, respectively.

While these differences between the two models exist in the outer disk regions
owing to differences in relative disk wind strength, the inner 10 AU of the
two disks are relatively comparable in terms of their surface densities. This
is the region where most planet formation occurs due to low solid accretion
rates at larger radii (with core accretion rate scaling with $\Sigma$). The
$\Sigma$ profiles of the two models are quite comparable in the inner 10 AU,
particularly for the first 2 Myr. The more rapid evolution of the combined
model becomes apparent on the surface density profile at 3 Myr, where $\Sigma$
is noticeably smaller in the combined model than in the turbulent scenario. On
the basis of comparing the two models’ (somewhat similar) surface densities,
and noting that planetary growth via solid accretion is a relatively fast
formation regime ($<$ 1 Myr), we do not expect significant differences in core
accretion rates between the two models.

When plotting the disk midplane temperature profiles in figure 3, we include
the radiative equilibrium profile (equation 11) to indicate the regions where
viscous heating is effective in both models. The turbulent and combined models
are quite different in terms of their temperature structures. Viscous heating
is more pronounced in the turbulent model due to the higher value of
$\alpha_{\rm{turb}}$, resulting in the viscously heated region (where
$T>T_{\rm{req}}$) extending to larger radii (10-20 AU at early disk evolution
times). We also see that at the location of the dead zone, where
$\alpha_{\rm{turb}}$ decreases by a factor of 100, the midplane temperatures
decrease towards the radiative equilibrium profile. This decrease is due to
the lower turbulence strength within the disks’ dead zones, which results in
less effective viscous heating. We therefore find that there are substantial
differences between the turbulent and combined models in terms of their
midplane temperatures. Since two of the traps in our model (the ice line and
heat transition) depend on the disk midplane temperature, we expect the
locations of these traps to be quite different between the two models. The
difference in the traps’ locations will also cause a significant difference
between the models’ chemical abundance profiles, which will depend most
sensitively on the disks’ temperature structures.

###  3.1 Planet traps and trap crossings

In figure 4, we explore differences in the planet traps’ locations and
evolutions between the two models at the fiducial setting of $\Sigma_{0}$. We
indeed find that the temperature-dependent traps, the heat transition and ice
line, exist at larger radii in the turbulent model due to the increased
viscous heating. Conversely, the outer edge of the dead zone lies farther out
in the combined model than in the turbulent case, despite their surface
density profiles being similar. While the midplane X-ray ionization is
sensitively dependent on the $\Sigma$ profile, there are several temperature-
dependent factors in the dead zone model (recombination rates, Ohmic
diffusivity $\eta_{O}$, etc). These factors, combined with the differences in
the disks’ temperature structures and evolution, ultimately affect the
location of the outer edge of the dead zone where the Elsasser number
$\Lambda=1$.

Figure 4: Planet traps’ evolutions are plotted for the turbulent (left) and
combined (right) disk models, both considering a fiducial setting of
$\Sigma_{0}$.

We see from figure 4 that within a typical disk lifetime of 3 Myr, all the
traps in our model converge to within $\lesssim$ 3 AU. We therefore predict
that substantial solid accretion rates can exist on planetary cores undergoing
trapped type-I migration at these traps’ locations, as high surface densities
persist in both disk models’ inner regions even after 2-3 Myr of evolution. In
the combined model, all traps exist within 10 AU after only 1 Myr of disk
evolution. Conversely, in the turbulent model, this is only the case for the
ice line and dead zone, and the heat transition evolves within 10 AU after
$\sim$ 2 Myr of disk evolution. Planet formation will be most efficient within
10-20 AU based on the disks’ surface density profiles (which set the solid
accretion timescales) that decrease sharply outside of this radius range. At
the fiducial $\Sigma_{0}$, the traps’ evolution to within this region by 1 Myr
in the combined model indicates that planet formation will be effective in
each of the three traps. In the turbulent model, we expect this to be also
true for formation in the ice line and dead zone traps. Formation within the
heat transition in the turbulent model, however, will likely result in much
longer formation timescales, due to its large radius at times up to 2 Myr.

The evolution of planet traps in the turbulent-dominated model using the new
disk framework of Chambers (2019) can be compared with trap evolution
corresponding to the pure turbulence disk model (Chambers, 2009) we previously
investigated in Alessi & Pudritz (2018)222We refer the reader to figure 5 of
Alessi & Pudritz (2018), which shows the traps’ evolution for a fiducial set
of model parameters the purely viscous disk model of Chambers (2009).. In the
case of turbulent-dominated model treated in this paper using the Chambers
(2019) formalism, the initial ordering of traps is that the heat transition is
the outermost trap at 25 AU, followed by the dead zone at 15 AU, and then the
ice line as the innermost and closest to star at around 5 AU. On the other
hand, the pure turbulence case has a different initial ordering of the outer
two planet traps, with the X-ray dead zone being the farthest out in the disk,
followed by the heat transition. In this regard, we see that, despite using
the same $\alpha_{\rm{turb}}=10^{-3}$ value, the different disk models do
change the initial locations and relative ordering of the outer two planet
traps. The explanation is that early on in these different models, the pure
turbulence disk column density at 0.1 Myr at 10 AU as an example (about 300 g
cm-1 ) is larger than that of the turbulence dominated case (about 200 g cm-1)
. The greater column depth means that the disk will have lower disk ionization
by X-rays out to larger distances leading to a more extensive region of MRI
damping and dead zone. We attribute this difference to the fact that in our
turbulent-dominated model the Chambers (2019) 80% of disk accretion is driven
by disk winds and this has the effect of reducing the column density below
what the pure turbulence driven case can do. The two disk models are similar,
however, in that the traps rather quickly evolve into the inner disk
($\lesssim$ 10 AU). In both disk models, the dead zone does so faster
($\lesssim$ 1 Myr) than the heat transition, which shifts to within 10 AU
after about 2 Myr.

The relative evolution of the different types of traps shown in figure 4 plays
a particularly important role in the character of the populations that they
produce. There are three important points to notice. First consider the
evolution of the ice line in both models. In the high viscosity model, the ice
line trap starts at about 5 AU. With time we see that it drifts slowly down to
1 to 2 AU where it stalls. The reason here is that the viscosity heats the
disk to temperatures higher than would result from just pure radiative heating
by the host star. As the column density drops so does the viscous heating rate
and the ice line moves inward to smaller disk radii. In the low viscosity
model, the iceline starts at around 3 AU due to lower viscous heating and
moves relatively less. This marks the second characteristic of trap evolution
in that the ice line become static in disk radius once viscous heating and
heating by stellar radiation are comparable - which begins precisely when the
heat transition reaches the ice line position. That asymptotic limit is at
around 1-2 AU for a solar mass star and this marks an important feature in the
M-a diagram as we shall see.

The third, and equally important aspect is the crossing of the dead zone and
ice line trap in low and high viscosity models. From the two panels in figure
4 the high viscosity case, this takes place at about 0.7 Myr in our fiducial
model, whereas for the low viscosity case it occurs later at about 2 Myr. This
has possible implications for a planet that is moving inwards while trapped on
the dead zone because it could be intercepted by the ice line and remain at
1-2 AU. This suggests that the ice line can potentially act as a filter that
prevents a larger population of Hot Jupiters from forming. We investigate this
point further in the Discussion section after the the population synthesis
results are laid out.

Figure 5 shows a few representative planet evolution tracks in the M-a
diagram. Planet formation models are calculated for each of the two models
(turbulent and combined), and at 3 settings of $\Sigma_{0}$: fiducial, high,
and low (as described in section 2.4). While our planet formation models are
calculated in long-lived 10 Myr disks, we include time marks along each planet
formation track at 1 Myr intervals. This demonstrates the effect that disk
lifetime will have on these models, as shorter lifetimes will simply truncate
the formation tracks at positions in the M-a diagram as indicated by the time
marks. This grid of planet formation models therefore contains a significant
amount of information foreshadowing the outcomes of full population synthesis
calculations, as it shows the effect of two of the most crucial varied disk
parameters; namely the disk mass and lifetime.

Figure 5: We show a grid of planet formation tracks pertaining to the
turbulent (left column) and combined (right column) disk models. We
investigate the effect of disk surface density $\Sigma_{0}$ by considering
three settings: high $\Sigma_{0}$ = 3 $\Sigma_{0,\rm{fid}}$ (top row),
fiducial $\Sigma_{0,\rm{fid}}$ (middle row), and low $\Sigma_{0}$ = 0.33
$\Sigma_{0,\rm{fid}}$ (bottom row). All models assume a long disk lifetime of
10 Myr. Open circles along the formation tracks indicate planets’ positions at
1 Myr intervals, indicating their positions if a shorter disk lifetime were to
be considered.

We first examine the fiducial $\Sigma_{0}$ case for the turbulent model. We
see that formation within each of the 3 traps gives rise to gas giants across
a wide span of orbital radii when a 10 Myr disk lifetime is used. The ice line
produces a warm gas giant at $\sim$0.6 AU after 3 Myr. Its position in the
disk within 5 AU for the entirety of the disk’s evolution results in efficient
planet formation in this trap. The dead zone forms a hot Jupiter at $\sim$0.07
AU after only 2 Myr, due to its rapid inward migration into the inner, high
density regions of the disk. The heat transition also forms a hot Jupiter at
$\sim$0.06 AU, but requires the full 10 Myr disk lifetime to do so. Disk
lifetimes of 4-9 Myr would truncate the heat transition planet’s formation
within the super Earth region of the M-a diagram between 0.06-0.4 AU. This
result confirms that the heat transition’s large radius and slow inward
migration results in the trap being relatively inefficient for planet
formation in this model.

Changing the disk surface density has two counteracting effects on planet
formation rates. Increasing the disk surface density gives rise to higher
accretion rates at a given orbital radius $r$. However, the second effect that
opposes this higher accretion rate is that increasing the surface density also
results in the positions of planet traps being shifted radially outwards to
regions with lower $\Sigma$. This shift is as a result of the higher surface
densities throughout the disk, and correspondingly larger viscous heating.
Therefore, planets forming at a given trap in models with different
$\Sigma_{0}$ values will accrete at different radii. When comparing the effect
of disk surface density on planet formation, then, one should consider the net
effect on the _local_ accretion rate at the traps’ positions.

In the case of the high $\Sigma_{0}$ turbulent model, the ice line again
produces a warm gas giant, but its final position is outside 1 AU and it forms
in only 2 Myr. The dead zone planet’s formation is similar to the fiducial
model, resulting in a hot Jupiter at roughly 0.07 AU after 2 Myr. Even after
10 Myr of formation, the heat transition is only able to produce a super Earth
at roughly 0.06 AU. In this case, the outward shift of the trap with the
increase in $\Sigma_{0}$ was sufficient to reduce its accretion such that the
trap failed to produce a gas giant. The high $\Sigma_{0}$ model is interesting
in that it produces each of the three main classes of observed exoplanets: a
super Earth, a hot Jupiter, and a warm Jupiter.

Table 2: A summary of the final planet masses and orbital radii (at the end of each disk’s 10 Myr-lifetime) in each model shown in figure 5, where we investigated the effect of different $\Sigma_{0}$ values on planet formation results in both the turbulent-dominated and combined turbulence & winds models. | Ice Line | Heat Transition | Dead Zone
---|---|---|---
Model |  $a_{p}$, $M_{p}$ |  $a_{p}$, $M_{p}$ |  $a_{p}$, $M_{p}$
Turbulence-Dominated Model
High $\Sigma_{0}$ | 1.3 AU, 729 M⊕ | 0.05 AU, 4.6 M⊕ | 0.06 AU, 414 M⊕
Fiducial $\Sigma_{0}$ | 0.5 AU, 372 M⊕ | 0.05 AU, 89 M⊕ | 0.06 AU, 411 M⊕
Low $\Sigma_{0}$ | 0.86 AU, 303 M⊕ | 0.06 AU, 90 M⊕ | 0.06 AU, 418 M⊕
Combined Turbulence & Winds Model
High $\Sigma_{0}$ | 1.9 AU, 786 M⊕ | 0.6 AU, 256 M⊕ | 0.25 AU, 165 M⊕
Fiducial $\Sigma_{0}$ | 1.2 AU 446 M⊕ | 0.55 AU, 246 M⊕ | 0.25 AU, 166 M⊕
Low $\Sigma_{0}$ | 0.9 AU, 316 M⊕ | 0.5 AU, 229 M⊕ | 0.26 AU, 170 M⊕

When examining the low $\Sigma_{0}$ turbulent scenario, we find that this
model also results in the ice line producing a warm gas giant at roughly 1 AU.
However, its formation time increases to 5 Myr. Since the ice line’s location
in the disk is relatively insensitive to $\Sigma_{0}$, the formation timescale
at this trap intuitively scales with disk surface density. This formation
track shows that an appreciable range of disk lifetimes, 1-4 Myr, will
truncate this planet’s formation and result in a super Earth or Neptune
between 1-2 AU. The decrease in $\Sigma_{0}$ from the fiducial setting has a
significant effect on the locations of the dead zone and heat transition
traps, and ultimately results in their formation timescales decreasing. The
dead zone again forms a hot Jupiter at 0.07 AU, but in this model does so
within 1 Myr. The heat transition forms a hot Jupiter at 0.06 AU in 9 Myr.
This formation track also shows that a large range of disk lifetimes, in this
case 3-8 Myr, will result in the formation of a super Earth with orbital
radius between 0.06-0.3 AU.

Tracks for the turbulent-dominated model reproduce a main result from Alessi &
Pudritz (2018) for which this model is most directly comparable; i.e.
$\alpha_{\rm{turb}}=10^{-3}$, constant $f_{\rm{dtg}}=0.01$ \- even though we
are using a different disk model. The point is that when using an X-ray
ionized disk and its related dead zone, there is a clear separation between
warm gas giants formed near 1 AU in the ice line, and hot Jupiters formed
within 0.01 AU from the dead zone. Based on the results of figure 5, we expect
that this result will remain even when a full population is considered with
distributions of disk lifetimes and $\Sigma_{0}$. With that said, it should be
recalled that this is indeed a different model, with disk winds carrying a
significant amount ($\sim$ 80%) of angular momentum.

Figure 5 shows the planet formation tracks corresponding to the combined
turbulence and winds models. We find a much more compact configuration among
gas giants that are formed in each of the three traps at all settings of
$\Sigma_{0}$ than we found in the turbulent model. In this circumstance, the
formed gas giants have the most compact orbital configuration at low
$\Sigma_{0}$, and the largest span of orbital radii is encountered in the high
$\Sigma_{0}$ case. The ice line forms a gas giant within 1 Myr at each
investigated $\Sigma_{0}$. The ice line forms gas giants at the largest
separation from the host star among each of the traps. Their radii range from
1-2 AU, shift outwards slightly as $\Sigma_{0}$ is increased. The ice line is
thus extremely efficient in gas giant formation in the combined model, which
is a result of its location within 3 AU for the entire disks’ evolution.

In contrast to the turbulent model, the heat transition produces gas giants at
$\sim$ 0.6 AU in 2-3 Myr of disk evolution in the combined model, with
shortest formation times corresponding to the low $\Sigma_{0}$ case. The final
orbital radius of the resulting gas giant is insensitive to the setting of
$\Sigma_{0}$. In the combined model, the smaller amount of viscous heating
results in the heat transition existing at a smaller orbital radius where disk
surface densities are higher. This shift to smaller orbital radii increases
the efficiency of planet formation in this trap when compared to the turbulent
model.

In the combined model, formation in the dead zone trap is also quite
different. The trap produces gas giants between 0.2-0.3 AU whose orbital radii
are also quite insensitive to $\Sigma_{0}$. These planets’ formation
timescales are between 3-4 Myr, with shortest formation timescales pertaining
to the lower $\Sigma_{0}$ cases. Comparing to the turbulent model, the dead
zone’s location in the combined scenario is at a larger radius. The trap
rapidly evolves inwards in both cases. While in the turbulent model, this trap
exclusively produced hot Jupiters within 1-2 Myr, the larger initial radius of
the dead zone in the combined scenario results in gas giants forming with
longer formation times and larger final orbital radii.

In table 2, we summarize the results of our planet formation models shown in
figure 5. We find that a larger range of planet formation outcomes are
achieved in the turbulent case when varying the disk lifetime and surface
density. Not only are the final configurations of planets more compact in the
combined scenario, but their orbital radii and formation times are quite
insensitive to variation in $\Sigma_{0}$. For example, across the range of
$\Sigma_{0}$ settings investigated, the combined model does not produce any
hot Jupiters or short-period super Earths. Additionally, we generally found
shorter gas giant formation timescales at all settings of $\Sigma_{0}$ in the
combined disk models, indicating that the super Earth population will be
limited at the lower setting of $\alpha_{\rm{turb}}=10^{-4}$. As a result, gas
giants will form too efficiently to achieve a reasonable comparison with the
data.

Varying $\Sigma_{0}$ has a larger effect on planet formation models in the
turbulent scenario, which affects planets’ formation times and/or orbital
radii. We also find that the turbulent model readily shows that the formation
of each of the observed classes of planets can be achieved. The ice line
produces warm gas giants and super Earths near 1 AU; the dead zone forms hot
Jupiters, and the heat transition results in super Earths forming over a range
of orbital radii, as well as hot Jupiters in the longest-lived disks.

Our planet formation results thereby suggest that a full population model in
the turbulent scenario with $\alpha_{\rm{turb}}=10^{-3}$ may achieve a better
comparison with the observed M-a distribution than the “combined” turbulence
and disk winds model with $\alpha_{\rm{turb}}=10^{-4}$. To fully explore this
idea, we continue our investigation using full planet populations in the
following section 4.

Finally, we refer the reader to Appendix A333The appendix for this work is
included as online supplementary material, where we show disk evolution and
planet formation results corresponding to disk models with low levels of
turbulence ($\alpha_{\rm{turb}}=10^{-6}$), where the effects of the disk wind
outflow strength are analyzed. In this investigation, we found that low
settings of the $\alpha_{\rm{turb}}=10^{-6}$ give rise to planet formation
that is entirely insensitive to changes in disk parameters such as lifetime or
mass. This indicates that, in a full population run as shown in section 4,
minimal variability in final planet masses or orbital radii would be found
even if disk parameters were stochastically sampled. Additionally, we found
that high settings of the wind outflow strength ($K$ = 1) result in
interesting effects on disk evolution, namely a local maximum in $\Sigma$ at
$\sim$ 10 AU (which was also found in Chambers 2019) giving rise to an
additional planet trap. We emphasize, however, that the high setting of the
wind outflow parameter produce large, unconstrained outflow rates
$\dot{M}_{\rm{wind}}\gtrsim\dot{M}_{\rm{acc}}$ that are in contention with the
theoretical and observational constraint provided by equation 14.

##  4 Results II: Planet Populations & the Effect of Relative Turbulence and
Disk-Winds Strength

Figure 6: A sequence of populations’ M-a distributions is shown to compare the
effect of the $\alpha_{\rm{turb}}$ value. Within each population, a constant
$\alpha_{\rm{turb}}$ values is used, for which we consider 10-4 in the left
panel, $10^{-3}$ in the middle, and 0.007 in the right.

We now turn to examining data for whole synthetic populations of planets.
Within these populations, we account for variability in planet formation
environments by stochastically varying disk lifetime, initial surface density,
and metallicity. In all of the migration tracks that we compute for planets in
the M-a diagram, we truncate our computations at an inner radius of 0.01 AU.
In so doing, we neglect the effects of the inner disk radius cutoff or the
planet trap at the inner edge of the dead zone that are known to halt further
radial inward migration at such radii (Masset et al., 2006).

In figure 6, we begin by considering populations at three different
$\alpha_{\rm{turb}}$ values, held constant within each population. We show
populations corresponding to $\alpha_{\rm{turb}}=10^{-4}$ (corresponding to
$\alpha_{\rm{wind}}=3\times 10^{-4}$) and 10-3 ($\alpha_{\rm{wind}}=3\times
10^{-4}$)444the same values used in the previous section 3 to analyze
individual disks and planet formation models as well as a high setting of
$\alpha_{\rm{turb}}=0.007$ ($\alpha_{\rm{wind}}=1\times 10^{-4}$) . The high
setting was chosen to correspond with the upper limit found from observations
of the TW-Hya disk in Flaherty et al. (2018).

The M-a distributions corresponding to the $\alpha_{\rm{turb}}=10^{-4}$ and
$10^{-3}$ populations represent the collection of individual tracks shown at
these settings of $\alpha_{\rm{turb}}$ from the previous section 3, whose
scatter comes from sampling a distribution of disk lifetimes, surface
densities, and metallicities. At the low $\alpha_{\rm{turb}}=10^{-4}$ value,
each trap populates a somewhat restricted region of the M-a space, producing
no planets within 0.1 AU. The dead zone and heat transition traps have
comparable planet distributions, forming most planets at $\sim$ 0.1-0.4 AU
with a range of masses between super Earths and gas giants. Additionally, both
traps’ distributions have a tail of super Earths extending 2 AU, with a small
amount of failed cores at yet larger $a_{p}$ and lower masses. The ice line
forms mostly warm gas giants, with a smaller fraction of Neptune-mass planets
with $a_{p}$ mostly in the 1-2 AU range. These results are very much in line
with what we found in section 3 with individual formation tracks. The M-a
distribution in the $\alpha_{\rm{turb}}=10^{-4}$ case covers a limited region
of the M-a space due to planet formation being insensitive to disk parameters;
most importantly the disk surface density.

In the $\alpha_{\rm{turb}}=10^{-3}$ population, there is a significant
increase in the diversity of the planet formation outcomes, but it remains the
case that the synthetic population covers a limited region of the M-a space
and does not have sufficient scatter to compare reasonably with the observed
M-a distribution. Again, the resulting population is a direct extension of the
results found for individual formation tracks, covered in section 3, and each
traps’ output in the full population shows similar results. The ice line
mostly produces super Earths and Neptunes between 1-2 AU, as well as gas
giants with $a_{p}$ between roughly 0.5-2 AU. The dead zone shows rapid inward
migration and planet formation timescales, resulting in almost exclusively hot
Jupiter formation. Lastly, the heat transition forms super Earths at fairly
small $a_{p}$ between roughly 0.05-0.5 AU, with a tail of failed cores
extending to lower masses and large orbital radii.

The final panel of figure 6 shows the population resulting from a high setting
of $\alpha_{\rm{turb}}=0.007$, for which individual tracks were not shown in
the previous section 3. The ice line predominantly forms super Earths between
1-2 AU with a small number of gas giants formed at 0.1-0.5 AU. The dead zone
only forms hot Jupiters, which in this population are even more massive than
in the $\alpha_{\rm{turb}}=10^{-3}$ case. Lastly, we find that the heat
transition exclusively forms failed cores at low $M_{p}$ and large $a_{p}$,
failing to even form Earth-mass planets during planet formation.

To summarize, we find that as we transition from low to high
$\alpha_{\rm{turb}}$, planet formation in the ice line and heat transition
become less efficient, while the distribution of planets from the dead zone
shifts to smaller orbital radii and larger masses. To explain this sequence,
we recall that viscous heating is reduced at lower values of
$\alpha_{\rm{turb}}$, reducing the midplane temperature. Since all traps’
positions are dependent on disk temperature, their position in the disk
changes (i.e. recall figure 4). Both the ice line and heat transition’s
positions change to smaller $a_{p}$ at lower settings of $\alpha_{\rm{turb}}$,
leading to planets accreting from regions of the disk with higher solid
surface densities, thus resulting in more efficient growth. The position of
the dead zone depends sensitively on both the disk surface density and
temperature, and was shown in figure 4 to shift to larger radii at lower
settings of $\alpha_{\rm{turb}}$, although in both cases the trap quickly
evolves to $\sim$ 0.1 AU within 1-2 Myr. This rapid inward evolution leads to
the dead zone producing massive planets at all settings of
$\alpha_{\rm{turb}}$.

The trend of the dead zone planets to shift inwards at higher
$\alpha_{\rm{turb}}$ is natural and that may explain the tendency to form more
massive hot Jupiters at smaller $a_{p}$ with increased $\alpha_{\rm{turb}}$:
planets will accrete from regions of the disk with higher surface densities.
On the other hand, there are reasons to think that building hot Jupiters at
close in distances are difficult. We have already noted that depending on the
detailed physics of the opacities in these accreting atmospheres, it may prove
too difficult to cool the gas effectively. Another consideration is that the
forming planet has a smaller Hills radius at small disk radii limiting the
core growth there substantially (Coleman et al., 2017) . We see however that
in our planet tracks the cores of planets destined to form Hot Jupiters in our
models have already reached the SuperEarth mass scale well outside of the
inner disk regions, at around 1 AU. This can be clearly seen by examining the
tracks in figure 7, which we discuss next.

Overall, we find that any individual setting of $\alpha_{\rm{turb}}$ produces
a population whose M-a distribution does not compare reasonably with the
observed distribution. In all cases, there is insufficient scatter and the
synthetic populations cover a limited region of the M-a space. It is certainly
intriguing that the populations at different settings of $\alpha_{\rm{turb}}$
cover different regions of the M-a diagram.

Figure 7: We show a time series of the varied $\alpha_{\rm{turb}}$
population’s formation tracks at 4 intermediate times leading to their final
M-a distribution at 10 Myr. Tracks are coloured based on the zone the
resulting planet populates. Data points denote the termination of a planet
formation track at the end of its disk lifetime. In the bottom right panel, we
show the final population’s M-a distribution and tracks with colour now
indicating the trap each planet formed in during the trapped type-I migration
regime. A complete animation of this time-series is included as part of this
paper’s online supplementary material: [click here]

###  4.1 Populations from a Distribution of Relative Disk Turbulence and Wind
Strength

We now compute populations whose $\alpha_{\rm{turb}}$ parameter is sampled
from a distribution. As covered in section 2.5, we construct an
$\alpha_{\rm{turb}}$ log-normal distribution ranging from 10-4 to 0.007, with
a mean of 10-3. We recall that the distribution’s upper limit coincides with
the observationally-inferred upper limit of the turbulence strength in TW-Hya
(Flaherty et al., 2018). Our model results motivate a lower-limit of $10^{-4}$
to be reasonable, as planet formation results become increasingly less
sensitive to disk parameters as $\alpha_{\rm{turb}}$ is decreased, and linking
scatter in the M-a and M-R relations to variability in disk parameters remains
the underlying motivation of using a population synthesis approach. We have
continued our investigation of planet formation in disks with low levels of
turbulence in Appendix A by considering disks with
$\alpha_{\rm{turb}}=10^{-6}$. We found that this trend persisted, with planet
formation results becoming _entirely_ insensitive to changes in disk lifetimes
and surface densities. Our approach of varying $\alpha_{\rm{turb}}$ in our
populations is also motivated by the probable variability of disk $\alpha$
parameters between different systems, and we are reflecting that variability
in our models by stochastically varying this key input parameter.

In figure 7, we show a time-sequence of formation tracks resulting from
including $\alpha_{\rm{turb}}$ in our list of stochastically-varied
parameters555The time-sequence shown includes several snapshots from a full
animation of the formation of this population, included as part of this
paper’s online supplementary material. Tracks are coloured based upon the zone
of the M-a diagram the resulting planet populates. Data points indicate the
termination of a given planet formation track at a particular planet’s disk
lifetime. Most formation tracks terminate at times near 3 Myr which is the
average of the disk lifetime distribution. At the final time of 10 Myr, all
formation tracks have terminated revealing the final planet distribution. For
clarity, we have only shown planet tracks that populate the various zones of
the M-a space (omitting, for example, formation tracks of sub-Earth-mass
cores).

Figure 8: The final M-a distribution of the varied $\alpha_{\rm{turb}}$
population is shown directly after the disk phase (i.e. before
photoevaporation in considered) in the left panel, and after up to 2 Gyr of
photoevaporation in the right panel. We include the frequencies by which
planets populate various zones of the M-a space.

Lastly, at the final time in the sequence presented in figure 7, we also show
formation tracks with colour indicating the trap each planet formed in during
their trapped type-I migration phase. This reveals that hot Jupiters arise
from the dead zone, warm Jupiters from the ice line, and super Earths from
each of the three traps. However, the small $a_{p}\lesssim$ 1 AU portion of
the super Earth distribution is predominantly formed in the heat transition,
while super Earths at larger orbital radii have formed in the ice line.

An important result here is that many of the hot Jupiters form quite rapidly,
reaching the zone 1 region of the M-a space prior to 1 Myr. These tracks
result from planet formation in the dead zone, whose evolution brings planets
into the inner regions of disk quickly, allowing them to accrete from high-
density regions leading to rapid formation. Many warm gas giants near 1 AU
also form quickly, accreting substantial material within 1 Myr. These tracks
correspond to formation in the ice line, showing that the trap efficiently
forms planets due to its location in the inner, high surface density region of
the disk at early times in disk evolution ($\lesssim$ 3 AU). Super Earth
formation tracks terminate at a range of times from intermediate ($\sim$ 3
Myr) to long disk lifetimes near 10 Myr.

We have argued in earlier works that Super Earths are likely to be failed
cores in the sense that their formation timescales exceed their disk lifetime
(Alessi et al., 2017). Figure 7 reveals that this can arise by either short
disk lifetimes (in the case of terminating formation by 3 Myr), or notably
long formation times (for which the disk lifetime is long, over several Myr).
Many tracks that populate the super Earth region of the M-a space with long
formation timescales correspond to planet formation in the heat transition
trap. Since the heat transition resides outside of 5-10 AU for the first 1-2
Myr of disk evolution, solid accretion rates onto cores remain small at these
times leading to inefficient planet formation (not producing many gas giants).

###  4.2 General Properties of Planet Populations in the M-a Diagram

In figure 8, we show the varied $\alpha_{\rm{turb}}$ population’s M-a
distribution both before (immediately following formation) and after
photoevaporation. We see that the resulting populations range over a
substantial portion of the M-a space. When compared with figure 6 with its
constant settings of $\alpha_{\rm{turb}}$, we see significantly more scatter
and an M-a distribution with more correspondence to observations. It is also
clear that photoevaporation does not remove much mass from planets which is
why the two panels in figure 8 look so similar. It is well known that
photoevaporation contributes mostly to changing the radius of close in
SuperEarths, not their masses.

A striking feature in this M-a diagram is that the super Earth and Neptune
planets are of decidedly mixed origins, with contributions from each of the
three traps. We also see direct formation of short period super Earths within
$a_{p}\lesssim 0.03$ AU from the protoplanetary disk phase. Our previous
population synthesis models in the pure turbulent disk scenario was limited in
its ability to directly form short-period super Earths (Alessi et al., 2020b).
However, in the combined turbulence and disk winds scenario, we find super
Earths with quite small $a_{p}$ forming directly from the protoplanetary disk.

Figure 9: The super Earth and Neptune region (zone 5) of the M-a space is
shown for our varied $\alpha_{\rm{turb}}$ population before (left) and after
(right) photoevaporation. Planets’ colours denote their ice compositions, the
main indicator of their overall solid composition. Datapoint shapes indicate
the trap each planet formed in, or in the after photoevaporation case, the
small number of planets that formed as zone 1 sub-Saturns and lost atmospheric
mass to produce a short-period Neptune.

The varied $\alpha_{\rm{turb}}$ population’s gas giant distribution consists
of a large number of hot Jupiters formed within the dead zone, and a large
number of warm gas giants formed in the ice line between 0.5-2 AU, with a
reduction in gas giant frequency between roughly 0.1-0.5 AU. This reduction in
gas giant frequency at these orbital radii corresponds with the reduced gas
giant frequency found near orbital periods of 10 days in occurrence rate
studies (Santerne et al., 2016; Petigura et al., 2018). This also reproduces
the result we found in Alessi & Pudritz (2018) where our populations showed a
clear separation between the hot and warm Jupiter populations. However,
occurrence rate studies also suggest that the gas giant frequency should have
an overall increasing trend with increasing orbital radii, reaching a maximum
between 3-10 AU (Cumming et al., 2008). There is tension between our varied
$\alpha_{\rm{turb}}$ population’s M-a distribution and the observed trend in
gas giant frequency in two aspects: 1) the population forms too many hot
Jupiters due to efficient formation in the dead zone trap; and 2) the warm
Jupiter distribution formed in the ice line sharply falls off at $a_{p}>$ 2
AU, instead of extending to larger radii in the 3-10 AU range.

Each trap leaves different and somewhat distinct imprints on the M-a diagram.
The dead zone efficiently forms hot Jupiters with $a_{p}\lesssim$ 0.1 AU as
well as some Neptunes with $a_{p}\lesssim$ 0.5 AU. The trap quickly evolves
into the high density, inner regions of the disk leading to rapid formation of
massive planets. The ice line forms many warm gas giants, as well as super
Earths and Neptunes mostly in the 0.5-2 AU range. This trap is situated within
5 AU for the entirety of the disk’s evolution, leading to rapid growth of
cores as they accrete from a high surface density region of the disk. The heat
transition has the most inefficient planet formation as it resides in the
outer disk for the majority of evolution. It contributes a substantial amount
of super Earths between 0.1-0.6 AU, with a tail of sub-Earth mass cores
extending to larger orbital radii that fall below observational limits.

Photoevaporation has the greatest effect on planets at smaller orbital radii
(as it is a flux-driven effect) and also for planets with lower surface
gravity as they are more susceptible to mass loss. Comparing the zone 5 M-a
distribution before and after photoevaporation, we indeed see that planet
masses reduce within $a_{p}\lesssim$ 0.2 AU. In particular, planets that
formed in the heat transition show the most noticeable mass reduction. This is
because the heat transition planets form with lower core masses due to their
delayed growth compared to planets in the ice line and dead zone traps. This
leaves more atmospheric mass available for stripping, leading to a greater
overall mass difference after photoevaporation. Many of the short-period super
Earths that form in the dead zone are mostly/entirely stripped, but their mass
difference is not as noticeable due to their higher core masses (and lower
atmospheric mass fractions). Our analysis of this population’s M-R
distribution in what follows will provide more insight into this argument.

Figure 10: We show the varied $\alpha_{\rm{turb}}$ population’s M-R
distribution before (left) and after photoevaporation (right). Data point
colours indicate planets’ orbital radii, and the shapes denote the trap they
formed in. We also indicate, in the right panel, planets that populate zone 1
immediately following the disk phase, but experience mass loss leading to
Neptune-mass planets. Photoevaporation improves our comparison with the
observed data (plotted in grey with error bars), particularly for planets
whose radii are heavily increased due to accreted gas ($M_{p}\gtrsim$ 2 M⊕).
Finally, we show M-R contours corresponding to three different uniform solid
compositions: ice rich planets (taking ice line core compositions from Alessi
et al. 2020b), dry, Earth-like compositions (using dead zone core
compositions), and pure iron cores.

###  4.3 Planetary Compositions and the M-R diagram for Populations

Figure 9 shows the ice fractions of planets throughout zone 5 as computed via
disk equilibrium chemistry models. As we found in Alessi et al. (2020b), a
planet’s ice composition is the main indicator of its overall solid
composition, with only small variations in the iron-to-silicate ratio existing
throughout the disk. We find a large variety in ice abundances on super Earths
and Neptunes as is expected given their varied origins. Planets that form in
the ice line have roughly one third of their solid mass comprised of ice. Both
the dead zone and heat transition traps form a mix of dry and rocky super
Earths (Earth-like compositions), and ice-rich planets depending on where the
planet accretes with respect to the ice line. Many of the planets forming in
these two traps accrete solids inside the ice line leading to dry planets,
however a portion accrete outside and can achieve up to $\sim$ 50% of their
solid mass in ice. The position of the traps relative to the ice line at the
timing of solid accretion is dependent on the exact disk parameters in the
population.

Our composition results found in zone 5 planets differ from those found in our
previous work that considered the pure turbulent disk scenario, Alessi et al.
(2020b), in three ways. First, the dead zone planets here show a range of
solid compositions from dry and rocky to a smaller portion with substantial
ice contents, whereas in our previous work the trap exclusively formed planets
with nearly zero ice content. We attribute the difference in results to a
difference in the trap’s evolution in the two disk models. Second, we find a
larger variety of super Earth solid compositions in this work’s population
since zone 5 has planets formed in each of the three traps. In our previous
work, depending on the initial disk radius considered, at most two of the
three traps contributed a substantial portion of super Earths. Lastly, while
there is an overall trend for dry planets to exist at small $a_{p}\lesssim$
0.5 AU, and ice-rich planets at $a_{p}\sim$ 0.5-2AU, there is not a clear
division in planet composition at any $a_{p}$ since the heat transition and
dead zone form planets with a range of compositions, and the traps overlapping
regions of zone 5. In our previous work, we found clear divisions between
super Earths that formed in different traps, and therefore having markedly
different compositions.

The M-R distribution for the varied $\alpha_{\rm{turb}}$ population in figure
10, compares the distributions before and after atmospheric photoevaporation.
Planets are coloured to show their orbital radii, a key parameter to consider
when examining photoevaporation’s effect. We overlay our computed populations
with the observed exoplanet M-R distribution. The data set we have used is the
same as that used in Alessi et al. (2020b) which is for planets around solar
mass stars and extracted from the NASA Exoplanet Data Base in March 2020. To
guide our analysis, we also plot three M-R relations corresponding to three
different uniform solid compositions: an average solid composition found in
the ice line, which is an ice-rich planet having roughly a third of its solid
mass in ice, a dry dead zone planet composition which corresponds to a rocky
Earth-like planet, and lastly a pure iron core. These contours show the range
of planet radii a planet can achieve from variations in solid compositions
alone. They are useful in our analysis of our M-R distributions as they
indicate planets with little/no atmospheres which will lie near or within the
contours. Planets significantly above the M-R relations indicate planets with
some atmospheric mass. For this point, we note that since gas is the lightest
constituent material out of which planets are composed, only a small amount of
gas can significantly increase a planet’s radius (Lopez & Fortney, 2013).

Our computed populations’ M-R distribution before photoevaporation shows a
division between core-dominated planets and planets whose radii are
significantly influenced by gas at roughly 2 M⊕. Planets at higher masses have
accreted sufficient gas to lie significantly above the solid M-R contours. At
$M_{p}\gtrsim$ 2 M⊕, a large portion of our computed M-R distribution compares
reasonably with the upper portion of the observed M-R distribution. However, a
fraction of these planets, many of which formed in the heat transition, have
radii that lie significantly above the observed distribution. These planets
have particularly high atmospheric mass fractions, and lower core masses than
planets formed in the other two traps at comparable total planet masses. Our
M-R distribution before photoevaporation has many planets that are at large
radii for a given mass in comparison to the observed distribution, which is a
reproduced result from Alessi et al. (2020b). We note that all of these
planets form as planets whose formation was terminated during their slow gas
accretion phase, so it follows that their radii are somewhat influenced by
atmospheres. Since we are considering our M-R distribution as formed from the
disk, we see that fairly massive atmospheres can be accreted in the absence of
photoevaporatiive losses. As gas accretion is extremely sensitive to the
setting of atmospheric opacity during formation (Mordasini et al., 2014;
Hasegawa & Pudritz, 2014), our M-R distribution may differ were we to explore
envelope opacity values other than those we found in Alessi et al. (2017).

Photoevaporation will reduce the radii determined by transits for close in
planets. We quantify this by plotting the final M-R distribution in the right
panel of figure 10. Photoevaporation indeed improves our comparison to the
observed distribution in several ways. First, the planets that formed in the
heat transition whose radii lied well above the observed distribution have
their radii reduced to compare reasonably with the upper portion of the
observed M-R distribution. These planets have rather low core masses and are
more susceptible to atmospheric mass-loss due to their lower surface
gravities. Planets with orbital radii $\lesssim$ 0.05-0.1 AU are significantly
or entirely stripped due to their proximities to their host stars and large
XUV-flux received, while planets at larger $a_{p}$ have unaffected radii. We
see that for all masses up to 30 M⊕, there are planets at small $a_{p}$ that
are entirely stripped as shown by their radii coinciding with the solid-core
M-R contours. At the higher $M_{p}$ end of the M-R distribution, there are a
small number of stripped cores that originally formed as zone 1 planets (sub-
Saturns). We note that many of these planets whose atmospheres are entirely
stripped have quite large core masses, explaining why they show significant
radius differences due to photoevaporation while only showing modest
differences in planet masses (i.e. figure 9).

Photoevaporation is therefore important in understanding the fate of gas
accreted by planets that migrate too close to their host stars. However, its
results are largely bimodal in the sense that planets for the most part either
retain their atmospheres or are fully stripped. This is seen in our
population’s distribution as planets distribute either along the solid M-R
contours (stripped planets) or along the high radius distribution that
coincides with the upper portion of the observed distribution. However, our
models are limited in their ability to compare with the portion of the M-R
distribution that exists at intermediate radii for a given mass. We find very
few planets that, after photoevaporation, achieve radii in this intermediate
state (between fully stripped or retained atmospheres). Finally, the effect of
planets’ solid compositions on the M-R diagram can be seen among fully
stripped planets, whose radii slightly differ depending upon their
compositions. However, these radius differences are small and are within
typical radius observational uncertainties.

##  5 Discussion & Conclusions

We have investigated the effects of MRI-turbulence and MHD disk winds as well
as a range of initial disk conditions on disk evolution, and the properties
and compositions of the resulting planetary populations. The relative strength
of turbulence and disk winds is a critical factor in these models. We first
considered various constant settings of turbulence and wind strengths
initially, focusing simulations on just three characteristic values for
$\alpha_{\rm{turb}}$, ($10^{-4}$, $10^{-3}$, and 0.007) that cover the range
of values quoted in the observational literature. Disk winds in our model are
essential to carry off the angular momentum from the inner dead zone region of
the disk, as recent simulations have clearly shown (Bai & Stone, 2013; Gressel
et al., 2015). Even at a setting of $\alpha_{\rm{turb}}=10^{-3}$, over 80% of
the disk accretion is being driven by winds as opposed to turbulence in our
models (since typical settings of $f_{w}\gtrsim 0.8$). We found that each of
these constant $\alpha_{\rm{turb}}$ models produced only a rather limited
variety planets populating the M-a diagram even though we included the full
observationally-constrained distributions of initial disk surface densities
and lifetimes in accordance with previous population synthesis efforts (i.e.
Ida & Lin 2004b; Mordasini et al. 2009; Hasegawa & Pudritz 2012). While all of
the constant $\alpha_{\rm{turb}}$ populations had insufficient scatter to
compare reasonably to the observed M-a distribution, the
$\alpha_{\rm{turb}}=10^{-3}$ population showed the most promising results as
it produced the widest variety of outcomes in terms of planet classes formed.
At this setting, super Earths were formed over a range of orbital radii
between $\sim$0.05-2 AU, and there was a separation between warm Jupiters at
$a_{p}\sim$ 0.6-2 AU formed in the ice line, and hot Jupiters formed in the
dead zone. Both higher and lower constant settings of $\alpha_{\rm{turb}}$
resulted in populations that were much more confined in terms of their M-a
distributions.

This result reinforces our argument that if the scatter in populations is
significantly linked to initial conditions for the disk population, then a
single constant values of turbulence strengths do not well represent the
physics of turbulence across a whole population of initial disks. We therefore
carried out simulations by adding a turbulence-strength parameter,
$\alpha_{\rm{turb}}$, to our set of stochastically-varied disk parameters in
our population synthesis calculations. We chose an average value of
$\alpha_{\rm{turb}}=10^{-3}$ for the distribution that was sampled in the