Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras
(USC Thesis Other)
Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
Controlling Electronic Properties of Two-Dimensional Quantum Materials:
Simulation at the Nexus of the Classical and Quantum Computing Eras
by
Lindsay Elizabeth Bassman
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
PHYSICS
May 2020
Copyright 2020 Lindsay Elizabeth Bassman
ii
Acknowledgements
I am immensely grateful to my advisor, Professor Aiichiro Nakano, for his teaching,
guidance, patience, and above all, his mentorship. I did not truly understand what a mentor was
until I joined his group. Co-PIs of the CACS group, Professors Rajiv Kalia and Priya Vashista,
have also provided invaluable training in becoming a successful scientist. I would like to thank
my PhD committee members including Profs. Nakano, Kalia, and Vashishta, along with
Professor Moh El-Naggar and Professor Paulo Branicio for their thought-provoking questions
and comments throughout my qualifying exam and PhD defense.
I would like to thank my peers at USC, in particular, Zoe, Nicky, Aravind, Pankaj, and
Subodh, for their collaborative efforts, advice on navigating life in and after PhD program, and
friendship. I would also like to thank all my collaborators for their contributions to my published
work. I would especially like to thank my collaborators in Japan, Professor Shimojo, Hiroyuki,
and Shogo, for their hospitality during my visit to Japan, and the development of and assistance
with a large code base upon which a portion of my dissertation is based.
Finally, I am extremely thankful to my family for their unconditional support, constant
encouragement, and enduring curiosity and interest in my scientific endeavors. I am especially
grateful to my Mom for giving me the confidence to follow my scientific aspirations; to my Dad
for advising me to maintain a practical skillset; to my partner, Nate, who was there every day to
celebrate my successes and put my setbacks into perspective; and to my dog, Tiki, for her
company throughout it all.
iii
Table of Contents
Acknowledgements ........................................................................................................................ ii
List of Tables ................................................................................................................................... v
List of Figures ................................................................................................................................. vi
Abstract ............................................................................................................................................ x
Chapter 1: Introduction ............................................................................................................. 1
1.1 Non-Adiabatic Quantum Molecular Dynamics Simulations ....................................................... 2
1.2 Machine-Learning-Aided Simulations ......................................................................................... 4
1.3 Dynamic Simulations on Quantum Computers ............................................................................ 6
Chapter 2: Two-Dimensional Transition Metal Dichalcogenides ........................................... 10
Chapter 3: Non-Adiabatic Quantum Molecular Dynamics Simulations ................................. 15
3.1 Introduction ................................................................................................................................ 15
3.2 Methods ...................................................................................................................................... 16
3.2.1 Non-Adiabatic Quantum Molecular Dynamics .................................................................................... 16
3.2.2 Dynamic Simulations ........................................................................................................................... 20
3.2.3 Static Simulations ................................................................................................................................. 22
3.3 Results ........................................................................................................................................ 23
3.3.1 Electronic Structure of Monolayer ....................................................................................................... 23
3.3.2 Light-Driven Lattice Dynamics ........................................................................................................... 24
3.3.3 Evidence of Strong Electron-Phonon Coupling ................................................................................... 26
3.3.4 Dynamics of Monolayer in a Heat Bath ............................................................................................... 28
3.3.5 Simulation Size Effect .......................................................................................................................... 29
3.4 Discussion .................................................................................................................................. 31
Chapter 4: Machine-Learning-Aided Simulations .................................................................. 33
4.1 Introduction ................................................................................................................................ 33
4.2 Methods ...................................................................................................................................... 36
4.2.1 Data Collection and Preparation .......................................................................................................... 37
4.2.2 Feature Vector Selection ...................................................................................................................... 39
4.2.3 Reference Data Preparation .................................................................................................................. 44
4.2.4 Gaussian Process Regression ............................................................................................................... 46
4.2.5 Bayesian Optimization ......................................................................................................................... 46
4.3 Results ........................................................................................................................................ 49
iv
4.3.1 Gaussian Process Regression ............................................................................................................... 49
4.3.2 Bayesian Optimization ......................................................................................................................... 54
4.4 Discussion .................................................................................................................................. 60
Chapter 5: Dynamic Simulations on Quantum Computers ..................................................... 62
5.1 Introduction ................................................................................................................................ 62
5.2 Methods ...................................................................................................................................... 66
5.2.1 Quantum Circuits for Dynamic Simulation ......................................................................................... 69
5.2.2 Wavefunction Simulator ...................................................................................................................... 71
5.2.3 Simulated Noisy Qubits ....................................................................................................................... 71
5.2.4 Currently Available NISQ Quantum Computers ................................................................................. 72
5.2.5 First Principles Calculation of Exchange Interaction Term ................................................................. 74
5.2.6 Compilation Methods ........................................................................................................................... 76
5.3 Results ........................................................................................................................................ 78
5.3.1 Material Simulation on Quantum Computers ...................................................................................... 78
5.3.2 Domain-Specific Quantum Circuit Compiler ...................................................................................... 83
5.4 Discussion .................................................................................................................................. 88
Chapter 6: Conclusion ............................................................................................................. 91
References ..................................................................................................................................... 94
v
List of Tables
Table 1. Electronegativity (EN), first ionization potential (IP) and atomic radius (AR) values for the five
elements of which the hetero-structures in this study are comprised, along with their normalized values
given in parentheses ..................................................................................................................................................... 40
Table 2. Components and dimension of the various feature vectors tested. ................................................................ 41
Table 3. Mean square errors on test data for various feature vectors, with varying sizes of training data set
created from the 126 unique 3-layered structures. In bold is the selected feature vector for our study. ................... 43
Table 4. Percentages of the 500 BO runs in which the optimal structure is correctly identified for models
built using various feature vectors shown in column 1, for the various target properties given in columns 2-
4. The best accuracy for each target property is shown in boldface font. ................................................................... 44
Table 5. Pseudocode for the Bayesian optimization algorithm and details on use of acquisition function. ............... 47
Table 6. Pseudocode for dynamic simulation of qubits. For each simulated time-step, a different circuit is
created. Next, in a user-define number of independent trials, the qubits are initialized to their t=0 state, the
particular circuit for the given time step is executed, and the qubits are measured. Measurements from all
trial runs are averaged together to give an estimate for the final quantum state of the system at each time-
step. .............................................................................................................................................................................. 70
Table 7. Decoherence and fidelity statistics for the qubits used on the IBM quantum processor. .............................. 73
Table 8. Decoherence and fidelity statistics for the lattice used on the Rigetti quantum processor. .......................... 74
Table 9. Algebraic identities for the most used gate sets in TFIM circuits. ................................................................ 85
Table 10. Pseudocode detailing the order in which algebraic identities are applied for compilation of
TFIM circuits for execution on IBM’s quantum computer. ......................................................................................... 86
vi
List of Figures
Figure 1. (a) Top and (b) side views of the simulated MoSe2 monolayer. (c) Schematic of the simulated
photoexcitation (left), and time evolution of the resulting electronic energies and occupations (right). Here,
eight electrons from the valence-band maximum are promoted to the conduction-band minimum. Red, blue,
and black lines denote electronic levels occupied by two electrons, single electron, and no electron,
respectively. Electronic levels with nontrivial dynamics are shown in thick lines. ..................................................... 21
Figure 2. Electronic structure of monolayer MoSe2. (a) Band structure plot denoting by color the character
of the Mo d-orbitals that make up the majority of the state at each point, as well as a partial density of
states of the different d-orbitals of Mo. The top valence (conduction) band is denoted with a green (red),
dashed circle, indicating the state is predominantly composed of 𝑑𝑧2 (𝑑𝑥𝑦 and 𝑑𝑥2−𝑦2) orbitals. (b)
Charge densities for the top two valence bands and the bottom two conduction bands, showing how the
electronic charge density changes when we excite electrons from their ground state. ............................................... 24
Figure 3. Time evolution of Debye−Waller factor along the (110) and (300) planes, with (red) and without
(blue) photoexcitation, averaged over three simulation runs beginning from identical positions but different
initial velocities in the micro-canonical ensemble. The shading shows the standard deviation across the
three independent simulations. ..................................................................................................................................... 25
Figure 4. Debye Waller Factors of the (110) and (300) planes in the monolayer MoSe2 crystal calculated
at 300K. This DWF value differs noticeably from 1 at this temperature, in contrast to the near-unity value
observed at 10K. The shading shows the standard deviation across the three independent simulations. ................... 26
Figure 5. (a) Electronic density of states of MoSe2 monolayer near the band edge (blue line). Red line
shows integrated density of states, which is used to determine the Fermi level for a given value of
electron−hole pair density, n(e−h). (b) Fermi surface for excited carrier density, n(e−h), values ranging
from 0 to 2×10
14
cm
−2
. The Fermi surface is localized at the K-points at minimal excitation (red contours),
while exposing Σ-pockets at higher excited electron−hole densities (black and blue contours). The three
nesting vectors 𝑞𝑛1, 𝑞𝑛2, and 𝑞𝑛3. correspond to reciprocal vectors ΓM, ΓK, and ΓΣ, respectively. (c)
Phonon dispersion curves in the ground state and two increasingly excited states. ................................................... 27
Figure 6. Debye-Waller factor plotter versus time along the (110) (a) and (300) (b) planes, at two different
electron-hole densities, in the canonical ensemble. ..................................................................................................... 29
Figure 7. Debye-Waller factor plotted versus time along the (110) and (300) planes, for a simulation cell
with 48 atoms (red) and 108 atoms (blue) upon photoexcitation, each averaged over three simulation runs
beginning from identical position but different initial velocities in the canonical ensemble. ..................................... 31
Figure 8. Workflow for optimal structure and property prediction. First, structure files for a family of
three-layered materials are created and uploaded to the Materials Project (MP) database. Second, the MP
infrastructure performs all DFT calculations, and subsequently, transport calculations using BoltzTraP
code are performed. A snapshot of the material property data computed by MP database is pictured, along
with the thermoelectric parameters computed by BoltzTraP. Third, a numerical feature vector is assigned
vii
to uniquely represent each structure. Fourth, machine learning models are designed based on the data to
make predictions for either a material property or an optimal structure. ................................................................... 37
Figure 9. (a) Electronic band structure, showing the CBM dispersion curve highlighted in red, the VBM
dispersion curve highlighted in green, and the band gap denoted by the magenta arrow. (b)
Thermoelectric EFF versus carrier-dopant concentration for both n-type (green) and p-type (purple). ................... 38
Figure 10. Histogram of the number of iterations required to discover one of the five structures with the
highest band gap in the 500 independent BO runs. ..................................................................................................... 48
Figure 11. Histogram of number of iterations required to discover one of the five structures with the
highest peak thermoelectric EFF value in the 500 independent BO runs. .................................................................. 49
Figure 12. Effect of training set size on predicted band gap and its confidence interval on test data set for
two different training data set sizes. ............................................................................................................................. 50
Figure 13. Effect of training data set size on the predicted CBM and VBM for one of the three-layered
hetero-structures. (a) With 40% training set size, mean square error (MSE) between the predicted and
ground truth is 0.177 and 0.153 for CBM and VBM, respectively. (b) With 70% training set size, MSE
decreases to 0.085 and 0.067 for CBM and VBM, respectively. .................................................................................. 51
Figure 14. Effect of training dataset size on the predicted EFF function for one of the p-doped three-
layered hetero-structure. (a) Total training data set size used here is 40% of 126 structures. Mean square
error between the ground truth and predicted EFF is 2.78. (b) The same as (a) but with 70% training set
size, mean square error on the same structure decreases to 1.22. .............................................................................. 51
Figure 15. (a) Predicted (blue) and ground truth (red) band gap values of three-layer hetero-structures in
the test data set. Yellow shading represents a 95% confidence interval (CI) of the predicted results. (b)
Predicted (dashed lines) and ground truth (solid lines) values for the VBM and CBM dispersion curves in
momentum space for a sample three-layered hetero-structure. (c, d) Predicted (dashed lines) and ground
truth (solid lines) thermoelectric EFF curves versus carrier concentration for a sample n(p)-type-doped
three-layered hetero-structure. In all models, 60% of the three-layered structures were randomly selected
to comprise the training data set. ................................................................................................................................. 53
Figure 16. (a) Band gap values for all three-layered hetero-structures, where hetero-structures on the x-
axis are ordered by increasing band gap value. Highlighted points denote hetero-structures returned most
frequently as the optimal structure by a BO model searching for the maximum band gap (red) and a
desired value of 1.1 eV (green). (b) Pie chart showing the distribution of optimal band gap values returned
in 500 independent BO searches for the maximum band gap. (c) Pie chart showing the distribution of
optimal band gap found in 500 independent BO searches for a desired band gap value of 1.1 eV. MoSe2-
WSe2-WSe2, with a band gap of 1.05 eV, has the closest band gap value to 1.1 eV and was returned as the
optimal structure in 91% of the 500 independent BO searches. WSe2-WSe2-WSe2 and MoS2-MoS2-WS2, with
band gap values of 1.23 and 1.26 eV, respectively, were returned in 7% of the BO searches, while WSe2-
MoSe2-WSe2 and MoSe2-WSe2-MoSe2, with band gap values of 1.01 and 1.04 eV, respectively, were
returned in the remaining 2% of BO searches. ............................................................................................................ 55
viii
Figure 17. Initial and final maximum band gap values for ten different BO runs on the class of four-layered
hetero-structures. Optimal structures, along with their band gap value are labeled .................................................. 56
Figure 18. (a, b) Thermoelectric EFF values for all p-type-doped (a) and n-type-doped (b) three-layered
hetero-structures, where hetero-structures on the x-axis are ordered by increasing EFF value. Points
highlighted in red and green denote hetero-structures returned most frequently as the optimal structure by
a BO model searching for the p-type and n-type-doped structures, respectively, with maximum EFF value.
(c, d) Pie chart showing the distribution of optimal thermoelectric structures returned in 500 independent
BO searches for the p-type-doped (c) and n- type-doped (d) materials with maximum EFF value. ........................... 57
Figure 19. Band structures along lines of high symmetry in the first Brillouin zone for (a) MoSe2-WS2-WS2
and (b) WSe2-WTe2-WSe2, which are the two best n-type materials, and (c) WTe2-MoTe2-WTe2 and (d)
MoSe2-WSe2-WSe2, which are the two best p-type materials based on global optimization of the maximum
electronic fitness function. The red dots represent conduction band minima, while the green dots denote
valence band maxima in all four panels. ...................................................................................................................... 58
Figure 20. Schematic showing how a simplified model of Re-doped monolayer MoSe2 is mapped onto the
qubits of a quantum computer. (a) Top-down view of the Re-doped MoSe2 monolayer where Mo, Se, and
Re atoms are depicted by pink, yellow, and purple spheres, respectively. Grey arrows are superposed on
the Re atoms, representing the spin of the extra, unpaired electron each Re atom possesses. Inset shows a
side view of the material and a representation of the excited 𝐸′′ phonon mode. (b) Spins of the extra,
unpaired electron of each Re atom (grey arrows), with exchange interactions between neighboring spins of
strength Jz, in the presence of an external magnetic field with frequency 𝜔𝑝ℎ and amplitude 𝜀𝑝ℎ are
mapped onto the qubits of a quantum computer, shown in their Bloch sphere representation. .................................. 67
Figure 21. Quantum circuit diagram simulating time evolution to the first time-step for a three-qubit
system. .......................................................................................................................................................................... 69
Figure 22. Topology of the IBM Q16 Melbourne quantum processor. Blue circles represent qubits, while
black lines represent physical connections between pairs of qubits. ........................................................................... 73
Figure 23. Topology of the Rigetti Aspen quantum processor. Blue circles represent qubits, while black
lines represent physical connections between pairs of qubits. ..................................................................................... 74
Figure 24. Schematic describing the calculation of the exchange interaction term Jz between unpaired
electron spins in neighboring Re atoms. Re atoms are depicted by cyan spheres, while Se atoms are
depicted by yellow spheres. The red block arrows indicate the unpaired electron of the Re atom is spin-up,
while the blue block arrows denote spin-down. ........................................................................................................... 75
Figure 25. Time evolution of the average magnetization of 2- (red), 3- (green), and 4-qubit (blue) systems
with electron-phonon coupling strengths 𝜀𝑝ℎ = 0.2𝐽𝑧 (a),0.5𝐽𝑧 (b), 𝐽𝑧 (c) and 5𝐽𝑧 (d). The black dotted
line shows zero magnetization. ..................................................................................................................................... 79
Figure 26. Simulation results for a 2-qubits system (red dots) on the IBM quantum processor (a-d) and the
Rigetti Aspen quantum processor (e-h), compared to theoretical results from simulated noisy qubits (black,
dashed lines) and the wavefunction simulator (black, solid lines). The black dotted lines show zero average
ix
magnetization. Results are shown for varying electron-phonon coupling strengths 𝜀𝑝ℎ = 0.2𝐽𝑧 (a,e),0.5𝐽𝑧
(b,f), 𝐽𝑧 (c,g) and 5𝐽𝑧 (d,h). .......................................................................................................................................... 80
Figure 27. Results of a TFIM simulation run on Rigetti’s quantum computer, based on circuits produced by
a naïve compiler (blue) and an optimized compiler (red) for a two-qubit system. The black line represents
numerical results from a noiseless, simulated quantum computer. ............................................................................. 82
Figure 28. General structure of circuits for dynamic simulation. (a) This unit of unitary gates performs
evolution over one time-step of size ∆𝑡 of the spins under the TFIM Hamiltonian. (b) The unit of gates
shown in (a) is repeated a total of n times to evolve the system through n time steps. Notice how, for a
time-dependent Hamiltonian, parameters within each unit change from time-step to time-step. ............................... 84
Figure 30. Performance comparison of the domain-specific (DS) compiler to IBM’s general-purpose
compiler. (a) Absolute gate count reduction using the DS compiler over IBM’s compiler. (b) Gate count
reduction using the DS compiler as a percentage of the gate count of IBM’s compiler. (c) Wall-clock
compiler time reduction using the DS compiler as a percentage of the wall-clock compiler time of IBM’s
compiler. ....................................................................................................................................................................... 87
x
Abstract
Monolayers of semiconducting transition metal dichalcogenides (TMDC) are emerging as
strong candidate materials for next-generation electronic devices. Interest in these two-
dimensional (2D) quantum materials stems from the ability to control their structural and
electronic properties via myriad schemes including mechanical strain, intercalation, carrier
doping, optical and phonon excitation, combination into hetero-structures, and alloying. Since
the integration of new materials into electronic devices often requires specific values for
numerous material properties, broad control over 2D TMDCs makes them highly amenable to
adaptation for new devices. For practical implementation, a deeper understanding of the
mechanisms and capabilities of control over different material properties at the atomic level is
essential. Achieving this level of accuracy with experimental study is difficult, often impossible.
Therefore, computer simulation of the control over material properties in quantum materials has
proven to be an invaluable tool. This thesis examines three methods of controlling different
material properties in 2D TMDCs using three different simulation techniques.
First, we investigate control over structural phase in 2D TMDCs via laser excitation using
non-adiabatic quantum molecular dynamics simulations. Simulation of a monolayer after optical
excitation showed a sub-picosecond disordering of the lattice, as measured by the Debye−Waller
factor, as well as increasing disorder for higher densities of photogenerated electron−hole pairs.
Furthermore, electronic structure analysis showed that the ultrafast, photo-induced lattice
dynamics arise from strong electron-phonon coupling resulting from electronic Fermi surface
nesting and its associated phonon-mode softening. Such strong lattice response to optical
excitation indicates that optical control over structural phase in 2D TMDCs is a possibility. The
xi
mechanistic understanding uncovered with these simulations helps guide optical control for
emerging applications such as phase change memories and neuromorphic computing.
Second, we investigate control over thermoelectric fitness in hetero-structures, comprised of
vertically stacked TMDC monolayers, via stacking order using machine-learning-aided
simulations. We developed a machine learning algorithm based on Bayesian optimization that
enables discovery of the optimal stacking order for thermoelectric fitness with high probability
using minimal quantum simulations. Given the remarkable success of our algorithm for such a
complex material property, we expect our algorithm to show similar success in finding optimal
materials for various other applications. This ability to guide expensive quantum simulations
with machine learning paves the way for discovering 2D TMDC hetero-structures for use in
future devices, which may have vastly different material requirements from current norms.
Finally, we investigate the control over emergent magnetism in 2D TMDCs via phonon
excitation using quantum simulations performed on quantum computers. We demonstrate
successful simulation of the control over magnetism by THz radiation of 2D TMDCs on IBM’s
Q16 Melbourne quantum processor and Rigetti’s Aspen quantum processor. Remarkably good
quantitative agreement was found between the results from the quantum computer and those
from simulated noisy qubits for a two-qubit system. The quality of the results can be attributed
to a special technique used by general-purpose quantum circuit compilers for two-qubit systems
to create short, constant-depth circuits for all time-steps of a simulation. Since this compilation
technique does not scale with system size, we developed a scalable, domain-specific compiler to
aid in circuit size reduction for simulation of larger systems. We anticipate the domain-specific
compiler will enable larger simulations (i.e. longer simulation times or larger system sizes) on
near-future quantum computers that would not otherwise be possible. Together, our early proof-
xii
of-concept simulations and domain-specific compiler lay the groundwork for simulations
performed on near-future quantum computers. Such simulations will deepen our understanding
of controlling magnetism, as well as other controllable electronic properties, in 2D TMDCs for
use in myriad future technologies.
The last half-century has witnessed exponential improvements in electronic devices, due
in main part to the miniaturization of their components. As we begin to reach the limits of such
miniaturization the search for alternate technologies and materials for next-generation electronic
devices becomes essential. The myriad controllable material properties of 2D TMDCs make
them excellent candidate materials for new quantum technologies. This thesis demonstrates the
viability of three techniques for controlling material properties in 2D TMDCs, paving the way
towards their integration into post-silicon quantum technologies.
1
Chapter 1: Introduction
Quantum materials, or materials that experience phenomena which can only be explained
by quantum mechanics, promise to bring about myriad new technologies due to the unusual
properties they exhibit, not found in their classical counterparts. Harnessing the potential of
these uniquely quantum properties, however, requires an understanding of the material at the
electronic level. Probing such materials at this level of accuracy in the lab is difficult at best,
impossible at worst. Computer simulation of both static and dynamic properties of quantum
materials, therefore, has proven to be an invaluable tool, enabling the theoretical study of how
electrons behave in these materials, and importantly, how they give rise to useful phenomena.
Two-dimensional (2D) quantum materials are excellent candidates for study with
simulation as they tend to require fewer atoms in the simulation cell than their 3D counterparts,
thus better fitting within computational resource constraints. Furthermore, 2D materials exhibit
interesting quantum properties arising from their spatial constraint in one dimension. Since the
advent of 2D graphene [1], a wide variety of 2D materials have been isolated with a suite of
interesting properties and applications [2]. In particular, 2D monolayers of semi-conducting
transition metal dichalcogenides (TMDCs) hold great promise for use in new quantum
technologies due to their facile synthesis and unique intrinsic properties [3], as well as highly
tunable electronic and structural properties through various strategies, including mechanical
strain [4], intercalation [5], carrier doping [6], optical excitation [7], combination into hetero-
structures [8], and alloying [9]. Novel devices often rely on precise values for numerous
material properties, which means the wide tunability of various attributes of 2D TMDCs make
them highly attractive as designer materials. A deeper understanding of the control of material
2
properties is therefore essential for integrating TMDCs into next-generation quantum
technologies.
In this thesis, we explore the optical control of structural phase in monolayer TMDCs, the
control of electronic properties such a band gap and thermoelectric fitness via the stacking order
of TMDC monolayers into vertical heterostructures, and the control of quantum magnetism in
doped TMDCs monolayers via phonon excitation. To showcase the variety of state-of-the-art
methods in simulation, we investigate each of these tunable properties using one of three
different simulation techniques: (i) non-adiabatic quantum molecular dynamics simulations on
classical supercomputers, (ii) machine-learning-aided static simulations on classical
supercomputers, and (iii) quantum dynamics simulations on currently available quantum
computers. We briefly elaborate on each of these simulations below, providing a basic
explanation of the method, its advantages and disadvantages, and a description of the particular
simulation performed with the method.
1.1 Non-Adiabatic Quantum Molecular Dynamics Simulations
Quantum molecular dynamics (QMD) simulations follow the trajectories of each atom in
the simulation cell, computing the interatomic forces quantum mechanically. Generally, these
are ground state calculations, in which the electrons adiabatically remain in their ground state as
the positions of the atomic nuclei are updated each time step. In non-adiabatic QMD (NAQMD),
electrons can occupy excited states and stochastically hop between states. The advantage of
QMD simulations is that they allow for observation of the movement of the electrons and nuclei
without disturbing the system. Moreover, such observations can be recorded on arbitrarily small
spatial and temporal scales, not achievable in experiment. Furthermore, the initial conditions of
3
the system as well as its external environment can be precisely defined and controlled. The main
advantage of NAQMD simulations in particular is that they allow for the simulation of the
dynamics of non-equilibrium systems, such as those that have been photo-excited. The
disadvantage of this method is its inability to scale to large system sizes (due to computational
expense) or systems with strong electron correlations (due to approximations used for exchange
and correlation interactions amongst electrons). While a full-machine run on one of the world’s
largest supercomputers can accommodate a simulation of 𝒪(10
!
) atoms [10], more reasonable
simulations generally only consist of 𝒪(10
"
) atoms [7,11]. Still, we show that NAQMD
simulations are able to give insight into the dynamics of relevant, optically excited materials.
In this thesis, NAQMD simulations are performed to study the lattice dynamics of a
model TMDC monolayer of MoSe2 after optical excitation. Using optical excitation to modulate
crystal structure is a promising technique for integrating 2D semi-conductors into new quantum
technologies in cases where methods such as chemical doping or strain engineering are
ineffective. Capitalizing on the ability of the NAQMD method to simulate a material with
varying levels of electronic excitation, we explore the dependence of lattice dynamics on the
photo-generated electron-hole density. The simulations show sub-picosecond disordering of the
lattice upon photo-excitation, as measured by the Debye-Waller factor, with increasing disorder
for higher densities of photo-generated electron-hole pairs. Further electronic structure analysis
shows that the rapid, photo-induced lattice dynamics are due to strong electron-phonon coupling
arising from electronic Fermi-surface nesting. These findings suggest the potential for optical
control of the structural phase of TMDC monolayers. Understanding the sub-picosecond atomic
dynamics after optical excitation is important for the realization of optical control of structure in
monolayer TMDCs to enable emerging applications such as phase change memories and
4
neuromorphic computing [12]. Chapter 3 provides more details about the NAQMD simulation
method, the simulation of photo-excited monolayer MoSe2, as well as an in-depth presentation
and discussion of the simulation results.
1.2 Machine-Learning-Aided Simulations
Over the last decade, machine learning methods have shown phenomenal success in the
materials science domain for materials screening and property prediction [13-15]. In the general
approach, material simulations are performed for a subset of a class of materials to create a set of
data used to train a machine learning model to predict the material properties of the remainder of
the class. A problem for which machine learning techniques are particularly useful is the
discovery of the optimal material (from within a class of materials) with respect to some
property. Without machine learning, exhaustive screening of the material search space is
usually the only means to a solution, which quickly goes beyond the bounds of existing
computational resources, even when performed in parallel on large supercomputers.
An alternative approach to brute force search involves performing property calculations
for only a subset of the materials in the search space to populate a data set used to train a
machine learning model to predict the property values for the remaining materials. In this way,
one could either simulate or predict the value of a material property for every material in the
search space to determine the optimal one, without exhaustively performing simulations for
every material. However, this can be problematic for complex material properties, since accurate
prediction can require computation of properties for a large percentage of the materials for the
training data set. For large classes of materials, performing simulations for a large percentage of
the class may still be too costly.
5
To limit the number of simulations that need to be performed, a branch of machine
learning, known as active learning, can be used. In this approach, simulations are performed for
a very small number of materials from the materials class, which are randomly selected to act as
the initial training data set. The active learning algorithm then uses a utility function, derived
from the initial training data set, to select the next best material to add to the training data set.
This material is the “best” to simulate in the sense that it either (i) exists in an area of the
material search space about which the model has little information (i.e. an exploratory choice), or
(ii) is a candidate for the optimal material following prediction trends developed by the model
(i.e. an exploitative choice). Upon performing a simulation of the selected “best” material and
adding it to the training data set, the utility function is updated, and the next best material to
simulate is determined. A small number of iterations of this process is performed, until there is a
high probability that the optimal material has been discovered using minimal material
simulations. The advantage of this method is that it greatly reduces the number of
computationally expensive simulations of quantum materials that need to be performed. The
disadvantage is that the returned optimal structure only has a high probability of being such, it is
not guaranteed.
In this thesis, we present an active learning algorithm used to expedite the discovery of
the optimal layered TMDC hetero-structure with respect to thermoelectric fitness and electronic
structure properties. These hetero-structures, comprised of vertically stacked, heterogenous
monolayers of TMDCs, are promising candidates for next generation optoelectronic and
thermoelectric devices. We show that the active learning method significantly reduces the
computational cost of discovering the optimal structure. Chapter 4 provides a description of the
6
TMDC hetero-structures, the estimator for their thermoelectric figure of merit, a technical
description of the machine learning methods used, as well as a detailed analysis of the results.
1.3 Dynamic Simulations on Quantum Computers
The newest of the three simulation techniques examined in this thesis, dynamic
simulations on quantum computers, is still in its nascent stages. The present era of quantum
computing has been dubbed the Noisy Intermediate-Scale Quantum (NISQ) [16] era, as the
currently available quantum computers have too few qubits to be able to address noise and
decoherence with error-correcting schemes. As such, materials simulations on quantum
computers have not led to any new discoveries to date, but the potential for quantum computers
to perform simulations that are effectively impossible on the largest classical supercomputers is
highly anticipated in the next few decades.
Quantum computers perform computation with two-state, quantum-mechanical systems
that serve as quantum bits, or qubits. Qubits take advantage of purely quantum mechanical
properties, such as superposition and entanglement, to outperform the most advanced classical
supercomputers for certain classes of problems. Over the past few decades, they have quickly
transformed from “imaginary devices” [17] to one of the most heavily funded technologies, with
newly engineered hardware and software emerging from Fortune 500 companies to start-ups to
government-backed centers across the globe. The National Quantum Initiative Act recently
signed into law by the 115
th
United States Congress, which seeks to “ensure the continued
leadership of the United States in quantum information science and its technology applications”
[18], exemplifies the importance of this novel technology to both science and national security.
Despite the recent experimental demonstration of quantum supremacy for a computational task
7
using random circuits with no known physical applications [19], the challenge of using a
quantum computer to solve a physically relevant scientific problem still remains.
The last decade has witnessed the growing success of quantum computers for simulating
static properties of quantum systems, i.e., the ground state energy of small molecules. However,
it remains a challenge to simulate quantum many-body dynamics on current-to-near-future NISQ
computers. Such dynamic simulations pose a challenge as current algorithms produce quantum
circuits that grow in size with each subsequent time-step of the simulation. Due to their high
gate-error rates and short decoherence times, however, NISQ computers can only produce high-
fidelity results for those quantum circuits smaller than some given circuit size. Thus, dynamic
simulations on NISQ computers have an upper limit on the number of time steps they can
simulate. This underscores the crucial role of quantum circuit compilers to produce executable
quantum circuits of minimal size, thereby maximizing the range of physical phenomena that can
be studied within the current NISQ fidelity budget.
The main advantage of performing simulations of quantum materials on quantum
computers is that they can efficiently store the wavefunction of a quantum many-body system.
Defining the many-body wavefunction of a quantum system entails storing a complex number
for each basis state of the system, which is associated with the probability of measuring the
system to be in that state. The number of basis states grows exponentially with the number of
particles in the system, so a classical computer will need a number of bits that increases
exponentially with the number of particles to store the wavefunction. However, a quantum
computer only requires a number of qubits that grows linearly with the number of particles in the
system. For example, a system with 50 spin-
#
"
particles will have 2
$%
≈1.13×10
#$
basis states,
which, using single precision floating point numbers to store wavefunction amplitudes would
8
require over 9 petabytes of memory on a classical computer. Adding just one more particle to
the system, for a total of 51 spins, would double the classical memory requirement to 18
petabytes, making the simulation impossible to run on even the world’s most advanced
supercomputer. On the other hand, the many-body wavefunction of a 50- and 51-spin system
can be stored on 50 and 51 qubits, respectively.
The disadvantage of performing simulations on quantum computers in the NISQ era is
the fact that current quantum computers are small and noisy, making it difficult to perform
simulations of systems large enough to discover new physics. A disadvantage of performing
simulations on quantum computers in general is the measurement problem; it is not possible to
measure the full quantum state of a system, instead, measurement causes the system to collapse
to one of the eigenstates of the measurement operator. Therefore, reconstructing the full many-
body wavefunction of a system requires a large number of measurements on an ensemble of
identical systems.
In this thesis, we demonstrate a proof-of-concept for a successful simulation of nontrivial
quantum dynamics on IBM’s Q16 Melbourne quantum processor and Rigetti’s Aspen quantum
processor; namely, the control of quantum magnetism via phonon excitation induced by THz
radiation in a 2D TMDC monolayer. Furthermore, we present a domain-specific quantum circuit
compiler, specifically designed for circuits simulating spin dynamics under a special class of
Hamiltonians, that outperforms the state-of-the-art general-purpose compiler in terms of circuit
size minimization as well as wall-clock compilation time. We find increasing gate reductions for
circuits simulating larger numbers of timesteps as well as for larger number of particles in the
simulation system, showing that our compilers scale with number of simulation time-steps and
system size.
9
As new quantum hardware continues to see improvements in decoherence times and gate
error-rates, hence allowing for larger circuit sizes, we anticipate that this domain-specific
compiler will enable simulations of larger systems with more timesteps not otherwise achievable
with general-purpose compilers. As such, this work lays a foundation for the promising study of
a wide variety of quantum dynamics on near-future quantum computers, including dynamic
localization of Floquet states and topological protection of qubits in noisy environments.
Chapter 5 provides a detailed explanation of how dynamic simulations are generally carried out
on current NISQ computers, as well as descriptions of the specific quantum simulations
performed and the domain-specific quantum circuit compiler developed.
10
Chapter 2: Two-Dimensional Transition Metal Dichalcogenides
Bulk elemental and compound semiconductors have provided the basis for an astounding
array of information technology devices over the last half century. However, improvements in
miniaturization and efficiency have begun to plateau as we reach the limit at which three-
dimensional structures cannot get any smaller without beginning to experience quantum
confinement in one or more degrees of freedom. As a result, smaller, faster, lower-power,
flexible, or transparent devices may require a paradigm shift from three-dimensional to 2D
components.
Of particular interest for 2D components of electronic, optoelectronic, and thermoelectric
devices are monolayers of semi-conducting TMDCs, due to their chemically stability, high
carrier mobilities [20,21], tunable electronic and transport properties [22-25], and their facile
synthesis [3]. Because of their 2D quantum confinement, these monolayers possess electronic
structures and properties not observed in their bulk counterparts. Examples include emergent
valley- and spin-dependent optical and electrical properties [26-30], making these monolayers
suitable for valleytronics; tunable semiconducting and metallic structural phases [31-34],
potentially useful for phase change memories; and strong light-matter interactions [35,36], which
qualify 2D TMDCs for photonics and photocatalytic applications [37,38]. Finally, their
resilience to large elastic strain [39-42] speaks to their potential in transparent and flexible
electronic devices.
Another important aspect of these monolayers is that they can exist in two distinct crystal
structures which possess markedly differently physical and electronic properties. The first, more
ubiquitous crystalline phase, labeled the H polytype, is characterized by a trigonal prismatic
arrangement of the chalcogen atoms around the transition metal atoms. In this atomic formation,
11
the monolayer exhibits semi-conducting behavior with a direct band gap. Furthermore, the lack
of inversion symmetry in this configuration leads to valley-dependent electronic properties. In
the second polytype, labeled T’, the chalcogen atoms have a distorted octahedral coordination
around the transition metal atoms. In this configuration, the material has semi-metallic
properties with complex topological states [3] and regains inversion symmetry. The precision
associated with the optical control of transformation between these two structural phases would
enable patterning in monolayers leading to the formation of defect-free semi-conductor-to-metal
hetero-phase homojunctions. Such junctions would ensure ohmic contact in 2D semi-conducting
devices, allowing for the integration of monolayers of TMDCs into new 2D technologies.
In this thesis, we examine three methods to control properties of 2D TMDCs to prime
them for integration into next-generation technologies based on 2D devices. The first method
involves optical modulation of the crystal structure of TMDC monolayers. Prior studies have
demonstrated unusually strong light-matter interactions in 2D TMDCs, which bolsters the
prospect of using optical excitation to control their structural phase. However, the nature of
electronic and structural dynamics in TMDCs in response to electronic excitation has not yet
been fully elucidated. Electronic excitation has previously been shown to induce bond
dissociation and atomic rearrangement on ultra-fast timescales in other material systems like
graphene [43], carbon nanotubes [44], transition metal oxides [45] and polymers [46]. Recent
experimental studies on atomically thin TMDCs have shown strong structural response to optical
and electronic excitation [47,48], often characterized by the activation of specific phonon modes.
In particular, ultrafast pump-probe experiments on bilayers of MoSe2 have shown ultrafast
conversion of electronic excitation to thermal energy of the lattice with near-unity quantum yield
12
[49]. Such ultrafast energy conversion is necessary for optical control of structural phase for
applications like phase change memories and neuromorphic computing [12].
The second method involves vertically stacking heterogeneous TMDC monolayers to
tune properties for various for next-generation optoelectronic and thermoelectric devices.
Stacking monolayers in this fashion opens the door to a new class of hybrid materials where one
can, in principle, engineer electronic, transport, optical and other properties in well-defined,
controllable material structures. An important aspect of these layered materials is the discrete
nature of the van der Waals (vdW) forces that bond the layers. This weak vdW bonding
facilitates synthesis of 2D monolayers from bulk via mechanical or chemical exfoliation. It also
allows for stacking of lattice-mismatched monolayers from different species of TMDCs, thereby
enabling the formation of unlimited combination of multi-layered hetero-structures. The
interlayer interactions are nonetheless sizable enough to strongly affect the electronic behavior.
Therefore, the electronic properties of these multi-layered hetero-structures can vary in a highly
nontrivial manner depending on the total number of layers, the specific layer ordering, and each
layer’s composition. Thus, these hetero-structures can, in principle, be engineered for any
desired combination of material properties important for specific electronic and thermoelectric
applications.
Good thermoelectrics have highly complex band structures, markedly different from
standard isotropic parabolic bands, on which the thermoelectric performance depends.
Importantly, the concept of reduced dimensionality to achieve high thermoelectric performance
has been very influential in thermoelectrics [50], but has proven difficult to achieve in practical
materials. Vertically stacked TMDC monolayers, which exhibit reduced dimensionality in each
13
layer and highly complex band structures due to stacking different species of TMDCs, therefore
have the potential to produce designer materials optimized for thermoelectrics.
The third and final method of control over material properties in 2D TMDCs presented in
this thesis is the control of emergent magnetism in substitutionally-doped TMDC monolayers.
While defect-free, ideal TMDC crystals are non-magnetic due to the lack of unsaturated bonds
and non-magnetic elements, there exists a wealth of research on the emergence of magnetism in
real 2D TMDC crystals arising from crystalline imperfections from edges, grain boundaries,
magnetic and non-magnetic dopants, and defects as well as from strain [51]. Of particular
interest for ultrafast control of magnetism is the recent theoretical prediction of emergent
magnetization in a TMDC monolayer induced by the excitation of a specific phonon mode [52].
Since THz photoexcitation of such systems has been shown to excite specific phonon modes on
sub-picosecond time scales [7,49], in principle, this technique could be used for optical control
of magnetism in 2D TMDCs on ultrafast timescales. Such control paves the way to the creation
of 2D magnetic semi-conductors, which, by integrating low-dimensionality, ferromagnetism, and
semi-conductivity, provide a cornerstone for novel high-speed, spintronic devices.
In the following chapters we perform both static and dynamic simulations of 2D TMDCs
to study the control of material properties in 2D TMDCs using three different simulation
techniques. We find that optical excitation of these materials has a tunable effect on atomic
displacement and phonon excitation, which in turn can lead to control over the structural phase
and emergent magnetism of the material, respectively. We also find that electronic and
thermoelectric properties can be tuned by stacking heterogenous monolayers of TMDCs in
different permutations. Such simulations help to guide experimentalists in their attempts to
optically control TMDCs and in their search for new hybrid- or meta-materials. All simulations
14
show that 2D TMDCs are promising candidate materials for next-generation quantum
technologies.
15
Chapter 3: Non-Adiabatic Quantum Molecular Dynamics
Simulations
3.1 Introduction
QMD is a well-established method for simulating the dynamic behavior of quantum
materials. After providing the species and initial positions of the atoms in the simulation cell,
QMD simulations step through time, tracking the trajectories of each atom while incrementally
updating their positions according to quantum mechanically derived inter-atomic forces. QMD
is an ab initio, or first principles, method, meaning that no empirical data is used to fit
parameters, rather only the equations of quantum and classical mechanics are used. Such
simulations provide two main advantages over experiment. First, they record the movement of
each atom in the system across every time step with arbitrary spatial and temporal precision,
without disturbing the system. Capturing such detail with experiment can be extremely difficult
and comes with limited precision. Furthermore, the measuring device, by definition, must
interact with the system, altering the future dynamics of the system. Second, the parameters of
both the system and its external environment are precisely defined and can be finely controlled in
QMD simulations. Settings such as initial positions of the atoms, the temperature of the heat
bath, the initial excitation of the electrons, etc. can be exactly defined and varied across runs to
study how turning different knobs of the system affects behavior. In experiment, many of these
factors must be approximated at best, and cannot be guaranteed to be constant across an
ensemble of experiments.
Conventionally, QMD simulations allow for adiabatic dynamics that assume the electrons
remain in their ground state throughout the simulation. In NAQMD simulations, electrons can
occupy excited states, as well as stochastically hop between them. This allows for simulation of
16
the dynamics of photo-excited materials. In this thesis, we use the open source software package
QXMD [53] to perform NAQMD simulations to study the electronic and lattice dynamics of a
model TMDC monolayer, MoSe2, upon electronic excitation. Monolayer TMDCs can exist in
two structural phases: a stable semi-conducting phase, and a metastable semi-metallic phase.
This makes these monolayers viable candidates for phase-change materials for use in phase-
change memories and neuromorphic computing. For integration of these monolayers into such
applications, optical control of the structural phase would be highly beneficial. Studying the
structural dynamics of the monolayer upon photo-excitation, therefore, is an important first step
towards optically controlled phase in monolayer TMDCs.
3.2 Methods
3.2.1 Non-Adiabatic Quantum Molecular Dynamics
NAQMD simulations [54] record the dynamics of a coupled electron-ion system, where
ionic dynamics are derived from the Hellman-Feynman theorem (i.e. ionic trajectories are
propagated using Newtonian mechanics based on quantum mechanically derived forces),
electronic dynamics are modeled within a time-dependent DFT framework [55], and non-
adiabatic transitions of electrons between energy bands are computed using the surface hopping
approach [56]. By alternating steps of (1) updating excited electron transition probabilities due
to new ionic positions with steps of (2) updating ionic configuration due to new electronic
wavefunctions, we are able to calculate the non-adiabatic electron-ion dynamics in a simulated
photoexcitation.
For each time-step of the simulation, the Hamiltonian of the electronic system is
computed, based in part on the positions of the atomic nuclei. The Hamiltonian of the system is
17
then plugged into the Schrodinger equation, from which the wavefunctions for the electronic
system can be solved using the self-consistent field method. The Hamiltonian for the entire
system of N electrons with positions {𝑟
&
| 𝑖 =1,…,𝑁} and Na nuclei with positions {𝑅
'
| 𝐼 =
1,…,𝑁
(
} is given by
𝐻
8
= −
ℏ
!
"*
"
∑ ∇
&
"
&
+∑
+
#
,
!
|.
$
/0
%
|
&,'
+
#
"
∑
,
!
2.
$
/.
&
2
&34
−∑
ℏ
!
"6
%
∇
'
"
&
+
#
"
∑
+
%
+
#
,
!
20
%
/0
#
2
'37
(1)
where ℏ is Planck’s constant, MI and ZI are the mass and charge of the I
th
nucleus; me and e are
the electron mass and charge, respectively. Lower-case indices span the electrons, while upper-
case indices span the nuclei. The first term is the kinetic energy of the electrons, the second term
accounts for the electron-nuclei interactions, the third term accounts for the electron-electron
interactions, the fourth term is the kinetic energy of the nuclei, and the last term accounts for the
nuclei-nuclei interactions.
The most commonly applied approximation is the Born-Oppenheimer approximation,
which separates the motion of the nuclei from the motion of the electrons. This is based on the
fact that the masses of the nuclei are many orders of magnitude larger than the electron mass, and
therefore the nuclei move on much slower time scales. The motion of the atomic nuclei is thus
neglected when solving the Schrodinger equation for the motion of the electrons. In this picture,
the nuclei are “frozen”, and we approximate the system as interacting electrons in an external
potential created by the nuclei. Thus, the kinetic energy of the nuclei (i.e. the 4
th
tern) in
equation (1) is dropped, and the Hamiltonian can be re-written as
𝐻
8
= 𝑇
>
+ 𝑉
>
,89
+𝑉
>
&:9
+𝐸
''
(2)
where, adopting Hartree atomic units in which ℏ=𝑚
,
=𝑒 =
!;
<
'
=1, the kinetic energy
operator for the electrons, 𝑇
>
, is given as 𝑇
>
= ∑ −
#
"
∇
&
"
&
, the external potential of the nuclei
18
acting on the electrons, 𝑉
>
,89
, is given as 𝑉
>
,89
= ∑ 𝑉
'
(|𝑟
&
−𝑅
'
|)
&,'
, the electron-electron
interactions, 𝑉
>
&:9
, are given as 𝑉
>
&:9
=
#
"
∑
#
2.
$
/.
&
2
&34
, and 𝐸
''
gives the interaction of the nuclei
with one another.
The dynamics of the nuclei are then determined with classical mechanics based on
quantum mechanically derived forces, derived from the Hellman-Feynman theorem as
𝑀
'
=
!
=9
!
𝑅
'
= −
>
>0
%
〈Ψ|𝐻
8
|Ψ 〉 (3)
where the Schrodinger equation is used to solve for the electronic wavefunction given by
𝑖ℏ
=?({B
$
;9})
=9
= 𝐻
8
Ψ({r
&
;𝑡}). (4)
Due to the difficulty of solving for the many-body wavefunction, DFT, along with the Kohn-
Sham ansatz, were developed to aid in the computation of the electronic wavefunction.
Based on the theoretical work of Hohenberg and Kohn [57,58] (extended to time-
dependent systems by Runge and Gross [55]), DFT posits that (1) there is a one-to-one
correspondence between the (time-dependent) external potential of an electronic system and
(time-dependent) electron charge density and (2) all properties of the ground state, many-
electron system are functionals of the electron charge density alone. For a system in a given
external potential (i.e. configuration of atomic nuclei), this means there is a unique electron
charge density corresponding to this potential. Furthermore, all properties of the system can be
derived from this charge density.
Since all properties of the system are functionals of the electron charge density, Kohn and
Sham proposed [57] working with an auxiliary system of N non-interacting electrons, each of
which only feels the average external potential of all the other electrons, with the same charge
density as the fully interacting system. It is much simpler to solve for the N single-electron
19
wavefunctions of the non-interacting system, rather than the N-body wavefunction in equation 4,
though in principle, these single-electron wavefunctions can be used to determine all properties
of the full many-body system. In practice, however, the average external potential that the non-
interacting electrons feel can only be approximated, which is the limiting factor in the accuracy
of this theory. The potential is derived using a mean-field approximation with additional
corrections for the exchange and correlation interactions between electrons.
The wavefunctions of the electrons are represented with a plane-wave basis. In order to
reduce the number of electrons for which wavefunctions need to be computed, a second
approximation using pseudopotentials (PPs) is applied. This approximation rests on the fact that
valence electrons are the ones that mainly participate in chemical interactions between atoms,
while core electrons remain tightly bound to their respective nuclei. Thus, the PP method wraps
the core electrons of an atom into the nucleus to create a positively charged ion, making it only
necessary to compute wavefunctions for the valence electrons. A separate PP must be calculated
for each atomic species, for which there are a few methods, including norm-conserving [59],
ultrasoft [60], and projector-augmented wave method [61].
For each time step in the NAQMD simulation, the wavefunction for each valence
electron in the system is calculated. The fewest switches surface hopping method [56] is then
used to stochastically determine if any electronic transitions between energy states occur. The
transition probability between two energy bands is based on the energy difference between the
two bands and the amount of overlap of the wavefunctions of the two energy states. A random
number generator decides if any electrons transition between states based on the transition
probabilities, and the electronic occupations numbers are then updated accordingly.
20
The excited electronic states are constructed from a basis set of ground state Kohn-Sham
(KS) orbitals within DFT. The electronic ground state orbitals of the system are computed using
the projector augmented wave PP method with projector functions generated for the 4d, 5s, and
5p states of Mo and 3d, 4s, and 4p states of Se. The generalized gradient approximation [62] is
used for the exchange-correlation energy with nonlinear core corrections. Electronic pseudo-
wavefunctions and the pseudo-charge density are expanded in planewaves with 408 and 3401 eV
(30 and 250 Rydberg) cutoffs, respectively. The energy functional is minimized with respect to
KS orbitals in an iterative fashion, using a preconditioned conjugate gradient method.
3.2.2 Dynamic Simulations
The simulated system is a monolayer of MoSe2 with 108 atoms, comprised of
6 × 6 × 1 unit cells in the a-, b- and c-directions, respectively. As shown in Figure 1, each
monolayer unit cell consists of one Mo atom sandwiched by two Se atoms. Figure 1(a) shows a
top-view of the material, while Figure 1(b) shows a side-view, where the c-axis is normal to the
monolayer surface. During simulations, periodic boundary conditions are applied in all
directions. To prevent periodic images of the monolayer from interacting with one another, we
add a vacuum layer of thickness 10 Å in the z-direction.
We perform NAQMD simulations in both the canonical and micro-canonical ensembles
with a unit time step of 50 a.u. (~1.2 fs) and 8 a.u (~0.2 fs), respectively. Optical excitation is
modeled by manually promoting a fixed number of electrons from the top of the valence band to
the bottom of the conduction band. This is illustrated in in left-hand side of Figure 1(c), where
eight electrons are promoted from the four bands at the valence-band maximum to the four bands
at the conduction-band minimum at time t = 0 ps. The right-hand side of Figure 1(c) shows the
21
energies and occupation numbers of the electronic bands as a function of time. Bands in red
represent those which are occupied by two electrons, bands in blue are singly-occupied, while
black bands are unoccupied. At time t ~ 0.1 ps an electron from a doubly occupied band hops to
an unoccupied band, resulting in these bands each becoming singly occupied.
Figure 1. (a) Top and (b) side views of the simulated MoSe2 monolayer. (c) Schematic of the simulated
photoexcitation (left), and time evolution of the resulting electronic energies and occupations (right). Here, eight
electrons from the valence-band maximum are promoted to the conduction-band minimum. Red, blue, and black
lines denote electronic levels occupied by two electrons, single electron, and no electron, respectively. Electronic
levels with nontrivial dynamics are shown in thick lines.
To predict the photo-excitation-derived dynamics of suspended monolayers, in the
absence of energy dissipation to or interaction with a substrate or heat bath, NAQMD is first
carried out within the micro-canonical ensemble. In such simulations, the number of electrons,
the volume of the simulation cell, and the total energy of the system are conserved, while the
temperature of the system is allowed to fluctuate. We excite electrons from the top two valence
bands, which we denote as the valence band maximum (VBM) to the bottom two conduction
bands, which we denote as the conduction band minimum (CBM). In total, four electrons are
excited, corresponding to an excited charge carrier density of 1.19 × 10
#!
cm
-2
. Analogous
simulations were performed in the canonical ensemble to simulate the monolayer laying on a
22
substrate. In these simulations, the number of electrons, the volume of the simulation cell, and
the temperature of the system are conserved.
3.2.3 Static Simulations
In addition to ionic dynamics, the electronic structure of several excited states was
characterized through the electronic band structure, wavefunction iso-surfaces, charge densities,
and phonon dispersion curves. Excited-state Fermi surfaces were plotted as reciprocal-space iso-
surfaces at energy levels corresponding to representative charge carrier concentrations of 0.2, 1,
and 2 ⨉ 10
14
cm
-2
using the Wannier90 package [63]. Phonon dispersion curves were produced
using density functional theory (DFT) implemented in the Vienna Ab initio Simulation Package
(VASP) [64-66] to corroborate the prediction that the emergence of Fermi surface nesting leads
to phonon softening.
Static calculations were performed on a 6 × 6 × 1 supercell of a monolayer MoSe2
crystal structure, containing 108 atoms, in both its electronic ground state, as well as at two
different electronically excited states. The wavefunctions of valence electrons were constructed
with a plane wave basis set with a kinetic energy cutoff of 520 eV. Reciprocal space was
sampled using a 3 × 3 × 1 Monkhorst-Pack k-point mesh with a Gaussian smearing of orbital
occupancies of 0.1 eV. DFT perturbation theory within VASP was used to calculate the Hessian
matrix, which was fed into the open-source phonopy [67] package to generate the dispersion
relations for the normal vibration modes of the MoSe2 unit cell. The electronically excited states
were modeled using the ∆𝑆𝐶𝐹 method, where electrons are promoted from the highest-lying,
occupied KS energy levels in the ground state, to the lowest-lying, unoccupied energy levels. We
modeled the two different electronic excitations by moving 1 or 2 electrons from the valence
23
band to the conduction band, equating to an electron-hole concentration of 0.297 × 10
#!
cm
-2
and 1.19 × 10
#!
cm
-2
, respectively.
3.3 Results
3.3.1 Electronic Structure of Monolayer
Figure 2 shows the electronic band structure and partial density of states (DOS) for the d-
orbitals of Mo, which are the predominant orbitals making up the VBM and CBM, as well as
charge densities for the VBM and CBM. The band structure plot is color coded to show the
dominant d-orbital characters along the high symmetry lines in the first Brillouin zone.
We observe that the highest point in the VBM is the K-point (marked in Figure 2(a) with
a green, dashed circle), with majority 𝑑
8F
and 𝑑
8
!
/F
! character. The lowest point in the CBM is
also the K-point (marked in Figure 2(a) with a red, dashed circle), but is mainly of 𝑑
G
! character.
The iso-surfaces of the charge densities of the VBM and CBM, shown in Figure 2(b), agree well
with this orbital character assessment, with the CBM charge density showing the signature 𝑑
G
!
shape, and VBM showing the signature 𝑑
8F
and 𝑑
8
!
/F
! shapes. It can be seen that the excited
electrons, which are removed from the VBM and placed in the CBM, remain centered around the
Mo atoms, however, while they are formerly localized within the x-y plane, excitation causes
them to spread out along the z-direction. Since both VBM and CBM are at the K-point in
reciprocal space, the monolayer exhibits a direct band gap, in agreement with experiments [68].
It should be noted that MoSe2 in its bulk form has an indirect band gap, which highlights one
way in which these 2D monolayers of TMDCs differ from their 3D counterparts.
24
Figure 2. Electronic structure of monolayer MoSe2. (a) Band structure plot denoting by color the character of the
Mo d-orbitals that make up the majority of the state at each point, as well as a partial density of states of the
different d-orbitals of Mo. The top valence (conduction) band is denoted with a green (red), dashed circle, indicating
the state is predominantly composed of 𝑑
(
! (𝑑
)*
and 𝑑
)
!
+*
!) orbitals. (b) Charge densities for the top two valence
bands and the bottom two conduction bands, showing how the electronic charge density changes when we excite
electrons from their ground state.
3.3.2 Light-Driven Lattice Dynamics
To quantify lattice disorder along a given reciprocal vector q U⃗, we compute the Debye
Waller Factor (DWF) as a function of time t, given by the equation
𝐷𝑊𝐹(𝑡)= 𝑒
/〈|I J⃗∙M JJ⃗(9)|
!
〉
(5)
where 𝑢 U⃗(𝑡) is the displacement vector of an atom from its original position at t = 0, and the
brackets áñ denote an average over all atoms in the simulation cell. Figure 3 shows the time
evolution of the DWF for the monolayer in the micro-canonical ensemble, along two
25
representative directions in hexagonal 2D crystals – the zigzag direction (300) and the armchair
direction (110). The red line shows the DWF for the monolayer after electronic excitation, while
the blue line shows the DWF for an unexcited monolayer at 10K. The relative lack of atomic
motion at 10K results in a near-unity value for the DWF. DWF values for the room temperature
crystal are shown in Figure 4.
Figure 3. Time evolution of Debye−Waller factor along the (110) and (300) planes, with (red) and without (blue)
photoexcitation, averaged over three simulation runs beginning from identical positions but different initial
velocities in the micro-canonical ensemble. The shading shows the standard deviation across the three independent
simulations.
The simulation results show a sub-picosecond, exponential decay of the DWF along both
the (110) and (300) planes upon electron excitation, a result that is validated by experiment [49].
Furthermore, we observe a high-frequency modulation in the DWF for times t = 0.5 to 1.5 ps
after excitation. These high-frequency features, with a time period of 0.2 ps, correspond to the
emission of coherent phonons with ω = 2.5 THz, which falls within the range of softened
vibrational frequencies for acoustic modes at the M- and K-points in reciprocal space (see Figure
5(c)). We also observe that these high-frequency vibrations are quenched after ~1 ps which is
26
consistent with recent experimental work on photoexcited atomically-thin MoSe2, which
demonstrated ultrafast lattice thermalization [49].
Figure 4. Debye Waller Factors of the (110) and (300) planes in the monolayer MoSe2 crystal calculated at 300K.
This DWF value differs noticeably from 1 at this temperature, in contrast to the near-unity value observed at 10K.
The shading shows the standard deviation across the three independent simulations.
These results suggest that the underlying cause for such rapid onset of electronic-excitation-
induced lattice disorder is strong electron-phonon coupling. To investigate further, we compute
the electronic structure, specifically the electronic Fermi surface, of the photo-excited MoSe2
monolayer.
3.3.3 Evidence of Strong Electron-Phonon Coupling
Figure 5(a) shows the total DOS (blue line) for electrons in the monolayer TMDC
system, as well as the integrated DOS (red line) which is used to compute the electronic Fermi
surface at different electron-hole concentrations. Contours of these Fermi surfaces at varying
27
levels of electronic excitation are shown in Figure 5(b). We see that at minimal excitation, the
Fermi surface is localized around the K (1/3, 1/3, 0) point in the Brillouin zone. However, upon
stronger photoexcitation the Fermi surface develops pockets between the Γ (0, 0, 0) and K (1/3,
1/3, 0) points, usually denoted as the Σ point.
Figure 5. (a) Electronic density of states of MoSe2 monolayer near the band edge (blue line). Red line shows
integrated density of states, which is used to determine the Fermi level for a given value of electron−hole pair
density, n(e−h). (b) Fermi surface for excited carrier density, n(e−h), values ranging from 0 to 2×10
14
cm
−2
. The
Fermi surface is localized at the K-points at minimal excitation (red contours), while exposing Σ-pockets at higher
excited electron−hole densities (black and blue contours). The three nesting vectors 𝑞
,
8888⃗
-
, 𝑞
,
8888⃗
.
, and 𝑞
,
8888⃗
/
. correspond
to reciprocal vectors ΓM, ΓK, and ΓΣ, respectively. (c) Phonon dispersion curves in the ground state and two
increasingly excited states.
Similar behavior is observed for other monolayer materials in the TMDC family upon
electronic excitation or charge doping [69,70]. With these pockets come parallel (nested) sheets
of Fermi surface, and the Fermi surface nesting vectors which connect them. Such nesting
28
greatly increases the number of possible electronic transitions at the nesting wave vectors and
causes phonon softening at corresponding points in momentum space. Indeed, calculated
phonon dispersion curves for the monolayer of MoSe2 show phonon softening upon
photoexcitation.
Figure 5(c) shows the phonon dispersion curves for the ground state and two different
excited states along lines of high symmetry in reciprocal space. The phonon dispersion curve for
the ground state exhibits all positive frequencies, indicating the stability of the structure in this
state. However, with mild excitation, we begin to observe zone-edge phonon softening at the M
point in the Brillouin zone. Upon higher excitation, we see additional phonon softening at the K
point as well as a uniform decrease in atomic vibration frequencies throughout the Brillouin
zone, which reflects a weakening of interatomic bonds and leads to larger atomic displacements
(i.e. decreasing DWF values as observed in the NAQMD simulations). The onset of softened
phonon modes can be attributed to strong electron-phonon coupling in the material [71,72],
which in turn provides a pathway for extremely efficient, non-radiative relaxation of excited
carrier energy into lattice disorder. This provides an explanation for the ultrafast conversion of
electronic energy into atomic motion as observed by the sub-picosecond decay in the DWF.
3.3.4 Dynamics of Monolayer in a Heat Bath
While the NAQMD simulations in the micro-canonical ensemble most closely model
suspended monolayers, supported monolayers, which would be found in integrated devices are
more closely simulated with under the canonical ensemble. This is because the substrate on
which a monolayer would be placed in a device can affect the molecular dynamics through
modified electron-hole screening as well as acting as a weak heat bath for the excited monolayer.
29
Canonical ensemble NAQMD simulations, where the excited monolayer is coupled to a Nosé-
Hoover heat bath, can model the latter effect. Figure 6 shows the DWF versus time, in the (110)
and (300) directions at two different electron-hole concentrations in the canonical ensemble.
We observe that the large decay in the DWF, reflecting significant lattice disorder,
persists despite the coupling to a thermostat at 10K indicating that atomic motion in the
deformed electronic potential energy surface dominates lattice order on picosecond timescales
after photoexcitation. Furthermore, we see an increase in the amount of lattice disorder with
increased electron-hole density. This finding is validated by recent experimental work which
found similar dynamics [49].
Figure 6. Debye-Waller factor plotter versus time along the (110) (a) and (300) (b) planes, at two different electron-
hole densities, in the canonical ensemble.
3.3.5 Simulation Size Effect
It is important to note that dynamics resulting from softening of both the M-point and K-
point vibration modes can be observed only in supercells commensurate with the periodicity of
both vibration modes. The 6 ⨉ 6 ⨉ 1 supercell considered in this study is the smallest system
30
that is compatible with both these modes. Smaller simulation cells, like 4 ⨉ 4 ⨉ 1, show
qualitatively different atomic displacements upon electronic excitation due to their inability to
accommodate all softened vibration modes. This arises from the fact that in computer
simulation, continuous wavefunctions must be approximated by a collection of discrete points,
and thus, the size of the supercell can play an important role in the resultant observed physics.
The number of repetitions, N, of a unit cell with a lattice vector of length a in a given direction,
will result in k-point values that are multiples of 2π/Na being sampled. Thus, as we change N,
we change which k-points are sampled. If an important k-point is not sampled due to an
unsuitable choice of simulation cell size, important physics can be missed. In our particular
system of monolayer MoSe2, we find that the M-point (0.5, 0.0, 0.0), at the edge of the Brillouin
zone, plays a crucial role in electron dynamics due to phonon softening at that point. Whereas
the simulation cell with 108 atoms (6´6´1 unit cells) does sample the M-point, the simulation
cell with 48 atoms (4´4´1 unit cells) does not.
Figure 7 shows the DWF as a function of time, upon electronic excitation with a carrier
concentration of approximately 2.5 × 10
#!
cm
-2
for these two different sized simulation cells.
While both systems show an initial sub-picosecond decay of the DWF, the smaller simulation
cell experiences far less disorder, and even begins to restore order after about a half picosecond.
As this is inconsistent with experimental results, we attribute this phenomenon as an artifact of
simulation, which can perhaps be attributed to important k-points not being sampled.
31
Figure 7. Debye-Waller factor plotted versus time along the (110) and (300) planes, for a simulation cell with 48
atoms (red) and 108 atoms (blue) upon photoexcitation, each averaged over three simulation runs beginning from
identical position but different initial velocities in the canonical ensemble.
3.4 Discussion
In conclusion, NAQMD simulations were performed to model photo-excitation of
monolayer MoSe2 crystals in the microcanonical and canonical ensembles. Electronic structure
characterization revealed that electronic excitation leads to the formation of nested Fermi surface
states, which are connected by nesting vectors that correspond closely to wave vectors of
softened phonon modes, indicating strong electron-phonon coupling. Sub-picosecond decay of
the DWF was observed in both ensembles, indicating ultrafast transfer of optical energy to
thermal motion.
The rapid disordering of the MoSe2 lattice upon photo-excitation found in the NAQMD
simulations (and validated by experiment) has important implications for optical control of
structural phase of TMDC monolayers. Importantly, the NAQMD simulations provide the
32
motion of individual atoms in the monolayer upon photo-excitation, illuminating how probable
an optically triggered structural phase change is, if at all. While our initial simulations did not
observe a photo-excitation-driven change in structural phase, finely tuning the excitation energy
in an ensemble of simulations may reveal what level excitation is necessary to achieve such a
transformation. Furthermore, other simulations parameters may be adjusted, such as
temperature, strain, electric field, etc. to discover what conditions are necessary to induce a phase
change. In this way, NAQMD simulations provide a powerful tool for determining how to
optically control structural phase in TMDC monolayers for use in new phase-dependent
technologies.
33
Chapter 4: Machine-Learning-Aided Simulations
4.1 Introduction
Over the last decade, machine learning techniques applied to the material science domain
have shown great success in material screening and property prediction. For property prediction,
the property of interest is measured or computed for some percentage of a class of materials, and
a machine learning model is trained to predict the value of a property for the remaining
materials. When applied to quantum materials, the training data will usually come from ab initio
simulations based on DFT, which quickly become computationally expensive for even modestly
sized systems. Of the materials for which properties are known, a majority portion is collected
into a training data set, while those remaining make up test data set. A statistical model is built
based on the training data set and then used to predict the material properties of the structures in
the test data set.
Since the true values of the given material property have been computed for the test data
set, it is possible to quantify the accuracy the machine leaning model by comparing the true
values to the predicted values. In fact, quantifying this accuracy is essential for estimating how
the model will generalize to new structures. Parameters of the machine learning model can be
adjusted to achieve a high level of prediction accuracy, although fitting too many parameters too
closely can constrain the generalizability of the model. With careful selection of model
parameters and penalization of overfitting of the statistical model, the model can be used to make
high accuracy predictions of material properties of the remaining materials in the class of
structures for which no property calculations have been performed, thus bypassing expensive ab
initio computations. Success has been shown in predicting material properties such as band gap
[73-76], dielectric breakdown strength [77,78], melting point [74], and defects [79] using various
34
machine learning methods, such as support vector regression [73,74], neural networks [80,81],
and kernel ridge regression [76,82].
We focus on the class of materials consisting hetero-structures formed from monolayers
of semi-conducting TMDCs. Such hetero-structures are promising candidate materials for
components in optoelectric and thermoelectric devices. Good thermoelectric materials have
highly complex band structures on which their thermoelectric performance critically depends.
Stacking TMDC monolayers of different species creates such complex band structures.
Furthermore, the concept of reduced dimensionality to achieve high thermoelectric performance
has been very influential in thermoelectrics [50], but has proven difficult to achieve in practical
materials. These TMDC heterostructures, composed of 2D layers with reduced dimensionality
by definition, may therefore be a good fit for thermoelectric devices. The ability to engineer
stacks may provide a unique opportunity to finally realize the benefits of reduced dimensionality
to achieve high thermoelectric performance.
In this thesis, we use Gaussian Process Regression (GPR) in the space of vertically
stacked TMDCs for prediction of two distinct types of properties derived from electronic
structure. The first property type includes attributes of the band structure, which are critical
parameters for optoelectronic materials. Specifically, we use GPR to predict the band gap value
of structures as well as the conduction band minimum (CBM) and valence band maximum
(VBM) dispersion curves in momentum space. The second property type is a simplified proxy
for the thermoelectric figure of merit, a complex combination of electrical and thermal transport
properties that is challenging to optimize [83]. Specifically, we focus on the electronic transport
component of the thermoelectric figure of merit and use GPR to predict a recently proposed,
band structure dependent, electronic fitness function [84] (EFF) as a function of dopant
35
concentration for both n- and p-type doped hetero-structures. With a complex structure
dependence, this material property provides an excellent test case for our machine learning
algorithms.
In many instances, instead of the ability to predict a structure’s material properties, we
only need to screen a class of materials to find the structure that has an optimal value for a given
property, defined by the specific application. Here, optimality does not necessarily refer to the
maximum or minimum value of the property; it can also refer to a property value that is closest
to some desired value. For example, finding a structure with optimal band gap could refer to the
structure with the maximum band gap, or the structure with a band gap closest to, say, the
Shockley-Queisser limit for efficiency of solar cells [85]. Building a regression model to predict
the material property value of a structure, and then using this model to find the optimal material
by predicting the material property values for all structures, is neither an efficient nor scalable
process for solving this problem. This is because building the model can require many material
computations to be performed to create a sufficiently large training data set for accurate
prediction.
Instead, a machine learning model based on active learning is well suited for optimal
material discovery. In active learning, each material has an associated reward, computed by an
auxiliary reward function, which signifies the decrease in uncertainty associated with the model
if a particular data point is selected. The material property is calculated for the material with the
highest reward, and the data point is appended to the training data set in each iteration during
training. Here, we use a type of active learning called Bayesian Optimization (BO), a method
that attempts to discover the optimal point in black box functions [86-89] with minimal function
evaluations, to discover the structure with an optimal property value.
36
In this thesis, we use BO to search for either the structure with maximum band gap, or a
band gap closest to 1.1 eV, the Shockley-Queisser limit for efficiency of solar cells [85]. We
also apply BO to find the best thermoelectric material. Since every structure has a convex EFF
curve versus dopant concentration, each structure will have a peak EFF value at some dopant
concentration, which can will from structure to structure. We use BO to find the structure that
has the maximum peak EFF value, which we presume means the structure has optimal
thermoelectric performance. While we focus on band gap and thermoelectric EFF values as
representative material properties, these machine learning methods can be used to predict or find
optimal structures with respect to any other material property for which data is available.
4.2 Methods
The class of materials we consider in this thesis consists of hetero-structures formed from
monolayers of group VI TMDCs: MoS2, MoSe2, MoTe2, WS2, WSe2, and WTe2. We take this
set to be our alphabet Σ = {MoS2, MoSe2, MoTe2, WS2, WSe2, WTe2}, where each species of
TMDC becomes a distinct symbol. Hetero-structures are then represented by strings, 𝑤 =
𝑎
#
𝑎
"
...𝑎
O
, where the string length |𝑤| corresponds to the number of layers N in the stacked
hetero-structure. For example, a three-layer hetero-structure may be written as 𝑤 = MoSe2-
WTe2-MoS2, where |𝑤| =3. We define the set of all N-layered hetero-structures as 𝐻
O
=
{𝑎
O
= 𝑎
#
𝑎
"
...𝑎
O
| 𝑎
&
∈Σ; i=1,...,N}. As the number of layers N increases, the size of 𝐻
O
increases exponentially. Furthermore, the computation time required to calculate the electronic
properties for each hetero-structure using standard DFT calculations increases rapidly as 𝒪(𝑛
P
),
where n is the number of electrons in the hetero-structure. Thus, performing exhaustive
exploration of 𝐻
O
quickly becomes intractable for large N. We therefore turn to techniques from
37
machine learning to reduce the number of materials calculations necessary for discovery of
optimal materials.
Figure 8 shows the four main steps in our property and optimal-structure prediction: data
preparation, data computation, determination of numerical feature vectors to represent the
structures, and machine learning model development. Details of each step are described in the
subsections below.
Figure 8. Workflow for optimal structure and property prediction. First, structure files for a family of three-layered
materials are created and uploaded to the Materials Project (MP) database. Second, the MP infrastructure performs
all DFT calculations, and subsequently, transport calculations using BoltzTraP code are performed. A snapshot of
the material property data computed by MP database is pictured, along with the thermoelectric parameters computed
by BoltzTraP. Third, a numerical feature vector is assigned to uniquely represent each structure. Fourth, machine
learning models are designed based on the data to make predictions for either a material property or an optimal
structure.
4.2.1 Data Collection and Preparation
While the advantage of using machine learning is to eliminate the need to perform
exhaustive material property calculations to find the optimal structure for a given application,
here we perform exhaustive DFT calculations for the entire class of three-layered hetero-
structures, 𝐻
P
, in order to validate the our machine learning methods against the ground truth.
First, we create structure files for all unique three-layered hetero-structures and upload them to
the Materials Project [90] (MP) database. The MP framework then performs electronic structure
38
calculations based on DFT to obtain band structures for each material, which we use an input for
transport calculations using BoltzTraP [91] to get the thermoelectric EFF for each structure as a
function of carrier-dopant concentration. This function reflects complex band structures as they
relate to electronic transport functions relevant to thermoelectric performance. Specifically, the
thermoelectric EFF is large for electronic structures that decouple the normally inverse
relationship between electrical conductivity and thermopower. The electronic band structure and
the thermoelectric EFF for a sample structure, MoSe2-WS2-WS2, are given in Figure 9.
Figure 9. (a) Electronic band structure, showing the CBM dispersion curve highlighted in red, the VBM dispersion
curve highlighted in green, and the band gap denoted by the magenta arrow. (b) Thermoelectric EFF versus carrier-
dopant concentration for both n-type (green) and p-type (purple).
Figure 9(a) shows the electronic band structure, where the conduction band minimum
(CBM) and valence band maximum (VBM) dispersion curves are shown by red and green lines,
respectively, and the magenta double arrow indicates the band gap. Figure 9(b) shows the
thermoelectric EFF as a function of carrier-dopant concentration for both n- and p-type doping.
The curves are always convex, with a maximum thermoelectric EFF value occurring at different
carrier-dopant concentrations for the various structures.
39
4.2.2 Feature Vector Selection
For the machine learning models to learn how the electronic property data relate to their
corresponding structures, a numerical representation of the structure is required, as opposed to
the character strings we, as humans, use to identify the different structures (e.g. MoSe2-WTe2-
MoS2). Since the material design space is limited to stacked layers of TMDCs, many atomic
properties often used in feature vectors are either irrelevant or too similar across our materials to
be useful. A literature review of attributes specifically used for band gap prediction showed
those most commonly used are: electronegativity (EN), first ionization potential (IP), atomic
radius (AR), atomic number, atomic mass, period in the periodic table, valency, orbital radii,
formal ionic charge, melting point, and lattice constant [73-75,80,92,93]. Furthermore, a recent
study concluded the following set of elemental attributes were “broad enough to capture a
sufficiently-diverse range of physical/chemical properties to be used to create accurate models
for many materials problems” [94]: atomic number; Mendeleev number; atomic weight; melting
temperature; column and row in periodic table; covalent radius; EN, number of ‘s’, ‘p’, ‘d’, ‘f’,
and total valence electrons; number of ‘s’, ‘p’, ‘d’, ‘f’, and total unfilled states; specific volume;
band gap energy; magnetic moment; and space group number.
Since too many attributes in a feature vector can be detrimental in smaller-sized data sets,
we focus on a few attributes which are most relevant and show the most variability among the
elements. Attributes like melting temperature, magnetic moment, and band gap (the attribute we
were trying to predict!) are not relevant to band gap prediction, while attributes related to
statistics of valence electrons, position in periodic table, ionic charge, and space group are too
similar amongst our materials to be of any use in differentiating amongst them. Atomic mass,
atomic number, and AR are highly correlated attributes among our elements, so we chose to use
40
just AR to best represent the size of the elements. EN, i.e. the tendency of the atom to attract an
electron based on energy, and IP, i.e. the energy required to remove an electron from the atom,
are highly relevant attributes in predicting band gaps as they are related to the valence and
conduction bands, respectively, which in turn determine the band gap. Thus EN, IP, and AR
were chosen as a minimal set of most relevant attributes to use in our feature vectors.
Based on these three atomic attributes, feature vectors can be composed to represent each
hetero-structure by combining any subset of these atomic properties for each atomic species in
each layer into a vector. Values for the EN, IP, and AR for each of the five atomic species found
in the hetero-structures are given in Table 1. We normalize EN, IP, and AR values so that each
property has a distribution with zero mean and unit variance, a step which aids in the training of
accurate machine learning models. The numbers in the parentheses in Table 1 show the
normalized values. For three-layered hetero-structures, a feature vector using one of these
atomic properties, would be a six-dimensional vector since each layer contains two atomic
species. A feature vector using two of these atomic properties per layer would be a twelve-
dimensional vector.
Element EN [Pauling] IP [kJ/mol] AR [pm]
Mo 2.16 (-0.971) 684.3 (-1.479) 190 (1.152)
W 2.36 (0.051) 770.0 (-0.727) 193 (1.220)
S 2.58 (1.175) 999.6 (1.288) 88 (-1.159)
Se 2.55 (1.022) 941.0 (0.774) 103 (-0.834)
Te 2.10 (-1.277) 869.3 (0.145) 123 (-0.378)
Table 1. Electronegativity (EN), first ionization potential (IP) and atomic radius (AR) values for the five
elements of which the hetero-structures in this study are comprised, along with their normalized values
given in parentheses
To determine the best feature vector for our GPR and BO models for predicting band gap
and thermoelectric EFF values with high accuracy and minimum model complexity, various
models are built using different subsets of the three atomic properties (EN, IP, AR) to represent
41
the three-layered hetero-structures. In total, we tested the performance of 19 feature vectors,
whose specifications and dimensions are listed in Table 2.
Model
Number
Atomic Properties for Each
Atomic Species in Each Layer
Dimension of Feature Vector for a
Three-Layered Hetero-Structure
1 EN 6
2 IP 6
3 AR 6
4 EN,IP 12
5 IP,AR 12
6 EN,AR 12
7 EN,IP,AR 18
8 {Mo, W}=EN; {S, Se, Te}=IP 6
9 {Mo, W}=IP; {S, Se, Te}=EN 6
10 {Mo, W}=EN; {S, Se, Te}=AR 6
11 {Mo, W}=AR; {S, Se, Te}=EN 6
12 {Mo, W}=IP; {S, Se, Te}=AR 6
13 {Mo, W}=AR; {S, Se, Te}=IP 6
14 {layer 1,layer 3}=IP; layer 2=EN 6
15 {layer 1,layer 3}=EN; layer 2=IP 6
16 {layer 1,layer 3}=EN; layer 2=AR 6
17 {layer 1,layer 3}=AR; layer 2=EN 6
18 {layer 1,layer 3}=IP; layer 2=AR 6
19 {layer 1,layer 3}=AR; layer 2=IP 6
Table 2. Components and dimension of the various feature vectors tested.
For some feature vectors, different atomic properties are used for different classes of
atomic species. For example, “{Mo, W}=EN; {S, Se, Te}=IP” in Table 2 denotes that EN is
used to represent transition metals, while IP is used to represent chalcogens. We also introduce
“stack-dependent” feature vectors, in which the atomic property used to represent a given layer
depends on that layer’s position in the stack. For three-layered hetero-structures, layers 1 and 3
(i.e. outermost layers) are indistinguishable from each other due to mirror symmetry, whereas
layer 2 (inner layer) is distinct from them. Thus, a stack-dependent feature vector can assign one
atomic property to represent layers 1 and 3, while a different atomic property to represent layer
2. For example, “{layer 1,layer 3}=IP; layer 2=EN” in Table 2 denotes that IP is used to
represent all atomic species in the outermost layers (layers 1 and 3), while EN is used to
42
represent atoms in the internal layer (layer 2). We tested all 19 feature vectors for each machine
learning model we developed: (i) GPR for band gap prediction, (ii) BO for prediction of optimal
band gap material, and (iii) BO for prediction of optimal thermoelectric material. In general, we
found that “stack-dependent” feature vectors provided the best prediction accuracy for all target
material properties, where the atomic property used to represent a given layer depends on that
layer’s position in the stack
To determine the best feature vector for the GPR model for band gap prediction, we
created 100 independent machine learning models for each feature vector using training data set
sizes ranging from 40% to 70% of all three-layered hetero-structures. The trained model then
predicted the band gap for the remaining materials which did not appear in the training data set.
One figure of merit for the accuracy of the predictions is the mean-squared error (MSE), which is
given by 𝑀𝑆𝐸 =
#
O
∑ b𝑌
&,9.M,
−𝑌
&,Q.,=&R9,=
d
"
O
0
&S#
, where Ns is the number of structures for which
predictions were made. Defining “true” band gap values to be those computed by DFT, the
MSE for the predicted band gap values can be calculated. The average MSEs over the 100
independent runs for each feature vector and training data set size are reported in Table 3.
We observe that GPR models built using IP and AR as atomic features generally show
higher accuracy than models using EN, particularly when EN is used for the chalcogens. Since
the combination of IP and AR in the feature vector does not increase the model accuracy
significantly, simpler feature vectors, comprised solely of either IP or AR, are considered the
best candidates for a GPR model for band-gap prediction. In our final GPR model, each three-
layered hetero-structure is represented by a six-dimensional vector, in which atoms of layers 1
and 3 are represented by AR and atoms in layer 2 by IP.
43
Feature Vector
Training Data
(40%)
50 structures
Training Data
(60%)
75 structures
Training Data
(70%)
88 structures
EN 0.0353 0.0237 0.0204
IP 0.0302 0.0217 0.0175
AR 0.0296 0.0214 0.0175
EN,IP 0.0369 0.0237 0.0202
IP,AR 0.0300 0.0216 0.0175
EN,AR 0.0363 0.0235 0.0201
EN,IP,AR 0.0370 0.0240 0.0201
{Mo, W}=EN; S, Se, Te=IP 0.0302 0.0215 0.0174
{Mo, W}=IP; {S, Se, Te}=EN 0.0466 0.0402 0.0382
{Mo, W}=EN; {S, Se, Te}=AR 0.0295 0.0214 0.0175
{Mo, W}=AR; {S, Se, Te}=EN 0.0406 0.0258 0.0213
{Mo, W}=IP; {S, Se, Te}=AR 0.0295 0.0214 0.0176
{Mo, W}=AR; {S, Se, Te}=IP 0.0302 0.0217 0.0175
{layer 1,layer 3}=IP; layer 2= EN 0.0333 0.0220 0.0174
{layer 1,layer 3}=EN; layer 2=IP 0.0321 0.0229 0.0202
{layer 1,layer 3}=EN; layer 2=AR 0.0309 0.0230 0.0203
{layer 1,layer 3}=AR; layer 2=EN 0.0329 0.0219 0.0175
{layer 1,layer 3}=IP; layer 2=AR 0.0300 0.0216 0.0175
{layer 1,layer 3}=AR; layer 2=IP 0.0296 0.0215 0.0175
Table 3. Mean square errors on test data for various feature vectors, with varying sizes of training data set created
from the 126 unique 3-layered structures. In bold is the selected feature vector for our study.
To determine the best feature vector for our BO models, we performed 500 independent
BO runs for each target material property. For each run, 5% of structures were randomly
selected for the initial training data set, followed by 30 iterations of BO. Performance of the
different models for various target variables are given in Table 4. For maximum band-gap
prediction, the accuracy was highest when atomic species in each layer are represented solely by
IP or AR, or by combination of IP and AR. On the other hand, models using EN to represent
chalcogens tended to have lower accuracy. For maximum thermoelectric EFF value for p-type
or n-type doped hetero-structures, the model built using stack-dependent feature vectors
performed best. In the case of n-type doping, maximum accuracy for optimal structure search
was observed when atoms in layers 1 and 3 were represented by AR and those in layer 2 by EN.
For p-type doping, BO was best able to find the optimal structure when layer 1 and 3 atoms were
44
represented by EN and layer 2 atoms by IP. Finally, for n-type doped hetero-structures, models
using EN for chalcogens tended to have higher accuracy. The best feature vector for each target
property was used in the corresponding prediction models presented below.
Feature Vector Max. Band Gap (%) N-type EFF (%) P-type EFF (%)
EN 58.2 41.0 87.4
IP 77.0 48.4 82.4
AR 78.0 48.2 80.8
EN,IP 56.0 39.2 87.5
IP,AR 79.0 44.4 82.4
EN,AR 58.6 45.2 83.6
EN,IP,AR 57.8 38.4 86.4
{Mo, W}=EN; {S, Se, Te}=IP 74.0 44.8 81.8
{Mo, W}=IP; {S, Se, Te}=EN 58.8 42.4 85.0
{Mo, W}=EN; {S, Se, Te}=AR 76.4 47.8 82.6
{Mo, W}=AR; {S, Se, Te}=EN 58.2 45.2 86.0
{Mo, W}=IP; {S, Se, Te}=AR 79.0 48.8 81.4
{Mo, W}=AR; {S, Se, Te}=IP 78.4 48.0 79.8
{layer 1,layer 3}=IP; layer 2=EN 68.8 49.2 82.2
{layer 1,layer 3}=EN; layer 2=IP 70.8 41.0 88.0
{layer 1,layer 3}=EN; layer 2=AR 69.0 45.6 85.4
{layer 1,layer 3}=AR; layer 2=EN 73.4 54.2 83.6
{layer 1,layer 3}=IP; layer 2=AR 79.0 42.8 81.2
{layer 1,layer 3}=AR; layer 2=IP 77.2 43.6 81.6
Table 4. Percentages of the 500 BO runs in which the optimal structure is correctly identified for models built using
various feature vectors shown in column 1, for the various target properties given in columns 2-4. The best accuracy
for each target property is shown in boldface font.
4.2.3 Reference Data Preparation
Structure files were prepared for all unique three-layered hetero-structures 𝑤 ∈𝐻
P
.
When preparing configurations for three-layered hetero-structures, one may naively think there
would be 6
P
=216 configurations, since each layer can be one of six different TMDC species.
However, a structure with stacking sequence ABC (ABB) is the same as the structure with
stacking sequence CBA (BBA), as these two structures are simply 180° rotations of one another.
For all such pairs, we prepared a structure file only for the first one in lexicographical order but
not the other (i.e. MoSe2-WSe2-WS2 was kept in the data set, while WS2-WSe2-MoSe2 was
discarded). This resulted in 126 unique three-layered hetero-structures.
45
In order to construct a simulation cell for each three-layered hetero-structure 𝑤, we use a
supercell of WTe2 (as determined by experiment [95], found in the Inorganic Crystal Structure
Database [96]) as a template for atomic position and substitute atomic species as necessary to
create each specific configuration for all unique three-layered hetero-structures. WTe2 has the
largest unit cell of all the Group VI TMDCs and thus was chosen as a template to avoid strain in
the layered systems. A unit cell of WTe2 has two layers, so in order to produce a three-layered
structure, 1.5 unit cells of WTe2, containing 9 atoms, were merged into one supercell. Periodic
boundary conditions were applied in all directions, with a 10 Å band of vacuum added in the z-
direction (normal to the structure surface) to prevent multiple images of the structure from
artificially interacting. All structure files were generated automatically and then uploaded to the
MP database [90] using the pymatgen [97] python library.
After submission, all calculations, job-scheduling, parallelization, and workflow
management was handled on supercomputers by MP’s infrastructure. Electronic structure data
was calculated based on DFT [57,58], using a plane-wave basis set and the projector augmented
wave (PAW) method [61] as implemented in the Vienna Ab initio simulation package [64-
66,98]. The exchange and correlation energies were approximated with the Perdew-Burke-
Ernzerhof functional [62]. First, both the unit cell and the atoms were allowed to relax to the
lowest energy configuration. Next, self-consistent field iterations were performed to obtain the
single-electron wave functions, and the electronic band structure was computed from the
resulting eigen-energies of the wave functions. The wave functions were constructed from a sum
over a plane wave basis set, which consisted of plane waves with kinetic energies up to 520 eV.
The electronic band structures computed by MP database were then fed as input to BoltzTraP
[91] to compute transport properties used to compute the thermoelectric EFF.
46
4.2.4 Gaussian Process Regression
A Gaussian process (GP) is a collection of random variables, in which any finite number
of these variables has a joint Gaussian distribution. GP’s are used to describe a distribution over
functions and are completely specified by their mean function, 𝑚(𝑥), and covariance function
(or kernel) 𝑘(𝑥,𝑥
T
). Thus, a GP may be written as 𝑓(𝑥) ~ 𝐺𝑃b𝑚(𝑥),𝑘(𝑥,𝑥
T
)d. Gaussian
process regression (GPR) is a non-parametric regression technique, which models a distribution
of functions that are consistent with a given set of N training observations (𝑥
O
,𝑦
O
). In our case,
x is the n-dimensional input feature vector for each structure and y is either the structure’s band
gap, CBM or VBM dispersion curve, or thermoelectric EFF curve. We take an n-dimensional
squared exponential kernel as the covariance kernel, shown in equation 5.
𝑘(𝑥,𝑥
T
)=expo−
∑ V8
$
/8
1
$
V
!
2
$34
W
$
!
p 𝑖 =1,…,𝑛 (5)
Here, 𝜎
&
are hyper-parameters associated with each dimension of the feature vector and are
estimated using the maximum likelihood estimate. After training the model, we use the model to
make predictions on the test data set, which contains structures the machine learning model has
not seen.
4.2.5 Bayesian Optimization
Bayesian optimization is an optimization technique, which is generally used for hyper-
parameter tuning of deep neural networks and optimization of black box functions. Pseudocode
for the algorithm is shown in Table 5. In the BO process, a GPR model is first built using a
randomly selected, small percentage of structures as the initial training data set (much smaller
47
than the training data set size in the GPR models presented above). Since the true functional
form of the relationship between structure and material property is unknown, and since the GPR
model only provides crude predictions due to the small size of the training set, the procedure
optimizes a surrogate function, called the acquisition function [87], to determine which structure
to evaluate next.
Bayesian Optimization Algorithm
1. Data set: 𝐷 = {𝑥
&
,𝑦
&
} 𝑖 =1 to 𝑛
2. Build GPR model: 𝑦 =𝑓
:
(𝑥) ~ 𝐺𝑃(𝑚(𝑥),𝑘(𝑥,𝑥
T
))
3. Bayesian optimization () {
for t =1 to tmax
a. Select next xt via acquisition function (see below)
b. Compute the value yt for new xt with GPR model
c. Add (𝑥
9
,𝑦
9
) into data set 𝐷
d. Update the GPR model
}
Selecting next xt via Acquisition Function
1. Compute Acquisition Function 𝑢(𝑥)
𝑢(𝑥)= o
(𝜇(𝑥)−𝑓
*(8
)Φ(𝑍)+𝜎(𝑥)𝜙(𝑍) 𝑖𝑓 𝜎(𝑥) >0
0 𝑖𝑓 𝜎(𝑥)=0
where 𝑍 =
X(8)/Y
567
W(8)
; 𝜇(𝑥) and 𝜎(𝑥) are the mean and standard deviation
values for x predicted by the GPR model 𝑓
:
(𝑥), which takes maximum value
𝑓
*(8
; 𝜙(𝑍) and Φ(𝑍) are the probability density function and cumulative
density function, respectively, of the standard normal distribution.
2. Return x that maximizes 𝑢(𝑥)
Table 5. Pseudocode for the Bayesian optimization algorithm and details on use of acquisition function.
The acquisition function selects the next structure based on a trade-off between
exploration (to diversify the search) and exploitation (to follow the trend found by the current
model). Among the available acquisition functions, such as probability of improvement, upper
confidence bounds, and expected improvement, we used expected improvement (EI), the most
widely used acquisition function due to its simple analytical expression (see Table 5). The EI
value for the remaining uncalculated structures is computed, and the material properties of the
structure with the maximum EI are calculated next. This completes one iteration of BO. Each
48
newly computed structure is added to the training data set, and the next iteration begins with a
new GPR model built from the augmented training data set. The total number of iterations is
defined by the algorithm designer and constrained by how expensive it is to calculate new
structure data. Here, we use 30 iterations for each BO run, as this was sufficiently many to
predict the optimal structure with high probability, but few enough to remain computationally
feasible.
For the development of all our BO models, five structures were chosen at random to act
as the initial training data set. A total of 500 independent BO runs were carried out for each
target material property to gather statistics on how frequently the true optimal structure was
found, results of which are given in section 4.3. Distributions of how many iterations it took to
find one of the top five optimal structures for each target material property are shown in Figure
10 and Figure 11 .
Figure 10. Histogram of the number of iterations required to discover one of the five structures with the highest
band gap in the 500 independent BO runs.
49
A number of iterations equal to zero corresponds to the case where the initial five data
points contained one of the five best structures. In the remaining cases, we see that at most thirty
iterations are required to determine one of the five best structures.
Figure 11. Histogram of number of iterations required to discover one of the five structures with the highest peak
thermoelectric EFF value in the 500 independent BO runs.
4.3 Results
4.3.1 Gaussian Process Regression
We built GPR models to predict (i) the band gap, (ii) the VBM and CBM dispersion
curves, and (iii) the thermoelectric EFF curve as a function of carrier-dopant concentration for
both n-doped and p-doped hetero-structures. In the case of predicting the band gap, the target
variable Y is a scalar, whereas in the case of predicting the VBM, CBM and EFF curves, the
target variable 𝑌
U⃗
is a vector created by discretizing a continuous curve into discrete points. A
different model is trained for each point in order to create a predicted curve. Training each
model involves randomly selecting a percentage of the class of three-layered hetero-structures to
serve as each regression model’s training data set. The remaining structures, which the
regression model has never encountered, are used as the test data set. Each model is then tasked
50
with predicting its respective target variable (band gap, discrete points in the VBM, CBM, and
EFF curves) for all structures in the test data set.
To determine the appropriate training data set size, regression models were built using
different percentages of the total data set (126 unique structures) as the training data set, ranging
from 40% to 70%. Figure 12 shows the band gap predictions for the test data set (all remaining
structures not in the training data set) for two models built using 40% and 70% of the data in the
training data set. The model built using 40% has a larger confidence interval (CI) and the
predicted band gap for many structures lies outside the 95% CI. In contrast, the model made
from 70% is more robust, with smaller CI.
Figure 12. Effect of training set size on predicted band gap and its confidence interval on test data set for two
different training data set sizes.
To further test the robustness of the model, all three-layered structures were randomly
split into training and test data sets, with either 40%, 60%, or 70% training data, 100 times.
After building the model, the MSE for each model was used to calculate the prediction accuracy.
The average prediction MSE for each set of 100 runs for models built using 40%, 60%, and 70%
of structures in the training data set was 0.0295, 0.0214, and 0.0175, respectively. While the
MSE only slightly decreases, the width of the CI becomes much smaller with the larger training
51
data set. We also observe that increasing the training data set beyond 60% shows little
improvement in results. Hence, a training data set comprised of 60% of the total number of
structures is sufficient to build a model to predict band gap both with low MSE and narrow CI,
and the results from this model are shown in Figure 15(a).
Figure 13. Effect of training data set size on the predicted CBM and VBM for one of the three-layered hetero-
structures. (a) With 40% training set size, mean square error (MSE) between the predicted and ground truth is 0.177
and 0.153 for CBM and VBM, respectively. (b) With 70% training set size, MSE decreases to 0.085 and 0.067 for
CBM and VBM, respectively.
Figure 14. Effect of training dataset size on the predicted EFF function for one of the p-doped three-layered hetero-
structure. (a) Total training data set size used here is 40% of 126 structures. Mean square error between the ground
truth and predicted EFF is 2.78. (b) The same as (a) but with 70% training set size, mean square error on the same
structure decreases to 1.22.
52
Analogously we tested optimal training data set sizes for prediction of CBM, VBM, and
thermoelectric EFF curves. For each of these material properties, the target variable is a vector.
The vector consists of 30 discretized pairs of data points, where each pair includes either (i) a
particular wave vector and its corresponding energy value of the VBM or CBM or (ii) a
particular carrier concentration and its corresponding EFF value. A regression model is built
using 40%, 60%, or 70% of the data in the training data set at each of these 30 points and the
predicted output from these 30 models are used to construct the shape of the predicted curve.
Figure 13 shows the predicted VBM and CBM dispersion curves for one of the structures in test
data where the training set consisted of 40% or 70% of the structures. Figure 14 shows the
predicted thermoelectric EFF for one of the structures in test data where the training set consisted
of 40% or 70% of the structures for p-doped three-layered hetero-structures. As in the case of
band gap prediction, a training data set comprised of 60% of the hetero-structures proved optimal
for CBM, VBM, and EFF curve prediction, with the results shown in Figure 15(b-d).
In summary, for all predicted material properties, we found that training data sets with
fewer than 60% of structures did not produce reliable predictions, while training sets with more
than 60% did not show additional improvement. Resulting predicted values versus their ground
truth values are shown in Figure 15 for one instance of a GPR model for each target variable,
where models were trained with 60% of the structures in their training data sets. The ground
truth values were obtained using DFT and BoltzTraP code.
Figure 15(a) shows predicted and ground truth band gap values of the three-layered
hetero-structures in the test data set from one instance of the band gap prediction regression
model with a particular, randomly selected training data set. Figure 15(b) shows the predicted
53
and ground truth values for the VBM and CBM dispersion curves in momentum space for a
sample three-layered hetero-structure.
Figure 15. (a) Predicted (blue) and ground truth (red) band gap values of three-layer hetero-structures in the test
data set. Yellow shading represents a 95% confidence interval (CI) of the predicted results. (b) Predicted (dashed
lines) and ground truth (solid lines) values for the VBM and CBM dispersion curves in momentum space for a
sample three-layered hetero-structure. (c, d) Predicted (dashed lines) and ground truth (solid lines) thermoelectric
EFF curves versus carrier concentration for a sample n(p)-type-doped three-layered hetero-structure. In all models,
60% of the three-layered structures were randomly selected to comprise the training data set.
Figure 15(c) shows the predicted and ground truth EFF curves versus carrier
concentration for a sample n-doped three-layered hetero-structures, while Figure 15(d) shows
curves for a sample p-doped hetero-structure. As can be seen, it is possible to build different
GPR models to successfully predict a wide variety of material properties.
54
4.3.2 Bayesian Optimization
For applications where we need only find the structure with a desired value of some
property (especially in cases where computation of each structure’s material property is
expensive), BO can be used for efficient (i.e. with minimal structure calculations) discovery of
the optimal structure. In this thesis, BO models were first created to find either the structure with
the maximum band gap or the structure with a band gap value closest to 1.1 eV (Shockley-
Queisser limit for efficiency of solar cells [85]). Since each BO run begins by randomly
selecting a small number of structures to compute for the initial training data set, we performed
500 independent BO runs and collected statistics on the optimal structure returned by each of the
independent runs to average out any effects from the particular, randomly selected initial training
data set.
Figure 16 shows a scatter plot of the band gap values for the three-layered hetero-
structures, where structures are ordered along the x-axis by increasing band gap value. Structures
that were most frequently returned as the optimal structure are highlighted in red (for maximum
band gap) and green (for desired band gap). The corresponding pie charts show the percentages
of the 500 independent BO runs that the most frequently found optimal structures were returned.
For example, in the BO search for the maximum band gap, Figure 16(b) shows that 79% of BO
runs correctly found the structure with the maximum band gap of 1.7 eV, and 99% of BO runs
found one of the top three maximum band gap structures. Similarly good performance was
observed when searching for a desired band gap. While no structure existed with a band gap of
exactly 1.1 eV, 100% of BO runs found a structure with a band gap within 0.1 eV of the desired
band gap value.
55
Figure 16. (a) Band gap values for all three-layered hetero-structures, where hetero-structures on the x-axis are
ordered by increasing band gap value. Highlighted points denote hetero-structures returned most frequently as the
optimal structure by a BO model searching for the maximum band gap (red) and a desired value of 1.1 eV (green).
(b) Pie chart showing the distribution of optimal band gap values returned in 500 independent BO searches for the
maximum band gap. (c) Pie chart showing the distribution of optimal band gap found in 500 independent BO
searches for a desired band gap value of 1.1 eV. MoSe2-WSe2-WSe2, with a band gap of 1.05 eV, has the closest
band gap value to 1.1 eV and was returned as the optimal structure in 91% of the 500 independent BO searches.
WSe2-WSe2-WSe2 and MoS2-MoS2-WS2, with band gap values of 1.23 and 1.26 eV, respectively, were returned in
7% of the BO searches, while WSe2-MoSe2-WSe2 and MoSe2-WSe2-MoSe2, with band gap values of 1.01 and 1.04
eV, respectively, were returned in the remaining 2% of BO searches.
BO was also applied to the class of four-layered hetero-structures to test the method on a
class of layered structures for which it is difficult to compute materials properties for the full set.
Since the full set of four-layered materials was not computed, we cannot guarantee that the BO
process found the optimal band gap structure (in this case, the structure with maximum band
gap), but we did find that the model was consistently able to find a material with a higher band
gap value than the highest band gap value computed in the initial training data set, as shown in
56
Figure 17. Here, the initial training data set was comprised of 30 randomly selected structures
and 30 iterations were performed in each BO run. Since a stack-dependent feature vector with
atoms in the outermost layers represented by IP and atoms in the innermost layers represented by
AR had the highest accuracy during BO search of maximum band gap, we used an analogous
feature vector to represent four-layer structures: atoms in layers 1 and 4 were represented by IP
and those in layers 2 and 3 by AR.
Figure 17. Initial and final maximum band gap values for ten different BO runs on the class of four-layered hetero-
structures. Optimal structures, along with their band gap value are labeled
Finally, BO was used to discover the three-layered hetero-structure with the highest peak
EFF value for both n-type and p-type doping. Figure 18 shows the peak EFF values for p-type
and n-type doped three-layered hetero-structures. Here, the structures are sorted in ascending
order of peak EFF value along the x-axis, with structures most frequently returned as the optimal
structure in a set of 500 independent BO runs highlighted in different colors. Corresponding pie
charts in Figure 18, show the percentages of the 500 independent BO runs in which the most
57
frequently found optimal structures were returned. The two materials found as the best
candidates for n-type (p-type) thermoelectric devices were MoSe2-WS2-WS2 and WSe2-WTe2-
WSe2 (WTe2-MoTe2-WTe2 and MoSe2-WSe2-WSe2).
Figure 18. (a, b) Thermoelectric EFF values for all p-type-doped (a) and n-type-doped (b) three-layered hetero-
structures, where hetero-structures on the x-axis are ordered by increasing EFF value. Points highlighted in red and
green denote hetero-structures returned most frequently as the optimal structure by a BO model searching for the p-
type and n-type-doped structures, respectively, with maximum EFF value. (c, d) Pie chart showing the distribution
of optimal thermoelectric structures returned in 500 independent BO searches for the p-type-doped (c) and n- type-
doped (d) materials with maximum EFF value.
As mentioned previously, the thermoelectric EFF curve depends nontrivially on the band
structure of the material; specifically, it identifies band structures that decouple the Seebeck
coefficient, S, from the electrical conductivity, σ, in a way that favors high thermoelectric
performance. Simple parabolic bands do not have good peak thermoelectric EFF values, while
complex band features may enhance or degrade the thermoelectric EFF, depending on finer
58
details of the band structure. The band structures of the optimal thermoelectric materials are
shown in Figure 19.
Figure 19. Band structures along lines of high symmetry in the first Brillouin zone for (a) MoSe2-WS2-WS2 and (b)
WSe2-WTe2-WSe2, which are the two best n-type materials, and (c) WTe2-MoTe2-WTe2 and (d) MoSe2-WSe2-
WSe2, which are the two best p-type materials based on global optimization of the maximum electronic fitness
function. The red dots represent conduction band minima, while the green dots denote valence band maxima in all
four panels.
As can be seen in Figure 19(a), MoSe2-WS2-WS2 shows six carrier pockets at the
conduction band minimum, which occurs along the Γ-K direction. These pockets are
anisotropic, and additionally the second conduction band has a minimum relatively close in
energy, which can become active at higher temperature and doping levels. This is a highly
favorable situation for thermoelectric performance. Figure 19(b) shows the band structure of
WSe2-WTe2-WSe2 is qualitatively similar, except that the second conduction band minimum is
very close in energy to the first conduction band. This type of feature is generally favorable for
59
thermoelectric performance especially at lower temperature. Thus, the two materials identified
by our BO algorithm indeed show band-structure features favorable for n-type thermoelectric
performance.
Turning to p-type thermoelectrics, the band structure of WTe2-MoTe2-WTe2, shown in
Figure 19(c), has three bands that are very close in energy at the valence band maximum, which
is at the K point of the Brillouin zone, a favorable feature for thermoelectrics. The band
structure of MoSe2-WSe2-WSe2, shown in Figure 19(d), shows a similarly favorable feature at
the valence band maximum. Thus, in addition to their high values of the fitness function, the
identified compounds also show band structures qualitatively consistent with high thermoelectric
performance.
It is noteworthy that these optimal materials are not trivially similar to each other and that
the detailed band features that lead to the high thermoelectric EFF values are not the same in all
cases. For example, one of the n-type materials has two WS2 layers out of three, while the other
has no sulfide layers at all. The two high-performance p-type compounds have very different
band gaps. This points to a highly nontrivial dependence of peak EFF value on material
composition. Nonetheless, our BO approach was able to find materials with high thermoelectric
fitness in such a complex search space.
In order to test the effectiveness and generalizability of BO for optimal property
prediction, we tested our algorithm on a recently published data set [94] of adsorption energies
for eight different adsorbates (H, O, N, N2, CH, CO, HO, and NO) on twenty-five different
transition-metal surfaces. The data set is complete in the sense the adsorption energies are given
for all permutations for adsorbate/surface material pairs, for a total of two-hundred adsorption
energies. A recent study [94] examined twelve different descriptors for a large number of
60
different metals in Cu-alloy surface materials for adsorption energy prediction of CH4-related
species. Group number, surface energy and melting point of the metal were found to be the three
most important features. Therefore, we used these three descriptors as part of our feature vector
to describe the transition-metal surface material. To describe the adsorbate, we chose to use
standard entropy of the molecule as this was the most readily available thermo-chemical
descriptor available for all the adsorbate molecules used in this data set. In a search for the
minimum adsorption energy among the adsorbate/surface pairs, 82% of 500 independent BO
runs were able to successfully find the optimal pair after evaluating only 20% of the data set.
We therefore conclude that the BO method can be successfully used for the discovery of a
maximum, minimum, or desired value of a range of material properties.
4.4 Discussion
In conclusion, we showed that BO can successfully find a (nearly) optimal structure with
high probability and minimal materials calculations, for a wide range of material properties.
Within the class of three-layered TMDC hetero-structures, the BO method required the smallest
number of materials calculations in order to discovery the optimal structure: compared to the
GPR prediction method, BO required 30% fewer materials for the training data set; and
compared to exhaustive brute force search, BO required 70% fewer materials calculations to find
the optimal material. It is probable that for hetero-structures with more than three layers, these
percentages increase even more, since the increase in the number of BO iterations required for N
layers will increase more slowly than the size of the material class of N-layers, 𝐻
O
, for 𝑁 >3.
Our results show that classes of quantum materials can be screened for an optimal structure when
machine learning is used to guide which materials simulations should be performed. For some
61
classes of quantum materials, it would be impossible to find the optimal structure without the aid
of machine learning, since the number of simulations required would be intractable, particularly
for complex target material properties. By reducing the number of simulations required,
machine learning enables simulations to effectively screen quantum materials for optimal
structures. Such machine learning methods will be exceedingly useful as more classes of
quantum materials are discovered and developed, as they can help guide experimentalists to the
particular quantum material in a class that has a high probability of being most effective in new
technologies.
62
Chapter 5: Dynamic Simulations on Quantum Computers
5.1 Introduction
Quantum computers promise to solve certain classes of problems that are intractable on
the most advanced classical computers by massively reducing either time-to-solution,
computational resource requirements, or both. By storing and processing information on
quantum bits, or qubits, quantum computers take advantage of quantum mechanical phenomena,
such as superposition and entanglement, to beat performance of their classical counterparts.
Nearly forty years after their theoretical conception as universal quantum simulators [99-101],
quantum computers are beginning to produce early successes in simulating quantum systems.
Current and near-future quantum computers, dubbed Noisy Intermediate-Scale Quantum
(NISQ) [16] computers, comprised of 𝒪(10
0
-10
2
)
qubits, are beginning to be put to use in
scientific research. Pioneering work showed that simulating systems of fermionic particles on
quantum computers could be carried out with polynomial complexity [102] (as opposed to the
exponential complexity on classical computers). Subsequent work demonstrated how chemical
properties of such systems can be gleaned from the quantum computer [103,104], culminating in
the first experimental implementation of such a simulation on a quantum computer in 2010
[105]. However, NISQ computers suffer from high gate- and measurement-error rates, as well as
qubit decoherence [16]. Furthermore, due to their small numbers of qubits, NISQ computers are
unable to exploit robust error-correcting schemes which are expected to give rise to fault-tolerant
quantum computers in the future. As a result, the only physical systems that have been simulated
on quantum computers to date [105-109] are too small to see a quantum advantage (i.e. these
systems can still be simulated on classical computers).
63
The vast majority of simulations carried out on NISQ computers have been static
calculations [106,108-111] with time-independent Hamiltonians. However, much stands to be
learned from the dynamic simulation of systems governed by time-dependent Hamiltonians,
though studies of this kind are still in their infancy [112-114]. While such simulations have
provided encouraging proofs-of-concept, the next great challenge is to learn new physics by
simulating a quantum system on a quantum computer, which cannot feasibly be simulated on any
classical computer.
A particularly promising route to discovering new physics with NISQ computers is the
dynamic simulation of strongly correlated quantum systems. On the one hand, the complexity of
such simulations on classical computers grows exponentially, quickly making simulations of
even modestly sized systems impossible to run on state-of-the-art classical supercomputers. On
the other hand, only modestly sized systems are required to study the dynamics of various
physical models, making it possible to fit such simulations on NISQ computers. Together, this
means there may exist a “goldilocks” quantum system that is too large for classical computers to
handle, yet small enough for NISQ computers simulate, thereby achieving physical quantum
supremacy.
One class of promising quantum many-body dynamics problems to be addressed on
NISQ computers is the ultrafast control of emergent properties by electromagnetic radiation in
atomically thin layered materials (LMs). Functional LMs will dominate materials science in this
century [115]. The attractiveness of LMs lies not only in their outstanding electronic, optical,
magnetic and chemical properties, but also in the possibility of easily tuning these properties on
demand by an external stimulus like electromagnetic radiation [49,116,117]. Especially
promising is using terahertz (THz) radiation to directly excite specific phonon modes in the
64
material, which in turn modify atomic arrangement, amounting to ultrafast control of electronic
properties in the material [118].
Of particular interest is controlling magnetism in LMs [9,52]. In a recent experimental
study, small levels of emergent magnetism were observed in single-layer, Re-doped MoSe2, a
prototypical LM comprised of all nonmagnetic elements [9]. Separately, a theoretical study
predicted net magnetization in a similar nonmagnetic LM created by the controlled excitation of
a specific phonon mode [52]. Since THz photoexcitation of such LMs has been shown to excite
specific phonon modes within a picosecond [7,49], in principle, it could be used to control
magnetism in LMs on ultrafast timescales. The study of such dynamic magnetism in LMs is
inherently a quantum many-body problem for which near-future NISQ computers might provide
valuable insights.
Such dynamic simulations are carried out on quantum computers by creating a different
quantum circuit, or sequence of quantum logic gates carried out on each qubit, for each time-step
of the simulation. Unfortunately, current algorithms generally produce quantum circuits that
grow in size with each time step [113]. Since NISQ computers can only execute circuits up to a
certain size with high fidelity, due to high gate-error rates and low qubit decoherence times, this
puts a limit on the number of time-steps a NISQ computer can feasibly simulate. It is therefore
crucial to develop quantum circuit compilers that optimize circuits (i.e. minimize circuit size) for
dynamic simulations on NISQ computers.
Still in their infancy, quantum circuit compilers are an integral part of the quantum
computing software stack. Their function is to transform a high-level quantum circuit, usually
the output of an algorithm, into an executable quantum circuit that can be carried out on the
quantum computer. High-level quantum circuits are defined by a set of arbitrary quantum logic
65
gates, each of which can be represented by a unitary matrix, acting on different subsets of qubits
in the system. Quantum computers, however, are only designed with the ability to perform a
small, universal set of one- and two-qubits gates, called their native gate set. Products of these
native gates can be designed to be equivalent to any N-qubit gate, and thus can be used to carry
out any high-level quantum circuit. A complicating factor in the NISQ-era is that different
quantum computers have different native gate sets. Therefore, the executable quantum circuit
produced by the quantum circuit compiler must only consist of those native gates specific to the
quantum machine upon which it will be executed. On top of mapping high-level circuits to
native gate circuits, the quantum circuit compiler is further expected to optimize the native gate
circuit, which in the NISQ-era equates to circuit size minimization.
In this thesis, we present two sets of results related to dynamic simulations of transverse
field Ising models (TFIMs) on quantum computers. The TFIM is considered a quintessential
model for studying quantum phase transitions [119], as well as myriad condensed matter
systems, such as ferroelectrics [120] and magnetic spin glasses [121]. While dynamics of the
TFIM have an analytical solution in one dimension [122], studying the dynamics of the TFIM in
higher dimensions generally requires numerical methods that are extremely computationally
expensive. Simulation of a TFIM spin system on a NISQ computer, is therefore an excellent
candidate for the discovery of new physics.
In the first set of results, we use IBM’s Q16 Melbourne quantum computer and Rigetti’s
Aspen quantum computer to simulate a TFIM system, specifically a simplified model of Re-
doped, monolayer MoSe2 with a specific phonon mode excited. By measuring the average
magnetization of the system as a function of time, we demonstrate an early proof-of-concept of
quantum dynamic simulation on NISQ computers. In the second set of results, we present a
66
domain-specific quantum circuit compiler, specifically designed for compiling circuits that
simulate time evolution under a TFIM Hamiltonian. The special-purpose compiler is able to
reduce quantum circuit sizes by ~25% and wall-clock compilation times by ~40% compared to
the state-of-the-art general-purpose compiler (values vary depending on system size and
simulation time-step). We anticipate that this domain-specific compiler will allow for successful
simulations of larger TFIM systems on near-future NISQ computers that could not otherwise be
achieved with general-purpose compilers, due to the shorter circuit sizes it produces.
5.2 Methods
In Re-doped MoSe2 monolayer, Re atoms substituting Mo atoms have been
experimentally shown to form clusters [9]. Therefore, we consider a 1D cluster of Re atoms in
the monolayer as a simple yet nontrivial quantum many-body testbed, which is amenable for
study on currently available NISQ computers. A schematic of the material with a four-Re-atom
cluster is shown from the top view in Figure 20(a). Since each Re atom has one additional,
unpaired electron compared to the Mo atom it replaces, we map the spin of each of those
electrons to the spin of a qubit on the quantum computer, as depicted in Figure 20(b). We
describe the magnetism using the Ising model, where the exchange interaction strength Jz is
computed from first principles (see section 5.3 for details). We simulate excitation of the E²
phonon mode in monolayer MoSe2 (shown in the inset of Figure 20(a)), as it has been shown to
be the only mode that appreciably couples to the spin motion, i.e. affects the magnetism in the
monolayer [52]. The E² phonon mode gives rise to an effective magnetic field through spin-
67
orbit coupling, which can be incorporated into the Ising model as an oscillatory transverse
magnetic field with the E²-phonon frequency of 𝑓
Z[
=
\
9:
";
=4.8 THz [52].
Figure 20. Schematic showing how a simplified model of Re-doped monolayer MoSe2 is mapped onto the qubits of
a quantum computer. (a) Top-down view of the Re-doped MoSe2 monolayer where Mo, Se, and Re atoms are
depicted by pink, yellow, and purple spheres, respectively. Grey arrows are superposed on the Re atoms,
representing the spin of the extra, unpaired electron each Re atom possesses. Inset shows a side view of the material
and a representation of the excited 𝐸′′ phonon mode. (b) Spins of the extra, unpaired electron of each Re atom (grey
arrows), with exchange interactions between neighboring spins of strength Jz, in the presence of an external
magnetic field with frequency 𝜔
;<
and amplitude 𝜀
;<
are mapped onto the qubits of a quantum computer, shown in
their Bloch sphere representation.
The resulting time-dependent Hamiltonian is given by
𝐻(𝑡)= −𝐽
G
∑ 𝜎
&
G
𝜎
&]#
G O/#
&S#
−𝜀
Z[
cosb𝜔
Z[
𝑡d∑ 𝜎
&
8 O
&S#
(6)
where 𝜀
Z[
is the amplitude of the effective magnetic field produced by the excited E² phonon
mode, which is controlled by the fluence of incident electromagnetic radiation in the THz range.
In Eq. (6), 𝜎
&
^
is the 𝛼-Pauli matrix acting on qubit 𝑖. In our simulations, all parameters of
Hamiltonian (6) are held fixed, except for 𝜀
Z[
, which we vary over physically reasonable
strengths compared to the coupling strength Jz.
The simulated Hamiltonian is parameterized by three constants, namely 𝐽
G
,𝜀
Q_
,𝜔
Q_
, which
each have associated frequencies. In order to find a suitable ∆𝑡, we must find the fastest frequency
associated with the Hamiltonian and choose a ∆𝑡 an order-of-magnitude smaller than one cycle of
68
this frequency so as not to miss any dynamics. The frequency associated with 𝐽
G
is 𝐽
G
ℎ ⁄ =
2.86265 THz (where h is Planck’s constant). The frequency associated with the largest 𝜀
Q_
used
in our simulations is (𝜀
Q_
=5𝐽
G
) ℎ ⁄ = 14.3133 THz. Finally, the frequency associated with the
phonon is 𝑓 =𝜔
Q_
2𝜋 ⁄ = 4.8 THz. Thus, the highest frequency in our simulation Hamiltonian is
14.3133 THz, which has a period of 69.865 fs. We have chosen a time-step of 3 fs to adequately
resolve this fast oscillation. As a validation, we ran our simulations with a time-step of 3 fs, 1.5
fs and 0.75 fs and obtained identical results.
We simulate time evolution of our model under the time-dependent Hamiltonian (6) by
acting on the qubits the time-ordered exponential form of the unitary operator 𝑈(𝑡) ≡𝑈(0,𝑡)=
𝒯𝑒𝑥𝑝−𝑖∫ 𝐻(𝑡)𝑑𝑡
9
%
in the atomic unit. In order to map this unitary operator into a set of gates
in a quantum circuit, we first discretize time with a small time step of Δ𝑡, during which 𝐻(𝑡) can
be regarded constant [123]. We then apply Trotter decomposition [124] by splitting the
Hamiltonian into components that are each easily diagonalizable on their own: 𝐻(𝑡)=𝐻
8
(𝑡)+
𝐻
G
, where 𝐻
8
(𝑡)= −𝜀
Z[
cos(𝜔
Z[
𝑡)∑ 𝜎
&
8 O
&S#
and 𝐻
G
= −𝐽
G
∑ 𝜎
&
G
𝜎
&]#
G O/#
&S#
. Thus, the time
evolution operator is approximated as
𝑈(𝑛∆𝑡)=∏ 𝑒
/&`
7
(a4]
4
!
b∆9)∆9 :/#
4S%
𝑒
/&`
=
∆9
+ 𝒪(∆𝑡). (7)
We note that Hamiltonian (6) is in the form of a 1D Ising model with an oscillating, transverse
magnetic field. In the case of a static magnetic field, this model was solved exactly by Pfeuty
[122] following the methods of Lieb, Shultz and Mattis [125]. More recently, Ref. [126]
proposed applying this method to quantum simulations on quantum computers to efficiently
create strongly correlated quantum states and simulate their dynamical evolution for arbitrary
times, while Ref. [107] carried out this proposal on IBM’s quantum computer. While these
transformations only apply to strictly 1D chains, it is useful to be able to simulate the spin
69
dynamics of arbitrary clusters, for which no such transformation is known. We have thus
employed a more general solution method, Eq. (7), which applies to arbitrary clusters. The
following subsections provide details of the underlying techniques.
5.2.1 Quantum Circuits for Dynamic Simulation
The basic quantum circuit for simulating the dynamics of the system involves (i)
initialization of the qubits, (ii) application of the time-evolution operator to the qubits, and (iii)
measurement of each qubit in the computational basis, which is assumed to be the z-basis. An
illustration of a sample quantum circuit simulating evolution to the first time-step for a three-
qubit system is shown in Figure 21.
Figure 21. Quantum circuit diagram simulating time evolution to the first time-step for a three-qubit system.
The top three horizontal lines in Figure 21 represent the three qubits, while the bottom
line represents the three-bit classical register that stores the results of qubit measurement.
Moving from left to right in the diagram represents forward motion in time. Colored boxes on
top of the qubit wires represent different quantum gates acting on their respective qubit. The
blue box labeled H is the Hadamard gate; the green box labeled Rz is a rotation of the qubit about
the z-axis by a given angle; blue two-qubit gates represented by circles with a cross are CNOT
gates; and finally, magenta boxes represent measurements of qubits into the classical register.
70
Initializing qubits to a desired initial state is a non-trivial task, for which a number of
different methods have been proposed [17,102,103,127-132]. We use a ferromagnetic
configuration (all spins up) as our initial state, since it is the ground state of the Hamiltonian (1)
in the absence of THz radiation. The goal of the simulations is to study the dynamics of
magnetism by switching on the THz field at time t = 0. Fortunately, all qubits are initialized in
the spin-up position by default on both IBM’s and Rigetti’s quantum processors, and therefore
qubit initialization is trivial, i.e. no quantum gates are required. After qubit initialization, the
time-evolution operator given in Eq. (7) is translated into a set of quantum gates and applied to
the qubits, followed by measurement of all qubits. In our case, the measurement of interest is the
time-dependent, average magnetization of the N qubits along the z-direction, given by 〈𝑚
G
(𝑡)〉 ≡
#
O
∑ 𝜎
&
G
&
(𝑡).
Since measurement of the qubits destroys their quantum state, in order to simulate
dynamic evolution of the qubits through time, the qubits must always be initialized to their time
𝑡 =0 values before applying the appropriate time-evolution operator 𝑈(𝑛∆𝑡) (a separate
quantum circuit) for each time step n.
Set num_trials
For n in range(T):
Create circuit 𝐶
:
to simulate evolution to time step n with 𝑈(𝑛∆𝑡)
For i in range(num_trials):
Initialize qubits to t=0 state
Run circuit 𝐶
:
Measure qubits
Table 6. Pseudocode for dynamic simulation of qubits. For each simulated time-step, a different circuit is created.
Next, in a user-define number of independent trials, the qubits are initialized to their t=0 state, the particular circuit
for the given time step is executed, and the qubits are measured. Measurements from all trial runs are averaged
together to give an estimate for the final quantum state of the system at each time-step.
71
Furthermore, measurement does not give the full quantum state of the qubits; instead
each qubit will only return a 0 or 1. Therefore, the circuit for each time step must be run a large
number of times to reconstruct an estimate for the full quantum state. Table 6 shows pseudocode
for a dynamic simulation of a quantum system for a total time 𝑇∆𝑡.
5.2.2 Wavefunction Simulator
A wavefunction simulator uses a classical computer to simulate qubit evolution by
evaluating the many-body wavefunction of a system of N qubits (stored as a 2
O
× 2
O
matrix)
after a set of gates (unitary matrices) has been applied to the initial wavefunction of the system.
As such, unlike a real quantum computer, the wavefunction simulator provides access to the
entire many-body quantum state (quantum computers only allow measurement of qubits, which
collapses the full quantum state to a set of measured 0’s and 1’s). Furthermore, the wavefunction
simulator does not experience any effects from noise or decoherence. Due to the exponential
growth of the matrix storing the many-body wavefunction with the number of qubits in the
system, wavefunction simulators cannot be efficiently used to simulate large systems. Results
for small systems are nevertheless useful for validation purposes.
5.2.3 Simulated Noisy Qubits
Both IBM and Rigetti provide simulators of noisy qubits, which run on classical
computers, and like the wavefunction simulator, simulate qubit evolution on the quantum
computer by computing the many-body wavefunction. The main difference between the two is
how measurement is performed. The wavefunction simulator measures the expectation value of
an observable by calculating the inner product of the many-body wavefunction and the
72
observable operator. Measurement with simulated noisy qubits is performed by stochastically
collapsing each qubit to a 0 or a 1, according to their respective probabilities defined by the
quantum state of the qubit. Therefore, simulating noisy qubits, like the quantum computer,
requires averaging over an ensemble of measurements to approximate the expectation value of
an observable at each given time-step. Furthermore, noisy qubit simulators allow for the user to
set a custom noise model for the simulated qubits. The four main noise parameters that can be
set are the decoherence times 𝑇1 and 𝑇2; the readout error, defined by 𝑝(0|0) (the probability
that a qubit is read out as a ‘0’ given that it is in the ‘0’ state) and 𝑝(1|1) (the probability that a
qubit is read out as a ‘1’ given that it is in the ‘1’ state); and on IBM’s simulator one can define
the single-gate error rates. Readout errors are generally presumed to be symmetric, meaning the
𝑝(0|0)=𝑝(1|1). We set the noise parameters of our simulated noisy qubits to match those
given for the qubits at the time of running our simulations on the quantum processors.
5.2.4 Currently Available NISQ Quantum Computers
One of the many currently available quantum computers from IBM, the Q16 Melbourne
quantum processor has 14 qubits implemented with superconducting circuits. This quantum
processor can natively perform three quantum gates on the qubits: (1) a two-qubit CNOT (or CX)
gate, (2) a single-qubit RX(
;
"
) gate, which performs a rotation of the qubit about the x-axis by an
angle of
;
"
, and (3) a single-qubit RZ(theta) gate, which performs a rotation of the qubit about the
z-axis by an angle theta (theta can take on any value in the range [0,2𝜋]). Figure 22 shows the
topology of the quantum processor.
73
Figure 22. Topology of the IBM Q16 Melbourne quantum processor. Blue circles represent qubits, while black
lines represent physical connections between pairs of qubits.
The quantum processor can be accessed over the cloud with the use of IBM’s python
library Qiskit. Quantum circuits for the processor can be developed using high-level functions
provided by the Qiskit library. Qubits 6 and 8 were used for all simulations. All simulations
were run on September 28
th
, 2019. The decoherence and fidelity statistics for the quantum
processor at the time of running are shown in Table 7.
T1 [𝜇𝑠] T2 [𝜇𝑠] Readout Error One-Qubit Gate Error Two-Qubit Gate Error
Qubit 6 59.385 77.588 0.0351 0.0026776 0.03695
Qubit 8 44.600 66.480 0.0641 0.0058133 0.03695
Table 7. Decoherence and fidelity statistics for the qubits used on the IBM quantum processor.
The Rigetti Aspen quantum processor is comprised of 16 qubits implemented with
superconducting qubits. The quantum processor can natively perform three gates on the qubits:
(1) a two-qubit CZ gate, (2) a single-qubit RX(theta) gate, which performs a rotation of the qubit
about the x-axis by an angle theta (theta can only take on values of ±𝜋 and ±
;
"
), and (3) a
single-qubit RZ(theta) gate, which performs a rotation of the qubit about the z-axis by an angle
theta (theta can take on any value in the range [0,2𝜋]). Figure 23 shows the topology of the
quantum processor.
74
Figure 23. Topology of the Rigetti Aspen quantum processor. Blue circles represent qubits, while black lines
represent physical connections between pairs of qubits.
The quantum processor can be accessed over the cloud with the use of Rigetti’s python
library PyQuil. Quantum circuits for the processor can also be developed using high-level
functions provided by the PyQuil library. Qubits 0 and 1 were used for all simulations, which
were performed on June 11
th
, 2019. The decoherence and fidelity statistics for the quantum
processor at the time of running are shown in Table 8.
T1 [𝜇𝑠] 27.07
T2 [𝜇𝑠] 21.43
Readout Fidelity 95.05%
One-Qubit Gate Fidelity 91.93%
Two-Qubit Gate Fidelity 94.52% ± 0.23%
Table 8. Decoherence and fidelity statistics for the lattice used on the Rigetti quantum processor.
5.2.5 First Principles Calculation of Exchange Interaction Term
First principles calculations based on spin-polarized density functional theory (DFT) are
used to compute the exchange interaction term Jz in Hamiltonian (6). The system used is a
monolayer of ReSe2 and includes 108 atoms, corresponding to 6 × 6 × 1 unit cells in a box of
dimensions 19.728 × 19.728 × 16.2 Å
3
with periodic boundary conditions applied in all
directions, as shown in Figure 24. This box includes a vacuum of 12 Å along the perpendicular
75
direction to avoid spurious interactions between neighboring images of the monolayer. The
electronic states are calculated using the projector augmented-wave (PAW) method [61,98].
Projector functions are generated for 5d, 6s, and 6p states of rhenium (Re) and for the 4s, 4p, and
4d states of selenium (Se). The generalized gradient approximation (GGA) is used for the
exchange-correlation energy [62].
Figure 24. Schematic describing the calculation of the exchange interaction term Jz between unpaired electron spins
in neighboring Re atoms. Re atoms are depicted by cyan spheres, while Se atoms are depicted by yellow spheres.
The red block arrows indicate the unpaired electron of the Re atom is spin-up, while the blue block arrows denote
spin-down.
To more accurately assess the effect of on-site Coulomb repulsion among the localized
electrons in the partially filled d subshells, the DFT+U method [133] was employed with the
parameters 𝑈
dee
=1.5 eV [134] for Re 5d electrons. The van der Waals interaction between
atoms is described with the DFT-D approach [135]. The momentum-space formalism was
utilized, where the plane-wave cutoff energies are 30 and 250 Ry for the electronic pseudo-wave
functions and pseudo-charge density, respectively. The Γ point is used for Brillouin zone
sampling for electronic-structure calculations. We use our own scalable parallel QMD simulation
software, called QXMD [53], for the first principles calculation.
76
The local magnetic moments were obtained as |M
f
|=0.68 and 0.048 𝜇
g
/atom for Re
and Se, respectively. The exchange interaction parameter 𝐽
G
is obtained according to the
equation
𝐽
G
=
h
>?@
/h
?@
!X
AB
!
(8)
where 𝜇
ij
is a net magnetic moment per unit cell, 𝐸
kl
and 𝐸
mkl
are the energy per unit cell in
ferromagnetic and anti-ferromagnetic alignment, respectively [136]. Equation 8 gives a value of
𝐽
G
=0.01183898 eV for this material. The unit cells for the ferromagnetic system and
antiferromagnetic system are different, as shown in Figure 24, and also differ from the usual
smallest unit cell of ReSe2 when ignoring spin.
5.2.6 Compilation Methods
There are two main methods for circuit optimization. The first method, which we call the
minimal universal circuit (MUC) method, requires the construction of a MUC for each integer
number of qubits. The circuit is universal in the sense that it can implement arbitrary unitary
evolution of the qubits and is minimal in the sense that it uses the fewest number of gates to
implement the unitary evolution. Each MUC has an associated set of parameters whose values
must be solved for based on the action the corresponding high-level circuit is meant to execute.
The method works by representing all gates of the high-level circuit in their matrix
representation and multiplying them all into one unitary matrix. This matrix is then used to
determine the values of the parameters in the MUC. Note that while the values of the parameters
in the MUC change from time-step to time-step, the size of the MUC remains the same,
irrespective of how many gates were present in the high-level circuit.
77
Despite the advantage of constant-size circuits, the MUC method has two substantial
drawbacks. First, the minimal universal circuit itself must be created for each qubit number.
Thus far, this has only been studied for two and three qubits [137-141]. Second, for an N-qubit
circuit, the compiler must compute the 2
O
×2
O
unitary matrix representing the action of all
gates on the qubits from which to derive the parameters for the universal circuit, causing the
memory resources of this method to grow exponentially with system size. Together, these issues
render the MUC method unsuitable for compilation of the larger-scale simulations required for
new discoveries on NISQ computers.
The second method of compilation, which we call the transform and reduce (TR) method,
involves using algebraic identities to substitute subsets of gates in the circuit with alternate (e.g.
smaller, native) sets of gates. The method works by scanning through a circuit, searching for
sets of gates can be contracted together or transformed to a set of reduced size. While this
method does shrink circuit size for each time-step, the executable circuits still grow with time-
step. The advantage is that this method can be made to scale efficiently with system size.
However, the large number of identities that can be applied, along with the enormous number of
permutations in which they can be applied, makes optimizing this method an NP-hard problem
[142,143].
Since as few as O(1) gates can make a significant difference in the fidelity of a quantum
circuit executed on a NISQ computer, the importance of quantum circuit optimization is
paramount. Successfully discovering new physics with NISQ computers will require larger
numbers of qubits, making the minimal universal circuit compilation method untenable due to its
poor scaling with system size. Therefore, near-future NISQ simulations require the scalability of
the TR method. To address the difficulty of finding good circuit minimization, we propose
78
developing compilers for a specific native gate set and a specific class of quantum circuits.
Creating such domain-specific compilers allows the designer to take advantage of the structure
of the specific problem to develop heuristics that enable better and faster optimization over a
general-purpose compiler.
5.3 Results
5.3.1 Material Simulation on Quantum Computers
To establish a ground truth for validation of our quantum computer results, we first
calculate the dynamics of our system using a wavefunction simulator. A wavefunction simulator
simulates a noiseless quantum computer on a classical computer, and unlike the quantum
computer, provides access to the full quantum state of the system. We perform simulations for
2-, 3-, and 4-qubit systems with various values of 𝜀
Q_
, keeping all other values in Hamiltonian
(6) constant. The initial state is defined as all spins up, giving an initial average magnetization of
1. Since the exchange interaction strength 𝐽
G
is positive, the system is ferromagnetic and thus,
the exchange term (first term) of the Hamiltonian will tend to keep the qubits aligned. The
phonon-induced effective magnetic field term (second term) of the Hamiltonian, which acts in a
direction perpendicular to the initial orientation of the qubits, will tend to push them out of
alignment, reducing the average magnetization of the system. Thus, as we increase the ratio
𝜀
Z[
𝐽
G
⁄ , we expect to see larger drops in the average magnetization at the start of the simulations.
Figure 25 shows the resultant time-dependent average magnetization for 2- (red), 3-
(green), and 4-qubit (blue) systems, with varying values of 𝜀
Q_
, using the wavefunction
simulator.
79
Figure 25. Time evolution of the average magnetization of 2- (red), 3- (green), and 4-qubit (blue) systems with
electron-phonon coupling strengths 𝜀
;<
= 0.2𝐽
(
(a),0.5𝐽
(
(b), 𝐽
(
(c) and 5𝐽
(
(d). The black dotted line shows zero
magnetization.
Indeed, we see that as the strength of the electron-phonon coupling constant, 𝜀
Q_
, is increased
from 0.2𝐽
G
(Figure 25(a)), to 0.5𝐽
G
(Figure 25(b)), to 𝐽
G
(Figure 25(c)) the average magnetization
decreases by greater amounts. Another notable observation is the size effect across all values of
𝜀
Z[
, in which the more qubits the system has, the smaller the phonon-induced change in average
magnetization. We can attribute this to the fact that as the number of qubits in the system
increases, the ratio of bulk qubits to edge qubits also increases. In the simulated finite 1D
cluster, edge qubits (those which only have one nearest neighbor) only contribute one exchange
interaction term to the Hamiltonian, while center qubits (those which have two nearest
neighbors) contribute two exchange terms. As the ratio of bulk qubits to edge qubits increases,
the ratio of exchange interaction terms to transverse magnetic field terms in Hamiltonian also
increases. Therefore, increasing numbers of qubits reduce the effects of the transverse magnetic
field, and thus the average magnetization is reduced by smaller amounts in systems with more
qubits. When the magnetic field becomes much stronger than the coupling strength, 5𝐽
G
(Figure
25(d)), the qubits are initially all rapidly flipped, but then proceed to cycle at different rates
depending on how many qubits are present in the system.
We perform the same simulations on IBM’s Q16 Melbourne quantum processor and
Rigetti’s Aspen quantum processor. Figure 26 shows the results for a 2-qubit system with
80
varying ratios of strengths between the exchange interaction term and the transverse magnetic
field term, on the IBM (Figure 26(a-d)) and Rigetti (Figure 26(e-h)) quantum computers. The
number of trials for each time step in the simulations was 1,024 on the IBM quantum computer
and 10,000 on the Rigetti quantum computer. The entire simulation (all time-steps with 1,024 or
10,000 trials per time-step) was run five independent times on each quantum processor, with the
average results shown in red circles (a line is included between red circles to guide the eye). The
results from the quantum computers are compared with those from simulated noisy qubits (black,
dashed lines) as well as the results from the wavefunction simulator (black, solid lines).
Figure 26. Simulation results for a 2-qubits system (red dots) on the IBM quantum processor (a-d) and the Rigetti
Aspen quantum processor (e-h), compared to theoretical results from simulated noisy qubits (black, dashed lines)
and the wavefunction simulator (black, solid lines). The black dotted lines show zero average magnetization.
Results are shown for varying electron-phonon coupling strengths 𝜀
;<
= 0.2𝐽
(
(a,e),0.5𝐽
(
(b,f), 𝐽
(
(c,g) and 5𝐽
(
(d,h).
The first thing to note in Figure 26 is that while the qubits all start in the spin-up position,
the average magnetization at time 𝑡 =0 is not measured to be 1, as expected. This can be
attributed to readout noise, in which qubits that are in the state ‘0’ have a small probability of
being measured to be in the ‘1’ state, and vice versa. Readout noise will plague results in a
81
systematic way throughout the entire simulation and is one of the sources of error that can be
included in noise models for simulated qubits.
Looking at results across the rest of the simulation time, Figure 26 shows good
correspondence between those from the quantum processor and the wavefunction simulator,
which essentially simulates a noiseless quantum computer. However, NISQ computers will
always contain elements of noise, so it is more appropriate to compare quantum computer results
to those from simulated noisy qubits. Indeed, near overlap of results is found between the
quantum processor and simulated noisy qubits. The close correspondence of these results
provides compelling evidence that our current noise models are capturing the largest sources of
error on currently available NISQ computers.
The main reason such good results were produced for two-qubit systems is due to a
quantum circuit compilation technique that IBM’s and Rigetti’s general-purpose compilers use
specifically for two-qubit systems. For this small system, the general-purpose compilers of IBM
and Rigetti (with optimization turned on) use the MUC method of compilation (see Section
5.3.6). In this way, a short, constant-depth circuit can be defined for every time-step in the
simulation. When optimization is not turned on for IBM’s or Rigetti’s compiler, the high-level
circuits are naïvely translated to native gate circuits with no circuit size optimization, and thus
the executable circuits grow with each time step. The stark difference in results from simulations
run on Rigetti’s quantum computer using the naïve general-purpose two-qubit compiler (blue)
versus the optimized general-purpose two-qubit compiler (red) is shown in Figure 27.
Since the optimized, constant-depth circuits are short enough to produce high-fidelity
results, the NISQ computer is able to simulate the system out to arbitrary time steps, as can be
seen by the close correspondence of the red line to the black line, which represents numerical
82
results from a noiseless, simulated quantum computer. The naïve compiler initially produces
circuits that are short enough and thus can produce results with close correspondence to
numerical results. However, after a certain number of time-steps, it produces circuits that are too
long and thus generate too much error, leading to results that drastically diverge from the ground
truth. As all other factors are held constant between the two simulations except the compiler,
these results emphasize the importance of quantum circuit compilers to achieving high-fidelity
results on NISQ computers.
Figure 27. Results of a TFIM simulation run on Rigetti’s quantum computer, based on circuits produced by a naïve
compiler (blue) and an optimized compiler (red) for a two-qubit system. The black line represents numerical results
from a noiseless, simulated quantum computer.
Unfortunately, the MUC method used for two-qubit systems is not a scalable approach.
For larger systems, which are required to obtain physical quantum supremacy, a different, more
scalable method of quantum circuit compilation must be applied. Therefore, our domain-specific
compiler, which scales with number of simulation time-steps and system size, targets circuits
simulating evolution of systems with four or more qubits.
83
5.3.2 Domain-Specific Quantum Circuit Compiler
The domain-specific compiler was developed to compile quantum circuits that perform
evolution under the TFIM Hamiltonian, defined by
𝐻(𝑡) = −𝐽
G
∑ 𝜎
&
G
𝜎
&]#
G O/#
&S#
−𝐵(𝑡)∑ 𝜎
&
8 O
&S#
(9)
This Hamiltonian models a system of spins with nearest neighbor exchange interactions of
strength 𝐽
G
, in the presence of a transverse magnetic field with a time-dependent amplitude
defined by 𝐵(𝑡). The model can be mapped to the quantum computer by representing the spins
of the TFIM with the spins of the qubits. The dynamics of the system is simulated by applying a
sequence of quantum logic gates to the qubits such that time-evolution of the qubits emulates the
time-evolution of the spins under their Hamiltonian, given by the unitary operator 𝑈(𝑡) ≡
𝑈(0,𝑡) =𝒯𝑒𝑥𝑝−𝑖∫ 𝐻(𝑡)𝑑𝑡
9
%
in the atomic unit. Deriving this sequence of quantum logic
gates involves discretizing time and applying the Trotter approximation [124] to arrive at the
following approximate time-evolution operator:
𝑈(𝑛∆𝑡)=∏ 𝑒
/&`
7
(a4]
4
!
b∆9)∆9 :/#
4S%
𝑒
/&`
=
∆9
+ 𝒪(∆𝑡), (10)
where 𝐻
8
(𝑡)=−𝐵(𝑡)∑ 𝜎
&
8 O
&S#
and 𝐻
G
=−𝐽
G
∑ 𝜎
&
G
𝜎
&]#
G O/#
&S#
.
Dynamic simulations are performed by creating a different circuit for each time-step of
the total simulation, where the circuit for time-step n involves simulating evolution from time
𝑡 =0 to 𝑡 =𝑛∆𝑡. Each circuit, therefore, consists of the product of a number of operators
proportional to the number of time steps n. The general structure of such TFIM circuits is shown
in Figure 28.
Figure 28a shows a block of gates enacting time-evolution of the system by one time-step
∆𝑡. The blue gates are all single-qubit operators that represent applying the transverse magnetic
field to each qubit for a time ∆𝑡, while the green gates are two-qubit operators that represent the
84
exchange interaction amongst nearest neighbor qubits acting for a time ∆𝑡. Figure 28b shows a
high-level circuit diagram for simulating evolution of the system for a total time 𝑛∆𝑡. Note how
the block of gates shown in Figure 28a is repeated n times in Figure 28b, though the parameters
of each unit vary per time-step as our Hamiltonian is time-dependent. As such, a regular
structure appears in the final circuits for TFIM simulations, which can be exploited in the design
of heuristics for domain-specific quantum circuit compilers.
Figure 28. General structure of circuits for dynamic simulation. (a) This unit of unitary gates performs evolution
over one time-step of size ∆𝑡 of the spins under the TFIM Hamiltonian. (b) The unit of gates shown in (a) is
repeated a total of n times to evolve the system through n time steps. Notice how, for a time-dependent
Hamiltonian, parameters within each unit change from time-step to time-step.
Our domain-specific compiler, based on the scalable TR approach, uses gate substitutions
and compressions based on algebraic identities to compile TFIM circuits into native gate circuits
executable on IBM’s quantum computer. It scans through each circuit of the total simulation,
searching for where such identities can be applied and performs the identities in an order that has
85
been designed with knowledge of their general structure of the TFIM circuits. All algebraic
identities used for gate transformation and reduction are detailed in Error! Reference source
not found..
The left-hand column numbers the identities for identification. The middle column gives
sets of gates commonly found in uncompiled TFIM circuits, while the right-hand column shows
the corresponding equivalent to that set of gates. Identity 1 is used for converting the Hadamard
gate, a high-level gate, into a set of native gates. Note that while this does not reduce gate count,
it is essential for producing circuits that are executable on the quantum processors. Identity 2 is
used to prime circuits for further reductions. Finally, identities 3-5 are simplifying identities that
are used purely to reduce gate count. Note that if 𝜃
&
=−𝜃
4
in identities 4 and 5, then a rotation
by zero radians results, which can be completely removed from the circuit.
Table 9. Algebraic identities for the most used gate sets in TFIM circuits.
Not only is it important to implement specific identities useful for a given circuit type and
native gate set, it is also important to use knowledge of the circuit structure to develop heuristics
for optimal ordering of identity application. Ordering identity application in different
permutations can lead to varied resultant circuits.
86
Table 10 shows pseudocode detailing the order in which the identities from Error!
Reference source not found. are applied.
Table 10. Pseudocode detailing the order in which algebraic identities are applied for compilation of TFIM circuits
for execution on IBM’s quantum computer.
The first loop applies identity 2 to swap the control and target qubits of each CNOT gate, which
primes the circuit for further reductions. The second loop uses identity 3 to remove consecutive
pairs of Hadamard gates acting on the same qubit. The third and fourth loops makes use of
identities 1, 4, and 5 to transform all Hadamard gates into native gates and contract RZ gates on
opposite side of a control qubit of a CNOT gates. The fifth and sixth loops also make use to
identities 1 and 4 to convert remaining Hadamard gates to native gates and concatenate
87
consecutive rotation gates about the same axis on the same qubit. Finally, not shown in the
pseudocode, is the removal of any trailing gates at the very end of the circuit that only alter the
phase of the system, and thus do not affect a measurement gate that immediately follows.
To test the performance of the compiler, we constructed an algorithm to create a
sequence of high-level circuits that simulate evolution through time of a system of spins under a
TFIM Hamiltonian, as given in equation (9). The high-level circuits were fed into our domain-
specific compiler as well as IBM’s general-purpose compiler, and we compared performance
based on the number of native gates in the compiled circuits as well as the wall-clock time for
compilation for each time-step for systems with varying numbers of qubits. Results are
summarized in Figure 29. Figure 29a shows the absolute native gate count reduction when using
the domain-specific compiler over IBM's compiler. As shown, the gate reduction increases with
increasing time step for all system sizes, which implies that the compiler scales with simulation
time-step. Furthermore, for each time-step the number of reduced gates increases with additional
qubits, indicating our compiler scales with growing system size.
Figure 29. Performance comparison of the domain-specific (DS) compiler to IBM’s general-purpose compiler. (a)
Absolute gate count reduction using the DS compiler over IBM’s compiler. (b) Gate count reduction using the DS
compiler as a percentage of the gate count of IBM’s compiler. (c) Wall-clock compiler time reduction using the DS
compiler as a percentage of the wall-clock compiler time of IBM’s compiler.
Figure 29b shows the number of reduced gates as a percentage of the circuit size
produced by IBM’s compiler. The percent reduction in gate count is shown to asymptote with
88
growing time-step for all system sizes in a range around 25%. Figure 29c shows the percent
reduction in wall-clock compilation time of the domain-specific compiler compared to IBM's
compiler, which asymptote to values ranging from 20-60% depending on time-step and system
size.
5.4 Discussion
In this work, we presented results of dynamic simulations of emergent magnetism in a
simplified model of an LM performed on two state-of-the-art quantum computers. Remarkably
good quantitative agreement was found between the results from the quantum computer and
those from simulated noisy qubits, and to a slightly lesser extent with those from the
wavefunction simulator. Incorporating error models to account for the decoherence times of the
qubits, readout noise, and gate error-rates brought the results from simulated noisy qubits into
closer agreement with the quantum computer results, compared to those from the noise-free
wavefunction simulator. This indicates that there is a good understanding of the largest sources
of error we currently face on available NISQ computers. This early proof-of-concept indicates
that near-future NISQ computers, capable of simulating larger systems, may soon be able to
deepen understandings of controlled magnetism, as well as other controllable electronic
properties in LMs, for use in myriad new technologies.
We also presented a domain-specific compiler developed to compile circuits built to
simulate time evolution of spin-systems under a TFIM Hamiltonian. When comparing our
compiler to the general-purpose compiler provided by IBM, we found that it produced circuits
with fewer native gates and took far less time to complete compilation. The compiler scaled
with both number of simulation time-steps as well as number of particles in the system. This
89
implies that when higher-performance NISQ computers arrive, our compiler may able to create
short enough executable circuits to get high-fidelity results for simulations that would not
otherwise be possible with general-purpose compilers.
A few points should be mentioned about our compiler. First, note that the compiler takes
advantage of the fact that for these dynamic simulations, only measurement probabilities of
different states need to be conserved when applying transformation and compression circuit
identities. This means that the overall global phase of the system, which has no observable effect
on measurement probabilities, need not be conserved, allowing for a larger number of useful
circuit identities to be used. Second, this compiler assumes all two-qubit gates in the high-level
circuit are performed on two qubits that are physically connected on the quantum computer.
This leverages that fact that TFIM circuits only contain two-qubit gates that act on nearest
neighbor qubits, which makes it relatively simple to create high-level circuits that do not require
complicated qubit (re)mapping. For small TFIM systems, particularly in one-dimension, one can
usually find a lattice (set of qubits) on the quantum computer with a topology that is compatible
with two-qubit gates in the high-level circuit. Of course, this may not always be possible, in
which case, our compiler can be run first to reduce gate count, and a secondary compiler can be
performed to map the resultant circuit to the specific topology of a given lattice of qubits.
Finally, we note that the compiler is deterministic in its approach to translation and compression;
a given high-level circuit will be compiled to the same executable circuit every time. This is as
opposed to other available compilers which return different compiled circuits, of varying size,
for the same high-level circuit on different runs.
In the immediate future, NISQ computers with O(10
2
) qubits and optimized compilers
will allow straightforward extensions of the quantum simulations presented in this Chapter. One
90
possible avenue for future simulation could be the dynamic study of Floquet systems. Here, one
could simulate a time-periodic Hamiltonian 𝐻(𝑡) =𝐻
%
+𝐻
#
(𝑡) where 𝐻
#
(𝑡)= 𝐻
#
(𝑡+𝑇)
harmonically oscillates with a period of T, such as the Hamiltonian presented in this work. The
Floquet formalism [144] can be used to deal with such time-periodic Hamiltonians, which casts
them into effective, quasi-time-independent Hamiltonians, with associated quasi-stationary
states, known as Floquet states. Using the periodic perturbations as a tuning parameter, various
Floquet states can be created, leading to novel phases of matter and engineered materials with
properties not normally seen in equilibrium. Powered by quantum computers, dynamic
simulations of so called “Floquet engineering” [145], could help guide experiment to achieve
myriad kinds of on-demand functionality in novel metamaterials. Note that tuning of the 𝜀
Z[
parameter in our simulations could be considered simulated Floquet engineering of the system;
however, larger simulations requiring more stable qubits would be required to obtain novel
discoveries. An exciting possibility for near-future NISQ computers could therefore be to study
dynamic localization of electronic Floquet states [146].
The second immediate extension of this work could be the study of topological protection
of qubits in dissipative environments. A recent theoretical study suggested the use of nontrivial
topological phases associated with edge states in 1D spin-chains for protecting the coherence of
qubits in a dissipative environment [147]. Near-future NISQ computers, with the ability to
extend our system, will provide an ideal platform to explore such possibilities.
91
Chapter 6: Conclusion
This thesis has explored three mechanisms for controlling electronic properties in 2D
quantum materials. The 2D materials of interest are monolayers of semi-conducting TMDCs,
which possess numerous favorable electronic and material properties including strong chemical
stability; strong light-matter interactions; and tunable electronic, thermoelectric, magnetic, and
structural properties. Harnessing this tunability for integration of 2D TMDCs into new
electronic devices, however, requires an atomistic-level understanding of how these materials
respond to various tuning parameters. Tracking this response with the necessary spatial and
temporal accuracy is most readily achieved through computer simulation. In this thesis, three
state-of-the-art computer simulation techniques were applied for studying control of electronic
properties in 2D TMDCs.
First, we investigated the viability of optical control over the structural phase of 2D TMDCs.
2D TMDCs can exist in two different crystalline structural phases: one semi-conducting, and the
other semi-metallic. The ability to use ultrafast laser pulses to transform between these two
polytypes could allow for easier integration of 2D TMDCs into next-generation technology. The
first step of assessing such optical control over structure requires examining the atomic-level
dynamics of the monolayers upon optical excitation. Therefore, we performed NAQMD
simulations of a photoexcited TMDC monolayer. The NAQMD simulations allowed for the
observation of the time-dependent Debye-Waller factor, a quantitative measure of atomic
displacement along a certain direction. Results from the simulation, corroborated by experiment
[49], showed a sub-picosecond transfer of energy from the excited electrons into atomic lattice
motion. To determine the origin of the ultrafast transfer of electronic energy to the lattice, the
electronic Fermi surface and phonon dispersion curves were computed for varying levels of
92
electron excitation. Together, these provide evidence for strong electron-phonon coupling in the
system, which we concluded to be the origin of the ultrafast transfer of energy from excited
electrons to lattice disorder. Such strong lattice response to optical excitation indicates that light-
induced, ultrafast semiconductor-to-metal structural phase transitions is a possible route for
integrating 2D TMDCs into new quantum technologies.
Second, we investigated the control of thermoelectric properties in 2D TMDCs via stacking
order. In particular, we were interested in layered hetero-structures comprised of monolayers of
Group 6 TMDCs, which due to the 2D nature of their layers, are believed to be good candidate
materials for thermoelectric applications. The primary goal was to efficiently determine the
optimal thermoelectric hetero-structure using minimal quantum material simulations. As such,
we developed a machine learning algorithm based on Bayesian optimization that allowed us to
find the optimal material with high probability using as few quantum simulations as possible.
For validation purposes, we tested our algorithm on the class of three-layered hetero-structures.
After running the algorithm 500 independent times, we found that it was able to discovery one of
the top four optimal p-type doped materials in 98% of those runs, and one of the top five optimal
n-type doped materials in 97% of those runs. The success of the Bayesian optimization method
in discovering the best thermoelectric materials within the class of three-layered materials
indicates its potential for discovering optimal materials with respect to other material properties
from amongst different material classes. This ability to guide expensive quantum simulations
through machine learning opens the door for their use in materials screening and discovery.
Finally, we investigated the control of magnetism in 2D TMDCs via phonon excitation. We
simulated a simplified model of Re-doped monolayer MoSe2 on a quantum computer, showing
how varying the electron-phonon coupling in the monolayer could produce varied magnetization
93
of the system. Remarkably good simulation results were achieved for a two-qubit system, which
were attributed mainly to an optimized two-qubit quantum circuit compiler. Since this optimized
compilation method does not scale with larger system sizes, we developed a domain-specific
quantum circuit compiler, which produces circuits with fewer gates than those produced by a
general-purpose compiler. We anticipate that the shorter circuits created by the domain-specific
compiler can lead to better performance of simulations of control over magnetism in 2D TMDCs
on near-future NISQ computers. A better understanding of control over magnetism in such
materials paves the way for their integration into novel high-speed, spintronic devices.
Over a half-century ago, Moore’s law predicted exponential improvements in computing
power for the foreseeable future. The law held true for many decades as the miniaturization of
transistors allowed for exponentially more to be packed onto silicon chips every couple years.
However, such miniaturization has begun to reach its limits as transistors are becoming small
enough for quantum effects to become meaningful. A race has therefore begun to find alternate
technologies and materials that, instead of being impaired by quantum effects, use them to their
advantage. As 2D materials, monolayers of TMDCs have many favorable electronic properties
due to their quantum character. Furthermore, many of these properties can be controlled,
enabling material engineering of the monolayers for specific applications. This thesis explored
three such techniques for material property control, providing a fundamental understanding, and
thus, the eventual capability to engineer monolayers of TMDCs for integration into next-
generation quantum technology devices.
94
References
[1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V.
Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
[2] A. Gupta, T. Sakthivel, and S. Seal, Progress in Materials Science 73, 44 (2015).
[3] S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, Nature Reviews
Materials 2, 17033 (2017).
[4] G. H. Ahn et al., Nature Communications 8, 1 (2017).
[5] N. Lu, H. Guo, L. Wang, X. Wu, and X. C. Zeng, Nanoscale 6, 4566 (2014).
[6] M.-K. Joo, B. H. Moon, H. Ji, G. H. Han, H. Kim, G. Lee, S. C. Lim, D. Suh, and Y. H.
Lee, Nano Letters 16, 6383 (2016).
[7] L. Bassman, A. Krishnamoorthy, H. Kumazoe, M. Misawa, F. Shimojo, R. K. Kalia, A.
Nakano, and P. Vashishta, Nano Letters 18, 4653 (2018).
[8] L. Bassman et al., npj Computational Materials 4, 1 (2018).
[9] V. Kochat et al., Adv Mat 29, 1703754 (2017).
[10] K. Shimamura, F. Shimojo, R. K. Kalia, A. Nakano, K.-i. Nomura, and P. Vashishta,
Nano Letters 14, 4090 (2014).
[11] H. Kumazoe, A. Krishnamoorthy, L. Bassman, R. K. Kalia, A. Nakano, F. Shimojo, and
P. Vashishta, Journal of Physics: Condensed Matter 30, 32LT02 (2018).
[12] V. K. Sangwan, H.-S. Lee, H. Bergeron, I. Balla, M. E. Beck, K.-S. Chen, and M. C.
Hersam, Nature 554, 500 (2018).
[13] K. Rajan, Materials Today 8, 38 (2005).
[14] R. LeSar, Statistical Analysis and Data Mining: The ASA Data Science Journal 1, 372
(2009).
[15] T. Mueller, A. G. Kusne, and R. Ramprasad, Reviews in Computational Chemistry
(2015).
[16] J. Preskill, Quantum 2, 79 (2018).
[17] C. Zalka, in Proceedings of the Royal Society of London A: Mathematical, Physical and
Engineering Sciences (The Royal Society, 1998), pp. 313.
[18] N. Q. I. Act, in H.R. 6227, edited by HR2018).
95
[19] F. Arute et al., Nature 574, 505 (2019).
[20] V. Podzorov, M. Gershenson, C. Kloc, R. Zeis, and E. Bucher, Applied Physics Letters
84, 3301 (2004).
[21] H. Schmidt et al., Nano Letters 14, 1909 (2014).
[22] S. J. McDonnell and R. M. Wallace, Thin Solid Films 616, 482 (2016).
[23] X. Congxin and L. Jingbo, Journal of Semiconductors 37, 051001 (2016).
[24] Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman, and M. S. Strano, Nature
Nanotechnology 7, 699 (2012).
[25] M. Chhowalla, H. S. Shin, G. Eda, L.-J. Li, K. P. Loh, and H. Zhang, Nat Chem 5, 263
(2013).
[26] T. Cao et al., Nature Communications 3, 887 (2012).
[27] K. F. Mak, K. He, J. Shan, and T. F. Heinz, Nature Nanotechnology 7, 494 (2012).
[28] H. Zeng, J. Dai, W. Yao, D. Xiao, and X. Cui, Nature Nanotechnology 7, 490 (2012).
[29] E. M. Mannebach et al., ACS Nano 8, 10734 (2014).
[30] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Physical Review Letters 108, 196802
(2012).
[31] D. Voiry et al., Nat Mater 12, 850 (2013).
[32] G. Eda, T. Fujita, H. Yamaguchi, D. Voiry, M. Chen, and M. Chhowalla, ACS Nano 6,
7311 (2012).
[33] Y.-C. Lin, D. O. Dumcenco, Y.-S. Huang, and K. Suenaga, Nature Nanotechnology 9,
391 (2014).
[34] K.-A. N. Duerloo, Y. Li, and E. J. Reed, Nature Communications 5 (2014).
[35] L. Britnell et al., Science 340, 1311 (2013).
[36] M. Bernardi, M. Palummo, and J. C. Grossman, Nano letters 13, 3664 (2013).
[37] K. F. Mak and J. Shan, Nat Photon 10, 216 (2016).
[38] L. Cao, MRS Bulletin 40, 592 (2015).
[39] S. Bertolazzi, J. Brivio, and A. Kis, ACS Nano 5, 9703 (2011).
[40] J. Feng, X. Qian, C.-W. Huang, and J. Li, Nat Photon 6, 866 (2012).
96
[41] D. Yu, J. Feng, and J. Hone, MRS Bulletin 39, 157 (2014).
[42] A. Van Der Zande and J. Hone, Nat Photon 6, 804 (2012).
[43] J. Bang, S. Meng, Y.-Y. Sun, D. West, Z. Wang, F. Gao, and S. B. Zhang, Proceedings of
the National Academy of Sciences 110, 908 (2013).
[44] J. Zhang et al., Nature Communications 8, 683 (2017).
[45] G. Kolesov, D. Vinichenko, G. A. Tritsaris, C. M. Friend, and E. Kaxiras, The Journal of
Physical Chemistry Letters 6, 1624 (2015).
[46] C. R. Bealing and R. Ramprasad, The Journal of Chemical Physics 139, 174904 (2013).
[47] E. M. Mannebach et al., Nano Letters 15, 6889 (2015).
[48] S. Cho et al., Science 349, 625 (2015).
[49] M. F. Lin et al., Nature Communications 8, 1745 (2017).
[50] L. Hicks and M. S. Dresselhaus, Physical Review B 47, 12727 (1993).
[51] A. V. Kolobov and J. Tominaga, in Two-Dimensional Transition-Metal Dichalcogenides
(Springer, 2016), pp. 365.
[52] D. Shin, H. Hübener, U. De Giovannini, H. Jin, A. Rubio, and N. Park, Nature
Communications 9, 638 (2018).
[53] F. Shimojo et al., SoftwareX 10, 100307 (2019).
[54] F. Shimojo et al., The Journal of Chemical Physics 140 (2014).
[55] E. Runge and E. K. U. Gross, Physical Review Letters 52, 997 (1984).
[56] J. C. Tully, The Journal of Chemical Physics 93, 1061 (1990).
[57] W. Kohn and L. J. Sham, Phys Rev 140, A1133 (1965).
[58] P. Hohenberg and W. Kohn, Phys Rev 136, B864 (1964).
[59] N. Troullier and J. L. Martins, Physical Review B 43, 1993 (1991).
[60] D. Vanderbilt, Physical Review B 41, 7892 (1990).
[61] P. E. Blöchl, Physical Review B 50, 17953 (1994).
[62] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
97
[63] A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari,
Computer Physics Communications 185, 2309 (2014).
[64] G. Kresse and J. Hafner, Physical Review B 47, 558 (1993).
[65] G. Kresse and J. Furthmüller, Computational Materials Science 6, 15 (1996).
[66] G. Kresse and J. Furthmüller, Physical Review B - Condensed Matter and Materials
Physics 54, 11169 (1996).
[67] A. Togo and I. Tanaka, Scripta Materialia 108, 1 (2015).
[68] Y. Zhang et al., Nature Nanotechnology 9, 111 (2014).
[69] S. Zeng, Y. Zhao, G. Li, and J. Ni, Physical Review B 94, 024501 (2016).
[70] A. Krishnamoorthy, L. Bassman, R. Kalia, A. Nakano, F. Shimojo, and P. Vashishta,
Nanoscale (2017).
[71] P. B. Allen and M. L. Cohen, Physical Review Letters 29, 1593 (1972).
[72] J. E. Moussa and M. L. Cohen, Physical Review B 74, 094520 (2006).
[73] J. Lee, A. Seko, K. Shitara, K. Nakayama, and I. Tanaka, Physical Review B 93, 115104
(2016).
[74] T. Gu, W. Lu, X. Bao, and N. Chen, Solid state sciences 8, 129 (2006).
[75] G. Pilania, J. E. Gubernatis, and T. Lookman, Computational Materials Science 129, 156
(2017).
[76] A. Mannodi-Kanakkithodi, G. Pilania, T. D. Huan, T. Lookman, and R. Ramprasad,
Scientific Reports 6, 20952 (2016).
[77] C. Kim, G. Pilania, and R. Ramprasad, Chemistry of Materials 28, 1304 (2016).
[78] C. Kim, G. Pilania, and R. Ramprasad, The Journal of Physical Chemistry C 120, 14575
(2016).
[79] E. D. Cubuk, S. S. Schoenholz, J. M. Rieser, B. D. Malone, J. Rottler, D. J. Durian, E.
Kaxiras, and A. J. Liu, Physical Review Letters 114, 108001 (2015).
[80] Z. Zhaochun, P. Ruiwu, and C. Nianyi, Materials Science and Engineering: B 54, 149
(1998).
[81] E. D. Cubuk, B. D. Malone, B. Onat, A. Waterland, and E. Kaxiras, The Journal of
Chemical Physics 147, 024104 (2017).
98
[82] T. D. Huan, A. Mannodi-Kanakkithodi, and R. Ramprasad, Physical Review B 92,
014106 (2015).
[83] J. Yang et al., npj Computational Materials 2, 15015 (2016).
[84] G. Xing, J. Sun, Y. Li, X. Fan, W. Zheng, and D. J. Singh, Physical Review Materials 1,
065405 (2017).
[85] W. Shockley and H. J. Queisser, Journal of Applied Physics 32, 510 (1961).
[86] J. Mockus, Bayesian approach to global optimization: theory and applications (Springer
Science & Business Media, 2012), Vol. 37.
[87] J. Snoek, H. Larochelle, and R. P. Adams, in Advances in Neural Information Processing
Systems2012), pp. 2951.
[88] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. de Freitas, Proceedings of the
IEEE 104, 148 (2016).
[89] H. J. Kushner, Journal of Basic Engineering 86, 97 (1964).
[90] A. Jain et al., APL Materials 1, 011002 (2013).
[91] G. K. Madsen and D. J. Singh, Computer Physics Communications 175, 67 (2006).
[92] G. Pilania, A. Mannodi-Kanakkithodi, B. Uberuaga, R. Ramprasad, J. Gubernatis, and T.
Lookman, Scientific Reports 6, 19375 (2016).
[93] P. Dey, J. Bible, S. Datta, S. Broderick, J. Jasinski, M. Sunkara, M. Menon, and K. Rajan,
Computational Materials Science 83, 185 (2014).
[94] T. Toyao, K. Suzuki, S. Kikuchi, S. Takakusagi, K.-i. Shimizu, and I. Takigawa, The
Journal of Physical Chemistry C 122, 8315 (2018).
[95] A. Yanaki and V. Obolonchik, Inorganic Materials 9, 1855 (1973).
[96] G. Bergerhoff, R. Hundt, R. Sievers, and I. Brown, Journal of Chemical Information and
Computer Sciences 23, 66 (1983).
[97] S. P. Ong et al., Computational Materials Science 68, 314 (2013).
[98] G. Kresse and D. Joubert, Physical Review B 59, 1758 (1999).
[99] R. P. Feynman, Int J Theor Phys 21, 467 (1982).
[100] S. Lloyd, Science 273, 1073 (1996).
[101] D. S. Abrams and S. Lloyd, Physical Review Letters 79, 2586 (1997).
99
[102] G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, Physical Review A 64, 022319
(2001).
[103] A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Science 309, 1704
(2005).
[104] D. A. Lidar and H. Wang, Phys Rev E 59, 2429 (1999).
[105] B. P. Lanyon et al., Nat Chem 2, 106 (2010).
[106] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M.
Gambetta, Nature 549, 242 (2017).
[107] A. Cervera-Lierta, Quantum 2, 114 (2018).
[108] J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. E. Kimchi-Schwartz, J. R.
McClean, J. Carter, W. A. De Jong, and I. Siddiqi, Physical Review X 8, 011021 (2018).
[109] C. Hempel et al., arXiv preprint arXiv:1803.10238 (2018).
[110] P. J. J. O’malley et al., Physical Review X 6, 031007 (2016).
[111] J. Du, N. Xu, X. Peng, P. Wang, S. Wu, and D. Lu, Physical Review Letters 104, 030502
(2010).
[112] H. Lamm and S. Lawrence, Physical Review Letters 121, 170501 (2018).
[113] N. Wiebe, D. W. Berry, P. Høyer, and B. C. Sanders, Journal of Physics A: Mathematical
and Theoretical 44, 445308 (2011).
[114] A. Smith, M. S. Kim, F. Pollmann, and J. Knolle, npj Quantum Information 5, 1 (2019).
[115] A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013).
[116] D. N. Basov, R. D. Averitt, and D. Hsieh, Nat Mater 16, 1077 (2017).
[117] I. Tung et al., Nat Photon 13, 425 (2019).
[118] X. Li, T. Qiu, J. H. Zhang, E. Baldini, J. Lu, A. M. Rappe, and K. A. Nelson, Science
364, 1079 (2019).
[119] S. Suzuki, J.-i. Inoue, and B. K. Chakrabarti, Quantum Ising phases and transitions in
transverse Ising models (Springer, 2012), Vol. 862.
[120] R. Blinc, B. Žekš, J. F. Sampaio, A. S. T. Pires, and F. C. S. Barreto, Physical Review B
20, 1991 (1979).
[121] W. Wu, B. Ellman, T. F. Rosenbaum, G. Aeppli, and D. H. Reich, Physical Review
Letters 67, 2076 (1991).
100
[122] P. Pfeuty, Ann Phys 57, 79 (1970).
[123] D. Poulin, A. Qarry, R. Somma, and F. Verstraete, Physical Review Letters 106, 170501
(2011).
[124] H. F. Trotter, Proceedings of the American Mathematical Society 10, 545 (1959).
[125] E. Lieb, T. Schultz, and D. Mattis, Ann Phys 16, 407 (1961).
[126] F. Verstraete, J. I. Cirac, and J. I. Latorre, Physical Review A 79, 032316 (2009).
[127] P. Kaye and M. Mosca, in Optical Fiber Communication Conference and International
Conference on Quantum Information (Optical Society of America, Rochester, New York,
2001), p. PB28.
[128] A. N. Soklakov and R. Schack, Physical Review A 73, 012307 (2006).
[129] A. Kitaev and W. A. Webb, arXiv preprint arXiv:0801.0342 (2008).
[130] H. Wang, S. Ashhab, and F. Nori, Physical Review A 79, 042335 (2009).
[131] L. Veis and J. Pittner, J Chem Phys 140, 214111 (2014).
[132] N. M. Tubman et al., arXiv preprint arXiv:1809.05523 (2018).
[133] A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Physical Review B 52, R5467 (1995).
[134] A. L. Ivanovskii, T. I. Chupakhina, V. G. Zubkov, A. P. Tyutyunnik, V. N. Krasilnikov,
G. V. Bazuev, S. V. Okatov, and A. I. Lichtenstein, Physics Letters A 348, 66 (2005).
[135] S. Grimme, Journal of Computational Chemistry 27, 1787 (2006).
[136] M. Kan, J. Zhou, Q. Sun, Y. Kawazoe, and P. Jena, The Journal of Physical Chemistry
Letters 4, 3382 (2013).
[137] S. S. Bullock and I. L. Markov, pp. 324.
[138] V. V. Shende, I. L. Markov, and S. S. Bullock, Physical Review A 69, 062321 (2004).
[139] F. Vatan and C. Williams, Physical Review A 69, 032315 (2004).
[140] F. Vatan and C. P. Williams, arXiv preprint quant-ph/0401178 (2004).
[141] G. Vidal and C. M. Dawson, Physical Review A 69, 010301 (2004).
[142] A. Botea, A. Kishimoto, and R. Marinescu.
[143] D. Herr, F. Nori, and S. J. Devitt, npj Quantum Information 3, 35 (2017).
101
[144] G. Floquet, in Ann de l’Ecole Norm Sup1883), pp. 47.
[145] M. Bukov, L. D'Alessio, and A. Polkovnikov, Advances in Physics 64, 139 (2015).
[146] M. Holthaus, Physical Review Letters 69, 351 (1992).
[147] L. C. Venuti, Z. Z. Ma, H. Saleur, and S. Haas, Physical Review A 96, 053858 (2017).
Abstract (if available)
Abstract
Monolayers of semiconducting transition metal dichalcogenides (TMDC) are emerging as strong candidate materials for next-generation electronic devices. Interest in these two-dimensional (2D) quantum materials stems from the ability to control their structural and electronic properties via myriad schemes including mechanical strain, intercalation, carrier doping, optical and phonon excitation, combination into hetero-structures, and alloying. Since the integration of new materials into electronic devices often requires specific values for numerous material properties, broad control over 2D TMDCs makes them highly amenable to adaptation for new devices. For practical implementation, a deeper understanding of the mechanisms and capabilities of control over different material properties at the atomic level is essential. Achieving this level of accuracy with experimental study is difficult, often impossible. Therefore, computer simulation of the control over material properties in quantum materials has proven to be an invaluable tool. This thesis examines three methods of controlling different material properties in 2D TMDCs using three different simulation techniques. ❧ First, we investigate control over structural phase in 2D TMDCs via laser excitation using non-adiabatic quantum molecular dynamics simulations. Simulation of a monolayer after optical excitation showed a sub-picosecond disordering of the lattice, as measured by the Debye−Waller factor, as well as increasing disorder for higher densities of photogenerated electron−hole pairs. Furthermore, electronic structure analysis showed that the ultrafast, photo-induced lattice dynamics arise from strong electron-phonon coupling resulting from electronic Fermi surface nesting and its associated phonon-mode softening. Such strong lattice response to optical excitation indicates that optical control over structural phase in 2D TMDCs is a possibility. The mechanistic understanding uncovered with these simulations helps guide optical control for emerging applications such as phase change memories and neuromorphic computing. ❧ Second, we investigate control over thermoelectric fitness in hetero-structures, comprised of vertically stacked TMDC monolayers, via stacking order using machine-learning-aided simulations. We developed a machine learning algorithm based on Bayesian optimization that enables discovery of the optimal stacking order for thermoelectric fitness with high probability using minimal quantum simulations. Given the remarkable success of our algorithm for such a complex material property, we expect our algorithm to show similar success in finding optimal materials for various other applications. This ability to guide expensive quantum simulations with machine learning paves the way for discovering 2D TMDC hetero-structures for use in future devices, which may have vastly different material requirements from current norms. ❧ Finally, we investigate the control over emergent magnetism in 2D TMDCs via phonon excitation using quantum simulations performed on quantum computers. We demonstrate successful simulation of the control over magnetism by THz radiation of 2D TMDCs on IBM’s Q16 Melbourne quantum processor and Rigetti’s Aspen quantum processor. Remarkably good quantitative agreement was found between the results from the quantum computer and those from simulated noisy qubits for a two-qubit system. The quality of the results can be attributed to a special technique used by general-purpose quantum circuit compilers for two-qubit systems to create short, constant-depth circuits for all time-steps of a simulation. Since this compilation technique does not scale with system size, we developed a scalable, domain-specific compiler to aid in circuit size reduction for simulation of larger systems. We anticipate the domain-specific compiler will enable larger simulations (i.e. longer simulation times or larger system sizes) on near-future quantum computers that would not otherwise be possible. Together, our early proof-of-concept simulations and domain-specific compiler lay the groundwork for simulations performed on near-future quantum computers. Such simulations will deepen our understanding of controlling magnetism, as well as other controllable electronic properties, in 2D TMDCs for use in myriad future technologies. ❧ The last half-century has witnessed exponential improvements in electronic devices, due in main part to the miniaturization of their components. As we begin to reach the limits of such miniaturization the search for alternate technologies and materials for next-generation electronic devices becomes essential. The myriad controllable material properties of 2D TMDCs make them excellent candidate materials for new quantum technologies. This thesis demonstrates the viability of three techniques for controlling material properties in 2D TMDCs, paving the way towards their integration into post-silicon quantum technologies.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials
PDF
Multi-scale quantum dynamics and machine learning for next generation energy applications
PDF
Neural network for molecular dynamics simulation and design of 2D materials
PDF
Large-scale molecular dynamics simulations of nano-structured materials
PDF
Atomistic modeling of the mechanical properties of metallic glasses, nanoglasses, and their composites
PDF
Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations
PDF
Simulation and machine learning at exascale
PDF
Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations
PDF
Reactive and quantum molecular dynamics study of materials: from oxidation and amorphization to interfacial charge transfer
PDF
Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties
PDF
Molecular dynamics simulations of nanostructures
PDF
High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks
PDF
Synthesis, characterization, and device application of two-dimensional materials beyond graphene
PDF
In2O3 COVID-19 biosensors and two-dimensional materials: synthesis and applications
PDF
Phase change heterostructures for electronic and photonic applications
PDF
Thermal properties of silicon carbide and combustion mechanisms of aluminum nanoparticle
PDF
Optoelectronic, thermoelectric, and photocatalytic properties of low dimensional materials
PDF
From first principles to machine intelligence: explaining charge carrier behavior in inorganic materials
PDF
Molecular dynamics simulation study of initial protein unfolding induced by the photo-responsive surfactants, azoTAB
PDF
Theoretical studies of lipid bilayer electroporation using molecular dynamics simulations
Asset Metadata
Creator
Bassman, Lindsay Elizabeth
(author)
Core Title
Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras
School
College of Letters, Arts and Sciences
Degree
Doctor of Philosophy
Degree Program
Physics
Publication Date
04/01/2020
Defense Date
01/29/2020
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
Bayesian optimization,machine learning,materials discovery,materials simulations,non-adiabatic quantum molecular dynamics,OAI-PMH Harvest,quantum compiler,quantum computing,quantum materials,quantum simulations
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Nakano, Aiichiro (
committee chair
), Branicio, Paulo (
committee member
), El-Naggar, Mohamed (
committee member
), Kalia, Rajiv (
committee member
), Vashishta, Priya (
committee member
)
Creator Email
bassman@usc.edu,LINDSAYBASSMAN@GMAIL.COM
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c89-279403
Unique identifier
UC11673082
Identifier
etd-BassmanLin-8240.pdf (filename),usctheses-c89-279403 (legacy record id)
Legacy Identifier
etd-BassmanLin-8240.pdf
Dmrecord
279403
Document Type
Dissertation
Rights
Bassman, Lindsay Elizabeth
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
Bayesian optimization
machine learning
materials discovery
materials simulations
non-adiabatic quantum molecular dynamics
quantum compiler
quantum computing
quantum materials
quantum simulations