varied $\alpha_{\rm{turb}}$ population. It is interesting that this optimal
value of $\alpha_{\rm{turb}}$ is in agreement with that used in our previous
population synthesis works approached within a pure viscous disk framework.
Perhaps this is a consequence of our normalization of the fiducial disk
surface density $\Sigma_{0}$ and initial accretion rate $\dot{M}_{0}$,
recalling that we set these initial disk parameters in accordance with our
fiducial viscous disk from Alessi et al. (2020a) in order to directly compare
the disk models’ results.

We found a drastic improvement in the final planetary M-a distribution when a
distribution of.$\alpha_{\rm{turb}}$ was used. As different individual
settings of $\alpha_{\rm{turb}}$ were shown to populate nearly distinct
regions of the M-a diagram (i.e. figure 6), sampling this parameter from a
distribution introduced a significant source of scatter in the final planetary
distribution. The “varied $\alpha_{\rm{turb}}$” population formed a
substantial amount of each of the observed planetary classes: warm Jupiters,
hot Jupiters, and super Earths. The super Earth population now extends to
small orbital radii $\lesssim 0.03$ AU, a region of the M-a diagram that our
previous models in the pure viscous framework did not populate. Additionally,
the varied $\alpha_{\rm{turb}}$ population readily produces a separation
between the hot and warm Jupiter populations, in line with occurrence rate
studies (Petigura et al., 2018).

###  5.1 Comparison with observed populations in M-a diagram

These results demonstrate that the relative strength of turbulence and disk
winds could be one of the key features that shape the outcomes of planet
formation and planetary populations. This is especially noteworthy as disk
models frequently prescribe a single value for wind or viscous strengths. Here
we address the successes and caveats regarding comparison of our models with
the data in the M-a diagram.

Although the current data is still somewhat limited, it is to be expected that
disks will exhibit a _range_ of turbulence and disk winds strengths. One
reason for this is that more massive disks with their higher column densities,
are better able screen ionizing X-rays over greater portions of the disk,
which in turn will reduce MRI turbulence. Thus there should be a correlation
between MRI driven turbulence amplitudes and disk mass, as has been argued by
Speedie et al. (2022). Our results indicate that the distribution of their
relative strengths will have crucial effects on planet formation. Even within
an individual disk, there is no theoretical basis for having $\alpha$
parameters that are constant with disk radius or time. In fact, numerical MHD
results show that the disk $\alpha$ parameters can vary within disks’ radial
and vertical extents in a non-trivial manner (i.e. Bai & Stone 2011; Lesur et
al. 2014; Gressel et al. 2020).

Since all of our previous population synthesis investigations within the
viscous regime adopted a constant value of $\alpha_{\rm{turb}}=10^{-3}$
(Alessi & Pudritz, 2018; Alessi et al., 2020a), this raises the important
question of why the pure viscous models resulted in scatter in the final M-a
distribution reasonably in line with the data, while constant
$\alpha_{\rm{turb}}$ in this combined winds and turbulence model do not. We
attribute this to the traps in this combined model being less sensitive to the
initial disk parameters, most notably the initial disk surface density
$\Sigma_{0}$. Since the location of the early stages of planet formation is
set by the position of traps in the type-I migration regime, the traps’
locations in large part shapes the radius distribution of the final planet
population in our models.

The main tension between our results and the M-a diagram distributions is the
overpopulation of the Hot Jupiter region (by a factor of 2) and an
underreprsentation of SuperEarths. Specifically comparing populations shown in
Figure 8b in this paper with the observational statistics quoted in Alessi et
al. (2020b) we find: Hot Jupiters 23.6 % vs 11.5 %; Warm Jupiters 19.4 % vs.
14.9 %; and SuperEarths and Mini Neptunes 47.3 % vs. 62.2 %. Clearly the
overproduction of Hot Jupiters in our simulations is at the cost of the
SuperEarth sector, as we have already indicated. We attribute this to
efficient planet formation in the dead zone which is a result of the trapid
migration of the trap to the inner disk where surface densities are high. What
are the caveats to this picture?

One technical point is our neglect of dust evolution via radial drift, as we
did in our pure viscous simulations in Alessi et al. (2020a). In the absence
of dust traps this process rapidly depletes solids outside the ice line with
the consequence that solid accretion outside the ice line occurred at a very
slow rate. Formation in the dead zone trap is particularly interesting as
planets experience delayed growth: solid accretion only begins at an
appreciable rate after the trap migrates inside the ice line. We have shown
that this occurs in the new disk model after $\sim$ 1 Myr (recall figure 4).
Radial drift would therefore have the effect of reducing the frequency of the
dead zone trap’s formation of hot Jupiters, with more planets ending up as
super Earths and Neptunes due to this delayed growth effect. Otherwise, we do
not expect radial drift to greatly affect our results of the varied
$\alpha_{\rm{turb}}$ population. As an example, when dust evolution was
included, the ice line was a prime location of efficient planet formation,
resulting in a mix of super Earths and warm Jupiters. These are precisely the
planet classes that are formed in this trap when radial drift is not included.

A plausible reason for the noted differences in the relative size of our
computed hot Jupiter and SuperEarth populations the observations is related to
our neglect of possible planet transfers during trap crossings - most notably
where the dead zone trap crosses the ice line trap on its way to the inner
part of the disk as its column density falls. As previously noted, the ice
line could intercept and capture the lower mass planets moving inwards with
the dead zone trap. The reason is that the corotation mass scale (which is the
characteristic mass at which the outward co-rotation torque peaks) is higher
in the active region beyond the dead zone - which the ice line enters for the
first time. While it is beyond this already overly long paper to compute this
effect, we suggest that since the co-rotation mass scale at which the co-
rotation torque peaks scales as $M_{corot}\propto\alpha^{2/3}h^{7/3}$ (Speedie
et al., 2022), lower mass planets moving with the dead zone would be trapped
by the ice line as this crossing occurs. Only the most massive planets on the
dead zone trap would move inside the ice line. This is likely to reduce the
Hot Jupiter population substantially, because we have seen that the dead zone
is the primary source for the hot Jupiter population. At the same time, this
would effectively increase the SuperEarth population particularly in the
region out towards 1 AU where the ice line trap lingers. We will need to
compute this effect in another paper but suggest that this may be the most
likely correction that would restore the correct balance between the
SuperEarth population statistics.

Planet formation near the heat transition trap is very inefficient due to its
location outside of the ice line for the majority of disks’ lifetimes. Out
there radial drift has depleted solids. This causes the trap to form only
failed cores, with the most massive planets only becoming super Earths. This
is also quite similar to the heat transition’s output in the varied
$\alpha_{\rm{turb}}$ population, although the heat transition does produce an
appreciable number of super Earths between 0.1-1 AU. We expect the formation
frequency of super Earths from the heat transition would be reduced if dust
evolution was included in our models.

The relative amount of radiative versus viscous heating is important in
understanding disk and planet evolution. We recall that the position of all
three of the traps depends on the midplane temperature structure - either
directly in the case of the ice line and heat transition, or indirectly in the
case of the dead zone outer edge. Viscous heating of the disk is only
contributed by MRI-turbulence, while the radiative heating profile is constant
and independent of disk $\Sigma$ where radiative heating dominates. As a
result, only the fraction of disk evolution driven by turbulence translates to
changes in viscous heating with $\Sigma$. Ultimately this results in the trap
positions being less sensitive to changes in $\Sigma$ in a population than
they otherwise would be in a pure viscous model and hence to less scatter in
the final M-a distribution. This is is another physical reason why the
variation of the $\alpha_{\rm{turb}}$ parameter is important and therefore
plays an significant role in the properties of populations in the combined
winds and turbulence disk picture.

The use of just single planet evolution in any disk leads to potential
limitations in understanding observed populations. Multiple planet can affect
migration and accretion while still in the disk phase and planet-planet
scattering effects will also play a role in the post disk phase of the
evolution of a planetary system. Multiple planet population forming in gas
disks may enter into stable near mean-motion resonances with closer in planets
which effectively reduces their rapid inward migration (Kley & Nelson, 2012).
Recent work has shown that such effects take place, but that even with
multiple planet effects included, differences are not typically greater than
factors of two in final orbital radii (Emsenhuber et al., 2021). After the gas
is dispersed, we first note that a single planet-planet scattering of equal
mass planets in which one is removed from the system, would reduce the orbital
radius by a factor of two or less (Nagasawa et al., 2008). Outward scattering
resulting from planet-planet interactions in the interior of the planetary
system (within 10 AU) to distances of the order 100 AU would be difficult to
achieve. These and other dynamical effects arising from the inward scattering
of planets by the Kozai effect need to be considered (Bryan et al., 2016) need
to be considered, but our point remains that variable viscosities are
important.

In extending our current model into multiplanet simulations as performed by
Emsenhuber et al. (2021), we note that planets on traps could anchor resonant
chains with other forming planets. The fact that the dead zone trap ultimately
becomes becomes the innermost of the three traps that arise in disk interiors
and that it would only trap the lowest mass planets may have important
implications for the origin of densely packed planetary systems.

###  5.2 Comparison with oberved populations in M- R diagram

The varied $\alpha_{\rm{turb}}$ population also gave rise to a M-R
distribution that, after accounting for post-disk atmospheric mass loss driven
by photoevaporation, corresponds reasonably to the observed distribution.
Atmospheric photoevaporation leads to a clear improvement in our population’s
M-R distribution over that which arises directly after the disk phase. In
particular, a fraction of super Earths form in our models with quite large
radii due to accreted gas, and lie above the observed data. Atmospheric
photoevaporation reduces all of these planets’ atmospheres to better
correspond with the data.

We find that planets at masses $\lesssim$ 3 M⊕ have radii that compare well
with the data, while super Earths and Neptunes at higher masses compare with
the high-radius end of the observed data at a given mass. This is similar to
our previous work, Alessi et al. (2020b), that investigated populations’ M-R
distributions in the pure viscous framework. However, we do achieve a somewhat
better M-R distribution than what we found in the pure viscous scenario due to
the varied $\alpha_{\rm{turb}}$ population forming more super Earths at small
orbital radii $\lesssim 0.1$ AU, which in turn leads to more atmospheric mass
loss from photoevaporation. We therefore find a larger fraction of planets
with masses $\gtrsim$ 3 M⊕ being partially or completely stripped, evolving
their positions on the M-R diagram that better compare with average observed
exoplanet radii for a given mass.

However, a large fraction of super Earths do form at yet larger orbital radii
where photoevaporation has a less-severe effect, and these planets remain at
somewhat large radii and at the upper envelope of the observational data in
the M-R diagram. As we found in Alessi et al. (2020b), photoevaporation only
has significant effects over a somewhat limited range of orbital radii
$\lesssim$ 0.1 AU. This limits the mechanism’s ability to improve our
populations’ M-R distribution, particularly at larger super Earth masses where
planets in our model tend to accrete substantial atmospheres, leading to
larger planetary radii than the bulk of the data. An additional means of
atmospheric mass loss that has significant effects over a larger range of
orbital radii would be advantageous for this purpose. This perhaps could be
the case for core-powered atmospheric mass loss, however both photoevaporation
and core-powered mass loss have only been investigated for planets at short
orbital periods of less than 1 year (Owen & Wu, 2013; Gupta & Schlichting,
2019).

###  5.3 Comparison with other disk wind results

One of the important aspects of disk winds for planet formation is that their
efficient transport of disk angular momentum drives the bulk of disk accretion
flow so that high levels of turbulence are not required in contrast with the
standard viscous disk picture. Low levels of turbulence that provide minimal
resistance to efficient dust settling to the disk midplane (Hasegawa et al.,
2017). However, as we have noted, our population results suggest that some
levels of turbulence are likely present and contribute to the diversity of
planetary populations. A number of models of disk winds and planet formation
have explored various regimes including a mix of turbulence and disk winds,
heavily mass loaded vs light disk winds, much lower levels of turbulence
leading,, or no turbulence at all. It is interesting to understand the
differences that these regimes might have in the kind of populations that
would arise.

First we explore in the Appendix A our own experiments with very low levels of
$\alpha_{\rm{turb}}$ as well as cases in which very inefficient, heavy disk
winds dominate; $K\simeq 1$. Using a value of $\alpha_{\rm{turb}}=10^{-6}$
(see Appendix) as typical of very low levels of turbulence, we find that the
disparity between the synthetic populations at constant $\alpha_{\rm{turb}}$
and the observed M-a distribution was exacerbated. In this regime of extremely
low turbulence, the planet formation tracks were entirely insensitive to disk
parameters. This result persisted regardless of whether we used the
constrained disk outflow ($\dot{M}_{\rm{wind}}/\dot{M}_{\rm{acc}}\simeq$ 0.1;
the fiducial setting), or relaxed this criterion to allow for a strong wind-
driven outflow as is often prescribed for winds-driven disk models. High winds
mass-loss rates resulted in a rapidly-decreasing disk surface density on
timescales shorter than typical core accretion timescales, producing failed
cores irrespective of disk surface density or lifetime. Our results indicate
that winds-dominated scenarios alone have difficulties in explaining the
diversity of planetary masses and orbital radii observed on the M-a diagram -
if indeed this is one of the major sources of scatter. However, this pertains
only to a trapped type-I migration scenario, and it may be the case that
laminar disks do not require trapping as we have prescribed.

The idea that planets can be prevented from undergoing rapid Type I migration
has been the focus of many investigations for more than a decade (Hasegawa &
Pudritz, 2011; Dittkrist et al., 2014; Bitsch et al., 2015; Coleman & Nelson,
2016; Bitsch et al., 2019; Speedie et al., 2022). In this type of approach
where type-I migration torques are applied to directly determine migration of
forming cores, the presence of traps (null torque locations) throughout the
disk can cause cores to converge to the traps’ locations. Type I migration can
be significantly reduced if a planet opens a gap in the disk, leading to Type
II migration. Recent work shows that the gap opening mass is given by the
traditional analysis (Lin & Papaloizou, 1986; Rafikov, 2002) wherein the
angular momentum flux of spiral waves excited by the planet is balanced by the
viscous torque. This has led to a new gap opening criterion for the ratio of
planet mass to stellar mass, q that opens a gap: $q\geq 5\alpha^{1/2}h^{5/2}$
(Kanagawa et al., 2018). For a low viscosity disk with
$\alpha_{\rm{turb}}=10^{-4}$ and an aspect ratio of $h=0.05$ this means that a
Neptunian mass can open a gap. Moreover TypeI migration is slowed because it
is driven by the gas surface density in the gap rather than by the unperturbed
disk. Population synthesis studies suggest that this can be accomplished by a
two - alpha disk model in which the bulk of the accretion flow is driven by
the wind, leaving a low value of the turbulence that greatly limits Type I
migration - (Ida et al., 2018; Matsumura et al., 2021) . This leads to good
agreement with the population of warm Jupiters.

Our simulations assume that planets as they become more massive in their traps
will enter conditions in which they begin to open a gap and undergo Type II
migrations. The co-rotation torque for sufficiently massive planets will
saturate somewhat below the gap opening mass scale resulting in the release of
the planets from such traps (Coleman & Nelson, 2016; Speedie et al., 2022).
These earlier studies were not informed by the recent results of Kanagawa et
al. (2018) in which the mass scale at which gaps are opened and slower type II
migration occurs is much lower - in the Neptune mass range for low viscosity
disks. This is about the mass scale at which planets are released from planet
traps. We find therefore that while our higher viscosity planet tracks may
lead to greater incursions of planets to smaller disk radii (perhaps by
factors of 2) that the lower viscosity tracks should remain accurate given
these new results.

It is possible that planets may be able to form in a winds-dominated disk
without rapid migration into the host star, making the trapped type-I
migration regime unnecessary in some systems. A number of studies (Ogihara et
al., 2018; Suzuki et al., 2016; Chambers, 2019) have shown that for the case
of strong disk winds, the inner regions of the disk are essentially evacuated.
Due to this hollowing out of the inner disk, the peak of the disk column
density distribution can occur out at 10 AU. This peak acts essentially as a
planet trap, and it effectively reduces rapid inwards Type I migration from
occurring. As noted in the Introduction, it is far from clear that disk wind
mass loss rates are comparable to accretion rates in general. Were this to be
generally true, one must wonder how star formation would occur in such disks
when at least half of the material moving through the disk would be ejected.
Star formation is inefficient since only about a third of the gas in a gaseous
core actually accretes onto a star (André et al., 2014). Models have shown
that the protostellar outflow can drive off the bulk of the accreting envelope
(Matzner & McKee, 2000).

This brings us to the case of disks which evolve solely by disk wind torques.
In the limit of vanishing turbulence flows within disks are laminar. Dynamical
instabilities such as the formation of vortices by Rossby wave instabilities
are no longer suppressed by viscous damping effects and can therefore play a
role in co-rotation torques. Detailed treatments of co-rotation torque in
wind-evolving disks (McNally et al. 2018, 2020; Kimmig et al. 2020) have
focused on the effects of winds in the immediate co-rotation region near the
planet. Disk winds with high enough mass loss rates can drive fast accretion
flows which result in outward directed co-rotation torques (Kimmig et al.,
2020; Speedie et al., 2022). As one moves from 2D to 3D, simulations reveal
the surprising importance of hydrodynamical buoyancy instabilities that do not
occur in 2D disks. These lead to a strong inward directed torque on the planet
(McNally et al., 2020). Detailed 3D MHD simulations are still needed to
further explore this effect and what consequences it might have for planetary
populations. We showed that planet migration is inwards at all disk radii in
the standard, constrained outflow case using the standard Paardekooper et al.
(2010) type-I torque formula (figure 12), so it is unclear how cores could
avoid encountering traps in our framework. The unconstrained outflow case does
show that a null torque exists and the direction of migration is outwards in
the inner disk. This is counteracted by the rapidly decreasing surface density
of the disk under the action of the assumed very strong wind, offering little
time for core accretion to take place.

###  5.4 Conclusions

This investigation has focused on the relative importance of disk turbulence
and winds as well as the initial distributions of disk properties and
turbulence strengths, in defining planetary populations and their chemical
compositions. All of our thousands of evolution tracks are based on single
planet formation within disks and we have discussed the likely effects of
multi-planet interactions. Our principle conclusions are as follows:

  * •

Our best comparison with the observed M-a distribution arises from populations
that incorporate a range of relative strengths of turbulence and disk winds.
When varying the $\alpha_{\rm{turb}}$ parameter that controls the strength of
turbulence in a population, more scatter is introduced into the final M-a
distribution that significantly improves correspondence with the data compared
to populations that consider constant settings of this disk parameter.

  * •

Most notably, this varied $\alpha_{\rm{turb}}$ population forms a large amount
of super Earths across a large orbital radius extent: 0.01 - 2 AU. This result
of the combined winds and turbulence disk model is therefore an improvement
over populations arising from the pure viscous scenario that we investigated
in previous works, as the combined model directly produces super Earths at
small $a_{p}<0.1$ AU from the disk phase, whereas the pure turbulent models
rarely formed super Earths at these small orbital radii. Additionally, the
super Earth region of the M-a distribution was populated by each of the three
planet traps in our model, which leads to a diverse range of planetary
compositions from dry-rocky super Earths to ice rich planets.

  * •

Using this varied $\alpha_{\rm{turb}}$ model within the combined winds and
turbulent disk framework, we obtain a reasonable comparison to the observed
M-R distribution within the super Earth-Neptune mass regime. Solid-dominated
planets with masses $\lesssim$ 2-3 M⊕ compare well with the data. Planets at
larger masses either retain their accreted atmospheres and correspond with
observed planets with largest radii at a given mass, or are stripped via
photoevaporation, reducing their radii to compare with average observed planet
radii at the corresponding mass.

  * •

Post-disk atmospheric photoevaporation has a critical role in reducing super
Earth radii and improving our comparison with the observed M-R distribution.
Compared to the population arising from the pure turbulent disk of Alessi et
al. (2020b), a larger number of super Earths form at smaller orbital radii in
this model, and therefore more planets are affected by photoevaporation. This
increases the number of stripped planets at masses greater than 3 M⊕, thereby
improving our comparison with the M-R data compared to that of the pure
turbulent disk.

  * •

The varied $\alpha_{\rm{turb}}$ population also produces a separation between
warm gas giants near 1 AU formed in the ice line and hot Jupiters formed in
the dead zone. 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 ice line forms many warm gas giants, as well as super Earths and Neptunes
mostly in the 0.5-2 AU range. The heat transition 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.

  * •

When considering populations with constant settings of the disk turbulence
parameter, a setting of $\alpha_{\rm{turb}}=10^{-3}$ resulted in a population
with the largest range of final planet properties and scatter on the M-a
diagram compared to lower (10-4) and higher (0.007) settings of this
parameter, that resulted in populations that covered very limited regions of
the M-a space. Planet formation in disks with $\alpha_{\rm{turb}}\simeq
10^{-3}$ is therefore most sensitive to disk parameters compared to lower and
higher settings. However, all of the constant $\alpha_{\rm{turb}}$ models
failed to produce populations with sufficient scatter in their M-a
distribution to compare reasonably with the data, and populations the varied
the turbulent strength were necessary to improve this comparison.

  * •

Winds-dominated models did not show variation in planet formation results.
They fail to reproduce the observed variation in the M-a distribution.

Our results clearly indicate that disk evolution and planet formation are
sensitive to the strengths of these two disk evolution mechanisms. They
strongly reinforce the idea that disk winds provide the dominant means of
angular momentum loss for most of the initial disk populations. Turbulence is
still expected to be present and important and is likely produced by weaker
dynamical instabilities that are left after MRI is suppressed. A distribution
of turbulence strengths naturally creates a greater variety of outcomes for
planet populations. We emphasize that planet traps continue to play a vital
role in defining planet populations. Disk winds should in general be highly
efficient in carrying off disk angular momentum so that the hollowing out of
the inner disk is not necessary to produce a population of warm Jupiters
trapped at a resulting column density peak out at several AU. The major caveat
in our calculations is that we overproduce hot Jupiters at the expense of
SuperEarths. This can be remedied by including the stripping of lower mass
planets migrating with the inward moving dead zone as it crosses the ice line
at a few AU. These results motivate further phyiscal modeling of disk
evolution processes, particularly in 3D MHD simulation of the co-rotation
region very near to accreting and migrating forming planets. This is a major
topic of current research to which we intend to make further contributions.

## Acknowledgements

We thank the anonymous referee for an insightful and very useful report that
improved the manuscript. We also thank Jessica Speedie, Alex Cridland, Nienke
van der Marel, and James Wadsley for many stimulating discussions. M.A.
acknowledges funding from the National Sciences and engineering Research
Council (NSERC) through the Alexander Graham Bell CGS/PGS Doctoral Scholarship
and from an Ontario graduate scholarship. R.E.P. is supported by an NSERC
Discovery Grant. This work made use of Compute/Caclul Canada.

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * ALMA Partnership et al. (2015) ALMA Partnership et al., 2015, ApJ, 808, L3
  * Alessi & Pudritz (2018) Alessi M., Pudritz R. E., 2018, MNRAS, 478, 2599
  * Alessi et al. (2017) Alessi M., Pudritz R. E., Cridland A. J., 2017, MNRAS, 464, 428
  * Alessi et al. (2020a) Alessi M., Pudritz R. E., Cridland A. J., 2020a, MNRAS, 493, 1013
  * Alessi et al. (2020b) Alessi M., Inglis J., Pudritz R. E., 2020b, MNRAS, 497, 4814
  * Alibert et al. (2011) Alibert Y., Mordasini C., Benz W., 2011, A&A, 526, A63
  * André et al. (2014) André P., Di Francesco J., Ward-Thompson D., Inutsuka S. I., Pudritz R. E., Pineda J. E., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 27 (arXiv:1312.6232), doi:10.2458/azu_uapress_9780816531240-ch002
  * Andrews (2020) Andrews S. M., 2020, ARA&A, 58, 483
  * Andrews et al. (2018) Andrews S. M., et al., 2018, ApJ, 869, L41
  * Ansdell et al. (2016) Ansdell M., et al., 2016, ApJ, 828, 46
  * Ansdell et al. (2018) Ansdell M., et al., 2018, ApJ, 859, 21
  * Armitage (2011) Armitage P. J., 2011, ARA&A, 49, 195
  * Armitage (2018) Armitage P. J., 2018, in Deeg H. J., Belmonte J. A., eds, , Handbook of Exoplanets. p. 135, doi:10.1007/978-3-319-55333-7_135
  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  * Bai (2015) Bai X.-N., 2015, ApJ, 798, 84
  * Bai (2016) Bai X.-N., 2016, ApJ, 821, 80
  * Bai & Stone (2011) Bai X.-N., Stone J. M., 2011, ApJ, 736, 144
  * Bai & Stone (2013) Bai X.-N., Stone J. M., 2013, ApJ, 769, 76
  * Bai et al. (2016) Bai X.-N., Ye J., Goodman J., Yuan F., 2016, ApJ, 818, 152
  * Baillié et al. (2016) Baillié K., Charnoz S., Pantin E., 2016, A&A, 590, A60
  * Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
  * Batalha et al. (2013) Batalha N. M., et al., 2013, ApJS, 204, 24
  * Bell & Lin (1994) Bell K. R., Lin D. N. C., 1994, ApJ, 427, 987
  * Benz et al. (2014) Benz W., Ida S., Alibert Y., Lin D., Mordasini C., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI, p. 691
  * Birnstiel et al. (2012) Birnstiel T., Klahr H., Ercolano B., 2012, A&A, 539, A148
  * Bitsch et al. (2015) Bitsch B., Lambrechts M., Johansen A., 2015, A&A, 582, A112
  * Bitsch et al. (2019) Bitsch B., Izidoro A., Johansen A., Raymond S. N., Morbidelli A., Lambrechts M., Jacobson S. A., 2019, A&A, 623, A88
  * Bjerkeli et al. (2016) Bjerkeli P., van der Wiel M. H. D., Harsono D., Ramsey J. P., Jørgensen J. K., 2016, Nature, 540, 406
  * Blandford & Payne (1982) Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883
  * Bryan et al. (2016) Bryan M. L., et al., 2016, ApJ, 821, 89
  * Chambers (2009) Chambers J. E., 2009, ApJ, 705, 1206
  * Chambers (2019) Chambers J., 2019, ApJ, 879, 98
  * Chatterjee et al. (2008) Chatterjee S., Ford E. B., Matsumura S., Rasio F. A., 2008, ApJ, 686, 580
  * Chen et al. (2016) Chen X., Arce H. G., Zhang Q., Launhardt R., Henning T., 2016, ApJ, 824, 72
  * Chiang & Goldreich (1997) Chiang E. I., Goldreich P., 1997, ApJ, 490, 368
  * Cleeves et al. (2013) Cleeves L. I., Adams F. C., Bergin E. A., 2013, ApJ, 772, 5
  * Cleeves et al. (2015) Cleeves L. I., Bergin E. A., Qi C., Adams F. C., Öberg K. I., 2015, ApJ, 799, 204
  * Coleman & Nelson (2014) Coleman G. A. L., Nelson R. P., 2014, MNRAS, 445, 479
  * Coleman & Nelson (2016) Coleman G. A. L., Nelson R. P., 2016, MNRAS, 457, 2480
  * Coleman et al. (2017) Coleman G. A. L., Papaloizou J. C. B., Nelson R. P., 2017, MNRAS, 470, 3206
  * Cridland et al. (2016) Cridland A. J., Pudritz R. E., Alessi M., 2016, MNRAS, 461, 3274
  * Cridland et al. (2019a) Cridland A. J., Eistrup C., van Dishoeck E. F., 2019a, A&A, 627, A127
  * Cridland et al. (2019b) Cridland A. J., van Dishoeck E. F., Alessi M., Pudritz R. E., 2019b, A&A, 632, A63
  * Cumming et al. (2008) Cumming A., Butler R. P., Marcy G. W., Vogt S. S., Wright J. T., Fischer D. A., 2008, PASP, 120, 531
  * D’Alessio et al. (1999) D’Alessio P., Calvet N., Hartmann L., Lizano S., Cantó J., 1999, ApJ, 527, 893
  * Dittkrist et al. (2014) Dittkrist K.-M., Mordasini C., Klahr H., Alibert Y., Henning T., 2014, A&A, 567, A121
  * Dong et al. (2015) Dong R., Zhu Z., Whitney B., 2015, ApJ, 809, 93
  * Dullemond et al. (2018) Dullemond C. P., et al., 2018, ApJ, 869, L46
  * Emsenhuber et al. (2021) Emsenhuber A., Mordasini C., Burn R., Alibert Y., Benz W., Asphaug E., 2021, A&A, 656, A70
  * Ferreira & Pelletier (1995) Ferreira J., Pelletier G., 1995, A&A, 295, 807
  * Fischer et al. (2016) Fischer D. A., et al., 2016, PASP, 128, 066001
  * Flaherty et al. (2018) Flaherty K. M., Hughes A. M., Teague R., Simon J. B., Andrews S. M., Wilner D. J., 2018, ApJ, 856, 117
  * Ford & Rasio (2008) Ford E. B., Rasio F. A., 2008, ApJ, 686, 621
  * Fung & Chiang (2016) Fung J., Chiang E., 2016, ApJ, 832, 105
  * Goldreich & Tremaine (1980) Goldreich P., Tremaine S., 1980, ApJ, 241, 425
  * Gressel et al. (2015) Gressel O., Turner N. J., Nelson R. P., McNally C. P., 2015, ApJ, 801, 84
  * Gressel et al. (2020) Gressel O., Ramsey J. P., Brinch C., Nelson R. P., Turner N. J., Bruderer S., 2020, ApJ, 896, 126
  * Gupta & Schlichting (2019) Gupta A., Schlichting H. E., 2019, MNRAS, 487, 24
  * Haffert et al. (2019) Haffert S. Y., Bohn A. J., de Boer J., Snellen I. A. G., Brinchmann J., Girard J. H., Keller C. U., Bacon R., 2019, Nature Astronomy, 3, 749
  * Hartmann et al. (1998) Hartmann L., Calvet N., Gullbring E., D’Alessio P., 1998, ApJ, 495, 385
  * Hasegawa (2016) Hasegawa Y., 2016, ApJ, 832, 83
  * Hasegawa & Pudritz (2010) Hasegawa Y., Pudritz R. E., 2010, ApJ, 710, L167
  * Hasegawa & Pudritz (2011) Hasegawa Y., Pudritz R. E., 2011, MNRAS, 417, 1236
  * Hasegawa & Pudritz (2012) Hasegawa Y., Pudritz R. E., 2012, ApJ, 760, 117
  * Hasegawa & Pudritz (2014) Hasegawa Y., Pudritz R. E., 2014, ApJ, 794, 25
  * Hasegawa et al. (2017) Hasegawa Y., Okuzumi S., Flock M., Turner N. J., 2017, ApJ, 845, 31
  * Hellary & Nelson (2012) Hellary P., Nelson R. P., 2012, MNRAS, 419, 2737
  * Hernández et al. (2007) Hernández J., et al., 2007, ApJ, 662, 1067
  * Ida & Lin (2004a) Ida S., Lin D. N. C., 2004a, ApJ, 604, 388
  * Ida & Lin (2004b) Ida S., Lin D. N. C., 2004b, ApJ, 616, 567
  * Ida & Lin (2008) Ida S., Lin D. N. C., 2008, ApJ, 673, 487
  * Ida et al. (2018) Ida S., Tanaka H., Johansen A., Kanagawa K. D., Tanigawa T., 2018, ApJ, 864, 77
  * Ikoma et al. (2000) Ikoma M., Nakazawa K., Emori H., 2000, ApJ, 537, 1013
  * Jackson et al. (2012) Jackson A. P., Davis T. A., Wheatley P. J., 2012, Monthly Notices of the Royal Astronomical Society, 422, 2024
  * Jin et al. (2014) Jin S., Mordasini C., Parmentier V., van Boekel R., Henning T., Ji J., 2014, ApJ, 795, 65
  * Johansen et al. (2007) Johansen A., Oishi J. S., Mac Low M.-M., Klahr H., Henning T., Youdin A., 2007, Nature, 448, 1022
  * Kanagawa et al. (2018) Kanagawa K. D., Tanaka H., Szuszkiewicz E., 2018, ApJ, 861, 140
  * Keppler et al. (2018) Keppler M., et al., 2018, A&A, 617, A44
  * Kimmig et al. (2020) Kimmig C. N., Dullemond C. P., Kley W., 2020, A&A, 633, A4
  * Klahr et al. (2018) Klahr H., Pfeil T., Schreiber A., 2018, in Deeg H. J., Belmonte J. A., eds, , Handbook of Exoplanets. p. 138, doi:10.1007/978-3-319-55333-7_138
  * Kley & Nelson (2012) Kley W., Nelson R. P., 2012, ARA&A, 50, 211
  * Kokubo & Ida (2002) Kokubo E., Ida S., 2002, ApJ, 581, 666
  * Lesur et al. (2014) Lesur G., Kunz M. W., Fromang S., 2014, A&A, 566, A56
  * Lesur et al. (2022) Lesur G., et al., 2022, arXiv e-prints, p. arXiv:2203.09821
  * Lin & Papaloizou (1986) Lin D. N. C., Papaloizou J., 1986, ApJ, 309, 846
  * Long et al. (2022) Long F., et al., 2022, ApJ, 931, 6
  * Lopez & Fortney (2013) Lopez E. D., Fortney J. J., 2013, ApJ, 776, 2
  * Louvet et al. (2018) Louvet F., Dougados C., Cabrit S., Mardones D., Ménard F., Tabone B., Pinte C., Dent W. R. F., 2018, A&A, 618, A120
  * Lyra et al. (2010) Lyra W., Paardekooper S.-J., Mac Low M.-M., 2010, ApJ, 715, L68
  * Madhusudhan (2018) Madhusudhan N., 2018, in Deeg H. J., Belmonte J. A., eds, , Handbook of Exoplanets. p. 104, doi:10.1007/978-3-319-55333-7_104
  * Manara et al. (2019) Manara C. F., Mordasini C., Testi L., Williams J. P., Miotello A., Lodato G., Emsenhuber A., 2019, A&A, 631, L2
  * Masset et al. (2006) Masset F. S., Morbidelli A., Crida A., Ferreira J., 2006, ApJ, 642, 478
  * Matsumura & Pudritz (2003) Matsumura S., Pudritz R. E., 2003, ApJ, 598, 645
  * Matsumura et al. (2021) Matsumura S., Brasser R., Ida S., 2021, A&A, 650, A116
  * Matzner & McKee (2000) Matzner C. D., McKee C. F., 2000, ApJ, 545, 364
  * McNally et al. (2017) McNally C. P., Nelson R. P., Paardekooper S.-J., Gressel O., Lyra W., 2017, MNRAS, 472, 1565
  * McNally et al. (2018) McNally C. P., Nelson R. P., Paardekooper S.-J., 2018, MNRAS, 477, 4596
  * McNally et al. (2020) McNally C. P., Nelson R. P., Paardekooper S.-J., Benítez-Llambay P., Gressel O., 2020, MNRAS, 493, 4382
  * Mordasini et al. (2009) Mordasini C., Alibert Y., Benz W., 2009, A&A, 501, 1139
  * Mordasini et al. (2012) Mordasini C., Alibert Y., Klahr H., Henning T., 2012, A&A, 547, A111
  * Mordasini et al. (2014) Mordasini C., Klahr H., Alibert Y., Miller N., Henning T., 2014, A&A, 566, A141
  * Mordasini et al. (2015) Mordasini C., Mollière P., Dittkrist K.-M., Jin S., Alibert Y., 2015, International Journal of Astrobiology, 14, 201
  * Mori et al. (2019) Mori S., Bai X.-N., Okuzumi S., 2019, ApJ, 872, 98
  * Morton et al. (2016) Morton T. D., Bryson S. T., Coughlin J. L., Rowe J. F., Ravichandran G., Petigura E. A., Haas M. R., Batalha N. M., 2016, ApJ, 822, 86
  * Mulders et al. (2017) Mulders G. D., Pascucci I., Manara C. F., Testi L., Herczeg G. J., Henning T., Mohanty S., Lodato G., 2017, ApJ, 847, 31
  * Müller et al. (2018) Müller A., et al., 2018, A&A, 617, L2
  * Murray-Clay et al. (2009) Murray-Clay R. A., Chiang E. I., Murray N., 2009, The Astrophysical Journal, 693, 23–42
  * Nagasawa et al. (2008) Nagasawa M., Ida S., Bessho T., 2008, ApJ, 678, 498
  * Nelson (2018) Nelson R. P., 2018, in Deeg H. J., Belmonte J. A., eds, , Handbook of Exoplanets. p. 139, doi:10.1007/978-3-319-55333-7_139
  * Nisini et al. (2018) Nisini B., Antoniucci S., Alcalá J. M., Giannini T., Manara C. F., Natta A., Fedele D., Biazzo K., 2018, A&A, 609, A87
  * Öberg et al. (2011) Öberg K. I., Boogert A. C. A., Pontoppidan K. M., van den Broek S., van Dishoeck E. F., Bottinelli S., Blake G. A., Evans II N. J., 2011, ApJ, 740, 109
  * Ogihara et al. (2018) Ogihara M., Kokubo E., Suzuki T. K., Morbidelli A., 2018, A&A, 615, A63
  * Ouyed & Pudritz (1997) Ouyed R., Pudritz R. E., 1997, ApJ, 482, 712
  * Owen & Jackson (2012) Owen J. E., Jackson A. P., 2012, Monthly Notices of the Royal Astronomical Society, 425, 2931
  * Owen & Wu (2013) Owen J. E., Wu Y., 2013, ApJ, 775, 105
  * Paardekooper et al. (2010) Paardekooper S. J., Baruteau C., Crida A., Kley W., 2010, MNRAS, 401, 1950
  * Paardekooper et al. (2011) Paardekooper S.-J., Baruteau C., Kley W., 2011, MNRAS, 410, 293
  * Paardekooper et al. (2022) Paardekooper S.-J., Dong R., Duffell P., Fung J., Masset F. S., Ogilvie G., Tanaka H., 2022, arXiv e-prints, p. arXiv:2203.09595
  * Pascucci et al. (2016) Pascucci I., et al., 2016, ApJ, 831, 125
  * Pascucci et al. (2022) Pascucci I., Cabrit S., Edwards S., Gorti U., Gressel O., Suzuki T., 2022, arXiv e-prints, p. arXiv:2203.10068
  * Pelletier & Pudritz (1992) Pelletier G., Pudritz R. E., 1992, ApJ, 394, 117
  * Petigura et al. (2018) Petigura E. A., et al., 2018, AJ, 155, 89
  * Pinilla et al. (2020) Pinilla P., Pascucci I., Marino S., 2020, A&A, 635, A105
  * Pinte et al. (2016) Pinte C., Dent W. R. F., Ménard F., Hales A., Hill T., Cortes P., de Gregorio-Monsalvo I., 2016, ApJ, 816, 25
  * Pinte et al. (2022) Pinte C., Teague R., Flaherty K., Hall C., Facchini S., Casassus S., 2022, arXiv e-prints, p. arXiv:2203.09528
  * Pudritz & Norman (1986) Pudritz R. E., Norman C. A., 1986, ApJ, 301, 571
  * Pudritz & Ray (2019) Pudritz R. E., Ray T. P., 2019, Frontiers in Astronomy and Space Sciences, 6, 54
  * Pudritz et al. (2018) Pudritz R. E., Cridland A. J., Alessi M., 2018, in Deeg H. J., Belmonte J. A., eds, , Handbook of Exoplanets. p. 144, doi:10.1007/978-3-319-55333-7_144
  * Rafikov (2002) Rafikov R. R., 2002, ApJ, 572, 566
  * Reiners et al. (2018) Reiners A., et al., 2018, A&A, 612, A49
  * Ribas et al. (2005) Ribas I., Guinan E. F., Güdel M., Audard M., 2005, ApJ, 622, 680
  * Rodenkirch et al. (2020) Rodenkirch P. J., Klahr H., Fendt C., Dullemond C. P., 2020, A&A, 633, A21
  * Rosotti et al. (2020) Rosotti G. P., Teague R., Dullemond C., Booth R. A., Clarke C. J., 2020, MNRAS, 495, 173
  * Rowe et al. (2014) Rowe J. F., et al., 2014, ApJ, 784, 45
  * Ruden & Pollack (1991) Ruden S. P., Pollack J. B., 1991, ApJ, 375, 740
  * Sándor et al. (2011) Sándor Z., Lyra W., Dullemond C. P., 2011, ApJ, 728, L9
  * Santerne et al. (2016) Santerne A., et al., 2016, A&A, 587, A64
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Simon et al. (2013) Simon J. B., Bai X.-N., Stone J. M., Armitage P. J., Beckwith K., 2013, ApJ, 764, 66
  * Speedie et al. (2022) Speedie J., Pudritz R. E., Cridland A. J., Meru F., Booth R. A., 2022, MNRAS, 510, 6059
  * Spruit (2010) Spruit H. C., 2010, in Belloni T., ed., , Vol. 794, Lecture Notes in Physics, Berlin Springer Verlag. p. 233, doi:10.1007/978-3-540-76937-8_9
  * Stepinski (1998) Stepinski T. F., 1998, Icarus, 132, 100
  * Suzuki et al. (2016) Suzuki T. K., Ogihara M., Morbidelli A. r., Crida A., Guillot T., 2016, A&A, 596, A74
  * Tabone et al. (2020) Tabone B., et al., 2020, A&A, 640, A82
  * Tazzari et al. (2017) Tazzari M., et al., 2017, A&A, 606, A88
  * Thompson et al. (2018) Thompson S. E., et al., 2018, ApJS, 235, 38
  * Trapman et al. (2019) Trapman L., Facchini S., Hogerheijde M. R., van Dishoeck E. F., Bruderer S., 2019, arXiv e-prints, p. arXiv:1903.06190
  * Valencia et al. (2007) Valencia D., Sasselov D. D., O’Connell R. J., 2007, ApJ, 665, 1413
  * Watson et al. (2016) Watson D. M., et al., 2016, ApJ, 828, 52
  * Winn & Fabrycky (2015) Winn J. N., Fabrycky D. C., 2015, ARA&A, 53, 409
  * Youdin & Goodman (2005) Youdin A. N., Goodman J., 2005, ApJ, 620, 459
  * Zhang et al. (2018) Zhang Y., et al., 2018, ApJ, 864, 76
  * Zhu & Stone (2018) Zhu Z., Stone J. M., 2018, ApJ, 857, 34
  * de Valon et al. (2020) de Valon A., Dougados C., Cabrit S., Louvet F., Zapata L. A., Mardones D., 2020, A&A, 634, L12
  * van Boekel et al. (2017) van Boekel R., et al., 2017, ApJ, 837, 132
  * van der Marel & Mulders (2021) van der Marel N., Mulders G., 2021, arXiv e-prints, p. arXiv:2104.06838

##  Appendix A Wind-dominated models & the effect of outflow strength

We now consider the case of a winds-dominated disk model with
$\alpha_{\rm{turb}}=10^{-6}$, which is same value used within the combined
model’s dead zone. We will examine the effect of the strength of the wind-
driven outflow by comparing two models: the constrained case where $K$ is
solved for using equation 14, and the unconstrained case where we use a high
setting of $K=1$. In this disk model we do not consider a dead zone as the
turbulence strength is quite low throughout the entire disk’s extent. The
related trap at the Ohmic dead zone’s outer edge is therefore not present in
these models.

Figure 11: We compare disk evolution in winds-dominated disks
($\alpha_{\rm{turb}}$ = 10-6) at two settings of disk outflow strength. The
left column pertains to the ‘constrained’ outflow where $K$ is determined
according to equation 14. The right column considers a high outflow strength
with $K$=1. The top row shows evolution of disk accretion rate, and middle and
bottom rows show radial profiles of surface density and midplane temperature,
respectively.

In figure 11, we plot the resulting disk evolution corresponding to each
scenario. We find that, even though the strength of turbulence has been
reduced by a factor of 100 compared to the combined model of the previous
section, the “constrained outflow” model’s accretion rate and surface density
evolution is quite comparable to the $\alpha_{\rm{turb}}=10^{-4}$ case. The
accretion rate decreases from its initial value of 6$\times 10^{-8}$ M⊙ yr-1
to roughly $10^{-9}$ M⊙ yr-1 after 3 Myr of disk evolution. The surface
density profiles also show the disk to contract with time in terms of its
outer radius; a result of winds-dominated evolution.

Figure 12: Profiles of the normalized type-I migration torque
$\Gamma_{\rm{I}}/\Gamma_{0}$ are shown for both the constrained and
unconstrained outflow models, calculated using equation 19. Inward type-I
migration will persist across the entire ‘constrained outflow’ disk’s extent.
Conversely, we find a planet trap near 10 AU related to the surface density
maximum in the ‘unconstrained outflow’ model.

The constrained outflow model’s temperature profile is, however, significantly
different from the $\alpha_{\rm{turb}}=10^{-4}$ case, as viscous heating has
been substantially decreased due to the reduced strength of turbulence. This
decrease results in the disk’s viscous region being limited to the innermost
region within 1 AU and early times in disk evolution $\lesssim$ 1 Myr. Between
1 and 2 Myr, viscous heating has decreased to the point where the entire
disk’s heating is dominated by host-stellar radiation. Once this has occurred,
the midplane temperature profiles simply correspond to the radiative
equilibrium profile (equation 11). As a consequence, the midplane temperature
profile will become completely static, since the radiative temperature profile
has no time-dependence. This will result in the ice line’s position being
stationary in the disk throughout its evolution. Furthermore, if disk
chemistry were to be computed in this disk, the profiles would not evolve with
time after the viscously-heated region has disappeared. We note that these
implications of a static, radiative-equilibrium temperature profile are a
result of the common assumption that we have used throughout this thesis that
the host-star’s luminosity is constant.

Now considering the unconstrained outflow scenario in the right column of
figure 11, we can immediately see the interesting effect of a strong, winds-
driven outflow on disk evolution. In this model, a larger fraction of the disk
winds’ stress is carried away in the outflow as opposed to contributing to
disk accretion. This results in the disk accretion rate $\dot{M}$ being
reduced by roughly two orders of magnitude compared to the constrained outflow
scenario.

Another compelling effect of a strong outflow can be readily seen in the disk
surface density profiles. Rather than being a decreasing function of disk
radius over the entire disk’s extent, we see that $\Sigma$ increases until
roughly 10 AU where it achieves a maximum value before decreasing in the outer
disk. This is a result of of the outflow surface density rate being a
decreasing function of radius $r$ (see equation 5), removing material more
efficiently in the inner disk than the outer disk. This feature is also
present in many of the surface profiles shown in the numerical treatment of
Suzuki et al. (2016).

Since, in addition to having a low setting of turbulence, the unconstrained
outflow model’s accretion is substantially lower, there is no viscously-heated
region in the disk. At all times, the midplane temperature is a result of
heating through radiation alone, and the midplane temperature profiles are
equal to $T_{\rm{req}}$ (equation 11) for the entirety of disk evolution.
There are two consequences of this result. First, the heat transition and its
related trap are not present in this model, as there is no transition into a
viscously-heated region present in the disk. Second, as we have described
above, the midplane temperature profile is completely static, which will
result in an ice line radius that does not change as the disk evolves.

To summarize, since we are considering an extremely low setting of
$\alpha_{\rm{turb}}$, there is no longer a dead zone or its related trap in
either winds-dominated model. Additionally, in the unconstrained outflow case,
there will be no heat transition, with only a static ice line remaining from
our standard set of three planet traps. Both the heat transition and ice line
will remain in the constrained outflow scenario.

Figure 13: We plot evolution of planet traps (top row) and planet formation
tracks (bottom row) for both the constrained outflow (left column) and
unconstrained outflow (right column) models. In the planet formation plots,
open circles designate the planets’ positions at 1 Myr intervals throughout
their disks’ 10 Myr-lifetimes.

However, the surface density maximum that is encountered in the unconstrained
outflow model has significant implications for planet migration. Type-I
migration models have shown that forming cores will experience zero net torque
at locations of surface density maxima (i.e. Sándor et al. (2011)). In fact,
this is our main motivation for analyzing the unconstrained outflow model,
despite the fact that it produces large outflow mass-loss rates that are in
contention with disk winds theory and observations (i.e. equation 14).

We adopt the standard type-I migration torque formula of Paardekooper et al.
(2010) to determine if a null torque location (i.e. a planet trap) will exist
in either model;

$\gamma\frac{\Gamma_{\rm{I}}}{\Gamma_{0}}=-2.5-1.7\beta+0.1\alpha+1.1(1.5-\alpha)+7.9\frac{\xi}{\gamma}\;,$
(19)

where $\Gamma_{\rm{I}}$ is the normalized type-I migration torque whose sign
indicates the direction of planet migration ($\Gamma_{\rm{I}}<0$ indicating
inward migration). $\alpha$ and $\beta$ correspond to the magnitudes of the
local power-law indices of the surface density, and temperature profiles,
respectively (i.e. $\Sigma\sim r^{-\alpha}$, $T\sim r^{-\beta}$). This torque
formula illustrates that the type-I migration torque is sensitive to the local
disk surface density and temperature profiles. We use an adiabatic index of
$\gamma=5/3$ corresponding to a monatomic gas. $\xi$ is the disk’s entropy
gradient; $\xi=\beta-(\gamma-1)\alpha$. Lastly, the torque normalization
factor $\Gamma_{0}$ is,

$\Gamma_{0}=\left(\frac{q}{h}\right)^{2}\Sigma_{p}r_{p}^{4}\Omega_{p}^{2}\;,$
(20)

where $q$ is the planet to star mass ratio, $h$ is the disk aspect ratio, and
all quantities are computed at the planet’s location.

We will, however, only calculate the _normalized_ type-I migration torque
$\Gamma_{\rm{I}}/\Gamma_{0}$ using equation 19, to determine if any planet
traps exist in either the constrained or unconstrained outflow models. These
will correspond to locations where $\Gamma_{\rm{I}}$ = 0. We note that, since
we are calculating the _normalized_ torque, we do not need to prescribe a core
mass (which would be needed to determine the torque normalization, equation
20).

The result of this calculation is shown in figure 12, where radial profiles of
the normalized torque are shown throughout both disks’ evolutions. We see that
in the constrained outflow model, the type-I migration direction will be
inward across the entire disk’s extent, and for the entire disk’s evolution.
The constrained model’s surface density profiles are similar to the ‘combined’
turbulence and winds model of the previous section, so we expect this result
to apply to that disk model as well.

In the case of the unconstrained outflow model, we do indeed find a planet
trap near 6 AU where $\Gamma_{\rm{I}}/\Gamma_{0}$ = 0, related to the maximum
in the disk surface density profile. Interior to the trap, the direction of
the type-I torque is outward, while exterior to the trap it is inward. The
directions of the type-I migration torque are therefore appropriated to
migrate planets towards the trap. We also find that the torque profiles in the
unconstrained outflow scenario are time-independent. In addition to the static
temperature profile, the radial profile of the surface density’s power-law
index $\alpha$ is also time-independent, which causes the static torque
profiles. This can be seen in figure 11, where the unconstrained outflow
model’s $\Sigma$ profile decreases while maintaining its radial structure.

Based on this result of figure 12, we obtain a second planet trap in the
unconstrained outflow case, in addition to the ice line. We refer to this trap
simply as the “null torque” for this model. However, since the torque profiles
and null torque radius are static, this trap will also be stationary in the
disk throughout its evolution.

In figure 13, we plot the evolution of the traps’ radii and planet formation
tracks in both the constrained and unconstrained outflow models. Here, we only
consider the fiducial setting of $\Sigma_{0}$. In the constrained outflow
model, we see that the heat transition is located at a much smaller radius in
the disk than either the turbulent or combined models of the previous section.
This is as a result of the lower $\alpha_{\rm{turb}}=10^{-6}$ and reduced
viscous heating. As previously outlined, the ice line’s position is stationary
in the disk near 1 AU, due to the static radiative equilibrium temperature
profile. Since the traps exist at small disk radii, and the constrained
outflow model maintains large surface densities within 10-20 AU, both traps
efficiently produce warm Jupiters near 1 AU within only 1 Myr of formation.

In the unconstrained outflow model, the radial locations of both the null
torque and ice line traps are time-independent. In contrast to the constrained
outflow model, planet formation in each of the two traps is extremely
inefficient due to the rapid decrease in disk surface density caused by the
strong outflow with $K=1$. Even considering a long disk lifetime of 10 Myr
does not result in massive planets forming, since the disk surface density has
decreased significantly by 1 Myr due to the outflow (see figure 11). The null
torque only forms a sub-M⊕ core in this disk model, while the ice line forms a
$\sim$ 3 M⊕ super Earth. In both cases, the planets form in-situ due to the
time-independence of their traps’ locations.

We have also separately investigated how the positions of the traps change in
both models with different settings of $\Sigma_{0}$. We find that the trap
locations are extremely/entirely insensitive to the setting of $\Sigma_{0}$.
This is mainly the case because the radiative equilibrium temperature profile,
which dominates midplane heating, is independent of disk surface density (only
depending on the radiative flux received by the disk). The resulting planet
formation tracks are also very similar to those shown in figure 13. Based on
these results that show limited planet formation and no variance with disk
parameters, we do not expect a population synthesis model to be particularly
interesting in either the constrained or unconstrained outflow models in terms
of the planet traps framework. However, the low settings of
$\alpha_{\rm{turb}}$ in a winds-dominated model may allow planets to form
under a general type-I migration regime, as calculated with equation 19, while
not being restricted to the locations of traps. We have discussed this idea
further in section 5.


# Reconstruction of manifold embeddings into Euclidean spaces via intrinsic
distances

Nikita Puchkin  HSE University and Institute for Information Transmission
Problems, Moscow<EMAIL_ADDRESS>,  Vladimir Spokoiny  WIAS and Humboldt
University, Berlin, and Institute for Information Transmission Problems and
Skoltech, Moscow<EMAIL_ADDRESS>,  Eugene Stepanov  St.Petersburg
Branch of the Steklov Mathematical Institute of the Russian Academy of
Sciences, Fontanka 27, 191023 St.Petersburg, Russia and ITMO University, St.
Petersburg, and HSE University, Moscow<EMAIL_ADDRESS>and  Dario
Trevisan  Dario Trevisan, Dipartimento di Matematica, Università di Pisa
Largo Bruno Pontecorvo 5
I-56127, Pisa<EMAIL_ADDRESS>
(Date: August 27, 2024)

###### Abstract.

We consider one of the classical manifold learning problems, that of
reconstructing up to an almost isometry an embedding of a compact connected
Riemannian manifold in a Euclidean space given the information on intrinsic
distances between points from its almost dense subset. It will be shown that
the most popular methods in data science to deal with such a problem, the
classical Multidimensional scaling (MDS) and the Maximum variance unfolding
(MVU) actually miss the point and may provide results very far from an
isometry (and even may give no biLipshitz embedding). We will then provide an
easy variational formulation of this problem which leads to an algorithm
always providing an almost isometric imbedding with given controlled small
distortion of original distances.

The article was prepared within the framework of the HSE University Basic
Research Program. The second author was supported by the German Ministry for
Education via the Berlin Center for Machine Learning (01IS18037I). The work of
the third author has been also partially financed by the RFBR grant
#20-01-00630A. The fourth author was partially supported by GNAMPA-INdAM 2020
project “Problemi di ottimizzazione con vincoli via trasporto ottimo e
incertezza” and University of Pisa, Project PRA 2018-49.

##  1\. Introduction

Let $M\subset\mathbb{R}^{n}$ be a smooth connected compact Riemannian manifold
which we will always think about as embedded in some Euclidean space
$\mathbb{R}^{n}$. By $d_{M}$ we denote the intrinsic (geodesic) distance on
$M$. Let $\\{y_{i}\\}\subset M$ and denote for the sake of brevity

$d_{ij}:=d_{M}(y_{i},y_{j}).$

In applications one often needs to reconstruct the embedding of the original
manifold in $\mathbb{R}^{n}$ (or at least discover some information about the
latter) knowing only the information about intrinsic distances. Namely, the
problem is as follows: the manifold $M$ is unknown, unlike the matrix
$\\{d_{ij}\\}$. Can one reconstruct the locations of all the points
$y_{j}\in\mathbb{R}^{n}$ (which, if the set $\\{y_{i}\\}\subset M$ is dense in
$M$, will give us the manifold $M$ itself)? Of course, such a problem cannot
be solved uniquely, since if $\tilde{M}\subset\mathbb{R}^{n}$ is an arbitrary
manifold isometric to $M$ via an isometry $f\colon M\to\tilde{M}$, then the
set $\\{f(y_{i})\\}\subset\tilde{M}$ will produce exactly the same intrinsic
distances.

In data science applications the set of points $\\{y_{i}\\}\subset M$ is of
course finite, consisting of quite a large (sometimes even huge) number of
points covering $M$ “almost densely”. In other words we might think that in
fact one deals with a finite subcollection $\\{y_{i}\\}_{i=1}^{N}$ of a dense
set $\\{y_{i}\\}\subset M$, with $N\in\mathbb{N}$ very large. An algorithm for
finding manifold embedding must produce then a set
$\\{x_{i}\\}_{i=1}^{N}\subset\Sigma$, with $\Sigma\subset\mathbb{R}^{n}$ some
other embedded manifold equipped with intrinsic distance $d_{\Sigma}$, so that
$d_{\Sigma}(x_{i},x_{j})\approx d_{ij}$ (ideally, $\approx$ standing for the
approximate equality with a small but quantifiable error).

Note that we are interested in reconstructing the embedding of the original
manifold $M$ into an Euclidean space (e.g. for the purpose of visualizing the
data on the latter), as opposed to the problem of reconstructing an abstract
manifold $M$ (e.g. determined by its metric tensor). The latter is solved in
[6], but its solution does not provide any explicit finite dimensonal
embedding. However, one might think that modern machine learning is abundant
of methods developed precisely for reconstructing embeddings of data
manifolds; they are mostly modifications of either multidimensional scaling
(MDS) or maximum variance unfolding (MVU). Some of them are also extensively
used as a background for dimension reduction methods in data science (for
instance, MDS is in fact used as a basic part of Isomap dimension reduction
method and its close relatives). Unfortunately it is not quite difficult to
observe that in general they both miss this point, that is, are unable to
reconstruct the embeddings of generic data up to some (almost) isometry, i.e.
with small distrortion of distances. This will be shown for both MDS and MVU
in section 3; moreover, based on the result from [1], we will show that the
most widely used classical MDS in fact, in a proper limit, may infinitely
distort the original distances when applied, say, to just a unit
circumference, by reconstructing in this case a snowflake-like closed curve
which is just Hölder continuous. Note that some newer heuristic methods of
dimension reduction with steadily growing popularity like SNE or t-SNE [7, 9]
also contain some core parts reconstructing data points just from the
distances; up to now they do not have any proven estimates on the distortion
of pairwise distances, but anyhow one cannot expect estimates of uniform
nature as it is clear from the methods themselves (see e.g. [3]).

We will further propose in section 5.1 a sightly different method (however
ideologically quite close to both MDS/Isomap and MVU) that reconstructs
isometric embeddings up to some small controllable error uniform over the
data. This will be based on a variational setting proposed in section 4 and on
a simple $\Gamma$-convergence result (Theorem 4.1). As an application in
section 6 we show that this method can be used also for topological data
reconstruction, i.e. for computing Čech cohomologies, and provide explicit
estimates on the input parameters for this purpose. Finally, in section 7 we
provide some numerical experiments on the proposed method.

##  2\. Notation and preliminaries

For a metric space $E$ equipped with distance $d$ and a curve
$\theta\colon[0,1]\to E$ we denote $|\dot{\theta}|$ its _metric derivative_
and by

$\ell(\theta):=\int_{0}^{1}|\dot{\theta}|(t)\,dt$

its _parametric length_. The notation $B_{r}(x)\subset E$ stands for the open
ball of radius $r>0$ with center $x\in E$. The Euclidean norm is denoted by
$|\cdot|$.

For a set $M\subset\mathbb{R}^{n}$ and $\varepsilon>0$ let
$(M)_{\varepsilon}\subset\mathbb{R}^{n}$ to be its open
$\varepsilon$-neighborhood, i.e. $(M)_{\varepsilon}:=\cup_{x\in
N}B_{\varepsilon}(x)$. We recall the notion of reach of $M$ introduced by
Federer in [5] and defined by

$\mathrm{Reach}\,(M):=\sup\left\\{\varepsilon>0\colon\mbox{every point of
$(M)_{\varepsilon}$ has a unique projection on $M$}\right\\}.$

For an $m\times n$ matrix $X$ we denote as usual by $X^{T}$ its transpose.
Vectors are silently identified with columns. The notation
$\mathrm{diag}\,(\alpha_{1},\ldots,\alpha_{r})$ stands for the diagonal matrix
with entries $(\alpha_{1},\ldots,\alpha_{r})$ over the diagonal. By $x\cdot y$
we denote the usual scalar product of vectors $x\in\mathbb{R}^{n}$ and
$y\in\mathbb{R}^{n}$.

##  3\. Main existing methods

The existing algorithms in manifold learning aimed at manifold reconstruction
from intrinsic distances are quite numerous, but can be essentially viewed as
close parents to just two basic ones, multidimensional scaling (MDS) and
maximum variance unfolding (MVU), which are aimed at reconstructing the
locations of points $y_{j}$ up to an isometric (or almost isometric, up to a
small error) embedding. This would be the case if the algorithm with input
$\\{d_{ij}\\}_{i,j=1}^{N}$, would produce the set of points
$\\{x_{i}^{N}\\}_{i=1}^{N}\subset\mathbb{R}^{n}$, and setting
$f_{N}(y_{i}):=x_{i}^{N}$, when $N\to\infty$, the functions $f_{N}$ tend to
some $f\colon M\to\mathbb{R}^{n}$ with $\Sigma:=f(M)$ (almost) isometric to
$M$. Unfortunately, as we show below, the existing methods in general miss
this point (though being applicable sometimes to particular situations).

###  3.1. Multidimensional scaling (MDS)

The classical multidimensional scaling (MDS) introduced by Torgerson and
further developed by many authors (see chapter 6 of [10] and references
therein), has been formulated for the situation when the distance $d_{M}$ is
Euclidean (which happens e.g. when $M$ is a convex subset of
$\mathbb{R}^{m}$). In this case, denote by $A$ the $N\times N$ matrix of
squared distances (called sometimes also proximity matrix) with entries
$A_{ij}:=d_{ij}^{2}$, and set

$B:=-\frac{1}{2}HAH,$

where $H$ is the $N\times N$ centering matrix
$H:=\mathrm{Id}_{N}-{\frac{1}{N}}\mathbf{1}\mathbf{1}^{T}$,
$\mathbf{1}\in\mathbb{R}^{N}$ standing for the column vector with $N$ unit
components, $\mathrm{Id}_{N}$ standing for the $N\times N$ identity matrix. It
is quite easy to show that $B$ is positive semidefinite with
$r:=\mathrm{rank}\,B\leq n$, and the decomposition $B=X^{T}X$ for an $r\times
N$ matrix $X$ allows to identify the points $x_{i}^{N}\in\mathbb{R}^{n}$,
$i=1,\ldots,N$, with the columns of $X$ [10, lemma 6.3]. Namely, denoting by
$\lambda_{i}$, $i=1,\ldots,r$, the nonzero eigenvalues of $B$ (numbered for
convenience in decreasing order, i.e.
$\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{r}>0$), and by
$\phi_{i}\in\mathbb{R}^{N}$ the corresponding eigenvectors, we may set
$X:=\Phi\,\mathrm{diag}\,(\sqrt{\lambda_{1}},\ldots,\sqrt{\lambda_{r}})$,
where $\Phi$ is the $r\times N$ matrix with columns $\phi_{i}$. In the recent
paper [2] one can find also the fine perturbation analysis of MDS as well as
on Isomap-type methods related to the latter, but only for the situation when
$M$ is isometric to a convex subset of a Euclidean space.

In practice however MDS method is quite frequently applied to the situation
when the distance $d_{M}$ is not necessarily Euclidean. For instance, it is
used as an essential part of the Isomap method proposed by Tenenbaum, de
Silva, and Langford for dimension reduction of a given set of points (see
chapter 8 of [10] and references therein), where it is applied to geodesic
distances in a specially constructed graph (in our setting, when the points
cover almost densely some regular submanifold of $\mathbb{R}^{n}$, these
distances in fact approximate the usual geodesic distances over this
submanifold). In the general case of not necessarily Euclidean distances one
just does exactly the same as in the Euclidean one but taking for
$\lambda_{i}$ the largest $r$ positive eigenvalues of $B$. This corresponds to
finding the points $x_{i}^{N}\in\mathbb{R}^{n}$ minimizing the stress
functional

$\mathrm{stress}(x_{1},x_{2},...,x_{N}):=\sum_{i,j=1}^{N}\left(d_{ij}^{2}-|x_{i}-x_{j}|^{2}\right),$

see theorem 6.3 in [10] (although there it is formally stated only in the case
all distances $d_{ij}$ are Euclidean). This however in general does not allow
to reconstruct the embedding of the original manifold $M$ up to an (almost)
isometry, as the following example shows.

######  Example 3.1.

Taking $M$ to be a unit circumference $S^{1}$, and $\\{y_{i}\\}_{i=1}^{N}$ to
be a set of equally spaced points in $M$, proposition 7.2.6 of [1] shows that
the MDS embedding of $\\{y_{i}\\}_{i=1}^{N}$ in $\mathbb{R}^{n}$ lies, up to a
rigid motion, on the closed curve $\gamma_{N}\colon[0,2\pi]\to\mathbb{R}^{n}$
defined by

$\displaystyle\gamma_{N}$ $\displaystyle(t):=$
$\displaystyle(a_{1}^{N}\cos(t),a_{1}^{N}\sin(t),\ldots,a_{2k+1}^{N}\cos((2k+1)t),a_{2k+1}^{N}\sin((2k+1)t),\ldots)\in\mathbb{R}^{n},$

where $\lim_{N}a_{j}^{N}=a_{j}:=\sqrt{2}/j$ (with $j$ odd). Clearly, in the
limit $N\to\infty$ and $n\to\infty$ this gives a closed snowflake-like curve
$\gamma$ homeomorphic but not isometric (nor even bilipschitz) to $M$; in
fact,

(3.1) $|\gamma(t)-\gamma(s)|=2\sqrt{\pi}|t-s|^{1/2}.$

To prove (3.1), we calculate

(3.2) $\displaystyle|\gamma(t)-\gamma(s)|^{2}$
$\displaystyle=\sum_{k=0}^{\infty}\frac{4}{(2k+1)^{2}}\left(\sin(2k+1)t-\sin(2k+1)s\right)^{2}$
$\displaystyle\qquad+\sum_{k=0}^{\infty}\frac{4}{(2k+1)^{2}}\left(\cos(2k+1)t-\cos(2k+1)s\right)^{2}$
$\displaystyle=16\sum_{k=0}^{\infty}\frac{1}{(2k+1)^{2}}-16\sum_{k=0}^{\infty}\frac{\cos\left((2k+1)(t-s)\right)}{(2k+1)^{2}}$
$\displaystyle=2\pi^{2}-16\sum_{k=0}^{\infty}\frac{\cos\left((2k+1)(t-s)\right)}{(2k+1)^{2}}.$

But from the Fourier expansion

$\displaystyle|x|=\frac{\pi}{2}-\frac{4}{\pi}\sum_{k=0}^{\infty}\frac{\cos(2k+1)x}{(2k+1)^{2}}$

for $x\in[-\pi,\pi]$, we get

$\displaystyle\sum_{k=0}^{\infty}\frac{\cos(2k+1)(t-s)}{(2k+1)^{2}}=\frac{\pi^{2}}{8}-\frac{\pi|t-s|}{4},$

which plugged in (3.2) gives (3.1) as claimed.

###  3.2. Maximum variance unfolding (MVU)

The method of maximum variance unfolding (MVU) (alternatively called also
semidefinite embedding (SDE)), has been introduced by Weinberger and Saul, see
chapter 9.1 of [10], and amounts to finding the points
$x_{i}^{N}\in\mathbb{R}^{n}$, $i=1,\ldots,N$ given the distance matrix
$\\{d_{ij}\\}_{i,j=1}^{N}$, by maximizing the total variance functional

$\mathrm{var}(x_{1},x_{2},...,x_{N}):=\sum_{i,j=1}^{N}|x_{i}-x_{j}|^{2}$

subject to the set of constraints

(3.3) $|x_{i}-x_{j}|^{2}=d_{ij}^{2}\quad\mbox{whenever $y_{i}$ is close to
$y_{j}$}.$

The condition of “$y_{i}$ close to $y_{j}$” is understood differently in
different versions of MVU, but most commonly as $d_{ij}\leq\varepsilon$ for
some fixed $\varepsilon>0$, so that (3.3) becomes

(3.4) $|x_{i}-x_{j}|^{2}=d_{ij}^{2}\quad\mbox{whenever
$d_{ij}\leq\varepsilon$}.$

The constraints (3.3) (or in particular (3.3)) are reformulated in an
equivalent way in terms of the Gram matrix $K$ with entries
$K_{ij}:=x_{i}\cdot x_{j}$ so that the above maximization becomes semidefinite
programming problem.

Similarly to Example 3.1 it is easy to show that MVU, and even more, any
method trying to preserve locally the distances as Euclidean ones, in general
not only does not allow to reconstruct the embedding of the original manifold
$M$ up to an (almost) isometry, but even worse, the constraints (3.4) may not
allow to reconstruct even something vaguely similar to $M$, as the example
below shows.

######  Example 3.2.

Taking again $M$ to be a unit circumference $S^{1}$, and
$\\{y_{i}\\}_{i=1}^{N}$ to be a set of equally spaced points in $M$, suppose
that $\\{x_{i}^{N}\\}_{i=1}^{N}\subset\mathbb{R}^{n}$ satisfy (3.4), i.e

$|x_{i}^{N}-x_{j}^{N}|^{2}=d_{ij}^{2}\quad\mbox{whenever
$d_{ij}\leq\varepsilon$}$

for some fixed $\varepsilon>0$, and that continuous functions $f_{N}\colon
M\to\mathbb{R}^{n}$ satisfying $f_{N}(y_{j})=x_{j}^{N}$ for all $j=1,\ldots,N$
converge as $N\to\infty$ to some continuous function $f\colon
M\to\mathbb{R}^{n}$. Then one has

$|f(u)-f(v)|=d_{M}(u,v)\quad\mbox{whenever $d_{M}(u,v)<\varepsilon$}.$

Parameterizing $M$ in a natural way over $[0,2\pi]$ by a curve
$\theta\colon[0,2\pi]\to M$, $\theta(t):=(\cos t,\sin t)$, for the curve
$\gamma\colon[0,2\pi]\to\mathbb{R}^{n}$ defined by $\gamma(t):=f(\theta(t))$
we have therefore that

$|\gamma(t)-\gamma(s)|=|t-s|$

whenever $|t-s|<\varepsilon$. Thus $\gamma(I)$ is a line segment for every
interval $I\subset[0,2\pi]$ of length $\ell(I)<2\pi$, which implies that
$\gamma([0,2\pi))$ is a nondegenerate line segment. But on the other hand one
must have $\gamma(0)=\gamma(2\pi)$, which is a contradiction.

##  4\. Variational setting

From now on we assume $M\subset\mathbb{R}^{n}$ to be a smooth connected
compact Riemannian submanifold of $\mathbb{R}^{n}$, with
$\mathrm{Reach}\,(M)=\alpha>0$ (which is of course automatically guaranteed
when $M$ is at least $C^{2}$-smooth). Let $\Sigma_{k}\subset M$ be a sequence
of closed sets such that the lower Kuratowski limit of $\Sigma_{k}$ is $M$,
i.e. for every $x\in M$ there is a sequence $x_{k}\in\Sigma_{k}$ such that
$\lim_{k}x_{k}=x$ (in particular, this holds when $\Sigma_{k}\to M$ as
$k\to\infty$ in the Hausdorff distance $d_{H}$).

Given an $\varepsilon>0$ and a $k\in\mathbb{N}$, we define the functionals

$F_{\varepsilon,k}\colon
C(M;\mathbb{R}^{n})\to\mathbb{R},F_{\varepsilon}\colon
C(M;\mathbb{R}^{n})\to\mathbb{R},$

by the formulae

(4.1) $\displaystyle F_{\varepsilon,k}(f)$
$\displaystyle:=\sup\left\\{\left|\dfrac{|f(x)-f(y)|^{2}}{d_{M}^{2}(x,y)}-1\right|\colon\\{x,y\\}\subset\Sigma_{k},0<d_{M}(x,y)\leq\varepsilon\right\\}.$
$\displaystyle F_{\varepsilon}(f)$
$\displaystyle:=\sup\left\\{\left|\dfrac{|f(x)-f(y)|^{2}}{d_{M}^{2}(x,y)}-1\right|\colon\\{x,y\\}\subset
M,0<d_{M}(x,y)\leq\varepsilon\right\\}.$

The scope of this section is to prove the following easy result.

######  Theorem 4.1.

Let $M$ be compact and connected, $\varepsilon_{0}:=\alpha$, $C_{1}$ and
$C_{2}$ be as in Lemmata A.1, A.2, $x_{0}\in M$, and denote

$\displaystyle\mathcal{C}_{k}$ $\displaystyle:=\left\\{f\in
C(M;\mathbb{R}^{n}),f(x_{0})=0,|f(x)-f(y)|\geq C_{2}d_{M}(x,y)\quad\mbox{for
all}\\{x,y\\}\subset\Sigma_{k}\right\\}.$ $\displaystyle\mathcal{C}$
$\displaystyle:=\left\\{f\in C(M;\mathbb{R}^{n}),f(x_{0})=0,|f(x)-f(y)|\geq
C_{2}d_{M}(x,y)\right\\}.$

Then the variational problems

$None$ $\min\left\\{F_{\varepsilon,k}(f)\colon f\in\mathcal{C}_{k}\right\\}$

have solutions for every $k\in\mathbb{N}$. If $f_{k}$ is a solution to
$(P_{k})$, then up to a subsequence (not relabeled) one has $\lim_{k}f_{k}=f$,
where $f$ solves

$None$ $\min\left\\{F_{\varepsilon}(f)\colon f\in\mathcal{C}\right\\}$

and satisfies

(4.2) $d_{M}(x,y)(1-3C_{1}\varepsilon)\leq|f(x)-f(y)|\leq
d_{M}(x,y)(1+3C_{1}\varepsilon),$

if $d_{M}(x,y)<\varepsilon$, and

(4.3) $d_{M}(x,y)(1-3C_{1}\varepsilon)\leq d_{\Sigma}(f(x),f(y))\leq
d_{M}(x,y)(1+3C_{1}\varepsilon),$

where $\Sigma:=f(M)$ and $d_{\Sigma}$ stands for the geodesic distance in
$\Sigma$, i.e.

$d_{\Sigma}(u,v):=\inf\left\\{\ell(\theta)\colon\theta\colon[0,1]\to\Sigma,\theta(0)=u,\theta(1)=v\right\\},$

$\ell(\theta)$ standing for the Euclidean length of $\theta$. Note that
$f_{k}(M)\to\Sigma$ in the sense of Hausdorff distance as $k\to+\infty$.

######  Remark 4.2.

Clearly, the map $f$ solving $(P)$ satisfies also

(4.4) $|f(x)-f(y)|\geq C_{2}d_{M}(x,y)$

for generic $(x,y)\subset M\times M$.

###### Proof.

By Lemma 4.4 for every $\varepsilon>0$ there is a $\bar{k}\in\mathbb{N}$
depending only on $\varepsilon$ such that for every $k\geq\bar{k}$ all $f\in
C(M;\mathbb{R}^{n})$ satisfying $F_{\varepsilon,k}(f)\leq K$ have equibounded
Lipschitz constants over $\Sigma_{k}$, and hence, up to redefining each $f$
over $M\setminus\Sigma_{k}$ as an extension with minimum Lipschitz constant
from $\Sigma_{k}$ to $M$, one has that the Lipschitz constants
$\mathrm{Lip}\,f$ over $M$ are equibounded for all such $f$, so that the set
$\cup_{k\geq\bar{k}}\\{f\in\mathcal{C}\colon F_{\varepsilon,k}(f)\leq K\\}$ is
compact by Ascoli-Arzelà theorem.

The class $\mathcal{C}$ is closed and by Lemma A.2 contains the translation
map $x\mapsto x-x_{0}$. Also the classes $\mathcal{C}_{k}$ are closed and,
moreover, $f\in\mathcal{C}$ once $f=\lim_{k}f_{k}$ with
$f_{k}\in\mathcal{C}_{k}$. Thus from Lemma 4.3 and the main property of
$\Gamma$-convergence minimizers $f_{k}$ solving $(P_{k})$ (which exist because
of coercivity and continuity of the functionals $F_{\varepsilon,k}$), up to a
subsequence (not relabeled) converge to a minimizer $f$ of $(P)$. By Lemma 4.5
one has (4.2) for $d_{M}(x,y)<\varepsilon$, and by Lemma 4.6 (with
$C:=1+3C_{1}\varepsilon$, $c:=1-3C_{1}\varepsilon$) one has (4.3). Finally,
note that when $\lim_{k}x_{k}=x$ in $M$, then $\lim_{k}f_{k}(x_{k})=f(x)$ in
view of the uniform convergence of $f_{k}$, which implies that
$f_{k}(M)\to\Sigma$ in the sense of Hausdorff distance as $k\to+\infty$. ∎

The following technical assertions have been used in the above proof.

######  Lemma 4.3.

Let $\mathcal{C}\subset C(M;\mathbb{R}^{n})$, $\mathcal{C}_{k}\subset
C(M;\mathbb{R}^{n})$ be closed sets of maps satisfying $f\in\mathcal{C}$ once
$f=\lim_{k}f_{k}$ with $f_{k}\in\mathcal{C}_{k}$, and

$\chi_{\mathcal{C}}(f):=\left\\{\begin{array}[]{rl}0,&f\in\mathcal{C},\\\
+\infty,&f\not\in\mathcal{C}.\end{array}\right.$

One has
$\Gamma^{-}\lim_{k}(F_{\varepsilon,k}+\chi_{\mathcal{C}_{k}})=F_{\varepsilon}+\chi_{\mathcal{C}}$,
where

$F_{\varepsilon}\colon C(M;\mathbb{R}^{n})\to[0,+\infty]$

is defined by (4.1).

###### Proof.

One has

(4.5) $\displaystyle\liminf_{k}$ $\displaystyle F_{\varepsilon,k}(f_{k})$
$\displaystyle=\liminf_{k}\sup\left\\{\left|\dfrac{|f_{k}(x)-f_{k}(y)|^{2}}{d_{M}^{2}(x,y)}-1\right|\colon\\{x,y\\}\subset\Sigma_{k},0<d_{M}(x,y)\leq\varepsilon\right\\}$
$\displaystyle\geq\liminf_{k}\sup\left\\{\left|\dfrac{|f_{k}(x)-f_{k}(y)|^{2}}{d_{M}^{2}(x,y)}-1\right|\colon\\{x,y\\}\subset\Sigma_{k},\rho\leq
d_{M}(x,y)\leq\varepsilon\right\\}$

for every $\rho>0$. Let $\\{\tilde{x},\tilde{y}\\}\in M$ be such that
$\rho\leq d_{M}(\tilde{x},\tilde{y})\leq\varepsilon$ and

$\left|\dfrac{|f(\tilde{x})-f(\tilde{y})|^{2}}{d_{M}^{2}(\tilde{x},\tilde{y})}-1\right|=\sup\left\\{\left|\dfrac{|f(x)-f(y)|^{2}}{d_{M}^{2}(x,y)}-1\right|\colon\\{x,y\\}\subset
M,\rho\leq d_{M}(x,y)\leq\varepsilon\right\\}.$

Since one has the existence of $\\{x_{k},y_{k}\\}\in\Sigma_{k}$ with
$\lim_{k}x_{k}=\tilde{x}$, $\lim_{k}y_{k}=\tilde{y}$, then

(4.6) $\displaystyle\liminf_{k}$
$\displaystyle\sup\left\\{\left|\dfrac{|f_{k}(x)-f_{k}(y)|^{2}}{d_{M}^{2}(x,y)}-1\right|\colon\\{x,y\\}\subset\Sigma_{k},\rho\leq
d_{M}(x,y)\leq\varepsilon\right\\}$
$\displaystyle\geq\lim_{k}\left|\dfrac{|f_{k}(x_{k})-f_{k}(y_{k})|^{2}}{d_{M}^{2}(x_{k},y_{k})}-1\right|$
$\displaystyle=\sup\left\\{\left|\dfrac{|f(x)-f(y)|^{2}}{d_{M}^{2}(x,y)}-1\right|\colon\\{x,y\\}\subset
M,\rho\leq d_{M}(x,y)\leq\varepsilon\right\\},$

Thus combining (4.5) and (4.6), we get

$\liminf_{k}F_{\varepsilon,k}(f_{k})\geq\sup\left\\{\left|\dfrac{|f(x)-f(y)|^{2}}{d_{M}^{2}(x,y)}-1\right|\colon\\{x,y\\}\subset
M,\rho\leq d_{M}(x,y)\leq\varepsilon\right\\}$

for every $\rho>0$. Hence taking $\sup$ in the above estimate over
$\rho\in(0,\varepsilon)$, we get

$\displaystyle\liminf_{k}F_{\varepsilon,k}(f_{k})$
$\displaystyle\geq\sup\left\\{\left|\dfrac{|f(x)-f(y)|^{2}}{d_{M}^{2}(x,y)}-1\right|\colon\\{x,y\\}\subset
M,0<d_{M}(x,y)\leq\varepsilon\right\\}$ $\displaystyle=F_{\varepsilon}(f),$

which means

$\Gamma^{-}\liminf_{k}F_{\varepsilon,k}\geq F_{\varepsilon},$

Thus the inequality

$\Gamma^{-}\liminf_{k}(F_{\varepsilon,k}+\chi_{\mathcal{C}_{k}})\geq
F_{\varepsilon}+\chi_{\mathcal{C}}$

follows because of the requirement on $\mathcal{C}_{k}$.

The inequality

$\Gamma^{-}\limsup_{k}(F_{\varepsilon,k}+\chi_{\mathcal{C}_{k}})\leq
F_{\varepsilon}+\chi_{\mathcal{C}}$

is immediate since $F_{\varepsilon,k}(f)\leq F_{\varepsilon}(f)$ for every
$f\in C(M;\mathbb{R}^{n})$, and $\mathcal{C}\subset\mathcal{C}_{k}$ which
concludes the proof. ∎

######  Lemma 4.4.

If $M$ is compact and connected, then for every $\varepsilon>0$ there is a
$\bar{k}\in\mathbb{N}$ (depending only on $\varepsilon$) such that for every
$k\geq\bar{k}$ and for every $f\in C(M;\mathbb{R}^{n})$ satisfying
$F_{\varepsilon,k}(f)\leq K$ one has

$|f(x)-f(y)|\leq Cd_{M}(x,y)$

for all $\\{x,y\\}\in\Sigma_{k}$, $x\neq y$, where $C>0$ depends only on $K$.

###### Proof.

If $f\in C(M;\mathbb{R}^{n})$ satisfies $F_{\varepsilon,k}(f)\leq K$, then

(4.7) $\displaystyle\frac{|f(x)-f(y)|}{d_{M}(x,y)}\leq C$

for all $\\{x,y\\}\in\Sigma_{k}$, $0<d_{M}(x,y)<\varepsilon$ (with
$C:=\sqrt{K+1}$).

Consider a finite $\varepsilon$-net $\\{x_{j}\\}_{j=1}^{N}$ of $M$. For each
pair of indices $i,j=1,\ldots,N$ fix a “polygonal line” $\theta_{ij}$ in $M$
connecting $x_{i}$ to $x_{j}$ made of a finite number of geodesic segments of
lengths at most $\varepsilon/2$ with vertices on a geodesic curve between
$x_{i}$ and $x_{j}$ and such that $\ell(\theta_{ij})<2d_{M}(x_{i},x_{j})$.
Note that the set of all vertices $\mathcal{V}$ of all $\theta_{ij}$ is
clearly finite.

Choose a $\bar{k}\in\mathbb{N}$ such that for every $k\geq\bar{k}$ for every
pair $\\{u,v\\}\in\mathcal{V}$ of consecutive vertices of some $\theta_{ij}$
there are points $\\{u_{k},v_{k}\\}\in\Sigma_{k}$ with

$d_{M}(u_{k},u)<\frac{\varepsilon}{2},\quad
d_{M}(v_{k},v)<\frac{\varepsilon}{2},\quad d_{M}(u_{k},v_{k})<2d_{M}(u,v)$

(the latter inequality implies in particular that
$d_{M}(u_{k},v_{k})<\varepsilon$). Therefore, for every $\theta_{ij}$ there is
a polygonal line $\sigma_{ij}$ with the same number of vertices as in
$\theta_{ij}$, all vertices in $\Sigma_{k}$, each of its vertices being at
most $\varepsilon/2$ close (in $d_{M}$) to the respective vertex of
$\theta_{ij}$, geodesic distance between consecutive vertices at most
$\varepsilon$, and, finally,

$\ell(\sigma_{ij})<2\ell(\theta_{ij})<4d_{M}(x_{i},x_{j}).$

Denoting the consequtive vertices of $\sigma_{ij}$ by $y_{0},\ldots,y_{m}$ (of
course, $m$ depends on $i$ and $j$), we get

(4.8) $\displaystyle|f(\sigma_{ij}(0))-f(\sigma_{ij}(1))|$
$\displaystyle=|f(y_{0})-f(y_{m})|\leq\sum_{j=0}^{m-1}|f(y_{j})-f(y_{j+1})|$
$\displaystyle\leq
C\sum_{j=0}^{m-1}d_{M}(y_{j},y_{j+1})=C\ell(\sigma_{ij})<4Cd_{M}(x_{i},x_{j}).$

Finally, for arbitrary $\\{x,y\\}\in\Sigma_{k}$, $d_{M}(x,y)\geq\varepsilon$,
we find a couple $\\{x_{i},x_{j}\\}$ such that

$d_{M}(x_{i},x)<\varepsilon,\quad d_{M}(x_{j},y)<\varepsilon,$

and estimate

(4.9) $\displaystyle|f(x)-f(y)|$
$\displaystyle\leq|f(x)-f(\sigma_{ij}(0))|+|f(\sigma_{ij}(0))-f(\sigma_{ij}(1))|+|f(\sigma_{ij}(1))-f(y)|$
$\displaystyle<Cd_{M}(x,\sigma_{ij}(0))+4Cd_{M}(x_{i},x_{j})+Cd_{M}(\sigma_{ij}(1),y)$

in view of (4.8). But

(4.10) $\displaystyle d_{M}(x,\sigma_{ij}(0))$ $\displaystyle\leq
d_{M}(x,x_{i})+d_{M}(x_{i},\sigma_{ij}(0))\leq 3\varepsilon/2,$

and analogously

(4.11) $\displaystyle d_{M}(y,\sigma_{ij}(1))$ $\displaystyle\leq
3\varepsilon/2,$

while

(4.12) $\displaystyle d_{M}(x_{i},x_{j})$ $\displaystyle\leq
d_{M}(x_{i},x)+d_{M}(x,y)+d_{M}(y,x_{j})\leq 2\varepsilon+d_{M}(x,y).$

Plugging (4.10), (4.11) and (4.12) into (4.9), we get

(4.13) $\displaystyle|f(x)-f(y)|$
$\displaystyle<C(11\varepsilon+4d_{M}(x,y))\leq 15Cd_{M}(x,y),$

because $d_{M}(x,y)\geq\varepsilon$. Together with (4.7) the estimate (4.13)
proves the claim. ∎

######  Lemma 4.5.

Let $\varepsilon_{0}$ and $C_{1}$ be as in Lemma A.1. If $f\in
C(M;\mathbb{R}^{n})$ is a minimizer of $F_{\varepsilon}$ with
$\varepsilon<\varepsilon_{0}\wedge 1/3C_{1}$ over some class containing some
rigid translation, then

$d_{M}(x,y)(1-3C_{1}\varepsilon)\leq|f(x)-f(y)|\leq
d_{M}(x,y)(1+3C_{1}\varepsilon),$

if $d_{M}(x,y)<\varepsilon$. The same holds for $\\{x,y\\}\subset\Sigma_{k}$
when $f\in C(M;\mathbb{R}^{n})$ is a minimizer of $F_{\varepsilon,k}$.

###### Proof.

Note that $F_{\varepsilon}(f)$ is invariant with respect to the compositions
of $f$ with rigid translations, and in particular the value of
$F_{\varepsilon}$ over any rigid translation is equal to that over the
identity map $\mathrm{id}$. One has then

$F_{\varepsilon}(f)\leq
F_{\varepsilon}(\mathrm{id})\leq(1+C_{1}\varepsilon)^{2}-1\leq
3C_{1}\varepsilon,$

the second inequality being valid due to Lemma A.1, which implies

$\left|\dfrac{|x-y|^{2}}{d_{M}^{2}(x,y)}-1\right|\leq(1+C_{1}\varepsilon)^{2}-1$

whenever $d_{M}(x,y)<\varepsilon$, showing the statement for
$F_{\varepsilon}$. The proof for $F_{\varepsilon,k}$ is identical. ∎

######  Lemma 4.6.

Suppose that for some $C>0$ and $f\in C(M;\mathbb{R}^{n})$ one has

(4.14) $|f(x)-f(y)|\leq Cd_{M}(x,y),$

if $d_{M}(x,y)<\varepsilon$. Then

(4.15) $d_{\Sigma}(f(x),f(y))\leq Cd_{M}(x,y),$

where $\Sigma:=f(M)$. If, moreover, $f$ is injective with the inverse function
$f^{-1}\colon\Sigma\to M$ continuous (which is the case, e.g. when $f$ is
proper), and

(4.16) $cd_{M}(x,y)\leq|f(x)-f(y)|$

for some $c>0$, if $d_{M}(x,y)<\varepsilon$, then also

(4.17) $cd_{M}(x,y)\leq d_{\Sigma}(f(x),f(y)).$

######  Remark 4.7.

The estimate (4.17) cannot hold for $f$ just satisfying (4.14) and (4.16) only
for $d_{M}(x,y)<\varepsilon$, unless $f$ is injective, as can be seen from the
example of $M$ a line segment of length $2\pi$ (identified with $[0,2\pi]$)
and $f\colon[0,2\pi]\to\mathbb{R}^{2}$ defined by $f(t):=(\cos t,\sin t)$.

###### Proof.

If $\theta\colon[0,1]\to M$ is a Lipschitz curve, then (4.14) implies that so
is $\sigma:=f\circ\theta\colon[0,1]\to\mathbb{R}^{n}$ and its metric
derivative $|\dot{\sigma}|$ satisfies

(4.18) $|\dot{\sigma}|(t)\leq C|\dot{\theta}|(t)$

for a.e. $t\in[0,1]$. Then (4.18) gives (4.15). In fact, if $\theta$ is a
geodesic curve connecting $x=\theta(0)$ to $x=\theta(0)$, then

$d_{\Sigma}(f(x),f(y))\leq\int_{0}^{1}|\dot{\sigma}|(t)\,dt\leq
C\int_{0}^{1}|\dot{\theta}|(t)\,dt=Cd_{M}(x,y).$

If $f$ is injective, then take an arbitrary $\delta>0$, and consider a
rectifiable curve $\sigma\colon[0,1]\to f(M)\subset\mathbb{R}^{n}$ such that
$\sigma(0)=f(x)$, $\sigma(1)=f(y)$ with

$\int_{0}^{1}|\dot{\sigma}|(t)\,dt\leq d_{\Sigma}(f(x),f(y))+\delta.$

Denote $\theta(t):=f^{-1}(\sigma(t))$ for every $t\in[0,1]$. If $f^{-1}$ is
continuous, then so is $\theta$, and hence for every $t\in[0,1]$ one has
$d_{M}(\theta(t),\theta(t+s))<\varepsilon$ once $s$ is sufficiently small.
Thus from (4.15) we get

(4.19) $c|\dot{\theta}|(t)\leq|\dot{\sigma}|(t)$

for a.e. $t\in[0,1]$. Therefore,

$cd_{M}(x,y)\leq c\int_{0}^{1}|\dot{\theta}|(t)\,dt\leq
d_{\Sigma}(f(x),f(y))+\delta,$

and taking the limit in the above inequality as $\delta\to 0^{+}$, we arrive
at (4.17). ∎

##  5\. Discrete variational setting

Let $\\{y_{i}\\}\subset M$ be a dense set in $M$, and denote for the sake of
brevity

$d_{ij}:=d_{M}(y_{i},y_{j}).$

Given an $\varepsilon>0$ and a $k\in\mathbb{N}$, we define the functional
$F_{\varepsilon,k}\colon(\mathbb{R}^{n})^{k}\to\mathbb{R}$ by the formula

$F_{\varepsilon,k}(x_{1},\ldots,x_{k}):=\max\left\\{\left|\dfrac{|x_{i}-x_{j}|^{2}}{d_{ij}^{2}}-1\right|\colon
i,j=1,\ldots,k,i\neq j,d_{ij}<\varepsilon\right\\}.$

The following statement is just a “translation” of Theorem 4.1 to this
sequence.

######  Proposition 5.1.

Let $\varepsilon_{0}:=\alpha$, $C_{1}$ and $C_{2}$ be as in Lemmata A.1, A.2.
Assume that $(x_{i}^{k})_{i=1}^{k}\in(\mathbb{R}^{n})^{k}$ be a minimizer of
$F_{\varepsilon,k}$ with $\varepsilon<\varepsilon_{0}\wedge 1/(C_{1}+1)$ over
the set $X^{k}\subset(\mathbb{R}^{n})^{k}$ defined by

$X^{k}:=\left\\{((x_{i})_{i=1}^{k}\in(\mathbb{R}^{n})^{k}\colon|x_{i}-x_{j}|\geq
C_{2}d_{ij},i,j=1,\ldots,k\right\\}.$

Then up to a subsequence $x_{i}^{k}\to x_{i}$ as $k\to\infty$, and

(5.1) $d_{ij}(1-C_{1}\varepsilon)\leq|x_{i}-x_{j}|\leq
d_{ij}(1+C_{1}\varepsilon),$

whenever $d_{ij}<\varepsilon$. Further, there is a $C>0$ depending only on
$\alpha$ and $D$ such that

$d_{ij}(1-C\varepsilon)\leq d_{\Sigma}(x_{i},x_{j})\leq
d_{ij}(1+C\varepsilon)$

for all $\\{i,j\\}\subset\mathbb{N}$.

###  5.1. Reduction to semidefinite programming problem

The problem of minimizing $F_{\varepsilon,k}$ over the set
$X^{k}\subset(\mathbb{R}^{n})^{k}$ is written as minimizing the convex
function $G_{\varepsilon,k}$ of a matrix defined by

$\displaystyle
G_{\varepsilon,k}(x_{1},\ldots,x_{k}):=\max\left\\{\left|\dfrac{K_{ii}+K_{jj}-2K_{ij}}{d_{ij}^{2}}-1\right|\colon
i,j=1,\ldots,k,i\neq j,d_{ij}<\varepsilon\right\\}.$

over the set of positive semidefinite matrices $K$ satisying the set of convex
constraints

$K_{ii}+K_{jj}-2K_{ij}\geq C_{2}^{2}d_{ij}^{2},i,j=1,\ldots,k.$

The solution $K$ of the latter problem is the Gram matrix of a set of vectors
$\\{x_{i}\\}$, i.e. $K_{ij}=\ x_{i}\cdot x_{j}$, $i,j=1,\ldots,k$, which
minimize $F_{\varepsilon,k}$ over $X^{k}$.

Adding a new scalar variable $t\in\mathbb{R}$ one reduces the above problem to
the following _semidefinite programming_ problem (i.e. a problem of
minimization of a linear function with linear constraints over the cone of
positive semidefinite matrices), namely

(5.2) minimize $\displaystyle\quad t\quad\mbox{over the
pairs}\quad(t,K)\quad\mbox{subject to}$ $\displaystyle-td_{ij}^{2}\leq
K_{ii}+K_{jj}-2K_{ij}-d_{ij}^{2}\leq td_{ij}^{2},$
$\displaystyle\qquad\mbox{for all }i,j=1,\ldots,k,i\neq j,d_{ij}<\varepsilon,$
$\displaystyle K_{ii}+K_{jj}-2K_{ij}\geq C_{2}^{2}d_{ij}^{2},$
$\displaystyle\qquad\mbox{for all }i,j=1,\ldots,k,i\neq j,$ $\displaystyle
K\quad\mbox{positive semidefinite $k\times k$ matrix}.$

##  6\. How to compute Čech cohomologies

Since $M$ and $\Sigma:=f(M)$ are homeomorphic (even biLipschitz equivalent) by
Theorem 4.1, they have the same homologies and cohomologies for every
reasonable (co)homology theory. Here as it is quite usual in topological data
analysis we are interested in Čech cohomologies of $M$ and provide a method to
compute the latter (recall that $M$ is unknown).

######  Lemma 6.1.

Assume that the conditions of Proposition 5.1 be satisfied. Then, if
$X\in\Sigma\cap B_{\sigma}(\bar{X})$, $\bar{X}\in\Sigma$, and $X=f(Y)$, $Y\in
M$, $\sigma\leq C_{2}\alpha$, then

$Y\in B_{\sigma^{\prime}}(\bar{Y}),\quad\mbox{where
}\sigma^{\prime}:=\left(1+C_{1}\frac{\sigma}{C_{2}}\right)\frac{\sigma}{C_{2}}.$

Vice versa, if $Y\in M\cap B_{\sigma^{\prime\prime}}(\bar{Y})$, $\bar{Y}\in
M$, where $\sigma^{\prime\prime}>0$ is defined by the equation

(6.1)
$\left(1+C_{1}\frac{\sigma^{\prime\prime}}{C_{2}}\right)\frac{\sigma^{\prime\prime}}{C_{2}}=\sigma,$

and $\sigma^{\prime\prime}\leq C_{2}\varepsilon$, then for $X=f(Y)$ one has
$X\in\Sigma\cap B_{\sigma}(\bar{X})$.

######  Remark 6.2.

It is worth noting that $\sigma^{\prime\prime}<\sigma<\sigma^{\prime}$.

###### Proof.

If $|X-\bar{X}|<\sigma$, then

$d_{M}(Y,\bar{Y})\leq\frac{\sigma}{C_{2}}$

in view of (4.4), so that the the first claim follows from Lemma A.1.

To prove the second claim note that if $|Y-\bar{Y}|<\sigma^{\prime\prime}$,
then by Lemma A.2 one has

$d_{M}(Y,\bar{Y})\leq\frac{\sigma^{\prime\prime}}{C_{2}},$

and hence it suffices to apply (4.2) minding the definition of
$\sigma^{\prime\prime}$. ∎

Let $\Lambda\subset\mathbb{N}$ be a finite set of indices such that for all
$Y\in M$ there is a $\lambda\in\Lambda$ and an $\bar{Y}_{\lambda}\in M$
satisfying

$d_{M}(Y,\bar{Y}_{\lambda})\leq\delta$

for some $\delta>0$. In other words, $\\{\bar{Y}_{\lambda}\\}$ is a finite
$\delta$-net of of $M$ (equipped with $d_{M}$).

Denote now by $C_{\Sigma}(r)$ the Čech complex built on the Euclidean balls
$B_{r}(\bar{X}_{\lambda})$, where $\bar{X}_{\lambda}:=f(\bar{Y}_{\lambda})$,
and by $C_{M}(r)$) the Čech complex built on the Euclidean balls
$B_{r}(\bar{Y}_{\lambda})$. We note that the vertices of all these complexes
may be considered the same (namely, the set of vertices of all them may be
identified with the index set $\Lambda$).

######  Lemma 6.3.

Under the conditions of Proposition 5.1 one has

$C_{M}(\sigma^{\prime\prime})\subset C_{\Sigma}(\sigma)\subset
C_{M}(\sigma^{\prime})$

when $0<\sigma\leq C_{2}\varepsilon_{0}$ and $\sigma^{\prime\prime}$ is
defined by (6.1) with $0<\sigma^{\prime\prime}\leq C_{2}\varepsilon_{0}$.

###### Proof.

Follows immediately from Lemma 6.1. ∎

######  Proposition 6.4.

Under the conditions of Proposition 5.1 assume $0<\sigma\leq C_{2}\alpha$ to
be so small that

$\sigma^{\prime}:=\left(1+C_{1}\frac{\sigma}{C_{2}}\right)\frac{\sigma}{C_{2}}<\frac{1}{2}\sqrt{\frac{3}{5}}\alpha,$

$\sigma^{\prime\prime}$ defined by (6.1) satisfy

$0<\sigma^{\prime\prime}<C_{2}\varepsilon$

and let $\delta\in(0,\varepsilon_{0})$ be such that

$\rho(\delta):=(1+C_{1}\delta)\delta\leq\sigma^{\prime}.$

Then $H^{*}(C_{\Sigma}(\sigma);\mathbb{R})\simeq H^{*}(M;\mathbb{R})$, $H^{*}$
standing for the Čech cohomology.

######  Remark 6.5.

One may take $\bar{Y}_{\lambda}$ to be drawn by sampling $M$ in i.i.d. way
according to the volume measure on $M$. In fact by proposition 3.2 of [8] one
has then that if $\\#\Lambda>n(M,\rho,p)$, then

$M\subset\bigcup_{\lambda}B_{\rho}(\bar{Y}_{\lambda})$

with probability at least $1-p$ and the number $n(M,\rho,p)$ depends
explicitly, besides $\rho$ and $p$, also on the total volume and the dimension
of $M$.

###### Proof.

From Lemma A.1 one has

$M\subset\bigcup_{\lambda}B_{\rho(\delta)}(\bar{Y}_{\lambda}),$

since $\delta<\varepsilon_{0}$. Therefore,

$M\subset\bigcup_{\lambda}B_{\sigma^{\prime}}(\bar{Y}_{\lambda}),$

since $\sigma^{\prime}>\rho(\delta)$. But

$2\sigma^{\prime\prime}<2\sigma^{\prime}<\sqrt{\frac{3}{5}}\lambda$

implies $H^{*}(C_{M}(\sigma^{\prime});\mathbb{R})\simeq
H^{*}((M)_{\sigma^{\prime}};\mathbb{R})$ and
$H^{*}(C_{M}(\sigma^{\prime\prime});\mathbb{R})\simeq
H^{*}((M)_{\sigma^{\prime\prime}};\mathbb{R})$ by proposition 3.1 of [8]. But
$(M)_{\sigma^{\prime\prime}}\subset(M)_{\sigma^{\prime}}\subset(M)_{\alpha}$,
then both $(M)_{\sigma^{\prime\prime}}$ and $(M)_{\sigma^{\prime}}$ are
clearly retractible to $M$, so that
$H^{*}((M)_{\sigma^{\prime\prime}};\mathbb{R})\simeq H^{*}(M;\mathbb{R})$ and
$H^{*}((M)_{\sigma^{\prime}};\mathbb{R})\simeq H^{*}(M;\mathbb{R})$.
Therefore,

$H^{*}(C_{M}(\sigma^{\prime});\mathbb{R})\simeq H^{*}(M;\mathbb{R})\simeq
H^{*}(C_{M}(\sigma^{\prime});\mathbb{R}),$

and it suffices now to apply Lemma 6.3. ∎

##  7\. Numerical experiments

In this section we present numerical experiments illustrating the performance
of our procedure for three sample datasets: a line segment, a two-dimensional
sphere and the Swiss Roll. In all the experiments, the constant $C_{2}$ from
(5.2) is set to $2/\pi$. For quantitative measure of the performance, we
compute the error

(7.1)
$Err=\frac{1}{k}\sqrt{\sum\limits_{i,j=1}^{k}\frac{\left(\widehat{d}_{ij}-d_{ij}\right)^{2}}{d_{ij}^{2}}},$

which reflects the average relative error in pairwise distances. Here
$\widehat{d}_{ij}$, $1\leq i,j\leq k$, stands for the pairwise distance
between the recovered i-th and j-th elements of the sample.

We start with the example of a line segment. We took $k=100$ equidistant
points on the unit interval $[0,1]$ and embedded them into $\mathbb{R}^{2}$
using our algorithm with $n=2$ and $\varepsilon=0.2$. The result is shown in
Figure 7.1. Though the recovered points do not lie on a segment, the error
(7.1) is equal to $Err=9\cdot 10^{-4}$, which is quite small.

Figure 7.1. Reconstruction of a line segment from pairwise distances. The
average distance error is $9\cdot 10^{-4}$.

In the example of a two-dimensional sphere we have two different setups. In
the first one, we took $k=100$ points on a grid on unit sphere $S^{2}$,
computed exact geodesic distances and applied the procedure with parameters
$n=3$ and $\varepsilon=0.6$. After that, we computed approximate pairwise
geodesic distances over the resulting point cloud as it is done is Isomap. As
a result, we obtained $Err=0.13$. In the second setup, we had $k=100$ points
drawn independently from uniform distribution on the sphere and computed exact
geodesic distances between them. After that, we performed the embedding into
$\mathbb{R}^{3}$ using our procedure with parameters $n=3$ and
$\varepsilon=0.6$, and computed approximate pairwise geodesic distances
between the embedded points again using the same method as in Isomap. As a
result, we obtained $Err=0.16$. The results of the sphere embedding are
displayed in Figure 7.2.

Figure 7.2. Reconstruction of a unit sphere from pairwise distances. Top line,
columns 1 and 2: points on a grid on the sphere. Top line, columns 3 and 4:
the recovered points of the unit sphere. Bottom line, columns 1 and 2: points
drawn from the uniform distribution on unit sphere. Bottom line, columns 3 and
4: the recovered points from approximate geodesic distances.

Finally, we provide the experiments on the widely known synthetic Swiss Roll
dataset from the Scikit-learn library in Python. Here we also have two
different setups. In the first one, we generated $k=100$ points, computed
pairwise Euclidean distances and applied the procedure with parameters $n=3$
and $\varepsilon=3$. The results are shown in Figure 7.3. After that we
computed pairwise Euclidean distances between the recovered points. The
resulting average relative error (7.1) was equal to $Err=0.4$. In the second
setup, we computed exact geodesic distances between points and applied the
procedure with parameters $n=2$ and $\varepsilon=3$. After that, we computed
the Euclidean distances between the embedded points. The resulting average
relative error (7.1) was equal to $Err=0.28$.

Figure 7.3. Reconstruction of the Swiss Roll from pairwise distances. Top
line: sample points in $\mathbb{R}^{3}$ (left), sample points, projection on
the XOZ plane (center), reconstructed Swiss Roll (right). Bottom line: sample
points in $\mathbb{R}^{3}$ (left), sample points, projection on the XOZ plane
(center), embedding into $\mathbb{R}^{2}$ (right).

##  Appendix A Auxiliary lemmata on embedded submanifolds of positive reach

Throughout this section $M\subset\mathbb{R}^{n}$ is a smooth connected
embedded compact Riemannian manifold, with $\mathrm{Reach}\,(M)=\alpha>0$, the
intrinsic geodesic distance on $M$ being denoted by $d_{M}$.

######  Lemma A.1.

Let $M$ be as above and let $u,v$ be any two points on $M$, such that
$d_{M}(u,v)<\alpha$. Then

$0\leq\dfrac{d_{M}(u,v)}{|u-v|}-1\leq C_{1}|u-v|$

for all $\\{u,v\\}\subset M$ with $C_{1}:=\left(3\alpha\sqrt{3}\right)^{-1}$.

###### Proof.

Since $|u-v|\leq d_{M}(u,v)<\alpha$, then by theorem 2.1 of [4] we have

$\frac{d_{M}(u,v)}{2\alpha}\leq\arcsin\frac{|u-v|}{2\alpha}.$

Let $f(t):=\arcsin t$, $t<0.5$. Since $f(0)=0$, $f^{\prime}(0)=1$ and

$f^{\prime\prime}(t)=\frac{x}{(1-x^{2})^{3/2}}\leq\frac{4}{3\sqrt{3}},$

we have

$f(t)\leq t+\frac{2t^{2}}{3\sqrt{3}}.$

Thus,

$\frac{d_{M}(u,v)-|u-v|}{2\alpha}\leq\frac{|u-v|^{2}}{6\alpha^{2}\sqrt{3}},$

which yields

$\frac{d_{M}(u,v)}{|u-v|}-1\leq\frac{|u-v|}{3\alpha\sqrt{3}}$

as claimed. ∎

######  Lemma A.2.

Let $M$ be as above. Then

$\dfrac{|u-v|}{d_{M}(u,v)}\geq C_{2}$

for all $\\{u,v\\}\subset M$ and for some $C_{2}>0$. In particular, one can
take $C_{2}:=(2/\pi)\wedge(2\alpha/D)$, where $D$ stands for the intrinsic
diameter of $M$.

###### Proof.

If $|u-v|<2\alpha$, then by theorem 2.1 from [4], one has

$d_{M}(u,v)\leq 2\alpha\arcsin\dfrac{|u-v|}{2\alpha},$

or, equivalently,

$\dfrac{|u-v|}{d_{M}(u,v)}\geq\dfrac{\sin t}{t},\quad\mbox{where
}t:=\arcsin\dfrac{|u-v|}{2\alpha}.$

Since $0\leq t<\pi/2$, then

$\dfrac{|u-v|}{d_{M}(u,v)}\geq\dfrac{2}{\pi}.$

In the remaining case $|u-v|\geq 2\alpha$, one has

$\dfrac{|u-v|}{d_{M}(u,v)}\geq\dfrac{2\alpha}{D},$

since $d_{M}(u,v)\leq D$. ∎

## References

  * [1] Henry Adams, Mark Blumstein, and Lara Kassab.  Multidimensional scaling on metric measure spaces.  Rocky Mountain Journal of Mathematics, 50(2):397–413, 2020.
  * [2] Ery Arias-Castro, Adel Javanmard, and Bruno Pelletier.  Perturbation bounds for procrustes, classical scaling, and trilateration, with applications to manifold learning.  Journal of Machine Learning Research, 21(15):1–37, 2020.
  * [3] Sanjeev Arora, Wei Hu, and Pravesh K. Kothari.  An analysis of the t-SNE algorithm for data visualization, 2018.
  * [4] Jean-Daniel Boissonnat, André Lieutier, and Mathijs Wintraecken.  The reach, metric distortion, geodesic convexity and the variation of tangent spaces.  J. Appl. Comput. Topol., 3(1-2):29–58, 2019.
  * [5] Herbert Federer.  Curvature measures.  Transactions of the American Mathematical Society, 93(3):418–491, 1959.
  * [6] Charles Fefferman, Sergei Ivanov, Matti Lassas, and Hariharan Narayanan.  Reconstruction of a Riemannian manifold from noisy intrinsic distances.  arXiv preprint arXiv:1905.07182, 2019.
  * [7] Geoffrey Hinton and Sam Roweis.  Stochastic neighbor embedding.  In Proceedings of the 15th International Conference on Neural Information Processing Systems, NIPS’02, page 857–864, Cambridge, MA, USA, 2002\. MIT Press.
  * [8] Partha Niyogi, Stephen Smale, and Shmuel Weinberger.  Finding the homology of submanifolds with high confidence from random samples.  Discrete Comput. Geom., 39(1-3):419–441, 2008.
  * [9] Laurens van der Maaten and Geoffrey Hinton.  Visualizing data using t-SNE.  Journal of Machine Learning Research, 9:2579–2605, 2008.
  * [10] Jianzhong Wang.  Geometric structure of high-dimensional data and dimensionality reduction.  Springer, 2012.