# TAP: Accelerating Large-Scale DNN Training Through
Tensor Automatic Parallelisation

Ziji Shi<EMAIL_ADDRESS>National University of Singapore  Alibaba Group
Le Jiang  Alibaba Group  Ang Wang  Alibaba Group  Jie Zhang  Alibaba Group
Xianyan Jia  Alibaba Group  Yong Li  Alibaba Group  Chencan Wu  Alibaba Group
Jialin Li  National University of Singapore  Wei Lin  Alibaba Group

###### Abstract

Model parallelism has become necessary to train large neural networks.
However, finding a suitable model parallel schedule for an arbitrary neural
network is a non-trivial task due to the exploding search space. In this work,
we present a model parallelism framework TAP that automatically searches for
the best data and tensor parallel schedules. Leveraging the key insight that a
neural network can be represented as a directed acyclic graph, within which
may only exist a limited set of frequent subgraphs, we design a graph pruning
algorithm to fold the search space efficiently. TAP runs at sub-linear
complexity concerning the neural network size. Experiments show that TAP is
$20\times-160\times$ faster than the state-of-the-art automatic parallelism
framework, and the performance of its discovered schedules is competitive with
the expert-engineered ones.

##  1 Introduction

111This paper was completed in October 2022.

Recent years have witnessed a burgeoning of large deep neural networks (DNNs)
that deliver unprecedented accuracy across a wide range of AI tasks. In fact,
the rate of DNN model size increase has far surpassed the growth in
accelerator memory capacity. As a result, model parallelism has been proposed,
where model weights are sharded onto multiple machines during distributed DNN
training, to address this memory capacity issue.

There are two main paradigms in model parallelism: pipeline parallelism and
tensor parallelism. Pipeline parallelism splits the model by layers. Only
activations are communicated during the forward pass, while gradient tensors
are exchanged in the backward phase. Pipeline parallelism has drawn much
attention recently; many proposed algorithms attempt to find the optimal
pipeline schedule that minimizes the pipeline idle time ("bubble size").
However, pipeline parallelism suffers from two significant drawbacks: 1) each
layer has to be small enough to fit into a single accelerator’s memory, and 2)
if the model has an imbalanced architecture, interleaving different layers can
be challenging. Tensor parallelism (or tensor sharding) is therefore proposed
as an alternative. Tensor parallelism partitions the model weights and
distributes them to multiple devices, lifting the restriction on single layer
size.

Unfortunately, manual specification of tensor parallelism can be challenging
and error-prone. The optimal partitioning scheme usually depends on both the
neural network architecture and the hardware system. Prior approaches tightly
couple user code, parallel strategies, and accelerator hardware, narrowing the
adaptability of model parallelism in practice (e.g., cloud environment).

Can an optimal tensor parallel plan be found automatically? Existing work on
automating model parallel training either rely on user hints or a brute-force
search over the entire space. [14, 26] incorporate user annotations as hints
to derive the pipeline or tensor sharding schedule, where users are required
to specify the number of pipeline stages or the mapping between layer to
physical device mesh. By inferring the full schedule based on expert
knowledge, it was made easier to write parallel plans than manual
specification [20, 24]. But the downside is that the manual specification may
not generalize across different hardware topologies. [15, 33, 30] propose to
search for the parallel schedule over a pre-defined search space. They
leverage Dynamic Programming to find the best device placement for the shards
of a neural network. However, the brute force search can be very slow in
reality, as they need to deal with exploding search space. Modern neural
networks can contain hundreds of thousands of operations, each having multiple
possibilities of sharding. The search for the optimal tensor parallel plan can
be reduced to the makespan problem, which has been proven to be NP-hard[16].

We argue that a brute-force search over the entire space is not necessary. We
observe that 1) most neural networks contain shared subgraphs, which can be
utilized to reduce search efforts when choosing parallel strategies, and 2)
communication is the main bottleneck when training in tensor parallelism, and
usually, it is impossible to overlap the contiguous partitions in a block
using tensor parallelism. Therefore, we can drastically accelerate the search
for better strategies by searching over unique neural network sub-modules
only, and we should evaluate the candidate strategies with communication cost.

Based on those observations, we designed a system drastically reducing search
space. Our contributions consist of the following:

  * •

A set of computational graph intermediate representations (IRs) that abstract
away from the implementation details of low-level operations;

  * •

A graph pruning algorithm that leverages the shared substructure to facilitate
efficient searching;

  * •

A communication-based cost model that accurately captures the amount of
communication necessary in tensor-parallel training.

We present TAP, an deep learning framework that automatically derives a
tensor-parallel plan for arbitrary neural networks without requiring expert
annotations. TAP drastically reduces the search time using shared subgraphs,
achieving $20\times-160\times$ speedup in finding a competitive sharding plan.
Evaluations show that our approach can find a comparable tensor parallel plan
similar to an expert-designed solution.

##  2 Background

Figure 1: Parallel strategies for distributed training

###  2.1 Data Parallelism

Data Parallelism is a commonly used parallel strategy that scales a model from
a single worker to multiple workers by replicating the weights, as shown in
Fig. 1b. Each worker maintains a full replica of the same model but trains on
different data slices. In the forward pass, the workers compute the parameter
updates independently on their own data slice, and different workers average
the gradients through AllReduce to ensure that the model parameters of all
workers are consistent during the backward pass.

###  2.2 Model Parallelism

Model parallelism is proposed to solve the downside of data parallelism: model
weight has to be able to fit into the memory of a single accelerator. Model
parallelism distributes the model weight into different devices and
synchronizes the full model through collective communication[7]. Model
parallelism can be divided into pipeline parallelism and tensor parallelism,
depending on the point of view.

####  2.2.1 Pipeline parallelism

Pipeline parallelism divides the model across the layers and distributes it to
different devices[13, 19, 10], shown in Fig. 1c. In the forward pass, the
training begins at GPU2, the data flows from GPU2 to GPU1 and finally reaches
GPU0; the backward pass starts on GPU0, calculates the gradients, and updates
the model parameters in the opposite direction to the forward pass. Pipeline
parallelism aims to minimize device idle time. Ideally, the size of each layer
should be similar so that each device’s workload is relatively balanced during
the training process. But in practice, the model architecture could be very
heterogeneous, and the inter-device communication speed could differ
significantly from inter-rack to intra-rack, thus greatly hindering the
training speed.

####  2.2.2 Tensor parallelism

Tensor parallelism splits the model layer and distributes it to different
devices to disperse the computational overhead of the layer[31, 27, 20], shown
in Fig. 1d. For each device, only part of the input tensors are stored in its
local memory, therefore the full result needs to be aggregated from partial
results on other devices through collective communication. Tensor parallelism
can alleviate the problem of training heterogeneous models using pipeline
parallelism, and achieve better performance.

###  2.3 Automatic Parallelism

Automatic parallelism is a recent line of research on automatically
distributing the local model from one device to multiple devices using the
data and model parallel strategies. Existing works on automatic parallelism
rely on user hints or brute-force searches over the entire space.

####  2.3.1 User hint

User-hint-based automatic parallelism is used to help users scale single-
device programs to multi-device. For example, GSPMD[31] infers the operator
partitioning scheme based on user annotations to scale single-device programs.
Whale[14] allows for incorporating user hints to perform semi-auto
parallelisation for large models and introduces a hardware-ware load balance
algorithm. However, user-hint-based automatic parallelism approaches require
users to have a deep understanding of both system and the model, and the hard-
coded user hints may not be transferable when the model or system changes.

####  2.3.2 Search algorithm

Recent work has proposed fully automatic approaches based on search algorithms
to optimize distributed DNN training. For example, Tofu[30] uses a recursive
search algorithm based on dynamic programming and DNN-specific heuristics to
minimize communication for the entire dataflow graph. Flexflow[15] uses
randomized search to find the best parallel strategy in the SOAP (Sample,
Operator, Attribute, and Parameter) space. Alpa[33] optimizes large DL models
through two-level optimizations: inter-operator and intra-operator. It
automates inter-operator parallelism by dynamic programming and intra-operator
parallelism by integer linear programming. Unity[28] represents
parallelisation and algebraic transformations as substitutions on a unified
graph representation, uses a novel hierarchical search algorithm to identify
an optimized sequence of substitutions, and scales to large numbers of GPUs
and complex DNNs.

####  2.3.3 Challenge of exploding search space

The search-based approaches face the challenge of exploding search space as
model size scales, resulting in a huge time cost. Concretely, each tensor
(assuming 2D) presents three possible sharding options: not sharding, sharding
on the first dimension (row-wise), or sharding on the second dimension
(column-wise). Given a neural network $G(E,V)$ with $V$ weight tensors, there
exists $3^{V}$ possible sharding plans. Therefore, there exists no polynomial
time solution to find an optimal sharding plan.

##  3 Approach

In this section, we formulate the problem of searching for an optimal tensor
parallel schedule, followed by our observation of the common presence of
shared sub-structures in a large neural network, leading to the motivation of
our design.

###  3.1 Problem Formulation

A neural network can be represented as a directed acyclic graph $G(E,V)$
comprised of $L$ layers. The set of vertices $V$ represents the operators, and
the set of edges $E$ represents the data flow from producer to consumer
operators. Operators can optionally carry a weight tensor. During the forward
pass, an edge represents an activation tensor, while in the backward phase, it
represents a gradient tensor. A layer $L_{i}\in L$ is either a layer or a
cluster of operators with a similar composition. Let the physical training
system be $S(m,n)$ where $m$ is the number of worker nodes, and $n$ is the
number of accelerators per worker node. A parallel plan $p$ is a new graph
mathematically equivalent to $G$, but the order of nodes may change, and the
communication node may be inserted into $p$. The cost function, $Cost(p,S)$,
measures training latency for a given plan and training system. The goal is to
find an optimal parallel plan $p^{*}$ where:

$\displaystyle\underset{p}{\text{minimize}}$ $\displaystyle Cost(p,S)$ subject
to $\displaystyle p(X)=G(X)\forall X$

How to find such a plan in an automated manner? Figure 2 illustrates the
typical workflow of an auto-parallel system. The system first attempts to
constrain the search space for splitting a model. With a more manageable
search space, a search algorithm then produces one or more candidate plans for
evaluation. All candidate plans are evaluated by a cost model that chooses the
best plan with the lowest cost based on its evaluation criteria.

Figure 2: General recipe of automatic model parallelism frameworks.

The end-to-end duration to produce an optimal schedule is a critical metric
for an auto-parallel system. We identify three main factors that contribute to
this overall completion time: the size of the search space, the time
complexity of the searching algorithm, and the speed of the evaluation method.

###  3.2 Challenges and Observations

As we see earlier, a major challenge faced by auto-parallel systems is the
search space explosion problem. This exponential increase in candidate space
has led to impractical search time to derive a parallel schedule for modern
large models [33] (§ 6.3). It creates a dilemma: auto-parallel systems aim to
accelerate large model training; however if the derivation step itself is too
slow, it may offset the benefit of using an auto-parallel system.

| Scaling
---
Technique
Task | Model | # Params | Shared Subgraph (SS) | # of SS
By width | Vision | ResNet50[12] | 23M | Conv | $50\times$
Vision + Language | CLIP-Base[21] | 63M | Transformer | $12\times$
Language Model | WideNet[32] | 63M | MoE layer | $32\times$
Vision | ViT-Huge[9] | 632M | Transformer | $32\times$
Vision | V-MoE[25] | 15B | MoE layer | $24\times$
By depth | Speech | wav2vec 2.0[4] | 317M | Conv, Transformer | $7\times,24\times$
Language Model | BERT[8] | 340M | Transformer | $24\times$
Language Model | T5-Large[22] | 770M | Transformer | $24\times$
Language Model | GPT-3[5] | 175B | Transformer | $96\times$
Language Model | Switch Transformer[11] | 1571B | MoE layer | $15\times$
Table 1: Shared subgraphs exist on many neural network models. "Conv" means
convolutional layer, "MoE" means Mixture-of-Expert layer.

How to effectively reduce this large candidate search space? To answer this
question, we studied common scaling techniques for popular DNN models and
summarized our findings in Table 1. We conclude that these techniques can be
grouped into two major categories: scaling on the width by increasing the
dimension of layers (e.g., adding the number of classes, adding attention
heads, increasing the convolutional channels), or scaling on the depth by
increasing the number of layers. Our key observation is that both techniques
start with a base subgraph, i.e., a group of layers or operators, and expand
from it. For instance, large-scale pre-trained language models such as BERT[8]
and T5[22] consist of tens of transformer layers, and multi-class object
classification networks like ResNet-50[12] are made of convolutional layers.

Furthermore, by analyzing expert-engineered parallel schedules ([20, 23, 24]),
we observe that parallel schedules are primarily identical for the same type
of layers. The underlying reason is that similar layers share a similar amount
of computation and memory consumption. This has motivated us to explore the
possibilities of reusing the parallel schedules discovered for the same layer
to save search effort.

###  3.3 Motivating Examples

Figure 3: Motivating examples.

We are motivated by two commonly encountered scenarios as depicted in Fig. 3:
scaling along the width (a), or scaling along the depth (b). In an e-commerce
setting, there exist hundreds of thousands to millions of types of
merchandise. Therefore, to classify a product image, a classification model
like ResNet[12] must have a very wide fully connected (FC) layer. As shown in
Fig. 3a, the size of the FC layer (205M floating point numbers) can be
disproportionately larger than the feature extraction module (24M). In fact,
the single gigantically wide layer can be too large to fit into an
accelerator, therefore, preventing the use of pipeline parallelism.

Another scenario is the family of models built using dense transformer
architecture as shown in Fig. 3b. Transformer[29] has been the building block
for many large language models (Table 1). A typical transformer layer
comprises a self-attention and a multi-level perception (MLP) layer. The self-
attention layer has 4 weight tensors Q, K, V, W, and the MLP layer has 2
weight tensors: one for the intermediate layer and one for the output layer. A
dense transformer model usually scales up by stacking transformer layers [8,
22, 5, 9]. Due to the similarity in model architecture, we may reuse the
sharding plan found for one layer on all transformer layers[20].

Under both scenarios, tensor parallelism proves to be a more general solution,
which has motivated us to formulate the problem under a unified view.

###  3.4 Split-Replica-Communication (SRC) Abstraction

Inspired by the observation that heavy bulky matrix multiplications of tensors
drive model parallelism, TAP uses a unified view to represent a data and
tensor parallelism schedule using the SRC abstraction.

##### Split(axis).

Split means sharding the tensor on a target axis, and different device stores
different partitions. Under this view, data parallelism is a special case for
tensor parallelism where the tensor shards on the first dimension (batch
dimension).

##### Replica.

Replica means to replicate the tensor on different devices without sharding.
For example, in data parallelism, the model tensors are replicated while the
input tensors are sharded.

##### Communication.

Splitting results in partial values. Also, the back-propagation requires
aggregating different copies of gradients. Therefore, additional communication
operators may be required to combine the partial results to ensure
mathematical equivalence.

It is worth noting that we are not required to use all S&R&C. We may have
patterns expressed as R-only, S&C, R&C, or S&R&C. Under the SRC abstraction,
the original operator expression

$Y=Op(A,B)$

will be converted to a distributed version by

$SR(Y)=Comm(Op(SR(A),SR(B))$

With the SRC abstraction, we can define general sharding rules for each
operator. They are implemented as ShardingPatterns in TAP . If there is no
viable way to split, we can always fall back to replicating the tensors (data
parallel).

Fig. 4 illustrate a matrix multiplication (MatMul) pattern in TAP under SRC
abstraction. $X$ and $W$ are 2D matrices, representing input tensor and weight
tensor respectively. $X$ is sharded column-wise while $W$ is sharded row-wise.
$X_{0}$ and $W_{0}$ are stored on device 0. Notice that standard MatMul still
works for the sharded tensors, hence $Y_{0}$ shares the same shape as $Y$.
However, $Y_{0}$ is just a partial result, and it needs to be added with
$Y_{1}$ from device 1 to get $Y$. As such, an AllReduceSum communication is
required to sum them up for mathematical equivalence.

Figure 4: Example using SRC abstraction to perform MatMul on two sharded
tensors. Only S&C were used.

##  4 Design and Implementation

Figure 5: Overview of the TAP system.

###  4.1 Overview

As depicted in Fig. 5, given a neural network represented as a graph, TAP
first converts the graph into its intermediate representation(§ 4.2) called
GraphNode and removes auxiliary nodes. Then TAP performs graph pruning(§ 4.3),
restricting the search space from the complete graph to the subgraphs instead.
After pruning, TAP explores the possible sharding opportunities using pre-
defined sharding patterns(§ 4.4) and validates the candidate plans(§ 4.5). If
a valid plan is found, it is then evaluated using the cost model(§ 4.6). In
the end, TAP takes the best plan, performs additional optimizations on
communication, and rewrites the model into a parallel version(§ 4.7.1).

Example 1: Auto-parallel code with TAP on 2 workers each with 8 GPUs

⬇

1import tensor_auto_parallel as tap

2mesh = [2, 8]

3tap.auto_parallel(tap.split(mesh))

4model_def()

###  4.2 Intermediate Representation

TAP defines a family of high-level Intermediate Representations (IRs) to
facilitate the derivation of parallel schedules. Compared to MLIR HLO[17], TAP
IRs operate on a coarser granularity while preserving the necessary
information for sharding.

Upon obtaining the original neural network graph, TAP first trims the graph by
deleting the auxiliary operators (Step 1 in Fig. 5). This will remove the
initialization and checkpoint-related operators, which will be recovered when
converted back to a neural network graph later. As a result, the remaining
graph will consist of only computing and communication operators.

TAP IRs consists of:

##### GraphNode.

A GraphNode represents a group of computing or communication operators. It can
be a layer or a logical group of operators, which is the basic unit for
deriving the sharding schedule. The TAP graph is made of GraphNode while
preserving the directed edges from the original DAG. Using the GraphNode IR,
we reduce the number of nodes in the T5-large model from 60k to 1015 weight
variables.

##### Sharding pattern.

A GraphNode could have multiple ways of sharding. For instance, a 2D matrix
weight can be split on either dimension or replicated. TAP defines each
sharding pattern using the SRC abstraction. TAP also establishes the cost of
each sharding pattern based on communication cost.

##### Sharding plan.

A sharding plan is a set of subgraphs (blocks of GraphNodes) with sharding
patterns connecting them.

###  4.3 Pruning using Shared Subgraph

It is common for DNN models to contain shared subgraphs. If we could identify
the shared subgraphs, we could prune the search space by searching only within
the subgraph. We propose a graph pruning algorithm to compress the search
space into a shared structure (Step 2):

Algorithm 1 Graph Pruning

1:procedure PruneGraph($modelDef,minDuplicate$)

2: $nodeTree\leftarrow\emptyset$

3: $maxDepth\leftarrow modelDef.depth$

4: for all $depth\in maxDepth\cdots 1$ do

5: $nodeTree[depth]\leftarrow longestCommonPrefix(modelDef.nodes.name)$

6: $opCount=findSimilarBlk(nodeTree[depth])$

7: if $opCount\geq minDuplicate$ then

8: $subgraphs.append(nodeTree[depth])$

9: else

10: $break$

11: end if

12: end for

13: $\textbf{return}subgraphs$

14:end procedure

In deep learning frameworks like TensorFlow[3], each variable is referred to
by the operator that produces it. As such, variables under the same layer
share the same name scope because they receive input from the same operator.
Therefore, it is possible to cluster operators that fall under the same name
scope.

Algorithm 1 starts by constructing a nodeTree, which identifies and groups the
GraphNodes on each level by using the longest common prefix algorithm on the
GraphNodes names (line2-5). After that, it finds the blocks of GraphNodes with
a similar composition of operators and compares the number of operators with
the minimum duplicate threshold (line 7). As the depth decreases, we will see
a larger subgraph with less homogeneous compositions. Notice that multiple
shared subgraphs may exist since a neural network may have multiple leaf
nodes.

###  4.4 Sharding Plan Generator

Algorithm 2 Derivation of Optimal Plan

1:procedure DerivePlan($modelDef,shardingPatterns$)

2: $subgraphs\leftarrow PruneGraph(modelDef)$

3: $candidatePlans\leftarrow enumerateAllPlans(subgraphs)$

4: $validPlans\leftarrow\\{\\}$

5: for all $p\in candidatePlans$ do

6: $validated\leftarrow PatternRouting(p)$

7: if $validated$ then

8: $validPlans.insert(p)$

9: end if

10: end for

11: $bestPlan\leftarrow min(QueryCost(validPlans))$

12: $\textbf{return}bestPlan$

13:end procedure

A sharding pattern, defining the way a GraphNode can be sharded, also serves
as the directed edge between nodes. According to the SRC abstractions, the
communication pattern is determined once the split/replica decision is made.
Under the hood, the sharding patterns connect to each other like a chain.

After pruning, TAP proceeds to derive the optimal plan (Step 3 and 4) using
Algorithm 2. In the first phase, TAP enumerates all possible sharding plans
given the subgraphs. TAP only needs to work on hundreds of plans thanks to
pruning. However, not every plan is valid because we only have weekly
connected subgraphs. Therefore, the candidate plans need to be validated by
checking the connectivity (line 5-10). Upon checking, TAP evaluates the
performance of each plan using a cost model and selects the best plan.

###  4.5 Pattern Routing

Algorithm 3 Plan Validation

1:procedure PatternRouting($currPlan$)

2: $TopoSort(currPlan)$

3: $nodesQ\leftarrow currPlan.root$

4: while $nodesQ\neq\emptyset$ do

5: $currNode\leftarrow nodesQ.dequeue()$

6: for all $childNode\in currNode.next()$ do

7: $sp\leftarrow lookUpShrdPatn(currNode,childNode)$

8: if $sp\neq\emptyset$ then

9: if $childNode==currPlan.leaf$ then

10: return TRUE

11: else

12: $nodeQ.enqueue(childNode)$

13: end if

14: end if

15: end for

16: end while

17: return FALSE

18:end procedure

In the pattern routing step (Algorithm 3), TAP tries to assemble the weakly
connected GraphNodes into a valid sharding plan by checking the
connectivities. This is to ensure the success of graph rewriting (Step 5). TAP
does so using breadth-first-search (BFS) starting from the root node, and the
goal is to make sure there exists at least a connected path from the root to
the leaf chained using the sharding patterns.

One challenge is that a pair of contracting sharding patterns may have
different input and output tensors, and a consumer operator’s input is not
ready until its producer is ready. In other words, dependencies exist between
GraphNodes, but the information was kept in the original edges and could be
lost when we perform pruning.

To solve it, we perform a topological search for the GraphNode based on the
readiness of the input tensor. We leverage that neural networks can be
represented using a directed acyclic graph, and reconstruct the edges based on
the order of the producer-consumer relationship. This way, TAP avoids checking
the order for every pair of GraphNodes.

###  4.6 Cost Model

To build a cost model, we first profile different tensor parallel plans to
understand the bottleneck. Fig. 6 summarizes the result. Data were collected
from two nodes interconnected by 32 Gbps Ethernet, each equipped with 8 GPUs.
We observe that inter-node communication is the main bottleneck for tensor
parallelism, and the best plan is not necessarily the one that splits every
weight tensor, in line with [7].

Figure 6: Time breakdown for tensor parallel plans on T5-large model on 8 and
16 GPUs (8w/16w). DP means data parallel, MHA means sharding the multi-head
attention, FFN means sharding the feed-forward layer, and Megatron refers to
the tensor sharding plan described in [20].

As the number of devices increases from $8\times$ to $16\times$, the
difference between communication time and computation time is further
pronounced. This is because the bottleneck has shifted from high-speed intra-
node communication (PCI-e) to slower inter-node communication (Ethernet).

Furthermore, the best tensor parallel plan for 16 GPUs (16w-FFN) only shards
the weight in the feed-forward layer. We conjecture that with more tensors
split instead of replicated, there are fewer FLOPs per device and the
computation time is lower. However, this comes at the cost of having more
communication. In the case of training in the data center where nodes are
interconnected by Ethernet, the speed bottleneck may shift from computation to
communication instead. Therefore, communication cost is the main consideration
when we design the cost model.

Motivated by the two observations above, TAP adopts the amount of
communication as the cost of each plan. But calculating the correct amount of
communication is not straightforward, as there exists a few practical
challenges as below:

##### Counting communicated parameter.

Although the amount of FLOPs is mostly symmetrical between the forward and
backward phases, their communication patterns differ. During the forward
phase, all activation tensors will be communicated, yet DNN frameworks will
only communicate non-constant parameters during the backward phase. Therefore,
the TAP cost model should be able to identify the weight that requires
communication.

##### Gradient overlap/aggregation.

The communication/computation patterns are different in the forward and
backward phases. During the forward phase, the computation of the current
layer is blocked until the input arrives from the previous layer, creating the
temporal dependency. However, the communication during the backward phase can
be overlapped with computation because the weight update stage is independent
of gradient communication.

Furthermore, a large fraction of weights in a neural network are small tensors
of tens to hundreds of bytes, but collectively they cause a large number of
communications. To reduce the overhead of sending many small gradient packets,
TAP fuses multiple gradients into one and overlaps with computation during the
graph rewriting phase(§ 4.7.1). As a result, the amount of gradient
communication may not fully translate into communication time in the backward
phase.

##### Efficiency of collective communications.

We observe that the communication efficiencies are different for different
collective communication. For instance, we observe that AlltoAll and AllGather
take more time to communicate the same amount of messages compared to the
heavily optimized AllReduce in NCCL.

TAP addresses these issues using an analytical cost model based on the
tensor’s communication method, shape, and data format. Each sharding pattern
is associated with a cost, and the total cost is calculated by summing all
pattern costs along the critical path.

###  4.7 Graph Rewriting

After evaluating the cost of each sharding plan, TAP assembles the parallel
plan. It does so by first restoring the original order of operators. Then, TAP
identifies optimization opportunities that can be performed through gradient
packing. In the end, TAP passes the resulting parallelized neural network plan
to the deep learning framework runtime.

####  4.7.1 Gradient packing.

TAP optimizes communication by packing multiple small gradient updates into a
larger one during the gradient synchronization stage, saving the communicator
setup overhead. During the backward pass, each layer produces gradients that
will be synchronized with other workers. This can be very time-consuming, as
the number of gradient packets sent equals the total number of trainable
parameters in a neural network.

As a communication optimization technique, TAP fuses packets smaller than a
threshold $\mu$ into one larger one. To prevent the aggregated packet from
growing too large and deferring the weight update stage, TAP segments the
packets into equally sized chunks. Therefore, the gradient synchronization and
weight update stages can be pipelined, allowing the former to transmit while
updating the weight parameters.

###  4.8 Limitation and Future Work

To further optimize the memory consumption, TAP could leverage other
orthogonal techniques such as Auto Mixed Precision (AMP) [1], recomputation
[6], and pipeline parallelism. Since both AMP and TAP optimize on the graph
representation of the neural network, they can be made into different passes.
Also, gradient checkpointing can be used to offload the selected GraphNode
onto the main memory. TAP may also be used with pipeline parallelism through
automatic[19, 18, 10, 14] or manual placements.

##  5 Complexity Analysis

This section presents an in-depth analysis of the arithmetic complexities of
known algorithms in auto-parallelism.

###  5.1 Complexity Analysis of Existing Solutions

Following the discussion on related work, we analyze the complexities of two
automatic model parallel frameworks and present it in Table 2. We define the
total complexity as:

$\begin{split}total\\_complexity={}&search\\_complexity\\\
&+num\\_plans*evaluation\\_complexity\\\ \end{split}$

Framework | Search Space | Search Algorithm | Evaluation | Total
---|---|---|---|---
FlexFlow | $N(4E,4V)$ | $B$ | $O(V+E)$ | $O(BV+BE)$
Alpa | $N(kE,kV)$ |  | Inter-Op: $O(V^{2}L)$
---
Intra-Op: $O(E(V+E))$
$O(V+E)$ | $O(V^{2}L(V+E^{2}))$
TAP | $N(\frac{E}{2CL},\frac{V}{2CL})$ | $O(\frac{E+V}{L})$ | $O(\frac{E}{L})$ | $O(\frac{E+V}{L})$
Table 2: Complexities of selected auto model parallel frameworks.

####  5.1.1 FlexFlow.

FlexFlow operates on four dimensions (S/O/A/P), and there was no space
reduction. Therefore, the search space is $N(4E,4V)$. As search complexity,
FlexFlow employs the Markov chain Monte Carlo (MCMC) algorithm. Thus, we use
$B$ to denote the number of trials in MCMC sampling. Furthermore, within each
trial, it needs to evaluate its performance by querying the cost model with
Depth-First-Search(DFS), hence its evaluation complexity is $O(V+E)$.

####  5.1.2 Alpa.

Alpa is formulated as a multi-level optimization problem: in the outer loop,
it searches for the inter-op plan using dynamic programming; in the inner
loop, it finds the intra-op parallel plan using integer linear programming.
First, since it operates at MLIR HLO, which is a finer IR than the TensorFlow
operator, we formulate the search space as $N(kE,kV)$ where $k\geq 1$. In the
outer loops, it uses a similar algorithm to [18] to search for pipeline slices
and map the slices to devise mesh. Optimization like operator clustering and
early pruning reduces the outer loop complexity to $(kV)^{2}L$. For the inner
loop, since the exact complexity of their ILP solver is unknown, we use a
lower bound by performing a BFS from each operator, and the complexity is
given as $kE(kV+kE)$. Finally, each trial needs to evaluate its performance by
querying the cost model, so the evaluation complexity is $kV+kE$.

####  5.1.3 TAP .

In TAP , we first reduce the search space by converting the TensorFlow graph
to TAP graph (by $C\times$, where $C\geq 1$). We then prune the tree by layer,
further reducing the complexity to $N(\frac{E}{2CL},\frac{V}{2CL})$. In the
searching stage, the result is derived by performing a BFS. Thus, the
complexity is $\frac{V+E}{2CL}$. For the evaluation stage, TAP needs to
evaluate the cost of each plan by querying the cost model, which depends only
on the size of the edges. Thus, the evaluation cost is $\frac{E}{2CL}$.

##  6 Evaluation

###  6.1 Setup

We first evaluate the pruning algorithm and the use of Just-In-Time
compilation for TAP . Then, for comparison with another auto-parallel
framework, we use Alpa version 0.7 running with JAX 0.3.5. Next, we use
Megatron running on PyTorch to compare against expert-engineered tensor
parallel plans. Finally, we present the training convergence running gigantic
neural networks.

The evaluation was performed on Company A’s public cloud node with 756GB main
memory, $2\times$ Intel 8163 CPUs at 24 cores each, and $8\times$ Nvidia V100
SXM2 32GB GPUs. Additionally, TAP builds on top of TensorFlow 1.12.

###  6.2 Micro Benchmarks

####  6.2.1 Pruning Algorithm

Figure 7: Tuning minDuplicates for Algorithm 1

One of the key contributions of TAP is the graph pruning algorithm, in which
minDuplicates determines the minimum size of subgraphs. If the threshold for
subgraphs is too low, we may still face exploding search spaces; if the
threshold is too high, we may see too few subgraphs, resulting in a longer
search time. Since the architecture of different neural networks may vary
significantly, it is desirable to have a robust threshold.

We explore a range of minDuplicates and report the number of unique subgraphs
found and graph pruning algorithm runtime in Fig. 7. We conclude that the
threshold is robust and does not require significant effort in tuning. Take
the T5-large model with 770M parameters as an example. When the threshold is
1, meaning the graph is unpruned, it contains 6561 nodes. After pruning, the
number has been drastically shrunk to just 5. As the threshold changes, the
number of identified unique blocks stays relatively stable, showing that our
graph pruning algorithm is insensitive to different thresholds.

Furthermore, we observe that the pruning algorithm is very efficient, taking
less than 12 seconds to find the subgraphs for T5-large, and less than a
second for the 152-layer 100K-class ResNet model, which proves the scalability
of TAP’s graph pruning algorithm.

####  6.2.2 XLA

Figure 8: Training time per iteration when XLA is enabled. Lower is faster.

XLA[2] is a JIT compiler for DNN frameworks, optimizing the training mainly
through fusing smaller kernels to reduce launch overhead. Like TAP , XLA
identifies the connected subgraphs and optimizes on the operator level. We
evaluate the time per iteration with and without XLA on the ResNet50 model
with varying numbers of classes. Fig. 8 shows that the improvement from XLA is
not consistent, and we observe a similar trend in T5 models, which has
performance improvement between $-9\%$ to $+1\%$. We believe the inconsistent
performance improvement results from new communication nodes being inserted
into the parallelised plans by TAP . Therefore, XLA may have difficulty
identifying the correct cluster of operators to fuse. Furthermore, XLA’s
operator clustering may hinder the degree of communication and computation
overlap, affecting the scaling efficiency. For this reason, we did not enable
XLA for the rest of the experiments.

###  6.3 End-to-End Evaluation

In this section, we compare with auto-parallel framework Alpa on search time
and performance of the discovered plan.

####  6.3.1 Search time.

As explained in § 5.1, TAP has a sub-linear time complexity, which is
desirable when the models’ size scales up. In the experiments with Alpa, we
present the end-to-end search time with respect to model scaling, defined by
the duration from the start of the experiment till the moment that the
training process begins. Due to time constraints, we shortlisted a search
space of 16 plans for T5 and 5 plans for ResNet, while we did not restrict the
search space for TAP .

Figure 9: End-to-end search time when scaling on the number of parameters for
dense transformer model. Figure 10: End-to-end search time when scaling on the
number of parameters for the large-scale classification model.

To scale the model along the depth, we increase the number of transformer
layers for T5, an encoder-decoder transformer architecture for language
modeling. Increasing the depth of dense transformer models is a common
practice to improve performance. Fig. 9 shows that, with rising parameters,
TAP can still find a plausible schedule in under 15 mins, which is
$21\times-67\times$ faster than Alpa.

To scale the model size along the width for the ResNet50 model, we choose to
increase the size of the classification layer. The original ResNet50 model has
1024 classes in the classification layer. As we increase the dimensions for
the classification layer, the total number of parameters also scales up. As
shown in Fig. 10, TAP is two orders of magnitude faster than Alpa in finding
the optimal solution. Our systems outperforms it by $103\times-162\times$.

We further analyze the time breakdown during the search. For example, for
24-layer T5-large (770M parameters), Alpa spent 5 mins profiling the operators
and 5 mins constructing the pipeline stages out of the operators. Instead, TAP
reduces the architecture to one transformer block and searches for shardable
parameters within that only, drastically reducing the search space. As a
result, Alpa takes 197 minutes to search for 16 candidate plans, while TAP
requires only 6 minutes to examine 729 candidate plans.

####  6.3.2 Training speed.

Figure 11: Training time per iteration for T5 (batch size=16). The blue band
represents the standard derivation.  Figure 12: Training time per iteration
for ResNet50 (batch size=1024).

We also evaluate the performance of the best plans produced by Alpa and TAP .
We observe that Alpa favors pipeline parallel schedules, while the optimal
schedule found by TAP is similar to the Megatron-style tensor parallel
schedule. Since the plans using pipeline parallelism require less
communication, the plans from Alpa have a higher throughput.

We also observe that as the width of the model increases, the performance of
TAP plans is better and more consistent. Fig. 12 shows the time to finish one
iteration of training for parallel plans of ResNet50. We first observe that
TAP consistently outperforms Alpa. Further, the variance (blue band) in plans
discovered by Alpa shows that it struggles to find consistently good plans.
222TAP only outputs the best plan out of all possible plans. Therefore it only
has one line.

###  6.4 Evaluation of Optimal Sharding Plan

We compare the best parallel plan found by TAP with an expert-engineered plan
described in [20]. Megatron runs on PyTorch.

Figure 13: Comparison between TAP and Megatron plan.

####  6.4.1 Memory and training speed.

We observe from Fig. 13 that the best parallel plan found by TAP is more
memory efficient than Megatron, making it more scalable for large neural
network models since the memory capacity is usually the bottleneck.
Furthermore, regarding training speed, TAP’s best plan is only within 2.3% to
14.8% slower than the Megatron plan.

####  6.4.2 Visualization of discovered plans.

Figure 14: Selected sharding plans discovered by TAP for T5. Each box
represents a trainable variable in the transformer layer.

We plot some of the sharding plans discovered by TAP in Fig. 14. For dense
transformer models, we observe that TAP usually shards the weight variable in
the multi-head attention layer while keeping the embedding and layernorm
layers replicated. Therefore, we conjecture that the attention and FFN layers
are usually more parameter-heavy, making them more suitable for weight
sharding.

We also found that TAP was not only able to discover Megatron-style fully
sharded plans or data-parallel plans but also can it find partially sharded
plans that only split multi-head attention (MHA-only) or feed-forward layers
(FFN-only). To our surprise, the best plan found by TAP on the experiment
system is the FFN-only plan, where the multi-head attention (green) gets
replicated, and the feed-forward layer (yellow) is sharded. Unlike fully
sharded plans like Megatron, the FFN-only plan will make more efficient use of
the available GPU memory and save commutation when GPU resource is abundant so
that it does not need to aggregate the partial activations in the forward
pass.

###  6.5 Scaling beyond Single Worker

To push the boundary of the scalability, we train the 100 billion parameters
M6-MoE-100B model with 128 NVIDIA V100 GPUs and 1 trillion parameter M6-MoE-1T
with 480 NVIDIA V100 GPUs. We scale model parameters by ten times while only
increasing GPU count by 3.75 times. Besides the resource saved per parameter,
M6-MoE-1T showed a significant model quality gain compared to M6-MoE-100B, as
shown in Fig. 15.

Figure 15: Training loss of M6-MoE-100B and M6-MoE-1T.

##  7 Conclusion

We present TAP, an automatic parallelism framework that efficiently discovers
data/tensor parallel plans for large models. Leveraging the observation that
shared subgraphs widely exist in neural networks, we design a pruning
algorithm and SRC abstraction that efficiently reduces the search space with a
sub-linear end-to-end complexity. The best plans found by TAP are comparable
with the state-of-the-art expert-engineered plans while only taking minutes to
discover. TAP will be open-sourced on GitHub.

## References

  * [1] Automatic mixed precision for deep learning.  https://developer.nvidia.com/automatic-mixed-precision.
  * [2] Xla: Optimizing compiler for machine learning.  https://www.tensorflow.org/xla.
  * [3] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dan Mane, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viegas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-Scale Machine Learning on Heterogeneous Distributed Systems.  2016\.
  * [4] Alexei Baevski, Yuhao Zhou, Abdelrahman Mohamed, and Michael Auli.  wav2vec 2.0: A framework for self-supervised learning of speech representations.  Advances in Neural Information Processing Systems, 33:12449–12460, 2020.
  * [5] Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al.  Language models are few-shot learners.  Advances in neural information processing systems, 33:1877–1901, 2020.
  * [6] Tianqi Chen, Bing Xu, Chiyuan Zhang, and Carlos Guestrin.  Training deep nets with sublinear memory cost.  arXiv preprint arXiv:1604.06174, 2016.
  * [7] Jeffrey Dean, Greg S. Corrado, Rajat Monga, Kai Chen, Matthieu Devin, Quoc V. Le, Mark Z. Mao, Marc Aurelio Ranzato, Andrew Senior, Paul Tucker, Ke Yang, and Andrew Y. Ng.  Large scale distributed deep networks.  Technical report, 2012.
  * [8] Jacob Devlin, Ming Wei Chang, Kenton Lee, and Kristina Toutanova.  BERT: Pre-training of deep bidirectional transformers for language understanding.  Technical report, 2019.
  * [9] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  arXiv preprint arXiv:2010.11929, 2020.
  * [10] Shiqing Fan, Yi Rong, Chen Meng, Zongyan Cao, Siyu Wang, Zhen Zheng, Chuan Wu, Guoping Long, Jun Yang, Lixue Xia, Lansong Diao, Xiaoyong Liu, and Wei Lin.  DAPPLE: A pipelined data parallel approach for training large models.  Proceedings of the ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming, PPOPP, 21:431–445, 2021.
  * [11] William Fedus, Barret Zoph, and Noam Shazeer.  Switch transformers: Scaling to trillion parameter models with simple and efficient sparsity, 2021.
  * [12] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, volume 2016-Decem, pages 770–778, 2016\.
  * [13] Yanping Huang, Youlong Cheng, Ankur Bapna, Orhan Firat, Mia Xu Chen, Dehao Chen, Hyouk Joong Lee, Jiquan Ngiam, Quoc V. Le, Yonghui Wu, and Zhifeng Chen.  GPipe: Efficient training of giant neural networks using pipeline parallelism.  Advances in Neural Information Processing Systems, 32, 2019.
  * [14] Xianyan Jia, Le Jiang, Ang Wang, Wencong Xiao, Ziji Shi, Jie Zhang, Xinyuan Li, Langshi Chen, Yong Li, Zhen Zheng, Xiaoyong Liu, and Wei Lin.  Whale: Efficient giant model training over heterogeneous gpus.  In USENIX Annual Technical Conference. USENIX, 2022.
  * [15] Zhihao Jia, Matei Zaharia, and Alex Aiken.  Beyond Data and Model Parallelism for Deep Neural Networks.  arXiv, 2018.
  * [16] Shui Lam and Ravi Sethi.  Worst case analysis of two scheduling algorithms.  SIAM Journal on Computing, 6(3):518–536, 1977.
  * [17] Chris Lattner, Mehdi Amini, Uday Bondhugula, Albert Cohen, Andy Davis, Jacques Pienaar, River Riddle, Tatiana Shpeisman, Nicolas Vasilache, and Oleksandr Zinenko.  MLIR: Scaling compiler infrastructure for domain specific computation.  In 2021 IEEE/ACM International Symposium on Code Generation and Optimization (CGO), pages 2–14, 2021.
  * [18] Zhuohan Li, Siyuan Zhuang, Shiyuan Guo, Danyang Zhuo, Hao Zhang, Dawn Song, and Ion Stoica.  TeraPipe: Token-Level Pipeline Parallelism for Training Large-Scale Language Models.  2021\.
  * [19] Deepak Narayanan, Aaron Harlap, Amar Phanishayee, Vivek Seshadri, Nikhil R. Devanur, Gregory R. Ganger, Phillip B. Gibbons, and Matei Zaharia.  Pipedream: Generalized pipeline parallelism for DNN training.  SOSP 2019 - Proceedings of the 27th ACM Symposium on Operating Systems Principles, pages 1–15, 2019.
  * [20] Deepak Narayanan, Mohammad Shoeybi, Jared Casper, Patrick LeGresley, Mostofa Patwary, Vijay Korthikanti, Dmitri Vainbrand, Prethvi Kashinkunti, Julie Bernauer, Bryan Catanzaro, Amar Phanishayee, and Matei Zaharia.  Efficient Large-Scale Language Model Training on GPU Clusters Using Megatron-LM.  International Conference for High Performance Computing, Networking, Storage and Analysis, SC, 2021.
  * [21] Alec Radford, Jong Wook Kim, Chris Hallacy, Aditya Ramesh, Gabriel Goh, Sandhini Agarwal, Girish Sastry, Amanda Askell, Pamela Mishkin, Jack Clark, et al.  Learning transferable visual models from natural language supervision.  In International Conference on Machine Learning, pages 8748–8763. PMLR, 2021.
  * [22] Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J. Liu.  Exploring the limits of transfer learning with a unified text-to-text transformer.  Technical report, 2020.
  * [23] Samyam Rajbhandari, Jeff Rasley, Olatunji Ruwase, and Yuxiong He.  Zero: Memory optimizations toward training trillion parameter models.  International Conference for High Performance Computing, Networking, Storage and Analysis, SC, 2020-Novem, 2020.
  * [24] Jie Ren, Samyam Rajbhandari, Reza Yazdani Aminabadi, Olatunji Ruwase, Shuangyan Yang, Minjia Zhang, Dong Li, and Yuxiong He.  ZeRO-offload: Democratizing billion-scale model training.  2021 USENIX Annual Technical Conference, pages 551–564, 2021.
  * [25] Carlos Riquelme, Joan Puigcerver, Basil Mustafa, Maxim Neumann, Rodolphe Jenatton, André Susano Pinto, Daniel Keysers, and Neil Houlsby.  Scaling vision with sparse mixture of experts.  Advances in Neural Information Processing Systems, 34:8583–8595, 2021.
  * [26] Noam Shazeer, Youlong Cheng, Niki Parmar, Dustin Tran, Ashish Vaswani, Penporn Koanantakool, Peter Hawkins, Hyouk Joong Lee, Mingsheng Hong, Cliff Young, Ryan Sepassi, and Blake Hechtman.  Mesh-tensorflow: Deep learning for supercomputers.  Advances in Neural Information Processing Systems, 2018-Decem:10414–10423, 2018.
  * [27] Mohammad Shoeybi, Mostofa Patwary, Raul Puri, Patrick LeGresley, Jared Casper, and Bryan Catanzaro.  Megatron-LM: Training Multi-Billion Parameter Language Models Using Model Parallelism.  2019\.
  * [28] Colin Unger, Zhihao Jia, Wei Wu, Sina Lin, Mandeep Baines, Carlos Efrain, Quintero Narvaez, Vinay Ramakrishnaiah, Nirmal Prajapati, Pat Mccormick, Jamaludin Mohd-yusof, Jongsoo Park, Misha Smelyanskiy, Alex Aiken, Pat Mccormick, Jamaludin Mohd-yusof Xi, and Luo Dheevatsa.  Unity : Accelerating DNN Training Through Joint Optimization of Algebraic Transformations and Parallelization This paper is included in the Proceedings of the.  2022\.
  * [29] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in Neural Information Processing Systems, 2017-Decem:5999–6009, 2017.
  * [30] Minjie Wang, Chien chin Huang, and Jinyang Li.  Supporting very large models using automatic dataflow graph partitioning.  Proceedings of the 14th EuroSys Conference 2019, 2019.
  * [31] Yuanzhong Xu, HyoukJoong Lee, Dehao Chen, Blake Hechtman, Yanping Huang, Rahul Joshi, Maxim Krikun, Dmitry Lepikhin, Andy Ly, Marcello Maggioni, Ruoming Pang, Noam Shazeer, Shibo Wang, Tao Wang, Yonghui Wu, and Zhifeng Chen.  GSPMD: General and Scalable Parallelization for ML Computation Graphs.  2021\.
  * [32] Fuzhao Xue, Ziji Shi, Futao Wei, Yuxuan Lou, Yong Liu, and Yang You.  Go wider instead of deeper.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 8779–8787, 2022.
  * [33] Lianmin Zheng, Zhuohan Li, Hao Zhang, Yonghao Zhuang, Zhifeng Chen, Yanping Huang, Yida Wang, Yuanzhong Xu, Danyang Zhuo, Joseph E. Gonzalez, and Ion Stoica.  Alpa: Automating Inter- and Intra-Operator Parallelism for Distributed Deep Learning.  2022\.


# Two molecular devices for superconducting spintronics

Cristina Mier  Centro de Física de Materiales CFM/MPC (CSIC-UPV/EHU), 20018
Donostia-San Sebastián, Spain  Donostia International Physics Center (DIPC),
20018 Donostia-San Sebastián, Spain  Alex Fétida  Institut de Phsyique et
Chime de Matériaux de Strasbourg, CNRS, Strasbourg, Cédex, France  Roberto
Robles  Centro de Física de Materiales CFM/MPC (CSIC-UPV/EHU), 20018 Donostia-
San Sebastián, Spain  Parmenio Boronat  Centro de Física de Materiales CFM/MPC
(CSIC-UPV/EHU), 20018 Donostia-San Sebastián, Spain  Institut de Ciència del
Materials, ICMUV, 46980 Paterna, València, Spain  Divya Jyoti  Centro de
Física de Materiales CFM/MPC (CSIC-UPV/EHU), 20018 Donostia-San Sebastián,
Spain  Donostia International Physics Center (DIPC), 20018 Donostia-San
Sebastián, Spain  Nicolás Lorente  Centro de Física de Materiales CFM/MPC
(CSIC-UPV/EHU), 20018 Donostia-San Sebastián, Spain  Donostia International
Physics Center (DIPC), 20018 Donostia-San Sebastián, Spain  Laurent Limot
Institut de Phsyique et Chime de Matériaux de Strasbourg, CNRS, Strasbourg,
Cédex, France  Deung-Jang Choi<EMAIL_ADDRESS>Centro de Física de Materiales
CFM/MPC (CSIC-UPV/EHU), 20018 Donostia-San Sebastián, Spain  Donostia
International Physics Center (DIPC), 20018 Donostia-San Sebastián, Spain
Ikerbasque, Basque Foundation for Science, 48013 Bilbao, Spain

###### Abstract

We create two molecular devices with superconducting junctions, using
nickelocene molecules, single Fe atoms, and Pb electrodes at low temperature.
We find contrasting behavior based on the coordination of the Fe atom: one
device shows low-bias features in its differential conductance due to the
superposition of multiple Andreev reflections (MAR) and Fe-induced in-gap
states. The other reveals interference between MAR and in-gap states,
showcasing the diversity achievable in atomically engineered devices with
identical components.

Superconductivity is a cornerstone in condensed matter physics and technology
both for its conceptual importance as well as for its potential to drive new
research Tinkham (2004). Recent developments are showing two promising trends
that integrate spintronics and superconductivity. The first direction focuses
on spin-polarized current injection Linder and Robinson (2015), while the
second explores non-reciprocal charge transport for rectification purposes
Ando  _et al._ (2020); Nadeem  _et al._ (2023). And indeed the impact of these
new trends can have important consequences in spin-based information
technology Cai  _et al._ (2023). Recently, it has been possible to probe the
atomic limit of non-reciprocal charge transport and rectification, resulting
in the first atomic-scale superconducting diode Trahms  _et al._ (2023). This
last example underscores the captivating strategy of combining atomic spins
and superconducting junctions for probing the atomic limits of radically new
quantum devices.

Among atomic quantum devices, molecule-based devices offer reproducibility and
functionality. In conjunction with superconductivity, interesting effects
appear. A notable example is the enhancement of spin excitation lifetimes, a
key parameter in quantum technologies. Recent experiments have revealed that
molecular magnets exhibit long lifetimes on a superconductor, attributed to
the gapped quasiparticle spectra characteristic of superconductors Heinrich
_et al._ (2013). The magnetic properties of these molecular junctions can vary
considerably, influenced by factors like molecular magnetic anisotropy and the
nature of the electrode coupling Hatter  _et al._ (2015); Farinacci  _et al._
(2018). Furthermore, the high density of states at the superconducting gap
edge amplifies inelastic conductivity signatures in transport measurements.
This enhancement becomes instrumental in distinguishing molecular spin
processes Kezilebieke  _et al._ (2019); Mier  _et al._ (2021a); Drost  _et
al._ (2023).

Figure 1: Two distinct molecular devices. Shown in (a), atomic scheme of a Fe
atom (red circle) underneath the first atomic layer of Pb (111) addressed by a
Nc-molecule-terminated Pb-coated tip. The Nc molecule contains two C5H5 rings
and a Ni atom (blue). (b) Constant-current image of three of these subsurface
Fe atoms close to an Ar-induced defect Song and Seo (2017) ($\mbox{I}=200$ pA,
$\mbox{V}=10$ mV). The second device (c) is a Pb-terminated tip positioned on
a Nc monolayer where the Fe impurity (red) is located at the molecule/Pb (111)
interface. (d) The constant-current image ($\mbox{I}=30$ pA, $\mbox{V}=-30$
mV) shows Nc dimers Bachellier  _et al._ (2016); Mier  _et al._ (2021a), where
two Fe atoms are located underneath two Nc molecules (brighter spots).

In this study, we demonstrate the effective integration of molecules into
superconducting molecular devices. Our approach involved the construction of
two distinct devices, each comprising a single magnetic molecule, nickelocene
(Nc, formed by two cyclo-pentadyenil rings, C5H5, and an inner nickel atom),
and an iron (Fe) atom, positioned between two lead (Pb) superconducting
electrodes in a scanning tunneling microscope (STM) setup. These junctions
exhibit weak-link characteristics Likharev (1979); Brand  _et al._ (2017) and
achieve stability through the inclusion of the molecule SI .

The weak link in these devices is characterized by distinct tunneling
phenomena: Cooper-pair tunneling manifesting as the Josephson effect at zero
bias and multiple Andreev reflection (MAR) processes at finite bias. Such
phenomena are not unique to our system but have been observed in other STM
superconducting configurations as well Naaman  _et al._ (2001); Naaman and
Dynes (2004); Rodrigo  _et al._ (2004); Ternes  _et al._ (2006); Randeria  _et
al._ (2016); Jäck  _et al._ (2015); Farinacci  _et al._ (2018); Karan  _et
al._ (2022). Notably, in our experiments, identical components were employed
to construct different atomic-scale devices, each distinguished by unique
atomic arrangements. This resulted in devices with markedly different
transport properties.

Figure 2:  Differential conductance of the two different devices. The schemes
(a) and (b) show the two molecular junctions together with the MAR. In (b)
$\varepsilon$ refers to the binding energy of the in-gap state. (c)
Differential conductance of the subsurface Fe atom where a molecule tip is
approached leading to set-point conductances $G=0.04,2.0,20,60$ $\mu$S
(purple, green, cyan, yellow). The purple dI/dV signal ($G=0.04$ $\mu$S) is
multiplied by the marked factor for clarity. (d) Differential conductance
obtained on the Fe below the molecular layer for set-point conductances
$G=10,20,30,40$ $\mu$S (purple, green, cyan, yellow). As in (c), purple
corresponds to the lower transparency and yellow to the higher transparency of
the junction. (e) and (f) are the calculated conductances (color code follows
the same trend as in the experiment), where the dashed lines correspond to
pure MAR processes given by $(\Delta+\Delta_{tip})/n\approx 2\Delta/n=2.7/n$
meV. In (e), we can observe in-gap states induced by the sub-surface Fe atom,
while in (f) we see split-off peaks for odd-numbered MAR showing the
composition of in-gap and MAR in the detected peaks. Due to the truncation in
Floquet terms and the absence of fluctuating phases, the very-low bias region
is not computed.

To conduct our measurements, we used an STM at a base temperature of 2.5 K
under ultra-high-vacuum conditions. We evaporated the Nc molecules onto the Pb
(111) surface resulting in self-assembled molecular layers SI . We used a Pb-
coated STM tip, leading to a gap of $\sim 2\Delta$ in the differential
conductance as a function of bias ($\Delta+\Delta_{tip}\approx 2\Delta$). The
procedure was completed with the in-situ Fe evaporation, reaching a coverage
of less than 1% of individual Fe atoms on the Pb surface.

Utilizing the atomic-manipulation capabilities of our STM, we constructed the
first device depicted in Fig. 1 (a). This was achieved by first picking up one
of the Nc molecules from the surface using the Pb-coated tip and then
approaching it to an Fe atom. The incorporation of the molecule at the tip’s
apex Mier  _et al._ (2021a) allowed for atomic resolution imaging of the Pb
(111) surface, Fig. 1 (b). The close match between the experimental and
computed STM images using density functional theory (DFT) SI  confirms that
the Fe atoms are subsurface, in agreement with previous observations for other
transition-metal atoms on Pb surfaces Trahms  _et al._ (2023); Ruby  _et al._
(2015); Choi  _et al._ (2017).

The other device, Fig. 1 (c), is created by bringing the STM tip close to the
Nc ontop of a Fe atom, seen as bright spots in Fig. 1 (d). The new device
consists of two superconducting electrodes (tip and surface) and a Nc-Fe
junction where the Nc molecule is coupled to the tip and the Fe is coupled to
the surface. Furthermore, both the constant-current image of Fig. 1 (d) and
our DFT calculations show that the lower cyclo-pentadienyl ring increases its
overlap with the Fe atom creating a strong bond that prevents the Fe atom from
slipping beneath the Pb (111) surface SI . Consequently, the tunneling
electrons from the Pb-coated tip are able to interact with the superconducting
surface, which includes Fe-induced in-gap states Randeria  _et al._ (2016)
within the molecular junction.

In Figure 2, panels (a) and (b) provide schematic illustrations of the two
devices. Panel (a) depicts a subsurface Fe atom, denoted by a red dot, beneath
the superconducting substrate. This atom is probed by a Nc-terminated STM tip,
represented as a pair of cyclo-pentadienyl rings in yellow, connected by a Ni
atom, shown as a blue dot. On the right of this illustration, an energy versus
distance diagram presents a one-electron view of the electron-transfer
processes at the tunnel junction, in analogy with established representations
of electron tunneling in superconducting systems Tinkham (2004); Cuevas and
Belzig (2003). In this configuration, the magnetic impurity, being remote from
the surface, does not disrupt the MAR processes. Conversely, panel (b)
illustrates the scenario where an Fe atom is positioned between the molecule
and the substrate. Here, the presence of the magnetic impurity on the surface
modifies the in-gap states, which in turn, interacts with MAR processes as
depicted by the more intricate scheme. These unconventional MAR processes can
be pictured as excitations of the in-gap states by electrons or holes that are
Andreev reflected from the superconductor Andersen  _et al._ (2011); Randeria
_et al._ (2016); Farinacci  _et al._ (2018); Karan  _et al._ (2022); Villas
_et al._ (2020).

The experimental implications of these two distinct configurations are
captured in panels (c) and (d). Each panel exhibits four differential
conductance ($dI/dV$) curves as a function of the sample bias for various STM
set points. These set points range from a tunneling regime (illustrated in
purple), predominantly displaying the quasiparticle peak at $\sim 2\Delta$,
along with minor in-gap features, to higher differential conductance (shown in
green, cyan, and yellow as the conductance increases) where in-gap structures
become more pronounced.

The low-bias electron-transport behaviors of the two devices are starkly
different. Indeed, at large junction resistance (purple in Fig. 2 (c) and
(d)), both curves show a clear quasiparticle peak at $2\Delta$ and some
smaller features near $\Delta$. Increasing the junction conductance by
bringing the tip closer modifies these features differently for the two
devices: for (c), the $2\Delta$ peaks are mostly stable, getting broader; for
(d), however, the peak splits giving a shoulder at $2\Delta$ and a peak
shifting to lower bias. In case (c), the intensity of these peaks invariably
increases with conductance, whereas in case (d), the intensity of the split
peak initially increases but subsequently decreases, exhibiting an unexpected
non-monotonic behavior. The next peak at lower bias exhibits a consistent
growth in magnitude with increased conductance in both cases, a trend that is
also observed for smaller bias. Ultimately, a pronounced Josephson-like
feature is observed, which can be accounted for in terms of the thermal
fluctuations of the superconducting phases Naaman  _et al._ (2001); Naaman and
Dynes (2004); Steiner  _et al._ (2023); SI .

To gain a deeper understanding of the observed phenomena, we have formulated
two models comprising two BCS superconductors linked via a molecular bridge,
here the Nc molecule. The molecular bridge is accounted for by an electron
transmission, denoted as $\alpha$, and the Fe atom is simplified to a Kondo
Hamiltonian. In this notation, $J$ represents the Kondo exchange coupling
constant, while $S$ denotes the spin of the Fe atom, which is treated as a
classical variable. The dynamics of the system is solved using a Floquet
expansion of the two-time non-equilibrium Green’s functions within Nambu space
SI ; Cuevas  _et al._ (1996); Andersen  _et al._ (2011); Taberner (2017). The
difference between the two models is how the Fe atom relates to the Nc
transmission SI . For the first device, the Fe atom affects the transmission
in the bulk, while for the second device, the Fe affects the transmission in
the junction.

The theoretical results, depicted in Fig. 2 panels (e) and (f), show good
agreement with the experimental data if the process of moving the STM tip
closer not only increases the electron transmission coefficient $\alpha$ but
also the exchange interaction $J$. By varying these two parameters within our
model, we are able to understand the peak evolution of the experimental
graphs, Figs 2 (c) and (d).

Figure 2 (e) corresponds to a model of a tunneling junction characterized by
$\alpha$ that is further coupled with a second transmission to a bulk
impurity. The bulk impurity yields in-gap states via $JS$. At low $\alpha$,
the full electron transmission can be seen as the contribution of two coherent
electronic steps: firstly, the transmission to the substrate with MAR
depending on $\alpha$; secondly, the transmission through the substrate
modified by the in-gap states introduced by the bulk impurity SI .

As a consequence, Fig. 2 (e) shows the superposition of the clean-substrate
MAR (dashed lines for the first three MAR orders) together with the in-gap
states. In agreement with Fig. 2 (c), the peak at $2\Delta$ is stationary
because it corresponds to the MAR $n=1$ peak, independent of the subsurface Fe
atom. In contrast, the succeeding peak at a lower voltage originates from an
in-gap state solely produced by the Fe impurity. Experimentally, the peak at
$V\approx 1.15$ mV for lower conductance (purple line) corresponds to a
thermal replica of this Fe-induced in-gap state, easily confirmed by our
theoretical results. As conductance increases, this feature is overshadowed by
the more rapidly intensifying $n=3$ MAR peak. This analysis shows that the
present Nc molecular bridge at 2.5 K allows us to reveal three orders of MAR.

Figure 2 (f) illustrates a model in which MAR processes incorporate in-gap
states, as detailed in Refs. Randeria  _et al._ (2016); Villas  _et al._
(2020); SI . Notably, even at lower tunneling rates, the peak near $2\Delta$
exhibits a discernible shift and distortion. This peak splits into two
distinct components as the tunneling rate increases: one remains anchored at
$2\Delta$, while the other evolves in response to an increase in $JS$. To
reproduce our experimental observations, we increase $JS$ concomitant with the
enhancement in electron transmission $\alpha$. Specifically, the second MAR
peak (n=2) shows a pronounced dependency on $\alpha$, overcoming the in-gap
state split from $2\Delta$. The $n=2$ peak reveals a composite nature,
comprising in-gap state contributions. Similarly, the third MAR peak (n=3)
mirrors the behavior of the first (n=1), forming a sharp shoulder. This leads
us to identify an odd-even effect in the MAR peaks: odd-numbered peaks (n=1,
3, …) exhibit in-gap contributions causing additional broadening and
distortion of the low-bias peaks. In contrast, even-numbered peaks remain
stationary, with more intensity and broadening as the junction’s transparency
increases. The root of this odd-even effect lies in the alignment of the MAR’s
final state with an in-gap state for odd-ordered reflections.

Density Functional Theory (DFT) calculations enable us to elucidate the
differences between the two devices SI . Figure 3 (a) and (b) show the
projected density of states (PDOS) on the molecular states (represented in
green for Nc) and on the Fe atom (in red for Fe), corresponding to the
molecule adsorbed vertically as per the atomic schemes of the two devices in
(c) and (d) respectively.

The PDOS for the first device (a) distinctly reveals the alignment of majority
and minority spins between Fe and Nc, resulting in a ferromagnetic
configuration. This is exemplified in the spin density in (c), where red and
blue isosurfaces represent up-spin and down-spin densities, respectively.
Here, the Fe and Nc molecules interact indirectly through the first Pb layer,
as reflected in their mutual magnetic ordering. Compressing the molecule by 1
Å, a notable spin inversion occurs in the Fe atom (indicated by magenta dashed
lines for Fe (-1)).These results indicate that, despite the close proximity of
Fe and the molecule in the first device, significant structural compression is
required for direct interaction.

Figures 3 (b) and (d) pertain to the interface-Fe device. Here, Fe and the
molecule maintain an antiferromagnetic configuration, even under 1 Å
compression, showing that Fe and Nc are in direct interaction. Consequently,
increasing their mutual interaction via molecular compression does not alter
their mutual magnetic ordering but leads to significant broadening of the
spin-up frontier orbitals of Nc.

Figure 3:  Density-functional theory results. The projected density of states
(PDOS) on the molecular states (green, Nc) and on the (a) subsurface or (b)
interface Fe atom (red, Fe) corresponding to the atomic schemes in (c) and
(d), respectively. The PDOS for the first device (a) clearly shows the
alignment of majority and minority spins between Fe and Nc, in agreement with
the spin density in (c), where the red isosurface shows the up-spin density
and blue the down-spin one. When the top cyclopentadienyl ring of the molecule
is rigidly displaced 1 Å towards the surface (equivalent to reducing the
junction gap in 1 Å), the Fe atom changes its spin (magenta dashed lines, Fe
(-1)). In the second device (b) we see that the majority-spin PDOS have
opposite signs corresponding to the antiferromagnetic configuration shown in
(d).

In conclusion, the integration of Nc molecules, Fe adatoms, and Pb electrodes
yields distinctly different electronic transport behaviors, critically
dependent on the positioning of the Fe adatom relative to the Nc molecule. The
computed diffusion barrier for Fe adsorbed on the surface is about 50 meV,
leading to its subsurface localization, as corroborated by experimental
imagery SI . This results in the Fe atom being encapsulated by the Pb layer,
which attenuates the direct coupling to the Nc tip. Consequently, in this
configuration, the Nc molecule functions as a non-magnetic intermediary in the
superconducting electrodes, while the Fe atoms contribute in-gap states that
are extrinsic to the MAR processes.

In stark contrast, when the Nc molecule is adsorbed onto the surface prior to
Fe atom addition, the Fe atom tends to nestle between the Nc molecule and the
surface, forming a robust bond. This setup gives rise to a significantly
diminished total magnetic moment, with the Nc and Fe elements aligning in an
antiferromagnetic fashion. This pronounced interaction between the Fe and Nc
alters the weak link considerably, with the in-gap states induced by the Fe
atom directly influencing the MAR processes, thereby manifesting a distinctly
different signature from that observed in the prior configuration.

These findings underscore the feasibility of engineering molecular systems
that alter the properties of superconducting devices. Notably, molecules
embodying multiple magnetic centers can substantially expand the design
parameters essential for the advancement of superconducting spintronic
applications.

We thank projects PID2021-127917NB-I00 by MCIN/AEI/10.13039/501100011033,
QUAN-000021-01 by Gipuzkoa Provincial Council, IT-1527-22 by Basque
Government, 202260I187 by CSIC, ESiM project 101046364 by EU, Marie
Skłodowska-Curie grant 847471 by EU, the International Center for Frontier
Research in Chemistry (Strasbourg), and computational resources by Finisterrae
II (CESGA). Views and opinions expressed are however those of the author(s)
only and do not necessarily reflect those of the EU. Neither the EU nor the
granting authority can be held responsible for them.

## References

  * Tinkham (2004) M. Tinkham,  _Introduction to Superconductivity_, 2nd ed. (Dover Publications, 2004).
  * Linder and Robinson (2015) J. Linder and J. W. A. Robinson, Nat Phys 11, 307 (2015).
  * Ando  _et al._ (2020) F. Ando, Y. Miyasaka, T. Li, J. Ishizuka, T. Arakawa, Y. Shiota, T. Moriyama, Y. Yanase, and T. Ono, Nature 584, 373 (2020).
  * Nadeem  _et al._ (2023) M. Nadeem, M. S. Fuhrer, and X. Wang, Nature Reviews Physics 5, 558 (2023).
  * Cai  _et al._ (2023) R. Cai, I. Zutic, and W. Han, Adv. Quantum Tech. 6, 2200080 (2023).
  * Trahms  _et al._ (2023) M. Trahms, L. Melischek, J. F. Steiner, B. Mahendru, I. Tamir, N. Bogdanoff, O. Peters, G. Reecht, C. B. Winkelmann, F. von Oppen, and K. J. Franke, Nature 615, 628 (2023).
  * Heinrich  _et al._ (2013) B. W. Heinrich, L. Braun, and K. J. Pascual, J. I. Franke, Nature Physics 9, 765 (2013).
  * Hatter  _et al._ (2015) N. Hatter, B. W. Heinrich, M. Ruby, J. I. Pascual, and K. J. Franke, Nature Communications 6, 8988 (2015).
  * Farinacci  _et al._ (2018) L. Farinacci, G. Ahmadi, G. Reecht, M. Ruby, N. Bogdanoff, O. Peters, B. W. Heinrich, F. von Oppen, and K. J. Franke, Physical Review Letters 121, 196803 (2018).
  * Kezilebieke  _et al._ (2019) S. Kezilebieke, R. Žitko, M. Dvorak, T. Ojanen, and P. Liljeroth, Nano Letters 19, 4614 (2019).
  * Mier  _et al._ (2021a) C. Mier, B. Verlhac, L. Garnier, R. Robles, L. Limot, N. Lorente, and D.-J. Choi, The Journal of Physical Chemistry Letters 12, 2983 (2021a).
  * Drost  _et al._ (2023) R. Drost, S. Kezilebieke, J. L. Lado, and P. Liljeroth, Phys. Rev. Lett. 131, 086701 (2023).
  * Song and Seo (2017) S. Y. Song and J. Seo, Scientific Reports 7, 12177 (2017).
  * Bachellier  _et al._ (2016) N. Bachellier, M. Ormaza, M. Faraggi, B. Verlhac, M. Vérot, T. Le Bahers, M.-L. Bocquet, and L. Limot, Physical Review B 93, 195403 (2016).
  * Likharev (1979) K. K. Likharev, Rev. Mod. Phys. 51, 101 (1979).
  * Brand  _et al._ (2017) J. Brand, P. Ribeiro, N. Néel, S. Kirchner, and J. Kröger, Phys. Rev. Lett. 118, 107001 (2017).
  * (17) See the supplementary material.
  * Naaman  _et al._ (2001) O. Naaman, W. Teizer, and R. C. Dynes, Phys. Rev. Lett. 87, 097004 (2001).
  * Naaman and Dynes (2004) O. Naaman and R. C. Dynes, Solid State Com. 129, 299 (2004).
  * Rodrigo  _et al._ (2004) J. G. Rodrigo, H. Suderow, S. Vieira, E. Bascones, and F. Guinea, J Phys: Cond Mat 16, R1151 (2004).
  * Ternes  _et al._ (2006) M. Ternes, W.-D. Schneider, J.-C. Cuevas, C. P. Lutz, C. F. Hirjibehedin, and A. J. Heinrich, Phys. Rev. B 74, 132501 (2006).
  * Randeria  _et al._ (2016) M. T. Randeria, B. E. Feldman, I. K. Drozdov, and A. Yazdani, Physical Review B 93, 161115 (2016).
  * Jäck  _et al._ (2015) B. Jäck, M. Eltschka, M. Assig, A. Hardock, M. Etzkorn, C. R. Ast, and K. Kern, App. Phys. Lett. 106, 013109 (2015).
  * Karan  _et al._ (2022) S. Karan, H. Huang, C. Padurariu, B. Kubala, A. Theiler, A. M. Black-Schaffer, G. Morrás, A. L. Yeyati, J. C. Cuevas, J. Ankerhold, K. Kern, and C. R. Ast, Nature Physics 18, 893 (2022).
  * Ruby  _et al._ (2015) M. Ruby, F. Pientka, Y. Peng, F. von Oppen, B. W. Heinrich, and K. J. Franke, Physical Review Letters 115, 087001 (2015).
  * Choi  _et al._ (2017) D.-J. Choi, C. Rubio-Verdú, J. de Bruijckere, M. M. Ugeda, N. Lorente, and J. I. Pascual, Nature Communications 8, 15175 (2017).
  * Cuevas and Belzig (2003) J. C. Cuevas and W. Belzig, Phys. Rev. Lett. 91, 187001 (2003).
  * Andersen  _et al._ (2011) B. M. Andersen, K. Flensberg, V. Koerting, and J. Paaske, Phys. Rev. Lett. 107, 256802 (2011).
  * Villas  _et al._ (2020) A. Villas, R. L. Klees, H. Huang, C. R. Ast, G. Rastelli, W. Belzig, and J. C. Cuevas, Phys. Rev. B 101, 235445 (2020).
  * Steiner  _et al._ (2023) J. F. Steiner, L. Melischek, M. Trahms, K. J. Franke, and F. von Oppen, Phys. Rev. Lett. 130, 177002 (2023).
  * Cuevas  _et al._ (1996) J. C. Cuevas, A. Martín-Rodero, and A. L. Yeyati, Phys. Rev. B 54, 7366 (1996).
  * Taberner (2017) C. O. Taberner, Master’s thesis, U. Copenhagen (2017).
  * Ormaza  _et al._ (2017) M. Ormaza, P. Abufager, B. Verlhac, N. Bachellier, M.-L. Bocquet, N. Lorente, and L. Limot, Nat. Comm. 8, 1974 (2017).
  * Harada  _et al._ (1996) Y. Harada, H. Takayanagi, and A. A. Odintsov, Phys. Rev. B 54, 6608 (1996).
  * Johnson (1928) J. B. Johnson, Phys. Rev. 32, 97 (1928).
  * Nyquist (1928) H. Nyquist, Phys. Rev. 32, 110 (1928).
  * Mier  _et al._ (2021b) C. Mier, D.-J. Choi, and N. Lorente, Phys. Rev. B 104, 245415 (2021b).
  * Mier  _et al._ (2022) C. Mier, D.-J. Choi, and N. Lorente, Phys. Rev. Res. 4, L032010 (2022).
  * González  _et al._ (2020) S. A. González, L. Melischek, O. Peters, K. Flensberg, K. J. Franke, and F. von Oppen, Phys. Rev. B 102, 045413 (2020).
  * Rogovin and Scalapino (1974) D. Rogovin and D. Scalapino, Ann. Phys. 86, 1 (1974).
  * Gerlach  _et al._ (2021) E. Gerlach, F. Jiménez, and R. L. Ayala, “Lecture notes in advanced quantum mechanics,” (2021).
  * Salkola  _et al._ (1997) M. I. Salkola, A. V. Balatsky, and J. R. Schrieffer, Phys. Rev. B 55, 12648 (1997).
  * Shiba (1968) H. Shiba, Progress of Theoretical Physics 40, 435 (1968).
  * Moca  _et al._ (2008) C. Moca, E. Demler, B. Jankó, and G. Zaránd, Phys. Rev. B 77, 174516 (2008).
  * Levy Yeyati  _et al._ (1995) A. Levy Yeyati, A. Martín-Rodero, and F. J. García-Vidal, Phys. Rev. B 51, 3743 (1995).
  * Keldysh (1964) L. V. Keldysh, Zh. Eksp. Teor. Fiz. 47, 1515 (1964).
  * Kantorovich (2020) L. Kantorovich, Phys. Rev. B 101, 165408 (2020).
  * Jauho (2006) A. P. Jauho, “Introduction to the keldysh nonequilibrium green function technique,”  (2006).
  * Hyrkäs  _et al._ (2019) M. J. Hyrkäs, D. Karlsson, and R. van Leeuwen, J. of Phys. A: Math. and Theo. 52, 215303 (2019).
  * Dynes  _et al._ (1978) R. C. Dynes, V. Narayanamurti, and J. P. Garno, Phys. Rev. Lett. 41, 1509 (1978).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, The Journal of Chemical Physics 132, 154104 (2010).
  * Grimme  _et al._ (2011) S. Grimme, S. Ehrlich, and L. Goerigk, Journal of Computational Chemistry 32, 1456 (2011).
  * Tersoff and Hamann (1983) J. Tersoff and D. R. Hamann, Phys. Rev. Lett. 50, 1998 (1983).
  * Bocquet  _et al._ (2009) M.-L. Bocquet, H. Lesnard, S. Monturet, and N. Lorente, in  _Computational Methods in Catalysis and Materials Science_, edited by R. A. v. Santen and P. Sautet (Wiley-VCH Verlag GmbH & Co. KGaA, 2009) pp. 199–219.
  * Lorente and Robles (2019) N. Lorente and R. Robles, “STMpw (v1.0b2). Zenodo. https://doi.org/10.5281/zenodo.3581159,”  (2019).