# Unbiasing Enhanced Sampling on a High-dimensional Free Energy Surface with
Deep Generative Model

Yikai Liu  [<EMAIL_ADDRESS>Tushar K. Ghosh  [<EMAIL_ADDRESS>Guang
Lin  [<EMAIL_ADDRESS>Ming Chen<EMAIL_ADDRESS>[

###### Abstract

Biased enhanced sampling methods utilizing collective variables (CVs) are
powerful tools for sampling conformational ensembles. Due to high intrinsic
dimensions, efficiently generating conformational ensembles for complex
systems requires enhanced sampling on high-dimensional free energy surfaces.
While methods like temperature-accelerated molecular dynamics (TAMD) can adopt
many CVs in a simulation, unbiasing the simulation requires accurate modeling
of a high-dimensional CV probability distribution, which is challenging for
traditional density estimation techniques. Here we propose an unbiasing method
based on the score-based diffusion model, a deep generative learning method
that excels in density estimation across complex data landscapes. We test the
score-based diffusion unbiasing method on TAMD simulations. The results
demonstrate that this unbiasing approach significantly outperforms traditional
unbiasing methods, and can generate accurate unbiased conformational ensembles
for simulations with a number of CVs higher than usual ranges.

###### keywords:

Enhanced Sampling, Machine Learning, Score-based Diffusion Model, Molecular
Dynamics Simulation

Purdue ME] Department of Mechanical Engineering, Purdue University, West
Lafayette, IN, 47906 Purdue Chemistry] Department of Chemistry, Purdue
University, West Lafayette, IN, 47906 Purdue ME] Department of Mechanical
Engineering, Purdue University, West Lafayette, IN, 47906 Purdue Chemistry]
Department of Chemistry, Purdue University, West Lafayette, IN, 47906
IR,NMR,UV

Molecular Dynamics (MD) simulations have emerged as a primary computational
tool for studying the thermodynamic and kinetic properties of complex systems
in chemistry, biology, and material science 1. With the assistance of
supercomputers 2, 3, it is now possible to perform milliseconds of all-atom MD
simulations for medium-sized proteins. To extend MD simulations to broader
time and length scales, multiple enhanced sampling methods have been designed
to increase MD sampling efficiency. 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 Among
all enhanced sampling methods, CV-based enhanced sampling methods focus on
several important degrees of freedom that capture systems’ essential dynamics.
By biasing the probability distribution along CVs, CV-based enhanced sampling
methods encourage systems to cross high energy barriers and explore different
regions of the energy landscape more efficiently.

A critical assumption behind CV-based enhanced sampling methods is the
manifold hypothesis 15, 16, 17, which posits that high-dimensional all-atom
configurations often lie along a low-dimensional latent manifold, and such a
low-dimensional manifold can accurately describe the important features of the
high-dimensional systems. Traditionally, physics-based CVs are chosen from
experimentally measurable properties, geometric descriptors and order
parameters with important underlying physics, such as end-to-end distance 18,
radius of gyration 19, or backbone torsion angles of proteins 20. Recently,
machine-learning-based methods that utilize dimensionality reduction
techniques have been applied to design CVs 21, 22, 23, 24, 25, 26, 27, 28, 29,
30, 31, 32. Regardless of CV categories, it is challenging to fully describe a
complex system in biochemistry and material science with one or two
parameters, a typical number of CVs used in many enhanced sampling
simulations. For example, following Two Nearest Neighbors method 33, we
estimate the intrinsic dimension (minimum numbers of parameters) to accurately
describe Amyloid-beta 42 (A$\beta_{42})$34 to be 7.13, as shown in Fig.(1). A
number of CVs below 7 will not fully describe this system.

Figure 1: (a) The Two Nearest Neighbors method calculates the ratio $\mu$
between the nearest and the second nearest neighbour distance for each data
point (red dots). $F(\mu)$ is the accumulative distribution of $\mu$. The
slope (black fitting line) of $-\log(1-F(\mu))$ as a function of $\log(\mu)$
is an estimation of the intrinsic dimension of A$\beta_{42}$ (see the
Supporting Information for details of the Two Nearest Neighbors method). (b)
The molecular structure of A$\beta_{42}$.

For most enhanced sampling methods 35, 36, 37, 38 that utilize biasing
potentials to assist crossing of energy barriers, it is challenging to use
many CVs, primarily due to the difficulty of constructing accurate biasing
potential for a high-dimensional FES. Other methods like driven-adiabatic free
energy dynamics/temperature accelerated molecular dynamics (TAMD)6, 7, 35, 39,
which enhance MD sampling by rising the temperature of certain degrees of
freedom, are capable of adopting many CVs in enhanced sampling simulations,
since biasing potential is not required. However, when handling many CVs, the
task of unbiasing TAMD trajectories to produce unbiased conformational
ensembles becomes increasingly complex. This complexity arises from the need
for accurate modeling of the high-dimensional, often multimodal CV probability
distribution. Traditional density estimation techniques, such as histogram
methods, kernel Density Estimation (KDE) 40, Nearest Neighbor Density
Estimation 41, and Gaussian Mixture Model (GMM) 42, often struggle to
accurately capture the nuances of such intricate distributions. For instance,
KDE suffers in high-dimensional spaces as it may produce overly smooth or
distorted estimates due to the lack of data across the expansive high-
dimensional space. On the other hand, GMM suffers from scalability,
initialization sensitivity, and a trivial model selection process. In this
paper, we leverage the score-based diffusion model (SBDM)43 for accurate
unbiasing of enhanced sampling simulations with many CVs. In our study, we
evaluate the performance of the SBDM-based unbiasing method in TAMD
simulations. We will demonstrate that SBDM excels in constructing CV
probability distributions, and can adapt to non-Euclidean CV such as torsion
angles, with minor changes to the model architecture 44. These capacities
endow the SBDM-based unbiasing method with superior performance and
versatility.

We will first introduce the TAMD method and its unbiasing formula. For a
system of $N$ particles, we denote its Cartesian coordinates by
$\mathbf{r}\equiv(\mathbf{r}_{1},\mathbf{r}_{2},\ldots,\mathbf{r}_{N})$, and
$n$ collective variables by
$\mathbf{q}\equiv(q_{1}(\mathbf{r}),\ldots,q_{n}(\mathbf{r}))$. In TAMD,
$\mathbf{q}$ are coupled with extended variables
$\mathbf{z}\equiv(z_{1},\ldots,z_{n})$ with stiff harmonic potentials
$\sum_{i}\kappa_{i}/2(q_{i}(\mathbf{r})-z_{i})^{2}$. $\mathbf{z}$ typically
shares the same topology as $\mathbf{q}$. It has been proved that the free
energy surface $A(\mathbf{q})$ can be approximated with the free energy
surface of extended variables $A_{\kappa}(\mathbf{z})$ when
$\kappa_{i}\rightarrow\infty$ for all $\kappa_{i}$. TAMD introduces a high
temperature $T_{h}\gg T$ for $\mathbf{z}$, and maintains $\mathbf{r}$ at
desired temperature $T$. To keep the thermodynamic properties of the system,
$\mathbf{z}$ are adiabatically decoupled from $\mathbf{r}$ by assigning each
$z_{i}$ a fictitious mass $\mu_{i}\gg 1$. By defining
$\beta_{h}=1/k_{\mathrm{B}}T_{h}$ where $k_{\mathrm{B}}$ is the Boltzmann
constant, the joint probability distribution of $\mathbf{r}$ and $\mathbf{z}$
from a TAMD simulation satisfies 35

$P_{\mathrm{TAMD}}(\mathbf{r},\mathbf{z})\approx
P_{T_{h}}(\mathbf{z})P(\mathbf{r}|\mathbf{z})\;\;,$ (1)

where $P_{T_{h}}(\mathbf{z})\propto\exp\\{-\beta_{h}A_{\kappa}(\mathbf{z})\\}$
is the marginal probability distribution of $\mathbf{z}$ and
$P(\mathbf{r}|\mathbf{z})$ is the Boltzmann distribution at the physical
temperature $T$ conditional on $\mathbf{z}$. Eq.(1) is exact if all
$\mu_{i}\rightarrow\infty$. The free energy of collective variables
$A_{\kappa}(\mathbf{z})$ can be easily obtained from TAMD with
$A(\mathbf{z})\approx-k_{\mathrm{B}}T_{h}\log P_{T_{h}}(\mathbf{z})$.

We are often interested in intuitive

Figure 2:  The free energy profiles generated from the unbiasing method of
using Eq.(2) are shown as solid lines with an “inaccurate” estimation of the
high-temperature probability. The blue, orange, and green lines correspond to
unbiasing results from $k_{\mathrm{B}}T_{h}=$3, 6, and 9, with
$k_{\mathrm{B}}T=1$ in all three cases. The black dashed line is the ground
truth with accurate estimation of the high-temperature probability.

features that are different from CVs used in an enhanced sampling to
understand properties of the simulated system. Unbiasing enhanced sampling
trajectories is necessary to project biased simulation data onto intuitive
features. Assuming $\mathbf{Y}(\mathbf{r})$ is a set of low dimensional
intuitive features of interests, the equilibrium probability of
$\mathbf{Y}(\mathbf{r})=\mathbf{y}$ can be written as

$P(\mathbf{y})=\int\delta(\mathbf{y}-\mathbf{Y}(\mathbf{r}))\omega(\mathbf{z})P_{\mathrm{TAMD}}(\mathbf{r},\mathbf{z})\mathrm{d}\mathbf{r}\mathrm{d}\mathbf{z},$
(2)

where $\omega(\mathbf{z})=P_{T_{h}}(\mathbf{z})^{T_{h}/T-1}$ is the unbiasing
weight. If $\mathbf{Y}$ can be written as a function of CVs, Eq.(2) can be
reduced to
$P(\mathbf{y})=\int\delta(\mathbf{y}-\mathbf{Y}(\mathbf{z}))\omega(\mathbf{z})P_{T_{h}}(\mathbf{z})\mathrm{d}\mathbf{z}$.

A good estimation of $P_{T_{h}}(\mathbf{z})$ is crucial for obtaining an
accurate $\omega(\mathbf{z})$ in TAMD. Errors in estimating
$P_{T_{h}}(\mathbf{z})$ magnify errors in $\omega(\mathbf{z})$ at a high
$T_{h}$, leading to inaccurate $P(\mathbf{y})$. We design a toy problem to
demonstrate the importance of accurately modeling $P_{T_{h}}(\mathbf{z})$ in
unbiasing TAMD. In this example, a one dimensional probability density
$P(\mathbf{z})$ at $k_{\mathrm{B}}T=1$ is constructed with a mixture of
Gaussians. The probability is

Figure 3: (a) The mechanism of denoising SBDM. During the diffusion process,
the data (in this demonstration a picture) is gradually perturbed to an
isotropic Gaussian noise via a forward SDE. (b) The flow chat of unbiasing
TAMD with SBDM. We perform TAMD simulation at a high temperature $T_{h}$, and
construct SBDM for $\mathbf{z}$. The time-dependent score function, squared in
orange, is used to perform density estimation of CV, $P(\mathbf{z})\equiv
P_{T_{h}}(\mathbf{z})$, from a simple distribution, as shown above the red
line. The weight $\omega(\mathbf{z})$ of configurations is evaluated from the
estimated CV probability $P(\mathbf{z})$. The weight is used to compute the
unbiased distribution of features of interest $P(\mathbf{y})$.

scaled to different $T_{h}$ followed by convolution with a Gaussian kernel to
“perturb” the high temperature probability $\tilde{P}_{T_{h}}(\mathbf{z})$.
$\tilde{P}_{T_{h}}(\mathbf{z})$ is used to calculate weight
$\omega(\mathbf{z})$. The “unbiased” distribution of $\mathbf{z}$,
$P_{\mathrm{ub}}(\mathbf{z})$, is then calculated with Eq.(2) (see the
Supporting Information for details of the toy problem). Fig.(2) demonstrates
that the error in $P_{\mathrm{ub}}(\mathbf{z})$ is more sensitive to errors in
$\tilde{P}_{T_{h}}(\mathbf{z})$ when $T_{h}$ is higher.

We then briefly review the SBDM framework and its usage in density estimation.
The SBDM perturbs data to noise prior with a diffusion process over a unit
time by a linear stochastic differential equation (SDE):

$d\mathbf{z}=f(\mathbf{z},t)dt+g(t)d\mathbf{w},\;t\in[0,1],$ (3)

where $f(\mathbf{z},t),g(t)$ are user-defined drift and diffusion functions of
the SDE and $\mathbf{w}$ denotes a standard Wiener process. In this paper, an
SDE with the drift term $f(\mathbf{z},t)=0$ is used. With carefully designed
$g(t)$, the marginal probability of $\mathbf{z}$ at diffusion time $t$,
$P_{t}(z)$, changes from the data distribution at $t=0$ to approximately a
simple Gaussian distribution at $t=1$.

For any diffusion process in Eq.(3), it has a corresponding reverse-time SDE
45 :

$d\mathbf{z}=[f(\mathbf{z},t)-g^{2}(t)\nabla_{\mathbf{z}}\log
P_{t}(\mathbf{z})]dt+g(t)d\bar{\mathbf{w}},$ (4)

with $\bar{\mathbf{w}}$ a standard Wiener process in the reverse-time. The
trajectories of the reverse SDE have the same marginal densities as the
forward SDE. Thus, the reverse-time SDE can gradually convert noise to data.
The SBDM parameterizes the time-dependent score function
$\nabla_{\mathbf{z}}\log P_{t}(\mathbf{z})$ in the reverse SDE with a neural
network $\mathbf{s}_{\boldsymbol{\theta}}(\mathbf{z}(t),t)$. To estimate
$\nabla_{\mathbf{z}}\log P_{t}(\mathbf{z})$, a time-dependent score-based
model $\mathbf{s}_{\boldsymbol{\theta}}(\mathbf{z}(t),t)$ can be trained via
minimizing a denoising score matching loss:

$\displaystyle M(\theta)$
$\displaystyle=\left\|\mathbf{s}_{\boldsymbol{\theta}}(\mathbf{z}(t),t)-\nabla_{\mathbf{z}(t)}\log
P(\mathbf{z}(t)|\mathbf{z}(0))\right\|_{2}^{2}$ (5) $\displaystyle J(\theta)$
$\displaystyle=\underset{\boldsymbol{\theta}}{\arg\min}\mathbb{E}_{t}\left\\{\mathbb{E}_{\mathbf{z}(0)}\mathbb{E}_{\mathbf{z}(t)\mid\mathbf{z}(0)}M(\theta)\right\\},$

with $t$ uniformly sampled between $[0,1]$. We note that
$\nabla_{\mathbf{z}(t)}\log P(\mathbf{z}(t)|\mathbf{z}(0))$ is not explicitly
required in the score matching loss while samples following
$P(\mathbf{z}(t)|\mathbf{z}(0))$ are needed.

Finally, SBDM defines a deterministic way to compute the data distribution
$P_{0}(\mathbf{z})$ as follows, with $f(\mathbf{z},t)=0$:

$\displaystyle\log P_{0}(\mathbf{z}(0))$ $\displaystyle=\log
P_{1}(\mathbf{z}(1))$
$\displaystyle-\frac{1}{2}\int_{0}^{1}g^{2}(t)\nabla_{\mathbf{z}}\cdot\mathbf{s}(\mathbf{z}(t),t)dt.$
(6)

CVs can be defined in spaces with different topologies. For example, $n$
torsion angles are defined on a hypertorus space $\mathbb{T}^{n}$, and
quaternions which representing rigid-body rotation are defined on a 3D unit
sphere $\mathbb{S}^{3}$. Therefore, extending the SBDM to different topologies
is important. In this study we will focus on SBDM on a hypertorus space. The
theory behind SBDM holds for compact Riemannian manifolds, with subtle
modifications. For $\mathbf{z}\in M$, a Riemannian manifold (such as
hypertorus $\mathbb{T}^{n}$), $\mathbf{w}$ being the Brownian motion on the
manifold, and $f(\mathbf{z},t)\in T_{\mathbf{z}}M$, a tangent space, Eq.(4)
still holds 44.

As shown in Eq.(5), training a denoising score matching model requires
sampling from the perturbation kernel $P(\mathbf{z}(t)|\mathbf{z}(0))$ of the
forward diffusion defined by Eq.(3). We consider the perturbation kernel on
$\mathbb{T}^{n}$ with wrapped normal distribution:

$\displaystyle\mathbf{U}_{\mathbf{d}}(t)=\frac{\left\|\boldsymbol{\mathbf{z}}(0)-\boldsymbol{\mathbf{z}}(t)+2\pi\mathbf{d}\right\|^{2}}{2\sigma^{2}(t)},\;\mathbf{d}\in\mathbb{Z}^{n}$
(7) $\displaystyle
P(\mathbf{z}(t)|\mathbf{z}(0))\propto\sum_{\mathbf{d}}\exp\left(-\mathbf{U}_{\mathbf{d}}(t)\right),$

where $g(t)$ and $\sigma(t)$ are related by
$g(t)=\sqrt{\mathrm{d}\sigma^{2}(t)/\mathrm{d}t}$. The rest of the terms in
the loss function in Eq. (5) remain the same.

Figure 4: Molecular structure of alanine dipeptide is shown in (a). (b, c) FES
w.r.t backbone dihedral angles of alanine dipeptide are presented, obtained by
metadynamics (b) and SBDM-TAMD (c). The red/black points represent the
location of important minima and saddle points on the FES. (d-f) The absolute
free energy difference between metadynamics and SBDM-TAMD, KDE-TAMD, and GMM-
TAMD.

Note that density estimation in Eq.(6) is applicable to both variables in
Euclidean and hypertorus space. The architecture of an SBDM model is highly
flexible. For example, an SBDM model can use Residual neural network (ResNet)
46, U-Net 47, Graph Neural Network (GNN) 48, etc.

In the following section, we will demonstrate how incorporating SBDM can
fulfill the strict density estimation accuracy requirement of unbiasing TAMD,
thus allowing TAMD to generate correct unbiased ensembles. We tested the
efficiency and accuracy of SBDM to unbiasing TAMD simulations on three
systems: (1) alanine dipeptide (2) glutamine dipeptide, and (3) met-enkephalin
(see the Supporting Information for simulation details). In all three systems,
we conducted TAMD with torsion angles as collective variables, with
$T_{h}=1200K$ in the first example while $T_{h}=900K$ in the second and third
examples. The physical variables were maintained at temperature $T=300K$ in
all three experiments. We then unbiased TAMD with density estimation performed
by SBDM on hypertorus space (see the Supporting Information for details of
training SBDM models). Unbiasing TAMD with SBDM as density estimation method
is referred as SBDM-TAMD in the rest of the paper. The general framework of
SBDM-TAMD is summarized in Fig.(3). For a fair comparison, we unbiased the
same TAMD simulations with kernel density estimation (KDE-TAMD) and Gaussian
Mixture Model (GMM-TAMD) in all three systems, and Normalizing Flow on
hypertorus space (NF-TAMD) 49, 50 in the last system. We performed converged
metadynamics simulations at $300K$ as the baseline results of all three
systems. We will compare the unbiasing results from different density
estimation methods to the baseline results.

The first proof-of-concept example is an alanine dipeptide in the aqueous
solution with implicit solvent. This has been a benchmark system with well-
established FES in previous studies. 51, 52 In both TAMD and the baseline
metadynamics, we used the Ramachandran angles ($\phi$, $\psi$) as CVs. The
accuracy of SBDM-TAMD is demonstrated in Fig.(4). The SBDM capably captures
the intricate free energy delineated by metadynamics, as shown in Fig.(4)(c)
and Fig.(4)(d). Notably, the model accurately locates all six free energy
minima, with depth of the minima quantitatively accurate (<1 kJ/mol energy
difference). Furthermore, the saddle points - representative of transitional
between conformational basins, are also well reproduced. Energy errors in low
free energy saddle points are less than 2 kJ/mol compared to metadynamics
results. This match implies that the SBDM-TAMD is accurate for studying

Figure 5: Molecular structure of Glutamine dipeptide is shown in (a). FES of
backbone dihedral angles of Glutamine dipeptide are calculated by metadnyamics
(b) and SBDM-TAMD (C). The red/black points represent the location of
important minima and saddle points on the FES. (d)-(f) shows he absolute free
energy difference between metadynamics and SBDM-TAMD, KDE-TAMD, and GMM-TAMD.

both thermodyanmics and kinetics of alanine dipeptide conformational changes.
We notice that for this two-dimensional problem, KDE and GMM perform
relatively well, but still with less accuracy compared to SBDM. Both methods
can locate the free energy minima but with a larger energy difference of up to
2.5 kJ/mol from the baseline metadynamics. The free energy of saddle points
have even larger deviation, up to 3.5 kJ/mol. The clear disparities in free
energy differences, as shown in Fig. 4(e) and Fig. 4(f), especially in regions
corresponding to the annotated energy minima and saddle points, illustrate the
inadequacies of KDE and GMM for precise FES calculations.

The second system we studied is Glutamine dipeptide in the aqueous solution
with explicit solvent. A TAMD simulation was performed with five dihedral
angles on the backbone and the side chain
($\phi,\psi,\chi_{1},\chi_{2},\chi_{3}$) as CVs to enhance the sampling of
both backbone and side chain conformations. We also run a benchmark
metadynamics with two backbone diheral angles as CVs. The result is
demonstrated in Fig.(5). Projecting unbiased SBDM-TAMD trajectory onto
backbone dihedral angles quantitavely matches the bechmark FES from
metadynamics. Energy errors of minima on the projected 2D FES are within 2
kJ/mol while errors of low-free-energy saddles are approximately 5 kJ/mol. As
a comparison, both KDE and GMM struggle to uphold their precision with the
increased CV dimensions. Both of these two methods introduce larger errors at
$\alpha$ and $\alpha_{L}$. This example highlights the SBDM’s superior
adaptability and accuracy in high-dimensional CV compared to traditional
methods.

The final, more challenging system we studied is the oligopeptide met-
enkephalin (Tyr-Gly-Gly-Phe-Met) in the aqueous solution with explicit water,
which is a common test case for enhanced sampling techniques. 53, 35, 54 For
TAMD, we chose ten backbone dihedral angles
($\phi_{1},\psi_{1},...,\phi_{5},\psi_{5}$) as CVs. The baseline metadynamics
simulation was performed with a 2D stochastic kinetic embedding (StKE), a
manifold learning method which serves as a low-dimensional CV representation
that preserves kinetic information 25. We projected the unbiased SBDM-TAMD
trajectory to features like end-to-end distance $d_{\mathrm{ee}}$ and the
radius of gyration $R_{\mathrm{g}}$. As a comparison, we unbiased the
benchmark metadynamics simulation 55, 56, 57 and projected configurations from
metadynamics onto the same features. We want to emphasize that generating an
optimal machine-learning-based CVs for metadynamics (StKE in this work) is
non-trivial and requires extra simulation data. From Fig.(6), the SBDM-TAMD
exhibits high performance, demonstrating a FES that closely aligned with the
metadynamics baseline. The minima and barriers on each projected one-
dimensional FES predicted by SBDM is in

Figure 6: (a) Molecular structure of met-enkephalin. (b-c) Free energy of end-
to-end distance and radius of gyration, with red dots indicating stable and
meta-stable states.

good agreements with benchmark, suggesting SBDM-TAMD is capable of modeling
thermodynamics and kinetics of a polypeptide. For example, as shown in Fig.
6(b), both metadynamics and SBDM-TAMD identify conformations with small
$d_{\mathrm{ee}}$ as a meta-stable conformation with consistent meta-
stability, while KDE and GMM fail to predict this meta-stable conformation.
Also, KDE and GMM fail to predict the location of the global minimum on the
FES of $d_{\mathrm{ee}}$ with an error of $\sim$3.0 Å. Although NF correctly
predicts the location of the global minimum, it overestimates the stability of
this minimum. Similarly, KDE ,GMM, and even NF, predict the most possible
$R_{\mathrm{g}}\approx 5$ Å, while SBDM-TAMD predicts that the most stable
conformations have $R_{\mathrm{g}}<4$ Å, which is consistent with the
metadynamics’ result.

In this paper, we developed an unbiasing method based on score-based diffusion
model, a deep generative learning model, to generate unbiased conformational
ensembles from collective-variable-based biased enhanced sampling simulations.
Our method can adapt to simulations with collective variables of large amounts
and different topologies. We test the unbiasing method on temperature-
accelerated molecular dynamics, an enhanced sampling method that can utilize
many collective variables to efficiently explore a high-dimensional free
energy surface. Numerical experiments across three systems of 2,5 and 10
collective variables underscore our unbiasing method’s exceptional accuracy
and adaptability to high-dimensional collective variables. Although
biomolecules were used in the numerical experiments, the developed approach
can be used to other problems like modeling material phase transitions.
Looking ahead, we aim to explore the method’s potential in unbiasing
simulations with even larger amounts of collective variables, and form
accurate high-dimensional biasing potentials for wider ranges of enhanced
sampling methods.

Y.L., T.K.G. and M.C. acknowledge support from Purdue startup funding. G.L.
and Y.L. gratefully acknowledge the support of the National Science Foundation
(DMS-2053746, DMS-2134209, ECCS-2328241, and OAC-2311848), and U.S. Department
of Energy (DOE) Office of Science Advanced Scientific Computing Research
program DE-SC0021142, DE-SC0023161, and the Uncertainty Quantification for
Multifidelity Operator Learning (MOLUcQ) project (Project No. 81739), and DOE
– Fusion Energy Science, under grant number: DE-SC0024583.

Supporting information contains details of MD simulations, SBDM training, Two
Nearest Neighbors intrinsic dimension estimation, and toy model.

## References

  * Ciccotti et al. 2014 Ciccotti, G.; Ferrario, M.; Schuette, C., et al. Molecular dynamics simulation. _Entropy_ 2014, _16_ , 1
  * Shaw et al. 2014 Shaw, D. E.; Grossman, J.; Bank, J. A.; Batson, B.; Butts, J. A.; Chao, J. C.; Deneroff, M. M.; Dror, R. O.; Even, A.; Fenton, C. H. et al. Anton 2: raising the bar for performance and programmability in a special-purpose molecular dynamics supercomputer. SC’14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. 2014; pp 41–53
  * Shaw et al. 2021 Shaw, D. E.; Adams, P. J.; Azaria, A.; Bank, J. A.; Batson, B.; Bell, A.; Bergdorf, M.; Bhatt, J.; Butts, J. A.; Correia, T. et al. Anton 3: twenty microseconds of molecular dynamics simulation before lunch. Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. 2021; pp 1–11
  * 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
  * Maragliano and Vanden-Eijnden 2006 Maragliano, L.; Vanden-Eijnden, E. A temperature accelerated method for sampling free energy and determining reaction pathways in rare events simulations. _Chemical physics letters_ 2006, _426_ , 168–175
  * Abrams and Tuckerman 2008 Abrams, J. B.; Tuckerman, M. E. Efficient and direct generation of multidimensional free energy surfaces via adiabatic dynamics without coordinate transformations. _The Journal of Physical Chemistry B_ 2008, _112_ , 15742–15757
  * Darve et al. 2008 Darve, E.; Rodríguez-Gómez, D.; Pohorille, A. Adaptive biasing force method for scalar and vector free energy calculations. _The Journal of chemical physics_ 2008, _128_ , 144120
  * Barducci et al. 2011 Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2011, _1_ , 826–843
  * Abrams and Vanden-Eijnden 2010 Abrams, C. F.; Vanden-Eijnden, E. Large-scale conformational sampling of proteins using temperature-accelerated molecular dynamics. _Biophysical Journal_ 2010, _98_ , 26a
  * Yu and Tuckerman 2011 Yu, T.-Q.; Tuckerman, M. E. Temperature-accelerated method for exploring polymorphism in molecular crystals based on free energy. _Physical review letters_ 2011, _107_ , 015701
  * Yu et al. 2014 Yu, T.-Q.; Chen, P.-Y.; Chen, M.; Samanta, A.; Vanden-Eijnden, E.; Tuckerman, M. Order-parameter-aided temperature-accelerated sampling for the exploration of crystal polymorphism and solid-liquid phase transitions. _The Journal of chemical physics_ 2014, _140_ , 214109
  * Samanta et al. 2014 Samanta, A.; Tuckerman, M. E.; Yu, T.-Q.; E, W. Microscopic mechanisms of equilibrium melting of a solid. _Science_ 2014, _346_ , 729–732
  * Tiwary and Parrinello 2013 Tiwary, P.; Parrinello, M. From metadynamics to dynamics. _Physical review letters_ 2013, _111_ , 230602
  * Fefferman et al. 2016 Fefferman, C.; Mitter, S.; Narayanan, H. Testing the manifold hypothesis. _Journal of the American Mathematical Society_ 2016, _29_ , 983–1049
  * Narayanan and Mitter 2010 Narayanan, H.; Mitter, S. Sample complexity of testing the manifold hypothesis. _Advances in neural information processing systems_ 2010, _23_
  * Gorban and Tyukin 2018 Gorban, A. N.; Tyukin, I. Y. Blessing of dimensionality: mathematical foundations of the statistical physics of data. _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 2018, _376_ , 20170237
  * Hummer and Szabo 2001 Hummer, G.; Szabo, A. Free energy reconstruction from nonequilibrium single-molecule pulling experiments. _Proceedings of the National Academy of Sciences_ 2001, _98_ , 3658–3661
  * Bussi et al. 2006 Bussi, G.; Gervasio, F. L.; Laio, A.; Parrinello, M. Free-energy landscape for $\beta$ hairpin folding from combined parallel tempering and metadynamics. _Journal of the American Chemical Society_ 2006, _128_ , 13435–13441
  * Granata et al. 2013 Granata, D.; Camilloni, C.; Vendruscolo, M.; Laio, A. Characterization of the free-energy landscapes of proteins by NMR-guided metadynamics. _Proceedings of the National Academy of Sciences_ 2013, _110_ , 6817–6822
  * Das et al. 2006 Das, P.; Moll, M.; Stamati, H.; Kavraki, L. E.; Clementi, C. Low-dimensional, free-energy landscapes of protein-folding reactions by nonlinear dimensionality reduction. _Proceedings of the National Academy of Sciences_ 2006, _103_ , 9885–9890
  * Ferguson et al. 2011 Ferguson, A. L.; Panagiotopoulos, A. Z.; Debenedetti, P. G.; Kevrekidis, I. G. Integrating diffusion maps with umbrella sampling: Application to alanine dipeptide. _The Journal of chemical physics_ 2011, _134_ , 135103
  * Rohrdanz et al. 2011 Rohrdanz, M. A.; Zheng, W.; Maggioni, M.; Clementi, C. Determination of reaction coordinates via locally scaled diffusion map. _The Journal of chemical physics_ 2011, _134_ , 124116
  * Noé and Clementi 2015 Noé, F.; Clementi, C. Kinetic distance and kinetic maps from molecular dynamics simulation. _Journal of chemical theory and computation_ 2015, _11_ , 5002–5011
  * Zhang and Chen 2018 Zhang, J.; Chen, M. Unfolding hidden barriers by active enhanced sampling. _Physical review letters_ 2018, _121_ , 010601
  * Fu et al. 2022 Fu, X.; Xie, T.; Rebello, N. J.; Olsen, B. D.; Jaakkola, T. Simulate time-integrated coarse-grained molecular dynamics with geometric machine learning. _arXiv preprint arXiv:2204.10348_ 2022,
  * Ceriotti et al. 2011 Ceriotti, M.; Tribello, G. A.; Parrinello, M. Simplifying the representation of complex free-energy landscapes using sketch-map. _Proceedings of the National Academy of Sciences_ 2011, _108_ , 13023–13028
  * Tiwary and Berne 2016 Tiwary, P.; Berne, B. Spectral gap optimization of order parameters for sampling complex molecular systems. _Proceedings of the National Academy of Sciences_ 2016, _113_ , 2839–2844
  * 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
  * Ribeiro et al. 2018 Ribeiro, J. M. L.; Bravo, P.; Wang, Y.; Tiwary, P. Reweighted autoencoded variational Bayes for enhanced sampling (RAVE). _The Journal of chemical physics_ 2018, _149_ , 072301
  * Wang et al. 2019 Wang, Y.; Ribeiro, J. M. L.; Tiwary, P. Past–future information bottleneck for sampling molecular reaction coordinate simultaneously with thermodynamics and kinetics. _Nature communications_ 2019, _10_ , 3573
  * Nuske et al. 2014 Nuske, F.; Keller, B. G.; Pérez-Hernández, G.; Mey, A. S.; Noé, F. Variational approach to molecular kinetics. _Journal of chemical theory and computation_ 2014, _10_ , 1739–1752
  * Facco et al. 2017 Facco, E.; d’Errico, M.; Rodriguez, A.; Laio, A. Estimating the intrinsic dimension of datasets by a minimal neighborhood information. _Scientific reports_ 2017, _7_ , 12140
  * Findeis 2007 Findeis, M. A. The role of amyloid $\beta$ peptide 42 in Alzheimer’s disease. _Pharmacology & therapeutics_ 2007, _116_ , 266–286
  * Chen et al. 2012 Chen, M.; Cuendet, M. A.; Tuckerman, M. E. Heating and flooding: A unified approach for rapid generation of free energy surfaces. _The Journal of Chemical Physics_ 2012, _137_ , 024102
  * Debnath and Parrinello 2020 Debnath, J.; Parrinello, M. Gaussian Mixture-Based Enhanced Sampling for Statics and Dynamics. _The Journal of Physical Chemistry Letters_ 2020, _11_ , 5076–5080, PMID: 32510225
  * Zhang et al. 2018 Zhang, L.; Wang, H.; E, W. Reinforced dynamics for enhanced sampling in large atomic and molecular systems. _The Journal of Chemical Physics_ 2018, _148_ , 124113
  * Zhang et al. 2020 Zhang, J.; Lei, Y.-K.; Yang, Y. I.; Gao, Y. Q. Deep learning for variational multiscale molecular modeling. _The Journal of Chemical Physics_ 2020, _153_ , 174115
  * Cuendet and Tuckerman 2014 Cuendet, M. A.; Tuckerman, M. E. Free energy reconstruction from metadynamics or adiabatic free energy dynamics simulations. _Journal of chemical theory and computation_ 2014, _10_ , 2975–2986
  * Chen 2017 Chen, Y.-C. A tutorial on kernel density estimation and recent advances. _Biostatistics & Epidemiology_ 2017, _1_ , 161–187
  * Mack and Rosenblatt 1979 Mack, Y.; Rosenblatt, M. Multivariate k-nearest neighbor density estimates. _Journal of Multivariate Analysis_ 1979, _9_ , 1–15
  * Reynolds et al. 2009 Reynolds, D. A., et al. Gaussian mixture models. _Encyclopedia of biometrics_ 2009, _741_
  * Song et al. 2021 Song, Y.; Sohl-Dickstein, J.; Kingma, D. P.; Kumar, A.; Ermon, S.; Poole, B. Score-Based Generative Modeling through Stochastic Differential Equations. International Conference on Learning Representations. 2021
  * De Bortoli et al. 2022 De Bortoli, V.; Mathieu, E.; Hutchinson, M.; Thornton, J.; Teh, Y. W.; Doucet, A. Riemannian score-based generative modelling. _Advances in Neural Information Processing Systems_ 2022, _35_ , 2406–2422
  * Anderson 1982 Anderson, B. D. Reverse-time diffusion equation models. _Stochastic Processes and their Applications_ 1982, _12_ , 313–326
  * He et al. 2016 He, K.; Zhang, X.; Ren, S.; Sun, J. Deep residual learning for image recognition. Proceedings of the IEEE conference on computer vision and pattern recognition. 2016; pp 770–778
  * Ronneberger et al. 2015 Ronneberger, O.; Fischer, P.; Brox, T. U-Net: Convolutional Networks for Biomedical Image Segmentation. Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015. Cham, 2015; pp 234–241
  * Scarselli et al. 2008 Scarselli, F.; Gori, M.; Tsoi, A. C.; Hagenbuchner, M.; Monfardini, G. The graph neural network model. _IEEE transactions on neural networks_ 2008, _20_ , 61–80
  * Rezende et al. 2020 Rezende, D. J.; Papamakarios, G.; Racaniere, S.; Albergo, M.; Kanwar, G.; Shanahan, P.; Cranmer, K. Normalizing flows on tori and spheres. International Conference on Machine Learning. 2020; pp 8083–8092
  * Köhler et al. 2021 Köhler, J.; Krämer, A.; Noé, F. Smooth normalizing flows. _Advances in Neural Information Processing Systems_ 2021, _34_ , 2796–2809
  * Strodel and Wales 2008 Strodel, B.; Wales, D. J. Free energy surfaces from an extended harmonic superposition approach and kinetics for alanine dipeptide. _Chemical Physics Letters_ 2008, _466_ , 105–115
  * Smith 1999 Smith, P. E. The alanine dipeptide free energy surface in solution. _The Journal of chemical physics_ 1999, _111_ , 5568–5579
  * Hénin et al. 2010 Hénin, J.; Fiorin, G.; Chipot, C.; Klein, M. L. Exploring Multidimensional Free Energy Landscapes Using Time-Dependent Biases on Collective Variables. _Journal of Chemical Theory and Computation_ 2010, _6_ , 35–47, PMID: 26614317
  * Chen et al. 2015 Chen, M.; Yu, T.-Q.; Tuckerman, M. E. Locating landmarks on high-dimensional free energy surfaces. _Proceedings of the National Academy of Sciences_ 2015, _112_ , 3235–3240
  * Sutto et al. 2010 Sutto, L.; D’Abramo, M.; Gervasio, F. L. Comparing the efficiency of biased and unbiased molecular dynamics in reconstructing the free energy landscape of met-enkephalin. _Journal of Chemical Theory and Computation_ 2010, _6_ , 3640–3646
  * Sutto et al. 2012 Sutto, L.; Marsili, S.; Gervasio, F. L. New advances in metadynamics. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2012, _2_ , 771–779
  * Sicard and Senet 2013 Sicard, F.; Senet, P. Reconstructing the free-energy landscape of Met-enkephalin using dihedral principal component analysis and well-tempered metadynamics. _The Journal of Chemical Physics_ 2013, _138_ , 235101


# Anisotropy of Yu-Shiba-Rusinov states in NbSe2

Mateo Uldemolins<EMAIL_ADDRESS>Donostia International Physics
Center (DIPC), 20018, Donostia–San Sebastián, Spain  Université Paris-Saclay,
CNRS, Laboratoire de Physique des Solides, 91405, Orsay, France  Freek Massee
Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides, 91405,
Orsay, France  Tristan Cren  Sorbonne Université, CNRS, Institut des
Nanosciences de Paris, UMR7588, F-75252 Paris, France  Andrej Mesaros
Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides, 91405,
Orsay, France  Pascal Simon  Université Paris-Saclay, CNRS, Laboratoire de
Physique des Solides, 91405, Orsay, France

###### Abstract

The spatial structure of in-gap Yu-Shiba-Rusinov (YSR) bound states induced by
a magnetic impurity in a superconductor is the essential ingredient for the
possibility of engineering collective impurity states. Recently, a saddle-
point approximation [Phys. Rev. B 105, 144503] revealed how the spatial form
of a YSR state is controlled by an anisotropic exponential decay length, and
an anisotropic prefactor, which depends on the Fermi velocity and Fermi-
surface curvature. Here we analyze STM data on YSR states in NbSe2, focusing
on a key issue that the exponential decay length predicted theoretically from
the small superconducting gap is much larger than the observed extent of YSR
states. We confirm that the exponential decay can be neglected in the analysis
of the anisotropy. Instead, we extract the anisotropic prefactor directly from
the data, matching it to the theoretical prediction, and we establish that the
theoretical expression for the prefactor alone captures the characteristic
flower-like shape of the YSR state. Surprisingly, we find that up to linear
order in the superconducting gap the anisotropic prefactor that determines the
shape of YSR states is the same as the anisotropic response to the impurity in
the underlying normal metal. Our work points out the correct way to analyze
STM data on impurities in small-gap superconductors, and reveals the
importance of the normal band structure’s curvature and Fermi velocity in
designing multi-impurity in-gap states in superconductors.

##  I Introduction

The local response to single magnetic impurities in conventional
superconductors appears most directly in the form of pairs of in-gap
excitations localized on the impurity known as Yu-Shiba-Rusinov (YSR) states
[1, 2, 3]. YSR states are nowadays routinely probed by scanning tunneling
microscopy (STM) experiments [4], and while an old problem, they have
attracted a renewed wave of interest owing to their potential to extract
information on the properties of the host [5, 6, 7], and study magnetism [8,
9] and spin-transport phenomena [10, 11, 12, 13] at the atomic scale. At the
same time, much experimental and theoretical effort has been devoted to
designing complex structures in real space (e.g. impurity chains or
lattices)[14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27] which are
suggested to realize topological states of matter harboring Majorana zero
modes [28]. Understanding the spatial structure of YSR states is key for
developing this program.

In an isotropic $d$-dimensional system, the local density of states (LDOS) of
YSR states far away from the impurity behaves according to $\sim
e^{-r/\xi}/r^{d-1}$, with the superconducting coherence length controlling the
exponential decay,
$\xi\sim\frac{v_{\mathbf{F}}\hbar}{\Delta}\sim\xi_{\mathrm{SC}}$. The
dependence of the power-law on the dimensionality of the substrate has been
evidenced by STM measurements on two-dimensional superconductors ($d=2$), e.g.
NbSe2, where the in-gap LDOS extends over several nanometers [29, 30, 31],
thereby rendering this class of materials a promising platform for
constructing collective impurity states [32, 33, 34, 35, 36]. In general, the
measured in-gap LDOS oftentimes displays a marked anisotropy on length scales
largely surpassing the typical impurity size, and hence attributed to the
lower symmetry of the band structure [30, 37]. The theory of spatial structure
of YSR states is hence intricate, requiring a blend of information about the
substrate’s band structure, pairing function, and impurity coupling.

Recently, the theory of the quasiparticle focusing effect (QFE) [38, 39] was
extended to superconducting substrates by some of the authors [40], leading to
the following insights in the limit of weak pairing,
$\Delta/\varepsilon_{\mathrm{F}}\ll 1$, and at long distances
($rk_{\mathrm{F}}\gg 1$):

  1. 1.

The QFE persists, i.e. the response to the impurity in a given real-space
direction is controlled by momenta (denoted critical) on the Fermi surface at
which the Fermi velocity lies parallel or antiparallel to that real-space
direction. This selectivity in reciprocal space controls the spatial
anisotropy of the YSR states; in particular, YSR states are focused in
directions perpendicular to flatter sections of the Fermi surface.

  2. 2.

Specifically to gapped superconductors, the in-gap YSR LDOS has an additional
exponentially decaying factor with angularly dependent exponential decay
length, apart from a prefactor and oscillations, all related to the geometry
of the Fermi surface through a simple analytical expression. Importantly, the
prefactor and exponential decay length are essentially always in phase, i.e.
they agree on the direction along which the LDOS is focused.

This theory describes very accurately the anisotropy of the characteristic
exponential decaying length of YSR states simulated on a lattice [40] and
successfully predicts the main anisotropy directions of YSR states observed in
NbSe2. However, the tight-binding description of the Nb-derived bands at the
Fermi level [41] in combination with the small superconducting gap
($\Delta\sim 1$ meV [42, 43]), yields an in-gap LDOS with large exponential
decay length ($\xi_{\mathrm{TB}}\sim 100$ nm) [29, 44, 35, 36] that widely
surpasses the extent of the observed YSR states, as well as the
superconducting coherence length reported in this material
($\xi^{\mathrm{exp}}_{\mathrm{SC}}\sim 10$ nm [45, 46, 47]). Therefore, there
is a problem in the natural interpretation that the anisotropy of the
coherence length in NbSe2 fundamentally controls the shape of the YSR LDOS,
and theoretical radial fitting of an exponential decay is questionable.

In this work, we show that the angular dependence of the LDOS prefactor
accounts for the anisotropy of YSR states in this material, while the radial
fitting of the LDOS can be circumvented, i.e., the (anisotropic) exponential
decay can be neglected. This evokes the QFE in normal metals [38, 39], where
standard Friedel oscillations decay algebraically and lack the exponential
factor. In fact, our analytical expression of the prefactor has the same
functional form, in terms of the Fermi velocity and curvature of the Fermi
surface, as the LDOS in case of a normal metal. A difference in principle
arises when this expression is evaluated at the critical momenta, which vary
as the band structure of the metal changes into the one of the superconductor.
Surprisingly, we find that up to linear order in
$\Delta/\varepsilon_{\mathrm{F}}\ll 1$ the prefactor of the impurity-induced
LDOS in the superconductor exactly equals the one in the parent normal metal.
We show that this prefactor can be reliably extracted from STM data on NbSe2,
and that the characteristic flower-shaped YSR LDOS is faithfully modelled by
the prefactor alone.

The rest of the paper is organized as follows: we start by presenting the STM
data in Sec. II. Then, in Sec. II.1, we briefly recall the main analytical
expressions derived in [40] and we analyze the spatial anisotropy of the LDOS
prefactor. This analysis is complemented in Sec. II.2, where we study the
Fourier transform of the YSR states. We conclude with a short discussion in
Sec. III. Additional details concerning the experiment and the theory
expressions are presented in the Appendices.

##  II Results

Figure 1: Spatial data sets 1 and 2. (a) Atomically-resolved constant-current
image of the Se surface of 2H-NbSe2, $V=4$ mV, $I=160$ pA. (b) In-gap current
recorded at $V=-0.6$ mV showing a six-fold YSR state. The (c,d) The same as
(a,b) for another YSR state, $V=4.2$ mV, $I=100$ pA. The lines at the bottom-
left corner indicate the crystallographic axes of 2H–NbSe2. Their length
corresponds to 2 nm in both data sets.

We consider STM data on 2H-NbSe2 crystals containing a few tens of ppm of
magnetic atoms measured at 300 mK with a metallic tip (for additional details
concerning the fabrication and measurement conditions, see Appendix A). YSR
bound states manifest as an in-gap peak in the differential-conductance
spectrum at the impurity location (see Appendix A). In Fig. 1, we present
atomically-resolved data around two independent impurities on the same sample.
The constant-current image [Fig. 1 (a) and (c)] shows the triangular geometry
of the Se lattice without any visible imperfection, implying that the defect
giving rise to the YSR state is most likely an atom substitution in the Nb
layer. The corresponding current images recorded before the onset of the
superconducting coherence peak [Fig. 1 (b) and (d)] show a marked six-fold
anisotropy that constitutes an excellent example of the QFE on a
superconductor.

###  II.1 Characterization of the LDOS prefactor in real space

Let us begin by briefly recalling the theoretical framework to describe the
QFE on superconductors. We model the impurity-substrate system as a classical
spin on an $s$-wave superconductor.

The Hamiltonian reads,

$\displaystyle H=H_{0}+H_{\mathrm{imp}},$ (1a) $\displaystyle
H_{0}=\sum_{\bm{k}\sigma}\varepsilon_{\bm{k}\sigma}c^{\dagger}_{\bm{k}\sigma}c_{\bm{k}\sigma}+\Delta\sum_{\bm{k}}\;c^{\dagger}_{\bm{k}\uparrow}c^{\dagger}_{\bm{-k}\downarrow}+\mathrm{h.c.},$
(1b) $\displaystyle
H_{\mathrm{imp}}=-J\left(c^{\dagger}_{\bm{r}_{0}\uparrow}c_{\bm{r}_{0}\uparrow}-c^{\dagger}_{\bm{r}_{0}\downarrow}c_{\bm{r}_{0}\downarrow}\right),$
(1c)

where $\varepsilon_{\bm{k}}$ is the spin-degenerate normal energy dispersion
of the substrate, $\Delta$ is the BCS superconducting gap and $J$ is the
amplitude of the magnetic exchange coupling between the superconducting
electrons and the spin of the impurity at $\bm{r}_{0}$. Specifically,
$\varepsilon_{\bm{k}}$ is an effective fifth-nearest neighbours tight-binding
energy dispersion on a triangular lattice that describes one of the two
lowest-lying Nb 4d bands of 2H-NbSe2 [41] (see Appendix C). This description
neglects interlayer coupling of the bulk crystal and assumes that the magnetic
impurity couples to one band only [29, 44]. It also neglects spin-orbit
coupling [48], however, it provides a faithful description of the geometry of
the Fermi surface and band structure of the substrate, which is the crucial
element to characterize the QFE. In addition, Eq. (1c) assumes a fully local
and isotropic scatterer, thereby neglecting the magnetic ion’s $d$-orbital
structure and the adsorption site’s symmetry. These add to the YSR-state
wavefunction some degree of anisotropy [8, 49]; however, the observed six-fold
anisotropy of the LDOS [Fig. 1 (b) and (d)] occurs on length scales of several
nanometers that largely surpass the typical orbital size, therefore it can be
safely ascribed to the QFE, and justifies the model’s simplification.

The LDOS at a distance $\bm{r}$ from the impurity and at the bound-state
energy $E_{\mathrm{S}}$ is computed by means of a saddle-point approximation
on the bare propagator valid in the far-field limit ($rk_{F,\mathrm{min}}\gg
1$, with $k_{F,\mathrm{min}}$ the minimum Fermi wave vector of the Fermi
contour). By further taking the small-gap limit
($\Delta/\varepsilon_{\mathrm{F}}\ll 1$), the LDOS takes the approximate form

$\rho(\bm{r})\sim\frac{1}{r}\sum_{j,j^{\prime}}\Gamma_{j,j^{\prime}}(\theta_{\bm{r}})e^{-r/\xi_{j,j^{\prime}}(\theta_{\bm{r}})}f[\bm{k}^{\pm}_{j,j^{\prime}}(\theta_{\bm{r}})\cdot\bm{r}],$
(2)

where the angular-dependent prefactor reads

$\Gamma_{j,j^{\prime}}(\theta_{\bm{r}})=\frac{1}{|\nabla\varepsilon_{\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})}|\sqrt{\kappa_{\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})}}|\nabla\varepsilon_{\widetilde{\bm{k}}_{j^{\prime}}(\theta_{\bm{r}})}|\sqrt{\kappa_{\widetilde{\bm{k}}_{j^{\prime}}(\theta_{\bm{r}})}}},$
(3)

with $|\nabla\varepsilon_{\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})}|$ the Fermi
velocity and $\kappa_{\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})}$ the curvature
of the Fermi contour evaluated at $\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})$.
The critical momenta $\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})$ are solutions
to the saddle-point equations in the small-gap limit, namely, points on the
Fermi contour for which the Fermi velocity is parallel to $\bm{r}$,

$\frac{\nabla\varepsilon_{\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})}}{|\nabla\varepsilon_{\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})}|}=\bm{\hat{r}},$
(4)

hence their dependence on the polar angle $\theta_{\bm{r}}$ which controls the
spatial anisotropy of the LDOS. We recall that the gapped nature of the
superconductor yields complex solutions to the saddle-point equations [40].
However, up to and including the linear order in
$\Delta/\varepsilon_{\mathrm{F}}$, these solutions [i.e. the vectors
$\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})$] are real and positioned on the
Fermi contour of the underlying metal. Therefore, the prefactor
$\Gamma_{j,j^{\prime}}(\theta_{\bm{r}})$ in Eq. (3) in the small-gap limit is
simply determined by the band structure of the normal metal (see Appendix B).
The exponential decay length $\xi_{j,j^{\prime}}(\theta_{\bm{r}})$ depends on
the Fermi velocity at $\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})$ and
$\widetilde{\bm{k}}_{j^{\prime}}(\theta_{\bm{r}})$, and
$f[\bm{k}^{\pm}_{j,j^{\prime}}(\theta_{\bm{r}})\cdot\bm{r}]$ is a linear
combination of oscillating functions with typical characteristic inverse
lengths
$\bm{k}^{\pm}_{j,j^{\prime}}(\theta_{\bm{r}})=\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})\pm\widetilde{\bm{k}}_{j^{\prime}}(\theta_{\bm{r}})$.
The indices $j,j^{\prime}=1,\dots,N$ label the various solutions to Eq. (4).
The form in Eq. (2), derived in Appendix B, is more compact than the form
first derived in Ref. [40], as it pairs up the two contributions from the
wavevectors whose respective Fermi velocities are parallel and antiparallel to
the observation direction [the latter satisfying Eq. (4) with
$\hat{\bm{r}}\rightarrow-\hat{\bm{r}}$]. The Fermi contour of NbSe2 consists
of three disconnected pockets centered at $\Gamma$, $K$ and $K^{\prime}$,
therefore, $N=3$.

To analyze the spatial anisotropy of the STM data we consider radial averaging
in annuli around the impurity and hence introduce the following quantity:

$I_{\mathrm{int}}(\theta_{\bm{r}})=\frac{1}{M}\sum_{\bm{R}_{m}\in
M(\theta_{\bm{r}})}|\bm{R}_{m}|\cdot I(\bm{R}_{m}),$ (5)

where $I(\bm{R}_{m})$ is the measured current at point $\bm{R}_{m}$ sitting at
a distance $|\bm{R}_{m}|$ from the impurity, while $M(\theta_{\bm{r}})$ is the
set of $M$ points belonging to an angular sector centered on the impurity,
having infinitesimal angular width $d\theta$, an inner radius $R_{0}$, and
outer radius $R_{0}+\delta$. As long as $R_{0}$ is sufficiently large to make
the asymptotic approximation reasonable ($R_{0}\gg
1/k_{\mathrm{F}}^{\mathrm{min}}\sim 0.2\;\mathrm{nm}$ for the Fermi contour
under consideration), we can relate $I_{\mathrm{int}}(\theta_{\bm{r}})$ to the
radially-integrated LDOS,

$I_{\mathrm{int}}(\theta_{\bm{r}})\sim\rho_{\mathrm{int}}^{R_{0},\delta}(\theta_{\bm{r}})=\frac{1}{\delta}\int_{R_{0}}^{R_{0}+\delta}\textrm{d}rr\rho_{\mathrm{max}}(\bm{r}),$
(6)

where $\rho_{\mathrm{max}}(\bm{r})$ is the term in Eq. (2) with dominant
prefactor,
$\Gamma_{\mathrm{max}}(\theta_{\bm{r}})=\mathrm{max}\\{\Gamma_{j,j^{\prime}}(\theta_{\bm{r}})\\}$.
This simplification is legitimate as the various prefactors
$\Gamma_{j,j^{\prime}}(\theta_{\bm{r}})$ are essentially in phase (see
Appendix C for further details). We can further simplify Eq. (6) by making the
following assumptions:

  * •

The integration range $\delta$ is larger than the typical oscillation length,
therefore $f[\bm{k}^{\pm}_{j,j^{\prime}}(\theta_{\bm{r}})\cdot\bm{r}]$ is
averaged out and can be safely ignored.

  * •

The exponential decay length is sufficiently large so that
$e^{-r/\xi(\theta)}\sim 1$ over the integration range. This demands
$R_{0}+\delta\leq\xi^{\mathrm{exp}}_{\mathrm{SC}}\sim
10\;\mathrm{nm}\ll\mathrm{min}\\{\xi_{j,j^{\prime}}(\theta_{\bm{r}})\\}$. This
approximation also relies on the fact that the angular dependence of the
prefactor and the exponential decay length are in phase, and therefore there
are not any competition effects.

Under these assumptions $\rho_{\mathrm{int}}^{R_{0},\delta}(\theta_{\bm{r}})$
becomes independent of $\delta$, and further, we have

$\rho_{\mathrm{int}}^{R_{0},\delta}(\theta_{\bm{r}})\sim\Gamma_{\mathrm{max}}(\theta_{\bm{r}}).$
(7)

In Fig. 2 we present the result of our analysis in comparison with the model
curve

$\Gamma_{\mathrm{model}}(\theta_{\bm{r}})=a\Gamma_{\mathrm{max}}(\theta_{\bm{r}})+a_{0}.$
(8)

Here $a$ is a free parameter in the theory, and $a_{0}$ is constant shift on
the order of the background signal. We estimate the latter by taking the
signal average far away from the impurity times the typical integration length
$\overline{\bm{R}}_{n}\sim
R_{0}+\frac{\delta_{\mathrm{min}}+\delta_{\mathrm{max}}}{2}$.

Figure 2: Angular dependence of the LDOS prefactor. (a) and (b) panels
correspond to data set 1 and 2, respectively. Purple-to-yellow curves
represent experimental data points $|I_{\mathrm{int}}|(\theta)$ for
$\delta$-values specified in the legend and $d\theta=2\pi/100$. Legend inset
shows the real-space current data with a lattice vector aligned with the
horizontal axis ($\theta_{\bm{r}}=0$). We excluded a range of data points in
data set 1 due to interference with another YSR state (partially visible on
the bottom right of the image).

Experimental data points $\rho_{\mathrm{int}}^{R_{0},\delta}(\theta)$ above
the background signal collapse onto a universal curve, thereby supporting the
assumptions stated in the previous paragraph.

The relative angular position of the model curve with respect to the data
points is not a degree of freedom, but it is determined by the orientation of
the lattice vectors (see legend insets in Fig. 2), the latter being extracted
from the constant-current image of the atomic surface. We emphasize the
agreement between the positions of the maxima and minima of the model and the
experimental curves. While the six-fold nature of the YSR state may follow
trivially from the symmetry of the band structure, its orientation with
respect to the lattice does not. One could generally expect YSR states to be
oriented according to high-symmetry lines; however, it is a priori equally
reasonable for the “petals” of the LDOS to be along the lattice vectors or in-
between. The theory shows how this is, in fact, determined by the geometry of
the Fermi surface and Fermi velocity.

We note that the periodicity of the data points with $\theta_{\bm{r}}$
exhibits certain higher harmonics, namely, the signal becomes thinner at the
crests and wider at the valleys, in other words, the “petals” are narrow. This
non-trivial feature is also present in the model curve.

###  II.2 Fourier-space analysis

To further confirm the validity of our interpretation, we study the Fourier
transform (FT) of the in-gap recorded current [Fig. 3 (a)], and we compare it
against the FT of the asymptotic expression of the LDOS, which contains the
prefactor, the exponential decay, and the oscillations,

$\widetilde{\rho}(\bm{q})=\sum_{j,j^{\prime}}\int_{0}^{2\pi}d\theta_{\bm{r}}\Gamma_{j,j^{\prime}}(\theta_{\bm{r}})\widetilde{F}_{j,j^{\prime}}(\bm{q},\theta_{\bm{r}}),$
(9)

with $\widetilde{F}_{j,j^{\prime}}(\bm{q},\theta_{\bm{r}})\equiv\int
dre^{-r/\xi_{j,j^{\prime}}(\theta_{\bm{r}})}f[(\bm{k}^{\pm}_{j,j^{\prime}}(\theta_{\bm{r}})-\bm{q})\cdot\bm{r}]$
the integral over the radial coordinate $r$ (see Appendix D for details).

Figure 3: (a) Fourier transform’s modulus of the in-gap current shown in Fig.
1 (b). The blue hexagon indicates the FBZ and the blue cross marks one of the
Bragg peaks. (b) Fourier transform’s modulus of the analytical expression of
the LDOS [Eq. (9)] evaluated on the FBZ. (c) same as (b) with constant
prefactor.

The experimental FT of the YSR states exhibits two salient features that are
qualitatively captured by the model [cf. Fig. 3 (a) and (b)]. The first
comprises the linear segments parallel to the edge of the hexagonal First
Brillouin Zone (FBZ), corresponding to the short wavelength oscillations of
the YSR state. The second is the six-fold star at small $\bm{q}$ whose
orientation is rotated by $\pi/6$ with respect to the YSR state in real-space
(cf. Fig. 1) as expected from the fact that a focused signal in real-space
gives a focused feature in Fourier space that is in orthogonal direction. The
small star in momentum space is the manifestation of the QFE at large length
scales. Finally, we corroborate that the long-range anisotropy of the YSR
state, i.e., the small-$\bm{q}$ star feature in Fourier space, is controlled
by the prefactor in the asymptotic expression of the LDOS: we set
$\Gamma_{j,j^{\prime}}(\theta_{\bm{r}})=1$ in Eq. (9) and obtain a FT-map
which differs from the full $\widetilde{\rho}(\bm{q})$ exactly in that it
lacks the small-$\bm{q}$ six-fold star [Fig. 3 (c)].

##  III Conclusions

We show that the theory of the quasiparticle focusing effect in two-
dimensional superconductors [40] can be further simplified in application to
STM data of superconductors such as 2H-NbSe2, in which the coherence length
derived from the band structure exceeds the spatial extent of YSR states. Our
conclusions are established both by a practical analysis of annuli in real
space, and by analysis of features in Fourier space.

The key physical quantity for all practical purposes becomes the prefactor
$\Gamma(\theta_{\bm{r}})$ of the LDOS expression, which faithfully captures
the symmetry, orientation, and angular profile of YSR states observed in STM
experiments. The prefactor has a simple analytical expression depending on the
Fermi velocity and curvature of the Fermi surface only. In particular, the
LDOS around a magnetic impurity is preferentially focused in directions
perpendicular to flatter sections of the Fermi contour, and in the small-gap
limit it matches the QFE in the underlying normal metal [39, 38]. It would be
interesting to test this matching in future measurements of local response to
impurities in the normal metallic state of NbSe2.

In conclusion, this work demonstrates the applicability of asymptotic
approximations in the spirit of [39, 40, 38] for describing the local response
to point defects, and underlines the importance of the band structure geometry
for designing collective impurity states.

###### Acknowledgements.

M.U., A.M. and P.S. acknowledge the support of the French Agence Nationale de
la Recherche (ANR), under grant number ANR-22-CE30-0037. M.U. acknowledges
funding by the European Union NextGenerationEU/PRTR-C17.I1 and the Department
of Education of the Basque Government (IKUR Strategy). F.M. acknowledges
Laurent Cario for providing the samples and funding from H2020 Marie
Skłodowska-Curie Actions (grant number 659247) and the ANR
(ANR-16-ACHN-0018-01).

##  Appendix A Experimental details and additional data

2H-NbSe2 single crystals were grown using an iodine transport method and were
unintentionally doped by a few tens of ppm of magnetic species (Fe, Cr, Mn)
contained in the Niobium precursor, see Ref. [29] for more details. The
crystals were mechanically cleaved in cryogenic vacuum at T $\sim$ 20 K and
directly inserted into the STM head [50] at 4.2 K. An etched atomically sharp
and stable tungsten tip was used for all measurements, which were performed at
$T=0.3$ K. In Fig. 4 we present the $dI/dV$ spectra for the two impurities
considered in this work and in Fig. 5 (b,d) we analyze the anisotropy of the
signal recorded at positive bias to complement the data presented in the main
text.

Figure 4: Tunneling spectrum taken at the core of the YSR state in Fig. 1a, b
(a) and Fig. 1c, d (b) showing the strongest in-gap peak at
$V_{\mathrm{S}}=-0.1$ meV and $V_{\mathrm{S}}=0.25$ meV, respectively. Figure
5: Same analysis of the LDOS prefactor as in Fig. 2 of the main text for the
in-gap current recorded at positive bias, $V=0.6$ mV. (a) and (b) correspond
to data sets 1 and 2, respectively. The phenomenological parameters in the
model curve [Eq. (8)] for panels [Fig. 2 (a, b), Fig. 5 (a, b)] are:
$a_{0}=(-3,-3,-2.5,-2.1)\cdot 10^{-12}$ pA $\cdot$ nm and
$a/a_{\mathrm{c}}=(4.45,4.45,3.2,4.45)\cdot 10^{-7}$ pA $\cdot$ nm $\cdot$
meV1/2, with the model curve computed using Eqs. (25) to (27b) and “Band 2”
parameters in Table 1.

##  Appendix B Derivation of the LDOS expression

We derive the expression for the in-gap LDOS in Eq. (2), which is equal in
value to the main result from [40], but is further transformed to make
explicit that LDOS is real, and that the doubling of critical momenta in the
superconductor with respect to the normal metal does not essentially alter the
prefactor. We further show that the prefactor up to linear order in
$\Delta/\varepsilon_{\textrm{F}}$ matches the one of the normal metal.

The bare propagator at energy $E<\Delta$ from $\bm{r}_{a}=\bm{r}$ to
$\bm{r}_{b}=\bm{0}$ on an s-wave superconductor (described by the Hamiltonian
in Eq. (1b) of the main text) in Nambu space can be approximated in the far-
field limit ($rk_{F,\mathrm{min}}\gg 1$, with $k_{F,\mathrm{min}}$ the minimum
Fermi wave vector of the Fermi contour) as

$\hat{G}_{0}(\bm{r},\bm{0};E)\sim\frac{1}{\sqrt{r}}\sum_{j,\epsilon=\pm}\Gamma_{j,\epsilon}(\theta_{\bm{r}})\exp\left\\{-\frac{r}{\xi_{j,\epsilon}(\theta_{\bm{r}})}+i\left[\bm{r}\cdot\operatorname{Re}\left[\bm{k}_{j,\epsilon}(\theta_{\bm{r}})\right]+\epsilon\frac{\pi}{4}\right]\right\\}\begin{pmatrix}E+i\epsilon\omega&\Delta\\\
\Delta&E-i\epsilon\omega\end{pmatrix},$ (10)

where $\theta_{\bm{r}}$ is the polar angle on the plane defined by the vector
$\bm{r}$ from the impurity to the point where we seek to evaluate the LDOS,
$\omega=\sqrt{\Delta^{2}-E^{2}}$ and

$\displaystyle\Gamma_{j,\epsilon}(\theta_{\bm{r}})$
$\displaystyle=\frac{1}{|\nabla\varepsilon_{\bm{k}_{j,\epsilon}(\theta_{\bm{r}})}|\sqrt{\kappa_{\bm{k}_{j,\epsilon}(\theta_{\bm{r}})}}},$
(11a) $\displaystyle\xi_{j,\epsilon}(\theta_{\bm{r}})$
$\displaystyle=\frac{1}{\operatorname{Im}\left[\bm{k}_{j,\epsilon}(\theta_{\bm{r}})\right]\cdot\hat{\bm{r}}},$
(11b)

with $|\nabla\varepsilon_{\bm{k}_{j,\epsilon}(\theta_{\bm{r}})}|$ norm of the
gradient of the energy dispersion $\varepsilon_{\bm{k}}$, and
$\kappa_{\bm{k}_{j,\epsilon}(\theta_{\bm{r}})}$ the curvature of the Fermi
contour $\varepsilon_{\bm{k}}=0$ evaluated at
$\bm{k}_{j,\epsilon}(\theta_{\bm{r}})$. The points
$\bm{k}_{j,\epsilon}(\theta_{\bm{r}})$ are the solutions to the saddle-point
equations, namely complex momenta satisfying

$\displaystyle\varepsilon_{\bm{k}_{j,\epsilon}(\theta_{\bm{r}})}$
$\displaystyle=i\epsilon\omega,$ (12a)
$\displaystyle\frac{\nabla\varepsilon_{\bm{k}_{j,\epsilon}(\theta_{\bm{r}})}}{|\nabla\varepsilon_{\bm{k}_{j,\epsilon}(\theta_{\bm{r}})}|}$
$\displaystyle=\epsilon\bm{\hat{r}},$ (12b)

where $j=1,\dots,N$ denote the solutions from multiple disconnected pockets of
the Fermi contour, $\epsilon=\pm$ and $\bm{\hat{r}}$ the unit vector defined
by $\theta_{\bm{r}}$. For our present purpose of analyzing a small-gap
superconductor, we consider a perturbative expansion of the saddle-point
equations in $\Delta/\varepsilon_{\mathrm{F}}$, noting that the normal metal
is recovered as the pairing vanishes. This forces two simultaneous limits: (1)
$\Delta/\varepsilon_{\mathrm{F}}\rightarrow 0$, and (2)
$E_{S}\rightarrow\Delta$, which due to the only in-gap state being at
$E=E_{S}$ effectively means that in Eq. (12a) one has $\omega\rightarrow 0$.
In Ref. [40] it was shown that to linear order in
$\Delta/\varepsilon_{\mathrm{F}}$, the angle-dependent quantities in the
exponent of the propagator (10) take the following form:

$\displaystyle\operatorname{Re}\left[\bm{k}_{j,\epsilon}(\theta_{\bm{r}})\right]$
$\displaystyle=\widetilde{\bm{k}}_{j,\epsilon}(\theta_{\bm{r}}),$ (13a)
$\displaystyle\xi_{j,\epsilon}(\theta_{\bm{r}})$
$\displaystyle=\frac{|\nabla\varepsilon_{\widetilde{\bm{k}}_{j,\epsilon}(\theta_{\bm{r}})}|}{\omega},$
(13b)

with $\widetilde{\bm{k}}_{j,\epsilon}(\theta_{\bm{r}})$ the solutions to the
saddle-point equations in the $\omega\rightarrow 0$ limit,

$\displaystyle\varepsilon_{\widetilde{\bm{k}}_{j,\epsilon}(\theta_{\bm{r}})}$
$\displaystyle=0,$ (14a)
$\displaystyle\frac{\nabla\varepsilon_{\widetilde{\bm{k}}_{j,\epsilon}(\theta_{\bm{r}})}}{|\nabla\varepsilon_{\widetilde{\bm{k}}_{j,\epsilon}(\theta_{\bm{r}})}|}$
$\displaystyle=\epsilon\bm{\hat{r}},$ (14b)

that is, real momenta positioned on the Fermi contour. To zeroth order in
$\Delta/\varepsilon_{\mathrm{F}}$, it is obvious that the prefactor is real
and takes the same form as in the normal-metal propagator, i.e., Eq. (11a)
evaluated at $\widetilde{\bm{k}}_{j,\epsilon}(\theta_{\bm{r}})$,

$\Gamma_{j,\epsilon}(\theta_{\bm{r}})=\frac{1}{|\nabla\varepsilon_{\widetilde{\bm{k}}_{j,\epsilon}(\theta_{\bm{r}})}|\sqrt{\kappa_{\widetilde{\bm{k}}_{j,\epsilon}(\theta_{\bm{r}})}}}.$
(15)

The question arises whether there is any non-trivial correction to next-
leading order, i.e., linear in $\Delta/\varepsilon_{\mathrm{F}}$. To show that
there is no correction, we consider an expansion of the prefactor in terms of
a small complex quantity $z$ with
$\bm{k}_{j,\epsilon}(\theta_{\bm{r}})=\widetilde{\bm{k}}_{j,\epsilon}(\theta_{\bm{r}})+z$,
knowing that $z\rightarrow 0$ as $\Delta/\varepsilon_{\mathrm{F}}\rightarrow
0$. The key point now is that the angular profile of the measured YSR state
depends on the amplitude of the prefactor, so that we write
$\Gamma_{j,\epsilon}(\theta_{\bm{r}})\equiv|\Gamma|\exp(i\Phi)$ and consider
the expansion of the real function evaluated for a complex value of its
variable,
$|\Gamma(z)|\approx\sqrt{\left[\Gamma(k)+\frac{1}{2}\Gamma^{\prime\prime}(k)\textrm{Re}(z^{2})\right]^{2}+\left[\Gamma^{\prime}(k)\textrm{Im}(z)+\frac{1}{2}\Gamma^{\prime\prime}(k)\textrm{Im}(z^{2})\right]^{2}}\approx\Gamma(k)+O(\Delta^{2})$,
where $k\equiv\widetilde{\bm{k}}_{j,\epsilon}(\theta_{\bm{r}})$ is the real
solution in the metal [Eq. (14)]. We used that $\Gamma(k)$ is real, and
further assumed that the correction $\textrm{Im}(z)$ is at least linear in
$\Delta$, while the $\textrm{Re}(z)$ is at least quadratic, which is
reasonable due to the structure of Eq. (12a). We note that these properties of
$z$ can be confirmed in the particular example of an elliptic Fermi contour
worked out in Ref. [40]. For the second step in the expansion we also assumed
$\Gamma(k)\neq 0$, which is a condition for the entire saddle-point approach
to be valid anyway. Hence the correction of the amplitude of the prefactor is
at most second order in $\Delta/\varepsilon_{\mathrm{F}}$, so that up to and
including first order the prefactor is calculated as for the normal metal [Eq.
(15)]. We note that the prefactor’s phase $\Phi$ can in general be linear in
$\Delta/\varepsilon_{\mathrm{F}}$, thereby contributing an additional
$\theta_{\bm{r}}$-dependent phase to the bare propagator (and hence the LDOS)
whose analysis is beyond the scope of this work.

Continuing now the derivation of the LDOS, to all orders in perturbation
theory in the strength of the impurity potential $\hat{V}$, the modification
of the LDOS due to a point-like scatterer at $\bm{r}=\bm{0}$ is given by

$\rho(\bm{r},E)=-\frac{1}{\pi}\operatorname{Im}\operatorname{Tr}\left[\hat{M}\hat{G}_{0}(\bm{r},\bm{0};E)\hat{T}(E)\hat{G}_{0}(\bm{0},\bm{r};E)\right],$
(16)

where
$\hat{T}(E)=\hat{V}\left[1-\hat{G_{0}}(\bm{0},\bm{0};E)\hat{V}\right]^{-1}$ is
the T-matrix and $\hat{M}=(\tau_{0}+\tau_{z})/2$ projects out the electron-
electron component. Here and in the following $\tau_{0,x,y,z}$ denote the
Pauli matrices in Nambu space. Note that the counter-propagator satisfies
$\hat{G}_{0}(\bm{0},\bm{r};E)=\hat{G}_{0}^{*}(\bm{r},\bm{0};E)$. Since the
anisotropy of the LDOS is encoded in the bare propagator (10), we simplify by
considering a magnetic scatterer only, without non-magnetic scattering
amplitude, i.e. $\hat{V}=-J\tau_{0}$. Further, we neglect the anisotropy of
the band structure when evaluating $\hat{G_{0}}(\bm{0},\bm{0};E)$ in the
T-matrix,

$\hat{G_{0}}(\bm{0},\bm{0};E)\equiv\frac{1}{\mathcal{V}}\sum_{\bm{k}}\hat{G}_{0}(\bm{k},E)\sim-\frac{\pi\nu_{0}}{\sqrt{\Delta^{2}-E^{2}}}\begin{pmatrix}E&\Delta\\\
\Delta&E\end{pmatrix},$ (17)

where $\nu_{0}$ is the density of states at the Fermi level. Under these
assumptions and in the small-gap limit we have

$\rho(\bm{r})\sim-\frac{1}{\pi}\operatorname{Im}\frac{1}{r}\sum_{j,j^{\prime},\epsilon,\epsilon^{\prime}}\Gamma_{j,j^{\prime}}^{\epsilon,\epsilon^{\prime}}(\theta_{\bm{r}})\exp\left\\{-\frac{r}{\xi_{j,j^{\prime}}^{\epsilon,\epsilon^{\prime}}(\theta_{\bm{r}})}+i\left[\bm{r}\cdot\left(\widetilde{\bm{k}}_{j,\epsilon}(\theta_{\bm{r}})-\widetilde{\bm{k}}_{j^{\prime},\epsilon^{\prime}}(\theta_{\bm{r}})\right)+(\epsilon-\epsilon^{\prime})\frac{\pi}{4}\right]\right\\}\frac{N_{\epsilon,\epsilon^{\prime}}}{D},$
(18)

where

$\displaystyle\Gamma_{j,j^{\prime}}^{\epsilon,\epsilon^{\prime}}(\theta_{\bm{r}})$
$\displaystyle=\Gamma_{j,\epsilon}(\theta_{\bm{r}})\Gamma_{j^{\prime},\epsilon^{\prime}}(\theta_{\bm{r}}),$
(19a)
$\displaystyle\xi_{j,j^{\prime}}^{\epsilon,\epsilon^{\prime}}(\theta_{\bm{r}})$
$\displaystyle=\left(\frac{1}{\xi_{j,\varepsilon}(\theta_{\bm{r}})}+\frac{1}{\xi_{j^{\prime},\varepsilon^{\prime}}(\theta_{\bm{r}})}\right)^{-1},$
(19b) $\displaystyle N_{\epsilon,\epsilon^{\prime}}$
$\displaystyle=t_{0}(\Delta^{2}+E^{2})+2E\Delta
t_{x}+t_{0}\omega^{2}\epsilon\epsilon^{\prime}+$ (19c)
$\displaystyle+i\omega(Et_{0}+\Delta t_{x})(\epsilon-\epsilon^{\prime}),$
$\displaystyle D$ $\displaystyle=t_{0}^{2}-t_{x}^{2},$ (19d)

with $t_{0}=1-\alpha\frac{E}{\omega}$, $t_{x}=\alpha\frac{\Delta}{\omega}$ and
$\alpha=\nu_{0}\pi J$, and $\Gamma_{j,\epsilon}(\theta_{\bm{r}})$ and
$\xi_{j,\epsilon}(\theta_{\bm{r}})$ as defined in Eqs. (15) and (13b),
respectively. The energy of the in-gap bound states are given by the poles of
the T-matrix, i.e. $D=0$, and takes the well known form,

$E_{\mathrm{S}}^{\pm}=\pm\Delta\frac{1-\alpha^{2}}{1+\alpha^{2}}.$ (20)

To obtain the LDOS of the YSR state [51], we evaluate
$N_{\epsilon,\epsilon^{\prime}}$ at $E=E_{\mathrm{S}}^{+}$ and expand $D$ to
first order around $E=E_{\mathrm{S}}^{+}$ to obtain up to a constant
prefactor,

$\rho(\bm{r})\sim-\frac{1}{\pi}\operatorname{Im}\frac{1}{r}\sum_{j,j^{\prime},\epsilon,\epsilon^{\prime}}\Gamma_{j,j^{\prime}}^{\epsilon,\epsilon^{\prime}}(\theta_{\bm{r}})\exp\left\\{-\frac{r}{\xi_{j,j^{\prime}}^{\epsilon,\epsilon^{\prime}}(\theta_{\bm{r}})}+i\bm{r}\cdot\left(\widetilde{\bm{k}}_{j,\epsilon}(\theta_{\bm{r}})-\widetilde{\bm{k}}_{j^{\prime},\epsilon^{\prime}}(\theta_{\bm{r}})\right)\right\\}\frac{X_{\epsilon,\epsilon^{\prime}}}{E-E_{\mathrm{S}}},$
(21)

where
$X_{\epsilon,\epsilon^{\prime}}=(1+\alpha^{2})e^{i\varphi_{\epsilon,\epsilon^{\prime}}}$
with

$\varphi_{\epsilon,\epsilon^{\prime}}=\begin{cases}0&\mathrm{if}\;\epsilon=\epsilon^{\prime},\\\
\epsilon\left[\frac{\pi}{2}+\operatorname{arg}\left(1-\alpha^{2}+i2\alpha\right)\right]&\mathrm{if}\;\epsilon\neq\epsilon^{\prime}.\end{cases}$
(22)

We note that including a non-magnetic component to the scattering potential
yields a different expression for $X_{\epsilon,\epsilon^{\prime}}$, but
crucially, with the same dependence on $\epsilon,\epsilon^{\prime}$.

The expression in Eq. (21) can be written in a more compact form by noticing
that, owing to the evenness of the energy dispersion $\varepsilon_{\bm{k}}$,
for each point $\widetilde{\bm{k}}_{j,+}(\theta_{\bm{r}})$, there exists a
point $\widetilde{\bm{k}}_{l,-}(\theta_{\bm{r}})$ such that
$\widetilde{\bm{k}}_{l,-}(\theta_{\bm{r}})=-\widetilde{\bm{k}}_{j,+}(\theta_{\bm{r}})$.
If the $j$ pocket is centered around the $\Gamma$ point, then $j=l$, but this
is not the case in general. The prefactor and characteristic decay length
satisfy that
$\Gamma_{j,+}(\theta_{\bm{r}})=\Gamma_{l,-}(\theta_{\bm{r}})\equiv\Gamma_{j}(\theta_{\bm{r}})$
and
$\xi_{j,+}(\theta_{\bm{r}})=\xi_{l,-}(\theta_{\bm{r}})\equiv\xi_{j}(\theta_{\bm{r}})$,
therefore, by relabelling the critical points as
$\\{\epsilon\widetilde{\bm{k}_{j}}(\theta_{\bm{r}})\\}_{\begin{subarray}{l}j=1\dots
N\\\ \epsilon=\pm\end{subarray}}$, it becomes obvious that the only
contribution to the LDOS comes from the poles of the T-matrix,

$\rho(\bm{r})\sim\frac{1}{r}\sum_{j,j^{\prime}}\Gamma_{j,j^{\prime}}(\theta_{\bm{r}})e^{-r/\xi_{j,j^{\prime}}(\theta_{\bm{r}})}\sum_{\epsilon,\epsilon^{\prime}}\exp\left\\{i\left[\bm{r}\cdot\left(\epsilon\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})-\epsilon^{\prime}\widetilde{\bm{k}}_{j^{\prime}}(\theta_{\bm{r}})\right)+\varphi_{\epsilon,\epsilon^{\prime}}\right]\right\\}\delta(E-E_{\mathrm{S}}),$
(23)

where Eq. (2) in the main text is recovered by introducing

$f[\bm{k}^{\pm}_{j,j^{\prime}}(\theta_{\bm{r}})\cdot\bm{r}]\equiv\sum_{\epsilon,\epsilon^{\prime}}\exp\left\\{i\left[\bm{r}\cdot\left(\epsilon\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})-\epsilon^{\prime}\widetilde{\bm{k}}_{j^{\prime}}(\theta_{\bm{r}})\right)+\varphi_{\epsilon,\epsilon^{\prime}}\right]\right\\},$
(24)

that is a real-valued function. With the new notation the $\epsilon=-$
solutions to Eq. (28a) are constructed from the $\epsilon=+$ solutions,
therefore Eq. (4) in the main text considers the parallel solutions only.
Note, however, that the antiparallel solutions are encoded in
$f[\bm{k}^{\pm}_{j,j^{\prime}}(\theta_{\bm{r}})\cdot\bm{r}]$.

##  Appendix C Analysis of the tight-binding energy dispersion

The starting point to compute the model curve for the prefactor (Fig. 2 in the
main text) is a tight-binding energy dispersion that describes the two lowest-
lying Nb 4d bands of 2H-NbSe2,

$\varepsilon_{\bm{k}}=\mu+\sum_{i=1}^{5}t_{i}g_{i}(\bm{k}),$ (25)

where $t_{i}$ is the hopping amplitude to the $i$-th nearest neighbor and

$\displaystyle g_{1}(\bm{k})$ $\displaystyle=2\cos\zeta\cos\eta+\cos 2\zeta,$
(26a) $\displaystyle g_{2}(\bm{k})$ $\displaystyle=2\cos 3\zeta\cos\eta+\cos
2\eta,$ (26b) $\displaystyle g_{3}(\bm{k})$ $\displaystyle=2\cos 2\zeta\cos
2\eta+\cos 4\zeta,$ (26c) $\displaystyle g_{4}(\bm{k})$
$\displaystyle=\cos\zeta\cos 3\eta+\cos 5\zeta\cos\eta,+\cos 4\zeta\cos
2\eta,$ (26d) $\displaystyle g_{5}(\bm{k})$ $\displaystyle=2\cos 3\eta\cos
3\zeta+\cos 6\zeta,$ (26e)

with $\zeta=\frac{1}{2}k_{x}a_{\mathrm{c}}$ and
$\zeta=\frac{\sqrt{3}}{2}k_{y}a_{\mathrm{c}}$, and $a_{\mathrm{c}}$ is the
lattice constant. The hopping amplitudes are presented in Table 1. We follow
Refs. [29, 44] and assume that the adatom primarily couples to one of the
bands only. Specifically, we take Band 2, which we treat as a spin-degenerate
energy dispersion.

| $\mu$ | $t_{1}$ | $t_{2}$ | $t_{3}$ | $t_{4}$ | $t_{5}$
---|---|---|---|---|---|---
Band 1 | 10.9 | 86.8 | 139.9 | 29.6 | 3.5 | 3.3
Band 2 | 203.0 | 46.0 | 257.5 | 4.4 | -15.0 | 6.0
Table 1: Fitted tight-binding parameters for the two lowest-lying Nb 4d bands
of 2H-NbSe2. All parameters in meV, for $a_{\mathrm{c}}=3.444$ Å. Extracted
from Ref. [41]

The Fermi velocity and the curvature of the Fermi contour that control the
anisotropy of the LDOS are related to the energy dispersion as follows,

$\displaystyle|\nabla\varepsilon_{\bm{k}}|$
$\displaystyle=\sqrt{(\partial_{k_{x}}\varepsilon_{\bm{k}})^{2}+(\partial_{k_{y}}\varepsilon_{\bm{k}})^{2}},$
(27a) $\displaystyle\kappa_{\bm{k}}$
$\displaystyle=\frac{(\partial_{k_{x}}\varepsilon_{\bm{k}})^{2}\partial^{2}_{k_{y}}\varepsilon_{\bm{k}}+(\partial_{k_{y}}\varepsilon_{\bm{k}})^{2}\partial^{2}_{k_{x}}\varepsilon_{\bm{k}}-2\partial_{k_{y}}\varepsilon_{\bm{k}}\partial_{k_{y}}\varepsilon_{\bm{k}}\partial^{2}_{k_{x},k_{y}}\varepsilon_{\bm{k}}}{|\nabla\varepsilon_{\bm{k}}|^{3}}.$
(27b)

To compute the model curve for the prefactor we first obtain the critical
points in the FBZ by numerically solving

$\displaystyle\varepsilon_{\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})}$
$\displaystyle=0,$ (28a)
$\displaystyle\frac{\nabla\varepsilon_{\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})}}{|\nabla\varepsilon_{\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})}|}$
$\displaystyle=\bm{\hat{r}},$ (28b)

for a discrete set of $\theta_{\bm{r}}$ values. As the Fermi contour of the
energy dispersion in Eq. (25) has three non-equivalent pockets, we obtain six
critical points
$\\{\epsilon\widetilde{\bm{k}_{j}}(\theta_{\bm{r}})\\}_{\begin{subarray}{l}j=1,2,3\\\
\epsilon=\pm\end{subarray}}$ for a given direction in real space determined by
$\theta_{\bm{r}}$ (see Fig. 6). This yields a set of six distinct
$\Gamma_{j,j^{\prime}}(\theta_{\bm{r}})$ curves that are essentially in phase
[see Fig. 7 (a)]. Crucially, the prefactor curves are in phase with respect to
the characteristic decay length curves as well [Fig. 7 (b)], which justifies
the simplification of taking the dominant prefactor
$\Gamma_{\mathrm{max}}(\theta_{\bm{r}})=\mathrm{max}\\{\Gamma_{j,j^{\prime}}(\theta_{\bm{r}})\\}$
in the model curve.

Figure 6: Examples of critical points for various values of $\theta_{\bm{r}}$.
The blue lines indicate the Fermi contour for the chosen 2H-NbSe2 model, i.e.
Eq. (25) with parameters of Band 2. The black arrow represents $\hat{\bm{r}}$,
i.e., a vector in real space from the impurity location to an arbitrary
position where we evaluate the LDOS. The orange dots on the Fermi contour
indicate the corresponding critical points
$\epsilon\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})$, i.e., the solutions to Eqs.
(28a, 28b) ($\epsilon=+$) and their spatially-inverted counterparts
($\epsilon=-$). The orange arrows indicate the normalized gradient of the
energy dispersion at that point, parallel ($\epsilon=+$) and anti-parallel
($\epsilon=-$) to $\hat{\bm{r}}$. (a), (b) and (c) correspond to
$\theta_{\bm{r}}=\pi/6,\pi/3,3\pi/2$, respectively. Figure 7: (a) Angular
dependence of the LDOS prefactor. (b) Angular dependence of LDOS exponential
decay length with $\omega=1$ meV. Each curve denotes an interference term
$(j,j^{\prime})$ with $j,j^{\prime}=1,2,3$ for pockets centered at
$\Gamma=(0,0)$, $K=(0,4\pi/3a)$ and $K^{\prime}=-K$. Certain curves are three-
fold periodic only reflecting the reduced symmetry of the $K,K^{\prime}$
pockets [e.g. $(2,2)$], yet, the six-fold symmetry of the band structure is
recovered when taking into account all the interference terms. The curve
presented in Fig. 2 in the main text corresponds to curve $(1,1)$ in panel
(a). Compare with Fig. 6 to observe that both the prefactor and the
exponential decay length are maximal when $\hat{\bm{r}}$ is perpendicular to
the flattest sections of the Fermi contour.

##  Appendix D Derivation of the Fourier transform of the LDOS

We derive here the Fourier transform of the LDOS at the YSR-state energy:

$\begin{split}\widetilde{\rho}(\bm{q})&\equiv\int
d\bm{r}e^{-i\bm{q}\bm{r}}\rho(\bm{r})\\\
&=\sum_{j,j^{\prime}}\int_{0}^{\infty}dr\int_{0}^{2\pi}d\theta_{\bm{r}}\Gamma_{j,j^{\prime}}(\theta_{\bm{r}})\sum_{\epsilon,\epsilon^{\prime}}\exp\left\\{r\left[i\bm{p}_{\epsilon,\epsilon^{\prime}}(\bm{q},\theta_{\bm{r}})\cdot\hat{\bm{r}}-\frac{1}{\xi_{j,j^{\prime}}(\theta_{\bm{r}})}\right]\right\\}\\\
&=\sum_{j,j^{\prime}}\int_{0}^{2\pi}d\theta_{\bm{r}}\Gamma_{j,j^{\prime}}(\theta_{\bm{r}})\sum_{\epsilon,\epsilon^{\prime}}e^{i\varphi_{\epsilon,\epsilon^{\prime}}}\frac{1}{i\bm{p}_{\epsilon,\epsilon^{\prime}}(\bm{q},\theta_{\bm{r}})\cdot\hat{\bm{r}}-1/\xi_{j,j^{\prime}}(\theta_{\bm{r}})},\end{split}$
(29)

where
$\bm{p}_{\epsilon,\epsilon^{\prime}}(\bm{q},\theta_{\bm{r}})=\epsilon\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})-\epsilon^{\prime}\widetilde{\bm{k}}_{j^{\prime}}(\theta_{\bm{r}})-\bm{q}$,
and Eq. (9) in the main text is recovered by setting

$\widetilde{F}_{j,j^{\prime}}(\bm{q},\theta_{\bm{r}})\equiv\sum_{\epsilon,\epsilon^{\prime}}e^{i\varphi_{\epsilon,\epsilon^{\prime}}}\frac{1}{i\bm{p}_{\epsilon,\epsilon^{\prime}}(\bm{q},\theta_{\bm{r}})\cdot\hat{\bm{r}}-1/\xi_{j,j^{\prime}}(\theta_{\bm{r}})}.$
(30)

The maps in Fig. 3 (b) and (c) are obtained by replacing the angular integral
by a sum over a discrete set of $\\{\theta_{\bm{r}}\\}$ values, for each of
which the term
$\sum_{j,j^{\prime}}\Gamma_{j,j^{\prime}}(\theta_{\bm{r}})\widetilde{F}_{j,j^{\prime}}(\bm{q},\theta_{\bm{r}})$
is evaluated. Importantly, to ensure that $\widetilde{\rho}(\bm{q})$ has the
periodicity of the reciprocal lattice, the pairs of critical points
$\epsilon\widetilde{\bm{k}}_{j}(\theta_{\bm{r}})-\epsilon^{\prime}\widetilde{\bm{k}}_{j^{\prime}}(\theta_{\bm{r}})$
have to be mapped onto the FBZ. Further, as it was shown in Appendix B, the
$\varphi_{\epsilon,\epsilon^{\prime}}$ phase depends on the nature of the
scatterer. For simplicity, we set $\alpha=1$ in Eq. (22) which assumes the
strong-coupling regime, i.e. $E_{\mathrm{S}}=0$.

## References

  * Rusinov [1969] A. I. Rusinov, Sov. J. Exp. Theor. Phys. Lett 9, 85 (1969).
  * Shiba [1968] H. Shiba, Progr. Theor. Phys. 40, 435 (1968).
  * Yu [1965] L. Yu, Acta Phys. Sin. 21, 75 (1965).
  * Heinrich  _et al._ [2018] B. W. Heinrich, J. I. Pascual, and K. J. Franke, Progress in Surface Science 93, 1 (2018).
  * Kaladzhyan  _et al._ [2016] V. Kaladzhyan, C. Bena, and P. Simon, Phys. Rev. B 93, 214514 (2016).
  * Perrin  _et al._ [2020] V. Perrin, F. L. N. Santos, G. C. Ménard, C. Brun, T. Cren, M. Civelli, and P. Simon, Phys. Rev. Lett. 125, 117003 (2020).
  * Sukhachov  _et al._ [2023] P. O. Sukhachov, F. von Oppen, and L. I. Glazman, Phys. Rev. B 108, 024505 (2023).
  * Choi  _et al._ [2017] D.-J. Choi, C. Rubio-Verdú, J. De Bruijckere, M. M. Ugeda, N. Lorente, and J. I. Pascual, Nat. Commun. 8, 15175 (2017).
  * Choi  _et al._ [2018] D.-J. Choi, C. G. Fernández, E. Herrera, C. Rubio-Verdú, M. M. Ugeda, I. Guillamón, H. Suderow, J. I. Pascual, and N. Lorente, Phys. Rev. Lett. 120, 167001 (2018).
  * Huang  _et al._ [2020] H. Huang, C. Padurariu, J. Senkpiel, R. Drost, A. L. Yeyati, J. C. Cuevas, B. Kubala, J. Ankerhold, K. Kern, and C. R. Ast, Nat. Phys. 16, 1227 (2020).
  * Huang  _et al._ [2021] H. Huang, J. Senkpiel, C. Padurariu, R. Drost, A. Villas, R. L. Klees, A. L. Yeyati, J. C. Cuevas, B. Kubala, J. Ankerhold, K. Kern, and C. R. Ast, Phys. Rev. Research 3, L032008 (2021).
  * Schneider  _et al._ [2021] L. Schneider, P. Beck, J. Wiebe, and R. Wiesendanger, Science Advances 7, eabd7302 (2021).
  * Vaxevani  _et al._ [2022] K. Vaxevani, J. Li, S. Trivini, J. Ortuzar, D. Longo, D. Wang, and J. I. Pascual, Nano Lett. 22, 6075–6082 (2022).
  * Nadj-Perge  _et al._ [2013] S. Nadj-Perge, I. K. Drozdov, B. A. Bernevig, and A. Yazdani, Phys. Rev. B 88, 020407 (2013).
  * Pientka  _et al._ [2013] F. Pientka, L. I. Glazman, and F. von Oppen, Phys. Rev. B 88, 155420 (2013).
  * Braunecker and Simon [2013] B. Braunecker and P. Simon, Phys. Rev. Lett. 111, 147202 (2013).
  * Nakosai  _et al._ [2013] S. Nakosai, Y. Tanaka, and N. Nagaosa, Phys. Rev. B 88, 180503 (2013).
  * Nadj-Perge  _et al._ [2014] S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon, J. Seo, A. H. MacDonald, B. A. Bernevig, and A. Yazdani, Science 346, 602–607 (2014).
  * Röntynen and Ojanen [2015] J. Röntynen and T. Ojanen, Phys. Rev. Lett. 114, 236803 (2015).
  * Li  _et al._ [2016] J. Li, T. Neupert, Z. Wang, A. MacDonald, A. Yazdani, and B. A. Bernevig, Nat. Commun. 7, 12297 (2016).
  * Jeon  _et al._ [2017] S. Jeon, Y. Xie, J. Li, Z. Wang, B. A. Bernevig, and A. Yazdani, Science 358, 772–776 (2017).
  * Ruby  _et al._ [2017] M. Ruby, B. W. Heinrich, Y. Peng, F. von Oppen, and K. J. Franke, Nano Lett. 17, 4473 (2017).
  * Kim  _et al._ [2018] H. Kim, A. Palacio-Morales, T. Posske, L. Rózsa, K. Palotás, L. Szunyogh, M. Thorwart, and R. Wiesendanger, Science Advances 4, eaar5251 (2018).
  * Palacio-Morales  _et al._ [2019] A. Palacio-Morales, E. Mascot, S. Cocklin, H. Kim, S. Rachel, D. K. Morr, and R. Wiesendanger, Science Advances 5, eaav6600 (2019).
  * Liebhaber  _et al._ [2022] E. Liebhaber, L. M. Rütten, G. Reecht, J. F. Steiner, S. Rohlf, K. Rossnagel, F. von Oppen, and K. J. Franke, Nat. Commun. 13, 10.1038/s41467-022-29879-0 (2022).
  * Farinacci  _et al._ [2023] L. Farinacci, G. Reecht, F. von Oppen, and K. J. Franke, arXiv:2307.09993 10.48550/arXiv.2307.09993 (2023).
  * Soldini  _et al._ [2023] M. O. Soldini, F. Küster, G. Wagner, S. Das, A. Aldarawsheh, R. Thomale, S. Lounis, S. S. P. Parkin, P. Sessi, and T. Neupert, Nat. Phys. 19, 1848–1854 (2023).
  * Yazdani  _et al._ [2023] A. Yazdani, F. Von Oppen, B. I. Halperin, and A. Yacoby, Science 380, eade0850 (2023).
  * Ménard  _et al._ [2015] G. C. Ménard, S. Guissart, C. Brun, S. Pons, V. S. Stolyarov, F. Debontridder, M. V. Leclerc, E. Janod, L. Cario, D. Roditchev, P. Simon, and T. Cren, Nat. Phys. 11, 1013 (2015).
  * Kim  _et al._ [2020] H. Kim, L. Rózsa, D. Schreyer, E. Simon, and R. Wiesendanger, Nat. Commun. 11, 4573 (2020).
  * Thupakula  _et al._ [2022] U. Thupakula, V. Perrin, A. Palacio-Morales, L. Cario, M. Aprili, P. Simon, and F. Massee, Phys. Rev. Lett. 128, 247001 (2022).
  * Kezilebieke  _et al._ [2018] S. Kezilebieke, M. Dvorak, T. Ojanen, and P. Liljeroth, Nano Lett. 18, 2311–2315 (2018).
  * Liebhaber  _et al._ [2019] E. Liebhaber, S. Acero González, R. Baba, G. Reecht, B. W. Heinrich, S. Rohlf, K. Rossnagel, F. von Oppen, and K. J. Franke, Nano Lett. 20, 339–344 (2019).
  * Ptok  _et al._ [2017] A. Ptok, S. Głodzik, and T. Domański, Phys. Rev. B 96, 184425 (2017).
  * Sticlet and Morari [2019] D. Sticlet and C. Morari, Phys. Rev. B 100, 075420 (2019).
  * Rütten  _et al._ [2024] L. M. Rütten, H. Schmid, E. Liebhaber, G. Franceschi, A. Yazdani, G. Reecht, K. Rossnagel, F. von Oppen, and K. J. Franke, arXiv:2404.16513 10.48550/arXiv.2404.16513 (2024).
  * Ortuzar  _et al._ [2022] J. Ortuzar, S. Trivini, M. Alvarado, M. Rouco, J. Zaldivar, A. L. Yeyati, J. I. Pascual, and F. S. Bergeret, Phys. Rev. B 105, 245403 (2022).
  * Weismann  _et al._ [2009] A. Weismann, M. Wenderoth, S. Lounis, P. Zahn, N. Quaas, R. G. Ulbrich, P. H. Dederichs, and S. Blügel, Science 323, 1190 (2009).
  * Lounis  _et al._ [2011] S. Lounis, P. Zahn, A. Weismann, M. Wenderoth, R. G. Ulbrich, I. Mertig, P. H. Dederichs, and S. Blügel, Phys. Rev. B 83, 035427 (2011).
  * Uldemolins  _et al._ [2022] M. Uldemolins, A. Mesaros, and P. Simon, Phys. Rev. B 105, 144503 (2022).
  * Rahn  _et al._ [2012] D. J. Rahn, S. Hellmann, M. Kalläne, C. Sohrt, T. K. Kim, L. Kipp, and K. Rossnagel, Phys. Rev. B 85, 224532 (2012).
  * Rodrigo and Vieira [2004] J. G. Rodrigo and S. Vieira, Physica C 404, 306 (2004).
  * Kuzmanović  _et al._ [2022] M. Kuzmanović, T. Dvir, D. LeBoeuf, S. Ilić, M. Haim, D. Möckli, S. Kramer, M. Khodas, M. Houzet, J. S. Meyer, M. Aprili, H. Steinberg, and C. H. L. Quay, Phys. Rev. B 106, 184514 (2022).
  * Liebhaber  _et al._ [2020] E. Liebhaber, S. Acero González, R. Baba, G. Reecht, B. W. Heinrich, S. Rohlf, K. Rossnagel, F. von Oppen, and K. J. Franke, Nano Lett. 20, 339 (2020).
  * Sanchez  _et al._ [1995] D. Sanchez, A. Junod, J. Muller, H. Berger, and F. Lévy, Physica B 204, 167 (1995).
  * Nader and Monceau [2014] A. Nader and P. Monceau, SpringerPlus 3, 16 (2014).
  * Mourachkine [2004] A. Mourachkine, Journal of Superconductivity 17, 711 (2004).
  * Xi  _et al._ [2016] X. Xi, Z. Wang, W. Zhao, J.-H. Park, K. T. Law, H. Berger, L. Forró, J. Shan, and K. F. Mak, Nat. Phys. 12, 139 (2016).
  * Ruby  _et al._ [2016] M. Ruby, Y. Peng, F. von Oppen, B. W. Heinrich, and K. J. Franke, Phys. Rev. Lett. 117, 186801 (2016).
  * Massee  _et al._ [2018] F. Massee, Q. Dong, A. Cavanna, Y. Jin, and M. Aprili, Rev. Sci. Instrum. 89, 093708 (2018).
  * Ruby  _et al._ [2015] M. Ruby, F. Pientka, Y. Peng, F. Von Oppen, B. W. Heinrich, and K. J. Franke, Phys. Rev. Lett. 115, 087001 (2015).