Modeling and Simulation of an

Oxygen Exchange Induced Resistive Switching

in a Two Layer Valence Change Memory Cell
Master Thesis

M.Sc. in Nanotechnologies for ICTs

Author: Supervisors:
Aria Sharifian Rainer Waser
s289756 Carlo Ricciardi
 Abstract

Some performances of the resistive switching mechanisms in Redox-based Random Access Memory
(ReRAM) Valence Change Mechanism (VCM) devices are listed to make it one of the candidates for
the future innovation of non-volatile memory application. There is a need for detailed and precise
explanations of their microscopic and nanoscopic behaviors.
Hence, in this thesis, features of the complex dynamics of bipolar resistive switching of ReRAM
devices, more specifically VCM resistive memory devices with HfOx/TaOx bilayer stacks, are brought
to light. A series of one-dimensional simulations with the finite element method is carried out using
COMSOL Multiphysics software to investigate the impact of applied external biases on internal prop-
erties such as oxygen vacancy concentration and conduction band energies. These are done by two
different approaches direct implementation differential equations and using built-in physical models
from COMSOL Multiphysics software.
From the results, it can be seen that there is a great change in the internal properties of the device
with the application of the pulsed voltages, and the changes depend on the magnitude, polarity, and
duration of the applied bias. In other words, the effects of magnification of the device’s conductivity
on voltages and duration, respectively, bring distinct resistance states to the device which may be
either a high resistance state or a low resistance state. More specifically, this can be seen through the
characteristics of the plotted device resistance versus the oxygen vacancy concentration. It is possible
to observe that the increase in oxygen vacancy concentration in the conducting oxide layer is always
accompanied by an increase in device resistance.
In this work it is tried to reproduce these findings using the Semiconductor Module of COMSOL
Multiphysics. Yet, there are certain problems in modeling the ions as charge carriers and the phe-
nomenon of direct tunneling. In order to prevent these problems, the possibility of enhancing the
methodology of simulation through the development of new advanced physical models and modules
is worked out.
It goes on to examine the effect of different electrode materials on device behavior in such a manner
that symmetric compositions of electrodes would cause symmetrical oxygen vacancy distributions
within the device. The material selection and their corresponding dependence on electrode metal
work functions further underscores the importance of material selection in device design.
Overall, this thesis sheds light on the complicated nature of VCM ReRAM devices, in the context

i
 ii

of simulation techniques and research in general, in order to exploit their potential to the hilt with
rising memory technologies.
 Acknowledgments

I would like to express my gratitude to everyone at RWTH-Aachen University’s IWE-2 Institute who
supported me during my thesis. Special thanks to Professor Rainer Waser, Dr. Stephan Menzel for
letting me join their team, and my daily supervisor Nils Sommer for his invaluable guidance and
support.

I also want to acknowledge the contributions of former students Astrid Marchewka and Fenja Berg,
whose work on the implementing models greatly benefited my research.

To my professors and friends at Politecnico di Torino, thank you for your support throughout my
master’s studies, specially Professor Carlo Ricciardi who helped me to find this great opportunity to
conduct my thesis abroad.

And finally, heartfelt thanks to my family for their unwavering support throughout my academic
journey.

Aria Sharifian - May, 2024 Aachen

iii
 Table of Contents

1 Introduction 1

2 Theoretical principles 3
2.1 General characteristics of resistive switching devices . . . . . . . . . . . . . . . . . . . 3
2.2 Valence change memory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2.1 Understanding VCM phenomena in crystalline structures using SrTiO3 model . 4
2.3 Charge transport principles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3.1 Schottky contacts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3.2 Ionic transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.3.3 Tunneling mechanism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3.4 Thermo-ionic emission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.4 ReRAM characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.4.1 Electro-forming of filament . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.4.2 SET and RESET in VCM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.4.3 Complementary switching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4.4 Bi-layer ReRAM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.5 Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.5.1 Hafnium oxide - HfOx . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.5.2 Tantalum oxide - TaOx . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3 Physical model for the simulation 13

3.1 Poisson equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1.1 General model equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1.2 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2 Electron continuity equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.2.1 General equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.2.2 Flux term (tunneling) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.2.3 Direct tunneling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2.4 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

1
TABLE OF CONTENTS 2

3.3 Ion continuity equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.4 Heat conduction equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

4 Simulations 23
4.1 Phase I-a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.1.1 Oxygen vacancy concentration . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.1.2 Conduction band energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.2 Phase I-b . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.2.1 Current-voltage characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.3 Phase I-c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.3.1 Oxygen vacancy concentration . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.3.2 Conduction band energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.4 Phase I-d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.4.1 Reaching relaxed state . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.4.2 Oxygen vacancy in different states . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.4.3 Resistance states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.5 Phase II-a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.5.1 Effect of applied voltage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.6 Phase II-b . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.6.1 Current-voltage characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.7 Phase III-a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.7.1 Reaching relaxed state . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.7.2 Oxygen vacancy concentration and resistance state . . . . . . . . . . . . . . . . 54
4.7.3 Effect of electrode metal material . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

5 Summary and outlook 60

5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

A Temperature fluctuation results 62

List of Figures 63

List of Tables 68

References 69
 Chapter 1

Introduction

Non-Volatile Memory (NVM) stands as a prevalent category of computer memories, capable of pre-
serving stored data even in the absence of power. High-density of cells, low fabrication costs, fast
read and write speeds, low energy consumption, and high write cyclability (i.e. high performance and
endurance) are among the desired characteristics of NVMs. Nowadays Silicon-based flash memory
devices are the most used NVM due to their high density and economical advantages. Nevertheless,
they face challenges such as limited endurance, low write speeds, and the requirement of high voltages
for writing operations [1].
Memristors (memory-resistor) or memristive devices can fulfill these requirements. They are highly
scale-able, fast in writing operations and low energy consuming in switching mechanisms. Two-
terminal memristor devices can retain their internal resistance state according to previously applied
voltages or currents. They are special in a sense that their behavior cannot be replicated using
conventional electronic circuit elements [2].
Memristors were initially predicted as the fourth passive electrical circuit element, which can relate
the electric charge with the magnetic flux [3]. A better description of memristors is devices with
pinched-hysteresis loop in voltage-current characteristics, the loop size is frequency dependant. The
natural application for these kind of devices is Resistive Random Access Memory (RRAM), although
other applications for example in alternative logic architectures can be considered for them [4].
The main idea behind redox-based random access memories (ReRAMs), is the fact that they can be
switched between two resistance states, i.e., High Resistance State (HRS or OFF) or Low Resistance
State (LRS or ON). Based on the state of the ReRAM memory cell, a bit of information is stored
through changes in the the dielectric material situated between the two metal electrodes [5]. As it can
be seen in the Figure 1.1, based on the resistive switching mechanisms, ReRAMs can be classified.

1
Introduction 2

Figure 1.1: Classification of the resistive switching effects which are considered for non-volatile
memory applications. From: [1], license number: 5787050663774.

Switching between a crystalline phase (ON-state) and an amorphous phase (OFF-state) is used in
phase change memories (PCM) [6]. In electro-chemical metallization (ECM), an electro-chemical
deposition and dissolution of a certain metal inside the cell is utilized to achieve resistance switching
[7]. For example, Ag ions can drift inside a ’I’ layer which can result in the growth of a highly
conductive filament of Ag, so the On state is achieved [8]. When the polarity of the applied voltage
is switched, a dissolution of the filament happens, so the cell is reset into the OFF state. Also there
are thermo-chemical mechanisms (TCM) which are based on the changes in stochiometry due to
temperature changes inside the material [1].
The mechanism that is in the focus of this thesis is Valence Change Memory effect. The main principle
of VCM can be seen as change in the valence state of the material in a resistive layer which changes the
conductivity of the cell [9]. Lately, there has been significant interest in bilayer resistive random access
memories due to their suggested superior endurance and performance in comparison to single-layer
devices [1].
In this thesis, we are focusing on bi-layer VCM structures in which two layers of resistive materials
are used. Lately, there has been significant interest in bilayer resistive random access memories due
to their suggested superior endurance and performance in comparison to single-layer devices [10]. To
gain a theoretical comprehension of the physical processes which are involved in resistive switching
bi-layer VCM devices, numerical simulations of these ReRAMs are becoming increasingly crucial [11].
A detailed numerical model to analyze VCM resistive switching is done by Astrid Marchewka [12]. In
this thesis the model is extended such as it enables to simulate bi-layer structures and analyze their
behavior.
 Chapter 2

Theoretical principles

This chapter serves the purpose of introducing the theoretical foundations which are essential to this
thesis. It explicates the fundamental principles underlying the used device.

2.1 General characteristics of resistive switching devices

The resistive switching phenomenon pertains to the reversible alteration of a material’s resistance
state under external electrical stimuli, transitioning between two or more stable states. This property
is harnessed in the development of re-writable nonvolatile memory devices, commonly known as resis-
tance (switching) random-access memories. These memories comprise numerous resistively switching
memory cells arranged in a structured memory matrix [12].
Broadly, a Resistive Random Access Memory (RRAM) device is composed of a Metal-Insulator-Metal
(MIM) structure, wherein the insulating or resistive layer is positioned between two conductive metal
layers [1].
Typically, the resistance states in this context are designated as the High Resistance State (HRS)
and Low Resistance State (LRS). The HRS is often referred to as the OFF state, while the LRS is
designated as the ON state. The transition from the High Resistance State to Low Resistance State
is termed ”Set”, whereas the transition from LRS to HRS is referred to as ”Reset”. These transitions
can be achieved by applying either a voltage or an electric current. Frequently, an electroforming step
is required in the cell before activating the switching mechanism [1].
In the context of voltage-driven resistive switching, the operational characteristics can be categorized
into two types: unipolar and bipolar switching. Unipolar switching occurs independently of the applied
voltage’s polarity. In this mode, the Set process necessitates a higher voltage compared to the Reset
process. Within the bipolar switching mode, distinct voltage polarities are requisite for the Set and
Reset processes. This asymmetry in voltage response arises due to a certain degree of asymmetry in
the device. The operational dynamics of these modes are illustrated schematically in Figure 2.1.

3
Theoretical principles 4

Figure 2.1: a) Illustration of unipolar switching, where both Set and Reset processes are unaf-
fected by the voltage polarity. b) Diagram depicting bipolar switching, where the Set process
occurs at one voltage polarity, and the Reset process takes place at the reverse polarity. Adapted
from: [1], license number: 5787050663774.

2.2 Valence change memory

In VCM, the metal-insulator-metal structure comprises an active electrode (AE) and an Ohmic counter
electrode (OE). The insulating layer is positioned between these electrodes functions as a mixed-ionic-
electronic conducting layer [13].
While a detailed comprehension of the underlying processes is currently not available, there is a
general consensus regarding resistive switching in VCM cells. This consensus suggests a reliance on
oxygen-vacancy migration and associated redox reactions [14].
Consequently, reduction and oxidation processes precipitate valence changes within the cation sublat-
tice, thereby justifying the naming ”valence-change memory”. Certain attributes of VCM structures
include their bipolar switching nature and the necessity for an electroforming step to facilitate the
switching process.

2.2.1 Understanding VCM phenomena in crystalline structures using SrTiO3

model
Valence change mechanism can be explained using a SrTiO3 model. SrTiO3 belongs to the Perovskite
class and features a cubic lattice structure. Within the unit cell, a titanium (Ti) atom is positioned
at the center, encompassed by eight strontium (Sr) atoms situated at the corners of the cell. The four
oxygen (O) atoms are strategically positioned at the face center in proximity to the titanium cations
[15], as depicted schematically in Figure 2.2.
In the case of SrTiO3 , the electrons originate from the last s-orbital and the inner d-orbital of titanium
atom participating in chemical bondings with oxygen. Consequently, titanium can release four elec-
trons for bonding, transforming into a fourfold charged positive ion, Ti4+ . This release of electrons is
named as oxidation. This oxidation results in a deficiency of conduction electrons available for current
transport, leading to insulating behavior in the oxide. However, in the presence of oxygen vacancies
within the structure of SrTiO3 , not all electrons from the outer shells of titanium are required for
oxygen bonding. Titanium undergoes a local reduction, changing to a double or threefold charged
ion, Ti2+ or Ti3+ . The liberated electrons in the conduction band contribute to the current transport,
thereby reducing the resistivity of the oxide. This dynamic interplay of reduction and oxidation by
oxygen induces a change in the valence state of the oxide [1].
Theoretical principles 5

Figure 2.2: Sketch of the lattice cubic structure of SrTiO3 .

Oxygen vacancies in the crystal lattice constitute point defects. Generally, point defects include vacan-
cies, interstitials, and substitutional effects. in its original state, how much the lattice is disordered
depends on different equilibrium situations involving defects. The Schottky equilibrium represents
one such equilibrium, describing the formation of anions and cation vacancies at elevated tempera-
tures. Additionally, the equilibria involving oxygen exchange between the lattice and the ambient
atmosphere, as well as the electronic equilibrium governing the recombination of electron-hole pairs,
significantly influence the lattice disorder [13]. Further elaboration is available in [13].

2.3 Charge transport principles

This section introduces and provides a brief explanation of the electronic and ionic transportation
fundamentals employed in the study of Valence Change Memory (VCM) within this thesis.

2.3.1 Schottky contacts

When a metal and a semiconductor come into contact and reach thermal equilibrium, their Fermi
energy levels become equal. Achieving this equilibrium involves adjusting the semiconductor’s energy
level by an amount known as the barrier height, which is equal to the difference between the two
work functions. The work function, denoted as qWM for metals and qχcs + qϕn for semiconductors,
represents the energy gap between the vacuum level and the Fermi energy level. Here, qχcs is the
electron affinity measured from the bottom of the conduction band (EC ) to the vacuum level, and
the energy difference between EC and the Fermi level is qϕn . In these expressions, q represents the
elementary positive charge. Therefore, the barrier height can be expressed as [16]:

qϕBn0 = q(WM − χsc ) (2.1)

A graphical representation (band diagram) of a Schottky barrier formed between a metal and a
semiconductor is depicted in Figure 2.3.
Theoretical principles 6

Figure 2.3: Energy-band diagram of a metal/n-type semiconductor contact.

Nevertheless, the accuracy of the Equation 2.1 is not always upheld in experiments. Real-world
scenarios often involve a physical separation between the two materials, potential presence of interface
states, and the possibility of image-force lowering affecting the barrier [16].

2.3.2 Ionic transport

Given that the switching behavior of Valence Change Memory devices is determined by oxygen vacan-
cies, understanding their transport characteristics is crucial. The ion transport mechanisms relevant
to VCM devices contains thermo-diffusion, drift, and diffusion. In essence, ion transport relies on
hopping within an energy landscape, therefore, overcoming an activation energy (WA ) is necessary.
In the absence of an external field, the ion hopping is due to isotropic diffusion movement. Conversely,
the presence of an external electrical field (E) arising from an applied voltage (V ) induces a reduction
in the energy barrier along the field direction by −e · a · E2 , where e represents the electron charge and
a is the hopping distance. This field-induced modification facilitates ion movement in the direction of
the field, leading to a drift motion. A sketch of the ionic transport in presences and the absence of
the external force is depicted in Figure 2.4.
Theoretical principles 7

Figure 2.4: Schematic of the potential energy landscape for ion hopping at a) absence of external
force b) when applying voltage . Redrawn form [17] © 2012 IEEE.

The comprehensive ion migration is understood through the superposition of all transport modes.
In the context of Valence Change Memory devices, this implies that the migration of oxygen vacan-
cies within the crystal lattice results in a localized alteration of the oxide stoichiometry. Regions
enriched with oxygen vacancies consequently exhibit an increased availability of conduction electrons
for electronic transport [17].

2.3.3 Tunneling mechanism

In classical physics, the transmission or reflection of electrons when encountering an energy barrier is
determined by the energy of the electron and the energy level of the barrier. If the electron possesses
greater energy, it will be transmitted, while lower energy results in reflection. Quantum mechanics
introduces the phenomenon of tunneling, wherein electrons may traverse a barrier with energy lower
than the barrier itself [11]. The extent of tunneling is influenced by various factors. Figure 2.5
illustrates two examples comparing classical mechanics and quantum mechanics for tunneling through
a barrier.

Figure 2.5: Transmission probability in different energy levels of the electrons in an energy
barrier based on the classical and quantum mechanics equations. On left and on the right the
amount of the transmission is different. This is due to the fact that thickness of the barrier is
changed.

A mathematical model for the tunneling process at the interface of the connection between two
different materials is provided by Tsu and Esaki [18]. In this model, the tunneling through the barrier
Theoretical principles 8

is expressed as:
Z Ec,max
4πmeff q
Jtunnel = T C(Ez )Nsup (Ez )dEz (2.2)
h3 Ec,min

In the provided expression, h represents Planck’s constant, meff denotes the effective mass of the
particle, and q stands for the elementary charge. The transmission coefficient (T C), corresponding to
the transmission probability, is determined by the electron’s energy levels (Ez ), and Nsup represents
the supply function for that specific energy level. The terms Ec,min and Ec,max denote the minimum
and maximum values of the energy barrier at the interface. Detailed equations for the transmission
coefficient and supply function are elaborated in Chapter 3.

2.3.4 Thermo-ionic emission

At different temperatures, electrons may occupy higher energy levels. Thermally excited electrons can
possess energy greater than the potential barrier they encounter, leading to a conduction mechanism
across the potential barrier known as thermionic emission. The distinction between thermionic emis-
sion and the tunneling process lies in the fact that in thermionic emission only the height (energy) of
the barrier is significant, not its shape [16]. The current density for thermionic emission can be cal-
culated similarly to tunneling current density. The key difference is that thermally induced electrons
have higher energy than the barrier, resulting in a transmission probability equal to 1. Therefore,
thermionic emission transport can be expressed as:
Z Ec,max
4πmeff q
Jthermo−ionic = Nsup (Ez )dEz (2.3)
h3 Ec,min

2.4 ReRAM characteristics

In this section, key attributes and properties essential for comprehending resistive switching in VCM
ReRAMs are briefly introduced and clarified.

2.4.1 Electro-forming of filament

Prior to any resistance switching in the Valence Change Memory device, an essential electroforming
step is required. The electroforming step involves applying a voltage to induce a local reduction of the
oxide, leading to the generation of oxygen vacancies. This process results in the creation of doubly
ionized oxygen vacancies that migrate toward the electrode [19]. As a consequence, a conductive
filament is formed inside the oxide layer, as illustrated in Figure 2.6.
Following the formation of the conductive filament, the Valence Change Memory device can be set
to either a High Resistance State or a Low Resistance State depending on the polarity of the applied
voltage and the characteristics of the oxide material. Typically, filament growth initiates from the
ohmic electrode, characterized by a lower potential barrier, and progresses towards the active electrode.
Titanium is commonly employed as the material for the ohmic electrode, while gold or platinum is
frequently used for the active electrode.
The forming voltage, which denotes the voltage which is required for the electroforming process, is
greater than the Set voltage and is depending on the thickness of the oxide (middle layer). During
the forming process, characterized by elevated voltages and currents, the localized temperatures are
expected to be very high [13].
Theoretical principles 9

Figure 2.6: By applying the forming voltage, a filament of oxygen vacancies is built starting
from the Ohomic electrode.

2.4.2 SET and RESET in VCM

The stable switching process in VCM devices transpires in the region between the plug and the
active electrode, with the plug representing the segment of the oxygen vacancy filament present in the
OFF state. The Set process involves applying a negative voltage to the active electrode, causing the
migration of oxygen vacancies toward the active electrode and a consequent reduction in the potential
barrier. As a result, the cell achieves higher conductivity, entering the ON (LRS) state. Reset, on
the other hand, is accomplished by applying a reversed-polarity voltage, causing oxygen vacancies to
move back into the filament (plug) and leading to re-oxidation. This action returns the cell to the
OFF (High Resistance State, HRS) state, characterized by low conductivity. A detailed depiction of
the Set and Reset processes for a bipolar switching VCM is presented in Figure 2.7 [13].
Theoretical principles 10

Figure 2.7: The schematic representation illustrates the Set and Reset characteristics of a Va-
lence Change Memory (VCM) device. In the OFF state (A), there is no conducting filament
between the electrodes. Upon applying a negative voltage, ions are attracted to the top electrode,
resulting in the creation of a filament during the Set process (B) and leading to the ON state (C)
with a conducting filament between both electrodes. Conversely, applying a positive voltage to
the active electrode during the Reset process (D) prompts the migration of oxygen ions back into
the plug, resulting in a measurable reset. The device then returns to the OFF state (A). This de-
piction is derived from [13] (https://marketplace.copyright.com/rs-ui-web/mp/license/4df5fdc4-
11c2-46d8-a901-d8de4401c120/1bc1ebd0-c59a-4b64-a722-a258b16fb766).

2.4.3 Complementary switching

Beside the bipolar switching, another common type of switching is called complementary switching
mechanism. The concept of the complementary switching is based on creating depletion areas near
top or bottom electrode. Hence, for each voltage polarity, a depletion layer is created between the
filament and one of the electrodes. In result, both ON and OFF state can be seen in positive and
negative applied voltages. An example of the complementary switching is displayed in Figure 2.8
In complementary switching, the initial state is OFF, in which a gap between the bottom electrode
and filament is present. By applying a positive voltage this gap is filled due to the migration of the
oxygen vacancies, hence, the cell will be in ON state (Set1). By increasing the voltage further, a
re-oxidation happens in near the top electrode, so the depletion gap is created between that and the
filament, this will turn the cell into OFF state (Reset1). By switching the polarity and applying a
negative voltage the gap near the top electrode can be closed and the cell will get to ON state again
(Set2). by increasing the value of the negative voltage, a gap is created near the bottom electrode
which put the cell once more in OFF state (Reset2) [20].
Theoretical principles 11

Figure 2.8: Sketch of complementary switching. In both voltage polarities HRS and LRS are
present. Applying a positive voltage results in Set 1 then Reset 1. Applying negative voltage
leads to Set 2 followed by Reset 2.

2.4.4 Bi-layer ReRAM

The bi-layer concept represents an innovative approach in Redox-based Random Access Memory
(ReRAM) devices, involving the usage of two oxide layers between the electrodes. One possible
configuration is the usage of homogeneous layers, where an oxygen-deficient layer is positioned at
the Ohmic electrode, and a stoichiometric layer is placed at the active electrode, as an example:
TiO2/TiO(x < 2). Another approach involves the usage of different oxide types with distinct band
gap energies and formation energies, resulting in a heterogeneous bi-layer structure. An example of
such a combination is Al2 O3 /TiO2 [13].
Bi-layer Resistive Random Access Memories (RRAMs) exhibit promising attributes in terms of resis-
tive switching behavior, demonstrating features essential for ReRAM applications. The incorporation
of two oxide layers between electrodes creates advantageous characteristics to bi-layer ReRAMs, in-
cluding robust endurance, accelerated switching speed, and prolonged data retention capabilities [21].
Integrating a bi-layer structure into ReRAMs not only enhances resistive switching behavior but also
plays an important role in improving the device’s retention and stability during prolonged use. The
utilization of two oxide layers between electrodes contributes to the overall robustness of the ReRAM,
ensuring sustained performance over extended periods [22].
The ReRAM investigated in this thesis features a heterogeneous bi-layer structure with a TaOx /HfOx
interface.

2.5 Materials
In this section the two oxides which are used in the stacked VCM cell under study are introduced and
briefly reviewed.

2.5.1 Hafnium oxide - HfOx

ReRAMs utilizing HfO2 have been extensively investigated. This is primarily because of its common
application in transistor gate stacks and its notable resistive switching properties. The discussions
Theoretical principles 12

encompass modifications to the Schottky barrier at electrodes, the migration of electrode ions, space-
charge-limited current, and the phenomenon of filamentation caused by oxygen deficiency [23].
Within the Valence Change Mechanism, materials like HfOx have demonstrated the generation of
oxygen vacancies induced by dielectric effects. in the presence of an oxygen vacancy filament, the LRS
is established when the oxygen-deficient path connects to the metal electrode. Conversely, the HRS
occurs when the filament is disrupted due to the re-migration of oxygen ions back to the vacancies.
Despite extensive studies on ReRAMs based on HfOx , the precise role of oxygen vacancies in electron
transport remains unclear [23].

2.5.2 Tantalum oxide - TaOx

Until now, resistance switches based on amorphous TaOx have been widely discussed, demonstrating
notable attributes such as high switching speed, long endurance, and low operation bias voltage when
compared to ReRAM devices utilizing alternative materials [24].
The switching mechanism in resistance switches relying on oxygen-deficient materials like TaOx is
commonly attributed to the drift-diffusion of oxygen ions or vacancies in the presence of an electric
field.
Recent experimental investigations into the filament components in TaOx -based resistance switches
have revealed the presence of significant continuous Ta-rich regions, denoted as TaO1−x , within the
LRS. These Ta-rich regions play a crucial role as the conduction path in ReRAM cells. However, the
origin of the switching mechanism in TaOx -based resistance switches remains a topic of controversy.
There is a strong need to elucidate the switching processes, particularly focusing on understanding
the detailed structure and electronic behavior of the oxygen-deficient filament [24].
 Chapter 3

Physical model for the simulation

The primary simulation model employed in this thesis is derived from the physical model outlined by
Astrid Marchewka [12]. To investigate the VCM ReRAM cell with stacked layers, a 1-Dimensional
model, as illustrated in Figure 3.1, is simulated using Finite Element Method (FEM) modeling within
the COMSOL Multiphysics software [25].

Figure 3.1: Schematic of 1D model used in simulation of bi-layer VCM cell

The simulation involves the numerical solution of four simultaneous differential equations:

• Poisson equation

• Electron continuity equation

• Ion continuity equation

• Heat conduction equation

These four equations and their boundary and initial conditions are described in this chapter.

13
Physical model for the simulation 14

3.1 Poisson equation

3.1.1 General model equation

The Poisson equation is solved within the both oxide regions, it is as:

∇ (−ϵ0 ϵr ∇ϕ) = ρq. (3.1)

Here, the dependent variable is the electrostatic potential ϕ. ϵ0 represents the free space permittivity,
and ϵr is the relative permittivity, which is material-dependent. q denotes the elementary charge,
and ρ signifies the charge density. Since donor dopants have dominant presence, the charge density
considers them and the electron concentration (n), leading to:

ρ = −n + ND . (3.2)

In this context, ND represents the donor dopants concentration, with n indicating the electron density.
For the studied VCM cell, only doubly ionized oxygen vacancies are considered as donors, resulting
in:

ND = 2VÖ . (3.3)

Here, VÖ represents the oxygen vacancy concentration, determined by solving the ion continuity
equation introduced in Section 3.3. The electron density (n) can be obtained using:
 
2 Efn − Ec
n = NC √ F0.5 . (3.4)
π kB T
where kB is the Boltzmann constant, T is the temperature (obtained using the heat conduction
equation in Section 3.4), and NC is the effective density of conduction band states, given by:
  32
2πmeff kB T
NC = 2 . (3.5)
h2
Here, meff is the electron effective density, and h represents Planck’s constant.
In Equation 3.4, F0.5 denotes the Fermi-Dirac integral of the order of 21 , Efn is the electrons quasi-
Fermi level (found by solving the electron continuity equation in the Section 3.2), and EC is the
conduction band edge energy, calculated as:

EC = −χ − qϕ. (3.6)
1
where χ denotes the electron affinity. The Fermi-Dirac integral of the order of 2 is defined as:
∞
t0.5
Z
1
F0.5 (x) = dt. (3.7)
Γ(0.5 + 1) 0 exp(t − x) + 1
Efn −EC E−EC
Here, Γ is the gamma function, and for this case, x = kB T and t = kB T .

3.1.2 Boundary conditions

To establish the boundary conditions for the Poisson equation, it is essential to determine the electric
potential (ϕ) at the interfaces of both the top and bottom electrodes with the oxide layers. In our
Physical model for the simulation 15

model, there exist two Schottky contacts: one at the interface of the top electrode and oxide 1, and
another at the interface between oxide 2 and the bottom electrode. The general equation describing
the Schottky contact potential is given by:

kB T ND
ϕSchottky = −[ϕBn0 + invF D0.5 ( )]. (3.8)
q NC
√
Here, invF D0.5 is defined as the inverse Fermi-Dirac integral multiplied by 2π , and ϕBn0 denotes
the total barrier height at the interfaces, obtainable from Equation 2.1. The term kBqT .invF D0.5 ( N
NC )
D

represents the disparity between the quasi-Fermi level and the conduction band energy. An illustrative
example of the Schottky contact band diagram is presented in Figure 3.2.

Figure 3.2: Band diagram at a Schottky contact, the Schottky potential of the interface is the
green dotted line.

Due to variations in donor concentration within the oxide region, a discrepancy may arise between
the intrinsic Fermi level and quasi-Fermi level. At each metal-semiconductor interface, this difference
can be expressed as:
 
Eg krmB T ND
ϕdiff = + invF D0.5 . (3.9)
q·2 q NC
Here, Eg denotes the band gap energy of the oxide, representing the energy gap between the valence
and conduction energy bands, expressed as:

Eg = EC − EV . (3.10)

Since the potentials are needed for both the top and bottom interfaces, one can determine the relative
change in potentials for the top and bottom electrodes and add this to the potential of the top
electrode:

ϕrelative = ϕdiff, top − ϕdiff, bottom . (3.11)

To find this difference, consider that the quasi-Fermi levels of two semiconductors (or oxides) tend
E
to align at their intersection. By utilizing their electron affinity, the difference in their q·2g can be
Physical model for the simulation 16

determined, equating to:

Eg,1 Eg,2 χ1 − χ2
− = . (3.12)
2q 2q q
in which χ1 and χ2 are the electron affinities
 ofthe oxide 1 and oxide 2.
Similarly, the difference in kBqT · invF D0.5 N
NC
D
can be straightforwardly calculated by:

kB Ttop ND,top kB Tbottom ND,bottom

.invF D0.5 ( )− .invF D0.5 ( ). (3.13)
q NC,top q NC,bottom
where ND , NC , and T should represent the values at the interfaces. Combining Equations 3.11, 3.12,
and 3.13, the relative potential difference at the top contact interface is given by:

χ1 − χ2 kB Ttop ND,top kB Tbottom ND,bottom

ϕreletive = + .invF D0.5 ( )− .invF D0.5 ( )). (3.14)
q q NC,top q NC,bottom
During the simulation, we apply a voltage to the upper contact to assess the properties and function-
ality of Valence Change Memory. It is essential to consider this applied voltage, denoted as Vapp , and
add it to the potential of the upper electrode.
In summary, the boundary conditions for the Poisson equation are established accordingly:

ϕtop = ϕSchottky−top + ϕrelative + Vapp . (3.15)

ϕbottom = ϕSchottky−bottom . (3.16)

However in the implementation, of this thesis the metal work functions of the two electrode metals
are directly used as the boundary condition for simplicity:

ϕtop = −WM,top + qVapp . (3.17)

ϕbottom = −WM,bottom . (3.18)

where WM,top and WM,bottom represent the metal work function of the material used as the top and
bottom electrodes. In Figure 3.3, a schematic illustrating the applied boundary condition for Poisson
equation is presented.
Physical model for the simulation 17

Figure 3.3: Position and value of the parameters used as boundary conditions in Poisson equa-
tion.

3.2 Electron continuity equation

3.2.1 General equation

The continuity equation for electrons is expressed in terms of the drift-diffusion equation:

1 1
− ∇jn = − ∇(−qµn n∇ϕ + qDn ∇n + qnDTn ∇T ) = −Rn . (3.19)
q q

Here, µn represents the electron mobility, Dn is the electron diffusion coefficient, Rn denotes the
electron recombination rate, and jn is the electron current density. The term DTn stands for the
thermal diffusion coefficient for electrons.
Assuming the temperature gradient is considerably smaller than the gradient of electric potential and
electron density, this term can be neglected, resulting in:

1 1
−Rn = − ∇jn ≈ − ∇(−q.µn .n∇ϕ + q.Dn ∇n). (3.20)
q q

The diffusion coefficient and electron mobility are related by the Einstein relation [16]:
q
µn = Dn . (3.21)
kB T

Utilizing the Boltzmann approximation [16], the electron density n can be expressed as:

Efn − EC
n = NC · exp . (3.22)
kB T

Applying equations 3.6, 3.21 and 3.22, the equation 3.20 can be rewritten as:

∇(nµn ∇Efn ) = qRn . (3.23)

This represents the simplified version of the electron drift-diffusion equation, where the dependent
variable is Efn quasi-Fermi level energy.
Physical model for the simulation 18

3.2.2 Flux term (tunneling)

The flux term in this simplified electron drift-diffusion equation accounts for the recombination rate
of tunneling in the z direction (where z signifies the direction along the length of the 1D model). The
electron recombination rate is given by:

4.meff .π dE
Rn = Nsup (Ez )T C(Ez ) . (3.24)
h3 dz
Here, dE
dz represents the energy derivative in the z direction, and Ez is the energy in the z direction.
The supply function Nsup is defined as the difference between Fermi-Dirac distributions on the two
sides of the interface:

Ez −Efn,1
1 + exp(− kB .T )
Nsup = kB .T ln( Ez −Efn,2
). (3.25)
1+ exp(− kB .T )
Here the Efn,1 and Efn,2 are Fermi levels at the two sides of the interface.
The transmission probability T C can be evaluated as:
Z z0
2 p
T C(Ez ) = exp(− 2.meff (Ez0 − Ez )dz). (3.26)
ℏ 0

where ℏ is the reduced Planck’s constant, and z0 is the point at which the electron’s energy is to
Ez . A sketch of the tunneling mechanism through the potential barrier at the metal-oxide interface
is shown in Figure 3.4.

Figure 3.4: Sketch of tunneling mechanism through the potential barrier at the metal-oxide
interface.

Having Equation 3.23, the current density can be expressed as the gradient of the quasi-Fermi level:

J = µn n∇Efn . (3.27)

A more detailed explanation is provided at Gehring’s work [26].

Physical model for the simulation 19

3.2.3 Direct tunneling

One can determine the current resulting from direct tunneling, which occurs between two oxide layers.
When two oxide or semiconductor layers come into contact, similar to Schottky contacts, their Fermi-
level energies tend to align, causing a discontinuity in the conduction band. An illustration of this
discontinuity is provided in Figure 3.5.

Figure 3.5: Schematic of conduction band energy at oxide-oxide interface.

This discontinuity results in an energy barrier and electrons may tunnel through and by that overcome
this barrier and generate a current. This current density can be calculated as outlined in Subsection
2.3.3:
Z EC,max(oxide1)
4πmeff q
Jdirect = T C(Ez )Nsup (Ez )dEz . (3.28)
h3 EC,min(oxide2)

Here, EC,max and EC,min represent the maximum and minimum values of the conduction band energy
around the barrier, i.e., the conduction band energy in oxide 1 and oxide 2. The supply function and
transmission probability are determined using the same equations as discussed in Subsection 3.2.2.

3.2.4 Boundary conditions

To establish the boundary conditions, one can initialize the quasi-Fermi level with specific values. For
computational simplicity, setting the initial value of the quasi-Fermi level to zero is a common choice.
In this approach, the quasi-Fermi level at the top electrode is set to −qVapp (the negative sign is used
since the electron Fermi level is under consideration), while at the bottom electrode, it is set to 0:

Efn,top = −qVapp . (3.29)

Efn,bottom = 0. (3.30)

In Figure 3.6, a schematic illustrating the applied boundary condition for electron drift diffusion
equation is presented.
Physical model for the simulation 20

Figure 3.6: Position and value of the parameters used as boundary conditions in electron conti-
nuity equation.

3.3 Ion continuity equation

Ions, similar to electrons, are governed by a drift-diffusion equation. As previously discussed, consid-
ering only the doubly ionized oxygen vacancies, the drift-diffusion equation for ions can be expressed
as:

∂VÖ 1
− ∇(2q.µVÖ .VÖ ∇ϕ + 2q.DVÖ ∇VÖ − 2q.VÖ .DT VÖ ∇T ) = −RVÖ ,2 . (3.31)
∂t 2q
In this equation, VÖ is the concentration of the doubly ionized oxygen vacancy and is the dependent
∂V
variable, while ∂tÖ is the oxygen vacancy flux. Also µVÖ represent the oxygen vacancy mobility and
DT VÖ is the thermal diffusion coefficient of the oxygen vacancy. RVÖ ,2 is the total recombination rate
of oxygen vacancy and it can be out equal to zero.
Again, with the neglect of the temperature gradient concerning other gradients and using the rela-
tionship between mobility and diffusion coefficient, Equation 3.31 will be transformed into:

DVÖ q ∂VÖ
∇(−DVÖ ∇VÖ − (− (−∇ϕ)VÖ )) + = 0. (3.32)
kB T ∂t
This is the simplified oxygen vacancy continuity equation. Here, −∇Φ represents the electric field,
and DVÖ is the oxygen vacancy diffusion coefficient, which can be determined as:

WA
DVÖ = D0 exp(− ). (3.33)
kB T
Here, D0 is the pre-factor for the maximum diffusion coefficient (at infinite temperature), and WA
is the migration barrier height for oxygen ionization, with both parameters dependent on the oxide
material. For a more detailed explanation of oxygen ionization, refer to [12].
Implementing boundary conditions for the oxygen vacancy continuity equation involves ensuring that
there is no ion exchange at the metal-oxide interface, resulting in the flux terms being set to zero at
interface towards the metal material.

DVÖ q
n̂top .(−DVÖ ∇VÖ − (− (−∇ϕ)VÖ )) = 0. (3.34)
kB T
Physical model for the simulation 21

DVÖ q
n̂bottom .(−DVÖ ∇VÖ − (− (−∇ϕ)VÖ )) = 0. (3.35)
kB T
Here n̂top and n̂bottom are normal vectors at top electrode and bottom interface pointing towards the
electrodes. Hence, they have the coordinates as:

n̂top = (0, 0, −1). (3.36)

n̂bottom = (0, 0, 1). (3.37)

DV q
Also the terms (−DVÖ ∇VÖ − (− kBÖT (−∇ϕ)VÖ ) represent the flux term of the oxygen vacancy con-
tinuity equation and can be represented by JˆVÖ . Therefore, the boundary conditions for the oxygen
vacancy continuity equation can be written as:

n̂top .JˆVÖ = 0. (3.38)

n̂bottom .JˆVÖ = 0. (3.39)

In Figure 3.7, a schematic illustrating the applied boundary condition for oxygen vacancy continuity
equation is presented.

Figure 3.7: Position and value of the parameters used as boundary conditions in ion drift-
diffusion equation.

3.4 Heat conduction equation

The fourth differential equation that can be numerically solved is the heat conduction equation, derived
from the principles of Joule heating:

∇(−κ∇T ) = j(−∇ϕ). (3.40)

Joule heating, also referred to as the heating resulting from the flow of electric current, is characterized
by the thermal conductivity κ, the total current density j, and the electric field represented by
−∇pP hi. In this equation, the temperature T serves as the dependent variable, and the focus of the
simulation is on the study of Joule heating in oxide regions.
Physical model for the simulation 22

The model employs the steady-state formulation of the heat conduction equation, as the equilibrium
states of electrons and the lattice are achieved more rapidly compared to the time scales of ionic drift
and diffusion.
Ambient temperature is employed as the boundary condition for the heat conduction equation at both
the top and bottom electrodes.
Ttop = Tenvironment (3.41)

Tbottom = Tenvironment (3.42)

In Figure 3.8, a schematic illustrating the applied boundary condition for Joule heating equation is
presented.

Figure 3.8: Position and value of the parameters used as boundary conditions in heat conduction
equation.
 Chapter 4

Simulations

In this chapter, simulations are conducted to investigate the behavior of VCM Resistive Random-
Access Memory (RRAM) devices. These simulations, based on a rigorous physical model, aimed to
uncover the complex mechanisms underlying resistive switching phenomena in VCM ReRAMs. By
systematically varying parameters such as applied voltages, pulse time, and others, the simulations
provided valuable insights into device performance. Results were analyzed to understand oxygen
vacancy dynamics, conduction band energy profiles, and current-voltage characteristics.

4.1 Phase I-a

In the initial phase of the simulation, pulsed voltages are applied to the VCM bi-layer contacts using
the parameters outlined in Table 4.1. The simulation utilizes the equations presented in Chapter 3,
and a COMSOL project is created based on the formulation of the four differential equations without
employing any built-in physics of the COMSOL software.
The applied voltages on the top electrode exhibit pulsed behavior, and the results are collected after
the pulses conclude. A schematic depicting of the applied voltage behavior used in this phase is
illustrated in Figure 4.1.
The pulsed signal is investigated for both positive and negative polarity, with amplitudes from −2 V
to +2 V.
The objective of this study is to comprehend the impact of the polarity and amplitude of the external
voltage on the distribution of oxygen vacancies along the oxide layers. Alterations in the distribution
of oxygen vacancies can potentially influence the conductivity of the VCM cell.
As reported in Appendix A, the temperature changes are insignificant, hence, in further studies results
regarding the temperature are not reported.

23
Simulations 24

Vapp

applied voltage Vapp [V]

0 1
time t [s]
Figure 4.1: Vapp (t) - Applied voltage function for simulation phase I-a. The increase and decrease
of the voltage happens in a short time (10−5 s).

Parameter Symbol Value Unit

oxide 1 thickness toxide1 2 nm
oxide 2 thickness toxide2 5 nm
electron affinity of oxide 1 χ1 2.4 eV
electron affinity of oxide 2 χ2 3.8 eV
band gap energy of oxide 1 Eg1 5.8 eV
band gap energy of oxide 2 Eg2 4.4 eV
relative permittivity of oxide 1 ϵr1 22
relative permittivity of oxide 2 ϵr2 20
thermal conductivity of oxide 1 κ1 1 W/(m·K)
thermal conductivity of oxide 2 κ2 1 W/(m·K)
electron mobility of oxide 1 µn1 5 cm2 /(Vs)
electron mobility of oxide 2 µn2 10 cm2 /(Vs)
effective electron mass for oxide 1 meff,n1 0.9 · me kg
effective electron mass for oxide 2 meff,n2 1.7 · me kg
initial density of donor dopants for oxide 1 ND1 1 · 1016 1/cm3
initial density of donor dopants for oxide 2 ND2 1 · 1020 1/cm3
diffusion coefficient pre-factor for oxide 1 D01 1 · 10−7 cm2 /s
−7
diffusion coefficient pre-factor for oxide 2 D02 1 · 10 cm2 /s
oxygen ion activation energy for oxide 1 WA1 0.5 eV
oxygen ion activation energy for oxide 2 WA2 0.5 eV
top electrode metal work function WM,top 5.2 eV
bottom electrode metal work function WM,bottom 4 eV
Table 4.1: Parameters used in simulation (me denotes the electron mass).
Simulations 25

4.1.1 Oxygen vacancy concentration

The initial outcome to analyze following the application of an external voltage for one second to the
VCM structure is the oxygen vacancy concentration (VÖ ). As detailed in Section 3.3, this variable is
the subject of the simplified continuity Equation 3.32.
Figure 4.2 illustrates the distribution of oxygen vacancies along the z axis in the two oxides after the
completion of the pulsed voltage application. Notably, even in the absence of an external voltage,
there is movement of the oxygen vacancies attributable to the diffusion components of the continuity
equation. The plot reveals that applying a positive voltage to the top electrode leads to a reduction in
the concentration of oxygen vacancies in oxide 1, with the distribution in oxide 2 remaining relatively
unchanged. Conversely, the application of a negative biased voltage increases the density of oxygen
vacancies in oxide 1 while leaving oxide 2 more or less unaffected. This observation indicates that the
polarity of the applied voltage dictates the direction of the shift in the oxygen vacancy distribution
profile.

10 26

10 25
cocentration VOB [1=m3 ]

10 24

10 23

+2 V pulse
+1 V pulse
10 22 0 V pulse
-1 V pulse
-2 V pulse
10 21
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.2: VÖ (z) - Oxygen vacancy concentration along the bi-layer ReRAM cell. The black
line represents the condition without applied voltage, the orange lines depict the scenario with
positive voltage application, and the blue lines illustrate the effect of negative voltage pulses.
Due to the simulation geometry, z < 2 nm corresponds to oxide 1 (HfOx ), and z > 2 nm
represents oxide 2 (TaOx ) areas.

Subsequently, the impact of the amplitude of the externally applied voltage pulse on the oxygen
vacancy concentration distribution can be examined. Figure 4.3 presents a comparison of the changes
in oxygen vacancy concentration in oxide 1 following the application of various positive voltage pulses.
Noticeably, an increase in the amplitude of the positive bias stimuli results in a progressively smaller
oxygen vacancy concentration compared to the equilibrium case.
Similarly, the anticipated effect is expected to be present in cases involving negative voltage pulses.
Figure 4.4 provides a comparison of different amplitudes of applied negative voltage pulses to the top
electrode against the equilibrium case. As evident from the figure, an increase in the amplitude of the
negative voltage pulse leads to a higher concentration of oxygen vacancy within oxide 1 compared to
Simulations 26

10 26

10 25
cocentration VOB [1=m3 ]

10 24

10 23

+2 V pulse
+1.5 V pulse
10 22 +1 V pulse
+0.5 V pulse
0 V pulse
10 21
0 0.5 1 1.5 2 2.5 3
position z [nm]
Figure 4.3: VÖ (z) - Oxygen vacancy concentration along the bi-layer ReRAM cell. The diagram
is zoomed in to observe details near and inside the oxide 1 regions. The black line represents
the equilibrium case (without external voltages), and the orange lines depict the oxygen vacancy
distribution after applying different amounts of positive voltage pulses.

the equilibrium state.

Hence, the influence of application of pulsed voltages on the oxygen vacancy concentration is observed
in both positive and negative bias scenarios. These alterations in ionized oxygen vacancies are expected
to influence the conductivity of the ReRAM structure.

4.1.2 Conduction band energy

The second parameter of interest in the pulsed input simulation is the conduction band edge energy
(EC ). As discussed in Section 3.1, the conduction band energy is determined after solving the Poisson
equation and finding the potential. It also depends on the electron affinity of the oxide, as indicated
by Equation 3.6. The conduction band energy profile plays a crucial role in the conductivity of the
VCM cell, as it influences the number of electrons (charge carriers) according to Equation 3.4.
Figure 4.5 illustrates the conduction band energy along the z axis after the application of different
positive voltage pulses. This graph reveals the relationship between the amplitude of the external
positive voltage applied to the top electrode and the conduction band energy in the two oxides. As
observed in the plot, only the application of +0.5 V leads to a meaningful change in the conduction
band energy, resulting in a decrease in the energy of the conduction band. No other amplitude of
positive voltage had a significant effect on the conduction band energy in the two oxides. This change
of behavior for a specific amplitude of voltage is not reasonable and the origins are unknown.
Figure 4.6 displays EC (z) after the application of negative voltage pulses in comparison to the situation
where no voltage is applied to the top electrode. As evident in the plot, applying different negative
voltage pulses to the VCM cell does not affect the energy of the conduction band. Consequently, the
distribution of EC remains largely unchanged after the pulse application.
Simulations 27

10 26

10 25
cocentration VOB [1=m3 ]

10 24

10 23

-2 V pulse
22
-1.5 V pulse
10 -1 V pulse
-0.5 V pulse
0 V pulse
10 21
0 0.5 1 1.5 2 2.5 3
position z [nm]
Figure 4.4: VÖ (z) - Oxygen vacancy concentration along the bi-layer ReRAM cell. The diagram
compares different amplitudes of applied negative voltage pulses (blue lines) with the equilibrium
case (black line). The diagram concentrates on areas of oxide 1.

3
+2 V pulse
+1.5 V pulse
conduction band energy EC [eV ]

2.5 +1 V pulse
+0.5V pulse
0 V pulse
2

1.5

0.5

0
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.5: EC (z) - Conduction band energy along the bi-layer VCM cell. The black line
represents the energy without the pulsed voltage, orange lines represent the situation after
applying positive voltage pulses with no profound change in EC , and the green line represents
the energy for the voltage pulse that has an effect on the conduction band energy level.
Simulations 28

3
-2 V pulse
-1.5 V pulse
conduction band energy EC [eV ]

2.5 -1 V pulse
-0.5 V pulse
0 V pulse
2

1.5

0.5

0
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.6: EC (z) - Conduction band energy along the bi-layer VCM cell. The black line
represents the energy without the pulsed voltage, blue lines represent the situation after applying
negative voltage pulses.
Simulations 29

4.2 Phase I-b

In this phase of the simulations, the oxygen vacancy continuity equation is excluded, and the extracted
oxygen vacancy distribution (VÖ (z)) from Section 4.1 is utilized. In this step, a range of positive and
negative voltages are applied to the previously pulsed VCM layers. Other three main differential
equations are still used in the COMSOL project in the same way implemented in Section 4.1 also used
parameters are again the ones mentioned in Table 4.1.
Figure 4.7 illustrates the applied voltage on the top electrode in phase I-b as the simulation progresses
(Vapp (t)).
The primary objective of this sweep voltage application after the pulses is to investigate the conduc-
tivity of the stack-layer VCM under different oxygen vacancy distributions. In other words, the effect
of the ionized oxygen vacancy distribution on the possible resistive switching of the device is studied.

3
applied vooltage Vapp [V ]

-3

0 3 6 9 12
time t [s]
Figure 4.7: Vapp (t) - Applied voltage on the top electrode in the time evolution for phase I-b of
the simulation.

4.2.1 Current-voltage characteristics

The most crucial outcome to study in this phase is the current density at different applied voltages,
i.e., the JV characteristics. The current density can be determined using Equation 3.27, subsequent
to finding the electron density and solving the quasi-Fermi level differential equation. This current
density (J) can be obtained at each geometrical point in the device. Typically, the maximum current
density is observed near the middle of oxide 2. Figure 4.8 shows examples of current density profile in
different voltages of the sweep, in all of them middle point of the oxide 2 can be seen as the maximum
point. Hence, this point is usually utilized as the reference. In the following results middle of oxide 2
(z = 4.5 nm) is chosen as the reference point to measure the current density.
Figure 4.9 illustrates the JV characteristics of the VCM cell subjected to pulses with positive, negative,
and zero volts. Firstly, it is evident that the amplitude of the current density is much lower in negative
biases compared to positive biases, possibly due to the asymmetric configuration of the stacked layer
structure. Secondly, by comparing the current density under different pulsed conditions, it can be
observed that although small, applying a negative voltage pulse increases the current density in both
positive and negative applied biases. In contrast, applying a positive voltage pulse decreases the
current density compared to the case when no pulse is applied to the bi-layer VCM. This difference
Simulations 30

10 10

10 0
current density J [A/m2]

10 -10

Vapp = -3V
Vapp = -1V
10 -20
Vapp = +1V
Vapp = +3V
middle of oxide 2
-30
10
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.8: J(z) - Current density distribution at different Vapp , where z = 4.5 nm (dashed line)
represents the middle of oxide 2.

in JV behavior can be attributed to the alteration in oxygen vacancy distribution resulting from the
pulsed signals applied to the top electrode in Section 4.1.
To gain further insights into the influence of pulse voltage on the conductivity of the structure, JV
characteristics of the positively pulsed cases can be examined, as illustrated in Figure 4.10. As evident
from the plots, further increasing the voltage of the positively pulsed bias leads to a decrease in current
density when the VCM is subjected to sweep voltage testing. In other words, applying higher voltages
as positive bias pulses results in lower conductance observed in the cell. However, it is important to
note that these changes are relatively small.
The behavior observed for negative voltage pulses mirrors that of positive pulses. As depicted in
Figure 4.11, increasing the amplitude of negative voltage pulses leads to higher current densities. In
essence, applying greater negative voltages during pulse steps enhances conductivity, aligning with
the changes in oxygen vacancy distribution following exposure to pulsed voltage in VCM devices.
Indeed, the observed changes in current density suggest a direct correlation between the applied
pulsed voltages and the distribution of oxygen vacancies, which serve as charge carriers within the
VCM structure. This behavior underscores the significance of oxygen vacancies in modulating the
conductivity of the device, as their migration across the oxide layers influences the flow of current
under different bias conditions.
Simulations 31

a b
10 2 10 10

10 0

10 5
-2
10
current density J [A/m2]

current density J [A/m2]

10 -4
10 0

10 -6

10 -5
10 -8

10 -10
10 -10
+1 V pulsed +1 V pulsed
10 -12
without pulse without pulse
-1 V pulsed -1 V pulsed
10 -14 10 -15
-3 -2 -1 0 0 1 2 3
applied voltage V app [V] applied voltage V app [V]
Figure 4.9: JV characteristics of the VCM cell under a) negative voltages and b) positive voltages.
The black lines represent the case where no voltage pulse is applied, the blue line shows the
current for the case where a −1 V pulse is applied, and the +1 V pulse case is depicted with the
orange line (all the currents are measured in z = 4.5 nm).
Simulations 32

a b

10 0

10 5

10 -2
current density J [A/m2]

current density J [A/m2]

10 -4 10 0

10 -6

10 -5

10 -8
+2 V pulsed +2 V pulsed
+1 V pulsed +1 V pulsed
without pulse without pulse
10 -10 10 -10
-3 -2 -1 0 0 1 2 3
applied voltage V app [V] applied voltage V app [V]
Figure 4.10: JV characteristics of the VCM cell under a) negative voltages and b) positive
voltages. The black lines represent the case where no voltage pulse is applied, the orange lines
show the current for the cases where a positive pulse is applied (all the currents are measured
in z = 4.5 nm).
Simulations 33

a b
10 2 10 10

10 0

10 5
-2
10
current density J [A/m2]

current density J [A/m2]

10 -4
10 0

10 -6

10 -5
10 -8

10 -10
10 -10
-2 V pulsed -2 V pulsed
10 -12
-1 V pulsed -1 V pulsed
without pulse without pulse
10 -14 10 -15
-3 -2 -1 0 0 1 2 3
applied voltage V app [V] applied voltage V app [V]
Figure 4.11: JV characteristics of the VCM cell under a) negative voltages and b) positive
voltages. The black lines represent the case where no voltage pulse is applied, the blue lines
show the current for the cases where a positive pulse is applied (all the currents are measured
in z = 4.5 nm).
Simulations 34

4.3 Phase I-c

In this phase of the study, the focus shifts to understanding the impact of the duration for which
voltage is applied to the VCM. Using the same methodology and parameters as in Phase I-a, the
duration of the applied pulses is extended to 5 s and 10 s, allowing for a comparison with the results
obtained from one second pulses in Section 4.1. Figure 4.12 shows the applied voltage as a function
of the simulation time for this phase. By varying the pulse duration, insights can be gained into how
the temporal aspect of the voltage application influences the behavior of the VCM device. The used
voltages in this phase are: -2 V, -1 V, 0 V, 1 V and 2 V.

Vapp
applied voltage Vapp [V]

0 tsim
time t [s]
Figure 4.12: Vapp (t) - Applied voltage function for simulation phase I-c. The increase and
decrease of the voltage happens in a short time (10−5 s).

4.3.1 Oxygen vacancy concentration

The initial property to examine following the alteration in pulse duration is the distribution of oxygen
vacancy concentration (VÖ (z)) across the two oxides. As detailed in Subsection 4.1.1, this distribution
was previously analyzed for pulse duration of one second. Now, the impact of extending the pulse
duration can be investigated to determine any variations in the oxygen vacancy distribution.
To establish a baseline, the impact of time in the absence of a pulsed voltage is examined first.
Figure 4.13 illustrates the variation in the oxygen vacancy distribution across the ReRAM structure
under the absence of external voltage application, at different simulation times. As depicted in the
plot, even without applying any voltage to the top electrode, there is noticeable movement of oxygen
vacancies. This movement is primarily driven by the diffusion component of the ion continuity equation
(Equation 3.32). With longer simulation times, this phenomenon becomes more pronounced, leading
to an increase in the concentration of oxygen vacancies within the oxide 1 region (the region with
initially lower concentration).
The subsequent analysis involves examining the influence of pulse voltage duration under positive
bias conditions. Figures 4.14 and 4.15 present the spatial distribution of oxygen vacancies at various
pulse durations for applied voltages of +1 V and +2 V respectively. The results for the case when the
applied pulsed voltage is +1 V reveal that while the positive bias exhibits a reverse effect compared
to the duration of the pulse, it does not entirely prevent the increase in oxygen vacancy concentration
Simulations 35

10 27

10 26
cocentration VOB [1=m3 ]

10 25

10 24

10 23

10 s pulse
10 22
5s pulse
1 s pulse
10 21
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.13: Variation of oxygen vacancy concentration (VÖ (z)) across different positions of the
VCM structure under an external voltage pulse of 0 V, applied for durations of 1 s (black line),
5 s (green line), and 10 s (purple line).

in the oxide 1 region. However, the increment is lower compared to the scenario without applied
bias. Conversely, in the plots corresponding to Vapp = +2 V, it is evident that a higher magnitude
of positive pulsed voltage can counteract the effect of diffusion. With increased pulse duration, the
concentration of oxygen vacancies decreases in the oxide 1 regions. This suggests that the magnitude
of the positive bias significantly influences the changes in oxygen vacancy distribution along the VCM
structure.
In the final comparison, we explore the impact of simulation duration on cases where a negative
pulse voltage is applied to the top electrode of the VCM. Figures 4.16 and 4.17 illustrate the oxygen
vacancy concentration across the structure in different pulse time durations for −1 V and −2 V pulses,
respectively. The results indicate that increasing the duration of the applied negative pulsed voltage
leads to a rise in the oxygen vacancy increment in the oxide 1 region. Notably, both negative bias
application and longer simulation times produce a similar effect on this property, accelerating the
increase in oxygen vacancy concentration compared to scenarios with no applied pulse. Furthermore,
with an increase in the magnitude of the negative bias, the reduction in oxygen vacancy concentration
in oxide 2 regions becomes more pronounced as the pulse duration increases.

4.3.2 Conduction band energy

The second parameter that warrants examination after applying pulses of varying durations is the
conduction band energy (EC (z)). This property, obtained using the methods described in Subsection
4.1.2, is crucial for understanding electrical conduction of the device.
To begin, we analyze the impact of time when no external voltage is applied to the top electrode.
Figure 4.18 illustrates the conduction band energy across different regions of the structure. As evident
in the plot, increasing the simulation time without any external bias applied to the top electrode leads
to a rise in the energy of the conduction band across all regions of the VCM. This change can be
Simulations 36

10 27

10 26
cocentration VOB [1=m3 ]

10 25

10 24

10 23

10 s pulse
10 22
5s pulse
1 s pulse
10 21
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.14: Variation of oxygen vacancy concentration (VÖ (z)) across different positions of the
VCM structure under an external voltage pulse of +1 V, applied for durations of 1 s (black line),
5 s (green line), and 10 s (purple line).

attributed to the diffusion contribution of the continuity equations utilized in the model.
One can then examine how the conduction band energy changes when an external voltage is applied
in the form of a pulse. While the results for a pulse duration of 1 s have been previously discussed
in Subsection 4.1.2, it is now valuable to observe the impact of prolonging the pulse duration on
the conduction band across the oxide layers. Figure 4.19 illustrates the conduction band at various
positions when voltage pulses of different polarities and magnitudes are applied for a duration of 5 s.
from the plot, it can be inferred that for an applied voltage magnitude of 1 V, regardless of polarity,
the conduction band energy remains relatively stable even as the pulse duration increases to 5 s.
However, with a higher voltage amplitude, noticeable changes in the conduction band energy become
apparent at this pulse duration. Specifically, for cases with a magnitude of 2 V, there is an increase in
the conduction band energy when the pulse polarity is negative, whereas a decrease in EC is observed
when a positively biased voltage is applied.
One can further explore the alterations in the conduction band energy with a longer simulation
time of ten seconds. Extending the simulation duration allows for a deeper understanding of how the
conduction band evolves across the oxide layers under prolonged voltage pulses. Figure 4.20 illustrates
the conduction band at various positions when voltage pulses of different polarities and magnitudes
are applied for a duration of 10 s. According to the results depicted in the plot, when the duration of
the pulse application is extended to 10 s, its effect on the conduction band energy becomes notably
evident. In cases of positive biases, the conduction band energy experiences a decrease, particularly
pronounced in oxide 2 regions, and this effect becomes more evident with increasing voltage amplitude.
Conversely, for negative biases, there is an increase in the energy level of the conduction band, with
this phenomenon being further enhanced when higher magnitudes of negative voltage are used as the
pulse.
Simulations 37

10 28

10 26
cocentration VOB [1=m3 ]

10 24

10 22

10 20 10 s pulse
5s pulse
1 s pulse
10 18
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.15: Variation of oxygen vacancy concentration (VÖ (z)) across different positions of the
VCM structure under an external voltage pulse of +2 V, applied for durations of 1 s (black line),
5 s (green line), and 10 s (purple line).

10 28
cocentration VOB [1=m3 ]

10 26

10 24

10 s pulse
5s pulse
1 s pulse
10 22
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.16: Variation of oxygen vacancy concentration (VÖ (z)) across different positions of the
VCM structure under an external voltage pulse of −1 V, applied for durations of 1 s (black line),
5 s (green line), and 10 s (purple line).
Simulations 38

10 28

10 26
cocentration VOB [1=m3 ]

10 24

10 22

10 20 10s pulse
5s pulse
1s pulse
10 18
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.17: Variation of oxygen vacancy concentration (VÖ (z)) across different positions of the
VCM structure under an external voltage pulse of −2 V, applied for durations of 1 s (black line),
5 s (green line), and 10 s (purple line).

3
t sim = 10 s
t sim = 5 s
2.5
t sim = 1 s

2
energy EC [eV]

1.5

0.5

0
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.18: EC (z) - Conduction band energy in various positions of the VCM bi-layer struc-
ture at different simulation durations without external bias. The black line represents the 1 s
simulation time, the green line corresponds to 5 s duration, and the purple line depicts the 10 s
case.
Simulations 39

3
Vapp = + 2 V
Vapp = +1 V
2.5
Vapp = 0 V
Vapp = -1 V

2 Vapp = -2 V
energy EC [eV]

1.5

0.5

0
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.19: EC (z) - Conduction band energy in various positions of the VCM bi-layer structure
at 5 s duration of pulse applied voltage. The black line represents the case without external
voltage, the orange lines correspond to positive voltages, and the blue lines depict the case when
the applied pulses have negative polarity.

3
Vapp = +2 V
Vapp = +1 V
2.5
Vapp = 0 V
Vapp = -1 V

2 Vapp = -2 V
energy EC [eV]

1.5

0.5

0
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.20: EC (z) - Conduction band energy in various positions of the VCM bi-layer structure
at 10 s duration of pulse applied voltage. The black line represents the case without external
voltage, the orange lines correspond to positive voltages, and the blue lines depict the case when
the applied pulses have negative polarity.
Simulations 40

4.4 Phase I-d

In Section 4.3, it was evident that increasing the simulation time without applying any external
electrical bias leads to a more equal distribution of oxygen vacancy concentration between oxide 1 and
oxide 2 regions. Therefore, it is anticipated that further increases in simulation time will not alter the
oxygen vacancy distribution, reaching a relaxed state.

4.4.1 Reaching relaxed state

In this phase of the simulation, the parameters used are detailed in Table 4.1, and the implementation
mirrors that of Section 4.3. However, longer simulation times are examined (up to 60 seconds).
Initially, the device undergoes simulation without any voltage pulses for varying durations. Figure
4.21 illustrates the oxygen vacancy distribution after simulating in different times up to one minute
without applied voltage. As depicted in the plot, the oxygen vacancy distribution profile continues to
evolve with increasing simulation time, even up to one minute. Consequently, a steady state is not
attained within the 60-second timeframe. Thus, extending the simulation duration is necessary to
reach a relaxed-state profile, in which the mean concentration would not be changed over time.
In Figure 4.22, the changes in oxygen vacancy distribution in the absence of an external voltage are
depicted for different durations up to five minutes. As depicted in the plot, after three minutes of
simulation time, the oxygen vacancy concentration distribution appears to stabilize, showing little
change thereafter. Consequently, this distribution can be considered as the steady-state configuration
for further analysis.

4.4.2 Oxygen vacancy in different states

After establishing the equilibrium distribution of oxygen vacancies (relaxed state), it serves as a
reference for subsequent steps. As depicted in Figure 4.15, the application of a positive voltage as
large as 2 V leads to a decrease in the concentration of oxygen vacancies in the oxide 1 region. To
further explore this phenomenon, the relaxed state can be utilized as the initial condition for the
oxygen vacancy distribution, followed by the application of a positive voltage to the top electrode.
The changes in distribution relative to the relaxed state are then monitored when an external force is
applied to the VCM cell.
Figure 4.23 illustrates the mean value of the oxygen vacancy distribution over time when a positive
bias of magnitude 2 V is applied, starting from the relaxed state. The plot indicates that applying
a positive bias of 2 V to the top electrode results in a decrease in the oxygen vacancy concentration
within the oxide 1 region compared to the steady-state concentration. However, it takes some time
for the external bias to overcome diffusion effects before the decrement becomes apparent. Also 10
s is not enough to reach the initial concentration for the simulations, before trying to reach relaxed
state (from Equation 3.3).
To further explore the alterations of the oxygen vacancy concentration in the oxide 1 region from the
steady state under the influence of positive voltage, additional simulation times and various positive
voltage magnitudes are investigated.
Figure 4.24 illustrates the fluctuations in the mean value of VÖ in the oxide 1 region for different posi-
tive voltage levels over a simulation period of 100 s, all initialized from the steady-state concentration.
As depicted in the plot, applying a +1 V bias does not alter the concentration of oxygen vacancies
in the first oxide from the steady-state configuration, indicating an inability to overcome the redis-
Simulations 41

10 28

10 26
cocentration VOB [1=m3 ]

10 24
t sim = 10 s
t sim = 20 s
t sim = 30 s
t sim = 40 s
10 22
t sim = 50 s
t sim = 60 s
t sim = 0 s

10 20
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.21: VÖ (z) - Distribution of oxygen vacancies across the VCM layers after no voltage is
applied to the top electrode. For 0 s, 10 s, 20 s, 30 s, 40 s, 50 s and 60 s of simulation.

tribution caused by diffusion. However, higher positive voltages result in a decreased concentration
in the oxide 1 regions. With greater bias magnitude, the decrement occurs more rapidly, and the
distribution approaches the initial value (derived from doping value) sooner.
To unravel the behavior behind the initial gradual decrease followed by a subsequent rapid linear
decrease in the oxygen vacancy distribution within oxide 1 when positive voltages are applied, an in-
depth analysis of the simulation results at various points in the geometry and different time intervals is
conducted. Figure 4.25 illustrates the oxygen vacancy distribution initiated from the steady state and
subjected to a positive voltage. The figure indicates that over time, the peak in the oxygen vacancy
concentration gradually shifts towards the bottom electrode. Specifically, around tsim = 6 s, the peak
reaches the interface between oxide 1 and oxide 2. This observation aligns with the findings from
Figure 4.24, wherein the rapid reduction of the mean value of VÖ in oxide 1 begins approximately 6 s
after the voltage is applied.

4.4.3 Resistance states

In the preceding subsection, various states arising from the application of positive bias are explored,
each potentially exhibiting distinct electrical properties. Given the relevance of characterizing resis-
tance states in ReRAM studies, the focus now shifts to examining the resistance per area in these
different states.
Figure 4.26 illustrates the resistivity of the device in various states of the oxygen vacancy concentration,
starting from the relaxed state and the subsequent application of a +2 V bias. This conductivity is
depicted over time to capture the evolution of resistance under the influence of the applied voltage.
As depicted in Figure 4.26, transitioning away from the relaxed state and driving the peak in the
oxygen vacancy concentration distribution towards the bottom electrode initially leads to a significant
Simulations 42

10 30

10 25
cocentration VOB [1=m3 ]

10 20

10 15
t sim = 1 min
t sim = 2 min
10 t sim = 3 min
10
t sim = 4 min
t sim = 5 min

10 5
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.22: VÖ (z) - Distribution of oxygen vacancies across the VCM layers after no voltage
is applied to the top electrode. For 1 min, 2 min, 3 min, 4 min, and 5 min of simulation. The
result is stabilized after 3 minutes.

decrease in resistance. However, over time, this change becomes more gradual. This observation
suggests that after a certain duration, most of the oxygen vacancies have migrated into the oxide 2
area, diminishing the time-dependent effect on resistance.
The findings obtained here are consistent with recent research conducted by Aussen et al. [27]. In
their study, aimed at investigating the potentiation (enhancement of conductivity) of an Al2 O3 /TiOx
device, a sequence of short-term pulses (1 ms) was applied. Subsequently, the conductivity of the
device was measured at each step. They observed a rise in conductivity with increasing duration of
positive bias application. Similar to the results observed in this subsection, the change in conductivity
exhibited a non-linear behavior, initially increasing rapidly before tending towards saturation. The
same behavior was also observed during depression (lowering the conductivity with negative voltage).
To validate the observed phenomena in reverse, the oxygen vacancy distribution resulting from ap-
plying a +2 V bias for 30 s is preserved as the starting point (representing the LRS), initiating a new
simulation. This simulation commences from the LRS distribution and applies a negative bias.
Figure 4.27 shows the resistivity of the VCM device and the mean oxygen vacancy concentration in
oxide 1 at various time intervals, originating from the LRS distribution and subjected to a constant
−2 V bias. It is evident from the figure that upon the application of a negative bias, the oxygen
vacancies tend to migrate towards the oxide 1 regions, consequently leading to an increase in the
device’s resistance. Furthermore, the change in resistance is more pronounced at the outset, gradually
attenuating as the oxygen vacancy concentration peak reaches the interface between the two oxide
layers.
Simulations 43

10 26

10 25
[1=m3 ]

10 24
mean cocentration VO;ave
B

Voltage applied +2 V
Steady-State
Initial V Ö (ND)

10 23

10 22

10 21
0 2 4 6 8 10
time t [s]
Figure 4.23: VÖ,ave (t) - Mean value of the oxygen vacancy concentration in the oxide 1 region
at different time intervals when a positive bias of +2 V is applied to the top electrode, with the
initial distribution set to the relaxed state.
Simulations 44

10 26

10 25
[1=m3 ]

10 24
mean cocentration VO;ave

Vapp = +3 V
B

Vapp = +2 V
Vapp = +1 V
Steady-State
Initial V Ö (ND)
10 23

10 22

10 21
1 2 5 10 20 50 100
time t [s]
Figure 4.24: VÖ,ave (t) - Mean value of the oxygen vacancy concentration in the oxide 1 region
at different times, when voltages of +3 V (orange), +2 V (blue), and +1 V (green) are applied
to the top electrode. The initial distribution corresponds to the relaxed state.
Simulations 45

10 30

10 25
cocentration VOB [1=m3 ]

10 20

10 15

Steady-state
10
t=2s
10 t=6s
t=8s
t = 10 s
10 5
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.25: VÖ (z) - Oxygen vacancy concentration in different simulation domains when +2
V is applied to the top electrode. The red lines represent cases where the drop in the mean
value of oxide 1 is minimal, the black line indicates the point where the maximum oxygen
vacancy concentration is located near the interface, and blue lines correspond to cases where the
decrement of the oxygen vacancy mean value is rapid for oxide 1 regions.

#10 9
9 10 26

7 [1=m3 ]
10 25
resistance R [+ . m 2]

6
cocentration VO;ave

5
B

10 24
4

2 10 23

0 10 22
0 5 10 15 20 25 30
time t [s]
Figure 4.26: R(t) vs. VÖ,ave (t)- Resistance times area measured for a +0.5 V bias in the middle of
oxide 2 (z = 4.5 nm) versus the mean value of the doubly ionized oxygen vacancy concentration
in oxide-1, in time evolution during which a +2 V bias is applied to the top electrode starting
from the relax-state distribution.
Simulations 46

#10 9
10 10 26

9
10 25

[1=m3 ]
8
resistance R [+ . m 2]

7
10 24

cocentration VO;ave
B
6

5 10 23

4
10 22
3

2 10 21
0 5 10 15 20 25 30
time t [s]
Figure 4.27: R(t) vs. VÖ,ave (t)- Resistance times area measured for a +0.5 V bias in the middle of
oxide 2 (z = 4.5 nm) versus the mean value of the doubly ionized oxygen vacancy concentration
in oxide-1, in time evolution during which a −2 V bias is applied to the top electrode starting
from the LRS distribution.
Simulations 47

4.5 Phase II-a

In this phase of the simulation, the focus shifts to utilizing the COMSOL built-in Semiconductor
Module. Here, the oxygen vacancy continuity equation, as described in Section 4.1, is implemented.
Given the small fluctuations in temperature observed in previous phases (Appendix A), the heat con-
duction equation is disregarded. Instead of implementing the Poisson equation and electron continuity
equation manually, the Semiconductor Module is employed, offering various methods for calculating
and determining the required parameters. Detailed descriptions of the equations and methodologies
can be found in [28]. The parameters used remain consistent with those outlined in Table 4.1. Due
to issues arising from the initial values in the Semiconductor Module, leading to non-convergence
during simulation, an additional step is introduced. This step involves conducting a semiconductor
equilibrium study before initiating the time dependent study.
To compare the utilization of the built-in module with the self-implemented simulation, a pulsed
voltage is once again applied to the top electrode of the VCM, akin to Sections 4.1 and 4.3. Figure
4.28 illustrates the behavior of the applied voltages in this phase, with pulse durations (simulation
time) set at one and five seconds, and applied biases ranging from −2 V to +2 V.

Vapp
applied voltage Vapp [V]

0 tsim
time t [s]
Figure 4.28: Vapp (t) - Applied voltage function for simulation phase II-a. The increase and
decrease of the voltage happens in a short time (10−5 s).

4.5.1 Effect of applied voltage

In the initial analysis, a comparison of results is conducted across different pulse durations. Oxygen
vacancy concentration and conduction band energy following the application of voltage pulses to the
top electrode are examined. The former is determined by solving the oxygen vacancy continuity
equation, while the latter is obtained through the Semiconductor Module of the COMSOL software.
Figure 4.29 illustrates the distribution of oxygen vacancies, while Figure 4.30 displays the conduction
band energy level after pulses with duration of 1 s. Interestingly, unlike the self-implemented projects
which were used in previous cases, the oxygen vacancy distribution (VÖ (z)) and conduction band
energy (EC (z)) remain consistent across different magnitudes and polarities of the applied pulse voltage
when using the Semiconductor module of COMSOL Multiphysics software.
Simulations 48

10 26

10 25
cocentration VOB [1=m3 ]

10 24

10 23

+2 V pulse
22 +1 V pulse
10
0 V pulse
-1 V pulse
- 2V pulse
10 21
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.29: VÖ (z) - Oxygen vacancy distribution across different positions of the VCM structure
under external voltage pulses applied for durations of one second. The black line represents the
case without external voltage, while orange and blue lines depict positive and negative pulse
voltages, respectively. Results are derived from solving the ion continuity equation and using the
built-in Semiconductor module in COMSOL. Since all the lines have the same values everywhere,
the orange and black lines are not visible.

After examining the results from one-second pulses, a similar approach is employed for a pulse duration
of 5 s to corroborate the findings. Figure 4.31 demonstrates the distribution of oxygen vacancies, while
Figure 4.32 shows the conduction band energy level following five-second duration pulses. Similar
to the previous case, the oxygen vacancy distribution (VÖ (z)) and conduction band energy (EC (z))
exhibit consistency across various magnitudes and polarities of the applied pulse voltage when utilizing
the Semiconductor Module of COMSOL Multiphysics software, contrasting with the results obtained
from the self-implemented project.
Regarding the obtained results in this phase, the difference observed in the behavior when utilizing
the built-in Semiconductor Module compared to the implementation of equations using simple math-
ematical modules in COMSOL Multiphysics is currently unclear. These unexpected differences raises
questions regarding the underlying algorithms and numerical methods employed in the Semiconductor
Module, as well as the accuracy of its simulations compared to custom implementations.
Simulations 49

3
Vapp = +2 V
2.5 Vapp = +1 V
Vapp = 0 V
2 Vapp = -1 V
energy EC [eV]

Vapp = -2 V
1.5

0.5

-0.5
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.30: EC (z) - Conduction band energy across various positions of the VCM bi-layer
structure after applying pulsed voltage for one second. The black line denotes the case without
external voltage, while orange and blue lines represent positive and negative pulse voltages,
respectively. Results are obtained from solving the ion continuity equation and employing the
built-in Semiconductor Module in COMSOL. Since all the lines have the same values everywhere,
the orange and black lines are not visible.
Simulations 50

10 26

10 25
cocentration VOB [1=m3 ]

10 24

10 23

+2 V pulse
22 +1 V pulse
10
0 V pulse
-1 V pulse
-2 V pulse
10 21
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.31: VÖ (z) - Oxygen vacancy distribution across different positions of the VCM structure
under external voltage pulses applied for durations of five seconds. The black line represents
the case without external voltage, while orange and blue lines depict positive and negative pulse
voltages, respectively. Results are derived from solving the ion continuity equation and using the
built-in Semiconductor module in COMSOL. Since all the lines have the same values everywhere,
the orange and black lines are not visible.
Simulations 51

3
Vapp = +2 V
2.5 Vapp = +1 V
Vapp = 0 V
2 Vapp = -1 V
energy EC [eV]

Vapp = -2 V
1.5

0.5

-0.5
0 1 2 3 4 5 6 7
position z [nm]
Figure 4.32: EC (z) - Conduction band energy across various positions of the VCM bi-layer
structure after applying pulsed voltage for five seconds. The black line denotes the case without
external voltage, while orange and blue lines represent positive and negative pulse voltages,
respectively. Results are obtained from solving the ion continuity equation and employing the
built-in Semiconductor Module in COMSOL. Since all the lines have the same values everywhere,
the orange and black lines are not visible.
Simulations 52

4.6 Phase II-b

After observing no significant change in the distribution of charge carriers following the coupling of
the oxygen vacancy equation and Semiconductor Module in the previous phase of the simulation, it
was decided to utilize the values obtained in phases I-a and I-c for VÖ (z) in Sections 4.1 and 4.3.
subsequently, the oxygen vacancy equation was removed from the simulation, akin to Section 4.2.
However, in this phase, the Semiconductor Module was employed to determine the JV-characteristics
after applying a sweep voltage to the top electrode. Figure 4.33 illustrates the applied voltage to the
VCM layers at different times during this phase of the simulation.

3
applied vooltage Vapp [V ]

-3

0 3 6 9 12
time t [s]
Figure 4.33: Vapp (t) - Applied voltage on the top electrode in the time evolution for phase II-b
of the simulation.

4.6.1 Current-voltage characteristics

The critical variable to examine after subjecting the VCM model to sweep voltage simulation is the
current density across different applied voltages, known as the JV characteristics. Although the
Semiconductor Module computes various current densities, it is beneficial to focus on observing the
drift electron current density along the z-axis, particularly at the midpoint of the second oxide layer.
By analyzing these specific current results, one can compare them with the outcomes obtained in
Section 4.2.
Despite numerous attempts, we were unable to successfully run the simulation using the built-in
Module of the COMSOL software. The primary challenge arose from incorporating direct tunneling
into the model, as convergence issues occurred consistently whenever attempts were made to include
it. Even efforts to reconstruct the model based on a previously tested two-dimensional model for
heterojunction tunneling [29] from the COMSOL website proved ineffective.
Furthermore, excluding direct tunneling from the model and applying sweep voltage to the top elec-
trode did not yield the expected behavior in terms of current inside the device. Figure 4.34 illustrates
an example of the current-voltage characteristics obtained from applying a sweep voltage to a device
previously subjected to a 1 s pulse bias. As depicted in the plot, the changes in current throughout
the voltage evolution do not align with expected outcomes.
Simulations 53

10 0
current density J [A/m2]

10 -20

10 -40

10 -60

10 -80
-3 -2 -1 0 1 2 3
applied voltage Vapp [V]
Figure 4.34: J(Vapp ) - Current density under different applied voltages, obtained using the
Semiconductor Module, for a device previously subjected to a 1 s pulse with a bias of +1 V.

4.7 Phase III-a

In this simulation phase, the objective is to replicate the experiments conducted in Section 4.4, with
the materials and geometry utilized by Aussen et al [27]. Consequently, while the methodologies and
physical model remain consistent with the previously implemented projects, adjustments are made to
the simulation parameters to align with the new materials (AlOx and TiOx ), as detailed in Table 4.2.

4.7.1 Reaching relaxed state

The initial step involves seeking the relaxed-state (or steady-state) distribution of doubly ionized
oxygen vacancies. This can be accomplished by removing the device from biasing and allowing the
simulation to run for an extended duration until a stable distribution is attained.
Figure 4.35 illustrates the evolution of oxygen vacancy distribution across various sections of the
VCM device over time. After five minutes of simulation, a noticeable change persists in the oxygen
vacancy concentration, indicating the necessity for longer simulation durations to achieve a steady-
state. Figure 4.36 illustrates the oxygen vacancy concentration distribution after allowing the device
to remain unbiased for up to 20 minutes. As depicted in the plot, with increasing time, the doubly
ionized vacancies tend to migrate from areas corresponding to the oxide 2 region towards regions near
the top electrode (oxide 1 areas). After a period (around 12 min), the distribution stabilizes, reaching
a steady-state where further increases in time do not significantly alter the vacancy concentration.
This distribution can be regarded as the relaxed-state oxygen vacancy concentration for subsequent
analyses.
Comparison with results obtained in Section 4.4 reveals a longer time required to reach the steady-
Simulations 54

Parameter Symbol Value Unit

oxide 1 thickness toxide1 1.2 nm
oxide 2 thickness toxide2 7 nm
electron affinity of oxide 1 χ1 2.58 eV
electron affinity of oxide 2 χ2 3.5 eV
relative permittivity of oxide 1 ϵr1 2.2
relative permittivity of oxide 2 ϵr2 8
thermal conductivity of oxide 1 κ1 1 W/(m ·K)
thermal conductivity of oxide 2 κ2 1 W/(m ·K)
electron mobility of oxide 1 µn1 5 cm2 /(V · s)
electron mobility of oxide 2 µn2 10 cm2 /(V · s)
effective electron mass for oxide 1 meff,n1 0.25 · me kg
effective electron mass for oxide 2 meff,n2 0.25 · me kg
initial density of donor dopants for oxide 1 ND1 1 ·1016 1/cm3
initial density of donor dopants for oxide 2 ND2 1 ·1020 1/cm3
−7
diffusion coefficient pre-factor for oxide 1 D01 1 ·10 cm2 /s
−7
diffusion coefficient pre-factor for oxide 2 D02 1 ·10 cm2 /s
oxygen ion activation energy for oxide 1 WA1 0.5 eV
oxygen ion activation energy for oxide 2 WA2 0.5 eV
top electrode metal work function (Pt) WM,top 5.84 eV
bottom electrode metal work function (Cr) WM,bottom 4.5 eV
Table 4.2: Parameters used in simulation from [27] (me denotes the electron mass).

state for the new device with the set parameters. this difference might be due to the fact that from
the current device having a longer geometry, resulting in a lengthier diffusion process required for
saturation.

4.7.2 Oxygen vacancy concentration and resistance state

After identifying the relaxed-state distribution of doubly ionized oxygen vacancies across the device,
we attempted to apply positive and negative voltages to the device to observe changes in both oxygen
vancancies and resistance of the VCM Bi-Layer structure. However, numerical errors (non convergings)
prevented us from obtaining usable results. These failures could be due to the high levels of oxygen
vacancy concentration near the top electrode, as shown in Figure 4.36. We could not find a solution
to address this issue while trying to achieve the relaxed-state in Subsection 4.7.1.

4.7.3 Effect of electrode metal material

In this phase of the simulations, all parameters listed in Table 4.2 remained constant. The only
alteration was the substitution of the material for the bottom electrode with Platinum, matching it
with the top electrode. To implement this change, the metal work function of the bottom electrode
was adjusted to 5.84 eV, aligning it with the top electrode.
The subsequent step involved letting the simulation run without applying any bias to the device while
monitoring changes in oxygen vacancy concentration to establish a relaxed-state distribution. Figure
4.37 illustrates the oxygen vacancy concentration across different segments of the device under the
absence of external bias. Over time, the vacancies tend to accumulate near the electrodes and move
away from the oxide regions.
Comparing the recent findings with those outlined in Subsection 4.7.1 yields two conclusions. Firstly,
employing two metals with identical work functions for the top and bottom electrodes leads to a
Simulations 55

10 30
t sim = 0 min
t sim = 1 min
t sim = 2 min
10 28 t sim = 3 min
cocentration VOB [1=m3 ]

t sim = 4 min
t sim = 5 min

10 26

10 24

10 22
0 1 2 3 4 5 6 7 8
position z [nm]
Figure 4.35: VÖ (z) - Distribution of oxygen vacancies across the VCM layers after no voltage is
applied to the top electrode. For 0 min, 1 min, 2 min, 3 min, 4 min, and 5 min of simulation.

symmetric distribution of oxygen vacancies in the relaxed-state. This contrasts with the scenario
where different materials are used, where vacancies tend to accumulate primarily around the electrode
with the higher metal work function. Secondly, in the simulations conducted in this subsection, the
relaxed-state is achieved more rapidly due to the reduced need for vacancy movement.
Simulations 56

10 30

10 25
cocentration VOB [1=m3 ]

10 20

10 15 t sim = 0 min
t sim = 4 min
t sim = 8 min

10 10 t sim = 12 min
t sim = 16 min
t sim = 20 min

10 5
0 1 2 3 4 5 6 7 8
position z [nm]
Figure 4.36: VÖ (z) - Distribution of oxygen vacancies across the VCM layers after no voltage is
applied to the top electrode. For 0 min, 4 min, 8 min, 12 min, 16 min, and 20 min of simulation.

10 28

10 27

10 26
cocentration VOB [1=m3 ]

10 25

10 24
t sim = 0 min
t sim = 1 min
10 23 t sim = 2 min
t sim = 3 min
22 t sim = 4 min
10
t sim = 5 min

10 21
0 1 2 3 4 5 6 7 8
position z [nm]
Figure 4.37: VÖ (z) - Distribution of oxygen vacancies across the VCM layers after no voltage is
applied to the top electrode. For 0 min, 1 min, 2 min, 3 min, 4 min, and 5 min of simulation.
When Platinum is used for bottom electrode material.
Simulations 57

4.8 Discussion
The results obtained from simulating the VCM model with the implemented model demonstrate the
expected changes in internal properties of the oxide layer upon the application of external voltages.
Specifically, when external voltages are applied in the form of pulses, the alterations within the
geometry are contingent upon the voltage magnitude, polarity, and duration of application. Of utmost
importance among these changes is the variation in oxygen vacancy concentration. It is observed that
positive voltages lead to a decrease in distribution within the oxide-1 regions, with a more pronounced
effect as the voltage magnitude increases. Conversely, negatively biased pulses exhibit a reverse effect,
increasing the concentration in oxide-1 areas, with the increment becoming more evident with larger
negative voltages. Subsequently, an analysis of the current characteristics reveals that alterations in
oxygen vacancy distribution can significantly impact the conductivity of the VCM cell. Configurations
with higher oxygen vacancy distribution in oxide-1 regions exhibit greater electrical conductivity
compared to those with lower vacancy levels, suggesting that changes in resistance state in a bi-layer
VCM cell can be attributed to oxygen migration between the layers. Table 4.3 reports the behavior
seen from applying 1 s pulses to the top electrode.

Table 4.3: Changes observed from applying 1 s pulses to the top electrode.
with positive bias with negative bias
VÖ in oxide 1 lower higher
VÖ in oxide 2 higher lower
EC in oxide 1 no meaningful change no meaningful change
EC in oxide 2 no meaningful change no meaningful change
J in positive Vapp lower higher
J in negative Vapp lower higher

Furthermore, the impact of pulse duration is investigated in the subsequent analysis. It is observed
that as the simulation time increases, the dominance of oxygen vacancy diffusion becomes prevalent
in most cases, resulting in an increase in VÖ distribution within oxide-1 regions. This trend persists
across different applied voltage values except for cases involving +2 V. Additionally, in longer pulse
durations, a change in conduction band energy is observed, wherein positive pulses lead to a decrease
in oxide-2 and negative pulses cause a rise in oxide-2 region, with the extent of change increasing with
simulation duration, also in 0 V the conduction band energy increases in oxide 2 areas with increasing
the duration. Table 4.4 reports the behavior seen from applying longer pulses to the top electrode. To
delve deeper into oxygen vacancy behavior, the simulation time is extended to achieve the equilibrium
distribution of vacancies, which is attained after three minutes. Subsequent simulations are then
conducted with the oxygen vacancy concentration set to the equilibrium state at the outset, followed
by the application of a positive voltage. Notably, applying a positive voltage to the top electrode
results in a reduction in oxygen vacancies within oxide-1 regions, shifting them towards the second
oxide. This effect is only feasible with positive voltages exceeding 2 V, with a faster decrease in the
oxygen vacancy concentration occurring as the magnitude of the positive bias increases. Table 4.5
presents the duration required to attain the initial value of the oxygen vacancy distribution (doping)
under the influence of a positive voltage applied to the top electrode, with the initial distribution set
to the steady-state configuration.
To investigate the non-linear change observed in the mean value of oxygen vacancy concentration in
oxide-1 over time, the evolution of the oxygen vacancy profile was examined at different time points
starting from the relaxed state. It was observed that after the peak value of the VÖ (z) distribution
Simulations 58

Table 4.4: Changes observed in different parameters when the duration of applied pulse is increased,
with respect to the cases with a pulse duration of 1 s.
5 s duration 10 s duration
EC in oxide 2 region higher higher
VÖ in oxide 1 in Vapp = 0 V higher higher
VÖ in oxide 1 in Vapp < 0 V higher higher
VÖ in oxide 1 in Vapp = +1 V higher higher
VÖ in oxide 1 in Vapp = +2 V lower lower

Table 4.5: Time taken to reach the initial VÖ mean value in oxide-1 for different positive voltages
when the starting concentration is the steady-state distribution.
Applied voltage [V] Time to reach initial value [s]
+1 did not reach after 100
+2 37
+3 6

reached the interface, the decrease in vacancy concentration mean value became more pronounced.
Subsequently, the resistance state corresponding to each distribution of oxygen vacancy concentration
after applying positive voltage to the steady-state device was measured. After the point at which
the peak value of vacancies reached the oxide 1-oxide 2 interface, the decrease in resistance slowed
down and tended to saturate. This suggests a relation between oxygen vacancy dynamics and electrical
conductivity within the device. The complementary effect becomes apparent upon applying a negative
bias to the top electrode, initiating from the previous low resistance distribution. This manifests as an
increase in the oxygen vacancy concentration within the oxide-1 regions, alongside a rise in the device’s
resistance. Table 4.6 provides a summary of the properties associated with the various resistance states
obtained from simulations conducted on the VCM bi-layer device.
Following the examination of pulse and steady applied biases using a simulation constructed with the
outlined physical model, an attempt was made to replicate the experiment using a built-in module
within the COMSOL Multiphysics software, namely the Semiconductor Module. However, unlike the
findings in the preceding phases, the alteration in the internal properties of the VCM layers, such as
oxygen vacancy concentration and conduction band energy, is not apparent after applying pulses with
varying polarity and magnitude. These simulations only reveal the effect of pulse duration, with the
reason for such behavior remaining unknown.
After exploring the effects of pulsed biases using the Semiconductor Module, attempts were made to
determine the resistance of the device in various states based on the previous distribution of oxygen
vacancies via applying sweep voltages. However, these attempts were unsuccessful as well.
In the final stages of the simulations, a different set of material parameters was tested using the same
methods as before. These materials were sourced from a previous study, with the goal of examining
how the distribution of oxygen vacancy concentration behaves under various biases. Once again, when
the device was left unbiased for a period, a relaxed-state was observed in the vacancy concentration,
with higher concentrations near the top electrode. The only notable difference was the longer time
required to reach this state, likely due to the device’s longer geometry with the new parameters.
However, when the device was subjected to external bias, the simulations did not converge. As a
diagnostic test, the material used for the bottom electrode was changed to observe its effect. The
results revealed that the distribution of oxygen vacancies in the relaxed-state is influenced by the
metal work functions of the device electrodes. Vacancies tend to migrate towards the electrode with
the higher work function at equilibrium. Additionally, if the same material is used for both the top
Simulations 59

Table 4.6: Resistance states in the VCM ReRAM model with corresponding applied biases, resistances,
and vacancy concentrations in oxide-1. Obtained from simulations with parameters reported in Table
4.1.
Resistance State Applied Bias Resistance [Ω.m2 ] Vacancy Concentration in Oxide-1
HRS (OFF) Negative ≈ 10 × 109 High
LRS (ON) Positive ≈ 1 × 109 Low

and bottom electrodes, the distribution of oxygen vacancies within the device becomes symmetric.
 Chapter 5

Summary and outlook

5.1 Summary
This thesis focuses on 1-D simulations conducted on a ReRAM VCM bi-layer structure composed
of Hafnium oxide and Tantalum oxide. These simulations were performed using the finite element
method with COMSOL Multiphysics software. The thesis outlines the physical model and equations
employed to construct the simulation, including boundary conditions and additional properties.
Two distinct approaches were utilized in the simulation process. Initially, the differential equations
were directly implemented, followed by the utilization of existing built-in physical models within the
software.
Across the implemented phases, the effects of applying pulsed voltages to the device were investigated.
The results indicate that external bias application can induce changes in internal properties of the
stacked layers, notably in the concentration of doubly ionized oxygen vacancies and conduction band
energies. These alterations are dependent on the magnitude, polarity, and duration of the applied
bias, with higher voltages and longer durations magnifying the effects. These changes in material
properties can subsequently influence the conductivity of the device, which is a key feature of ReRAM
technology.
To deepen our understanding of the relationship between oxygen exchange between the layers and
the electrical resistance of the device, efforts were made to achieve a steady-state (relaxed-state)
by applying zero voltage to the device for an extended period. Subsequently, starting from the
relaxed state and applying positive bias, changes in oxygen vacancy concentration and, consequently,
resistance were observed. The simulation was also done in reversed manner, applying negative voltage
to the previously positive-biased device. A distinct relationship was observed between the vacancy
concentration within the oxide 1 (conducting oxide) and the resistance state of the device. In essence,
a higher concentration of oxygen vacancies within the oxide 1 region corresponds to a high resistance
state, while a lower concentration in oxide 1 leads to a low resistance state. This temporal evolution

60
Summary and outlook 61

was evident in the time-dependent results, with initial rapid changes followed by saturation as most
of the oxygen vacancies migrated toward bottom electrode (oxide 2 regions). In summary, applying a
positive voltage of adequate magnitude can be likened to the SET process, while the application of a
sufficiently negative bias initiates the RESET process. Both processes exhibit time dependence.
In contrast, while utilizing the built-in module, attempts to replicate the same pulse measurements
did not yield anticipated results for unknown reasons. In a subsequent phase, efforts were made to
measure the conductivity of the device using the distributions obtained from the implemented phases,
again those simulation were not successful. The reasons behind the received errors are not clear.
In the final stages, we tested different material parameters and observed the behavior of oxygen va-
cancy concentration under various biases. Despite achieving a relaxed-state in vacancy concentration,
convergence issues arose during simulations under external bias. We experimented with altering the
material for the bottom electrode, revealing a dependency on metal work functions. This affected
vacancy migration towards the electrode with the higher work function and led to symmetric vacancy
distribution when using identical materials for both electrodes.

5.2 Outlook
Throughout this thesis, a limited set of parameters were investigated to understand their impact on
the internal and external properties of the bi-layer VCM ReRAM device. However, due to challenges
encountered in numerical modeling, further exploration with a more sophisticated physical model may
be necessary to overcome these simulation difficulties. This advanced model could facilitate the testing
of a wider range of materials and geometries, providing deeper insights into device behavior.
Moreover, the use of the built-in Semiconductor Module within the COMSOL Multiphysics software
did not prove effective for analyzing a 1-D model when ion charge carriers were included, and attempts
to implement direct tunneling were unsuccessful. Therefore, it is recommended that an improved
module be developed to address these issues and enhance the analysis capabilities for such devices.
 Appendix A

Temperature fluctuation results

Figure A.1 illustrates the temperature variations across the VCM layer at different time intervals
during phase I-a simulation with Vapp = 2 V. It is evident from the plot that the temperature changes
are negligible, indicating their insignificance in the simulation results.

300

290
temperature T [K]

t=0s
280
t = 0.1 s
t = 0.2 s
t = 0.5 s
270 t=1s

260

250
0 1 2 3 4 5 6 7
position z [nm]
Figure A.1: T (z) - Temperature in different parts of the VCM layers in different times of the
simulation when +2 V pulse is applied to the top electrode for one second.

62
List of Figures

1.1 Classification of the resistive switching effects which are considered for non-volatile
memory applications. From: [1], license number: 5787050663774. . . . . . . . . . . . . 2

2.1 a) Illustration of unipolar switching, where both Set and Reset processes are unaffected
by the voltage polarity. b) Diagram depicting bipolar switching, where the Set process
occurs at one voltage polarity, and the Reset process takes place at the reverse polarity.
Adapted from: [1], license number: 5787050663774. . . . . . . . . . . . . . . . . . . . . 4
2.2 Sketch of the lattice cubic structure of SrTiO3 . . . . . . . . . . . . . . . . . . . . . . . 5
2.3 Energy-band diagram of a metal/n-type semiconductor contact. . . . . . . . . . . . . . 6
2.4 Schematic of the potential energy landscape for ion hopping at a) absence of external
force b) when applying voltage . Redrawn form [17] © 2012 IEEE. . . . . . . . . . . . 7
2.5 Transmission probability in different energy levels of the electrons in an energy barrier
based on the classical and quantum mechanics equations. On left and on the right the
amount of the transmission is different. This is due to the fact that thickness of the
barrier is changed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.6 By applying the forming voltage, a filament of oxygen vacancies is built starting from
the Ohomic electrode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.7 The schematic representation illustrates the Set and Reset characteristics of a Valence
Change Memory (VCM) device. In the OFF state (A), there is no conducting fila-
ment between the electrodes. Upon applying a negative voltage, ions are attracted
to the top electrode, resulting in the creation of a filament during the Set process
(B) and leading to the ON state (C) with a conducting filament between both elec-
trodes. Conversely, applying a positive voltage to the active electrode during the Reset
process (D) prompts the migration of oxygen ions back into the plug, resulting in a
measurable reset. The device then returns to the OFF state (A). This depiction is
derived from [13] (https://marketplace.copyright.com/rs-ui-web/mp/license/4df5fdc4-
11c2-46d8-a901-d8de4401c120/1bc1ebd0-c59a-4b64-a722-a258b16fb766). . . . . . . . . 10

63
LIST OF FIGURES 64

2.8 Sketch of complementary switching. In both voltage polarities HRS and LRS are
present. Applying a positive voltage results in Set 1 then Reset 1. Applying nega-
tive voltage leads to Set 2 followed by Reset 2. . . . . . . . . . . . . . . . . . . . . . . 11

3.1 Schematic of 1D model used in simulation of bi-layer VCM cell . . . . . . . . . . . . . 13

3.2 Band diagram at a Schottky contact, the Schottky potential of the interface is the green
dotted line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.3 Position and value of the parameters used as boundary conditions in Poisson equation. 17
3.4 Sketch of tunneling mechanism through the potential barrier at the metal-oxide interface. 18
3.5 Schematic of conduction band energy at oxide-oxide interface. . . . . . . . . . . . . . . 19
3.6 Position and value of the parameters used as boundary conditions in electron continuity
equation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.7 Position and value of the parameters used as boundary conditions in ion drift-diffusion
equation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.8 Position and value of the parameters used as boundary conditions in heat conduction
equation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

4.1 Vapp (t) - Applied voltage function for simulation phase I-a. The increase and decrease
of the voltage happens in a short time (10−5 s). . . . . . . . . . . . . . . . . . . . . . . 24
4.2 VÖ (z) - Oxygen vacancy concentration along the bi-layer ReRAM cell. The black line
represents the condition without applied voltage, the orange lines depict the scenario
with positive voltage application, and the blue lines illustrate the effect of negative
voltage pulses. Due to the simulation geometry, z < 2 nm corresponds to oxide 1
(HfOx ), and z > 2 nm represents oxide 2 (TaOx ) areas. . . . . . . . . . . . . . . . . . . 25
4.3 VÖ (z) - Oxygen vacancy concentration along the bi-layer ReRAM cell. The diagram
is zoomed in to observe details near and inside the oxide 1 regions. The black line
represents the equilibrium case (without external voltages), and the orange lines depict
the oxygen vacancy distribution after applying different amounts of positive voltage
pulses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.4 VÖ (z) - Oxygen vacancy concentration along the bi-layer ReRAM cell. The diagram
compares different amplitudes of applied negative voltage pulses (blue lines) with the
equilibrium case (black line). The diagram concentrates on areas of oxide 1. . . . . . . 27
4.5 EC (z) - Conduction band energy along the bi-layer VCM cell. The black line represents
the energy without the pulsed voltage, orange lines represent the situation after applying
positive voltage pulses with no profound change in EC , and the green line represents
the energy for the voltage pulse that has an effect on the conduction band energy level. 27
4.6 EC (z) - Conduction band energy along the bi-layer VCM cell. The black line represents
the energy without the pulsed voltage, blue lines represent the situation after applying
negative voltage pulses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.7 Vapp (t) - Applied voltage on the top electrode in the time evolution for phase I-b of the
simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.8 J(z) - Current density distribution at different Vapp , where z = 4.5 nm (dashed line)
represents the middle of oxide 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
LIST OF FIGURES 65

4.9 JV characteristics of the VCM cell under a) negative voltages and b) positive voltages.
The black lines represent the case where no voltage pulse is applied, the blue line shows
the current for the case where a −1 V pulse is applied, and the +1 V pulse case is
depicted with the orange line (all the currents are measured in z = 4.5 nm). . . . . . . 31
4.10 JV characteristics of the VCM cell under a) negative voltages and b) positive voltages.
The black lines represent the case where no voltage pulse is applied, the orange lines
show the current for the cases where a positive pulse is applied (all the currents are
measured in z = 4.5 nm). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.11 JV characteristics of the VCM cell under a) negative voltages and b) positive voltages.
The black lines represent the case where no voltage pulse is applied, the blue lines show
the current for the cases where a positive pulse is applied (all the currents are measured
in z = 4.5 nm). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.12 Vapp (t) - Applied voltage function for simulation phase I-c. The increase and decrease
of the voltage happens in a short time (10−5 s). . . . . . . . . . . . . . . . . . . . . . . 34
4.13 Variation of oxygen vacancy concentration (VÖ (z)) across different positions of the
VCM structure under an external voltage pulse of 0 V, applied for durations of 1 s
(black line), 5 s (green line), and 10 s (purple line). . . . . . . . . . . . . . . . . . . . . 35
4.14 Variation of oxygen vacancy concentration (VÖ (z)) across different positions of the
VCM structure under an external voltage pulse of +1 V, applied for durations of 1 s
(black line), 5 s (green line), and 10 s (purple line). . . . . . . . . . . . . . . . . . . . . 36
4.15 Variation of oxygen vacancy concentration (VÖ (z)) across different positions of the
VCM structure under an external voltage pulse of +2 V, applied for durations of 1 s
(black line), 5 s (green line), and 10 s (purple line). . . . . . . . . . . . . . . . . . . . . 37
4.16 Variation of oxygen vacancy concentration (VÖ (z)) across different positions of the
VCM structure under an external voltage pulse of −1 V, applied for durations of 1 s
(black line), 5 s (green line), and 10 s (purple line). . . . . . . . . . . . . . . . . . . . . 37
4.17 Variation of oxygen vacancy concentration (VÖ (z)) across different positions of the
VCM structure under an external voltage pulse of −2 V, applied for durations of 1 s
(black line), 5 s (green line), and 10 s (purple line). . . . . . . . . . . . . . . . . . . . . 38
4.18 EC (z) - Conduction band energy in various positions of the VCM bi-layer structure at
different simulation durations without external bias. The black line represents the 1 s
simulation time, the green line corresponds to 5 s duration, and the purple line depicts
the 10 s case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.19 EC (z) - Conduction band energy in various positions of the VCM bi-layer structure at 5
s duration of pulse applied voltage. The black line represents the case without external
voltage, the orange lines correspond to positive voltages, and the blue lines depict the
case when the applied pulses have negative polarity. . . . . . . . . . . . . . . . . . . . 39
4.20 EC (z) - Conduction band energy in various positions of the VCM bi-layer structure
at 10 s duration of pulse applied voltage. The black line represents the case without
external voltage, the orange lines correspond to positive voltages, and the blue lines
depict the case when the applied pulses have negative polarity. . . . . . . . . . . . . . 39
4.21 VÖ (z) - Distribution of oxygen vacancies across the VCM layers after no voltage is
applied to the top electrode. For 0 s, 10 s, 20 s, 30 s, 40 s, 50 s and 60 s of simulation. 41
LIST OF FIGURES 66

4.22 VÖ (z) - Distribution of oxygen vacancies across the VCM layers after no voltage is
applied to the top electrode. For 1 min, 2 min, 3 min, 4 min, and 5 min of simulation.
The result is stabilized after 3 minutes. . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.23 VÖ,ave (t) - Mean value of the oxygen vacancy concentration in the oxide 1 region at
different time intervals when a positive bias of +2 V is applied to the top electrode,
with the initial distribution set to the relaxed state. . . . . . . . . . . . . . . . . . . . 43
4.24 VÖ,ave (t) - Mean value of the oxygen vacancy concentration in the oxide 1 region at
different times, when voltages of +3 V (orange), +2 V (blue), and +1 V (green) are
applied to the top electrode. The initial distribution corresponds to the relaxed state. 44
4.25 VÖ (z) - Oxygen vacancy concentration in different simulation domains when +2 V is
applied to the top electrode. The red lines represent cases where the drop in the mean
value of oxide 1 is minimal, the black line indicates the point where the maximum
oxygen vacancy concentration is located near the interface, and blue lines correspond
to cases where the decrement of the oxygen vacancy mean value is rapid for oxide 1
regions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.26 R(t) vs. VÖ,ave (t)- Resistance times area measured for a +0.5 V bias in the middle
of oxide 2 (z = 4.5 nm) versus the mean value of the doubly ionized oxygen vacancy
concentration in oxide-1, in time evolution during which a +2 V bias is applied to the
top electrode starting from the relax-state distribution. . . . . . . . . . . . . . . . . . . 45
4.27 R(t) vs. VÖ,ave (t)- Resistance times area measured for a +0.5 V bias in the middle
of oxide 2 (z = 4.5 nm) versus the mean value of the doubly ionized oxygen vacancy
concentration in oxide-1, in time evolution during which a −2 V bias is applied to the
top electrode starting from the LRS distribution. . . . . . . . . . . . . . . . . . . . . . 46
4.28 Vapp (t) - Applied voltage function for simulation phase II-a. The increase and decrease
of the voltage happens in a short time (10−5 s). . . . . . . . . . . . . . . . . . . . . . 47
4.29 VÖ (z) - Oxygen vacancy distribution across different positions of the VCM structure
under external voltage pulses applied for durations of one second. The black line rep-
resents the case without external voltage, while orange and blue lines depict positive
and negative pulse voltages, respectively. Results are derived from solving the ion con-
tinuity equation and using the built-in Semiconductor module in COMSOL. Since all
the lines have the same values everywhere, the orange and black lines are not visible. . 48
4.30 EC (z) - Conduction band energy across various positions of the VCM bi-layer structure
after applying pulsed voltage for one second. The black line denotes the case without
external voltage, while orange and blue lines represent positive and negative pulse
voltages, respectively. Results are obtained from solving the ion continuity equation
and employing the built-in Semiconductor Module in COMSOL. Since all the lines have
the same values everywhere, the orange and black lines are not visible. . . . . . . . . . 49
4.31 VÖ (z) - Oxygen vacancy distribution across different positions of the VCM structure
under external voltage pulses applied for durations of five seconds. The black line
represents the case without external voltage, while orange and blue lines depict positive
and negative pulse voltages, respectively. Results are derived from solving the ion
continuity equation and using the built-in Semiconductor module in COMSOL. Since
all the lines have the same values everywhere, the orange and black lines are not visible. 50
4.32 width=14cm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
LIST OF FIGURES 67

4.33 Vapp (t) - Applied voltage on the top electrode in the time evolution for phase II-b of
the simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.34 J(Vapp ) - Current density under different applied voltages, obtained using the Semi-
conductor Module, for a device previously subjected to a 1 s pulse with a bias of +1
V. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.35 VÖ (z) - Distribution of oxygen vacancies across the VCM layers after no voltage is
applied to the top electrode. For 0 min, 1 min, 2 min, 3 min, 4 min, and 5 min of
simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.36 VÖ (z) - Distribution of oxygen vacancies across the VCM layers after no voltage is
applied to the top electrode. For 0 min, 4 min, 8 min, 12 min, 16 min, and 20 min of
simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.37 VÖ (z) - Distribution of oxygen vacancies across the VCM layers after no voltage is
applied to the top electrode. For 0 min, 1 min, 2 min, 3 min, 4 min, and 5 min of
simulation. When Platinum is used for bottom electrode material. . . . . . . . . . . . 56

A.1 T (z) - Temperature in different parts of the VCM layers in different times of the simu-
lation when +2 V pulse is applied to the top electrode for one second. . . . . . . . . . 62
List of Tables

4.1 Parameters used in simulation (me denotes the electron mass). . . . . . . . . . . . . . 24

4.2 Parameters used in simulation from [27] (me denotes the electron mass). . . . . . . . . 54
4.3 Changes observed from applying 1 s pulses to the top electrode. . . . . . . . . . . . . . 57
4.4 Changes observed in different parameters when the duration of applied pulse is in-
creased, with respect to the cases with a pulse duration of 1 s. . . . . . . . . . . . . . 58
4.5 Time taken to reach the initial VÖ mean value in oxide-1 for different positive voltages
when the starting concentration is the steady-state distribution. . . . . . . . . . . . . . 58
4.6 Resistance states in the VCM ReRAM model with corresponding applied biases, resis-
tances, and vacancy concentrations in oxide-1. Obtained from simulations with param-
eters reported in Table 4.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

68
 Bibliography

[1] Rainer Waser, Regina Dittmann, Georgi Staikov, and Kristof Szot. Redox-Based Resistive Switch-
ing Memories – Nanoionic Mechanisms, Prospects, and Challenges. Advanced Materials, 21(25-
26):2632–2663, 2009.

[2] J. Joshua Yang, Dmitri B. Strukov, and Duncan R. Stewart. Memristive devices for computing.
Nature Nanotechnology, 8(1):13–24, 2012.

[3] L. Chua. Memristor-The missing circuit element. IEEE Transactions on Circuit Theory,
18(5):507–519, 1971.

[4] Leon Chua. Resistance Switching Memories Are Memristors, pages 21–51. Springer International
Publishing, 2014.

[5] A. V. Fadeev and K. V. Rudenko. To the Issue of the Memristor’s HRS and LRS States Degra-
dation and Data Retention Time. Russian Microelectron, 50(5):311–325, 2021.

[6] Matthias Wuttig and Noboru Yamada. Phase-change materials for rewriteable data storage.
Nature Material, 6(11):824–832, 2007.

[7] Michael N. Kozicki, Maria Mitkova, and Ilia Valov. Electrochemical Metallization Memories. In
Resistive Switching, pages 483–514. John Wiley Sons, Ltd, 2016.

[8] M. N. Kozicki, M. Yun, L. Hilt, and A. Singh. Applications of programmable resistance changes
in metal-doped chalcogenides. In Electrochemical Society, MEETING ABSTRACTS- ELEC-
TROCHEMICAL SOCIETY -ALL DIVISIONS-, number 1, page 849. Electrochemical Society;,
1999.

[9] Stephan Menzel, Matthias Waters, Astrid Marchewka, Ulrich Böttger, Regina Dittmann, and
Rainer Waser. Origin of the Ultra-nonlinear Switching Kinetics in Oxide-Based Resistive
Switches. Advanced Functional Materials, 21(23):4487–4492, 2011.

[10] Chun-Yang Huang, Chung-Yu Huang, Tsung-Ling Tsai, Chun-An Lin, and Tseung-Yuen Tseng.
Switching mechanism of double forming process phenomenon in ZrOx/HfOy bilayer resistive

69
BIBLIOGRAPHY 70

switching memory structure with large endurance. Applied Physics Letters, 104(6):062901, 02
2014.

[11] Fenja Berg. Simulation, fabrication and electrical characterisation of ReRAM bilayer structures:
Understanding the influence of stacked oxides on forming and switching behaviour of ReRAM.
Master’s thesis, RWTH Aachen, 2017.

[12] Astrid Marchewka. A numerical simulation model of valence-change-based resistive switching.

PhD thesis, RWTH Aachen, 2017.

[13] R. Waser. Nanoelectronics and Information Technology. Wiley-VCH, 2012.

[14] Benedikt Arndt, Francesco Borgatti, Francesco Offi, Monifa Phillips, Pedro Parreira, Thorsten
Meiners, Stephan Menzel, Katharina Skaja, Giancarlo Panaccione, Donald A. MacLaren, Rainer
Waser, and Regina Dittmann. Spectroscopic indications of tunnel barrier charging as the switch-
ing mechanism in memristive devices. Advanced Functional Materials, 27(45):1702282, 2017.

[15] Lesly E. Smart and Eliane A. Moore. Solid State Chemistry: An Introduction. Taylor & Francis,
1995.

[16] S.M. Sze and Kwok K. Ng. Metal-Semiconductor Contacts. In Physics of Semiconductor Devices,
pages 134–196. John Wiley Sons, Ltd, 2006.

[17] Stefano Larentis, Federico Nardi, Simone Balatti, David C. Gilmer, and Daniele Ielmini. Resis-
tive Switching by Voltage-Driven Ion Migration in Bipolar RRAM—Part II: Modeling. IEEE
Transactions on Electron Devices, 59(9):2468–2475, 2012.

[18] R. Tsu and L. Esaki. Tunneling in a finite superlattice. Applied Physics Letters, 22(11):562–564,
06 1973.

[19] Mario Lanza. A Review on Resistive Switching in High-k Dielectrics: A Nanoscale Point of View
Using Conductive Atomic Force Microscope. Materials, 7:2155–2182, 2014.

[20] F. Nardi, S. Balatti, S. Larentis, and D. Ielmini. Complementary switching in metal oxides: To-
ward diode-less crossbar RRAMs. In 2011 International Electron Devices Meeting, pages 31.1.1–
31.1.4, 2011.

[21] Myoung-Jae Lee, Chang Bum Lee, Dongsoo Lee, Seung Ryul Lee, Man Chang, Ji Hyun Hur,
Young-Bae Kim, Chang-Jung Kim, David H. Seo, Sunae Seo, U-In Chung, In-Kyeong Yoo, and
Kinam Kim. A fast, high-endurance and scalable non-volatile memory device made from asym-
metric Ta2O5x/TaO2x bilayer structures. Nature Materials, 10(8):625–630, 2011.

[22] Umesh Chand, Chun-Yang Huang, and Tseung Yuen Tseng. Mechanism of High Temperature
Retention Property (up to 200 °C) in ZrO2-Based Memory Device With Inserting a ZnO Thin
Layer. IEEE Electron Device Letters, 35:1019–1021, 2014.

[23] Zhongrui Wang, HongYu Yu, Xuan Anh Tran, Zheng Fang, Jinghao Wang, and Haibin Su.
Transport properties of HfO2−x based resistive-switching memories. Phys. Rev. B, 85:195322,
May 2012.

[24] Bo Xiao and Satoshi Watanabe. Oxygen vacancy effects on an amorphous-TaOx-based resistance
switch: a first principles study. Nanoscale, 6(17):10169–10178, 2014.
BIBLIOGRAPHY 71

[25] COMSOL AB, Stockholm, Sweden. Comsol multiphysics® v.6.1. https://www.comsol.com, 2024.

[26] A. Gehring. Simulation of Tunneling in Semiconductor Devices. PhD thesis, Technische Univer-
sität Wien, 2003.

[27] Stephan Aussen, Felix Cüppers, Carsten Funck, Janghyun Jo, Stephan Werner, Christoph
Pratsch, Stephan Menzel, Regina Dittmann, Rafal Dunin-Borkowski, Rainer Waser, and Susanne
Hoffmann-Eifert. Correlation between Electronic Structure, Microstructure, and Switching Mode
in Valence Change Mechanism Al2O3/TiOx-Based Memristive Devices. Advanced Electronic
Materials, 9(12):2300520, 2023.

[28] COMSOL AB, Stockholm, Sweden. Semiconductor module user’s guide. COMSOL Multi-
physics® v.6.1, 2023.

[29] COMSOL. Heterojunction tunneling. https://www.comsol.com/model/

heterojunction-tunneling-65761, 2023. Accessed on March, 2024.