DeepMaterials LLC, Silver Spring, MD 20906, USA

# Inverse Design of Next-generation Superconductors Using Data-driven Deep
Generative Models

Daniel Wines<EMAIL_ADDRESS>Material Measurement Laboratory, National
Institute of Standards and Technology, Gaithersburg, Maryland 20899, USA  Tian
Xie  Microsoft Research AI4Science, Cambridge, United Kingdom CB1 2FB  Kamal
Choudhary  Material Measurement Laboratory, National Institute of Standards
and Technology, Gaithersburg, Maryland 20899, USA

###### Abstract

Finding new superconductors with a high critical temperature ($T_{c}$) has
been a challenging task due to computational and experimental costs. We
present a diffusion model inspired by the computer vision community to
generate new superconductors with unique structures and chemical compositions.
Specifically, we used a crystal diffusion variational autoencoder (CDVAE)
along with atomistic line graph neural network (ALIGNN) pretrained models and
the Joint Automated Repository for Various Integrated Simulations (JARVIS)
superconducting database of density functional theory (DFT) calculations to
generate new superconductors with a high success rate. We started with a DFT
dataset of $\approx$1000 superconducting materials to train the diffusion
model. We used the model to generate 3000 new structures, which along with
pre-trained ALIGNN screening results in 61 candidates. For the top candidates,
we performed DFT calculations for validation. Such approaches go beyond the
funnel-like materials design approaches and allow for the inverse design of
next-generation materials.

Keywords: generative modeling; superconductivity; inverse design; density
functional theory; high-throughput; materials discovery; machine learning

After superconductivity was discovered in 1911 by Onnes 1, the efforts to
identify novel superconducting materials with high transition temperatures
($T_{c}$) has been an intense area of research in materials science and
condensed matter physics 2, 3. There have been systematic computational
efforts to identify Bardeen–Cooper–Schrieffer (BCS) conventional
superconductors 4, 5 with high-$T_{c}$ prior to costly experimental
investigation 6, 7, 8, where density functional theory-perturbation theory
(DFT-PT) calculations have been performed to obtain the electron-phonon
coupling (EPC) parameters. In addition, various machine learning approaches
have been utilized to accelerate the search for high-$T_{c}$ superconductors
6, 9, 10, 11, 12, 13, 14. However, these typical funnel-like funnel screening-
based approaches are not sufficient for inverse materials design, where
instead of engineering from structure to property, the goal is to engineer
from a target property to crystal structure.

The inverse design of novel materials becomes a difficult problem to solve due
to the fact that there are an infinite number of possible materials with
varying properties that are dependent on chemical composition and crystal
structure. In order for inverse design to be successful, an algorithm that can
successfully create new candidate structures based on a high quality and
diverse dataset of materials data is required. There currently exists abundant
sources of high quality calculated material properties in databases such as
Materials Project 15, Open Quantum Materials Database (OQMD) 16, 17, C2DB 18,
19 and our own JARVIS 20 database, which all contain millions of high-
throughput 21, 20 DFT calculations. The main obstacle in inverse design has
been the method or algorithm used to generate new candidate materials. To
circumvent this obstacle, generative machine learning algorithms can be used
to successfully design new candidate structures. With the popularity of
generative tools such as DALL-E 22, which uses deep learning models to
generate digital images from natural language descriptions (prompts), the
interest in using deep generative models for scientific applications has
immensely increased 23, 24, 25, 26. One of the main challenges with generative
models for periodic materials stems from the creation of representations that
are translationally and rotationally invariant 27, 28, 29, 30, 31, 32, 33, 34.

In a recent work, Xie et al. developed a crystal diffusion variational
autoencoder (CDVAE) model24 (https://github.com/txie-93/cdvae) for periodic
structure generation. The CDVAE consists of a variational autoencoder 35 and a
diffusion model 36, 37, 38 that works directly with the atomic coordinates of
the structures and uses an equivariant graph neural network 39 to ensure
invariance without the need for representations such as graphs or descriptors.
The CDVAE consists of three simultaneously trained networks: 1) the encoder,
which encodes onto the latent space, 2) the property predictor, which samples
the latent space and predicts a structure and composition, and 3) the decoder,
which is a diffusion model that denoises the randomly initialized atom types
into a material that is similar to the training set distribution. More details
of the CDVAE method can be found in Xie et al. 24. A recent work by Lyngby and
Thygesen 25 successfully applied the CDVAE model to discover new, stable 2D
materials and vastly expand the space of 2D materials (on the order of
thousands). Another recent work by Moustafa et al. 26 used the CDVAE model to
discover more than 500 new stable one-dimensional materials. In this work, we
trained a CDVAE model (optimizing $T_{c}$ in the latent space) with DFT
computed data of 1058 superconducting materials from the JARVIS database and
generated thousands of new candidate superconductors. We screened these
candidate structures further by predicting the properties using pretrained
deep learning models for fast predictions. After narrowing down the pool of
potential candidate superconductors, we performed DFT calculations to verify
our predictions and assessed the dynamical and thermodynamic stability of the
newly predicted materials.

We utilized the atomistic line graph neural network (ALIGNN) 40
(https://github.com/usnistgov/alignn) to make deep learning predictions of the
superconducting properties of each material. This pretrained model for
superconducting properties was specifically developed in ref. 6. ALIGNN is
implemented in the deep graph library 41 and PyTorch 42. In the ALIGNN
framework, a structure is represented as a graph, where the elements are nodes
and bonds are edges. Nine input features are assigned to each node in the
graph. These features include first ionization energy, electron affinity,
electronegativity, block, valence electrons, group number, covalent radius,
and atomic volume. The edge features are the bond distances, where the cutoff
for the radial basis function is 8 Å. A 12-nearest neighbor periodic graph
construction is used. The line graph is constructed from the atomistic graph
using bond distances as nodes and bond angles as edge features. Edge-gated
graph convolution is used for updating the nodes and edge features by use of a
propagation function. An edge-gated graph convolution on the bond graph with
an edge-gated convolution on the line graph is what composes one layer. Bond
messages are produced from the line graph convolution that propogate to the
atomistic graph, where the bond features and atom features are further
updated. With regards to predicting superconducting properties, we used a
batch size of 16, 90:5:5 split and training for 300 epochs, where the test set
was not used at all during training. We kept the hyperparameters of the model
the same as the original ALIGNN paper 40.

We utilize the publicly available JARVIS 20 infrastructure for our DFT and
deep learning goals mentioned above. JARVIS (Joint Automated Repository for
Various Integrated Simulations, https://jarvis.nist.gov/) is a collection of
databases and tools to automate materials design using classical force-field,
density functional theory, machine learning calculations and experiments.
JARVIS-DFT is a density functional theory based database of over 75000
materials with several material properties such as formation energy, band gap
with different level of theories 43, solar-cell efficiency 44, topological
spin-orbit coupling spillage 45, 46, 47, elastic tensors 48, dielectric
tensors, piezoelectric tensors, infrared and Raman spectrum 49, electric field
gradients 50, exfoliation energies 51, two-dimensional (2D) magnets 52, and
bulk 6 and 2D superconductors 7, all with stringent DFT-convergence setup 53.
Our ALIGNN and CDVAE models were trained on the 1058 DFT calculations of
superconducting properties presented in ref. 6.

To verify top candidate superconductors, we followed the workflow used to
generate the training data (in ref. 6), where we performed EPC calculations
using non-spin polarized DFT-PT 54, 55 (using the interpolated/Gaussian
broadening method 56) with the Quantum Espresso (QE) software package 57,
PBEsol functional 58, and the GBRV 59 pseudopotentials. The EPC parameter is
derived from spectral function ${\alpha}^{2}F(\omega)$ which is calculated as
follows:

${\alpha}^{2}F(\omega)=\frac{1}{2{\pi}N({\epsilon_{F}})}\sum_{qj}\frac{\gamma_{qj}}{\omega_{qj}}\delta(\omega-\omega_{qj})w(q)$
(1)

where $\omega_{qj}$ is the mode frequency, $N({\epsilon_{F}})$ is the DOS at
the Fermi level ${\epsilon_{F}}$, $\delta$ is the Dirac-delta function, $w(q)$
is the weight of the $q$ point, $\gamma_{qj}$ is the linewidth of a phonon
mode $j$ at wave vector $q$ and is given by:

$\gamma_{qj}=2\pi\omega_{qj}\sum_{nm}\int\frac{d^{3}k}{\Omega_{BZ}}|g_{kn,k+qm}^{j}|^{2}\delta(\epsilon_{kn}-\epsilon_{F})\delta(\epsilon_{k+qm}-\epsilon_{F})$
(2)

Here, the integral is over the first Brillouin zone, $\epsilon_{kn}$ and
$\epsilon_{k+qm}$ are the DFT eigenvalues with wavevector $k$ and $k+q$ within
the $n$th and $m$th bands respectively, $g_{kn,k+qm}^{j}$ is the electron-
phonon matrix element. $\gamma_{qj}$ is related to the mode EPC parameter
$\lambda_{qj}$ by:

$\lambda_{qj}=\frac{\gamma_{qj}}{\pi hN(\epsilon_{F})\omega_{qj}^{2}}$ (3)

Now, the EPC parameter is given by:

$\lambda=2\int\frac{\alpha^{2}F(\omega)}{\omega}d\omega=\sum_{qj}\lambda_{qj}w(q)$
(4)

with $w(q)$ as the weight of a $q$ point. The superconducting transition
temperature, $T_{c}$ can then be approximated using McMillan-Allen-Dynes 60
equation as follows:

$T_{c}=\frac{\omega_{log}}{1.2}\exp\bigg{[}-\frac{1.04(1+\lambda)}{\lambda-\mu^{*}(1+0.62\lambda)}\bigg{]}$
(5)

where

$\omega_{log}=\exp\bigg{[}\frac{\int
d\omega\frac{\alpha^{2}F(\omega)}{\omega}\ln\omega}{\int
d\omega\frac{\alpha^{2}F(\omega)}{\omega}}\bigg{]}$ (6)

In Eq. 5, the parameter $\mu^{*}$ is the effective Coulomb potential
parameter, which we take as 0.1. It is important to note that the robustness
of this workflow was heavily benchmarked against experimental data and higher
levels of theory in ref. 6 and 7, which indicates that the training data for
these deep learning models is of high quality, given the level of theory used
to produce the data.

Figure 1: a) The full inverse design workflow for new superconductors using
DFT, ALIGNN and the CDVAE generative model and b) the probability that
compounds containing a given element in the CDVAE generated structures have an
ALIGNN predicted $T_{c}$ $>$ 3 K.

The full inverse design workflow proposed in this work is depicted in Fig.
1a). The first step of this workflow involved training the CDVAE model on the
1058 DFT calculations in the JARVIS-SC (superconducting) database (from ref.
6). With regards to the CDVAE model for inverse design, the target property to
optimize was $T_{c}$ for new candidate superconductors. After training on the
JARVIS-SC data, we generated 3000 candidate materials using the CDVAE model,
where $T_{c}$ was optimized in the latent space. To ensure that the CDVAE
model optimized $T_{c}$ in the latent space and actually “learned”
superconductivity from the JARVIS-SC training data, we analyzed the property
loss function (which is depicted in Fig. S1). As seen in Fig. S1, the loss
function is sharply decreasing and converging around 1.2 K. This converged
value is comparable to the mean absolute error (MAE) achieved in the ALIGNN
model for $T_{c}$ reported in ref. 6.

With such a large amount of crystal structures to investigate, it is
impracticable to perform DFT calculations of the EPC workflow for all of the
candidates. For this reason, we used our deep learning property prediction
tool (ALIGNN) to screen all 3000 candidate materials. ALIGNN has demonstrated
success for predicting properties such as formation energy and band gap and
more recently superconducting $T_{c}$. The screening criteria we established
for further investigation of these candidates was an ALIGNN predicted: $T_{c}$
$>$ 5 K, $E_{\textrm{form}}$ $<$ 0 eV/atom, and $E_{\textrm{gap}}$ $<$ 0.05
eV. This criterion is based on the fact that we want to further investigate
superconducting materials which have a high $T_{c}$, are potentially stable,
and are metallic (high density of states at the Fermi level). Although
negative formation energy is a stringent requirement for stability,
investigations of dynamical and thermodynamic stability are still needed to
confirm stability (phonon spectrum and energy above the convex hull). After
performing this screening with ALIGNN, we found 61 materials which fit this
criterion. We then went on to perform the full DFT-EPC workflow for these 61
materials and computed the $T_{c}$ using the McMillan-Allen-Dynes equation.
After computing the $T_{c}$ of these materials with DFT, we found that 32
structures have a $T_{c}$ above 5 K. Upon investigation of the phonon density
of states and chemical composition, we find that 7 of these structures have
negative phonon frequencies, indicating dynamical instability. The remaining
25 candidate superconductors are presented in Table 1.

Figure 2: The chemical formula, distribution of densities, and distribution of
packing fractions for a) - c) the 3000 CDVAE generated structures and for d) -
f) the 1058 structures used for training from the JARVIS-DFT database.

Before discussing the more detailed results of the 25 new candidate
superconducting materials generated by our workflow depicted in Fig. 1a), we
will discuss the intermediate steps of the workflow in more detail. The CDVAE
model has the capability to produce new structures with chemical and
structural diversity. Fig. 1b) depicts the probability that compounds
containing a given element in the 3000 CDVAE generated structures have an
ALIGNN predicted $T_{c}$ above 3 K (displayed as a heat map overlayed on the
periodic table). We observe that compounds containing C, N, B, O, V, Mn, Nb,
Ru, and Ta are the most abundant with regards to structures that have a higher
predicted $T_{c}$ with ALIGNN. This is not surprising, due to the fact that
several of the materials in the DFT training set that have a high value of
$T_{c}$ contain such elements (see ref. 6 for details).

Although the abundance of elements in the CDVAE structures (with an ALIGNN
predicted $T_{c}$ $>$ 3 K) follows similar trends to the training data, the
CDVAE candidate structures are not restricted to stoichiometries and crystal
structures already present in the training data. This is illustrated by Fig.
2, where the chemical formula, distribution of densities, and distribution of
packing fractions for a) - c) the 3000 CDVAE generated structures and for d) -
f) the 1058 structures used for training from the JARVIS-DFT database is
depicted. A comparison of Fig. 2a) and 2d) emphasizes the difference in the
stoichiometry between the CDVAE generated structures and the training data
from JARVIS. For example, in the JARVIS training set, a majority of the data
contains only three different atomic species (A, B, and C) and the unit cells
have 5 formula units (see Fig. 2d)). In contrast, the CDVAE generated data
based on the JARVIS training set see Fig. 2a)) has a number of structures that
contain four atomic species (A, B, C, D) and unit cells that have a larger
number of formula units (i.e., AB7). In addition, the CDVAE generated data
contains a substantial number of monoelemental structures, while this is not a
prevalent component in the training set. When comparing Fig. 2b) and 2e) and
2c) and 2f), we observe very different distributions for densities and packing
fractions between the CDVAE structures and JARVIS training data. Specifically,
the CDVAE generated structures have a much larger distribution of structures
with a low density and packing fraction. This can be due in part to the low
symmetry of the CDVAE generated structures, which gives rise to larger crystal
volumes for the newly generated unit cells and therefore smaller densities and
packing fractions. In fact, the CDVAE structures are all of space group P1.
The likelihood of the CDVAE model generating low symmetry, chemically diverse
structures might be attributed to the CDVAE model sampling from a Gaussian
distribution to create new materials (to predict the number of atoms and
composition). Since the underlying distribution of the materials is non-
Gaussian, and the Gaussian distribution that the CDVAE model samples from is
not representative of the latent space, materials out of the distribution can
be generated. This is a limitation of CDVAE which could be addressed with a
better latent space encoding in the future.

Figure 3: The distribution of $T_{c}$ for the a) 1058 JARVIS-DFT structures
($T_{c}$ computed with DFT) and b) the 3000 CDVAE structures ($T_{c}$ computed
with ALIGNN). c) - d) the relation between EPC parameters for the CDVAE
candidate materials verified with DFT.

After discussing the structural and chemical diversity of the CDVAE generated
structures, we plan to discuss the steps of the workflow that involve the
prediction of $T_{c}$ with deep learning (ALIGNN) and DFT. The ALIGNN deep
learning property predictor for $T_{c}$ was previously pretrained on the
dataset in ref. 6 (1058 materials) and benchmarked. Using the pretrained
ALIGNN model allows us to filter the vast amount of CDVAE candidate
superconductors with instantaneous property prediction. Fig. 3a) and b) depict
the distribution of $T_{c}$ for the a) 1058 JARVIS-DFT structures ($T_{c}$
computed with DFT) and b) the 3000 CDVAE structures ($T_{c}$ computed with
ALIGNN). As seen in the figure, the distribution of $T_{c}$ is vastly
different. Most strikingly, the CDVAE generated data has a very low amount of
structures with a predicted $T_{c}$ close to zero (non-superconducting), as
opposed to a large number of materials in the JARVIS training set that have a
$T_{c}$ close to zero. In addition, the distribution of CDVAE+ALIGNN
temperatures resembles more of a Gaussian distribution, with a majority of the
$T_{c}$ values of the newly generated structures concentrated around 2.5 K - 6
K. This can be an indicator of the success of the CDVAE model for generating
new structures with an optimized target property, such as $T_{c}$, in the
latent space.

After screening the 3000 CDVAE generated structures with ALIGNN, we performed
DFT for the 61 candidates that met the ALIGNN screening criteria depicted in
Fig. 1a). At this stage, it is important to acknowledge the success of the
pretrained ALIGNN model with regards to predicting the $T_{c}$ of these newly
generated structures. Out of the 61 structures, 54 of them have a DFT computed
$T_{c}$ $>$ 0.5 K, which indicates that ALIGNN is 89 $\%$ successful with
filtering out superconducting materials from the CDVAE generated structures.
Out of the 61 structures, 32 of them have a DFT computed $T_{c}$ $>$ 5 K,
which indicates that ALIGNN has a 52 $\%$ success rate in filtering out
materials with a $T_{c}$ above 5 K from the CDVAE generated structures. More
details of the ALIGNN screening are given in the SI. It is important to note
that the ALIGNN predictions used for prescreening were performed on structures
generated from CDVAE prior to geometric relaxation with DFT. A full
description of the DFT and ALIGNN (on the relaxed and unrelaxed structures)
predictions are given in Table S1. To demonstrate how the predictive power of
ALIGNN extrapolates from unrelaxed structures to relaxed structures for
$T_{c}$ predictions, we plot the ALIGNN $T_{c}$ from the unrelaxed structures
versus the ALIGNN $T_{c}$ from the relaxed structures in Fig. S2, where we
observe that the predictions are relatively consistent between structures. We
perform similar analysis for formation energy and band gap in Fig. S3 and
Table S2. For $E_{\textrm{form}}$, we observe much closer agreement when
comparing DFT and ALIGNN predictions for the relaxed (MAE of 0.15 eV/atom)
structures compared to ALIGNN predictions for the unrelaxed (MAE of 0.42
eV/atom) structures (see Fig. S3). With regards to band gap screening, for 41
of the materials we ran DFT for, DFT and ALIGNN predicted them to be metallic
(in relaxed and unrelaxed configurations). For the remaining 20 structures, at
least one of the three methods predicted a band gap (either DFT, ALIGNN on
unrelaxed structures or ALIGNN on relaxed structures). A comparison of these
band gap results is given in Table S2. The success of the ALIGNN prediction of
$T_{c}$ is a direct consequence of the availability of training data. In terms
of deep learning models for material property prediction, 1058 structures is a
relatively small amount of data to train on. The reason for the smaller
training set size is due to the vast computational expense needed to perform
these DFT calculations of EPC properties. In principle, the ALIGNN prediction
of $T_{c}$ can be systematically improved by adding more DFT calculations to
the training, which is an ongoing effort of JARVIS. We show the relationship
between EPC parameters ($\lambda$, $\omega_{\textrm{log}}$, $T_{c}$) for the
61 materials we performed DFT calculations for in Fig. 3c) and d). Fig. 3c)
depicts an inverse relationship between $\lambda$ and $\omega_{\textrm{log}}$
and in Fig. 3d), we observe a somewhat positive relationship between $\lambda$
and $T_{c}$. These are typical behaviors of BCS superconductors and were
observed in our work on BCS bulk and 2D superconductors 6, 7. From the
colormap of Fig. 3c) and d), it is clear that a balance of high $\lambda$ and
$\omega_{\textrm{log}}$ is a necessary condition for a material to have a high
$T_{c}$. It is important to note that our focus on ambient condition
stoichiometric BCS superconductors significantly limits the value of $T_{c}$
that can be achieved. High-pressure superconductors, non-stoichiometric (i.e.
several cuprates and pnictides) superconductors, and unconventional
superconductors (not mediated by electron-phonon interactions) can achieve
much higher $T_{c}$ than BCS superconductors. We hope to extend our workflow
to these types of superconductors in the future, but for proof-of-concept and
validation of our methods, we focus on BCS superconductors in this work.

Figure 4: Top and side view of the top superconductor candidates (closest to the convex hull) generated with CDVAE and verified with DFT: a) VN2, b) MnN, c) VN, d) BN2Zr, e) NP2Sr, f) TaP2, g) AlN2Zr, h) NbRh, i) AlN2V (I), j) TiO2NbN, k) NBaP, l) ScO3Zr, m) Al2N, n) AlN2V (II), o) TiN. $T_{c}$, $E_{\textrm{form}}$ and $E_{\textrm{hull}}$ are also given for each material.  Table 1:  Chemical formula, JARVIS ID (JID), $T_{c}$, formation energy per atom and energy above the convex hull per atom of the 25candidate superconductors from CDVAE verified by DFT calculations. Structure | JID | $T_{c}$ | $E_{\textrm{form}}$ | $E_{\textrm{hull}}$
---|---|---|---|---
|  | (K) | (eV/atom) | (eV/atom)
VN2 | JVASP-161655 | 20.2 | -0.56 | 0.32
NTaB (I) | JVASP-161630 | 16.2 | -0.32 | 0.84
BORu | JVASP-161610 | 12.4 | -0.72 | 0.75
BTa2N | JVASP-161612 | 11.9 | -0.37 | 0.67
MnN | JVASP-161621 | 11.4 | -0.43 | 0.01
NTaB (II) | JVASP-161624 | 11.1 | -0.32 | 0.84
VN | JVASP-161656 | 9.6 | -0.98 | 0.34
BN2Zr | JVASP-161608 | 9.3 | -1.23 | 0.43
BTaNS | JVASP-161613 | 9.1 | -0.33 | 0.88
NP2Sr | JVASP-162663 | 8.4 | -0.54 | 0.22
TaP2 | JVASP-161649 | 8.1 | -0.45 | 0.25
AlN2Zr | JVASP-161600 | 7.8 | -1.40 | 0.39
NPdTi2 | JVASP-161629 | 7.7 | -0.78 | 0.47
PScSi2 | JVASP-161644 | 7.4 | -0.43 | 0.49
NbRh | JVASP-161638 | 7.4 | -0.41 | 0.00
AlN2V (I) | JVASP-161599 | 7.2 | -1.14 | 0.34
TiO2NbN | JVASP-161653 | 7.1 | -2.15 | 0.30
NBaP | JVASP-161626 | 6.6 | -0.59 | 0.22
NVBRu | JVASP-161631 | 6.3 | -0.14 | 0.74
ScO3Zr | JVASP-161647 | 6.3 | -3.43 | 0.15
Al2N | JVASP-161597 | 5.9 | -0.79 | 0.31
AlN2V (II) | JVASP-162662 | 5.4 | -1.14 | 0.34
TiN | JVASP-161652 | 5.3 | -1.56 | 0.35
ScBORuCa | JVASP-161646 | 5.2 | -1.12 | 0.58
B2TaS | JVASP-161604 | 5.2 | -0.07 | 0.58

The chemical composition, $T_{c}$, formation energy per atom and energy above
the convex hull of the top superconducting candidates are given in Table 1. As
seen in the table, all of the 25 structures have negative formation energy.
Although negative formation energy is a necessary prerequisite for
thermodynamic stability, it does not guarantee thermodynamic stability. For
this reason, we calculated the energy above the convex hull (shown in Table
1). We computed the convex hull of these structures from the formation energy
calculations of the different phases in JARVIS-DFT. From here, we observed 15
structures within an energy of 0.45 eV or less from the convex hull, and deem
that these structures are the most likely to be experimentally synthesized.
The atomic structures (all of P1 symmetry) of these 15 materials are given in
Fig. 4 (VN2, MnN, VN, BN2Zr, NP2Sr, TaP2, AlN2Zr, NbRh, AlN2V (I), TiO2NbN,
NBaP, ScO3Zr, Al2N, AlN2V (II), TiN). The properties and atomic structures of
all the materials we performed DFT verification for (61 materials) will be
available along with this manuscript and uploaded to the JARVIS-DFT database
(https://jarvis.nist.gov/). Although these 15 materials have an energy above
the convex hull between 0.0 eV/atom and 0.43 eV/atom, indicating that
thermodynamic stability is still not guaranteed, they still may be
synthesizable. In fact, Aykol et al. 61 demonstrated that nitride-based
polymorphs with a higher energy above the convex hull can be stable up to a
high energetic amorphous limit. Since a majority of these superconductors
contain nitrogen, this gives promise that they can eventually be
experimentally realized.

As a success metric of the CDVAE method for inverse design, it is important to
measure how diverse the chemical composition and crystal structure of the
generated materials are, not only with respect to the training data, but with
respect to larger areas of phase space. Out of the top 25 candidate
superconductors, 20 do not have a chemical composition in the JARVIS training
set used for CDVAE. Compounds with a stoichiometry of MnN, VN, AlN2Zr, NbRh
and TiN are part of the training set and were all found to have
superconducting phases in ref. 6. Although these five candidates have similar
chemistry to compounds in the training set, they all have entirely different
crystal structures, with the exception of NbRh (which has the same structure
as JVASP-20529). Due to the fact that the JARVIS database has less materials
(the primary focus of JARVIS is to expand the properties and accuracy of
materials), we checked other well-established materials databases for these
candidates. Specifically, we checked the Materials Project (over 150000
structures) and OQMD (over 1000000 structures). In addition, we searched the
literature to see if any of these materials have been synthesized. For VN2, we
found 3 compounds in OQMD with entirely different structures (1239835,
1233906, 1589664) with $E_{\textrm{hull}}$ ranging from 0.8 to 2.6 eV/atom.
For Al2N, we found 4 compounds in OQMD with entirely different structures
(1237731, 1415673, 1590028) with $E_{\textrm{hull}}$ ranging from (0.9 to 2.6)
eV/atom. We found 3 TaP2 compounds (1372133, 1590238, 1237725) in OQMD. The
$E_{\textrm{hull}}$ ranged from 0 to 1.2 eV/atom, but none of these phases
resemble the structure in this study, which has a similar crystal structure to
layered transition metal dichalcogenides (such as 2H-MoS2). Upon searching the
literature, we found that a semimetallic phase of TaP2 had been recently
synthesized 62 that possessed a different crystal structure than the
superconducting phase reported in this work. This experimentally synthesized
TaP2 was found to possess quantum oscillations and nontrivial topological
properties, but was not found to be superconducting 62. We found 1 AlN2V
compound in the Materials Project (mp-1247742), where the $E_{\textrm{hull}}$
was 0.37 eV/atom. When comparing this structure to the two AlN2V structures in
our study, mp-1247742 has a similar form to AlN2V (I) but a drastically
different structure than AlN2V (II). With regards to the five structures that
have the same stoichiometry as structures in the training set, MnN, VN, NbRh
and TiN have multiple entries in OQMD and Materials Project, while AlN2Zr has
just one entry in Materials Project. Our NbRh structure is the same as OQMD
30785, 17746, 339226, mp-1963 and JVASP-20529 (our $T_{c}$ is in agreement
with the value obtained in ref. 6) and our TiN structure is similar to the
OQMD 1230416 entry. Upon searching the experimental literature, certain phases
of TiN and VN have been reported to have superconducting properties 63, 64,
65, 66. We additionally searched the Supercon database (the Supercon database
is unavailable as of December 2021, but previous papers have made this dataset
available for further research 67, 68), which consists of over 16000
experimentally realized superconducting materials. We find no common materials
between the top 25 CDVAE candidates and the Supercon database. This further
proves that the CDVAE method can generate unique materials with a specific
desired property, covering previously undiscovered areas of phase space.

In this work, we used a multi step workflow, combining generative models, deep
learning property prediction and DFT to discover next-generation
superconducting materials. We demonstrated that the deep learning property
prediction using the ALIGNN model can accelerate the search for new
superconductors by instantaneously screening the properties of the newly
generated materials prior to DFT verification and experimental investigation.
Our search revealed 25 newly predicted candidate superconductors, with 15
structures being close to the convex hull and $T_{c}$ values as high as 20.2
K. Our approach goes beyond the standard funnel-like materials design workflow
and allows for inverse design of novel materials, populating previously
undiscovered areas of phase space.

##  1 Data Availability Statement

The data from the present work is available at
https://figshare.com/articles/dataset/Inverse_Design_of_Next-
Generation_Superconductors_Using_Data-Driven_Deep_Generative_Models/23681025.

##  2 Code Availability Statement

Software packages mentioned in the article can be found at
https://github.com/usnistgov/jarvis, https://github.com/knc6/cdvae and
https://github.com/txie-93/cdvae.

##  3 Competing interests

The authors declare no competing interests.

##  4 Supporting Information

CDVAE loss function, tabulated data for relaxed and unrelaxed structures
computed with ALIGNN and DFT, additional comparisons of ALIGNN for relaxed and
unrelaxed structures.

##  5 Acknowledgments

All authors thank the National Institute of Standards and Technology for
funding, computational, and data-management resources, specifically the NIST
Nisaba and Raritan HPC clusters. K.C. thanks the computational support from
XSEDE (Extreme Science and Engineering Discovery Environment) computational
resources under allocation number TG-DMR 190095. Contributions from K.C. were
supported by the financial assistance award 70NANB19H117 from the U.S.
Department of Commerce, National Institute of Standards and Technology.

## References

  * Kamerlingh Onnes 1911 Kamerlingh Onnes, H. The resistance of pure mercury at helium temperatures. _Commun. Phys. Lab. Univ. Leiden, b_ 1911, _120_
  * Poole et al. 2013 Poole, C. P.; Farach, H. A.; Creswick, R. J. _Superconductivity_ ; Academic press, 2013
  * Rogalla and Kes 2011 Rogalla, H.; Kes, P. H. _100 years of superconductivity_ ; Taylor & Francis, 2011
  * Cooper and Feldman 2010 Cooper, L. N.; Feldman, D. _BCS: 50 years_ ; World scientific, 2010
  * Giustino 2017 Giustino, F. Electron-phonon interactions from first principles. _Reviews of Modern Physics_ 2017, _89_ , 015003
  * Choudhary and Garrity 2022 Choudhary, K.; Garrity, K. Designing high-TC superconductors with BCS-inspired screening, density functional theory, and deep-learning. _npj Computational Materials_ 2022, _8_ , 244
  * Wines et al. 2023 Wines, D.; Choudhary, K.; Biacchi, A. J.; Garrity, K. F.; Tavazza, F. High-throughput DFT-based discovery of next generation two-dimensional (2D) superconductors. _Nano Letters_ 2023, _23_ , 969–978
  * Shipley et al. 2021 Shipley, A. M.; Hutcheon, M. J.; Needs, R. J.; Pickard, C. J. High-throughput discovery of high-temperature conventional superconductors. _Phys. Rev. B_ 2021, _104_ , 054501
  * García-Nieto et al. 2021 García-Nieto, P. J.; García-Gonzalo, E.; Paredes-Sánchez, J. Prediction of the critical temperature of a superconductor by using the WOA/MARS, Ridge, Lasso and Elastic-net machine learning techniques. _Neural Computing and Applications_ 2021, _33_ , 17131–17145
  * Stanev et al. 2018 Stanev, V.; Oses, C.; Kusne, A. G.; Rodriguez, E.; Paglione, J.; Curtarolo, S.; Takeuchi, I. Machine learning modeling of superconducting critical temperature. _npj Computational Materials_ 2018, _4_ , 29
  * Zhang et al. 2022 Zhang, J.; Zhu, Z.; Xiang, X. D.; Zhang, K.; Huang, S.; Zhong, C.; Qiu, H.-J.; Hu, K.; Lin, X. Machine learning prediction of superconducting critical temperature through the structural descriptor. _The Journal of Physical Chemistry C_ 2022, _126_ , 8922–8927
  * Meredig et al. 2018 Meredig, B.; Antono, E.; Church, C.; Hutchinson, M.; Ling, J.; Paradiso, S.; Blaiszik, B.; Foster, I.; Gibbons, B.; Hattrick-Simpers, J.; Mehta, A.; Ward, L. Can machine learning identify the next high-temperature superconductor? Examining extrapolation performance for materials discovery. _Mol. Syst. Des. Eng._ 2018, _3_ , 819–825
  * Roter and Dordevic 2020 Roter, B.; Dordevic, S. Predicting new superconductors and their critical temperatures using machine learning. _Physica C: Superconductivity and its Applications_ 2020, _575_ , 1353689
  * Menon and Ranganathan 2022 Menon, D.; Ranganathan, R. A generative approach to materials discovery, design, and optimization. _ACS Omega_ 2022, _7_ , 25958–25973
  * Jain et al. 2013 Jain, A.; Ong, S. P.; Hautier, G.; Chen, W.; Richards, W. D.; Dacek, S.; Cholia, S.; Gunter, D.; Skinner, D.; Ceder, G.; Persson, K. A. Commentary: The Materials Project: A materials genome approach to accelerating materials innovation. _APL Materials_ 2013, _1_ , 011002
  * Saal et al. 2013 Saal, J. E.; Kirklin, S.; Aykol, M.; Meredig, B.; Wolverton, C. Materials design and discovery with high-throughput density functional theory: The Open Quantum Materials Database (OQMD). _JOM_ 2013, _65_ , 1501–1509
  * Kirklin et al. 2015 Kirklin, S.; Saal, J. E.; Meredig, B.; Thompson, A.; Doak, J. W.; Aykol, M.; Rühl, S.; Wolverton, C. The Open Quantum Materials Database (OQMD): assessing the accuracy of DFT formation energies. _npj Computational Materials_ 2015, _1_ , 15010
  * Haastrup et al. 2018 Haastrup, S.; Strange, M.; Pandey, M.; Deilmann, T.; Schmidt, P. S.; Hinsche, N. F.; Gjerding, M. N.; Torelli, D.; Larsen, P. M.; Riis-Jensen, A. C.; Gath, J.; Jacobsen, K. W.; Mortensen, J. J.; Olsen, T.; Thygesen, K. S. The Computational 2D Materials Database: high-throughput modeling and discovery of atomically thin crystals. _2D Materials_ 2018, _5_ , 042002
  * Gjerding et al. 2021 Gjerding, M. N. et al. Recent progress of the Computational 2D Materials Database (C2DB). _2D Materials_ 2021, _8_ , 044002
  * Choudhary et al. 2020 Choudhary, K.; Garrity, K. F.; Reid, A. C.; DeCost, B.; Biacchi, A. J.; Hight Walker, A. R.; Trautt, Z.; Hattrick-Simpers, J.; Kusne, A. G.; Centrone, A., et al. The joint automated repository for various integrated simulations (JARVIS) for data-driven materials design. _npj Computational Materials_ 2020, _6_ , 1–13
  * Ong et al. 2013 Ong, S. P.; Richards, W. D.; Jain, A.; Hautier, G.; Kocher, M.; Cholia, S.; Gunter, D.; Chevrier, V. L.; Persson, K. A.; Ceder, G. Python Materials Genomics (pymatgen): A robust, open-source python library for materials analysis. _Computational Materials Science_ 2013, _68_ , 314–319
  * Ramesh et al. 2021 Ramesh, A.; Pavlov, M.; Goh, G.; Gray, S.; Voss, C.; Radford, A.; Chen, M.; Sutskever, I. Zero-shot text-to-image generation. 2021
  * Wu et al. 2022 Wu, K. E.; Yang, K. K.; van den Berg, R.; Zou, J. Y.; Lu, A. X.; Amini, A. P. Protein structure generation via folding diffusion. 2022
  * Xie et al. 2021 Xie, T.; Fu, X.; Ganea, O.-E.; Barzilay, R.; Jaakkola, T. Crystal diffusion variational autoencoder for periodic material generation. 2021; https://arxiv.org/abs/2110.06197
  * Lyngby and Thygesen 2022 Lyngby, P.; Thygesen, K. S. Data-driven discovery of 2D materials by deep generative models. _npj Computational Materials_ 2022, _8_ , 232
  * Moustafa et al. 2023 Moustafa, H.; Lyngby, P. M.; Mortensen, J. J.; Thygesen, K. S.; Jacobsen, K. W. Hundreds of new, stable, one-dimensional materials from a generative machine learning model. _Phys. Rev. Mater._ 2023, _7_ , 014007
  * Fung et al. 2022 Fung, V.; Jia, S.; Zhang, J.; Bi, S.; Yin, J.; Ganesh, P. Atomic structure generation from reconstructing structural fingerprints. _Machine Learning: Science and Technology_ 2022, _3_ , 045018
  * Noh et al. 2020 Noh, J.; Gu, G. H.; Kim, S.; Jung, Y. Machine-enabled inverse design of inorganic solid materials: promises and challenges. _Chem. Sci._ 2020, _11_ , 4871–4881
  * Kim et al. 2020 Kim, S.; Noh, J.; Gu, G. H.; Aspuru-Guzik, A.; Jung, Y. Generative adversarial networks for crystal structure prediction. _ACS Central Science_ 2020, _6_ , 1412–1420
  * Long et al. 2021 Long, T.; Fortunato, N. M.; Opahle, I.; Zhang, Y.; Samathrakis, I.; Shen, C.; Gutfleisch, O.; Zhang, H. Constrained crystals deep convolutional generative adversarial network for the inverse design of crystal structures. _npj Computational Materials_ 2021, _7_ , 66
  * Song et al. 2021 Song, Y.; Siriwardane, E. M. D.; Zhao, Y.; Hu, J. Computational discovery of new 2D materials using deep learning generative models. _ACS Applied Materials & Interfaces_ 2021, _13_ , 53303–53313
  * Noh et al. 2019 Noh, J.; Kim, J.; Stein, H. S.; Sanchez-Lengeling, B.; Gregoire, J. M.; Aspuru-Guzik, A.; Jung, Y. Inverse design of solid-state materials via a continuous representation. _Matter_ 2019, _1_ , 1370–1384
  * Zhao et al. 2021 Zhao, Y.; Al-Fahdi, M.; Hu, M.; Siriwardane, E. M. D.; Song, Y.; Nasiri, A.; Hu, J. High-throughput discovery of novel cubic crystal materials using deep generative neural networks. _Advanced Science_ 2021, _8_ , 2100566
  * Ren et al. 2022 Ren, Z. et al. An invertible crystallographic representation for general inverse design of inorganic crystals with targeted properties. _Matter_ 2022, _5_ , 314–335
  * Kingma and Welling 2013 Kingma, D. P.; Welling, M. Auto-encoding variational bayes. 2013; https://arxiv.org/abs/1312.6114
  * Sohl-Dickstein et al. 2015 Sohl-Dickstein, J.; Weiss, E.; Maheswaranathan, N.; Ganguli, S. Deep unsupervised learning using nonequilibrium thermodynamics. Proceedings of the 32nd International Conference on Machine Learning. Lille, France, 2015; pp 2256–2265
  * Song and Ermon 2019 Song, Y.; Ermon, S. Generative modeling by estimating gradients of the data distribution. 2019; https://arxiv.org/abs/1907.05600
  * Sohl-Dickstein et al. 2015 Sohl-Dickstein, J.; Weiss, E. A.; Maheswaranathan, N.; Ganguli, S. Deep unsupervised learning using nonequilibrium thermodynamics. 2015
  * Batzner et al. 2022 Batzner, S.; Musaelian, A.; Sun, L.; Geiger, M.; Mailoa, J. P.; Kornbluth, M.; Molinari, N.; Smidt, T. E.; Kozinsky, B. E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. _Nature Communications_ 2022, _13_ , 2453
  * Choudhary and DeCost 2021 Choudhary, K.; DeCost, B. Atomistic Line Graph Neural Network for improved materials property predictions. _npj Computational Materials_ 2021, _7_ , 1–8
  * Wang et al. 2019 Wang, M.; Yu, L.; Zheng, D.; Gan, Q.; Gai, Y.; Ye, Z.; Li, M.; Zhou, J.; Huang, Q.; Ma, C., et al. Deep Graph Library: Towards efficient and scalable deep learning on graphs. _arXiv e-prints_ 2019,
  * Paszke et al. 2019 Paszke, A.; Gross, S.; Massa, F.; Lerer, A.; Bradbury, J.; Chanan, G.; Killeen, T.; Lin, Z.; Gimelshein, N.; Antiga, L., et al. Pytorch: An imperative style, high-performance deep learning library. _Advances in neural information processing systems_ 2019, _32_ , 8026–8037
  * Choudhary et al. 2018 Choudhary, K.; Zhang, Q.; Reid, A. C.; Chowdhury, S.; Van Nguyen, N.; Trautt, Z.; Newrock, M. W.; Congo, F. Y.; Tavazza, F. Computational screening of high-performance optoelectronic materials using OptB88vdW and TB-mBJ formalisms. _Scientific data_ 2018, _5_ , 1–12
  * Choudhary et al. 2019 Choudhary, K.; Bercx, M.; Jiang, J.; Pachter, R.; Lamoen, D.; Tavazza, F. Accelerated discovery of efficient solar cell materials using quantum and machine-learning methods. _Chemistry of Materials_ 2019, _31_ , 5900–5908
  * Choudhary et al. 2021 Choudhary, K.; Garrity, K. F.; Ghimire, N. J.; Anand, N.; Tavazza, F. High-throughput search for magnetic topological materials using spin-orbit spillage, machine learning, and experiments. _Physical Review B_ 2021, _103_ , 155131
  * Choudhary et al. 2019 Choudhary, K.; Garrity, K. F.; Tavazza, F. High-throughput discovery of topologically non-trivial materials using spin-orbit spillage. _Scientific reports_ 2019, _9_ , 1–8
  * Choudhary et al. 2020 Choudhary, K.; Garrity, K. F.; Jiang, J.; Pachter, R.; Tavazza, F. Computational search for magnetic and non-magnetic 2D topological materials using unified spin–orbit spillage screening. _NPJ Computational Materials_ 2020, _6_ , 1–8
  * Choudhary et al. 2018 Choudhary, K.; Cheon, G.; Reed, E.; Tavazza, F. Elastic properties of bulk and low-dimensional materials using van der Waals density functional. _Physical Review B_ 2018, _98_ , 014107
  * Choudhary et al. 2020 Choudhary, K.; Garrity, K. F.; Sharma, V.; Biacchi, A. J.; Hight Walker, A. R.; Tavazza, F. High-throughput density functional perturbation theory and machine learning predictions of infrared, piezoelectric, and dielectric responses. _NPJ Computational Materials_ 2020, _6_ , 1–13
  * Choudhary et al. 2020 Choudhary, K.; Ansari, J. N.; Mazin, I. I.; Sauer, K. L. Density functional theory-based electric field gradient database. _Scientific Data_ 2020, _7_ , 1–10
  * Choudhary et al. 2017 Choudhary, K.; Kalish, I.; Beams, R.; Tavazza, F. High-throughput identification and characterization of two-dimensional materials using density functional theory. _Scientific reports_ 2017, _7_ , 1–16
  * Wines et al. 2023 Wines, D.; Choudhary, K.; Tavazza, F. Systematic DFT+U and Quantum Monte Carlo benchmark of magnetic two-dimensional (2D) CrX3 (X = I, Br, Cl, F). _The Journal of Physical Chemistry C_ 2023, _127_ , 1176–1188
  * Choudhary and Tavazza 2019 Choudhary, K.; Tavazza, F. Convergence and machine learning predictions of Monkhorst-Pack k-points and plane-wave cut-off in high-throughput DFT calculations. _Computational materials science_ 2019, _161_ , 300–308
  * Baroni et al. 1987 Baroni, S.; Giannozzi, P.; Testa, A. Green’s-function approach to linear response in solids. _Physical review letters_ 1987, _58_ , 1861
  * Gonze 1995 Gonze, X. Perturbation expansion of variational principles at arbitrary order. _Physical Review A_ 1995, _52_ , 1086
  * Wierzbowska et al. 2005 Wierzbowska, M.; de Gironcoli, S.; Giannozzi, P. Origins of low-and high-pressure discontinuities of $T\\_${$c$} in niobium. _arXiv preprint cond-mat/0504077_ 2005,
  * Giannozzi et al. 2009 Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I., et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. _Journal of physics: Condensed matter_ 2009, _21_ , 395502
  * Perdew et al. 2008 Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Constantin, L. A.; Zhou, X.; Burke, K. Restoring the density-gradient expansion for exchange in solids and surfaces. _Physical review letters_ 2008, _100_ , 136406
  * Garrity et al. 2014 Garrity, K. F.; Bennett, J. W.; Rabe, K. M.; Vanderbilt, D. Pseudopotentials for high-throughput DFT calculations. _Computational Materials Science_ 2014, _81_ , 446–452
  * McMillan 1968 McMillan, W. Transition temperature of strong-coupled superconductors. _Physical Review_ 1968, _167_ , 331
  * Aykol et al. 2018 Aykol, M.; Dwaraknath, S. S.; Sun, W.; Persson, K. A. Thermodynamic limit for synthesis of metastable inorganic materials. _Science Advances_ 2018, _4_ , eaaq0148
  * Wang et al. 2019 Wang, H.; Su, H.; Zhang, J.; Xia, W.; Lin, Y.; Liu, X.; Hou, X.; Yu, Z.; Yu, N.; Wang, X.; Zou, Z.; Wang, Y.; Liang, Q.; Zhen, Y.; Guo, Y. Quantum oscillations and nontrivial topological state in a compensated semimetal ${\mathrm{TaP}}_{2}$. _Phys. Rev. B_ 2019, _100_ , 115127
  * Pan et al. 1989 Pan, V.; Prokhorov, V.; Komashko, V.; Kaminsky, G.; Kousenetsov, M.; Tretiatchenko, C. Superconducting properties of TaN and VN films. _IEEE Transactions on Magnetics_ 1989, _25_ , 2000–2003
  * Romanov et al. 2018 Romanov, N. R.; Zolotov, P. I.; Vakhtomin, Y. B.; Divochiy, A. V.; Smirnov, K. V. Electron diffusivity measurements of VN superconducting single-photon detectors. _Journal of Physics: Conference Series_ 2018, _1124_ , 051032
  * Saveskul et al. 2019 Saveskul, N.; Titova, N.; Baeva, E.; Semenov, A.; Lubenchenko, A.; Saha, S.; Reddy, H.; Bogdanov, S.; Marinero, E.; Shalaev, V.; Boltasseva, A.; Khrapai, V.; Kardakova, A.; Goltsman, G. Superconductivity behavior in epitaxial $\mathrm{Ti}\mathrm{N}$ films points to surface magnetic disorder. _Phys. Rev. Appl._ 2019, _12_ , 054001
  * Jaim et al. 2015 Jaim, H. M. I.; Aguilar, J. A.; Sarabi, B.; Rosen, Y. J.; Ramanayaka, A. N.; Lock, E. H.; Richardson, C. J. K.; Osborn, K. D. Superconducting TiN films sputtered over a large range of substrate DC bias. _IEEE Transactions on Applied Superconductivity_ 2015, _25_ , 1–5
  * vstanev1 2021 vstanev1, Supercon. 2021; https://github.com/vstanev1/Supercon
  * Hamidieh 2018 Hamidieh, K. A data-driven statistical model for predicting the critical temperature of a superconductor. 2018


# Orbital relaxation length from first-principles scattering calculations

Max Rang  Paul J. Kelly  Faculty of Science and Technology and MESA+ Institute
for Nanotechnology, University of Twente, P.O. Box 217, 7500 AE Enschede, The
Netherlands

###### Abstract

The orbital Hall effect generates a current of orbital angular momentum
perpendicular to a charge current. Experiments suggest that this orbital
current decays on a long length scale that is of the order of the spin flip
diffusion length or longer. We examine this suggestion using first-principles
quantum mechanical scattering calculations to study the decay of orbital
currents injected from an orbitally-polarized lead into thermally disordered
bulk systems of selected transition metals. We find that the decay occurs over
only a few atomic layers. On this length scale the orbital current may be
converted into a spin current if the spin Hall angle is sufficiently large, as
for Pt. In Cu, Cr and V with small spin Hall angles, the conversion into a
spin current is negligible in the bulk and significant conversion only occurs
at interfaces.

Introduction. The orbital analogue of the spin Hall effect (SHE) [1, 2, 3, 4],
the orbital Hall effect (OHE) [5, 6, 7, 8], is predicted to have significantly
higher conversion rates than the SHE, in particular for 3$d$ transition metals
[9, 10]. While experimental results appear consistent with this theoretical
prediction [11, 12, 13, 14], long relaxation lengths need to be invoked for
the orbital currents. For Ti, for example, relaxation lengths of 50-$60\pm
15\,$nm [8] or $47\pm 11\,$nm [15] are extracted from experiment; for Cr,
magneto-optical Kerr effect (MOKE) experiments are consistent with a
relaxation length of $6.6\pm 0.6\,$nm [14] while spin-orbit torque (SOT)
experiments yield a value of $6.1\pm 1.7\,$nm [16]. Such large values are at
variance with the standard picture of orbital quenching in solids according to
which orbitally nondegenerate states contain equal amounts of $+m$ and $-m$
character because of time-reversal symmetry [17, 18]. Noting that all states
in a solid are orbitally nondegenerate at room temperature because thermal
disorder eliminates all rotational symmetry, we expect that an orbitally
polarized state $|e^{+im\phi}\rangle$ injected into a material at room
temperature can only hop to states containing equal amounts of
$|e^{+im\phi}\rangle$ and $|e^{-im\phi}\rangle$ character so the orbital
polarization is quenched on the length scale of the hopping.

Although it has been argued in a number of recent publications that $l_{\rm
of}$, the orbital analogue of the spin-flip diffusion length $l_{\rm sf}$, is
short [19, 20, 21], the methods used in these papers precluded quantitative
estimates of $l_{\rm of}$, in particular the effect of various types of
disorder. It is the purpose of this work to illuminate the region between
theoretical studies on perfectly crystalline and defect-free bulk materials on
the one hand and macroscopic experiments on inhomogeneous materials at room
temperature on the other, by determining the orbital relaxation length using
quantum mechanical scattering calculations for thermally disordered solids.

Figure 1: A spin-polarized current injected from a spin-polarized left lead
($\mathcal{L}$) shows a clear exponential decay in the scattering region
($\mathcal{S}$) of thermally disordered Pt at T=300K. The inset indicates that
the pink-shaded part of the scattering region is lead-like [22].

Method.—The calculations involve large scattering geometries (thousands of
atoms) and take temperature-induced disorder into account in the adiabatic
approximation [23]. The formalism implemented in the Twente Quantum Transport
(TQT) [24, 25, 26, 27] code considers the disordered material of interest as a
scattering region ($\mathcal{S}$) sandwiched between semi-infinite crystalline
left ($\mathcal{L}$) and right ($\mathcal{R}$) leads, see the inset to Fig. 1.
Solution of the time-independent Schrödinger equation by matching the
wavefunctions at the $\mathcal{L|S}$ and $\mathcal{S|R}$ interfaces takes the
form of a set of inhomogeneous linear equations [28]. The resulting wave
function contains all information on the linear response of the system to an
external potential [29, 25]. On computing the solution in a basis of tight-
binding linearized muffin-tin orbitals (TB-LMTO) [30, *Andersen:85,
*Andersen:prb86], we use the continuity equation to extract interatomic (spin
and orbital) currents from the wave function [26]. The flux $j$ of the
property $\hat{X}$ between atoms $P$ and $Q$ is

$\\!\\!\\!j^{PQ}_{X}=\frac{1}{i\hbar}\left[\left<\Psi_{P}\right|\hat{X}_{P}\hat{H}_{PQ}\left|\psi_{Q}\right>-\left<\Psi_{Q}\right|\hat{H}_{QP}\hat{X}_{Q}\left|\Psi_{Q}\right>\right],$
(1)

where $\hat{X}=1$ leads to an expression for the charge current,
$\hat{X}=\hat{S}_{\alpha}$ yields the spin current polarized in the $\alpha$
direction and $\hat{X}=\hat{L}_{\alpha}$ results in the orbital current
polarized in the $\alpha$ direction, a current of orbital angular momentum
[33]. If the operator $\hat{X}$ is such that
$\left<\Psi_{P}\right|\hat{X}\left|\Psi_{Q}\right>\neq 0$, this definition of
the flux is not well defined since the continuity equation requires local
conservation. In particular, we cannot use the definition of the orbital
angular momentum in the so-called “modern theory of orbital polarization” [34,
35] to define a current in the scattering approach; instead, we use the atom-
centered approximation of the orbital angular momentum operators to compute
orbital angular momentum currents. Interatomic currents were previously used
to compute the spin Hall angle (SHA), spin-flip diffusion length, spin memory
loss and the effect of interface disorder on those parameters [36, 26, 37, 38,
39, 40, 41].

To construct a spin-polarized current, the left lead is modified to shift the
energy of one of the spin channels so that the Fermi surface becomes
completely spin-polarized [26]. A charge current injected from the left lead
carries a spin current that decays exponentially in the scattering region of
thermally disordered bulk diffusive Pt, Fig. 1. At room temperature, a value
of the spin-flip diffusion length $l_{\mathrm{sf}}\approx 5.2$ nm [26, 38] can
be extracted directly. The setup for this calculation is illustrated in the
inset to Fig. 1. Constructing a spin-polarized system this way requires
introducing a Zeeman term in the Hamiltonian for the left lead, albeit only in
the spin degree of freedom

$\hat{H}=\hat{B}\cdot\hat{S},$ (2)

where the amplitude of $\hat{B}$ is chosen to completely spin-polarize the
Fermi surface.

To compute an orbital relaxation length, we want to inject an orbitally
polarized current. To do so, we introduce an orbital Zeeman term into the
Hamiltonian for the left lead

$\hat{H}=\hat{B}\cdot\hat{L},$ (3)

where the choice of amplitude of $\hat{B}$ is not so straightforward as in the
spin case. Unlike the spin degrees of freedom which are only coupled weakly
through the spin-orbit coupling (SOC), modelled in our case as

$\hat{H}_{\rm so}=\xi\hat{L}\cdot\hat{S},$ (4)

the orbital degrees of freedom are much more strongly coupled by the crystal
potential so that introducing a significant Zeeman orbital term completely
distorts the electronic band structure. The resulting computational setup is
sketched in the inset to Fig. 2 [22] and used to inject a current of orbital
angular momentum into thermally disordered, diffusive Pt. The most striking
feature of the figure is a very short orbital current decay length $l_{\rm
of}$.

Figure 2: An orbitally-polarized current injected from an orbitally-polarized
left lead decays on a length scale of only a few atomic layers in room-
temperature thermally disordered Pt [22].

Computational details.—The workflow for these calculations is as follows. We
begin with a standard density functional theory (DFT) calculation for each
transition metal system using the questaal code [42] with the von Barth-Hedin
LDA functional [43] in combination with a TB-LMTO $spd$ basis and a 21$\times
21\times$21 $\mathbf{k}$-point grid. The Kohn-Sham [44] atomic sphere [30,
*Andersen:85, *Andersen:prb86] potential is extracted and used to construct
the TB-LMTOs in the scattering calculation. The wave-function matching problem
[28], which has been shown to be equivalent to the more conventional non-
equilibrium Green’s function (NEGF) method [29], is set up and solved in the
TQT code. Periodic boundary conditions in the $x$ and $y$ directions ($z$ is
the transport direction) require sampling the two-dimensional Brillouin zone
with a sufficiently dense $\mathbf{k}$-point grid. From previous calculations
[26], a grid equivalent to 160$\times$160 $\mathbf{k}$-points for a 1$\times$1
unit cell is more than sufficient. Here we have used a 29$\times$29
$\mathbf{k}$-point grid with a 5$\times$5 supercell, which is equivalent to
145$\times$145 for a 1$\times$1\. Thermal disorder is modelled by randomly
displacing the atoms from their equilibrium positions in the scattering region
with a Gaussian distribution of displacements. We take multiple configurations
of the random disorder and use the average to simulate the thermodynamic mean.
In practice, 10 configurations is usually sufficient to obtain results
uniquely distinguishable from the noise inherent in the random displacements.
The root mean square disorder parameter is iteratively chosen so that the
resistivity of the system is as close to the experimental room temperature
resistivity as desired.

Results.—The orbital current injected into thermally disordered Pt decays very
rapidly as seen in Fig. 2. The orbital polarization is unity inside the left
lead (and the leftmost, pink shaded part of the scattering region [22]) but
decreases to 2.5% of this value within three atomic layers in diffusive Pt.
Fitting an exponential curve to these results is not very meaningful and would
result in an “orbital diffusion length” shorter than $1\,$nm. On plotting the
spin and orbital currents together in Fig. 3, it can be seen that the orbital
current is converted into a spin current within a few layers in disordered Pt.
In the few atomic layers where the orbital current decays rapidly, the spin
current increases from around 0.06 to 0.08. The polarization of this spin
current is not tremendously high, around 8% at its peak. It decays on a length
scale consistent with the spin-flip diffusion length $l_{\rm Pt}\equiv l_{\rm
sf}^{\rm Pt}\approx 5.2\,$nm.

Figure 3: A current $j_{l}$, orbitally-polarized in the $z$ direction
($j_{lz}$) and injected into thermally disordered Pt in the $z$ direction
($j_{lz}^{z}$), is converted into a spin current $j_{s}$ within a few atomic
layers after which $j_{s}$ decays on a length scale consistent with the spin-
flip diffusion length of Pt, $l_{\rm Pt}\equiv l_{\rm sf}^{\rm Pt}$. The inset
shows how the degree of orbital polarization depends on the size of the
orbital Zeeman splitting in the lead. Figure 4: The orbitally-polarized
current $j_{l}$ injected into thermally disordered Pt decays within a few
atomic layers.

The effect of the Pt SOC can be studied by switching it off in the scattering
region by setting $\xi$ to zero in (4), see Fig. 4. The short decay length of
the orbital current remains unchanged, of order a few atomic layers. The spin
current is zero throughout the scattering region since in the absence of SOC
no orbital current is converted into spin current. The sizable difference in
the orbital polarization in the lead in the absence of SOC ($\sim 14\%$
without versus $\sim 21\%$ with SOC) comes from the slight change in the band
structure when the SOC is turned off. This does not affect the principal
result that $l_{\rm of}$ is very short, of the order of a few atomic layers.

An unphysically large value of $B=40\,$Ry in (3) was used to achieve the high
degree of orbital polarization seen in Fig. 2 (where the Fermi energy in the
lead was shifted to maximize it) and was chosen for illustrative purposes. The
short length scale of the decay is unchanged on using a more reasonable value
of $B=0.05\,$Ry as in Fig. 3. To test the effect of the orbital Zeeman
splitting in the left lead, we calculated the orbital current injection for a
number of values of $B$ in Pt. The orbital polarization achievable with these
values of $|B|$ is shown in the inset to Fig. 3. For all polarization values,
the orbital current decays to less that $1/e$ times its original value within
three atomic layers.

Figure 5: The orbitally-polarized current injected from an orbitally-polarized
left lead decays on a length scale of only a few atomic layers in thermally
disordered Cr, V and Cu. For Cr, orbitally-polarizing fields have been used in
the transport direction ($j_{lz}^{z}$) and perpendicular to it ($j_{ly}^{z}$).
The orbital relaxation length is very short in all cases.

As shown in Fig. 5, we observe a short decay length not only for Pt, but also
for Cr, V and Cu. Although different amplitudes of $B$ can produce very
different values of the orbital polarization, the orbital current decays to
values indistinguishable from the fluctuations arising from disorder sampling
within a few atomic layers in all cases. This indicates that the length scale
$l_{\rm of}$ does not depend on the orbital Hall conductivity which we find to
be large for Cr and V but small for Cu [10].

Even a magnetic field in (3) as small as $B=0.007$ Ry ($\approx 100$ meV)
polarizes the lead sufficiently to allow the rapid decay of the orbital
current injected into V to be clearly registered. In this case, we expect the
electronic structure mismatch to be minimal. For the intermediate values of
$B=0.05$ Ry for Cu and $B=10$ Ry for Cr we find the same behaviour; in all
cases the orbital current is very short-ranged. As in Pt, the orbital current
injected into Cr and V generates a spin current through the spin-orbit
coupling. However, the conversion rates in Cr and V are so low that the spin
currents generated by the orbital-to-spin conversion are practically
indistinguishable from the noise inherent in the disordered calculations.

Fitting an exponential curve to extract values for $l_{\rm of}$ for the
different materials is not possible. We can obtain rough estimates by counting
the number of layers over which the polarization drops to $1/e$ times its
value in the lead. As shown in Table 1, all $l_{\rm of}$ values determined in
this way are $<1$ nm. The range shown for Cr follows from the slight
difference between $y$ and $z$ orbital polarizations; in both cases, the value
of $l_{\rm of}$ is very small.

Table 1: $l_{\mathrm{of}}$ values for Pt, Cu, V and Cr extracted from the decay of orbital polarization inside the scattering region. The number of atomic layers after which the polarization is below the value of $1/e$ times the value in the lead is counted and multiplied by the atomic interlayer distance to obtain an upper bound for $l_{\mathrm{of}}$. |  $1/e$ at layer # |  $d$ (nm) per atomic layer |  $l_{\mathrm{of}}$ (nm)
---|---|---|---
Pt | 3 | 0.23 | 0.6
Cu | 2 | 0.21 | 0.4
V | 3 | 0.14 | 0.5
Cr | 2–5 | 0.15 | 0.3–0.8

Discussion.— Experiments measuring the orbital Hall effect have been
interpreted in terms of long orbital diffusion lengths whereas the results
presented here show that the polarization of an orbital current injected into
a bulk metal drops off within a few atomic layers. The results for Pt show
that the orbital current is converted into a spin current on the same short
length scale. One possible explanation of the long length scale extracted from
experiment is that what is being observed is the long $l_{\rm sf}$ over which
the spin current resulting from the injected orbital current decays rather
than the orbital current itself.

Exactly how the orbital current gets converted into a spin current in systems
containing the light elements Cr and Ti is unclear. One possible explanation
is that there is some orbital-to-spin conversion taking place at an interface
or surface facilitated by the reduced symmetry of the interface.
Alternatively, a large spin-Hall current may be generated directly at the
interface where the $3d$ materials simply serve to facilitate symmetry
breaking by being different from the ferromagnetic material with which they
are interfaced [36]. Though the precise mechanism is probably dependent on the
experiment being performed, it is well known that the effect of SOC is greatly
enhanced by the symmetry lowering at an interface. The results presented here
show that the orbital current can be converted into a spin current within a
few atomic layers and this process could be compatible with an interface-
mediated mechanism. One conclusion we can draw from the results presented here
is that only a few atomic layers are needed to convert almost all of the
orbital current into a spin current and a thick slab of orbital-to-spin
conversion material may not be needed.

Our finding that $l_{\rm of}$ is very short is not necessarily bad news. The
origin of the torque exerted by injecting a spin or orbital current from a
nonmagnetic into a magnetic material is largely a matter of interpretation;
for example, no change has to be made to electronic structure based computer
codes that calculate SOT to include the OHE; it is included automatically. The
important insight from OHE studies is that the search for efficient SOT need
not be restricted to heavy metals with large spin-orbit coupling; the reduced
bandwidth and frequently large state density at the Fermi energy of 3$d$
elements being frequently more important [45, 46, 47]. One consequence of a
very short $l_{\rm of}$ is that the orbital-Hall conductivity is not a useful
figure of merit in the search for better SOT materials. Another is that
because the SOT is essentially a property of an A$|$B interface, it cannot be
factored into separate A and B parts.

Conclusion.—By computing the decay of orbital currents injected from an
orbitally polarized lead, we have shown that the orbital relaxation length in
transition metals is as short as might be expected from classical orbital
quenching arguments. In all cases studied here, the orbital polarization
decays to less than $1/e$ times its injected value within less than a
nanometer. It does not matter if the polarization is small or large. The decay
length does not depend on the spin Hall angle ($l_{\rm of}$ is short for both
Pt with a large SHA and for Cu with a negligible SHA) or the orbital Hall
angle ($l_{\rm of}$ is short for V and Cr as well as for Cu). The conversion
of the orbital current into a spin current does significantly depend on the
spin Hall angle, as expected. These results contradict the interpretation of
experimental findings in terms of significantly longer length scales. One
possible explanation is that the experiments are actually measuring the spin-
flip diffusion length instead. In the case of Pt, we have shown how the
injected orbital current is converted into a spin current within a few layers
after which it decays as we expect for the spin current.

##  I Acknowledgement

This work was sponsored by NWO Domain Science for the use of supercomputer
facilities.

## References

  * Dyakonov and Perel [1971] M. I. Dyakonov and V. I. Perel, Phys. Lett. A 35, 459 (1971).
  * D’yakonov and Perel [1971] M. I. D’yakonov and V. I. Perel, Zh. Eksp. Teor. Fiz. 13, 657 (1971), [JETP Letters-USSR 13, pp 467-469 (1971)].
  * Hirsch [1999] J. E. Hirsch, Phys. Rev. Lett. 83, 1834 (1999).
  * Valenzuela and Tinkham [2006] S. O. Valenzuela and M. Tinkham, Nature (London) 442, 176 (2006).
  * Tanaka  _et al._ [2008] T. Tanaka, H. Kontani, M. Naito, T. Naito, D. S. Hirashima, K. Yamada, and J. Inoue, Phys. Rev. B 77, 165117 (2008).
  * Kontani  _et al._ [2009] H. Kontani, T. Tanaka, D. S. Hirashima, K. Yamada, and J. Inoue, Phys. Rev. Lett. 102, 016601 (2009).
  * Go  _et al._ [2018] D. Go, D. Jo, C. Kim, and H.-W. Lee, Phys. Rev. Lett. 121, 086602 (2018).
  * Choi  _et al._ [2023] Y.-G. Choi, D. Jo, K.-H. Ko, D. Go, K.-H. Kim, H. G. Park, C. Kim, B.-C. Min, G.-M. Choi, and H.-W. Lee, Nature (London) 619, 52 (2023).
  * Jo  _et al._ [2018] D. Jo, D. Go, and H.-W. Lee, Phys. Rev. B 98, 214405 (2018).
  * Salemi and Oppeneer [2022] L. Salemi and P. M. Oppeneer, Phys. Rev. Materials 6, 095001 (2022).
  * Sala and Gambardella [2022] G. Sala and P. Gambardella, Phys. Rev. Research 4, 033037 (2022).
  * Rothschild  _et al._ [2022] A. Rothschild, N. Am-Shalom, N. Bernstein, M. Meron, T. David, B. Assouline, E. Frohlich, J. Xiao, B. Yan, and A. Capua, Phys. Rev. B 106, 144415 (2022).
  * Zhang  _et al._ [2022] J. Zhang, H. Xia, X. Zhang, Z. Yan, Y. Zhai, J. Chi, H. Xu, Y. Zuo, and L. Xi, Appl. Phys. Lett. 121, 172405 (2022).
  * Lyalin  _et al._ [2023] I. Lyalin, S. Alikhah, M. Berritta, P. M. Oppeneer, and R. K. Kawakami, Phys. Rev. Lett. 131, 156702 (2023).
  * Hayashi  _et al._ [2023] H. Hayashi, D. Jo, T. Gao, S. Haku, Y. Mokrousov, H.-W. Lee, and K. Ando, Communications Physics 6, 32 (2023).
  * Lee  _et al._ [2021] S. Lee, M.-G. Fang, D. Go, D. Kim, J.-H. Kang, T. Lee, G.-H. Lee, J. Hang, N. J. Lee, Y. Mokrousov, S. Kim, K.-J. Kim, K.-J. Lee, and B.-G. Park, Communications Physics 4, 234 (2021).
  * Mohn [2003] P. Mohn,  _Magnetism in the Solid State: An Introduction_ , edited by M. Cardona, P. Fulde, K. von Klitzing, and H.-J. Queisser, Springer Series in Solid-State Sciences, Vol. 134 (Springer, Berlin, 2003).
  * Tinkham [2003] M. Tinkham,  _Group Theory and Quantum Mechanics_ , Dover Books on Chemistry (Dover Publications, Minneola, New York, 2003).
  * Salemi  _et al._ [2021] L. Salemi, M. Berritta, and P. M. Oppeneer, Phys. Rev. Materials 5, 074407 (2021).
  * Belashchenko  _et al._ [2023] K. D. Belashchenko, G. G. Baez Flores, W. Fang, A. A. Kovalev, M. van Schilfgaarde, P. M. Haney, and M. D. Stiles, Phys. Rev. B 108, 144433 (2023).
  * Urazhdin [2023] S. Urazhdin, Phys. Rev. B 108, L180404 (2023).
  * foo [a] Because of how currents are extracted from the scattering wave functions [26], $j(z)$ is not calculated in the leads. In order to visualise $j(z)$ in the left-hand lead, the interface at $z=0$ is artificially shifted into the scattering region to a positive value of $z$. This is indicated in the figures by the pink shading.
  * Liu  _et al._ [2015] Y. Liu, Z. Yuan, R. J. H. Wesselink, A. A. Starikov, M. van Schilfgaarde, and P. J. Kelly, Phys. Rev. B 91, 220405(R) (2015).
  * Xia  _et al._ [2006] K. Xia, M. Zwierzycki, M. Talanana, P. J. Kelly, and G. E. W. Bauer, Phys. Rev. B 73, 064420 (2006).
  * Starikov  _et al._ [2018] A. A. Starikov, Y. Liu, Z. Yuan, and P. J. Kelly, Phys. Rev. B 97, 214415 (2018).
  * Wesselink  _et al._ [2019] R. J. H. Wesselink, K. Gupta, Z. Yuan, and P. J. Kelly, Phys. Rev. B 99, 144409 (2019).
  * Nair and Kelly [2021] R. S. Nair and P. J. Kelly, Phys. Rev. B 103, 195406 (2021).
  * Ando [1991] T. Ando, Phys. Rev. B 44, 8017 (1991).
  * Khomyakov  _et al._ [2005] P. A. Khomyakov, G. Brocks, V. Karpan, M. Zwierzycki, and P. J. Kelly, Phys. Rev. B 72, 035450 (2005).
  * Andersen and Jepsen [1984] O. K. Andersen and O. Jepsen, Phys. Rev. Lett. 53, 2571 (1984).
  * Andersen  _et al._ [1985] O. K. Andersen, O. Jepsen, and D. Glötzel, in  _Highlights of Condensed Matter Theory_ , International School of Physics ‘Enrico Fermi’, Varenna, Italy, edited by F. Bassani, F. Fumi, and M. P. Tosi (North-Holland, Amsterdam, 1985) pp. 59–176.
  * Andersen  _et al._ [1986] O. K. Andersen, Z. Pawlowska, and O. Jepsen, Phys. Rev. B 34, 5253 (1986).
  * foo [b] For convenience, we measure the charge density in units of the electron charge -$e$ (where $e$ is a positive quantity), the spin density in units of $\hbar/2$ and the orbital density in units of $\hbar$ so they all become particle densities and the fluxes become particle fluxes.
  * Thonhauser  _et al._ [2005] T. Thonhauser, D. Ceresoli, D. Vanderbilt, and R. Resta, Phys. Rev. Lett. 95, 137205 (2005).
  * Pezo  _et al._ [2022] A. Pezo, D. G. Ovalle, and A. Manchon, Phys. Rev. B 106, 104414 (2022).
  * Wang  _et al._ [2016] L. Wang, R. J. H. Wesselink, Y. Liu, Z. Yuan, K. Xia, and P. J. Kelly, Phys. Rev. Lett. 116, 196602 (2016).
  * Gupta  _et al._ [2020] K. Gupta, R. J. H. Wesselink, R. Liu, Z. Yuan, and P. J. Kelly, Phys. Rev. Lett. 124, 087702 (2020).
  * Nair  _et al._ [2021] R. S. Nair, E. Barati, K. Gupta, Z. Yuan, and P. J. Kelly, Phys. Rev. Lett. 126, 196601 (2021).
  * Gupta  _et al._ [2021] K. Gupta, R. J. H. Wesselink, Z. Yuan, and P. J. Kelly, Phys. Rev. B 104, 205426 (2021).
  * Gupta  _et al._ [2022] K. Gupta, R. Liu, R. J. H. Wesselink, Z. Yuan, and P. J. Kelly, Phys. Rev. B 106, 115425 (2022).
  * Liu  _et al._ [2022] R. Liu, K. Gupta, Z. Yuan, and P. J. Kelly, Phys. Rev. B 106, 014401 (2022).
  * Pashov  _et al._ [2020] D. Pashov, S. Acharya, W. R. Lambrecht, J. Jackson, K. D. Belashchenko, A. Chantis, F. Jamet, and M. van Schilfgaarde, Comput. Phys. Commun. 249, 107065 (2020).
  * von Barth and Hedin [1972] U. von Barth and L. Hedin, J. Phys. C: Sol. State Phys. 5, 1629 (1972).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Gupta [2019] K. Gupta,  _Disentangling interfaces and bulk in spin transport calculations_ , Ph.D. thesis, University of Twente, The Netherlands (2019).
  * Amin  _et al._ [2019] V. P. Amin, J. Li, M. D. Stiles, and P. M. Haney, Phys. Rev. B 99, 220405(R) (2019).
  * Nair [2021] R. S. Nair,  _Spin-orbit-coupling induced lateral spin transport from first principles_ , Ph.D. thesis, University of Twente, The Netherlands (2021).