# Deep-Learning Database of Density Functional Theory Hamiltonians for Twisted
Materials

Ting Bao  State Key Laboratory of Low Dimensional Quantum Physics and
Department of Physics, Tsinghua University, Beijing, 100084, China  These
authors contributed equally  Runzhang Xu<EMAIL_ADDRESS>State Key
Laboratory of Low Dimensional Quantum Physics and Department of Physics,
Tsinghua University, Beijing, 100084, China  These authors contributed equally
He Li  State Key Laboratory of Low Dimensional Quantum Physics and Department
of Physics, Tsinghua University, Beijing, 100084, China  Institute for
Advanced Study, Tsinghua University, Beijing 100084, China  Xiaoxun Gong
School of Physics, Peking University, Beijing 100871, China  Zechen Tang
State Key Laboratory of Low Dimensional Quantum Physics and Department of
Physics, Tsinghua University, Beijing, 100084, China  Jingheng Fu  State Key
Laboratory of Low Dimensional Quantum Physics and Department of Physics,
Tsinghua University, Beijing, 100084, China  Wenhui Duan
<EMAIL_ADDRESS>State Key Laboratory of Low Dimensional Quantum Physics
and Department of Physics, Tsinghua University, Beijing, 100084, China
Institute for Advanced Study, Tsinghua University, Beijing 100084, China
Frontier Science Center for Quantum Information, Beijing, China  Beijing
Academy of Quantum Information Sciences, Beijing 100193, China  Yong Xu
<EMAIL_ADDRESS>State Key Laboratory of Low Dimensional Quantum
Physics and Department of Physics, Tsinghua University, Beijing, 100084, China
Frontier Science Center for Quantum Information, Beijing, China  RIKEN Center
for Emergent Matter Science (CEMS), Wako, Saitama 351-0198, Japan

###### Abstract

Moiré-twisted materials have garnered significant research interest due to
their distinctive properties and intriguing physics. However, conducting
first-principles studies on such materials faces challenges, notably the
formidable computational cost associated with simulating ultra-large twisted
structures. This obstacle impedes the construction of a twisted materials
database crucial for data-driven materials discovery. Here, by using high-
throughput calculations and state-of-the-art neural network methods, we
construct a Deep-learning Database of density functional theory (DFT)
Hamiltonians for Twisted materials named DDHT. The DDHT database comprises
trained neural-network models of over a hundred homo-bilayer and hetero-
bilayer moiré-twisted materials. These models enable accurate prediction of
the DFT Hamiltonian for these materials across arbitrary twist angles, with an
averaged mean absolute error of approximately 1.0 meV or lower. The database
facilitates the exploration of flat bands and correlated materials platforms
within ultra-large twisted structures.

##  I Introduction

Figure 1: Schematic diagram of construction and application of DDHT database.
Uniformly distributed interlayer slidings together with random atomic
perturbations are introduced in the non-twisted bilayer supercells to
construct training datasets. Graph neural networks are employed here, using
materials structures as input and their DFT Hamiltonians as output. Deep
neural network models are trained by DFT calculation results of training
structures, and then generalized to study diverse moiré-twisted structures
with varying twist angles and constituent materials selected from C2DB. From
the DFT Hamiltonians predicted by DeepH, many physical properties can be
derived, including band structure, band topology, projected density of states,
etc.

The discovery of correlated states in magic-angle twisted bilayer graphene
(TBG) [1, 2] has sparked great research interest in the emergent physics and
exotic properties in moiré-twisted two-dimensional (2D) materials. In moiré
twist systems, the relative twist rotation between the two constituent van der
Waals (vdW) layers forms moiré patterns with tunable periodicity, generating a
moiré potential field that strongly influences the electronic structures and
introduces exotic quantum effects such as correlated states, band topology,
spin textures, and ferroelectric polarization [3]. These twist-induced
interesting underlying physics have been intensively studied on TBG, twisted
bilayer CrI3[4, 5], twisted bilayer BN (TBBN)[6, 7], twisted transition metal
dichalcogenides (TTMDs) [8] and so on. Other twist-induced properties, like
moiré excitons [9, 10, 11, 12], Wigner crystal states [13], and spin- and
bond-density waves [14, 15] further expand the knowledge of twisted materials.
These intriguing effects and physics, coupled with the rich external
modulations, position twisted materials not only as ideal quantum simulators
[16] in condensed-matter research but also as practical platforms for
twistronic applications.

Despite the growing interest in twisted materials and the related physics, the
field faces fundamental challenges from both experimental and theoretical
perspectives [17]. In experiments, the synthesis of high-quality twisted
samples and the precise control of twist angles are rather challenging. From
theoretical and calculation perspective, various kinds of numerical and
theoretical methods, such as DFT, tight-binding method, and continuum model,
are employed to make predictions for guiding experiments. However, for twisted
materials, all these theoretical methods encounter the vital efficiency-
accuracy dilemma. For instance, DFT methods scales as $\mathcal{O}(N^{3})$
with $N$ being number of atoms per periodic unit cell, which limits the sizes
of twisted structures to be studied. The empirical tight-binding methods or
continuum models are much more efficient, but suffers from the critical
problem of low accuracy. The revolutionary success of artificial intelligence,
especially deep learning neural networks, in many fields [18, 19, 20] brings
solution to this dilemma. In previous work DeepH [21], Li et.al. show that a
neural network (NN) model trained on small structures can be generalized to
tackle the DFT Hamiltonians of large twisted structures at significantly
reduced computational cost. However, only TBG, twisted bilayer bismuthene and
twisted MoS2 are studied using DeepH, which is not enough to validate the
effectiveness. A more systematic learning on DeepH methods is in urgent need,
which is not only to check the validness and universality of DeepH but also to
meet the demand of getting the DFT-level Hamiltonians for more materials and
boost further research.

Here, we report the construction of a deep-learning database of DFT
Hamiltonians for twisted materials named DDHT as shown in Fig. 1. DDHT
currently covers DeepH NN models for 124 types of 2D homo-bilayer and 5 types
of hetero-bilayer twist materials and can efficiently predict accurate DFT
Hamiltonians of their twisted structures across arbitrary twist angles. For
each type of material, the NN model is trained by the state-of-the-art deep-
learning representation of DFT Hamiltonians implemented in the E3-equivariant
DeepH method [22, 21]. We demonstrate the high reliability and DFT-comparable
accuracy of all NN models in our database by the very low validation loss
achieved in the NN training process and using the DFT-calculated Hamiltonians
at selected large twist angles as benchmark. Furthermore, taking 2D black
phosphorene (bP) as an exhibition, we use the corresponding NN model in DDHT
to obtain electronic band structure and find the ultra-flat valence bands in
twisted bilayer bP down to twist angles inaccessed by conventional DFT
methods, showing the effectiveness and efficiency of our database in getting
exotic electronic properties of ultra-large twisted materials. In our
exhibition of twisted bP and probably many other twisted materials systems,
the flat bands usually have bandwidth droping to only a few meV at small-angle
limit ($<$2.0 ∘) and their real space projection can be highly localized,
giving rise to rich exotic properties from doping modulation and strong
correlation. DDHT not only provides an optimal parameter strategy for NN
training and strongly supports the universality of DeepH method, but also
offers prior DFT-level knowledge of Hamiltonians and electronic structures for
researchers in the field of twisted materials and twistronics and thus greatly
accelerate the exploration in this field.

It is worth mentioning that unlike traditional material databases like
Materials Project [23], ICSD [24] or C2DB [25, 26], which includes specific
materials corresponding to discrete points in the materials parameter space,
our DDHT is applicable to twisted structures at arbitrary twist angles and
interlayer fluctuation of provided types of material, covering wide areas in
the materials parameter space with infinite single structure cases inside.
Furthermore, considering each model in the DDHT as an expert and adding a
simple classifier as the gate on types of material, out DDHT itself is
actually serving as a mixture of experts (MoE) model [27] with each NN item
experienced in one type of twisted material. Thus DDHT has dualism of database
and MoE model, we suppose this dualism is a new formunism and will be a
necessity of future databases in the AI era.

##  II Results

###  II.1 Construction of Database

Figure 2: Workflow of constructing deep-learning database. All the steps are
classified into three parts: candidate selection, dataset generation, and
neural network training and validation. The parameters used in construction
are explicitly displayed by typical values for key steps. The eventually
obtained DDHT database can be further used for unlimited cases.

Candidate 2D materials considered in the current database are adopted from
C2DB [25, 26]. They are further screened to construct twisted materials
structures and training datasets (left part of Fig. 2). In total, 124 types of
2D materials (Table S1) in the C2DB are selected following three criteria:
non-magnetic, thermally and dynamically stable, structurally simple ($\leq$ 8
atoms per unit cell, $\leq$ 3 elements). Their bilayer unit cells are
constructed and fully relaxed by DFT calculations using Vienna ab initio
simulation package (VASP) [28, 29] (see Note 1 in Supporting Information for
details of screening and relaxation). This allows us to obtain the optimal
structural parameters, including lattice constants, stacking orders, and
interlayer spacings. These structural parameters will be used to construct the
non-twisted supercell structures in the datasets and the corresponding 124 (5)
types of to-be-predicted moiré-twisted homo-bilayer (hetero-bilayer)
structures with varying twist angles $\theta$s. These candidates not only
cover most of the intensively studied 2D materials and their twisted
structures, such as graphene, black phosphorene (bP), and transition metal
dichalcogenides (TMDCs), but also include most light elements in the periodic
table (Fig. S1).

The schematic diagram of constructing and using our database is demonstrated
in Fig. 1 using TBG as an example material. Here, the construction involves
both the preparation of appropriate datasets and the training of deep neural
network models. The usage refers to the prediction of DFT Hamiltonians for
twisted materials with any given twist angle as well as various kinds of
derived physical properties. We first discuss the preparation of training
datasets. Moiré-twisted materials host rich structural diversity in terms of
twist angle and constituent materials, which should be captured by atomic
structures in the training datasets for accurate ${\mathcal{R}}\mapsto
H_{\text{DeepH-E3}}$ mapping by the trained models. This requirement is
satisfied by preparing atomic structures in datasets respectively (or
separately) for each type of twisted materials, and introducing diverse enough
variations to local atomic positions and interlayer stackings to cover all
possible local environments in twisted bilayer structures. For instance, TBG
and TBBN have different constituent 2D materials (graphene and BN), and the
non-twisted supercells in their datasets are generated using graphene and BN
bilayer unit cells respectively, while the twisted heterobilayer graphene/BN
have heterogeneity in layers and thus uses graphene/BN unit cells for
generating the dataset. The relative rotations induced by twisting between the
two vdW layers can be treated as locally approximated interlayer shear
slidings combined with perturbations to atomic positions in non-twisted
bilayer structures. Therefore, for each type of twisted material, the atomic
structures in its dataset are prepared (”Training data generation” part in
Fig. 2) by 1) making non-twisted supercells with sufficient large lateral size
from the most optimal bilayer unit cells to include all neighbors within cut-
off radius $R_{c}$ of local atomic basis and also to avoid interactions from
periodic images, and 2) introducing uniform distributed interlayer slidings
between the two vdW layers and random perturbations to the atom positions in
these non-twisted supercells to cover the local structural diversity in
twisted structures, with the mesh density of slidings and the upper limit of
perturbations determined based on balanced accuracy and efficiency. The local-
atomic-basis Hamiltonians $H_{\text{DFT}}$ of all these non-twisted,
perturbed, and shifted supercell structures in the dataset are calculated
self-consistently by first-principles DFT methods and further processed to
extract the decomposed hopping matrix segments $\mathcal{H}_{ij}$ that are in
one-to-one correspondence with local atomic structures $\mathcal{R}_{ij}$
within $R_{c}$. After representing $\mathcal{R}_{ij}$ by crystal graphs [30],
we finally prepared the complete training dataset for each twisted material by
incorporating these $\mathcal{H}_{ij}$–$\mathcal{R}_{ij}$ pairs.

Based on each training dataset, an E(3)-equivariant deep neural network model
is built and trained (”Neural network training and validation” part in Fig. 2)
using the state-of-the-art DeepH-E3 code [22, 21]. To ensure the high accuracy
and generalizability of the model, the training process runs for a sufficient
number of epochs and stops when the validation loss reaches a strict criterion
(see Sec. METHODS for details). Based on a self-made high-throughput
computation framework, the dataset preparation and the model training for all
$124+5$ types of twisted materials are carried out in high efficiency, which
in the end output $124+5$ E(3)-equivariant deep neural network models and
their corresponding non-twisted bilayer unit cells that both constitute the
database. The DFT calculations in preparation of training datasets and the
process of model training are both, to some extend, computationally expansive.
Typically, for each type of material, they will take 1-3 days on a 64-core CPU
node and 2-7 days on a NVIDIA V100 GPU card. However, the resulting well-
trained DeepH-E3 models in our database have strong generalizability and
universality over a wide range of $\theta$s (from 21.79∘ all the way down to
less than 1∘) and are capable of predicting DFT-accuracy $H_{\text{DeepH-E3}}$
of ultra-large (ultra small-angle) twisted structures beyond the reach of
common DFT codes with an acceptable cost. For example, the computational time
of predicting $H_{\text{DeepH-E3}}$ scales only linearly with the number of
atoms in lattice cells, allowing for the calculation of DFT-accuracy
Hamiltonian for magic-angle TBG in just about 100 seconds [21].

Other than the direct usage of our database in predicting
$H_{\text{DeepH-E3}}$ by direct neural-network mapping for twisted structures
with arbitrary $\theta$s and situation of interlayer fluctuation, several
important physical properties can be further calculated from the predicted
$H_{\text{DeepH-E3}}$ sparse matrices (dashed box on the right in Fig. 1),
which may include the eigen values and states, band structures, total and
projected density of states (DOS), and possible correlated flat bands near the
Fermi level and their real-space projections. It should be emphasized that our
database, serving as a MoE for DFT Hamiltonian prediction, offers an ultra-
efficient way (enhanced by orders of magnitude [21]) to bypass the most
efficiency-limiting DFT calculations in computational research flow in twisted
materials and still output Hamiltonians with similar accuracy.

Figure 3: The statistics over MAEs between the DeepH-E3-predicted
$H_{\text{DeepH-E3}}$ and the DFT-calculated $H_{\text{DFT}}$ of 124(5) types
of twisted homobilayer(heterobilayer) materials covered in our database. (a-d)
shows the MAEs of elementary substances, ternary compounds, binary compounds
(divided into four prototypes AB, $AB_{2}$, $A_{2}B_{3}$ and $AB_{3}$) and
hetero structures, respectively, by bar plot. The height of each bar displays
the MAE for one specific 2D material, which is obtained by averaging
Hamiltonian errors over all atoms and orbitals for its twisted structures with
the largest two commensurate twist angles $\theta=$ 21.79∘ (2-1) and 13.17∘
(3-2), except for $\theta=$ 10.87∘ (7-4) and 8.49∘ (9-5) in rectangular P4
(bP) and $\theta=$ 10.87∘ (7-4) and 8.49∘ (9-5) in square SnS2. The predicted
MAEs for most materials are in the order of sub meV, with only the twisted
heterobilayer Bi2Se3/Bi2Te3 and WS2/WSe2 being exceptional (about 1.2 meV).

###  II.2 Accuracy of the Database

To assess the accuracy and reliability of these DeepH-E3 models in our DDHT
database in predicting Hamiltonians and electronic structures of twisted
materials, we perform comprehensive statistics over mean absolute errors
(MAEs) between the DeepH-E3-predicted and the DFT-calculated Hamiltonians for
selected rigid-twisted structures of all 124 homo-bilayer materials.
Considering the heavy computational demand of DFT Hamiltonians of twisted
materials, we only cover the twisted supercells with either the largest two
commensurate twist angles (e.g. 21.79∘ and 13.17∘ corresponding to twist index
2-1 and 3-2 for hexagonal lattice) or the number of atoms less than 120 in the
MAE statistics. For each type of material, the absolute-error matrices are
calculated by taking the absolute difference between respective hoppings of
the predicted Hamiltonians of its twisted bilayer structures and that of the
DFT calculated ones. These absolute error values are further averaged with
respect to hopping pairs and containing chemical elements, resulting in the
material-wise MAEs as shown in Fig. 3. The overwhelming majority of hopping
MAEs are in the sub-meV order of magnitude, which are comparable to or even
less than commonly recognized errors in general DFT calculations. It’s safe to
say 1) the DeepH-E3 method can be efficiently generalized to much more
materials using the workflow and parameter strategy here in the DDHT 2)
Hamiltonian predictions by models in our DDHT database are indeed at DFT
accuracy. By expanding the MAE heatmaps with respect to orbitals of different
element pairs (see Note 5 in Supporting Information), we find no significant
large MAEs between these explicitly presented orbital-wise hoppings which may
qualitatively weaken the above conclusion.

It’s also important to note that the models in our database work well for
twisted materials both with and without spin-orbit coupling (SOC) in our
demonstrations as long as the datasets (thus the DeepH-E3 models) incorporate
SOC for respective (SOC-embedded) materials.

We further verify the accuracy and effectiveness of the models in our database
by calculating the electronic band structures from the above predicted
Hamiltonian matrices and presenting the close fit between the predicted energy
bands (energy dispersion relations) and the DFT-calculated ones. Here, two
representative materials, bP and 2H-TaS2, are chosen as demonstration based on
the general research interest in their unique features, namely, the potential
flat bands in twisted rectangular buckled lattice [15] of semiconducting bP
and the interplay among twisting, SOC, charge density wave, and
superconductivity in metallic TaS2 [31]. At large twist angles, the energy
bands obtained from the DeepH-E3 predictions (blue dots) and the DFT
calculations (red curves) are almost identical (left panels in Fig. 4a, b, d,
e) in the energy range near the Fermi level. Similar results are observed in
all other twisted homobilayer materials covered in our database, as shown in
the left panels of Fig. SX and also reported in previous studies [22, 21].

The twisted heterobilayers induce extra complexity in crystal structures and
interatomic hopping relations, which is more challenging for DeepH-E3 models
to give accurate predictions. We test the validity regarding the heterogeneity
in their two constituent layers in several selected twisted heterobilayer
materials, such as Bi2Se3/Bi2Te3, graphene/BN, MoS2/MoSe2, WS2/WSe2, and
MoTe2/WSe2, and include the DeepH-E3 models of these materials in our
database. Our database outputs accurate predictions for all these twisted
heterobilayers considered (Fig. 4c, f and Fig. SXXX-SXXX) and therefore holds
the accuracy and effectiveness toward heterogeneity in constituent layers.
Take twisted strong-SOC Bi2Se3/Bi2Te3 with $\theta$ being 21.79∘ and 13.17∘ as
an example, the predicted and DFT-calculated energy bands have negligible
difference (with MAE in the order of XXX meV), and more importantly, the non-
trivial Z2 topology of both twisted structures is also well reproduced by
prediction. Therefore, as a step forward, our database firstly generalized the
DeepH-E3 method and models to twisted heterobilayers with comparable accuracy
to homobilayer cases.

In addition to the DFT band structures, the predicted Hamiltonians by our
database can be used to calculate the electron density of states (DOS) through
identifying the basis projections of their eigenvectors. The right panels in
Fig. 4a-f display the total (grey curves) and the projected DOS (colored
curves) of twisted bP, 2H-TaS2, and Bi2Se3/Bi2Te3, with DOS of all other
large-angle twisted structures presented in the right panels of Fig. SX2. The
predicted total and projected DOS provides more detailed information on the
distribution of electronic states with respect to different energy levels,
chemical elements, atoms, and orbitals, supplementing the picture of
electronic structures other than dispersion relations demonstrated above.

At the current stage, the above demonstrations of our database in accurately
predicting DFT Hamiltonians and electronic structures cover all $124+5$ types
of twisted bilayer materials considered, which indicates its high
generalizability and reliability as a DFT-level deep-learning model database
that offers a wide set of DeepH-E3 models targeting the prediction
Hamiltonians of twisted materials, effectively bypassing the expensive first-
principles calculations. The research into twisted materials and related areas
can be greatly assisted and accelerated by our database, as well as the
predicted results, through offering guidance and reference prior to costly
investigations and screening for more exotic properties.

Figure 4: The DFT-calculated band structures (red curves) and predicted band
structures (blue dots) for twisted (a, b) bP with $\theta=$ 10.87∘ and 8.49∘,
and (c, d) 2H-phase TaS2, and (e, f) heterobilayer Bi2Se3/Bi2Te3 with
$\theta=$ 21.79∘ and 13.17∘. The twist angles of bP are different due to the
consideration of lattice compression in constructing twisted supercells from
rectangular lattices. The accuracy of both the DDHT models and predicted
electronic structures is demonstrated by almost identical bands from
predictions and DFT calculations.

###  II.3 Prediction of Flat Bands

While our current demonstrations of predicting DFT Hamiltonians and electronic
structures for twisted structures are limited to $\theta$ values as low as
13.17∘ for the majority of the considered materials, the remarkable
generalizability of our database extends the capability to explore twisted
structures with arbitrarily small $\theta$ angles—a domain beyond the
computational reach of most DFT codes while maintaining acceptable
computational costs. To illustrate this capability, we have chosen twisted bP
bilayer structures at various twist angles as exemplars. As shown in the Fig.
4a, as the twist angle of bP approaches smaller values, the time required for
the NN model to perform Hamiltonian inference, which is directly proportional
to the number of atoms in the unit cell, increases rapidly but remains within
the order of minutes. It is noteworthy that the widths of energy band near the
Fermi energy diminish as the twist angle decreases, as shown in Fig. 4b,
revealing the band width of the highest valence band and the average band
width of the top 10 valence bands with respect to the twist angle of bP. The
top two valence bands in Fig. 4a, d narrows their band width along the
X-$\Gamma$ high-symmetry path with $\theta$ decreased by about 1.93∘, and they
will be potentially flattened throughout the whole momentum space in twisted
supercells at smaller $\theta$s. In this study, we examine two bP twisted
supercells with $\theta=$ 2.36∘ (4756 atoms in unit cell) and 1.59∘ (10588
atoms in unit cell), which are relaxed by a machine-learning force-field model
(GAP+R6) for phosphorene [32] to eliminate the unreasonable band crossings
found in ion-clamped structures at small $\theta$s. Obtaining the band
structures of the two bP twisted supercells with thousands of atoms using ab
initio DFT methods consumes unacceptable computational resources, however,
they can be obtained using the DDHT model within days or even hours (mostly
consumed by matrix diagonalization). As expected, the predicted band
structures of the two small-angle bP twisted supercells, as shown in Fig. 5b,
c, clearly display the correlated flat-band feature in their top valence
states, offering a new opportunity to study flat-band physics and properties
in an ab initio manner. Obtaining the band structures of the two bP twisted
structures with thousands of atoms using ab initio DFT methods consumes
unacceptable computational resources, however, they can be obtained using the
DDHT model within days or even hours (mostly spent on matrix diagonalization).

For the twisted bilayer bP structure with twist angle $\theta=$ 2.36∘, our
DDHT model for bP predicts the presence of two stand-alone flat bands (Fig.
5b) separated by 8.4 meV at the valence band top. The band width of the first
(marked by red) and the second valence flat band (marked by blue) are 7.7 and
1.0 meV. Their small band width enables exotic correlated states when they are
partially occupied upon hole doping. Moreover, the closely located flat bands
with respect to energy make the correlated states from different flat bands
accessible through modulation of doping density in a small range. The real-
space projections of these two flat bands (5e) are quite localized and well
separated by half the moiré periodicity, which also support the strong
correlated states. The flat bands and their real-space projections of twisted
bP of smaller $\theta=$ 1.59∘ (Fig. 5c, f) are also predicted by our DDHT
model. In Fig. 5c, multiple stand-alone flat valence bands are identified in
energy range between $\pm$ 0.2 eV and have band width around 1.0 meV,
potentially giving rise to more pronounced correlated phases than in larger-
angle cases.

The accurate and efficient predictions of flat bands and their real-space
projections pave the way for discovering new flat-band materials and future
experimental and theoretical investigations into the correlated physics in
twisted bP and other twisted materials.

Figure 5: Prediction of flat bands in electronic structures of twisted bP
systems. (a) shows the number of atoms per unit cell of twisted structures
(blue bars) and the associated time of Hamiltonian inference using the trained
model in DDHT for different twist angles. The inference time increases almost
linearly with the number of atoms. (b) plots the bandwidth of highest occupied
band, average bandwidth and standard error of 10 bands, which contains 5
highest valance bands and 5 lowest conduction bands, against the twist angle.
The angle values on x-axis of (a, b) use logarithmic scale for clear
demonstration. Predicted band structures and real-space flat-band projections
of twisted bP structures with (c) $\theta=$ 2.36∘ and (d) 1.59∘ as examples.
The top two flat valence bands are outlined by red and blue and their
projections in real space are presented by transparent regions in the same
color. The twisted atomic structures are schematically displayed in the
background.

##  III Discussion

While all the previous discussions have mainly focused on ion-clamped atomic
structures of twisted materials, it should be noted that the DeepH-E3 models
in the DDHT are not limited to predicting only such structures. In fact, our
previous examples on the band-structure prediction of relaxed bP twisted
bilayers demonstrated that the prediction accuracy is not significantly
affected by structural relaxation. This is mainly because the atomic
displacements induced by relaxation are generally below or around 0.2 Å, the
resulting interatomic distances between neighboring atoms after relaxation are
still within the coverage of our training datasets (eg, Fig. S2). By using the
same model that predicts for ion-clamped TBG, the DeepH-E3-predicted energy
bands of relaxed magic-angle TBG fit closely to its DFT-calculated ones [22],
which further support this argument. The obtained MAEs in both cases are in
the sub-meV range, indicating that the models in the DDHT can capture the
structural relaxation accurately. Therefore, we can conclude that the DDHT
models are applicable to a wide range of twisted materials structures, whether
they are ion-clamped or relaxed.

Straightforwardly, the models in the DDHT are highly versatile in predicting
electronic structures of twisted multilayer structures, which can be treated
as combinations of twisted bilayers consist of the same or different
materials. For instance, in a twisted trilayer structure, where layers 1, 2,
and 3 are arranged in a top-to-bottom order, layer 1-2 and layer 2-3
constitute two twisted bilayers. As a result, the Hamiltonian and electronic
structure of this twisted trilayer can be accurately predicted using the same
models as in the bilayer cases. This argument can be extended to twisted
hetero-multilayer structures where the constituent layers are all different
materials. In this case, multiple models trained by heterobilayer training
datasets are required for accurate predictions. The accuracy of these models
has been previously demonstrated by the predictions of twisted heterobilayer
structures including graphene/BN, MoS2/MoSe2, WS2/WSe2, Bi2Se3/Bi2Te3, and
MoTe2/WSe2, and can potentially sustain in predicting twisted multi-
heterolayers.

The DDHT currently covers only a fraction of stable 2D materials (124 out of
about 2132 non-magnetic materials in the C2DB accessed at June 19, 2022), it
will cover more non-magnetic materials in the near-future development, by
using the optimal scheme of dataset preparation and DeepH-E3 model training
discussed in this article. Especially, we will take priority at offering
models involving materials which users of our DDHT concern most. Although the
models in the DDHT can predict arbitrary twist angles, the diagonalization to
obtain band structure, DOS or other properties will consume certain amount of
computational resource. Thus we allow users to upload their results help to
enrich the database. In conclusion, the DDHT will be updated with predictions
for more twist angles, especially for twisted materials and twist angles that
are currently of general research interest. Additionally, the magnetic
materials (about 464 stable in the C2DB) present a unique challenge due to
their spin degree of freedom, which requires the incorporation of both spin
moments and crystal structures as input in learning the neural network
representations of their Hamiltonians. This can be achieved through the
utilization of the recently developed xDeepH method [33], which can predict
DFT Hamiltonians of magnetic twisted materials. Therefore, in the future
development of the DDHT, xDeepH models and corresponding predictions of
magnetic materials can be included. In our vision, the DDHT has the potential
to cover all stable 2D materials in the future and serve as a handy tool for
researchers in the emerging field of twistronics.

The predicted results in the DDHT are not limited to Hamiltonians and
electronic structures, the identifications of topology in twisted materials
will also be included in future development. This will be achieved by
determining the irreducible representations of valence eigenstates at high-
symmetry points in momentum space under various symmetry operations and
comparing them with the elementary band representations (EBRs) in the
topological quantum chemistry [34, 35, 36]. If these irreducible
representations of a material are not linear combinations of EBRs, it would be
identified as having non-trivial topology. Based on the predicted eigenstates
by the DDHT models, theoretically, the band topology of twisted materials of
any chemical constitutions and twist angles can be accurately and efficiently
determined; meanwhile, the berry curvature may also be calculated to explore
exotic properties regarding topology. Note that the classification of topology
by topological quantum chemistry holds effective in the case of the generally-
existed flat bands in twisted materials [37]. In this way, our DDHT may
promisingly find all possible topological twisted materials and provide
opportunities in exploring the currently scarce twist angle-dependent
topological transitions.

##  IV Summary

In summary, by using the deep-learning method targeting at DFT Hamiltonians
implemented in the DeepH-E3 code, we train the deep neural network models for
124(5) types of homo(hetero)bilayer 2D materials and use them to predict DFT
Hamiltonians, as well as electronic structures, of their twisted structures at
DFT accuracy. Based on these well-trained models capable of doing accurate
predictions, we construct a deep-learning model database for Hamiltonian
prediction of twisted materials named DDHT, which consists of DeepH-E3 NN
models and related files like predicted DFT Hamiltonian matrices, eigenvalues
and eigenvectors, band structures, density of states, and identified flat
bands. Further, by screening the database, we list ultra flat bands in small-
angle twisted bP and other 19 types of large-angle twisted materials. The
accuracy of our DDHT models and predictions is demonstrated with comparable
MAEs and energy bands to respective DFT results. Our DDHT covers twisted
structures of a significant number of 2D materials and holds solid
universality over a wide range of twist angles down to about 1∘, showing the
reliability and universality of DeepH method and corresponding high-throughput
NN parameter strategy, offering DFT-level knowledge prior to deep
investigations and serving as a handy tool for researchers in the field of
twistronics. Another crucial contribution of this work if to perform the new
paradigm of building a material database, which is a collection of useful
models applied to series of materials instead of single material and can be
seen as a extensible MoE model. In the future development, our DDHT can be
expanded to cover all stable 2D materials and incorporate magnetism and
topology, satisfying the growing general interest toward exotic properties in
twisted materials.

##  V Methods

###  V.1 Preparation of Training Datasets

The bilayer unit cell for each type of 2D material in our database is fully
relaxed using Viena ab initio Simulation Package (VASP) [28, 38] to give the
most optimal lattice constants, stacking orders, and interlayer spacings. The
interactions between ions and electrons and the exchange correlations are
treated by the projected-augmented wave method [39, 40] and the Perdew-Burke-
Ernzerhof functional [41]. The interaction is corrected by DFT-D3 method [42,
43]. The convergence criteria for energy and force are $1\times 10^{-7}$ eV
and $1\times 10^{-3}$ eV/Å. The vacuum space between two periodic images along
the direction normal to the 2D plane is set to be more than 12 Å. The k-point
mesh is 11$\times$11$\times$1 using the Monkhorst-Pack scheme [44].

The atomic structures in the training dataset for each 2D material in the DDHT
are prepared as a series of $3\times 3$ non-twisted bilayer supercells (made
from the optimized primitive bilayer unit cells) with typically 24$\times$24
uniformly shifted interlayer stackings along two in-plane lattice vectors
followed by random perturbations to atomic positions (up to 0.1 Å) along three
respective Cartesian axes. In total, the training dataset for each material
has typically 576 perturbed and shifted $3\times 3$ supercell structures.

The Hamiltonian for each shifted and perturbed atomic structure in the
training datasets is calculated by first-principles density functional theory
(DFT) simulation code OpenMX [45, 46] using pseudo atomic orbital basis, norm-
conserving pseudo potential [47], and PBE exchange functional [41]. The energy
convergence criteria, cut-off energy, and $\Gamma$-centered k-point mesh are
set as $4\times 10^{-8}$ Hartree, 300 Ry, and no less than $3\times 3\times 1$
in the self-consistent calculations, respectively. The DFT Hamiltonians and
band structures used for evaluating the DDHT accuracy are calculated by OpenMX
code using similar settings. The training datasets, including the DFT
Hamiltonians and the corresponding structures, are processed prior to the
training. The DFT Hamiltonians are decomposed into small matrix blocks
$\mathcal{H}_{ij}$ that corresponds to local atomic structures
$\mathcal{R}_{ij}$ within the cut-off radius $R_{c}$, and these
$\mathcal{R}_{ij}$ are represented by crystal graphs where the on-site atoms
and the distances to their neighbors within $R_{c}$ are treated by vertex and
edge features.

###  V.2 Training and Prediction

The deep neural network models in the DDHT are trained by the DeepH-E3 code
[21, 22] with Euclidean symmetry of local coordinate transformation satisfied
by an E(3)-equivariant framework [48]. Details of the DeepH-E3 code used in
the DDHT are summarized as follows. The representation for the initial and the
intermediate vertex (edge) features is 64$\times$0e and 64$\times$0e +
32$\times$1o + 16$\times$2e + 8$\times$3o + 8$\times$4e, respectively, where
the first and the second number (separated by ”$\times$”) in each term
(connected by ”$+$”) stands respectively for the number of vectors and the
carried representation of spherical harmonics $l$, while the last letter e (o)
denotes even (odd) parity. Spherical harmonics with $l=0$ to 4 and $l=0$ to 5
are used for materials without and with spin-orbit coupling respectively.
During the model training, we employ an optimized learning scheme where the
training, validation, and testing set take respectively 80%, 15%, and 5% of
the total dataset, the batch size takes the value of 2, and the learning rate
starts with 0.002 and reduces by half when validation loss plateau (in this
case, decreasing less than 5% after every 120 epochs). The training takes at
least 1,000 epochs and is terminated upon the validation loss drops below
$1\times 10^{-7}$ for its averaged value or $1\times 10^{-5}$ for its maximum
value, ensuring the high accuracy and generalizability of DeepH-E3 models.
Similar loss values are given in evaluation of these models in testing sets
upon the stopping of model training.

The Hamiltonian matrices are predicted using only atomic structures of twisted
bilayer materials, corresponding overlap matrices, and trained DeepH-E3 models
as input. Based on the atomic structures and local-atomic basis, the overlap
matrices are directly obtained by a modified OpenMX code without conducting
DFT calculations. The $H_{\text{DeepH-E3}}$ matrices are sparse due to the
locality of local-atomic basis and the eigenvalues and eigenvectors are
calculated by efficient matrix diagonalization using the Pardiso package [49,
50] of Julia code.

##  VI Data Availability

The DeepH-E3 models of all 124 homobilayer and 5 heterobilayer materials and
their twisted structures of large $\theta$s will be available at the DDHT
website.

##  VII Acknowledgments

This work was supported by the Basic Science Center Project of NSFC (grant no.
2388201), the National Natural Science Foundation of China (grant no.
12334003), the National Science Fund for Distinguished Young Scholars (grant
no. 12025405) and the Ministry of Science and Technology of China (grant no.
2023YFA1406400). R.X. acknowledges the founding support from the Fellowship of
China Postdoctoral Science Foundation (Grant No. 2021TQ0187). to be updated

##  VIII Author Contributions

T.B. and R.X. contributed equally to this project.

##  IX Competing Interests

The authors declare no competing of interests.

## References

  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018b).
  * Sun  _et al._ [2024] X. Sun, M. Suriyage, A. R. Khan, M. Gao, J. Zhao, B. Liu, M. M. Hasan, S. Rahman, R.-s. Chen, P. K. Lam, _et al._ , Twisted van der waals quantum materials: Fundamentals, tunability, and applications, Chemical Reviews (2024).
  * Song  _et al._ [2021] T. Song, Q.-C. Sun, E. Anderson, C. Wang, J. Qian, T. Taniguchi, K. Watanabe, M. A. McGuire, R. Stöhr, D. Xiao, T. Cao, J. Wrachtrup, and X. Xu, Direct visualization of magnetic domains and moiré magnetism in twisted 2d magnets, Science 374, 1140 (2021).
  * Xu  _et al._ [2021] Y. Xu, A. Ray, Y.-T. Shao, S. Jiang, K. Lee, D. Weber, J. E. Goldberger, K. Watanabe, T. Taniguchi, D. A. Muller, K. F. Mak, and J. Shan, Coexisting ferromagnetic–antiferromagnetic state in twisted bilayer CrI3, Nature Nanotechnology 17, 143 (2021).
  * Stern  _et al._ [2021] M. V. Stern, Y. Waschitz, W. Cao, I. Nevo, K. Watanabe, T. Taniguchi, E. Sela, M. Urbakh, O. Hod, and M. B. Shalom, Interfacial ferroelectricity by van der waals sliding, Science 372, 1462 (2021).
  * Yasuda  _et al._ [2021] K. Yasuda, X. Wang, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Stacking-engineered ferroelectricity in bilayer boron nitride, Science 372, 1458 (2021).
  * Devakul  _et al._ [2021] T. Devakul, V. Crépel, Y. Zhang, and L. Fu, Magic in twisted transition metal dichalcogenide bilayers, Nature communications 12, 6730 (2021).
  * Jin  _et al._ [2019] C. Jin, E. C. Regan, A. Yan, M. I. B. Utama, D. Wang, S. Zhao, Y. Qin, S. Yang, Z. Zheng, S. Shi, K. Watanabe, T. Taniguchi, S. Tongay, A. Zettl, and F. Wang, Observation of moiré excitons in WSe2/WS2 heterostructure superlattices, Nature 567, 76 (2019).
  * Tran  _et al._ [2019] K. Tran, G. Moody, F. Wu, X. Lu, J. Choi, K. Kim, A. Rai, D. A. Sanchez, J. Quan, A. Singh, J. Embley, A. Zepeda, M. Campbell, T. Autry, T. Taniguchi, K. Watanabe, N. Lu, S. K. Banerjee, K. L. Silverman, S. Kim, E. Tutuc, L. Yang, A. H. MacDonald, and X. Li, Evidence for moiré excitons in van der waals heterostructures, Nature 567, 71 (2019).
  * Alexeev  _et al._ [2019] E. M. Alexeev, D. A. Ruiz-Tijerina, M. Danovich, M. J. Hamer, D. J. Terry, P. K. Nayak, S. Ahn, S. Pak, J. Lee, J. I. Sohn, M. R. Molas, M. Koperski, K. Watanabe, T. Taniguchi, K. S. Novoselov, R. V. Gorbachev, H. S. Shin, V. I. Fal’ko, and A. I. Tartakovskii, Resonantly hybridized excitons in moiré superlattices in van der waals heterostructures, Nature 567, 81 (2019).
  * Shimazaki  _et al._ [2020] Y. Shimazaki, I. Schwartz, K. Watanabe, T. Taniguchi, M. Kroner, and A. Imamoğlu, Strongly correlated electrons and hybrid excitons in a moiré heterostructure, Nature 580, 472 (2020).
  * Regan  _et al._ [2020] E. C. Regan, D. Wang, C. Jin, M. I. B. Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, M. Blei, J. D. Carlström, K. Watanabe, T. Taniguchi, S. Tongay, M. Crommie, A. Zettl, and F. Wang, Mott and generalized wigner crystal states in WSe2/WS2 moiré superlattices, Nature 579, 359 (2020).
  * Xian  _et al._ [2019] L. Xian, D. M. Kennes, N. Tancogne-Dejean, M. Altarelli, and A. Rubio, Multiflat bands and strong correlations in twisted bilayer boron nitride: Doping-induced correlated insulator and superconductor, Nano Letters 19, 4934 (2019).
  * Kennes  _et al._ [2020] D. M. Kennes, L. Xian, M. Claassen, and A. Rubio, One-dimensional flat bands in twisted bilayer germanium selenide, Nature Communications 11, 1124 (2020).
  * Kennes  _et al._ [2021] D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. N. Basov, A. N. Pasupathy, and A. Rubio, Moiré heterostructures as a condensed-matter quantum simulator, Nature Physics 17, 155 (2021).
  * Carr  _et al._ [2020] S. Carr, S. Fang, and E. Kaxiras, Electronic-structure methods for twisted moiré layers, Nature Reviews Materials 5, 748 (2020).
  * Senior  _et al._ [2020] A. W. Senior, R. Evans, J. Jumper, J. Kirkpatrick, L. Sifre, T. Green, C. Qin, A. Žídek, A. W. R. Nelson, A. Bridgland, H. Penedones, S. Petersen, K. Simonyan, S. Crossan, P. Kohli, D. T. Jones, D. Silver, K. Kavukcuoglu, and D. Hassabis, Improved protein structure prediction using potentials from deep learning, Nature 577, 706 (2020).
  * Zhang  _et al._ [2018] L. Zhang, J. Han, H. Wang, R. Car, and W. E, Deep potential molecular dynamics: A scalable model with the accuracy of quantum mechanics, Physical Review Letters 120, 143001 (2018).
  * Behler and Parrinello [2007] J. Behler and M. Parrinello, Generalized neural-network representation of high-dimensional potential-energy surfaces, Physical review letters 98, 146401 (2007).
  * Li  _et al._ [2022] H. Li, Z. Wang, N. Zou, M. Ye, R. Xu, X. Gong, W. Duan, and Y. Xu, Deep-learning density functional theory Hamiltonian for efficient ab initio electronic-structure calculation, Nature Computational Science 2, 367 (2022).
  * Gong  _et al._ [2023] X. Gong, H. Li, N. Zou, R. Xu, W. Duan, and Y. Xu, General framework for E(3)-equivariant neural network representation of density functional theory Hamiltonian, Nature Communications 14, 2848 (2023).
  * Jain  _et al._ [2013] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. A. Persson, Commentary: The materials project: A materials genome approach to accelerating materials innovation, APL Materials 1, 011002 (2013).
  * Hellenbrandt [2004] M. Hellenbrandt, The inorganic crystal structure database (icsd)—present and future, Crystallography Reviews 10, 17 (2004).
  * Haastrup  _et al._ [2018] S. Haastrup, M. Strange, M. Pandey, T. Deilmann, P. S. Schmidt, N. F. Hinsche, M. N. Gjerding, D. Torelli, P. M. Larsen, A. C. Riis-Jensen, J. Gath, K. W. Jacobsen, J. J. Mortensen, T. Olsen, and K. S. Thygesen, The Computational 2D Materials Database: High-throughput modeling and discovery of atomically thin crystals, 2D Materials 5, 042002 (2018).
  * Gjerding  _et al._ [2021] M. N. Gjerding, A. Taghizadeh, A. Rasmussen, S. Ali, F. Bertoldo, T. Deilmann, N. R. Knøsgaard, M. Kruse, A. H. Larsen, S. Manti, T. G. Pedersen, U. Petralanda, T. Skovhus, M. K. Svendsen, J. J. Mortensen, T. Olsen, and K. S. Thygesen, Recent progress of the Computational 2D Materials Database (C2DB), 2D Materials 8, 044002 (2021).
  * Jacobs  _et al._ [1991] R. A. Jacobs, M. I. Jordan, S. J. Nowlan, and G. E. Hinton, Adaptive mixtures of local experts, Neural computation 3, 79 (1991).
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational Materials Science 6, 15 (1996a).
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Physical review B 54, 11169 (1996b).
  * Xie and Grossman [2018] T. Xie and J. C. Grossman, Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties, Physical Review Letters 120, 145301 (2018).
  * Manzeli  _et al._ [2017] S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, 2d transition metal dichalcogenides, Nature Reviews Materials 2, 17033 (2017).
  * Deringer  _et al._ [2020] V. L. Deringer, M. A. Caro, and G. Csányi, A general-purpose machine-learning force field for bulk and nanostructured phosphorus, Nature Communications 11, 5461 (2020).
  * Li  _et al._ [2023] H. Li, Z. Tang, X. Gong, N. Zou, W. Duan, and Y. Xu, Deep-learning electronic-structure calculation of magnetic superstructures, Nature Computational Science 3, 321 (2023).
  * Bradlyn  _et al._ [2017] B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, and B. A. Bernevig, Topological quantum chemistry, Nature 547, 298 (2017).
  * Vergniory  _et al._ [2019] M. G. Vergniory, L. Elcoro, C. Felser, N. Regnault, B. A. Bernevig, and Z. Wang, A complete catalogue of high-quality topological materials, Nature 566, 480 (2019).
  * Vergniory  _et al._ [2022] M. G. Vergniory, B. J. Wieder, L. Elcoro, S. S. P. Parkin, C. Felser, B. A. Bernevig, and N. Regnault, All topological bands of all nonmagnetic stoichiometric materials, Science 376, eabg9094 (2022).
  * Călugăru  _et al._ [2021] D. Călugăru, A. Chew, L. Elcoro, Y. Xu, N. Regnault, Z.-D. Song, and B. A. Bernevig, General construction and topological classification of crystalline flat bands, Nature Physics 18, 185 (2021).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Physical Review B 54, 11169 (1996).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Physical Review B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Physical Review B 59, 1758 (1999).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Physical Review Letters 77, 3865 (1996).
  * Grimme  _et al._ [2011] S. Grimme, S. Ehrlich, and L. Goerigk, Effect of the damping function in dispersion corrected density functional theory, Journal of Computational Chemistry 32, 1456 (2011).
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (dft-d) for the 94 elements h-pu, The Journal of Chemical Physics 132, 154104 (2010).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Special points for brillouin-zone integrations, Physical Review B 13, 5188 (1976).
  * Ozaki [2003] T. Ozaki, Variationally optimized atomic orbitals for large-scale electronic structures, Phys. Rev. B 67, 155108 (2003).
  * Ozaki and Kino [2004] T. Ozaki and H. Kino, Numerical atomic basis orbitals from h to kr, Phys. Rev. B 69, 195113 (2004).
  * Morrison  _et al._ [1993] I. Morrison, D. M. Bylander, and L. Kleinman, Nonlocal hermitian norm-conserving vanderbilt pseudopotential, Phys. Rev. B 47, 6728 (1993).
  * Geiger and Smidt [2022] M. Geiger and T. Smidt, e3nn: Euclidean neural networks, arXiv preprint arXiv:2207.09453 10.48550/arXiv.2207.09453 (2022).
  * Alappat  _et al._ [2020] C. Alappat, A. Basermann, A. R. Bishop, H. Fehske, G. Hager, O. Schenk, J. Thies, and G. Wellein, A recursive algebraic coloring technique for hardware-efficient symmetric sparse matrix-vector multiplication, ACM Trans. Parallel Comput. 7, 37 (2020).
  * Bollhöfer  _et al._ [2019] M. Bollhöfer, A. Eftekhari, S. Scheidegger, and O. Schenk, Large-scale sparse inverse covariance matrix estimation, SIAM Journal on Scientific Computing 41, A380 (2019).


# Influence of imperfections on tunneling rate in $\delta$-layer junctions

Juan P. Mendez<EMAIL_ADDRESS>Shashank Misra  Denis Mamaluy
<EMAIL_ADDRESS>Sandia National Laboratories, Albuquerque, New Mexico,
87123

###### Abstract

The atomically precise placement of dopants in semiconductors using scanning
tunneling microscopes has been used to create planar dopant-based devices,
enabling the exploration of novel classical or quantum computing concepts,
which often require precise control over tunneling rates in their operation.
While the geometry of the dopants can be defined to sub-nanometer precision,
imperfections can still play a significant role in determining the tunneling
rates. Here, we investigate the influence of different imperfections in
phosphorous $\delta$-layer tunnel junctions in silicon: variations of
$\delta$-layer thickness and tunnel gap width, interface roughness, and
charged impurities. It is found that while most of the imperfections
moderately affect the tunneling rate, a single charged impurity in the tunnel
gap can alter the tunneling rate by more than an order of magnitude, even for
relatively large tunnel gaps. Moreover, it is also revealed that the tunneling
rate strongly depends on the electrical charge sign of the impurity.

##  I Introduction

Atomic precision advanced manufacturing (APAM) can be used to create 2D doped
regions (known as $\delta$-layers) in silicon that simultaneously have single-
atom precision [1, 2, 3, 4, 5, 6] and very high conductivity [7, 8, 9, 10, 11,
12]. APAM has application for exploring basic principles of novel electronic
devices, including nano-scale diodes and transistors for classical computing
and sensing systems [13, 14, 11, 15] (see e.g. Fig. 1 a and b). Primarily,
however, this technology is used to explore dopant-based qubits in silicon,
with recent advancements in understand exchange-based 2-qubit operations [16],
the limits to qubit fidelity from environmental noise [17], the advantages to
leveraging the number of dopants as a degree of freedom [18, 19], and the
exploration of many body [20] and topological [21] effects in dopant chains
(see e.g. Fig. 1 c). In principle, atomically precise fabrication imbues the
kind of control required by these applications, which have a high sensitivity
to tunnel rates.

However, in reality, APAM involves tradeoffs between a number of defect
mechanisms whose impact on tunnel rates have not been systematically studied,
and this work pursues. A general processing tradeoff exists where the point
defect density can be reduced by increasing the various processing
temperatures, at the expense of worse dopant placement uncertainty from
activating dopant diffusion [10]. More specifically, there is a well-known
intrinsic stochasticity from the underlying chemistry resulting in a dopant
placement uncertainty of $\pm 0.3$ nm [1, 4]. Moreover, after dopant
incorporation, $\delta$-layer devices must be capped with silicon at moderate
temperature to protect them, but adatom-mediated diffusion can lead to a loss
of out-of-plane sharpness that is on the order of 1 nm [22]. In contrast, the
low temperature capping layer growth also leads to charged point defects at a
density of $\sim 1$ defect in a ($10$ nm)3 volume [23]. Determining which of
these disorder mechanisms is most likely to create large variations in
tunneling rates in a simpler tunnel junction geometry will help inform how to
navigate these tradeoffs, and lays the groundwork to analyze the more
complicated case of qubits in the future.

Figure 1: Examples of atomic-precision advanced manufacturing nanodevices. a
$\delta$-layer tunnel junction field-effect transistor (FET) introduced by
[15]. b A $\delta$-layer tunnel junction is used for detecting electrons in
the quantum dot (QD) in [14]. c Two-qubit donor spin device from [16, 17].

Electron tunneling is understood in different terms at different length scales
– ranging from semi-empirical descriptions of resonant tunnel diodes in terms
of effective barrier heights [24] to atomistic hopping through single
molecules in break junctions [25]. Challenges arise in problems where the
macro-scale description of tunneling in terms of band structure cannot account
for atomic-scale details by simply renormalizing parameters. A direct quantum-
mechanical way to investigate tunneling requires an open-system charge
transport treatment [26]. In this work we will employ an efficient, charge
self-consistent, quantum transport implementation of open-system Non-
Equilibrium Green Function (NEGF) formalism, known as the Contact Block
Reduction (CBR) method [27, 28, 29, 30, 31, 32, 33, 34]. We combine it with an
effective mass description of the electronic structure, shown to be in a very
good agreement with tight-binding models for system sizes down to 3 nm [35,
36], to assess the effect of imperfections on the tunneling rate for
phosphorous $\delta$-layer tunnel junctions in silicon (Si: P $\delta$-layer
tunnel junction). The considered imperfections in this work include variations
of the $\delta$-layer thickness, small variations of the tunnel junction gap
length, roughness in the edges of the $\delta$-layers, and the presence of
charged impurities in the intrinsic gap.

Figure 2: Si:P $\delta$-layer tunnel junction (TJ) devices. a Ideal device,
which consists of a semi-infinite source and drain, in contact with the
channel. The channel is composed of a lightly doped Si body and Si cap and a
very thin, highly P doped-layer with an intrinsic gap of length $L_{gap}$. b
Device with roughness in $\delta$-layer edges; The edge roughness is modeled
as periodic protrusions of size $d_{1}\times d_{2}\times t$ with a periodicity
of $2\times d_{2}$. c Device with presence of a charged impurity, either
p-type or n-type, in the center of the intrinsic tunnel gap; The charged
impurity is represented as a green sphere in the figure.

##  II Methodology

To explore the impact of these defect mechanisms, we adopt a structure of a
$\delta$-layer tunnel junction, which consist of two highly-conductive
$\delta$-layers separated by an intrinsic semiconductor gap. In the open-
system NEGF framework, the computational device consists of a semi-infinite
source and drain, in contact with the channel of length $L$, which is composed
of a lightly doped Si body and Si cap and two very thin, highly P-doped layers
(referred to as left and right $\delta$-layers) separated by intrinsic gap of
length $L_{gap}$, as shown in Fig. 2 a. The channel length is chosen to be
$L=30~{}\text{nm}+L_{gap}$ to avoid the boundary effect between the source and
drain contacts, the device height is $H=12$ nm and the device width is chosen
to be $W=15$ nm to avoid size quantization effects on the conductive
properties of $\delta$-layer systems [37]. We consider three different
$\delta$-layer thicknesses: $t=0.2$ nm to approximate the true mono-atomic
$\delta$-layer, an intermediate $t=1.0$ nm and the ”thick” $\delta$-layer of
$t=5.0$ nm. The sheet doping density of $\delta$-layer is $N_{D}=1.0\times
10^{14}$ cm-2 ($N_{D}^{(2D)}=t\times N_{D}^{(3D)}$) and the doping densities
in the Si cap and Si body are $N_{A}=5.0\times 10^{17}$ cm-3 in all
simulations. Furthermore, all simulations are carried out at the cryogenic
temperature of $4$K, for which we can neglect inelastic scatterings [7, 38].

The simulations in this work are conducted using the open-system charge self-
consistent Non-Equilibrium Green Function (NEGF) Keldysh formalism [39],
together with the Contact Block Reduction (CBR) method [27, 28, 29, 30, 31,
32, 33, 34] and the effective mass theory. The CBR method allows a very
efficient calculation of the density matrix, transmission function, etc. of an
arbitrarily shaped, multiterminal two- or three-dimensoinal open device within
the NEGF formalism and scales linearly $O(N)$ with the system size $N$. As
validation, in our previous works [40, 41], we showed for Si: P $\delta$-layer
systems a very good agreement with experimental electrical measurements [7,
42, 43, 44], proving a excellent reliability of this framework to investigate
$\delta$-layer systems.

Within this framework, we solve self-consistently the open-system effective
mass Schrödinger equation and the non-linear Poisson equation [27, 30, 33]. We
employ a single-band effective mass approximation with a valley degeneracy of
$d_{val}=6$. For the charge self-consistent solution of the non-linear Poisson
equation we use a combination of the predictor-corrector approach and Anderson
mixing scheme [31, 33]. First, the Schrödinger equation is solved in a
specially defined (see the generalized Neumann BC in Section C of
Supplementary Material [45]) closed-system basis taking into account the
Hartree potential $\phi^{H}(\bm{r}_{i})$ and the exchange and correlation
potential $\phi^{XC}(\bm{r}_{i})$. Second, the LDOS of the open system,
$\rho(\bm{r}_{i},E)$, and the electron density, $n(\bm{r}_{i})$, are computed
using the CBR method for each iteration. Then the electrostatic potential and
the carrier density are used to calculate the residuum $F$ of the Poisson
equation

$\big{|}\big{|}\bm{F}[\bm{\phi}^{H}(\bm{r}_{i})]\big{|}\big{|}=\big{|}\big{|}\bm{A}\bm{\phi}^{H}(\bm{r}_{i})-(\bm{n}(\bm{r}_{i})-\bm{N}_{D}(\bm{r}_{i})+\bm{N}_{A}(\bm{r}_{i}))\big{|}\big{|},$
(1)

where $\bm{A}$ is the matrix derived from the discretization of the Poisson
equation and $\bm{N}_{D}$ and $\bm{N}_{A}$ are the total donor and acceptor
doping densities arrays, respectively. If the residuum is larger than a
predetermined threshold $\epsilon$, the Hartree potential is updated using the
predictor-corrector method, together with the Anderson mixing scheme. Using
the updated Hartree potential and the corresponding carrier density, the
exchange-correlation is computed again for the next step, and an iteration of
Schrodinger-Poisson is repeated until the convergence is reached with
$\big{|}\big{|}\bm{F}[\bm{\phi}^{H}(\bm{r}_{i})]\big{|}\big{|}<\epsilon=10^{-6}$
eV. Further details of the methodology are included in the Supplementary
Material [45]. In our simulations we have utilized a 3D real-space model, with
a discretization size of 0.2 nm along all directions, thus with about $10^{6}$
real-space grid points, and around 3,000 energy points were used. We have also
employed the standard values of electron effective masses, $m_{l}=0.98\times
m_{e}$, $m_{t}=0.19\times m_{e}$, the dielectric constant of Silicon,
$\epsilon_{Si}=11.7$, and the cryogenic temperature of T$=4.0$ K in all our
simulations.

##  III Results and discussion

###  III.1 Conductivity of ideal tunnel junctions

Figure 3:  Characteristic tunneling current curves. Total current $I$
(logarithmic scale) vs. tunnel gap length $L_{gap}$ for different values of
$\delta$-layer thickness $t$ and applied voltages ($1$ mV and $100$ mV). Black
dotted lines represent least-square fits to the exponential trend.
$N_{D}=1.0\times 10^{14}$ cm-2 and $N_{A}=5.0\times 10^{17}$ cm-3.  Figure 4:
Two conductivity regimes in $\delta$-layer tunnel junctions. Total current vs
voltage (blue curve, linear scale) and the corresponding differential
resistance $dV/dI$ (red curve, logarithmic scale) are shown for $L_{gap}=10$
nm, $N_{D}=1.0\times 10^{14}$ cm-2, $N_{A}=5.0\times 10^{17}$ cm-3 and $t=1$
nm. Figure 5: Local Density of States for $\delta$-layer tunnel junctions. The
$LDOS(E,x)$ for a tunnel junction of $L_{gap}=10$ nm is shown in a and b when
a voltage of 1 mV and 100 mV is applied to the drain contact, respectively.
The Fermi levels indicated in the figures correspond to the Fermi levels of
the source and drain contacts. In a and b, the corresponding effective 1D
potentials are also shown, calculated by integrating over the (y,z)-plane the
actual charge self-consistent 3D potentials weighted with the electron
density. $N_{D}=1.0\times 10^{14}$ cm-2, $N_{A}=5.0\times 10^{17}$ cm-3, and
$t=1$ nm.

In our previous work for infinite-width ($W\to\infty$) $\delta$-layer sheets
[41], we demonstrated that the distribution of dopants along the confinement
z-axis for a fixed sheet doping density ($N_{D}^{(2D)}$) significantly affects
the current. Conductivity decreases for sharper $\delta$-layer doping
profiles, which create strong transverse electric fields in their vicinity,
while it increases for broader $\delta$-layers doping profiles, which
conversely create weaker transverse electric fields. We also report in this
work that the same trend is observed for finite-width $\delta$-layer tunnel
junctions. In Fig. 3, the tunneling current $I$ vs the tunnel gap length
$L_{gap}$ for an ideal $\delta$-layer tunnel junction (see 2 a) is included
for different $\delta$-layer thicknesses and two voltages: 1mV and 100 mV. As
shown, the tunneling current decreases as the $\delta$-layer thickness
decreases for a fixed sheet doping density (i.e. the total charge density is
kept constant). We also find that the overall $I$ vs $L_{gap}$ trend is
practically exponential for all tunnel gap lengths $L_{gap}=0,...,12$ nm, i.e.
$\ln{I}\sim\ln{I_{L_{gap}=0}}-L_{gap}/B_{voltage}$, where $I_{L_{gap}=0}$ is
the current when $L_{gap}=0$, $L_{gap}$ is the tunnel gap length and
$B_{voltage}$ is a proportional constant related with the barrier height. As a
guide to the eyes, the exponential $I$ vs $L_{gap}$ trend is included in the
figure as black dotted lines for $t=1$ nm. However, a deviation from the
exponential trend can be noticed for large tunnel gaps $L_{gap}>7$ nm and an
applied voltage of $1$ mV (dashed lines in the figure). Conversely, for a
voltage of $100$ mV (continous lines in the figure), the deviation from the
ideal trend vanishes and the overall trend is exponential, thus shedding light
on two conductivity regimes for $\delta$-layer tunnel junctions: low- and
high-voltage regimes. As reported in [37], the deviation from the ideal trend
is the result of the quantized conduction band in both $\delta$-layers (left
and right) and a certain mismatch between them, which can only occurs for low
applied voltages.

Fig. 4 shows the I-V characteristic curve for an ideal tunnel junction of
length $L_{gap}=10$ nm and $t=$1 nm. We can discern two Ohmic behaviors,
elucidating again the existence of two distinct conductivity regimes
corresponding to the low and high voltage: the first one, between $0$ V and
$0.04$ V, with a resistance of approximately $50$ M$\Omega$; the second one,
above $0.08$ V, with a resistance of $2$ M$\Omega$. The resistance in the low-
voltage regime is one order of magnitude higher than in the high-voltage
regime. Between these two regimes, approximately between $0.04$ V and $0.08$
V, there is a transition region over which the resistance is reduced. The I-V
characteristics recently measured by [15] agree with the existence of the two
conductivity regimes.

To get a better understanding of the two conductivity regimes (see Fig. 4) and
the strong influence of the quantized conduction band on the tunneling current
for low voltages (see the oscillations for $L_{gap}>7$ nm in Fig. 3), we
examine the local density of states (LDOS), which represent the available
states that can be occupied by the free electrons in space-energy dimension.
For very low temperatures the states below Fermi level are occupied, whereas
the states above Fermi level are unoccupied. Fig. 5 shows the LDOS along
x-direction for a tunnel junction of $L_{gap}=10$ nm when a voltage of 1 mV in
a and 100 mV in b is applied to the drain contact. Additionally, the
corresponding effective 1D potential is included in the figure, exhibiting a
tunnel barrier height of approximately 55 meV for the equilibrium case, which
is in excellent agreement with experimental measurements in [15]. Firstly, as
shown in Fig. 5, the low-energy LDOS are strongly quantized, highlighted with
dashed lines in the figure. On the contrary, for high energies, the LDOS are
practically continuous in space-energy, thus these states are not quantized.
When a voltage is applied to the drain contact, the Fermi level corresponding
to the drain contact is reduced, resulting in lowering the energies of all
states in the right side as well. As a result, new unoccupied states in the
right $\delta$-layer will be available to be occupied by the tunneling
electrons coming from the left $\delta$-layer. When a low drain voltage is
applied, $<45-50~{}\text{mV}$, only the unoccupied quantized states in the
right $\delta$-layer will play a role in the tunneling process. If the
occupied quasi-discrete states near the Fermi level in the left side align
with the unoccupied quasi-discrete states in the right side, it will result in
a considerable increase of the tunneling current as shown in Fig. 3 for
$L_{gap}=10$ nm and $t=1$ nm. Conversely, if the overlap is minimum, as happen
for $L_{gap}=11$ nm and $t=1$ nm, the tunneling current will be reduced. For
low biases, this alternating mismatch can only exist for sufficiently large
tunnel gaps, $L_{gap}>7$ nm, because the coupling of the left and right
$\delta$-layer wave-functions for narrow tunnel gaps ($L_{gap}<7$ nm)
equalizes the electron states on both sides, increasing the overlap and thus
eliminating the mismatch. When a high bias is applied, as in Fig. 5 b, it
makes the continuous unoccupied high-energy states in the right side available
for tunneling from the left side, thus diminishing the influence of the
conduction band quantization on the current, as can be seen in $I$ vs
$L_{gap}$ plots in Fig. 3 for $100$ mV.

In the following we will evaluate the effect of diverse imperfections in
$\delta$-layer tunnel junctions on the tunneling current. We therefore remark
the importance of the evaluation of the tunneling rate for both conductivity
regimes: low- and high voltage regimes. For low voltages, strong quantization
effects on the tunneling current are expected, especially for large tunnel
gaps, because of the quantized low-energy conduction band, thus it will be
reflected in noisy characteristic curves. Conversely, in the high-bias regime,
no significant influence of the low-energy conduction band quantization is
expected on the tunneling current.

###  III.2 Effects of $\delta$-layer thickness deviations from mono-atomic
layer

Figure 6:  Effect of $\delta$-layer thickness variation. Current ratio,
$I_{t}/I_{t=0.2~{}\text{nm}}$, vs. tunnel gap length for different deviations
of the $\delta$-layer thickness from the ”ideal” mono-atomic layer.

We first investigate the effects of $\delta$-layer thickness variation on the
tunneling rate. Fig. 6 shows the current ratio between two different
$\delta$-layer thicknesses ($t=1,~{}5$ nm) with respect to the ideal mono-
atomic layer, which is approximately 0.2 nm, in dashed lines for a low-bias of
$1$ mV and in continuous lines for a high-bias of $100$ mV. Our results
suggest that the tunneling rate increases approximately up to two times for a
broadening of the $\delta$-layer of $1$ nm. Interestingly, the effect is
roughly constant for the whole considered tunnel gap $L_{gap}$ range. However,
for even wider $\delta$-layers, the tunneling rate further increases, between
two and seven times for a effective thickness of $5$ nm, developing now a
stronger dependence on the tunnel gap length. Furthermore, for large tunnel
gaps in the low-bias regime ($1$ mV), the strong non-monotonic dependence of
the current ratio on the gap length is the result of the space-energy
quantization of the conduction band and the mismatch in the overlapping
between the states in the left and right $\delta$-layers, discussed in Sect.
III.1. When a higher voltage is applied (e.g. for $100$ mV), it results in a
lower current ratio, and the quantization effects diminish as we discussed in
Sect III.1.

###  III.3 Effects of tunnel gap length variations and interface roughness

Figure 7:  Effect of tunnel gap length variation. Current ratio,
$I_{L_{gap}+\delta L_{gap}}/I_{L_{gap}}$ vs. tunnel gap length, $L_{gap}$, for
distinct applied voltages. $t=1.0$ nm, $N_{D}=1.0\times 10^{14}$ cm-2 and
$N_{A}=5.0\times 10^{17}$ cm-3.  Figure 8:  Effect of interface roughness.
Current ratio, $I_{ideal}/I_{roughness}$, vs. tunnel gap length, $L_{gap}$,
for an applied bias of $1.0$ mV. For each tunnel gap length, different
roughness sizes have been considered, ranging $d_{1}=0.8-2.0$ nm and
$d_{2}=0.6-3.4$ nm (see Fig. 2 b). The dots represent the average and the bars
represent the dispersion of the values. $t=1.0$ nm, $N_{D}=1.0\times 10^{14}$
cm-2 and $N_{A}=5.0\times 10^{17}$ cm-3.

Next we will assess how the variation of the tunnel gap length and the
roughness of the $\delta$-layer might affect the tunneling rate. As a first
approximation, the interface roughness can be modeled as a uniformed increase
or reduction of the average gap length $\langle L_{gap}\rangle$. We evaluate
in Fig. 7 the change in the tunneling current for small uniform variations,
such as $L_{gap}+\delta L_{gap}$ with $\delta L_{gap}=\pm 0.2$ nm. This small
perturbation is of the order of the stochasticity of APAM chemistry. Our
simulations suggest that a small variations of the tunnel gap length (of the
order of $\pm 0.2$ nm) can lead to a current change of around 20$\%$ with
respect to the ideal (or designed) length for the whole studied range. A
reduction of the tunnel gap length evidently results in an increase of the
tunneling rate, while an increase of the effective gap length leads to a
decrease of the tunneling rate by a similar magnitude. For large tunnel gaps
in the low-bias regime (see dashed lines for $1$ mV), our results similarly
present the expected non-monotonic behaviour due to the quantization effect,
but this effect vanishes when higher voltage is applied (see continuous lines
for $100$ mV).

A second-order analysis of the edge roughness can be performed assuming that
the average $\langle L_{gap}\rangle$ value does not change due to the
roughness. In this case, we model the roughness as periodic protrusion of size
$d_{1}\times d_{2}\times t$, with a periodicity of $2\times d_{2}$, as shown
in Fig. 2 b, instead of the uniformed variation evaluated above. To maintain
$\langle L_{gap}\rangle$ constant we have only considered in-phase roughness,
i.e. the protrusion of the left $\delta$-layer is exactly a mirror of the
right side. The analysis of out-phase roughness is the scope of this work and
it will be further investigated elsewhere. In our analysis, we have considered
different roughness sizes, varying the parameter $d_{1}$ from $0.8$ nm to
$2.0$ nm, and the parameter $d_{2}$ from $0.6$ nm to $3.4$ nm. Fig. 8
includes, for all evaluated roughness sizes, the current ratio between the
”non-ideal” tunnel junction device (with edge roughness) and the ”ideal”
device (without any roughness) for a voltage of $1$ mV. The simulations
indicate that edge roughness reduces the tunneling rate between $6\%$ and
$20\%$ for almost the whole considered gap range, predicting very similar
magnitude of the tunneling rate change as the uniformed variation of the
tunnel gap length. One can also notice that for large tunnel gaps, especially
for $L_{gap}=9-10$ nm, the tunneling rate is even further reduced up to 2.2
times due to the mentioned quantization of the conduction bands.

###  III.4 Effects of charged impurities in the tunnel gap

Figure 9: Local Density of States for $\delta$-layer tunnel junctions with
impurities. The $LDOS(E,x)$ for a tunnel junction of $L_{gap}=10$ nm is shown:
for a n-type impurity when a voltage of 1 mV and 100 mV is applied to the
drain contact in a and b, respectively; for a p-type impurity when a voltage
of 1 mV and 100 mV is applied in c and d, respectively. In all figures, the
corresponding effective 1D potentials are also shown in red color, calculated
by integrating over the (y,z)-plane the actual charge self-consistent 3D
potentials weighted with the electron density. The effective 1D potential for
the ideal tunnel junction is also included in orange dashed lines for
comparison purpose. $L_{gap}=10$ nm, $t=1.0$ nm, $N_{D}=1.0\times 10^{14}$
cm-2, and $N_{A}=5.0\times 10^{17}$ cm-3. Figure 10:  Localized states created
or depleted by a charged impurity. Shown is the $LDOS(E,x)$ difference for a
tunnel junction of $L_{gap}=10$ nm between: n-type impurity in the middle of
the tunnel gap and without impurity in a and b for an applied voltage of 1 mV
and 100 mV, respectively; p-type impurity and without impurity for an applied
voltage of 1mV and 100 mV in c and d, respectively. The corresponding
effective 1D potentials are also shown in red color, together with the
effective 1D potential for the ideal case in orange color. $N_{D}=1.0\times
10^{14}$ cm-2, $N_{A}=5.0\times 10^{17}$ cm-3, and $t=1$ nm. Figure 11:
Effect of charged impurities. Current ratio, $I_{non-ideal}/I_{ideal}$, vs.
tunnel gap length, $L_{gap}$, for tunnel junctions with a single n-type and
p-type impurities in the intrinsic gap: a t=0.2 nm; a t=1.0 nm; and c t=5.0
nm. The insets in a, b and c are a zoom of the result within the range between
3 nm and 8 nm. $N_{D}=1.0\times 10^{14}$ cm-2 and $N_{A}=5.0\times 10^{17}$
cm-3.

In the following, we evaluate the effects on the tunneling rate by the
presence of a single charged impurity in the tunnel gap, as shown in Fig. 2 c.
For the simulation set-up, we place in the center of the tunnel gap either a
n-type impurity (e.g. a phosphorus atom) or a p-type impurity (e.g. an
aluminium atom). While in this work we restrict our analysis to the center-gap
location, the influence of other locations may be also interesting to
investigate since the free electrons in $\delta$-layer systems form distinct
conducting layers perpendicular to the confinement direction, thus signaling a
highly non-homogeneous electron density distribution [41, 40].

Fig. 9 shows the LDOS of a $\delta$-layer tunnel junction with a single n-type
impurity for an applied voltage of 1 mV in a and for 100 mV in b, and with a
single p-type impurity for 1 mV in c and for 100 mV in d. In addition, Fig. 10
shows the LDOS difference for both applied voltages between the ideal tunnel
junction (i.e. the LDOS shown in Fig. 5) and the junction with the impurity
(i.e. the LDOS shown in Fig. 9), therefore it represents the localized states
created (in red color) or depleted (in blue color) by the impurity. As Fig. 9
a and b reveal, an n-type impurity in the middle of the tunnel gap creates
unoccupied states above the Fermi level within the tunnel gap, i.e. between
$x=15-25$ nm, as one might discern comparing panels a and b in Fig. 5 and Fig.
9. The new available states are clearly shown in red color in panels a and b
in Fig. 10. On the contrary, a p-type impurity depletes unoccupied available
states above the Fermi level as shown in Fig. 9 c and d within the tunnel gap
($x=15-25$ nm). The depletion of the states due a p-type impurity can be also
seen in Fig. 10 c and d in blue color since we are representing the LDOS
difference between the ideal tunnel junction and the one with the p-type
impurity. Interestingly, the presence of the impurity can be felt even far
from the impurity position as Fig. 10 also reveals for both impurity types.
The impurity creates a quantized perturbation in the unoccupied available
states (above the Fermi level) even far from the impurity location, shown as
ripples between $x=0-15$ nm and $x=15-40$ nm above the Fermi level in the
figure. However, the intensity of the perturbations vanish as we move away
from the impurity location.

In the semi-classical picture, the energy difference between the peak of the
effective electrostatic potential and the Fermi level is the energy barrier
which the electrons have to overcome to tunnel from one $\delta$-layer to the
other. Because of the presence of the impurity in the middle of a nano-scale
gap, the effective electrostatic potential is obviously affected. For an
n-type impurity in the middle of tunnel gap, the height of the electrostatic
barrier is reduced, as shown in Fig. 9 a, in which the height of the
electrostatic barrier without impurity (the orange dashed curve) is slightly
higher than the one with an n-type impurity (the red continuous curve). In
addition, the n-type impurity creates a dip in the electrostatic potential due
to the donor atom. For a p-type impurity, acceptor atom, the depletion of the
states above the Fermi level, in turn, increases the height of the energy
barrier in the electrostatic potential with respect to the ideal case, as
depicted in Fig. 9 c, where the height of the electrostatic potential for the
ideal tunnel junction (red continuous curve) is higher than the electrostatic
potential for the tunnel junction with the impurity (orange dashed curve).

Fig. 11 shows the current ratio versus the tunnel gap length $L_{gap}$ for
both impurity types and voltage regimes ($1$ mV and $100$ mV) in a, b and c
for a $\delta$-layer thicknesses of $t=0.2,~{}1.0,~{}5.0$ nm, respectively.
These simulations first reveal that a n-type impurity increases the tunneling
rate, whereas a p-type impurity decreases the tunneling rate. This result can
be explained in two different, but related ways. The first one, which
corresponds to a semi-classical picture, is by examining the electrostatic
potentials shown in Fig. 9 and discussed above: an n-type impurity decreases
the barrier height (i.e. the energy difference between the peak of the
electrostatic potential and the Fermi level), whereas a p-type impurity
increases the barrier height. Then the tunneling current increases or
decreases according to the change of the barrier height. The other way is to
examine the unoccupied and occupied states in the LDOS for the conduction
band. When we apply a positive voltage to the drain contact, the right Fermi
level and occupied and unoccupied states move down, falling below the Fermi
level of the source. As result, electrons injecting from the source can tunnel
into these available states. The presence of an n-type impurity creates
unoccupied states within the tunnel gap, just above the Fermi level in
equilibrium. Similarly, when a voltage is applied, these available states
might also fall within the source and drain Fermi energy levels and,
therefore, they became as intermediate states in which electrons can tunnel in
and out, reducing then the tunneling length (i.e. the width of the barrier in
the semiclassical picture) and, therefore, increasing the tunneling rate. On
the contrary, the presence of a p-type impurity creates a depletion of the
states within the tunnel gap, as shown panels c and d in Fig. 10 as a blue
cloud around the effective electrostatic potential, increasing then the
effective tunneling length and, therefore, reducing the tunneling current.

We can also observe from Fig. 11 two different behavior, corresponding to
narrow tunnel gaps $L_{gap}=3-7$ nm and large tunnel gaps $L_{gap}>7$ nm, in
which the deviation of the tunneling current from the ideal one behave
differently. For narrow tunnel gaps, the magnitude of the effect on the
tunneling rate is very similar for both impurity types and not very pronounced
(see the inset figures): up to $3.75$ times increase and $2.5$ times reduction
of the tunneling rate for a n-type and p-type impurity, respectively, when the
tunnel gap length is $7$ nm. Our simulations also exhibit that there is not
significant difference on tunneling rate change in both voltage regimes ($1$
mV for the low-voltage regime and $100$ mV for the high-voltage regime) for
narrow tunnel gaps: the change of the tunnel rate with respect to the ideal
tunneling current is only slightly higher for higher voltages. The deviation
in the magnitude of the effect on the tunneling current between both impurity
types starts approximately for $L_{gap}>7$ nm: the effect of a n-type impurity
becomes much more relevant than for a p-type impurity, especially for the
high-voltage regime ($100$ mV in the figure), in which it is a few times
higher than in the low-voltage regime ($1$ mV in the figure). Interestingly,
our simulations suggest that the tunneling current can be up to $60$ times
higher and $20$ times lower for an n-type and a p-type impurity, respectively,
for a tunnel gap length of $12$ nm. Additionally, in the low-voltage regime
for a p-type impurity, the simulations reveal that the maximum effect on the
tunneling current is reached for a tunnel gap length of approximately $10$ nm,
then it reaches a plateau regime in which the tunneling ratio roughly change
with the tunnel gap length. Finally, we note that the impurity is less felt by
thicker $\delta$-layers, as our result indicates by reducing the current ratio
as the thickness of the $\delta$-layer increases.

##  IV Conclusion

We have employed an open-system quantum transport analysis to investigate the
effect of diverse imperfections on the tunneling rate in Si: P $\delta$-layer
tunnel junctions. These imperfections span from geometry variations of the
$\delta$-layer thickness and junction gap length, to the presence of charged
impurities, either n-type or p-type, in the intrinsic gap. It is shown that
while most of the disorders moderately affect the tunneling rate, a single
charged impurity in the tunnel gap can alter the tunneling rate by more than
an order of magnitude. Contrary to predictions of semiclassical impurity
scattering (mobility-based) models, the electric sign of impurity plays an
important role in the tunneling rate: the change of current due to an n-type
impurity in Si: P $\delta$-layer systems is several times higher than for a
p-type impurity, especially for large tunnel gap lengths. Similarly, we can
extrapolate these findings to other systems such as Si: B $\delta$-layer
tunnel junctions, in which the influence of a p-type impurity in the intrinsic
gap, instead of an n-type, would result in a dramatic increase of the
tunneling current.

These results immediately suggest that the overall geometric fidelity of the
APAM device fabrication can be less important than mitigation of charged
impurities nearby the junction, which can lead to a strong change of tunneling
rates. APAM-based qubits in particular require having tunneling rates that are
tightly controlled between multiple closely-spaced objects. This includes
exchange coupling between a pair of donor-based qubits, initialization from a
tunnel-coupled single electron transistor, and spin-to-charge conversion from
the resonant tunneling of the single electron transistor to the leads in
Figure 1. Importantly, a change in tunneling rates from the fabricated
geometry due to impurities can be hard to compensate for using electrostatic
gates. The size of the gates (tens of nm) and the required spacing for them to
not leak to one another (tens of nm) is much larger than the length scale for
tunnel coupling (a few nm). In the APAM geometry, these gates are few in
number - enough to control electrostatics, but too few to control tunneling
rates. This leads to the conclusion that adopting APAM processing practices
that minimize charged defects is more important than those that preserve the
absolute device geometry.

Finally, the extreme sensitivity of $\delta$-layers tunnel junctions on the
tunneling current to the presence of charges in the vicinity of the tunnel gap
opens a great opportunity to use them for quantum FET-based sensors for
biological, chemical and radiation applications. The signal detection (either
due to radiation or chemicals or molecules) at the sensing area would be
strongly enhanced due to the conduction band quantization created by the
highly-confined $\delta$-layers. The sensing area can be placed above the
tunnel gap, replacing the traditional gate in a conventional geometry.
Contrary to traditional FET-based sensors, instead of needing to accumulate
enough charges at the sensing area to invert the full channel and detect the
signal [46], quantum FETs based on $\delta$-layers would allow to detect even
small signals that correspond to a single charge, thus significantly enhancing
the sensitivity with respect to traditional FET-based sensors.

## acknowledgement

The authors are thankful to FAIR DEAL’s team at Sandia National Laboratories
for the discussions during the project’s meetings. This work is funded under
Laboratory Directed Research and Development Grand Challenge (LDRD GC)
program, Project No. 213017, and under Laboratory Directed Research and
development (LDRD) program, Project No. 227155, at Sandia National
Laboratories. Sandia National Laboratories is a multimission laboratory
managed and operated by National Technology and Engineering Solutions of
Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for
the U.S. Department of Energy’s National Nuclear Security Administration under
contract DE-NA-0003525. This paper describes objective technical results and
analysis. Any subjective views or opinions that might be expressed in the
paper do not necessarily represent the views of the U.S. Department of Energy
or the United States Government.

## Author contributions

J.P. and D.M. performed the central calculations and analysis presented in
this work. S.M. initially guided the analysis of this work. The manuscript was
written by J.P. and D.M., and S.M contributed to write the introduction and
draw the conclusions of this manuscript.

## References

  * Wilson  _et al._ [2006] H. F. Wilson, O. Warschkow, N. A. Marks, N. J. Curson, S. R. Schofield, T. C. G. Reusch, M. W. Radny, P. V. Smith, D. R. McKenzie, and M. Y. Simmons, Phys. Rev. B 74, 195310 (2006).
  * Warschkow  _et al._ [2016] O. Warschkow, N. J. Curson, S. R. Schofield, N. A. Marks, H. F. Wilson, M. W. Radny, P. V. Smith, T. C. G. Reusch, D. R. McKenzie, and M. Y. Simmons, The Journal of Chemical Physics 144, 014705 (2016), https://doi.org/10.1063/1.4939124 .
  * Fuechsle  _et al._ [2012] M. Fuechsle, J. A. Miwa, S. Mahapatra, H. Ryu, S. Lee, O. Warschkow, L. C. L. Hollenberg, G. Klimeck, and M. Y. Simmons, Nature Nanotechnology 7, 242 (2012).
  * Ivie  _et al._ [2021] J. A. Ivie, Q. Campbell, J. C. Koepke, M. I. Brickson, P. A. Schultz, R. P. Muller, A. M. Mounce, D. R. Ward, M. S. Carroll, E. Bussmann, A. D. Baczewski, and S. Misra, Phys. Rev. Applied 16, 054037 (2021).
  * Wyrick  _et al._ [2022] J. Wyrick, X. Wang, P. Namboodiri, R. V. Kashid, F. Fei, J. Fox, and R. Silver, ACS Nano 16, 19114 (2022).
  * Campbell  _et al._ [2022] Q. Campbell, A. D. Baczewski, R. E. Butera, and S. Misra, AVS Quantum Science 4, 016801 (2022), https://doi.org/10.1116/5.0075467 .
  * Goh  _et al._ [2006] K. E. J. Goh, L. Oberbeck, M. Y. Simmons, A. R. Hamilton, and M. J. Butcher, Phys. Rev. B 73, 035401 (2006).
  * Weber  _et al._ [2012] B. Weber, S. Mahapatra, H. Ryu, S. Lee, A. Fuhrer, T. C. G. Reusch, D. L. Thompson, W. C. T. Lee, G. Klimeck, L. C. L. Hollenberg, and M. Y. Simmons, Science 335, 64 (2012), https://www.science.org/doi/pdf/10.1126/science.1214319 .
  * McKibbin  _et al._ [2013] S. R. McKibbin, G. Scappucci, W. Pok, and M. Y. Simmons, Nanotechnology 24, 045303 (2013).
  * Keizer  _et al._ [2015] J. G. Keizer, S. R. McKibbin, and M. Y. Simmons, ACS Nano 9, 7080 (2015), pMID: 26083628, https://doi.org/10.1021/acsnano.5b01638 .
  * Škereň  _et al._ [2020] T. Škereň, S. A. Koester, B. Douhard, C. Fleischmann, and A. Fuhrer, Nature Electronics 3, 524 (2020).
  * Dwyer  _et al._ [2021] K. J. Dwyer, S. Baek, A. Farzaneh, M. Dreyer, J. R. Williams, and R. E. Butera, ACS Applied Materials & Interfaces 13, 41275 (2021).
  * Mahapatra  _et al._ [2011] S. Mahapatra, H. Büch, and M. Y. Simmons, Nano Letters 11, 4376 (2011), pMID: 21919458, https://doi.org/10.1021/nl2025079 .
  * House  _et al._ [2014] M. G. House, E. Peretz, J. G. Keizer, S. J. Hile, and M. Y. Simmons, Applied Physics Letters 104, 113111 (2014), https://doi.org/10.1063/1.4869032 .
  * Donnelly  _et al._ [2021] M. B. Donnelly, J. G. Keizer, Y. Chung, and M. Y. Simmons, Nano Letters 21, 10092 (2021), pMID: 34797661, https://doi.org/10.1021/acs.nanolett.1c03879 .
  * He  _et al._ [2019] Y. He, S.K.Gorman, and D. Keith et al., Nature 571, 371 (2019).
  * Kranz  _et al._ [2020] L. Kranz, S. K. Gorman, B. Thorgrimsson, Y. He, D. Keith, J. G. Keizer, and M. Y. Simmons, Advanced Materials 32, 2003361 (2020), https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.202003361 .
  * Krauth  _et al._ [2022] F. Krauth, S. Gorman, Y. He, M. Jones, P. Macha, S. Kocsis, C. Chua, B. Voisin, S. Rogge, R. Rahman, Y. Chung, and M. Simmons, Phys. Rev. Appl. 17, 054006 (2022).
  * Fricke  _et al._ [2021] L. Fricke, S. J. Hile, L. Kranz, Y. Chung, Y. He, P. Pakkiam, M. G. House, J. G. Keizer, and M. Y. Simmons, Nature Communications 12, 3323 (2021).
  * Wang  _et al._ [2022] X. Wang, E. Khatami, F. Fei, J. Wyrick, P. Namboodiri, R. Kashid, A. F. Rigosi, G. Bryant, and R. Silver, Nature Communications 13, 6824 (2022).
  * Kiczynski  _et al._ [2022] M. Kiczynski, S. K. Gorman, H. Geng, M. B. Donnelly, Y. Chung, Y. He, J. G. Keizer, and M. Y. Simmons, Nature 606, 694 (2022).
  * Hagmann  _et al._ [2018] J. A. Hagmann, X. Wang, P. Namboodiri, J. Wyrick, R. Murray, M. D. Stewart, R. M. Silver, and C. A. Richter, Applied Physics Letters 112, 043102 (2018), https://doi.org/10.1063/1.4998712 .
  * Anderson  _et al._ [2020] E. M. Anderson, D. M. Campbell, L. N. Maurer, A. D. Baczewski, M. T. Marshall, T.-M. Lu, P. Lu, L. A. Tracy, S. W. Schmucker, D. R. Ward, and S. Misra, Journal of Physics: Materials 3, 035002 (2020).
  * Sze and Ng [2006] S. M. Sze and K. K. Ng,  _Physics of semiconductor devices_ (John Wiley & Sons, 2006).
  * Gehring  _et al._ [2019] P. Gehring, J. M. Thijssen, and H. S. J. van der Zant, Nature Reviews Physics 1, 381 (2019).
  * Ferry  _et al._ [2022] D. K. Ferry, J. Weinbub, M. Nedjalkov, and S. Selberherr, Semiconductor Science and Technology 37, 043001 (2022).
  * Mamaluy  _et al._ [2003] D. Mamaluy, M. Sabathil, and P. Vogl, J. Appl. Phys. 93, 4628 (2003), https://doi.org/10.1063/1.1560567 .
  * Mamaluy  _et al._ [2004] D. Mamaluy, A. Mannargudi, D. Vasileska, M. Sabathil, and P. Vogl, Semiconductor Science and Technology 19, S118 (2004).
  * Sabathil  _et al._ [2004] M. Sabathil, D. Mamaluy, and P. Vogl, Semiconductor Science and Technology 19, S137 (2004).
  * Mamaluy  _et al._ [2005] D. Mamaluy, D. Vasileska, M. Sabathil, T. Zibold, and P. Vogl, Phys. Rev. B 71, 245321 (2005).
  * Khan  _et al._ [2007] H. R. Khan, D. Mamaluy, and D. Vasileska, IEEE T. Electron Dev. 54, 784 (2007).
  * Khan  _et al._ [2008] H. R. Khan, D. Mamaluy, and D. Vasileska, Journal of Physics: Conference Series 107, 012007 (2008).
  * Gao  _et al._ [2014] X. Gao, D. Mamaluy, E. Nielsen, R. W. Young, A. Shirkhorshidian, M. P. Lilly, N. C. Bishop, M. S. Carroll, and R. P. Muller, J. Appl. Phys. 115, 133707 (2014).
  * Mendez  _et al._ [2021] J. P. Mendez, D. Mamaluy, X. Gao, and S. Misra, in  _2021 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD)_ (2021) pp. 210–214.
  * Wang  _et al._ [2005] J. Wang, A. Rahman, A. Ghosh, G. Klimeck, and M. Lundstrom, IEEE Transactions on Electron Devices 52, 1589 (2005).
  * Neophytou  _et al._ [2008] N. Neophytou, A. Paul, M. S. Lundstrom, and G. Klimeck, Journal of Computational Electronics 7, 363 (2008).
  * Mendez and Mamaluy [2022] J. P. Mendez and D. Mamaluy, Scientific Reports 10.1038/s41598-022-20105-x (2022).
  * Mazzola  _et al._ [2014] F. Mazzola, C. M. Polley, J. A. Miwa, M. Y. Simmons, and J. W. Wells, Applied Physics Letters 104, 173108 (2014), https://doi.org/10.1063/1.4874651 .
  * Keldysh [1965] L. V. Keldysh, Sov. Phys. J. Exp. Theor. Phys. 20, 1018 (1965).
  * Mendez  _et al._ [2020] J. P. Mendez, D. Mamaluy, X. Gao, E. M. Anderson, D. M. Campbell, J. A. Ivie, T.-M. Lu, S. W. Schmucker, and S. Misra, in  _2020 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD)_ (2020) pp. 181–184.
  * Mamaluy  _et al._ [2021] D. Mamaluy, J. P. Mendez, X. Gao, and S. Misra, Communications Physics 4, 205 (2021).
  * Goh and Simmons [2009] K. E. J. Goh and M. Y. Simmons, Appl. Phys. Lett. 95, 142104 (2009), https://doi.org/10.1063/1.3245313 .
  * Reusch  _et al._ [2008] T. C. G. Reusch, K. E. J. Goh, W. Pok, W.-C. N. Lo, S. R. McKibbin, and M. Y. Simmons, J. Appl. Phys. 104, 066104 (2008), https://doi.org/10.1063/1.2977750 .
  * McKibbin  _et al._ [2014] S. R. McKibbin, C. M. Polley, G. Scappucci, J. G. Keizer, and M. Y. Simmons, Appl. Phys. Lett. 104, 123502 (2014), https://doi.org/10.1063/1.4869111 .
  * [45] See Supplemental Material at [URL will be inserted by publisher] for further details of the methodology.
  * Shoorideh and Chui [2014] K. Shoorideh and C. O. Chui, Proceedings of the National Academy of Sciences 111, 5111 (2014), https://www.pnas.org/doi/pdf/10.1073/pnas.1315485111 .