# Collective Variables for the Study of Crystallization

Tarak Karmakar1,2,4 Michele Invernizzi1,3,4 Valerio Rizzi1,2,4 Michele
Parrinello1,2,4 1\. Institute of Computational Sciences, Faculty of
Informatics, Università della Svizzera italiana, Via Giuseppe Buffi 13, 6900
Lugano, Switzerland
2\. Department of Chemistry and Applied Biosciences, ETH Zurich, 8092 Zurich,
Switzerland
3\. Department of Physics, ETH Zurich, 8092 Zurich, Switzerland
4\. Italian Institute of Technology, Via Morego 30, 16163 Genova, Italy
<EMAIL_ADDRESS>
###### Abstract

The phenomenon of solidification of a substance from its liquid phase is of
the greatest practical and theoretical importance, and atomistic simulations
can provide precious information towards its understanding and control.
Unfortunately, the time scale for crystallization is much larger than what can
be explored in standard simulations. Enhanced sampling methods can overcome
this time scale hurdle. Here we employ the on-the-fly probability enhanced
sampling method that is a recent evolution of metadynamics. This method, like
many others, relies on the definition of appropriate collective variables able
to capture the slow degrees of freedom. To this effect we introduce collective
coordinates of general applicability to crystallization simulations. They are
based on the peaks of the three-dimensional structure factor that are combined
non-linearly via the Deep Linear Discriminant Analysis machine learning
method. We apply the method to the study of crystallization of a
multicomponent system, Sodium Chloride and a molecular system, Carbon Dioxide.

###### keywords:

American Chemical Society, LaTeX

+41 58 666 48 01

##  1 1\. Introduction

Understanding at the microscopic level the process of crystallization is of
great importance in many areas of physics, chemistry, and biology.1, 2, 3, 4,
5 From a more general point of view, crystallization provides a spectacular
example of symmetry-breaking first-order phase transition. In principle,
atomistic molecular dynamics simulations could offer invaluable insight, given
the fact that experimental studies are rather challenging. This possibility
did not escape the attention of the early simulators, and already in 1976,
Rahman studied the homogeneous crystallization of a Lennard-Jones liquid.6
However, since crystallization takes place in a macroscopic time scale, Rahman
was forced to apply deep and unrealistically fast temperature quenchings in
order to induce crystallization.6, 7 It became then clear that in order to
study crystallization, enhanced sampling methods are mandatory.8, 9, 10, 11,
12, 13, 14 Such methods allow overcoming the large free energy barriers
involved in the crystallization process, thus circumventing the presence of
kinetic bottlenecks.

The overwhelming majority of such studies have used methods that directly or
indirectly require the definition of a small set of collective variables
(CVs).11, 12, 13, 15, 16, 17, 18, 19, 20, 21 The CVs are meant to encode the
slow modes of the system such that, by applying a bias potential that is a
function of the CVs, it is possible to remove the kinetic bottlenecks.
Designing such CVs is often a non-trivial task and a successful choice of CVs
demands physical intuition and an in-depth understanding of the system
studied. Practical considerations suggest limiting the number of CVs used, and
much ingenuity has been invested into condensing several physical parameters
into a small number of CVs.22, 12, 23

The CVs can be categorized into three classes. In the first class, one could
include generic CVs such as enthalpy and an entropy surrogate whose aim is to
drive crystallization without prior knowledge of the crystal structure.15, 16
These CVs are useful when exploring polymorphism, but they are not very
efficient when it comes to calculating free energies. If the crystal structure
is known, a second class of CVs has been designed that describes the local
crystalline environment.24, 25, 22, 12, 26, 20 These CVs have been very useful
but it is difficult to extend their application to multicomponent systems, and
if not properly chosen, they can drive the system towards unwanted defective
structures. A third class of CVs, that also involves knowledge of the final
crystal structure, is obtained by analyzing the atomic density in the
reciprocal Fourier space and using as CVs the Bragg peak intensities that are
the experimental signature of crystallinity. Guided by this consideration, we
have recently used the peaks in the shperically-averaged structure factor as
possible CVs for crystallization.18, 19, 27 When applying this idea in the
practice, one usually needs to reduce the number of peaks considered e.g., by
linearly combining different peaks with a lengthy trial and error procedure,
as done in Refs. 23 and 28.

In this paper, instead of using the peaks of the spherically-averaged
structure factor following the Debye scattering equation as in Ref. 18, we
propose to use the full three-dimensional structure factor (SF). This has two
important advantages. On one hand, it provides a more accurate description by
resolving each scattering vector into its three-dimensional components. This
avoids the degeneracy intrinsic to the Debye peaks, which can lead to
mislabeling defected structures. On the other hand, the SF is also
computationally more efficient, since it scales linearly with the number of
atoms. However, an inevitable hurdle of this choice is that the 3D SF presents
a number of crystal peaks much greater than its spherically-averaged
counterpart. Thus it becomes crucial to avoid the lengthy procedure of Ref. 28
and replace it by a modern machine learning (ML) technique, in which several
descriptors are non-linearly combined to produce efficient CVs.29, 30 To this
effect, we use the Deep Linear Discrimination Analysis (Deep-LDA) method30 and
obtain a Deep-LDA CV that we then use in the On-the-fly Probability Enhanced
Sampling (OPES) method31, 32 to simulate the crystallization of multicomponent
and molecular systems. As prototypical examples, we have applied with success
our approach to study the crystallization of NaCl and CO2 from their liquid
phase. We have chosen these systems since they present a number of challenges
as discussed in the literature.11, 33, 34, 35

###  1.1 2\. OPES

In order to accelerate sampling, we use the recently developed OPES method.31,
32 In OPES, like in metadynamics36, 37, 38, 39, the bias is constructed on-
the-fly with the help of periodically added Gaussians $G(s,s_{n})$ that are
functions of the chosen CVs $s$ and are centered at the CVs current $s_{n}$
value, where $n$ labels the $n^{th}$ deposition. In metadynamics, the
Gaussians are directly added to the bias potential $V(s)$. In OPES instead,
the added Gaussians modify the equilibrium probability distribution $P(s)$
which in turn is related to the bias by:

$V(s)=-\frac{1}{\beta}\log\frac{p^{\text{tg}}(s)}{P(s)}$ (1)

in which $\beta^{-1}={k_{B}T}$ is the inverse temperature, and
$p^{\text{tg}}(s)$ is the target distribution we want the OPES simulation to
converge to. In the present application, we will aim at obtaining the so-
called well-tempered distribution
$p^{\text{tg}}(s)\propto[P(s)]^{\frac{1}{\gamma}}$, where $\gamma>1$ is called
bias factor. The merits of this choice have been amply discussed in the
metadynamics literature.37, 40, 41

In order to compute $P(s)$, one sets up an iterative scheme in which at the
$n^{th}$ step the probability distribution is estimated as:

$P_{n}(s)=\frac{\sum_{k}^{n}w_{k}G(s,s_{k})}{\sum_{k}^{n}w_{k}}$ (2)

and the weights $w_{k}=e^{\beta V_{k-1}(s_{k})}$ are given by the bias
potential previously deposited.

The relation between the $n^{th}$ estimate of $P_{n}(s)$ and that of
$V_{n}(s)$ is given by

$V_{n}(s)=(1-1/\gamma)\frac{1}{\beta}\log\left(\frac{P_{n}(s)}{Z_{n}}+\epsilon\right)$
(3)

where $Z_{n}$ is a normalization factor, and $\epsilon$ is a regularization
parameter that sets a maximum limit to the bias potential, in order to avoid
exploring too high free energy regions.

The $n^{th}$ free energy estimate then becomes

$F_{n}(s)=-\frac{1}{\beta}\log{P_{n}(s)}$ (4)

Similarly, any other ensemble average of an observable $O=O(\mathbf{R})$ that
is function of the coordinates $\mathbf{R}$, can be estimated via reweighting

$\langle O\rangle=\frac{\langle O\,e^{\beta V}\rangle_{V}}{\langle e^{\beta
V}\rangle_{V}}$ (5)

where $\langle\cdot\rangle_{V}$ denotes the ensemble average over the biased
simulation.

The discussion presented here is highly schematic and misses many of the
important technical details that make OPES efficient. These can be found in
Refs 31 and 32.

##  2 3\. LDA and Deep-LDA

In order to design efficient CVs, we follow the approach of Ref. 30 that
relies on the linear discriminant analysis (LDA) introduced a long time ago by
Fisher42. We specialize our discussion to the case in which only two
metastable states are present since here we are only interested in studying
transitions between solid (S) and liquid (L) phases. For a treatment of the
more general case in which several metastable states are present the
interested reader is referred to the literature.30

The transition between solid and liquid is a rare event occurring on a time
scale that far exceeds that of standard simulations. Thus, if we perform two
distinct simulations one in the solid and the other in the liquid, no
transition will be observed, and the configurations thus generated will be
separated in the atomic coordinate multidimensional space. We want to project
these data into a low dimensional set of CV that is still able to discriminate
between the two states. As stated earlier, here we shall apply a method
derived from LDA.

In LDA, one starts by selecting an exhaustive set of descriptors that are
$N_{d}$ functions of the atomic coordinates $\boldsymbol{R}$. For notational
convenience we arrange the $d_{i}(\boldsymbol{R})$ to form a vector
$\boldsymbol{d(R)}$, and in the following, we shall not make explicit the
$\boldsymbol{R}$ dependence of $\boldsymbol{d(R)}$. The descriptors could in
theory be all the $\boldsymbol{R}$, but it is better to preprocess the
$\boldsymbol{R}$ such that the result remains physically transparent, and the
natural symmetries of the system are respected. For instance, one could use as
descriptors distances or angles that have the property of being
translationally and rotationally invariant. In addition, for the purpose of
generating useful CVs, the descriptors should be able to correctly distinguish
between the various metastable and transition states. This requirement is not
necessarily satisfied by a CV that is generated in a discrimination procedure.
However, experience has shown that this is to a good extent the case.

The object of LDA is to find a linear combination of descriptors
$s=\mathbf{w^{T}d}$ such that when the data are projected on $s$ the variance
of S and L data is large while the total spread is low. In order to obtain the
optimal $\mathbf{w}$, one first runs two distinct simulations, one in the
solid and the other in the liquid. From the data thus collected one computes
the averages and covariances of the descriptors in each state and arrange them
to form the vectors $\boldsymbol{d_{S}}$ and $\boldsymbol{d_{L}}$ and matrices
$\boldsymbol{\Sigma_{S}}$ and $\boldsymbol{\Sigma_{L}}$, respectively. These
quantities are combined to form two $N_{d}\times N_{d}$ covariance matrices,
the within
$\boldsymbol{S}_{w}=\frac{1}{2}\left(\boldsymbol{\Sigma_{S}}+\boldsymbol{\Sigma_{L}}\right)$
and the between
$\boldsymbol{S}_{b}=\boldsymbol{(d_{S}-d_{L})(d_{S}-d_{L})^{T}}$ one. The
former measures the total spread of the projected distribution while the
latter informs on its total variance.

The condition of high separation in the $s$ projection between data pertaining
to the different states and low variance within each state is obtained by
maximizing the Fisher ratio:

$\displaystyle\frac{\boldsymbol{\mathbf{w}S_{b}\mathbf{w}^{T}}}{\boldsymbol{\mathbf{w}S_{w}\mathbf{w}^{T}}}$
(6)

with respect to $\mathbf{w}$. The value of $\mathbf{w}^{*}$ that maximizes the
Fisher ratio can be obtained as the eigenvector of the highest eigenvalue of
the generalized eigenvalue equation.

$\displaystyle\boldsymbol{S_{b}\mathbf{w_{i}}}=\nu_{i}\boldsymbol{S_{w}\mathbf{w_{i}}}$
(7)

The sought-after CV is then

$\displaystyle s=\boldsymbol{\mathbf{w}^{*T}d}$ (8)

In Deep-LDA, the $\boldsymbol{d}$ is first inputted into a feed-forward deep
neuronal network (NN), whose output $\boldsymbol{h}$ has dimension $N_{h}$.
One then starts with a guess for $\boldsymbol{h}$, calculates
$\boldsymbol{S_{b}}$ and $\boldsymbol{S_{w}}$ that, in the $\boldsymbol{h}$
basis, are now $N_{h}\times N_{h}$ matrices. The NN weights are optimized
using as loss function the eigenvector with lowest eigenvalue ($\nu$) of the
Fisher’s generalized equation Eq. 7 that is now a $N_{h}\times N_{h}$
equation, this leads to a new $\boldsymbol{h}$ and to a new optimization
procedure. The cycle is repeated until convergence leading to the Deep-LDA CV,

$\displaystyle s=\boldsymbol{\mathbf{w}^{T}h}$ (9)

This non-linear procedure has proven to lead to more efficient CVs than the
one based on LDA even when using its harmonic version HLDA.43

##  3 4\. Collective Variable

As discussed in the previous section we want to combine a set of descriptors
in a non-linear way so as to obtain efficient CVs for the study of
crystallization. Typically, crystal order parameters are defined locally, by
considering the environment of each atom up to a given cutoff radius. We
introduce instead a descriptor capable to take into account for the long-range
order of the crystal. Following Ref. 18 we propose to use the structure factor
as descriptor, but contrary to Ref. 18, here we will not use the isotropic
Debye formula, but rather the full three-dimensional structure factor, as
defined for a monoatomic substance

$S({\boldsymbol{k}})=\frac{1}{N}\left\lvert\sum_{i=1}^{N}e^{-i\boldsymbol{k}\cdot\boldsymbol{R_{i}}}\right\rvert^{2}$
(10)

where $\mathbf{k}$ is the three-dimensional scattering vector. For simplicity,
we consider a periodically repeating cubic simulation box of edge $L$, so
$\mathbf{k}$ can be written as

$\boldsymbol{k}=\frac{2\pi}{L}(l,m,n)$ (11)

with $l,m,n$ integer indexes. In the case of a crystal structure that is
commensurate to the MD cell, only a subset of $l,m,n$ will correspond to a
Bragg peak. Our choice is to set ourselves in this condition and choose also
the number of atoms such that perfect commensuration is possible. Out of
simplicity, we are considering here only the case of a cubic box, but
extension to more general type of periodic boundary conditions is
straightforward.

The advantage of Eq. 10 is that, contrary to the Debye formula, it scales
linearly with the number of atoms $N$ without the need to introduce a cutoff
radius. The main drawback is that the number of valid scattering vectors
$\mathbf{k}$ grows roughly as $N^{3}$. This however, does not constitute a
problem, because for any given crystal structure $S(\mathbf{k})$ will present
only few high intensity peaks at specific values of $\mathbf{k}$ while all the
others can be safely ignored. Notice that the peaks of the three-dimensional
SF should not be confused with the peaks of the spherically-averaged SF,
$S(|\mathbf{k}|)$, that is the one typically measured in powder X-ray
diffraction experiments.

In order to define the SF descriptors for a given system, we proceed as
follows. First we run a short unbiased simulation of the perfect crystal and
calculate the average $S(\mathbf{k})$. We then consider as descriptors only
those SF peaks that are greater than a threshold value $\varepsilon_{I}$ and
whose scattering vector is smaller in modulus than a preassigned cutoff
$|\mathbf{k}|<k_{c}$ where information on the X-ray crystal structure are
damped by the Debye-Waller factor. This typically still leaves us with a
sizable number of descriptors, thus they cannot be used directly as CVs, but
they can be combined via the Deep-LDA machine learning algorithm, in order to
obtain a single CV, as shown schematically in Fig. 1.

For practical reasons, before feeding them to the neural network, we take the
square root of the structure factor peaks, thus the descriptors are the
following:

$d_{\mathbf{k}}=\sqrt{S(\mathbf{k})}=\frac{1}{\sqrt{N}}\left\lvert\sum_{i=1}^{N}e^{-i\boldsymbol{k\cdot
R_{i}}}\right\rvert$ (12)

Figure 1: Schematic of the Deep-LDA protocol, left: A representative
spherically-averaged structure factor $S(|\boldsymbol{k}|)$ profile. Each of
the shown peaks is constituted of several SF peaks, corresponding to
scattering vectors $\mathbf{k}$ that have same modulus, but different
components. The corresponding descriptors ($\boldsymbol{d_{k}}$) are provided
as input for the Deep-LDA network.

An important characteristic of the proposed descriptors is that by selecting
only the $\mathbf{k}$ vectors corresponding to a given crystal configuration,
we obtain descriptors that are not rotationally invariant, but that favour the
formation of a crystal perfectly aligned with the simulation box. This, as
proposed in Ref. 26, is essential to avoid secondary defected states and
greatly simplify and speedup the calculation.

##  4 5\. Simulation details

The crystal structures of NaCl and CO2 were taken from the crystallographic
database. The unit cells were replicated to generate the supercells using the
AVOGADRO software.44 In the subsequent step, the supercells were minimized and
thermally equilibrated in constant temperature simulation. The corresponding
liquid phases were obtained by melting the crystal at a temperature higher
than their respective melting temperatures and afterward cooling down to the
desired temperature.

The systems were further equilibrated by carrying out short simulations in the
isothermal-isobaric ($NPT$) ensemble. The temperature and pressure of the
system were controlled by the stochastic velocity rescaling thermostat45 and
the Berendsen barostat,46 respectively. For the production $NPT$ simulations,
we retain the same thermostat but switch to an isotropic Parrinello-Rahman
barostat.47 A time step of 2 fs is used in all production simulations. The
coupling constants for the thermostat and barostat were 0.5 ps and 1 ps,
respectively. A cut-off of 0.7 and 0.85 was used for both the van der Waals
and the short-range Coulomb interactions for NaCl and CO2 systems,
respectively. Particle Mesh Ewald was utilized to treat the long-range
electrostatic interactions. The LINCS algorithm was used to constraint the
bonds. The Joung-Cheatham48 and EPM2-flexible potential model49, 50 parameters
were used to simulate NaCl and CO2, respectively.

The OPES bias uses adaptive kernel bandwidth and is updated every 1000 time
steps in the case of NaCl and every 500 for CO2, while the barrier parameter
is set to 600 and 200 kJ/mol respectively. All simulations were carried out in
GROMACS-2019.4 software51 patched with a private version of the PLUMED2.7-dev
plugin.52 The in-house built codes and input files that were used in this work
are available at PLUMED-NEST (…), a public repository of the PLUMED
consortium.53 The figures were prepared using the Visual Molecular Dynamics
(VMD) software.54

##  5 6\. Results and discussion

###  5.1 6.1. Sodium Chloride

While much attention has been paid to the nucleation of monatomic systems,
multicomponent systems have been much less studied. As a representative
example, we discuss here the case of NaCl.

Following the scheme described earlier we perform an equilibrium simulation of
NaCl in its rock salt structure at T=1250 K. For graphical convenience, we
report in Fig. 2(a) the spherically-averaged SF profile,
$S(|\boldsymbol{k}|)$, as obtained from this simulation. . We use the cutoff
values (See section Collective variable) $k_{c}$ = 3.4 and $\varepsilon_{I}$ =
7 to select the $\boldsymbol{k}$ whose intensity will form our set of
descriptors. The chosen $\boldsymbol{k}$’s correspond to the Miller indices
(111), (200), (220) and (311). Considering that we introduce for each
$\boldsymbol{k}$ orientation a different descriptor, this leads to 15
independent $d_{\mathbf{k}}$ (note that $\boldsymbol{k}$ and $-\boldsymbol{k}$
are equivalent). After having performed an unbiased simulations in the solid
(S) and liquid (L), we can compute the averages and fluctuations of the
$d_{\mathbf{k}}$’s in the S and L states and build from those the Deep-LDA CV
as prescribed in Sec. 3.

Figure 2: (a) The $S(|\mathbf{k}|)$ profiles for Na+ and Cl- ions. The two
profiles are identical, so for visibility the Na+ profile is shifted upward.
(b) The time series of NaCl CV $s$ obtained from the OPES simulation. The
color is based on the first descriptor that corresponds to the Miller index
(111) used in the Deep-LDA CV, red: crystal and blue: liquid. and (c) The free
energy profile as a function of $s$ for the crystallization process at 1250 K,
two representative configurations, (d) rock-salt crystal and the liquid
configurations are shown.

The fact that NaCl is a two-component system requires extending our scheme to
a multi-component system. In a strong ionic system like NaCl in which
establishing local neutrality is paramount, the two components Na+ and Cl-
will have a similar behavior. Thus, for each ionic species we determine its CV
$s_{Na}$ and $s_{Cl}$ as described previously. Finally, the CV used in the
OPES simulation is the linear combination, $s=0.5*(s_{Na}+s_{Cl})$. We perform
an OPES simulation with this CV at T 1250 K and P 1 bar.

A large number of reversible transitions between the solid and the liquid
phase manifests the effectiveness of the CV, Fig. 2(b), and allows us to
calculate the converged free energy profile for the crystallization process.
The free energy profile is shown in Fig. 2(c), and the free energy difference
($\Delta F$) between the two states converges to $\sim$10 kJ/mol in $\sim$ 160
ns simulation time (see Fig. S2). At T = 1250 K, both the crystal and liquid
basins have almost comparable free energies with a slightly deeper minimum for
the crystal basin. This agrees with the fact that the chosen temperature, T =
1250 K is a little below the melting temperature of the the Joung-Cheatham
potential model of NaCl which was found to be 1285 $\pm$ 5 K.48, 55

###  5.2 6.2. Carbon dioxide

We extend our approach to study the crystallization of molecular systems and
choose for this study, carbon dioxide. This seemingly simple system has a very
rich phase diagram consisting of gas, compressed fluid, and several
polymorphic crystal phases.56, 57, 58, 59, 60, 61, 62 At ambient conditions,
CO2 exists as a gas, and when cooled and compressed, it transforms into a
dense liquid which can further be converted to a crystal, popularly known as
dry ice. In this work, we restrict ourselves to studying the phase transition
from the liquid to its crystal polymorph-I (Fig. 3). We choose to simulate a
3$\times$3$\times$3 crystal super-cell containing 108 CO2 molecules at T = 400
K and P = 2 GPa. These T and P values are close to the liquid-solid phase
coexistence point.

Figure 3: (a) The $S(|\mathbf{k}|)$ profile of CO2 crystal. (b) Time evolution
of the CO2 CV $s$ obtained from the OPES simulation. The color is based on the
value of the first descriptor that corresponds to the Miller index (111) used
in the Deep-LDA CV, red: crystal and blue: liquid. (c) The free energy profile
for the crystallization process at T = 400 K and P = 2 GPa, and representative
(d) crystal and liquid configurations are shown.

In the CO2 structure, it suffices to specify the oxygen positions to determine
the full crystal structure. For this reason we build the CO2 CV from the
oxygen partial densities. Thus we calculate the $S_{OO}(\boldsymbol{k}$)
partial structure factor with respect to the oxygen atoms (Fig. 3). Here we
choose the peaks having $\varepsilon_{I}>6$ with Miller indices (111), (012),
(121), (302), (132), and (004) (Fig. 3(a)) are fed to the Deep-LDA network to
obtain the CV. Subsequently, we have carried out an OPES simulation with this
CV at T = 400 K and P = 2 GPa.

The system visits the S and L states several times, as can be seen in Fig.
3(b) and $\Delta F$ is converged within 300 ns simulation time (Fig. S4). The
free energy profile is shown in Fig. 3(c). A nucleation barrier of $\sim$140
kJ/mol and a free energy difference of $\sim$5 kJ/mol between the two states
is obtained. In agreement with the experiments, at the thermodynamic
conditions studied here, the liquid state is slightly more stable than the
crystal.

##  6 7\. Conclusions

In this paper, we introduce the use of the structure factor peaks as
descriptors for enhanced sampling simulations of crystallization events.
Contrary to the CVs previously proposed in Ref. 18, based on the Debye formula
for spherically-averaged SF, our descriptors are non-degenerate and scale more
efficiently with the system size. They also induce crystallization only along
the box directions, similarly to the CVs proposed in Ref. 26. These SF
descriptors, combined with the Deep-LDA ML method, provide a systematic and
efficient way for obtaining low dimensional CVs for a vast range of
crystallization problems.

The effectiveness of our approach in simulation of crystallization is
manifested by two representative examples, NaCl and CO2. In each case, the
Deep-LDA CV is found effective in reversibly transforming the system from its
liquid to the crystal state in relatively short simulation time. This allows
us to calculate the free energy of the liquid and crystal phase at the chosen
temperature and pressure. Our method opens up new avenues to efficiently study
the liquid-solid phase transition and estimate the stability of the crystal
phase relative to its molten state. We believe that this strategy can be
applied to study other phase transitions such solid-solid, amorphous-solid,
and liquid-liquid transitions.

The research was supported by the NCCR MARVEL (D&D1) funded by Swiss National
Science Foundation and the European Union Grant No.
ERC-2014-AdG-670227/VARMET. We thank CSCS, the Swiss National Supercomputing
Centre for providing the computational resources.

The Supporting Information (SI) contains free energy difference plots (Fig. S1
and S2) obtained from OPES simulations of NaCl and CO2 systems.

##  7 Code and data Availability

All input files, codes, scripts, and instructions to reproduce the results
presented in this manuscript can be found at the PLUMED-NEST repository as
plumID:XXX. Additionally, all simulations data are available from the authors
upon request.

## References

  * Hulliger 1994 Hulliger, J. Chemistry and crystal growth. _Angewandte Chemie International Edition in English_ 1994, _33_ , 143–162
  * Schüler and Frankel 1999 Schüler, D.; Frankel, R. B. Bacterial magnetosomes: microbiology, biomineralization and biotechnological applications. _Applied Microbiology and Biotechnology_ 1999, _52_ , 464–473
  * Sloan 2003 Sloan, E. D. Fundamental principles and applications of natural gas hydrates. _Nature_ 2003, _426_ , 353–359
  * Doherty 2006 Doherty, M. F. Crystal engineering: from molecules to products. _Chemical Engineering Education_ 2006, _40_ , 116–125
  * Erdemir et al. 2007 Erdemir, D.; Lee, A. Y.; Myerson, A. S. Polymorph selection: the role of nucleation, crystal growth and molecular modeling. _Current opinion in drug discovery & development_ 2007, _10_ , 746–755
  * Mandell et al. 1976 Mandell, M.; McTague, J.; Rahman, A. Crystal nucleation in a three-dimensional Lennard-Jones system: A molecular dynamics study. _The journal of chemical physics_ 1976, _64_ , 3699–3702
  * Mandell et al. 1977 Mandell, M.; McTague, J.; Rahman, A. Crystal nucleation in a three-dimensional Lennard-Jones system. II. Nucleation kinetics for 256 and 500 particles. _The Journal of Chemical Physics_ 1977, _66_ , 3070–3075
  * Auer and Frenkel 2001 Auer, S.; Frenkel, D. Prediction of absolute crystal-nucleation rate in hard-sphere colloids. _Nature_ 2001, _409_ , 1020–1023
  * ten Wolde and Frenkel 1997 ten Wolde, P. R.; Frenkel, D. Enhancement of protein crystal nucleation by critical density fluctuations. _Science_ 1997, _277_ , 1975–1978
  * Tribello et al. 2009 Tribello, G. A.; Bruneval, F.; Liew, C.; Parrinello, M. A molecular dynamics study of the early stages of calcium carbonate growth. _The Journal of Physical Chemistry B_ 2009, _113_ , 11680–11687
  * Giberti et al. 2013 Giberti, F.; Tribello, G. A.; Parrinello, M. Transient polymorphism in NaCl. _Journal of chemical theory and computation_ 2013, _9_ , 2526–2530
  * Giberti et al. 2015 Giberti, F.; Salvalaglio, M.; Mazzotti, M.; Parrinello, M. Insight into the nucleation of urea crystals from the melt. _Chemical Engineering Science_ 2015, _121_ , 51–59
  * Salvalaglio et al. 2013 Salvalaglio, M.; Vetter, T.; Mazzotti, M.; Parrinello, M. Controlling and predicting crystal shapes: The case of urea. _Angewandte Chemie_ 2013, _125_ , 13611–13614
  * Giberti et al. 2015 Giberti, F.; Salvalaglio, M.; Parrinello, M. Metadynamics studies of crystal nucleation. _IUCrJ_ 2015, _2_ , 256–266
  * Piaggi et al. 2017 Piaggi, P. M.; Valsson, O.; Parrinello, M. Enhancing entropy and enthalpy fluctuations to drive crystallization in atomistic simulations. _Physical Review Letters_ 2017, _119_ , 015701
  * Piaggi and Parrinello 2018 Piaggi, P. M.; Parrinello, M. Predicting polymorphism in molecular crystals using orientational entropy. _Proceedings of the National Academy of Sciences_ 2018, _115_ , 10251–10256
  * Gobbo et al. 2018 Gobbo, G.; Bellucci, M. A.; Tribello, G. A.; Ciccotti, G.; Trout, B. L. Nucleation of molecular crystals driven by relative information entropy. _Journal of Chemical Theory and Computation_ 2018, _14_ , 959–972
  * Niu et al. 2018 Niu, H.; Piaggi, P. M.; Invernizzi, M.; Parrinello, M. Molecular dynamics simulations of liquid silica crystallization. _Proceedings of the National Academy of Sciences_ 2018, _115_ , 5348–5352
  * Bonati and Parrinello 2018 Bonati, L.; Parrinello, M. Silicon liquid structure and crystal nucleation from ab initio deep metadynamics. _Physical review letters_ 2018, _121_ , 265701
  * Piaggi and Car 2020 Piaggi, P. M.; Car, R. Phase equilibrium of liquid water and hexagonal ice from enhanced sampling molecular dynamics simulations. _The Journal of Chemical Physics_ 2020, _152_ , 204116
  * Niu et al. 2020 Niu, H.; Bonati, L.; Piaggi, P. M.; Parrinello, M. Ab initio phase diagram and nucleation of gallium. _Nature Communications_ 2020, _11_ , 1–9
  * Santiso and Trout 2011 Santiso, E. E.; Trout, B. L. A general set of order parameters for molecular crystals. _The Journal of chemical physics_ 2011, _134_ , 064109
  * Zhang et al. 2019 Zhang, Y.-Y.; Niu, H.; Piccini, G.; Mendels, D.; Parrinello, M. Improving collective variables: The case of crystallization. _The Journal of chemical physics_ 2019, _150_ , 094509
  * Rein ten Wolde et al. 1996 Rein ten Wolde, P.; Ruiz-Montero, M. J.; Frenkel, D. Numerical calculation of the rate of crystal nucleation in a Lennard-Jones system at moderate undercooling. _The Journal of chemical physics_ 1996, _104_ , 9932–9947
  * Lechner and Dellago 2008 Lechner, W.; Dellago, C. Accurate determination of crystal structures based on averaged local bond order parameters. _The Journal of chemical physics_ 2008, _129_ , 114707
  * Piaggi and Parrinello 2019 Piaggi, P. M.; Parrinello, M. Calculation of phase diagrams in the multithermal-multibaric ensemble. _The Journal of chemical physics_ 2019, _150_ , 244119
  * Invernizzi and Parrinello 2019 Invernizzi, M.; Parrinello, M. Making the Best of a Bad Situation: A Multiscale Approach to Free Energy Calculation. _Journal of Chemical Theory and Computation_ 2019, _15_ , 2187–2194
  * Niu et al. 2019 Niu, H.; Yang, Y. I.; Parrinello, M. Temperature dependence of homogeneous nucleation in ice. _Physical review letters_ 2019, _122_ , 245501
  * Rogal et al. 2019 Rogal, J.; Schneider, E.; Tuckerman, M. E. Neural-Network-Based Path Collective Variables for Enhanced Sampling of Phase Transformations. _Physical Review Letters_ 2019, _123_ , 245701
  * Bonati et al. 2020 Bonati, L.; Rizzi, V.; Parrinello, M. Data-driven collective variables for enhanced sampling. _The Journal of Physical Chemistry Letters_ 2020, _11_ , 2998–3004
  * Invernizzi and Parrinello 2020 Invernizzi, M.; Parrinello, M. Rethinking Metadynamics: from bias potentials to probability distributions. _The Journal of Physical Chemistry Letters_ 2020, _11_ , 2731–2736
  * Invernizzi et al. 2020 Invernizzi, M.; Piaggi, P. M.; Parrinello, M. Unified Approach to Enhanced Sampling. _Physical Review X_ 2020, _10_ , 041034
  * Valeriani et al. 2005 Valeriani, C.; Sanz, E.; Frenkel, D. Rate of homogeneous crystal nucleation in molten NaCl. _The Journal of chemical physics_ 2005, _122_ , 194501
  * Gimondi and Salvalaglio 2017 Gimondi, I.; Salvalaglio, M. CO2 packing polymorphism under pressure: Mechanism and thermodynamics of the I-III polymorphic transition. _The Journal of Chemical Physics_ 2017, _147_ , 114502
  * Gimondi and Salvalaglio 2018 Gimondi, I.; Salvalaglio, M. CO 2 packing polymorphism under confinement in cylindrical nanopores. _Molecular Systems Design & Engineering_ 2018, _3_ , 243–252
  * Laio and Parrinello 2002 Laio, A.; Parrinello, M. Escaping free-energy minima. _Proceedings of the National Academy of Sciences_ 2002, _99_ , 12562–12566
  * Barducci et al. 2008 Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: a smoothly converging and tunable free-energy method. _Physical review letters_ 2008, _100_ , 020603
  * Barducci et al. 2011 Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2011, _1_ , 826–843
  * Valsson et al. 2016 Valsson, O.; Tiwary, P.; Parrinello, M. Enhancing important fluctuations: Rare events and metadynamics from a conceptual viewpoint. _Annual review of physical chemistry_ 2016, _67_ , 159–184
  * Barducci et al. 2010 Barducci, A.; Bonomi, M.; Parrinello, M. Linking well-tempered metadynamics simulations with experiments. _Biophysical journal_ 2010, _98_ , L44–L46
  * Dama et al. 2014 Dama, J. F.; Parrinello, M.; Voth, G. A. Well-tempered metadynamics converges asymptotically. _Physical review letters_ 2014, _112_ , 240602
  * Welling 2005 Welling, M. _Fisher linear discriminant analysis. department of computer science, university of toronto_ ; 2005
  * Mendels et al. 2018 Mendels, D.; Piccini, G.; Parrinello, M. Collective variables from local fluctuations. _The journal of physical chemistry letters_ 2018, _9_ , 2776–2781
  * Hanwell et al. 2012 Hanwell, M. D.; Curtis, D. E.; Lonie, D. C.; Vandermeersch, T.; Zurek, E.; Hutchison, G. R. Avogadro: an advanced semantic chemical editor, visualization, and analysis platform. _Journal of cheminformatics_ 2012, _4_ , 17
  * Bussi et al. 2007 Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. _J. Chem. Phys._ 2007, _126_ , 014101
  * Berendsen et al. 1984 Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. _The Journal of Chemical Physics_ 1984, _81_ , 3684–3690
  * Parrinello and Rahman 1981 Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. _J. Appl. Phys._ 1981, _52_ , 7182–7190
  * Joung and Cheatham III 2008 Joung, I. S.; Cheatham III, T. E. Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations. _The journal of physical chemistry B_ 2008, _112_ , 9020–9041
  * Nieto-Draghi et al. 2007 Nieto-Draghi, C.; de Bruin, T.; Pérez-Pellitero, J.; Bonet Avalos, J.; Mackie, A. D. Thermodynamic and transport properties of carbon dioxide from molecular simulation. _The Journal of Chemical Physics_ 2007, _126_ , 064509
  * Zhang and Duan 2005 Zhang, Z.; Duan, Z. An optimized molecular potential for carbon dioxide. _The Journal of Chemical Physics_ 2005, _122_ , 214507
  * Abraham et al. 2015 Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. _SoftwareX_ 2015, _1_ , 19–25
  * Tribello et al. 2014 Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. _Computer Physics Communications_ 2014, _185_ , 604–613
  * Bonomi et al. 2019 others,, et al. Promoting transparency and reproducibility in enhanced molecular simulations. _Nature methods_ 2019, _16_ , 670–673
  * Humphrey et al. 1996 others,, et al. VMD: visual molecular dynamics. _Journal of molecular graphics_ 1996, _14_ , 33–38
  * Aragones et al. 2012 Aragones, J.; Sanz, E.; Valeriani, C.; Vega, C. Calculation of the melting point of alkali halides by means of computer simulations. _The Journal of chemical physics_ 2012, _137_ , 104507
  * Datchi et al. 2009 Datchi, F.; Giordano, V. M.; Munsch, P.; Saitta, A. M. Structure of carbon dioxide phase IV: Breakdown of the intermediate bonding state scenario. _Physical review letters_ 2009, _103_ , 185701
  * Santoro and Gorelli 2006 Santoro, M.; Gorelli, F. High pressure solid state chemistry of carbon dioxide. _Chemical Society Reviews_ 2006, _35_ , 918–931
  * Yoo 2013 Yoo, C.-S. Physical and chemical transformations of highly compressed carbon dioxide at bond energies. _Physical Chemistry Chemical Physics_ 2013, _15_ , 7949–7966
  * Giordano and Datchi 2007 Giordano, V.; Datchi, F. Molecular carbon dioxide at high pressure and high temperature. _EPL (Europhysics Letters)_ 2007, _77_ , 46002
  * Iota and Yoo 2001 Iota, V.; Yoo, C.-S. Phase diagram of carbon dioxide: Evidence for a new associated phase. _Physical review letters_ 2001, _86_ , 5922
  * Bonev et al. 2003 Bonev, S.; Gygi, F.; Ogitsu, T.; Galli, G. High-pressure molecular phases of solid carbon dioxide. _Physical review letters_ 2003, _91_ , 065501
  * Etters and Kuchta 1989 Etters, R.; Kuchta, B. Static and dynamic properties of solid CO2 at various temperatures and pressures. _The Journal of chemical physics_ 1989, _90_ , 4537–4541


Current address: ]Max Planck Institut für Mikrostrukturphysik, Weinberg 2,
06120 Halle, Germany Current address: ]Department of Physics, Birla Institute
of Technology, Ranchi 835215 Jharkhand, India

Current address: ]Max Planck Institut für Mikrostrukturphysik, Weinberg 2,
06120 Halle, Germany Current address: ] Key Laboratory of Polar Materials and
Devices (MOE), School of Physics and Electronic Science, East China Normal
University, Shanghai 200241, China

# Hydrogen-impurity induced unconventional magnetism in semiconducting
molybdenum ditelluride

Jonas A. Krieger  [  Laboratory for Muon Spin Spectroscopy, Paul Scherrer
Institute, CH-5232 Villigen PSI, Switzerland  Swiss Light Source, Paul
Scherrer Institute, CH-5232 Villigen PSI, Switzerland  Laboratorium für
Festkörperphysik, ETH Zürich, CH-8093 Zürich, Switzerland  Daniel Tay
Laboratorium für Festkörperphysik, ETH Zürich, CH-8093 Zürich, Switzerland
Igor P. Rusinov  Tomsk State University, pr. Lenina 36, 634050 Tomsk, Russia
St. Petersburg State University, Universitetskaya nab. 7/9, 199034 St.
Petersburg, Russia  Sourabh Barua  [  Department of Physics, University of
Warwick, Coventry CV4 7AL, UK  Pabitra K. Biswas  ISIS Facility, Rutherford
Appleton Laboratory, Chilton, Didcot, Oxon OX110QX, United Kingdom  Lukas
Korosec  Laboratorium für Festkörperphysik, ETH Zürich, CH-8093 Zürich,
Switzerland  Thomas Prokscha  Laboratory for Muon Spin Spectroscopy, Paul
Scherrer Institute, CH-5232 Villigen PSI, Switzerland  Thorsten Schmitt  Swiss
Light Source, Paul Scherrer Institute, CH-5232 Villigen PSI, Switzerland
Niels B. M. Schröter  [  Swiss Light Source, Paul Scherrer Institute, CH-5232
Villigen PSI, Switzerland  Tian Shang  [  Laboratory for Multiscale Materials
Experiments, Paul Scherrer Institut,CH-5232 Villigen PSI, Switzerland  Toni
Shiroka  Laboratory for Muon Spin Spectroscopy, Paul Scherrer Institute,
CH-5232 Villigen PSI, Switzerland  Laboratorium für Festkörperphysik, ETH
Zürich, CH-8093 Zürich, Switzerland  Andreas Suter  Laboratory for Muon Spin
Spectroscopy, Paul Scherrer Institute, CH-5232 Villigen PSI, Switzerland
Geetha Balakrishnan  Department of Physics, University of Warwick, Coventry
CV4 7AL, UK  Evgueni V. Chulkov  Donostia International Physics Center, P.
Manuel de Lardizabal 4, San Sebastián, 20018 Basque Country, Spain
Departamento de Física de Materiales UPV/EHU, Centro de Física de Materiales
CFM - MPC and Centro Mixto CSIC-UPV/EHU, 20080 San Sebastián/Donostia, Spain
Tomsk State University, pr. Lenina 36, 634050 Tomsk, Russia  St. Petersburg
State University, Universitetskaya nab. 7/9, 199034 St. Petersburg, Russia
Vladimir N. Strocov<EMAIL_ADDRESS>Swiss Light Source, Paul Scherrer
Institute, CH-5232 Villigen PSI, Switzerland  Zaher Salman
<EMAIL_ADDRESS>Laboratory for Muon Spin Spectroscopy, Paul Scherrer
Institute, CH-5232 Villigen PSI, Switzerland

###### Abstract

Layered transition-metal dichalcogenides are proposed as building blocks for
van der Waals (vdW) heterostructures due to their graphene-like two
dimensional structure. For this purpose, a magnetic semiconductor could
represent an invaluable component for various spintronics and topotronics
devices. Here, we combine different local magnetic probe spectroscopies with
angle-resolved photoemission and density-functional theory calculations to
show that 2H-MoTe2 is on the verge of becoming magnetic. Our results present
clear evidence that the magnetism can be “switched on” by a hydrogen-like
impurity. We also show that this magnetic state survives up to the free
surface region, demonstrating the material’s potential applicability as a
magnetic component for thin-film heterostructures.

Layered transition-metal dichalcogenides (TMDs) are a promising class of
materials featuring interesting optoelectronic, superconducting, catalytic,
and topological properties. They can also be exfoliated down to monolayers,
making them suitable for potential applications in complex van der Waals (vdW)
heterostructures Geim and Grigorieva (2013); Rajamathi et al. (2017). Among
the TMDs, polymorphic MoTe2, exhibits a variety of electronic and magnetic
properties: Its metallic 1T’ phase turns into a Td-phase at
$\approx$250\text{\,}\mathrm{K}$$ Clarke et al. (1978). The latter is a
topological Weyl semimetal that features unconventional superconductivity at
low temperature Sun et al. (2015); Qi et al. (2016); Xu et al. (2018). The
2H-phase, the most stable one at ambient conditions, is an indirect band gap
semiconductor and has been successfully integrated in transistors and photo-
detectors Lezama et al. (2014); Xu et al. (2015); Ding et al. (2018).
Moreover, it has been demonstrated that the 2H-phase can be distorted into the
1T’-phase as a function of strain or gate voltage Song et al. (2016); Wang et
al. (2017). In addition, room temperature (RT) magnetic hysteresis has been
observed for small amounts of vanadium deposited on 2H-MoTe2 Coelho et al.
(2019). However, most interestingly, 2H-MoTe2 was shown by muon spin
spectroscopy ($\mu^{+}$SR) to be a magnetic semiconductor, exhibiting a long-
range magnetic order in the bulk at low temperatures, which was tentatively
attributed to Mo anti-site defects Guguchia et al. (2018). However, no other
technique has reported such long range magnetic order and the origin of this
magnetism remains still unclear. In fact, theoretical calculations also
predict that, apart from anti-site defects, hydrogen and transition-metal
dopants are able to induce spin polarization in MoTe2 in the monolayer limit
Ma et al. (2011); Kanoun (2018). Since a muon can be considered to be a light
hydrogen isotope, one possible explanation is that the magnetism observed in
$\mu^{+}$SR is not intrinsic, but rather induced by the muon itself. Another
important consideration is that in order to exploit these magnetic properties,
e.g. in thin-film heterostructures, one should first elucidate their origin
and evolution as a function of depth in the vicinity of an interface (with
vacuum or other materials).

Here, we address these questions by utilizing complementary local spin probe
measurements, with implanted muons, as well as with the host’s 125Te-NMR to
reveal the origin of the reported magnetism. We find that while some of the
implanted muons, which form a hydrogen-like state in 2H-MoTe2, detect
magnetism in this material, 125Te-NMR shows no evidence of any intrinsic
magnetic order. These findings were then confirmed by detailed density-
functional theory (DFT) calculations, reproducing the experimental 3D
electronic band structure. When including isolated hydrogen (i.e., muon like)
impurities, the DFT unveils an induced magnetic moment at the impurity.
Furthermore, we show that this hydrogen-induced magnetic state remains
unchanged towards the crystal surface, establishing 2H-MoTe2 as a switchable
magnetic semiconductor by hydrogen doping, thus paving the path towards new
possibilities for integrating this intriguing TMD into future applications in
thin-film heterostructures.

In order to gain a better understanding of the nature of the previously
observed magnetic signal in 2H-MoTe2, we combine different local magnetic
spectroscopies, (i) $\mu^{+}$SR which uses individually implanted, spin-
polarized muons, stopping at interstitial sites in the lattice, and (ii)
125Te-NMR on the tellurium atoms, intrinsic to the host compound. This is
illustrated in Fig. 1(a).

Figure 1:  (a) Schematic illustration of the different probing sites. The
muon’s stopping sites are interstitial, inside the layer and in the vdW gap,
whereas 125Te-NMR is sensitive to the local magnetism at the Te site. The
$\mu^{+}$SR asymmetry spectra measured in zero field at
$T=$$5.5\text{\,}\mathrm{K}$ with the initial muon spin polarization (b)
perpendicular ($\perp$) and (c) parallel ($\parallel$) to the $c$-axis of the
crystal. Coherent zero-field precession is observed only in (b), indicating
that the muons probe a magnetic field along the $c$-axis. (d) 125Te-NMR
spectral line shapes of a 2H-MoTe2 single crystal in a field of
$5.01\text{\,}\mathrm{T}$ perpendicular to the $c$-axis at different
temperatures. (e) 125Te-NMR powder line shapes at $15\text{\,}\mathrm{K}$ at
different applied fields. The lack of a pronounced temperature or field
dependence in the 125Te-NMR lines implies the absence of an intrinsic magnetic
order.

Both the muon and the 125Te nuclei have a spin of $1/2$ and are, therefore,
unaffected by quadrupolar interactions and thus act as purely magnetic probes
of the system. The 2H-MoTe2 single crystals used in this study were grown by
the Chemical Vapour Transport technique, using TeCl4 as the transporting
agent. In contrast, the previously studied samples from Ref. Guguchia et al.
(2018) were prepared using Te flux. Therefore, we expect the density and
nature of native defects in our samples to be different. This subtle
difference should affect their magnetic properties, if long-range magnetism
were induced by structural defects.

Indeed we observe a coherent precession of the muon spin in zero magnetic
field and at low temperature, see Fig. 1(b). However, our $\mu^{+}$SR
measurements show identical results to those published earlier, namely
yielding the same transition temperature and the same internal magnetic field
magnitude. The precession frequency at $5.5\text{\,}\mathrm{K}$ in Fig. 1(b)
corresponds to an internal magnetic field of
$199.1(3)\text{\,}\mathrm{m}\mathrm{T}$, consistent with the previously
reported value of $\mu_{0}H_{\mathrm{int}}=$200\text{\,}\mathrm{m}\mathrm{T}$$
Guguchia et al. (2018). Additional measurements as a function of temperature
and using low-energy muons stopping within
$10\text{\,}\mathrm{n}\mathrm{m}100\text{\,}\mathrm{n}\mathrm{m}$ from the
surface are presented in Supplementary Fig. S8 and show a similar behavior SI
. Therefore, the local magnetic properties observed with $\mu^{+}$SR in
2H-MoTe2 are robust and do not change in the near-surface region of the
crystals, nor do they depend on the details of the sample synthesis. This is
an indication that the previously reported structural defects Guguchia et al.
(2018) are unlikely to be the origin of long-range magnetic order.

Furthermore, measurements as a function of the angle between the $c$-axis and
initial muon spin polarization, as shown in Fig. 1(b,c), reveal that the
internal magnetic field sensed by the muon is oriented along the $c$-axis. A
maximal Larmor precession amplitude is detected when the initial muon spin
$\mathbf{S}\perp c$-axis (Fig. 1(b)) while it almost vanishes for
$\mathbf{S}\parallel c$-axis (Fig. 1(c)). This is a clear indication that the
local magnetic field sensed by the implanted muons points along the
crystallographic $c$-axis. A detailed analysis of the full angular dependence
between the initial muon spin and the $c$-axis is presented in Supplementary
Fig. S6 and supports this conclusion SI .

Contrary to the $\mu^{+}$SR results, we find no signatures of intrinsic
magnetism when we investigate 2H-MoTe2 with 125Te-NMR. Example spectra on
single crystalline 2H-MoTe2 in a field of $\approx$5.01\text{\,}\mathrm{T}$$
are shown in Fig. 1(d). The spectra neither exhibit a significant broadening
nor a resonance shift as a function of temperature, which is strong evidence
that the 125Te nuclear spin experiences no significant changes in the local
spin susceptibility. This is in stark contrast to the typical NMR hallmarks of
magnetic ordering, which consist of a temperature-dependent NMR line shift and
broadening, caused by the emergence of an internal magnetic field near the
magnetic transition, as observed for example in unenriched Fe Budnick et al.
(1961) or in the RAlGe material family Tay et al. (2020). Similar distinctive
NMR magnetic signatures are also expected for defect-induced magnetism, as has
been observed e.g. in the 2D system SiC Zhang et al. (2017).

Furthermore, we exclude the presence of a static magnetic field at the 125Te
site, by measuring powder line shapes at $15\text{\,}\mathrm{K}$ in different
applied fields, see Fig. 1(e). These spectra are significantly broader than
those of single crystals, a broadening which we attribute to an anisotropic
chemical shift. However, apart from a small change in the line shape, the line
width scales perfectly with the applied field. This lack of field dependence
is indicative of the non-magnetic origin of the line broadening, as it
excludes the presence of a static internal magnetic field at the 125Te site.
Therefore, these results clearly show that an intrinsic spin probe exhibits no
evidence of any long-range magnetic order in 2H-MoTe2. This is in clear
contrast to our observation with $\mu^{+}$SR and a strong indication that the
local static fields observed with muons originate from a muon-induced effect.

In order to reveal the origin of this effect, we use DFT to calculate the
ground state of 2H-MoTe2 and how it is affected by an individual implanted
muon. We start by refining our DFT calculation to reproduce the electronic
structure of 2H-MoTe2. For this purpose, we determine the three dimensional
electronic band dispersion with soft X-ray angle-resolved photoemission
spectroscopy (ARPES). Spectra along high symmetry directions at photon
energies corresponding to the $\Gamma$ and A planes are shown in Fig. 2(a-e).

Figure 2:  Comparison of soft X-ray ARPES spectra at
$\sim$$12\text{\,}\mathrm{K}$ with the calculated band structure (solid
lines). (a,b) Dispersion along $\Gamma-M$ at $\Gamma_{25}$, (c) along
$\Gamma-M$ at $\Gamma_{27}$, (d) along K-M at $\Gamma_{25}$, (e) along A-L
between $\Gamma_{26}$ and $\Gamma_{27}$ and (f) along $\Gamma$-A. (g) Constant
energy cuts at $\Gamma_{25}$. In (a) and (b) symmetric and antisymmetric bands
were probed by using $p$\- and $s$\- polarized light, respectively. There, the
width of the calculated lines represents the contribution of symmetric or
anti-symmetric orbitals to the band structure. All other spectra were probed
with circularly polarized light, which does not distinguish between symmetric
and anti-symmetric states. (h) Brillouin zone and high symmetry points of
2H-MoTe2.

These reveal sharp bands over the whole range of binding energies. We find
that the valence band maximum is formed at the $\Gamma$ point by a band
derived from states that are symmetric along $\Gamma-M$. In addition, we
measure the band dispersions in $k_{z}$ along $\Gamma-A$ as a function of
photon energy (Fig. 2(f)). Along this direction we expect that a two
dimensional state forms a flat band with no clear dispersion. In contrast, the
dispersive spectral weight in Fig. 2(f) reveals that most bands hybridize
across the vdW gap and gain a strong 3D character, despite the 2D layered
structure of 2H-MoTe2. Interestingly, the spectral weight exhibits modulations
with longer periodicity in $k_{z}$ than the period of the reciprocal lattice.
Such an effect has been previously observed in other 2H-TMDs and is connected
with the non-symmorphic symmetry group of the 2H structure of MoTe2 Weber et
al. (2018).

We compare the spectral weight in Fig. 2 to the calculated band structure.
Note that in general, ab-initio DFT of vdW materials with strong spin-orbit
coupling can be misleading, as exemplified by Td-MoTe2, where an on-site
correlation term ($U$) had to be introduced, to correctly reproduce the
experimental results Rhodes et al. (2017); Xu et al. (2018); Aryal and
Manousakis (2019). However, our DFT calculations, performed without a $U$
term, are in excellent agreement with the experimental bands. This indicates
that, unlike Td-MoTe2, in 2H-MoTe2 there are only weak correlations. We
further note that there are no signatures of intrinsic magnetism in the
calculations of bulk 2H-MoTe2.

By correctly reproducing the 3D band dispersions, we have identified a
reliable representation of 2H-MoTe2’s electronic structure and charge density
within DFT. We can now employ it to predict the behavior of implanted
$\mu^{+}$ spin probes, introduced as isolated impurities at interstitial
positions within a supercell. We consider two charge states of the muon by
adding a charge on the supercell ($+1$ for a diamagnetic $\mu^{+}$ and $+0$
for a neutral muon state) Möller et al. (2013); Bernardini et al. (2013). The
electrically neutral state of a muon corresponds to a muonium (Mu), a
hydrogen-like bound state between a muon and an electron. Figure 3 shows the
resulting stable stopping sites and the diffusion energy barrier along the
most energetically favorable paths between them.

Figure 3:  (a) Stopping sites and diffusion paths indicated in the
3$\times$3$\times$1 supercell. The muon ($\mu^{+}$) at sites b and g and
muonium (Mu) at sites b and a. Energy barrier for (b) $\mu^{+}$, and (c) Mu
diffusion along the different paths. The black lines show the ground-state
energy level of the particles at the local minima. (d) Calculated local spin
density of states for a muon stopping at site a, assuming a nonzero Hubbard
$U$.

In particular, we find that the diamagnetic $\mu^{+}$ has a stable site b
(denoted by its Wyckoff letter) inside the Mo-layer and a meta-stable site g
in the vdW gap. On the other hand, the neutral muonium state can be stable in
b or meta-stable in a.

Finally, we consider Mu- and muon-induced magnetism in these sites. Indeed,
for the Mu at the site a, and when employing a finite Hubbard term $U$ on the
hydrogen (i.e. Mu) potential, a significant spin density is induced in the
material, Fig. 3(d), while all other sites remain nonmagnetic (Supplementary
Fig. S5) SI . Such muon-induced magnetism is very similar to what was found
for hydrogen adsorbed on a monolayer of 2H-MoTe2 or on graphene Ma et al.
(2011); González-Herrero et al. (2016). As we discuss in more detail in the
Supplementary SI , such a magnetic muonium state is consistent with both the
zero-field as well as the transverse field behavior (Supplementary Fig. S7) of
the observed $\mu^{+}$SR signals. Therefore, we conclude that the 2H-MoTe2
does not exhibit an intrinsic magnetic order and that the magnetic signal
observed with $\mu^{+}$SR originates from a magnetic muonium state in site a.
As the muon can be considered a light hydrogen isotope, this induced magnetic
state should be the same for a hydrogen atom located at site a.

To conclude, using local magnetic probe measurements, ARPES and theoretical
DFT calculations, we are able to reveal that 2H-MoTe2 is a semiconductor on
the verge of becoming magnetic. The comparison of the intrinsic 125Te-NMR
probe with the implanted muon spin probe confirms that the observed local
field in $\mu^{+}$SR does not reflect an intrinsic long-range magnetic order.
Instead, we find that the muonium, being a hydrogen-like isolated impurity,
can be used to “switch-on” a local magnetic moment in this system.

Importantly, we find that the effect of muonium on the local magnetic
properties remains unchanged up to the surface of the crystals. This paves the
path towards the use of hydrogen impurities to transform 2H-MoTe2 into a
magnetic semiconductor building block for vdW heterostructures. We expect that
the use of such a material is not limited to applications in spintronics
devices Baltz et al. (2018) but may also include, e.g., Majorana
heterostructures Sau et al. (2010) or introduction of magnetic proximity at
the surface of topological materials Hesjedal and Chen (2017).

## Acknowledgements

This work is partially based on experiments performed at the Swiss Muon Source
(S$\mu$S) and Swiss Light Source (SLS), Paul Scherrer Institute, Villigen,
Switzerland. J.A.K. acknowledges support by the Swiss National Science
Foundation (SNF-Grant No. 200021_165910). I.P.R. acknowledges support from the
Ministry of Education and Science of the Russian Federation within State Task
Nr. 0721-2020-0033. The work at the University of Warwick was supported by
EPSRC, UK, through through Grants EP/M028771/1 and EP/T005963/1. The
magnetization measurements were supported by the Swiss National Science
Foundation (SNF-Grant No. 206021_139082). We are grateful to C. W. Schneider
for help with single crystal XRD measurements. We thank Z. Guguchia, and S.
Holenstein from PSI, CH, and A. Chatzichristos, D. Fujimoto, V. L. Karner, R.
M. L. McFadden, J. O. Ticknor, W. A. MacFarlane, and R. F. Kiefl from the
University of British Columbia, CA, for helpful discussions. Further, we note
that a preliminary version of these results was published as part of a PhD
thesis, Ref. Krieger (2020).

## References

  * Geim and Grigorieva (2013) A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013).
  * Rajamathi et al. (2017) C. R. Rajamathi, U. Gupta, N. Kumar, H. Yang, Y. Sun, V. Süß, C. Shekhar, M. Schmidt, H. Blumtritt, P. Werner, et al., Advanced Materials 29, 1606202 (2017).
  * Clarke et al. (1978) R. Clarke, E. Marseglia, and H. P. Hughes, Philosophical Magazine B 38, 121 (1978).
  * Sun et al. (2015) Y. Sun, S.-C. Wu, M. N. Ali, C. Felser, and B. Yan, Physical Review B 92, 161107(R) (2015).
  * Qi et al. (2016) Y. Qi, P. G. Naumov, M. N. Ali, C. R. Rajamathi, W. Schnelle, O. Barkalov, M. Hanfland, S.-C. Wu, C. Shekhar, Y. Sun, et al., Nature Communications 7, 11038 (2016).
  * Xu et al. (2018) N. Xu, Z. W. Wang, A. Magrez, P. Bugnon, H. Berger, C. E. Matt, V. N. Strocov, N. C. Plumb, M. Radovic, E. Pomjakushina, K. Conder, J. H. Dil, J. Mesot, R. Yu, H. Ding, and M. Shi, Physical Review Letters 121, 136401 (2018).
  * Lezama et al. (2014) I. G. Lezama, A. Ubaldini, M. Longobardi, E. Giannini, C. Renner, A. B. Kuzmenko, and A. F. Morpurgo, 2D Materials 1, 021002 (2014).
  * Xu et al. (2015) H. Xu, S. Fathipour, E. W. Kinder, A. C. Seabaugh, and S. K. Fullerton-Shirey, ACS Nano 9, 4900 (2015).
  * Ding et al. (2018) Y. Ding, N. Zhou, L. Gan, X. Yan, R. Wu, I. H. Abidi, A. Waleed, J. Pan, X. Ou, Q. Zhang, et al., Nano Energy 49, 200 (2018).
  * Song et al. (2016) S. Song, D. H. Keum, S. Cho, D. Perello, Y. Kim, and Y. H. Lee, Nano Letters 16, 188 (2016).
  * Wang et al. (2017) Y. Wang, J. Xiao, H. Zhu, Y. Li, Y. Alsaid, K. Y. Fong, Y. Zhou, S. Wang, W. Shi, Y. Wang, et al., Nature 550, 487 (2017).
  * Coelho et al. (2019) P. M. Coelho, H.-P. Komsa, K. Lasek, V. Kalappattil, J. Karthikeyan, M.-H. Phan, A. V. Krasheninnikov, and M. Batzill, Advanced Electronic Materials 0, 1900044 (2019).
  * Guguchia et al. (2018) Z. Guguchia, A. Kerelsky, D. Edelberg, S. Banerjee, F. v. Rohr, D. Scullion, M. Augustin, M. Scully, D. A. Rhodes, Z. Shermadini, et al., Science Advances 4, eaat3672 (2018).
  * Ma et al. (2011) Y. Ma, Y. Dai, M. Guo, C. Niu, J. Lu, and B. Huang, Physical Chemistry Chemical Physics 13, 15546 (2011).
  * Kanoun (2018) M. B. Kanoun, Journal of Alloys and Compounds 748, 938 (2018).
  * (16) Additional details are available in the Supplementary Information, containing Refs. Amato et al. (2017); Yaouanc and de Réotier (2011); Suter and Wojek (2012); Prokscha et al. (2008); Eckstein (1991); Harris et al. (2002); Strocov et al. (2010, 2014, 2012); Böker et al. (2001); Kresse and Furthmüller (1996); Kresse and Joubert (1999); Perdew et al. (1996); Möller et al. (2013); Bernardini et al. (2013); Henkelman and Jónsson (2000); Ozaki (2003); Momma and Izumi (2011); Guguchia et al. (2018); Tongay et al. (2012); Ma et al. (2011); González-Herrero et al. (2016); Herak et al. (2013); Patterson (1988); Cox (2003); Tian et al. (2020, 2021).
  * Budnick et al. (1961) J. I. Budnick, L. J. Bruner, R. J. Blume, and E. L. Boyd, Journal of Applied Physics 32, S120 (1961).
  * Tay et al. (2020) D. Tay, T. Shang, P. Puphal, E. Pomjakushina, H.-R. Ott, and T. Shiroka, Physical Review B 102, 241109(R) (2020).
  * Zhang et al. (2017) Z. T. Zhang, D. Dmytriieva, S. Molatta, J. Wosnitza, Y. Wang, M. Helm, S. Zhou, and H. Kühne, Physical Review B 95, 085203 (2017).
  * Weber et al. (2018) F. Weber, R. Hott, R. Heid, L. L. Lev, M. Caputo, T. Schmitt, and V. N. Strocov, Physical Review B 97, 235122 (2018).
  * Rhodes et al. (2017) D. Rhodes, R. Schönemann, N. Aryal, Q. Zhou, Q. R. Zhang, E. Kampert, Y.-C. Chiu, Y. Lai, Y. Shimura, G. T. McCandless, J. Y. Chan, D. W. Paley, J. Lee, A. D. Finke, J. P. C. Ruff, S. Das, E. Manousakis, and L. Balicas, Physical Review B 96, 165134 (2017).
  * Aryal and Manousakis (2019) N. Aryal and E. Manousakis, Physical Review B 99, 035123 (2019).
  * Möller et al. (2013) J. S. Möller, D. Ceresoli, T. Lancaster, N. Marzari, and S. J. Blundell, Physical Review B 87, 121108(R) (2013).
  * Bernardini et al. (2013) F. Bernardini, P. Bonfà, S. Massidda, and R. De Renzi, Physical Review B 87, 115148 (2013).
  * González-Herrero et al. (2016) H. González-Herrero, J. M. Gómez-Rodríguez, P. Mallet, M. Moaied, J. J. Palacios, C. Salgado, M. M. Ugeda, J.-Y. Veuillen, F. Yndurain, and I. Brihuega, Science 352, 437 (2016).
  * Baltz et al. (2018) V. Baltz, A. Manchon, M. Tsoi, T. Moriyama, T. Ono, and Y. Tserkovnyak, Reviews of Modern Physics 90, 015005 (2018).
  * Sau et al. (2010) J. D. Sau, S. Tewari, R. M. Lutchyn, T. D. Stanescu, and S. Das Sarma, Physical Review B 82, 214509 (2010).
  * Hesjedal and Chen (2017) T. Hesjedal and Y. Chen, Nature Materials 16, 3 (2017).
  * Krieger (2020) J. A. Krieger, Doctoral Thesis, ETH Zurich (2020), URL https://dx.doi.org/10.3929/ethz-b-000460760.
  * Amato et al. (2017) A. Amato, H. Luetkens, K. Sedlak, A. Stoykov, R. Scheuermann, M. Elender, A. Raselli, and D. Graf, Review of Scientific Instruments 88, 093301 (2017).
  * Yaouanc and de Réotier (2011) A. Yaouanc and P. de Réotier, _Muon Spin Rotation, Relaxation, and Resonance: Applications to Condensed Matter_ , International Series of Monographs on Physics (OUP Oxford, 2011), ISBN 978-0-19-959647-8.
  * Suter and Wojek (2012) A. Suter and B. M. Wojek, Physics Procedia 30, 69 (2012).
  * Prokscha et al. (2008) T. Prokscha, E. Morenzoni, K. Deiters, F. Foroughi, D. George, R. Kobler, A. Suter, and V. Vrankovic, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 595, 317 (2008).
  * Eckstein (1991) W. Eckstein, _Computer Simulation of Ion-Solid Interactions_ , Springer Series in Materials Science (Springer-Verlag, Berlin Heidelberg, 1991), ISBN 978-3-642-73515-8.
  * Harris et al. (2002) R. K. Harris, E. D. Becker, S. M. C. d. Menezes, R. Goodfellow, and P. Granger, Magnetic Resonance in Chemistry 40, 489 (2002).
  * Strocov et al. (2010) V. N. Strocov, T. Schmitt, U. Flechsig, T. Schmidt, A. Imhof, Q. Chen, J. Raabe, R. Betemps, D. Zimoch, J. Krempasky, et al., Journal of Synchrotron Radiation 17, 631 (2010).
  * Strocov et al. (2014) V. N. Strocov, X. Wang, M. Shi, M. Kobayashi, J. Krempasky, C. Hess, T. Schmitt, and L. Patthey, Journal of Synchrotron Radiation 21, 32 (2014).
  * Strocov et al. (2012) V. N. Strocov, M. Shi, M. Kobayashi, C. Monney, X. Wang, J. Krempasky, T. Schmitt, L. Patthey, H. Berger, and P. Blaha, Physical Review Letters 109, 086401 (2012).
  * Böker et al. (2001) T. Böker, R. Severin, A. Müller, C. Janowitz, R. Manzke, D. Voß, P. Krüger, A. Mazur, and J. Pollmann, Physical Review B 64, 235305 (2001).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Physical Review B 54, 11169 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Physical Review B 59, 1758 (1999).
  * Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
  * Henkelman and Jónsson (2000) G. Henkelman and H. Jónsson, The Journal of Chemical Physics 113, 9978 (2000).
  * Ozaki (2003) T. Ozaki, Physical Review B 67, 155108 (2003).
  * Momma and Izumi (2011) K. Momma and F. Izumi, Journal of Applied Crystallography 44, 1272 (2011).
  * Tongay et al. (2012) S. Tongay, S. S. Varnoosfaderani, B. R. Appleton, J. Wu, and A. F. Hebard, Applied Physics Letters 101, 123105 (2012).
  * Herak et al. (2013) M. Herak, A. Zorko, M. Pregelj, O. Zaharko, G. Posnjak, Z. Jagličić, A. Potočnik, H. Luetkens, J. van Tol, A. Ozarowski, H. Berger, and D. Arcon, Physical Review B 87, 104413 (2013).
  * Patterson (1988) B. D. Patterson, Reviews of Modern Physics 60, 69 (1988).
  * Cox (2003) S. F. J. Cox, Journal of Physics: Condensed Matter 15, R1727 (2003).
  * Tian et al. (2020) Y. Tian, N. Ghassemi, and J. H. Ross, Physical Review B 102, 165149 (2020).
  * Tian et al. (2021) Y. Tian, Y. Zhu, R. Li, Z. Mao, and J. H. Ross, Physical Review B 104, L041105 (2021).

## Supplemental Material: Hydrogen-impurity induced unconventional magnetism
in semiconducting molybdenum ditelluride

## Sample growth and characterization

High-quality single crystals of 2H-MoTe2 were grown by chemical vapor
transport method. Stoichiometric amounts of the starting materials, Mo and Te,
along with the transport agent TeCl4
($3\text{\,}\mathrm{m}\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}^{3}$), were
mixed together in a quartz ampule and sealed in vacuum. The quartz tube was
then placed for 3 weeks in a two-zone horizontal tube furnace with the charge
end at $800\text{\,}\mathrm{\SIUnitSymbolCelsius}$ and the end where the
crystals form at $750\text{\,}\mathrm{\SIUnitSymbolCelsius}$. The high quality
of the single crystals and the crystal structure was confirmed by x-ray
diffraction (XRD) and x-ray photoelectron spectroscopy (XPS), as shown in
Figure S1. The additional XPS peaks around $-460\text{\,}\mathrm{e}\mathrm{V}$
are the Te M4,5-N4,5N4,5 Auger lines.

Figure S1: (a) The Laue back reflection X-ray diffraction pattern for one of
the single crystals. (b) X-ray (0,0,n) diffraction peaks. The inset shows a
photograph of one of the single crystals. (c) XPS core level spectrum measured
at $h\nu=$950\text{\,}\mathrm{e}\mathrm{V}$$. The blue (red) lines show the
expected positions of the Mo (Te) peaks.

### Magnetization measurements

Magnetometry measurements have been performed with a Quantum Design
superconducting quantum interference device (SQUID) magnetic properties
measurement system (MPMS) and are shown in Fig. S2.

Figure S2:  Magnetization measurements on 2H-MoTe2. (a) Zero field cooled
(ZFC) magnetization in an applied field of $100\text{\,}\mathrm{m}\mathrm{T}$.
(b) Zoom into the low temperature region of (a). The solid line shows the fit
of a Curie-Weiss law ($T_{c}=$-0.9(5)\text{\,}\mathrm{K}$$). The arrow shows
the onset of the magnetic zero field oscillations observed in $\mu^{+}$SR. (c)
Field cooled (below $50\text{\,}\mathrm{K}$) and zero field cooled
magnetization for a $10\text{\,}\mathrm{m}\mathrm{T}$ out-of-plane field.

For an applied in-plane field of $100\text{\,}\mathrm{m}\mathrm{T}$, the
sample shows a Curie-Weiss like behavior. In contrast, if the field is applied
out-of-plane there is a signature of a transition around
$5\text{\,}\mathrm{K}$, see Fig. S2(b). Moreover, for an out-of-plane field of
$10\text{\,}\mathrm{m}\mathrm{T}$ this transition features a small hysteresis
in the magnetization signal, Fig. S2(c). Its magnitude is too small for an
intrinsically ferromagnetic sample and it cannot account for the observed
magnetic transition in $\mu^{+}$SR, which appears at much higher temperatures
($\sim$40\text{\,}\mathrm{K}$$). However, it could be attributed to various
origins, including magnetic defects Guguchia et al. (2018) or magnetic grain-
boundaries Tongay et al. (2012).

## ARPES

The soft X-ray ARPES experiments were performed on the X03MA beamline at the
Swiss Light Source of the Paul Scherrer Institute, Villigen, Switzerland
Strocov et al. (2010, 2014). The analyzer slit was oriented along the incoming
beam direction and the photon energy $h\nu$ was tuned in the range of
$320\text{\,}\mathrm{e}\mathrm{V}$ to $900\text{\,}\mathrm{e}\mathrm{V}$. The
high photon energies allow a better $k_{z}$ resolution Strocov et al. (2012)
than in previous investigations of 2H-MoTe2 Böker et al. (2001). The combined
beamline and analyzer resolution was better than
$0.15\text{\,}\mathrm{e}\mathrm{V}$ at
$h\nu=$$900\text{\,}\mathrm{e}\mathrm{V}$ and better than
$63\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ around
$h\nu\sim$$500\text{\,}\mathrm{e}\mathrm{V}$. The samples were cleaved in-situ
and kept at $\sim$12\text{\,}\mathrm{K}$$ in a vacuum lower than
${10}^{-10}\text{\,}\mathrm{m}\mathrm{b}\mathrm{a}\mathrm{r}$. We only present
the binding energy with respect to the valence band maximum, because 2H-MoTe2
tends to statically charge upon x-ray illumination, which causes a rigid shift
of the spectrum, see below.

### Out-of plane dispersion

In the soft X-rays energy range the photoelectrons have a larger escape depth
compared to standard ultraviolet ARPES. This reduces the $k_{z}$ broadening
and allows to determine the out-of plane dispersion. For a measured
photoelectron kinetic energy $E_{\mathrm{kin}}$, the perpendicular momentum is
given by:

$\mathbf{k}_{z}=\sqrt{\frac{2m_{\mathrm{e}}\left(E_{\mathrm{kin}}-V_{000}\right)}{\hbar^{2}}-\mathbf{k}_{\parallel}}+\mathbf{p}_{\gamma,z},$
(S1)

where $m_{\mathrm{e}}$ is the electron mass, $\mathbf{k}_{\parallel}$ the in-
plane momentum, $V_{000}$ the inner potential and $\mathbf{p}_{\gamma,z}$ the
component of the photon momentum ($h\nu/c$) perpendicular to the sample
surface 111The formulas for the specific measurement geometry can be found in
Ref. Strocov et al. (2014). The inner potential is a material and energy
dependent parameter. We have estimated it by comparing the measured band
structure with the expected periodicity of the Brillouin zone (cf. Fig. 1(f)
in the main paper). We used a value of
$V_{000}\approx$10\text{\,}\mathrm{e}\mathrm{V}$$, which gives a fair
agreement between theory and experiment, but slightly differs from the
previously used $V_{000}=$16(1)\text{\,}\mathrm{e}\mathrm{V}$$ Böker et al.
(2001) .

### Static charging

Because 2H-MoTe2 is a semiconductor, it tends to statically charge upon x-ray
irradiation at low temperature. This is shown in Fig. S3, which displays the
binding energy of several features in the band structure as a function of
photon flux.

Figure S3: (a) Binding energy of several features indicated in (b) as a
function of x-ray flux. (c) Angle integrated spectra corresponding to (a),
after rigidly aligning the valence band maximum.

The binding energies were calibrated by measuring the Fermi level of Au. For
some measurements the sample was additionally irradiated with a near UV laser.
This excites some charge carriers and helps to significantly reduce the
charging (Fig S3(a)). In both cases, the charging can be assumed to be a rigid
shift of the spectrum. This is evident from the linear dependence on the flux
in Fig S3(a) and from the collapse of angle integrated spectra after alignment
of the valence band maximum in Fig S3(c). Therefore, in the main paper we only
show the binding energy with respect to the valence band maximum.

### n-type samples

All measurements shown in the main text were performed on samples that were
freshly cleaved in-situ at $\sim$$12\text{\,}\mathrm{K}$ directly before the
measurement. However, after a couple of hours of ARPES measurements the
2H-MoTe2 samples usually degrade and the material becomes slightly $n$-type.
This is exemplified by the spectra in Fig. S4, where a small spectral weight
above the valence band was observed.

Figure S4: (a,c) Spectra measured on a sample that became n-type over time.
(b,d) Zooms into the region above the valence band maximum. The dashed line
indicates the estimated, effective position of the Fermi level.

These measurements clearly confirm that the band gap in 2H-MoTe2 is indirect.

## Density-functional theory calculations

The electronic structure was calculated with Density-Functional Theory (DFT)
and the Projected Augmented Wave (PAW) method from the VASP code Kresse and
Furthmüller (1996); Kresse and Joubert (1999). Spin-orbit coupling was treated
within the scalar-relativistic approach by the second variation method. The
structure was relaxed and the exchange-correlation functional was chosen to be
Perdew-Burke-Ernzerhof (PBE) Perdew et al. (1996). The parameter for the
relaxed cell used in the ARPES calculation were given by
$a=$3.491\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$,
$c=$13.647\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$, where the Mo and Te sit
at the $2c$ and $4f$ Wyckoff positions of space group 194, respectively, with
relative $z$-coordinates $0.116$ and $0.383$ for Te. Along $\Gamma$-M an
orbital decomposition into symmetric
($p_{z}$,$p_{y}$,$d_{z^{2}}$,$d_{x^{2}-y^{2}}$,$d_{zy}$) and antisymmetric
($p_{x}$,$d_{xy}$,$d_{xz}$) orbitals was preformed. The orbital character does
not fully separate into symmetric and antisymmetric bands because of spin-
orbit coupling. Nevertheless, when using $p-$ ($s-$) polarized light [Fig.
2(a) (2(b)) in the main text], only the symmetric (antisymmetric) part of the
electronic wave function, which primarily consists of these orbitals, will be
probed.

### Muon site determination

In order to simulate a single muon impurity, we used a hydrogen
pseudopotential in the relaxed $3\times 3\times 1$ supercell of 2H-MoTe2. The
assumed lattice constants in the $3\times 3\times 1$ supercell,
c$=$13.647\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$ and a corresponding
a$=$10.4745\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$, were the same as for the
band structure calculations. The charge state of the muon is accounted for by
assuming a neutral/charged supercell Möller et al. (2013); Bernardini et al.
(2013). Therein, we preformed nudged elastic band (NEB) Henkelman and Jónsson
(2000) as well as density of states calculations on the basis of the Linear
Combination of Pseudo Atomic Orbitals (LCPAO) approach as implemented in the
OPENMX code Ozaki (2003). The muon ground state energies were estimated by
approximating the NEB diffusion energy barriers with an anisotropic harmonic
oscillator potential. The crystal structures were visualized with VESTA Momma
and Izumi (2011). Some further details regarding the stable muon positions are
listed in Tab. 1.

Species | Wyckoff position | frac. coordinates |  $\Delta E$ $(eV)$ | ZPE $(eV)$
---|---|---|---|---
$\mu^{+}$ | $2b$ | 0 | 0 | 0.25 | 0 | 0.57
$\mu^{+}$ | $6g$ | 0.5 | 0 | 0.5 | 1.5 | 0.22
Mu | $2b$ | 0 | 0 | 0.25 | 0 | 0.54
Mu | $2a$ | 0 | 0 | 0.5 | 1.24 | 0.26
Table 1: Details of the identified stable and meta-stable site for several
implanted probe species. $\Delta E$ denotes the difference in the total energy
of the relaxed supercell compared to the most stable site and ZPE is the
approximate ground state energy of the muon.

We further show the local spin density of states at the different stopping
sites in Fig. S5.

Figure S5: Total density of states and local spin density of states calculated
for a muon (i.e. a hydrogen potential in a charged supercell) at site b (a,b),
site a (c,d), and site g (e,f), respectively. Note that site a is not a stable
stopping site for a diamagnetic muon.

It is evident from this figure, that there is no induced magnetism at the
stable muon stopping site b and g, see Fig. S5(b,f). For site b, the absence
of an induced spin polarization can be understood because of the strong
pairwise interaction of the muon with the neighboring Mo atoms. However, in
the case of a muonium (and also a muon) sitting at the hexagonal site inside
the vdW gap (site a) and when employing a finite Hubbard term $U$, a
significant spin density is induced, see Fig. S5(c,d). This is very similar to
what was found for hydrogen adsorbed on a monolayer 2H-MoTe2 or on graphene Ma
et al. (2011); González-Herrero et al. (2016). Note that site a is only a
stable stopping site of muonium but not of a muon. As argued in more detail in
the main text and below, such a muonium induced magnetic state can explain the
observed coherent muon spin precession, while being consistent with the
absence of a magnetic signal in magnetization and 125Te-NMR.

## Muon spin spectroscopy

The $\mu^{+}$SR experiments were performed on the general purpose spectrometer
(GPS) at the Swiss Muon Source, Paul Scherrer Institute, Switzerland Amato et
al. (2017). The temperature was controlled using a 4He flow cryostat. In order
to vary the orientation of the sample with respect to the initial muon spin
direction, the sample was mounted on a mechanical rotation-stage. Nearly
$100\text{\,}\mathrm{\char 37\relax}$ spin polarized muons were implanted. The
asymmetric emission of the muon’s decay positrons (lifetime
$\tau_{\mu}=$2.2\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$$) enables us
to use them as a local magnetic probe Yaouanc and de Réotier (2011). Note,
that at any given moment, there is only one muon in the sample. Therefore, we
can exclude that a muon probes magnetism induced by a nearby Mu or $\mu^{+}$.
The data were analyzed with musrfit software suite Suter and Wojek (2012).

### Zero field angle dependence

In order to determine the variation in the zero field (ZF) $\mu^{+}$SR spectra
as a function of sample orientation, the asymmetry spectra (such as the ones
shown in Fig. 1(b,c) of the main paper) were fitted to a sum of a Gaussian
damped oscillating part and two exponentially damped components:

$A(t)=A_{1}e^{-\frac{1}{2}\left({\sigma}t\right)^{2}}\cos\left(\gamma_{\mu}Bt+\varphi\right)+A_{2}e^{-\lambda_{\mathrm{slow}}t}+A_{3}e^{-\lambda_{\mathrm{fast}}t}.$
(S2)

We further assumed that the internal field $B$ and the depolarization rate
$\sigma$ of the oscillating part to be independent of the measured projection
of the polarization. Figure S6 shows the resulting oscillating ZF asymmetry
$A_{1}$ in different sets of detectors (corresponding to different projections
of the muon spin polarization) as a function of the sample orientation.

Figure S6: Oscillating asymmetry in different sets of detectors, measured for
different orientations of the sample. The inset depicts the measurement
geometry.

In order to understand this rotation dependence, one has to consider the
Larmor precession of the muon polarization $\mathbf{P}(t)$, defined as

$\frac{\rm{d}{\mathbf{P}}(t)}{\rm
dt}=\gamma_{\mu}\mathbf{P}(t)\wedge\mathbf{B}.$ (S3)

Here, $\mathbf{B}$ denotes the local magnetic field vector and
$\gamma_{\mu}=2\pi\times$135.5\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}\mathrm{/}\mathrm{T}$$
is the gyromagnetic ratio of the muon. The resulting precession can be written
as Herak et al. (2013); Yaouanc and de Réotier (2011):

$\displaystyle\mathbf{P}(t)$
$\displaystyle=\mathbf{p}^{\parallel}+\mathbf{p}^{\perp}_{1}\cos(\gamma_{\mu}Bt)+\mathbf{p}^{\perp}_{2}\sin(\gamma_{\mu}Bt),$
(S4) $\displaystyle\mathbf{p}^{\parallel}$
$\displaystyle=\frac{(\mathbf{P}_{0}\cdot\mathbf{B})\mathbf{B}}{B^{2}},\qquad\mathbf{p}^{\perp}_{1}=\mathbf{P}_{0}-\mathbf{p}^{\parallel},\qquad\mathbf{p}^{\perp}_{2}=\mathbf{p}^{\perp}_{1}\wedge\frac{\mathbf{B}}{B},$
(S5)

where $\mathbf{P}_{0}$ is initial polarization vector and $B=\|\mathbf{B}\|$
denotes the magnitude of the field.

We choose the coordinate system $(x,y,z)$ depicted in the inset of Fig. S6,
with ${\mathbf{P}_{0}\propto(0,1,1)}$, and the local field
$\mathbf{B}=B(0,\sin{(\theta+\alpha)},\cos{(\theta+\alpha)})$. The angle
$\theta$ denotes the rotation of the sample’s $c$-axis with respect to the
detectors and $\alpha$ is a hypothetical offset between the local magnetic
field seen by the muon and the $c$-axis of the crystal.

The asymmetry probed by a set of opposite detectors is proportional to the
projection of the polarization along that direction. Therefore, in the inset
of Fig. S6 the asymmetry in the up-down detectors is $A_{\text{up-
do}}\propto\mathbf{P}(t)\cdot\hat{\mathbf{e}}_{y}$, and the one in the forward
backward detectors $A_{\text{fw-
ba}}\propto\mathbf{P}(t)\cdot\hat{\mathbf{e}}_{z}$, where
$\hat{\mathbf{e}}_{y}$ and $\hat{\mathbf{e}}_{z}$ are the unit vectors in the
$y$ and $z$ directions.

Inserting these conventions into Eq. S4, we find that the oscillating part of
the asymmetry is given by $\mathbf{p}^{\perp}_{1}$ and that accordingly

$\displaystyle A_{\text{up-do}}$ $\displaystyle\propto
1-\sqrt{2}\sin{(\alpha+\theta+\pi/4)}\sin{(\theta+\alpha)},$ (S6)
$\displaystyle A_{\text{fw-ba}}$ $\displaystyle\propto
1-\sqrt{2}\sin{(\alpha+\theta+\pi/4)}\cos{(\theta+\alpha)}.$ (S7)

These curves are depicted as solid lines in Fig. S6, where
$\alpha=-$4(1)\text{\,}\mathrm{\SIUnitSymbolDegree}$$ was determined using
nonlinear least-squares optimization. As the sample was aligned with respect
to the detector frame by hand, this small value of $\alpha$ could easily
correspond to a small misalignment of the sample. Therefore, we conclude that
within our measurement precision, the local field points along the $c$-axis of
the crystal.

### Field dependence in transverse field

We have also determined the evolution of the muon oscillation frequencies as a
function of a transverse field applied along the crystallographic $c$-axis.
The fast Fourier transforms of representative spectra are shown in Fig. S7(a).

Figure S7: (a) Fourier transform of the $\mu^{+}$SR asymmetry spectra measured
in different applied fields. The curves have been offset for clarity. (b)
Evolution of the oscillation frequencies as a function of the applied field.

We find that there is one component in the spectra, which oscillates at a
frequency corresponding to the applied field, and which allows an accurate
calibration of applied field ($B_{\rm ext}$). This component is partially
attributed to a background contribution from a silver plate, onto which the
sample was mounted. Other contributions to this component could also include
some diamagnetic muons, which experience no strong magnetic field in 2H-MoTe2.
In addition, we observe other oscillating components which split up into peaks
at $\gamma_{\mu}|B_{loc}\pm B_{ext}|$, where
$B_{loc}\sim$200\text{\,}\mathrm{m}\mathrm{T}$$ is the field corresponding to
the ZF oscillation. The accurate values of these frequencies were extracted by
fitting the $\mu^{+}$SR time spectra to a sum of Gaussian damped oscillating
components plus a non-oscillating Gaussian damped contribution:

$A(t)=A_{0}e^{-\frac{1}{2}\left({\sigma_{0}}t\right)^{2}}+\sum_{n}{A_{n}e^{-\frac{1}{2}\left({\sigma_{n}}t\right)^{2}}\cos\left(\gamma_{\mu}B_{n}t+\varphi\right)}.$
(S8)

The resulting field values $B_{n}$ are shown in Fig. S7(b). We note that the
internal field $B_{loc}$ seems to point equally likely along the positive and
negative $c$-axis direction. At a first glance, this seems incompatible with a
muonium state, to which we attribute this signal in the main part of the
paper. Usually the behavior of muonium in a field is dominated by the
hyperfine coupling and described by the Hamiltonian

$\mathcal{H}=\hbar\mathbf{I}_{\mu}A\mathbf{S}_{e}+\hbar\gamma_{e}\mathbf{S}_{e}\cdot\mathbf{B}_{\mathrm{ext}}+\hbar\gamma_{\mu}\mathbf{I}_{\mu}\cdot\mathbf{B}_{\mathrm{ext}},$
(S9)

where $A$ is the hyperfine coupling tensor, $\mathbf{I}_{\mu}$, and
$\mathbf{S}_{e}$ are the spins of the muon and the bound electron and
$\gamma_{e}=2\pi\times$28\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}\mathrm{/}\mathrm{T}$$
is the electron’s gyromagnetic ratio. In a system without locally induced
magnetism, this results in the presence of four different oscillation
frequencies, which scale non-linearly as a function of field, see e.g. Refs.
Patterson (1988); Cox (2003); Yaouanc and de Réotier (2011).

However, in the case of muonium induced magnetism, we can assume that the
electron spin $\mathbf{S}_{e}$ is locked along its local easy axis, and
remains unaffected by the (relatively small) applied fields. Then Eq. S9 can
be simplified to,

$\mathcal{H}=\hbar\mathbf{I}_{\mu}A\mathbf{S}_{e}+\hbar\gamma_{\mu}\mathbf{I}_{\mu}\cdot\mathbf{B}_{\mathrm{ext}}=\hbar\gamma_{\mu}\mathbf{I}_{\mu}\cdot(\mathbf{B}_{loc}+\mathbf{B}_{\mathrm{ext}}),$
(S10)

where $A\mathbf{S}_{e}/\gamma_{\mu}=:\mathbf{B}_{loc}$ is an effective local
magnetic field characterizing the ZF oscillation. The field dependence of the
oscillation frequencies resulting from Eq. S10 is fully consistent with our
observation presented in Figs. S6, S7, with
$\mathbf{B}_{loc}=(\pm$200\text{\,}\mathrm{m}\mathrm{T}$)\mathbf{\hat{c}}$,
where $\mathbf{\hat{c}}$ is a unit vector along the $c$-axis.

### Low-energy $\mu$SR

We also investigate the evolution of the observed magnetism within
$\sim$100\text{\,}\mathrm{n}\mathrm{m}$$ from the surface using low energy
$\mu^{+}$SR at the $\mu$E4 beamline of the Swiss Muon Source, Paul Scherrer
Institute, Switzerland Prokscha et al. (2008). For these measurements,
multiple single crystals were glued on a Ni coated sample plate with the
$c$-axis parallel to the applied field, but perpendicular to the initial muon
spin. The implantation depth was controlled by varying the kinetic energy,
$E$, of the incident muon beam. The resulting stopping profiles were
calculated with the Trim.SP software Eckstein (1991) and are shown in Fig.
S8(b).

Figure S8: (a) $\mu^{+}$SR time spectrum measured in zero field at
$10\text{\,}\mathrm{K}$ with an implantation energy of
$22.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$. (b) Simulated implantation
profiles for different implantation energies. The average stopping depths
beneath the surface for $3\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ or
$22.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ are
$18\text{\,}\mathrm{n}\mathrm{m}$ and $99\text{\,}\mathrm{n}\mathrm{m}$,
respectively. (c) Representative weak transverse field asymmetry spectra at
$22.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ corresponding to (d). (d)
Fraction of diamagnetic muons as a function of temperature, normalized at RT.
The dashed line shows the bulk data from Ref. Guguchia et al. (2018). The
solid lines in (a,c) show fits to the data.

At $10\text{\,}\mathrm{K}$ and
$E=$$22.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ (corresponding to an average
implantation depth of $99\text{\,}\mathrm{n}\mathrm{m}$), the muon spins
exhibit a coherent precession in zero field (ZF) similar to what we observed
in the bulk, see Fig. S8(a). This spectrum was fitted to the following
function:

$A(t)=A_{1}e^{-\frac{1}{2}\left({\sigma_{\mathrm{osc}}}t\right)^{2}}\cos\left(\gamma_{\mu}B_{\mathrm{int}}t\right)+A_{2}e^{-\lambda_{\mathrm{slow}}t}+A_{3}e^{-\frac{1}{2}\left({\sigma_{\mathrm{fast}}}t\right)^{2}}.$
(S11)

The first term describes the muons that oscillate in a local magnetic field
$B_{\mathrm{int}}$. The second term represents the slowly depolarized
(nonmagnetic) signal and the third accounts for the strongly damped part due
to background contribution from muons stopping in the Ni coated sample plate.
This analysis is similar to that used in Ref. Guguchia et al. (2018), except
for the background contribution. We find that the oscillation frequency in
Fig. S8(a) corresponds to an internal static field
$B_{\mathrm{int}}=$$200(2)\text{\,}\mathrm{m}\mathrm{T}$. Within the
statistical error of the measurement, this is consistent with what is found in
the bulk Guguchia et al. (2018).

We then use a small magnetic field ($5\text{\,}\mathrm{m}\mathrm{T}$,
transverse to the initial muon spin) to identify the fraction of diamagnetic
muons, which is proportional to the non-magnetic fraction of the sample [Fig.
S8(c,d)]. These spectra were fitted with a single, Gaussian damped oscillating
polarization.

$A(t)=A_{0}e^{-\frac{1}{2}\left(\sigma
t\right)^{2}}\cos\left(\gamma_{\mu}Bt\right).$ (S12)

This function only describes the muons in a nonmagnetic environment
oscillating at a frequency close to the Larmor frequency in the applied field.
Any magnetic or strongly damped signal (cf. Eq. S11) leads to an effective
loss of $A_{0}$. The background contribution to $A_{0}$ was determined by
measuring an empty Ni plate with the same beamline settings. This temperature
independent value was subtracted from the measured $A_{0}$ and the resulting
curve is shown in Fig. S8(d). Independent of depth and within the experimental
uncertainties, all points in Fig. S8(d) closely follow the bulk behavior
reported in Ref. Guguchia et al. (2018). Therefore, we conclude that the bulk
magnetic properties observed with $\mu^{+}$SR in 2H-MoTe2 do not change in the
near-surface region of the crystals.

##  125Te-NMR

The 125Te-NMR signals were collected by means of spin-echo detection in a
$\pi/2$–$\pi$ sequence, using a
$62.5\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}125\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$
duplexer, with a $63\text{\,}\mathrm{d}\mathrm{B}$ preamplifier. The line
shapes were then extracted via a fast Fourier transformation of those signals.
Typical pulse widths of $5\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$ and
repetition times of $20\text{\,}\mathrm{s}$ were used, with each line shape
corresponding to more than 20’000 records. Measurements on a single crystal
were performed in a field of $5.007\text{\,}\mathrm{T}$ perpendicular to the
$c$-axis, while those on powders were performed in fields of
$7.067\text{\,}\mathrm{T}$ and $4.011\text{\,}\mathrm{T}$. The 125Te reference
frequencies were calibrated using 27Al-NMR on aluminum foil, which is known to
have a very small shift as a function of temperature Harris et al. (2002).

### Powder line shapes

Figure S9 shows the detailed temperature dependence of the powder line shape
of 2H-MoTe2 in $7.067\text{\,}\mathrm{T}$. Similar to the spectra of the
single crystal measurements shown in the main part of the paper, here, too, we
observe no significant peak-shift or broadening, that would be the two common
fingerprints of a magnetic transition.

Figure S9: 125Te-NMR powder line shapes of 2H-MoTe2 in a field of
$7.067\text{\,}\mathrm{T}$ at different temperatures. The spectra have been
offset for clarity. The black dashed line indicates the expected position of
the reference line.

We note that there are some small changes in the line shape as a function of
temperature. However, these are most likely due to a change in spectral weight
between different portions of the resonance line, owing to a difference in
their spin-lattice relaxation (SLR) temperature dependence.

### Spin-lattice relaxation

The spin-lattice relaxation rate ${T_{1}}^{-1}$ of the 125Te moments is shown
in Fig. S10.

Figure S10: 125Te-NMR spin lattice relaxation rate in 2H-MoTe2 in a field of
$7.067\text{\,}\mathrm{T}$ as a function of temperature. The dashed line
indicates the activated Arrhenius-like behavior at high temperature.

At high temperature, the SLR shows an exponential dependence on temperature
($1/T_{1}\propto e^{\Delta/T}$). This is a typical behavior observed in
semiconductors. However, we find that here the gap
$\Delta=$27.7(3)\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ is much smaller than
the spectroscopic band gap, which is on the order of
$\sim$1\text{\,}\mathrm{e}\mathrm{V}$$, cf. Fig. S4. This difference is likely
due to, e.g., the presence of a small amount of impurity-induced in-gap
states, which have too little spectral weight to be detected with ARPES, but
still can depolarize the 125Te nuclear spins. At present, we are unable to
determine the mechanism for the increase in relaxation rate at low
temperatures in 2H-MoTe2. However, we note that a similar increase in
relaxation was observed in the analogue compounds ZrTe2 and ZrSiTe Tian et al.
(2020, 2021). In both cases, the increase in relaxation rate is suggested to
arise from a contribution from Dirac electrons. However, 2H-MoTe2 is not a
Dirac system, so the fact that there is a similar increase in relaxation rate
suggests that the true origin of this effect might be due to some new
mechanism, yet to be understood. Finally, we note that the SLR has a minimum
around $40\text{\,}\mathrm{K}$, i.e., the onset temperature of the magnetic
oscillations observed by $\mu^{+}$SR. However, in the presence of an intrinsic
magnetic transition, the SLR should instead show a peak at this temperature
due to critical magnetic fluctuations. This finding confirms once more that
the observed local magnetic field in $\mu^{+}$SR is not intrinsic to 2H-MoTe2.

## References

  * Guguchia et al. (2018) Z. Guguchia, A. Kerelsky, D. Edelberg, S. Banerjee, F. v. Rohr, D. Scullion, M. Augustin, M. Scully, D. A. Rhodes, Z. Shermadini, et al., Science Advances 4, eaat3672 (2018).
  * Tongay et al. (2012) S. SITongay, S. S. Varnoosfaderani, B. R. Appleton, J. Wu, and A. F. Hebard, Applied Physics Letters 101, 123105 (2012).
  * Strocov et al. (2010) V. N. Strocov, T. Schmitt, U. Flechsig, T. Schmidt, A. Imhof, Q. Chen, J. Raabe, R. Betemps, D. Zimoch, J. Krempasky, et al., Journal of Synchrotron Radiation 17, 631 (2010).
  * Strocov et al. (2014) V. N. Strocov, X. Wang, M. Shi, M. Kobayashi, J. Krempasky, C. Hess, T. Schmitt, and L. Patthey, Journal of Synchrotron Radiation 21, 32 (2014).
  * Strocov et al. (2012) V. N. Strocov, M. Shi, M. Kobayashi, C. Monney, X. Wang, J. Krempasky, T. Schmitt, L. Patthey, H. Berger, and P. Blaha, Physical Review Letters 109, 086401 (2012).
  * Böker et al. (2001) T. Böker, R. Severin, A. Müller, C. Janowitz, R. Manzke, D. Voß, P. Krüger, A. Mazur, and J. Pollmann, Physical Review B 64, 235305 (2001).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Physical Review B 54, 11169 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Physical Review B 59, 1758 (1999).
  * Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
  * Möller et al. (2013) J. S. Möller, D. Ceresoli, T. Lancaster, N. Marzari, and S. J. Blundell, Physical Review B 87, 121108(R) (2013).
  * Bernardini et al. (2013) F. Bernardini, P. Bonfà, S. Massidda, and R. De Renzi, Physical Review B 87, 115148 (2013).
  * Henkelman and Jónsson (2000) G. Henkelman and H. Jónsson, The Journal of Chemical Physics 113, 9978 (2000).
  * Ozaki (2003) T. SIOzaki, Physical Review B 67, 155108 (2003).
  * Momma and Izumi (2011) K. Momma and F. Izumi, Journal of Applied Crystallography 44, 1272 (2011).
  * Ma et al. (2011) Y. Ma, Y. Dai, M. Guo, C. Niu, J. Lu, and B. Huang, Physical Chemistry Chemical Physics 13, 15546 (2011).
  * González-Herrero et al. (2016) H. González-Herrero, J. M. Gómez-Rodríguez, P. Mallet, M. Moaied, J. J. Palacios, C. Salgado, M. M. Ugeda, J.-Y. Veuillen, F. Yndurain, and I. Brihuega, Science 352, 437 (2016).
  * Amato et al. (2017) A. Amato, H. Luetkens, K. Sedlak, A. Stoykov, R. Scheuermann, M. Elender, A. Raselli, and D. Graf, Review of Scientific Instruments 88, 093301 (2017).
  * Yaouanc and de Réotier (2011) A. Yaouanc and P. de Réotier, _Muon Spin Rotation, Relaxation, and Resonance: Applications to Condensed Matter_ , International Series of Monographs on Physics (OUP Oxford, 2011), ISBN 978-0-19-959647-8.
  * Suter and Wojek (2012) A. Suter and B. M. Wojek, Physics Procedia 30, 69 (2012).
  * Herak et al. (2013) M. Herak, A. Zorko, M. Pregelj, O. Zaharko, G. Posnjak, Z. Jagličić, A. Potočnik, H. Luetkens, J. Van Tol, A. Ozarowski, et al., Physical Review B 87, 104413 (2013).
  * Patterson (1988) B. D. Patterson, Reviews of Modern Physics 60, 69 (1988).
  * Cox (2003) S. F. J. Cox, Journal of Physics: Condensed Matter 15, R1727 (2003).
  * Prokscha et al. (2008) T. Prokscha, E. Morenzoni, K. Deiters, F. Foroughi, D. George, R. Kobler, A. Suter, and V. Vrankovic, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 595, 317 (2008).
  * Eckstein (1991) W. Eckstein, _Computer Simulation of Ion-Solid Interactions_ , Springer Series in Materials Science (Springer-Verlag, Berlin Heidelberg, 1991), ISBN 978-3-642-73515-8.
  * Harris et al. (2002) R. K. Harris, E. D. Becker, S. M. C. d. Menezes, R. Goodfellow, and P. Granger, Magnetic Resonance in Chemistry 40, 489 (2002).
  * Tian et al. (2020) Y. Tian, N. Ghassemi, and J. H. Ross, Physical Review B 102, 165149 (2020).
  * Tian et al. (2021) Y. Tian, Y. Zhu, R. Li, Z. Mao, and J. H. Ross, Physical Review B 104, L041105 (2021).