Computer simulations of swift heavy ion effects in graphene and amorphous bulk materials

(1)

HU-P-D277

Computer simulations of swift heavy ion eﬀects in graphene and amorphous bulk

materials

Henrique Vázquez Muíños

Division of Materials Physics Department of Physics

Faculty of Science University of Helsinki

Helsinki, Finland

ACADEMIC DISSERTATION

To be presented for public examination with the permission of the Faculty of Science of the University of Helsinki, in the auditorium PIII of Porthania, on the 12th of March

2022 at 12 o’clock.

HELSINKI 2022

(2)

Helsinki 2022 Unknown

ISBN 978-951-51-7822-0 (PDF version) Helsinki 2022

(3)

Henrique Vázquez Muíños, Computer simulations of swift heavy ion eﬀects in graphene and amorphous bulk materials, University of Helsinki, 2022, 63 p. + appendices, University of Helsinki Report Series in Physics, (printed version), (PDF version)

ABSTRACT

Ion irradiation is capable of modifying the structure of materials. By changing the structure in a controlled fashion, it is possible to tune the material properties to create new devices. One of the most exciting applications of ion irradiation can be found in the semiconductor industry, where this technique is used to manufacture electronic components.

Swift Heavy Ions (SHI) are a speciﬁc type of ion irradiation characterized by their high mass, and elevated energies (> 1 MeV/amu). These ions interact mainly with the electrons in the material, and often leave traces of cylindrical defected regions known as ion tracks. These defects are few nanometers wide but can be up to several microns long.

SHI irradiation can be used to produce ion tracks of tunable size, making it a suitable technique for several promising applications. By analyzing the abundance of tracks one can estimate the age of geological and archaeological samples. The etched ion tracks in silica (a-SiO₂) serve as excellent templates to grow arrays of metallic nanowires. On the other hand, SHI irradiation of graphene allows also the engineering of high eﬃciency water desalinators.

Experimental measurements have shown that SHIs can produce defects in graphene, and that their size grows with the ion energies. Previous reports have identiﬁed a core-shell structure within the ion tracks of several amorphous materials: a-SiO₂, a-Si₃N₄ and a- Si. Nevertheless, current experimental techniques are not able to resolve the structure of these defects, nor the mechanism of their formation. In this study, we use Molecular Dynamics, augmented by the two-temperature model and a Monte Carlo-based electronic cascade model to study the formation mechanism and structure of defects at the atomic level, as well as to investigate the early electron dynamics after the ion impact in these materials.

Our simulations show that SHI irradiation can create pore-like defects in graphene; even when a considerable amount of the energy initially deposited by the ion is removed via electron emission during the development of the electron cascade. Moreover, we attribute the formation of the core-shell structure in a-SiO₂, a-Si₃N₄ and a-Si to the diﬀerent transient pressure levels at the time when the core and the shell of the track solidify; higher (lower) pressures leading to higher (lower) densities, respectively. This work constitutes a step forward at modeling the interaction of SHIs in 2D and amorphous materials, and sheds light on the mechanisms of defects formation in these type of materials under extreme conditions of high-energy ion irradiation.

(4)

1 INTRODUCTION

Ions are atoms that are electrically charged, i.e. they have different number of protons than electrons. Particle accelerators can generate artificially streams of high speed ions, known as ion beams. Ion beams have many applications, ranging from medicine, where they are used for cancer treatment [1–3], to imaging and material analysis, where they are employed to examine the structure of the material with great detail [4]. They are also used to implant ions in materials and modify their properties; this technique is employed by the semiconductor industry to produce most of the electric devices in our every day life. Ions can be found in space in form of cosmic rays or inside nuclear power plants as a byproduct of the fusion or fission reactions. Improved understanding of how these particles damage materials could enable the development of more resistant materials for their use in satellite components and power plants.

Ions interact diﬀerently with the material depending on their mass and energy. Low energy ions tend to interact mainly with the atomic nuclei of the target atoms, displacing them and creating defects in the material. High energy ions interact instead with the electrons in the target, exciting them to higher energy levels, and often causing disordered regions in the material. Low energy ions deposit the energy mainly to the atoms, whereas high energy ions deposit it to the electrons, therefore, these are often referred as nuclear and electronic energy deposition regimes. The nuclear regime has been studied to great detail since 1950s [5–7]. On the other hand, the electronic regime has been studied to a lesser extent, and remains only partially understood. In this dissertation we focus on a speciﬁc type of ions belonging to the latter regime known as Swift Heavy Ions (SHIs).

This term refers to ions with large masses (heavier than carbon) and high energies (in the MeV range).

SHIs are characterized by their capability of travelling long distances inside the material while interacting strongly with the electrons in it. In some materials, this interaction can result in the formation of cylindrical defects known as ion tracks, which tend to be several microns long but only few nanometers wide. SHI tracks were ﬁrst observed in insulating crystals in the late 1950s [8, 9]; however, with the development of new experimental techniques it was later discovered that they can form also in semiconductors and amorphous materials.

The small distances—thousand times smaller than a human cell—and short time scales—

trillion times faster than the blink of an eye—of the processes happening inside the ion track do not allow the exact mechanism of how ions damage the material to be discerned in experiments. In that respect, Molecular Dynamics (MD) simulations cover the desired space and time scales of this problem. This technique follows the trajectory of individual atoms in the material and enables the identiﬁcation of defects at the atomic level. We use this technique in this thesis to investigate the behaviour of the materials to the radiation.

(8)

In this dissertation we study the mechanisms of how the electron excitation induced by the SHIs causes modiﬁcations in the materials. We focus on two groups of materials:

amorphous and two-dimensional materials. SHI tracks in amorphous materials could have application in the manufacture of optical waveguides [10, 11]. Moreover, the damaged region in the track can be chemically etched, i.e. removed with a solvent that reacts chemically with the disordered structure only; track-etched material membranes could be used to build nanoﬂuid channels for directional transport of chemicals [12–16] or to manufacture nanowire matrices in conjunction with electrochemical deposition [17–19].

In two-dimensional materials, SHIs could be used to create pores of few nanometers in the layer. For example, the precise perforation of graphene would allow the engineering of DNA sequencers [20], water desalinators [21–24] of outstanding performance or nanoﬁlters of extraordinary sensitivity [25–28].

In publicationsI and II we focus on the defects created by SHIs in graphene, a two- dimensional material. In publicationIwe investigate with MD simulations the formation of pore-like defects in the graphene layer during SHI impacts and compare the results with experiments. In publicationII, we look deeper into the electron dynamics in this material during the ion impact in order to improve the modelling of two-dimensional materials. In that publication we employ a Monte Carlo (MC) model for transport of energetic electrons, Particle-in-Cell (PIC) and Time Dependent Density Functional Theory (TDDFT) techniques to study the ion-induced electron excitation in the material.

In publicationIII,IVwe investigate with MD simulations the structural transformations that lead to the formation of the core-shell ﬁne structure within the SHI ion tracks of a-SiO₂, a-Si₃N₄ and a-Si. In publicationIIIwe focus the analysis on the defects of the damaged regions and compare against experiments. In publicationIVwe investigate the mechanism responsible for the formation of the core-shell structure, linking the transformations with the transient temperature and pressure conditions inside the ion track. In publicationVwe study with MD simulations the formation mechanism of SHI tracks in a-C, and look at the structural changes in the damaged region within the track.

(9)

2 PURPOSE AND STRUCTURE

The purpose of this thesis is to provide insight at the atomic level into the mechanisms how SHIs damage materials. This study has a two-fold focus. First, simulate with state of the art models SHI impacts in graphene and determine if the observed experimental increase of Raman signal after irradiation is due to the formation of nanopores in the material. Second, resolve the formation mechanism and inner structure of the SHI tracks observed with small angle X-ray scattering in a-SiO₂, a-Si₃N₄, a-Si and a-C amorphous materials.

This dissertation is composed of ﬁve research articles, all of them published in peer- reviewed journals. The publications are listed in this chapter, together with a brief description. Chapter three introduces the key processes involved during the impact of a SHI. Chapter four provides an overview of previous irradiation studies of the materials in the scope of this thesis. Chapter ﬁve describes in detail the methods employed in the publications. Chapter six focuses on the early electron dynamics during the impact of SHIs in graphene, and on the resulting defects. Chapter seven investigates the structural changes within the ion tracks in a-SiO₂, a-Si₃N₄ and a-C materials, and the mechanism how these changes take place. Finally, the conclusions are given in chapter eight.

2.1 Summaries of the original publications

Publication I : Creating nanoporous graphene with swift heavy ions

H. Vázquez, E.H. Åhlgren, O. Ochedowski, A.A. Leino, R. Mirzayev, R. Kozubek, H.

Lebius, M. Karlu˜sic, M. Jak˜sic, A.V. Krasheninnikov and J. Kotakoski Carbon, 114, pp.511-518. (2017)

This publication focuses on the size and morphology of defects formed by SHIs in graphene.

SHI irradiation is a promising technique that could be used to create pores in graphene in a controlled manner for novel application. Up to date, however, these defects can only be measured indirectly via Raman spectroscopy, which does not give information on the type of defects. Our Molecular Dynamics simulations combined with the two-temperature model show the formation of nanopores in graphene. We report an increases of the pore size with the stopping power of the ion; this trend agrees well with the increase of the Raman signal in the experiments and suggests that SHIs can produce pore-like defects in the graphene layer.

Publication II : Electron cascades and secondary electron emission from graphene under ion irradiation

H. Vázquez, A. Kononov, A. Kyritsakis, N. Medvedev , A. Schleife, F. Djurabekova, Physical Review B 103.22 (2021): 224306

(10)

This work investigates the initial dynamics of the electron cascade induced by the ion impact in graphene. We focus our analysis on the electrons emitted during the ion impact as well as the electrons captured by the projectile. We employ a Monte Carlo based model and Time Dependent Density Functional Theory to simulate these processes. The simulations show a large number of electrons emitted during a single ion impact, the amount of energy removed by those electrons accounts between 15-70% of the energy initially deposited by the ion. These results reveal the importance of electron emission during SHI irradiation of 2D materials and suggests that the removal of energy by this process might lead to smaller defects in these materials compared to their bulk counterparts.

Publication III : Nanoscale density variations induced by high energy heavy ions in amorphous silicon nitride and silicon dioxide

P. Mota-Santiago, H. Vázquez, T. Bierschenk, F. Kremer, A. Nadzri, D. Schauries, F.

Djurabekova, K. Nordlund, C. Trautmann, S. Mudie and M.C. RidgwayNanotechnology, 29(14), p.144004 (2018).

This study focuses on the inner structure of SHI tracks in amorphous silicon dioxide and silicon nitride. SHI tracks in these materials are attractive for waveguide applications.

Small Angle X-ray Scattering experiments and Molecular Dynamics simulations show that the ion tracks have a clear inner structure, with an overdense shell and an underdense core.

Both Infrared Absorbtion Spectroscopy and the simulations show the presence of broken bonds in the core of the track. Silicon nitride exhibits three times larger density change than silicon dioxide, this diﬀerence was attributed to the higher thermal conductivity and larger drop in viscosity in the former material.

Publication IV : Ultrafast phase transitions in polyamorphic materials triggered by swift heavy ion impacts

H. Vázquez, F. Djurabekova,Physical Review Materials 5.6 (2021): 065603.

This publication continues the work of publication IV, but it focuses on the mechanism how the core-shell structure forms inside the SHI track. We study three amorphous materials which were shown to exhibit such density structures: silicon dioxide, silicon nitride and silicon. We simulate with Molecular Dynamics and the two-temperature model the ion impact and analyze the density, temperature and pressure inside the ion track. We show that the diﬀerent density values inside the track are consequence of elevated pressures levels inside the track at the moment when the track region solidiﬁes.

High and low pressures lead respectively to higher and lower densities. These results clarify the mechanism how the ﬁne structure of the track is formed in these amorphous materials.

(11)

Publication V : Graphitization of amorphous carbon by swift heavy ion impacts: Molecular dynamics simulation

K. Kupka, A.A. Leino, W. Ren, H. Vázquez, E.H. Åhlgren, K. Nordlund, M. Tomut, C. Trautmann, P. Kluth, M. Toulemonde and F. Djurabekova Diamond and Related Materials, 83, pp.134-140 (2018)

This work studies the structural changes inside the ion tracks in amorphous carbon.

Amorphous carbon is used in accelerator facilities as a charge stripper for ions. It is important to understand how the SHIs create defects in the material, as they aﬀect the eﬃciency and lifetime of the stripper foils. Simulations with Molecular Dynamics coupled with the two-temperature model show that the original sp³ -rich structure turns into sp² inside the ion track. The increase of graphitic carbon in the sample is consistent with experiments which reported an increase of electric conductivity after SHI irradiation.

2.2 Author’s contribution

The author performed the Molecular Dynamics calculations and analyzed the data in publicationsI,III andIV . In publicationII, the author carried out the Monte Carlo- based and Particle-in-Cell simulations and analyzed the corresponding data. The author wrote the draft of publicationsII andIV and contributed largely to the composition of the draft of publicationI . The author also contributed writing the Molecular Dynamics methods section in publicationIII. In publicationV, the author computed a small part of the simulations and took part in the ﬁnal corrections and submission of the manuscript.

Publications not included in the thesis

Unravelling the secrets of the resistance of GaN to strongly ionising radiation M.C. Sequeira, J.G. Mattei, H. Vázquez, F. Djurabekova, K. Nordlund, I. Monnet, P.

Mota-Santiago, P. Kluth, C. Grygiel, S. Zhang, E. AlvesCommunications physics. 2021 Mar 12;4(1):1-8.

Molecular dynamics simulation of the eﬀects of swift heavy ion irradiation on multilayer graphene and diamond-like carbon

J. Liu, H.V. Muíños, K. Nordlund and F. Djurabekova,2020. Applied Surface Science, 527, p.146495.

(12)

Highly active single-layer MoS2 catalysts synthesized by swift heavy ion irradiation

L. Madauß, I. Zegkinoglou, H.V. Muíños, Y.W. Choi, S. Kunze, M. Q. Zhao, C.H. Naylor, P. Ernst, E. Pollmann, O. Ochedowski and H. Lebius,2018. Nanoscale, 10(48), pp.22908- 22916.

Vaporlike phase of amorphous SiO2 is not a prerequisite for the core/shell ion tracks or ion shaping

H. Amekura, P. Kluth, P. Mota-Santiago, I. Sahlberg, V. Jantunen, A.A. Leino, H.

Vázquez, K. Nordlund, F. Djurabekova, N. Okubo, N. IshikawaPhysical Review Materi- als. 2018 Sep 4;2(9):096001.

Molecular Dynamics Simulations of Heavy Ion Induced Defects in SiC Schot- tky Diodes

A. Javanainen, H.V. Muíños, K. Nordlund, F. Djurabekova, K.F. Galloway, M. Tur- owski, R.D. SchrimpfIEEE Transactions on Device and Materials Reliability. 2018 May 31;18(3):481-3.

(13)

3 ION IRRADIATION OF MATERIALS

3.1 Types of ion irradiation

Ions travelling through a material interact with it mainly via two mechanisms: through repulsive forces between the nucleus of the projectile and the target atoms (nuclear interaction) or by exciting electrons of the target material and projectile (electronic interaction). The energy deposited in the material by the ion per unit length is often called stopping power. Nuclear and electronic stopping power correspond to the contributions to the stopping power of the ion associated with nuclear and electronic interactions, respectively. The relative importance of each of these contributions vary greatly with the energy of the projectile. Slow ions (less than 100 keV/amu) interact mainly via nuclear collisions and barely excite the electrons in the material; in this regime the nuclear stopping power is predominant. Projectiles in the intermediate energy regime can interact with the material via both electronic and nuclear mechanisms and both contributions to the stopping power are relevant. For fast ions, instead, electronic stopping power dominates as the ions deposit most of the energy via electron excitations.

In ﬁgure 1a we show the typical dependence of nuclear and electronic stopping power on projectile energy. In general, we see that heavier ions have much larger stopping powers than lighter ones. The ions lose electrons as they traverse the material at high speeds;

an increase in velocity leads consequently to an increase in charge state and electronic stopping power. The peak in electronic stopping power in ﬁgure 1a is called Bragg peak, and corresponds to the projectile velocity at which the ion becomes fully stripped of electrons. SHIs are characterized by their high velocity and span the region near to the Bragg peak all the way to higher velocities. In this regime, the ions interact mainly with the electrons in the material, leading to the formation of ion tracks (see Fig.1b).

Moreover, we shall mention also Highly Charged Ions (HCIs) [31–33]. These are usually heavy ions which can be stripped of electrons to very high charge states (> 10 e⁻) but have relatively low velocities. Despite their low velocities, their high charge state is able to excite many electrons and can produce signiﬁcant damage in the material. Unlike SHIs, which can travel several microns in the material, these ions only penetrate a few nanometers in the material and deposit all their energy in the surface.

3.2 Interaction of SHIs with the material

As a SHI travels through the material, it excites the electrons in the vicinity of the ion path to higher energy states. The electrons excited by the ion are called delta (δ) or primary electrons and often have large energies in the range of keV. Most of the excited electrons have velocities perpendicular to the ion trajectory, and only the most energetic ones have component in the ion direction [34]. Delta electrons travel balistically through

(14)

10⁴ 10³ 10² 10¹ 10⁰ 10¹ Projectile energy [MeV/u]

0 5 10 15 20 25

Stopping power [keV/nm]

He, electronic He, nuclear Au, electronic Au, nuclear

a b

Figure 1: (a) Nuclear and electronic stopping power as function of ion energy for He and Au projectiles in graphite according to SRIM database [29]. (b) Cross-section image with TEM of an ion track in quartz from Ref. [30]. Reprinted with permission from Elsevier.

the material and continue exciting more electron-hole pairs. The latter are referred as secondary electrons, and have considerably smaller energies than theδ-electrons. Despite their lower energy, secondary electrons do also scatter within the material, generating further excitation as they travel outwards from the track. Since each excited electron can excite more than one electron, the development of the electron excitations around the ion path is often referred as electron cascade.

The intensity of the ionization depends on the charge state [35–38] of the ion inside the material; generally, higher charge states induce larger ionization. The charge of the projectile changes throughout its trajectory. As the projectile enters the material, it starts to captures electrons from it; moreover, electrons from the projectile can be scattered oﬀ by the electrons from the bulk. The competition of both processes makes the charge saturate towards the equilibrium charge in that material. It takes approximately a few nanometers for an ion charge state to reach the equilibrium value, and after that, the charge state oscillates around it. The equilibrium charge of the ion in the material might not be a stable conﬁguration once the projectile exits. Some of the electrons captured by the ion are weakly bound once become detached outside the bulk; these electrons have the same velocity as the projectile and are often called convoy electrons [39].

Moreover, excited electrons reaching the surface can leave the material if they have enough energy to overcome the potential energy barrier. This energy barrier is composed of two contributions: work function, i.e. the energy necessary to excite an electron from the high- est occupied energy state into the vacuum, and electrostatic energy, the amount of energy that has to be spent to separate the charge of the electron from the positively charged layer. Emitted electrons are often referred as Secondary Electron Emission (SEE) [40–42].

Only electrons originated within a few nanometers from the surface are emitted into the

(15)

vacuum, therefore, the importance of this mechanism is negligible in bulk materials at perpendicular ion incidence. This eﬀect, however, can be critical under grazing incidence or in 2D materials, where the ion mainly interacts with the material within its surface.

Electrons in the material might be excited from the core shells or the valence bands. Holes in the core shell can be promoted into an upper shell or into the valence levels via two processes: Auger decay, where the decay is mediated by the emission of an outer electron, or radiative decays, where a photon instead is emitted. As the cascade evolves and the secondary electrons interact with each other, they become thermalized, i.e. their energy distribution follows the Thomas-Fermi distribution. This happens in the scale of tens of femtoseconds, however higher energy electrons can take longer times to equillibrate [34].

The development of the electron cascade, determines how the energy deposited by the ion is redistributed spatially around the ion path to be later transferred to the lattice.

3.3 Early electron dynamics after a SHI impact in 2D materials

Electron capture and SEE take place near the surface of the material. Two-dimensional materials by definition consist only of surface, and therefore these processes play a capital role during the evolution of the electron cascade induced by the ion. Two-dimensional materials comprise of a small number of atomic layers, and their thickness ranges from subnanometer to few nanometers. This short distance might not be sufficient for impacting ions to reach the equilibrium charge. Despite this, during the passage of the ion, the projectile can capture numerous electrons and deposit large amounts of energy in form of the potential energy of the holes created in the layer. The intensity of the electron capture process can contribute heftily to the deposition of energy; moreover, the charge of the projectile traversing the layer differs from the equilibrium charge in bulk, which can lead to different stopping power values than in bulk.

In 2D materials, the relative position of the ion impact point with respect to the lattice atoms can further affect the stopping power. In bulk materials, unless the ion trajectory is perfectly aligned with the positioning of the atoms in the crystal, the distance from the ion to the atoms on average is the same independently of where the ion enters the material. In 2D materials, there is only one layer of atoms, and the position where the ion crosses the layer can affect strongly the amount of energy the ion deposits in the material[43]. In figure 2a we showcase this effect for the impact of a H 25 keV ion in graphene.

The emission of electrons excited by the SHI also gains importance in 2D materials. As previously mentioned, in bulk, emission of excited electrons may only take place a few nanometers from the surface; in 2D materials, any electron is eligible to leave the surface if it has enough energy (schematics of electron emission shown in Fig. 2b). On the other hand, thin layers do not screen eﬃciently the charge of the holes left in the 2D material.

(16)

A B C F O impact point 0.40

0.45 0.50 0.55 0.60 0.65 0.70

number of electrons

e emission e capture

energy deposited [eV]

SRIM stopping [eV]

25 30 35 40 45

energy deposition [eV]

A COBF

a b

Figure 2: Figures from publication II. (a) Dependence of the deposited energy and number of emitted and captured electrons on the impact point of the ion in the graphene layer.

Results obtained with TDDFT for the ion H⁺80 keV. The inset shows the chosen impact points in the graphene layer. (b) Schematics of electron dynamics in the graphene layer during an ion impact.

The accumulation of charge exerts an attractive force on the electrons and may partially hinder the electron emission by increasing the barrier that electrons have to overcame to be emitted. These diﬀerences alter noticeably the dynamics of how the excitation induced by the ion evolves in the layer compared to bulk materials.

3.4 Energy transfer to the lattice

Electron cascades develop for few tens of femtoseconds, after which the electrons are thermalized. Electron-hole pairs around the track can recombine via several processes such as optical or Auger recombination, via excitons or by transferring energy to the lattice in form of collective vibrations. The strong excitation of the electron subsystem around the ion track leads to deformations in the lattice. However, the mechanism how the energy is transferred to the atoms and the material is damaged after the ion impact is not yet fully understood. There are several proposed mechanisms to explain the damage observed in SHI tracks: non-thermal melting,Coulomb explosion andthermal melting.

According toNon-thermal melting, the strong excitation of the electron subsystem alters the potential energy surface of the atoms, prompting a reorganization of the lattice while transferring large amounts of kinetic energy to the atoms. Nevertheless, for this reorganization to happen, the electron excitation has to be very high and it must be sustained enough time such that the atomic reorganizations takes place. This mechanism has been observed in some materials under intense electron excitation caused by laser pulses [44–

49]. Nevertheless, the excitations created by the ions are more localized and shorter-lived than in laser irradiation; this suggests that non-thermal eﬀects do not play a big role in SHI impacts [50].

(17)

TheCoulomb explosionmechanism proposes instead, that the lattice deformations appear as a consequence of the accumulation of positive charge in the track core [51–53]. The strong ionization produced by the ion would cause the electrons to move far from the ion track whereas the holes would remain positively charged in the center of the track.

The positively charged atoms would repel each other via Coulombic forces, transferring the potential energy into kinetic and melting the region around the ion path. For this mechanism to take place, the electrons leaving the core can not be replenished fast enough by the electronic subsystem in the material. This scenario is plausible in nanoscale isolated materials such as nanoclusters or molecules [54, 55]; however, in bulk is highly unlikely due to the abundance of electrons in the region around the track available to neutralize the charge in the core.

At last,thermal meltingsuggests that excited electrons deexcite with the help of atomic motions while promoting collective vibrations of atoms in the lattice, i.e. phonons. As the electrons deexcite, large number of phonons are generated, which translates into large kinetic energies and thermal motion—in other words, chaotic—of the atoms. The high temperatures in the track region prompt the melting of the lattice. This mechanism has been this far the most successful explaining most of the experimental observations in SHI tracks, and the scientiﬁc community has consensus that thermal eﬀects play a large role during SHI impacts [34, 56, 57].

(18)

4 ION-INDUCED DEFECTS

In the previous chapter we discussed the key processes taking place during the impact of a SHI and how the energy is deposited by the ion and later transferred into atomic motions. In this chapter we shortly review the defects that can form in 3D materials following the electron cascade generated by SHI impacts. We then proceed to the speciﬁcs of damage production in amorphous materials and graphene, the materials investigated in this dissertation.

4.1 Damage in bulk materials induced by SHI irradiation

SHIs can induce damage in a large variety of materials. Insulators have been shown to be particularly susceptible to SHI irradiation [34]; nevertheless, SHI tracks have also been observed in semiconductors [58–61] and metals [62–65]. Moreover, amorphous material are known to require smaller stopping powers to create tracks, and the size of the defects produced are larger than in their crystalline counterparts [34]. The composition of the ion track also varies depending on the material. Materials with stable amorphous phase (amorphizable) tend to show an amorphous core, whereas in non-amorphizable materials, the core is composed instead of a distinct crystal phase or contains large amount of defects [66, 67].

During the last decade, the development of experimental techniques has opened the possibility of exploring the fine structure within the ion track. Previous experiments have shown that in some crystals, the track exhibits a core-shell structure constituted of an amorphous core and a crystalline shell with defects [68, 69]. In some amorphous materials, a different core-shell structure was observed, composed of both amorphous core and shell regions but of different density [70, 71]. SHI impacts were also shown to be capable of annealing preexisting defects in SiC [72]. On the other hand, the opposite trend was observed in SrTiO3, where tracks are only observed if the material is pre-damaged with nuclear irradiation [73]. MD has proven to be an excellent tool to study the atomic dynamics during SHI impacts. This technique has been used to explain successfully the reaction of materials to SHI irradiation: recovery of defects [72], the elongation of metallic nanoclusters in amorphous matrices [74, 75] or the inner structure of the ion tracks [71, 76–78].

4.2 Ion tracks in amorphous materials

Amorphous materials exhibit clear short range order, i.e. the bond distances and bond angles of the atoms with their nearest neighbours show a strict ordering; however at further distances the arrangement of the atoms seems to form a random network (structure of a-SiO₂ shown in Fig. 3a). As a consequence, there is a vast number of atomic con- ﬁgurations with very similar energies in which the system can get trapped. The random

(19)

b c a

Figure 3: (a) Atomic bonds composing a-SiO₂ amorphous network. Oxygen and silicon atoms are represented as red light brown spheres respectively. (b) SAXS spectra of amorphous materials after SHI irradiation from publication III. (c) Schematics of core shell density structure (ﬁgure b) used to ﬁt the SAXS spectra from publication III.

nature of these materials does not allow the identification of changes in the structure of the network with imaging techniques as in the case of crystal lattices; nevertheless, network configurations can have different densities and this quantity can be easily measured.

The behaviour of amorphous materials under ion irradiation has been studied extensively in the last decades. Klaümunzer et al. showed that SHI irradiation induces dimensional changes in the material [79–82], and that the changes in dimensions are proportional to the ion fluence, i.e. the number of impacting ions per unit of area. This phenomenon is often referred asion hammering, and it has been attributed to the viscous relaxation of anisotropic SHI-induced stresses [83, 84]. In some materials such as a-SiO₂,ion hammer- ingis preceded by a compactification of the material at fluences lower than the incubation fluence [85]. Similar compactification has been observed under irradiation with ions in the nuclear regime [86, 87].

With the development of experimental techniques, it became possible to measure with high precision the density inside the tracks of amorphous materials. Small-angle X-ray scattering (SAXS) measurements showed that a-SiO₂ and a-Si exhibit a characteristic core-shell fine structure composed of an underdense core and an overdense shell [70, 71, 76] (see figures 3b and 3c). High-angle annular dark-field scanning transmission electron microscopy (HAADF-STEM) experiments also showed similar structures in a- Si₃N₄ thin layers [88]. The underdense core in the track has been associated with a disordered amorphous structure due to the fast quench of the track; this hypothesis was later supported by the observation of color centres inside the track cores of a-SiO₂ [89].

The overdense shell has been often attributed to the viscous ﬂow of defects initially proposed by Borodin et al. [90]. This model was introduced by the authors to explain

(20)

the mass redistribution inside the ion track during a SHI impact, however, it is not clear on what consists theﬂow of defects from an atomistic standpoint and how does it occur.

The core-shell structure has also been associated with the melt and boil radius in the track [88]. Moreover, Ref. [85] pointed at the density anomaly of a-SiO₂ in the molten phase to explain the densiﬁcation of a-SiO₂ samples at high SHI irradiation ﬂuences. In publicationIII, we study with atomistic simulations the inner composition of the core- shell structure of the track in SiO₂ and Si₃N₄. In publication IV we employ the same method to explain the mechanism how the core-shell forms in both those materials and why a-Si shows the opposite core-shell density structure.

Other amorphous materials show different features inside the track: a joint experimental and theoretical effort showed that SHI produce bow-tie voids in a-Ge, metallic glasses exhibit an underdense homogeneous core [91]. In some amorphous materials, the structural modifications in the ion tracks can lead to drastic changes in the macroscopic properties.

For example, several studies reported that sp³-rich amorphous carbon turns conductive after SHI irradiation [92–94]. In publicationV we employ MD simulations to examine the changes in the structure of a-C within the track, that are responsible for the property changes observed after irradiation.

4.3 Modiﬁcation of graphene with ion irradiation

Graphene is a material composed of a single layer of sp² hybridized carbon atoms. The atoms in the material are arranged in a honeycomb lattice forming a two-dimensional crystal (see Fig. 4a). Moreover, graphene is a semimetal of zero band gap, and exhibits a double cone band structure at theKpoint. Graphene has excellent tensile properties, and is capable of conducting electricity and heat extremely well; these unique characteristics, make it suitable for myriad of applications. Since it was ﬁrst discovered in 2004 [95], many studies have investigated the response of graphene to ion irradiation.

4.3.1 Low energy ion irradiation

Irradiation of graphene with low energy ions was shown to lead to the appearance of the D peak in Raman Spectroscopy measurements [96]. Raman Spectroscopy consists on shining an infra-red laser at the sample and measuring the shifts in frequency due to the inelastic scatterings of the photons with the lattice vibrations. The D peak shift in graphene and graphite is associated with symmetry-forbidden scatterings; defects in the lattice break this symmetry an prompt the appearance of the D peak in the Raman spectra [97]. Lucchese et al. reported an increase of the D peak with ion ﬂuence, and succeeded to model its behaviour by assuming a Raman active region around the defects;

the area of this region is defined by the characteristic relaxation length of the Raman scattering [96]. Despite the success of this approach, Eckmann et al. showed the difficulty of identifying specific defect types from the Raman spectra [98].

(21)

Ref. [99] studied with MD the production of vacancies and more complex defects in graphene by irradiation with keV ions. That work showed that faster ions can trigger in- plane nuclear cascades and produce larger defects. Focused Ion Beams (FIB) concentrate high ion ﬂuences in a narrow region. This technique was used for precise patterning of graphene [100]; nevertheless, it was shown to also produce undesirable outcomes in the layer such as amorphization [101].

4.3.2 Highly charged ion irradiation

HCI irradiation has recently been investigated in 2D materials. Studies in graphene showed that HCIs cause an increase of the Raman D peak [102]; however, up to date it has not been possible to link the Raman signal with the creation of pores, which suggests that the interaction with the defected substrate as the source of the D peak [103].

Theoretical work in Ref. [104] showed that graphene can provide tens of electrons to the HCI within few femtoseconds due its extreme mobility. Later research also pointed out that during the passage of a HCI up to 70 electrons can be emitted via Interatomic Coulombic Decay (ICD) [105]. All these evidences suggest that the D peak signal induced by HCIs appears due to sp bonding changes in the layer and that the superior electron transport in graphene prevent the formation of large defects. In other 2D materials, however, HCIs were shown to produce larger defects: in h-BN was observed increasing atomic sputtering with higher stopping power ions [106], whereas in MoS₂ and carbon nanomembranes HCIs were shown to produce nanometer size pores [107, 108].

4.3.3 Swift heavy ion irradiation

The controlled creation of nanopores in graphene has large potential for nanofiltering applications [103]; SHIs are deemed as a better alternative than low energy and HCI irradiation for the creation of such defects. Up to date there have been only few studies of SHI irradiation on 2D materials, most of them in graphene. The first studies of SHI irradiation were performed at grazing incidence irradiation angles. In this setup, the ions travel long distances parallel to the layer as they deposit energy via electron excitations in a line-like morphology. Experiments showed that a single ion can produce tears of several hundreds of nanometers [109, 110] (experimental images shown in Figs. 4b-d), and that the layers in the affected region can fold [111]. The length of the defect and the folding morphology were shown to depend strongly on the substrate [110]; moreover, it was observed for graphene that increasing number of layers reduced the probability of folding.

The defects produced under grazing incidence angles are too large for nanoﬁltering applications; for this purpose normal incidence is more promising. Measuring defects produced

(22)

b c

d a

Figure 4: (a) Schematic illustration of graphene lattice. (b-d) Experimental images from Ref. [110] of defects produced in supported (b, d) and suspended (c) graphene after irradiation with SHIs at grazing angle incidence. Figures relabeled and reprinted with permission from IOP publishing.

under perpendicular ion irradiation has proven, however, more challenging due to two fac- tors: lack of availability of in-situ STEM measurements in the irradiation facilities and impossibility to distinguish between ion-induced defects and preexistent ones from man- ufacturing or manipulation procedures [103]. These limitations are further aggravated in the case of graphene by its capability to heal itself [112], and its reactivity when exposed to air which leads to adsorption of hydrocarbons and impedes its imaging. Up to date, most studies of SHI perpendicular irradiation have employed the same Raman analysis as Lucheese et al. [96] to quantify the amount of defects induced in the graphene layer.

This approach, however, has two limitations: it provides no information regarding the types of defects, and the model depends on many ﬁtting parameters such as the eﬃciency for defect creation or the Raman relaxation length.

The high mobility of electrons in graphene [113] raises the question whether SHIs can produce large defects at all, or if the electrons in the layer are capable of dissipating the energy fast enough. Recent experiments investigated the irradiation of supported graphene followed by etching of the tracks in the substrate for nanofiltering application [114]. That work shows an increase of the permeability of device with fluence. This is only possible if pores in graphene coincide with the etched pores in PPMA—keep in mind that the etchant does not affect the carbon atoms—and supports the notion that pore-like defects are created in the layer.

Computational studies are capable of bypassing the experimental limitations allowing to observe the process of defect formation in graphene during a SHI impact as it unfolds.

Zhao et al. [115] made a ﬁrst attempt to model SHI impacts in graphene; the authors employed atomistic simulations, and modelled the electron subsystem of graphene with parameters ﬁtted to graphite. They predicted as stopping power threshold for defect formation 6.5 keV/nm for supported graphene and 8 keV/nm for suspended, whereas the

(23)

value for graphite was between 6.7-8 keV/nm [116, 117]. In publicationIII we employ Molecular Dynamics combined with the Two-temperature model to simulate SHI impacts in graphene, and correlate the defect size with the signal in Raman experiments to determine if SHIs form pore-like defects in the layer. In publicationII, we investigate the electron cascades generated by the SHI impacts in graphene, and how they may aﬀect the formation of defects.

(24)

5 METHODS

Simulations oﬀer the possibility to investigate a physical process at time and length scales not accessible to experimental techniques, and to test our current theoretical understanding of those process. In this chapter, we give an overview of the methods employed in the publications of the dissertation.

In publicationI,III,IVandVwe employed Molecular Dynamics to investigate the transformations at atomic resolution within the ion track during a SHI impact. This technique was combined with the Two-temperature model to describe how the energy deposited by the SHI to the electrons in the material is spatially redistributed and translated in into atomic vibrations. We applied these techniques in publicationI to determine whether SHIs can produce nanopores in graphene, in publicationIII, andIVto probe the inner structure of the core-shell pattern inside the ion track, and to reveal the mechanism how it forms, and in publicationVto investigate the transformations triggered by the SHI impact that lead to the change of macroscopic properties of sp³-rich carbon in experiments.

In section 3.3, we reviewed the electron cascades triggered by SHI ions, and the diﬀer- ent processes taking place in two-dimensional materials, such as electron emission and electron capture. In publication II we employed a Monte Carlo-based technique to investigate the electron cascades in graphene, and in particular, how the electron emission contributes to the dissipation of part of the energy deposited by the SHI. We also applied Particle-in-cell simulations to describe the eﬀect of the electrostatic forces during the process of electron emission. In the same publication, Time-dependent density functional theory was used to study the electron emission in graphene.

5.1 Atomistic Simulations

5.1.1 Molecular Dynamics

Molecular Dynamics was developed in the 1950s by Alder and Wainwright [118, 119]. This technique simulates the motion of the individual atoms by solving Newton’s equations of motion

fi = mir¨i, (1)

fi = −∂Ui

∂ri, (2)

where fi and ri are the force and position of the ith atom in the system, and U the potential energy. This equation can be derived from the Schrodinger equation—the equation that dictates the dynamics of particles at quantum level—applying the Born- Oppenheimer [120] approximation. This approximation assumes that, due to the diﬀer- ence in mass between electrons and atoms, the electrons dynamics are much faster than

(25)

the movement of the atomic nuclei, and therefore the electrons in the system relax to the ground state, i.e. the lowest energy level, in a faster timescale than that of the atomic motion. This allows the creation of a potential energy surface as a function of the nucleus positionUi({r}). The potential energy surface can be constructed with quantum mechanicalab initiomethods such as Density Functional Theory [121], however, the computational complexity of the simulations increases steeply with the number of atoms in the system. The use of interatomic potentials—mathematical functions ﬁtted to describe the interaction between atoms—allows this limitation to be overcome and enables the simulation of millions of atoms simultaneously.

In MD, we have a finite number of particles in the simulation, and therefore, at a given time we only have a reduced region of the phase space, i.e. all the possible velocity and position configurations for the atoms, of the thermodynamical ensemble. The connection between the atomistic simulations and the thermodynamical properties is guaranteed by the ergodic hypothesis, which states that the average of a given thermodynamical property over time is equivalent to its average over the ensemble. This allows to identification of the kinetic energy and stresses in the simulations with their equivalent thermodynamical properties of temperature and pressure.

The stress per atom is deﬁned with the virial stressesτik as follows:

τ_ij⁽^l⁾= 1 Ω⁽^l⁾

−m⁽^l⁾(u⁽_i^l⁾−u^avg₍_i₎)(u⁽_j^l⁾−u^avg_j )−1 2

N_n

m

(r_i⁽^l⁾−r⁽_i^m⁾)f_j⁽^lm⁾

, (3)

whereldenotes the atom index andΩ,m⁽^l⁾,r⁽^l⁾andu⁽^l⁾its corresponding atomic volume, mass, position and velocity respectively; the average velocity of the system is denoted as uâvg. Moreover, f⁽^lm⁾ represents the force acting between atoms l and m; this term is computed for all atoms interacting with atoml(in our case, all the atoms within a cutoff distance around the particle). Moreover, the temperature in a given volume is defined as

T =

N_V

l=0

m⁽^l⁾(u⁽^l⁾)²

3NVkB , (4)

wherekB is the Boltzmann constant andNV is the number of particles contained in a given volumeV.

The simulations described in this dissertation were performed with PARCAS [122] and LAMMPS [123] codes. The former code employs the GEAR5 [124] algorithm to solve equation 2, whereas the latter one uses the Verlet algorithm [125]. The choice of code was inﬂuenced mainly by the availability of implementation for the desired potentials, as most of them are implemented in either one or the other code. The atomic structures in our simulations were optimized applying the Berendsen [126] thermostat and barostat in the simulations performed with the PARCAS code [122]. LAMMPS [123] includes the implementation of the Nose-Hoover [127] thermostat and barostat, known to give a better

(26)

description of the thermodynamical ensemble. Therefore, we decided to use instead the Nose-Hoover thermostat and barostat in all simulations performed with the LAMMPS code.

5.1.2 Two-temperature model (TTM)

The two-temperature model was ﬁrst introduced by Lifshift for metals [128], and describes how the energy initially deposited by the ion in form of electron excitation is later transferred to the lattice atoms. This model pressumes that defects in the material are formed via theinelastic thermal spike (i-TS) mechanism and that the energy is transferred to the atoms in a purely thermal manner. This model assumes that the electrons and the lattice atoms constitute two separate thermodynamical subsystems, in other words, that the electrons thermalize with themselves much faster than the atoms. These assumptions allow the modelling of the energy transport with two separate diﬀusion equations—one for electrons and another for the lattice—and an exchange term between them

Ce∂Te

∂t =∇(ke∇Te)−G(Te−Tl) +A(r) (5) Cl∂Tl

∂t =∇(kl∇Tl) +G(Te−Tl) (6) whereci andki are the heat capacity and thermal conductivity andi= e, l denote the electronic and lattice subsystem, respectively. The termGis called electron-phonon coupling, and determines how fast the energy stored in form of excited electrons is transferred into atomic motions. The term A(r) is spatial distribution of the energy source in the electron subsystem. In this case, the source term corresponds to the energy distribution due to the electron cascade produced by the ion once the electrons have thermalized. This distribution can be calculated for example with Monte Carlo codes as in Ref. [129–131].

In the literature often an analytical distribution is employed. For example, the Waligorski distribution [131] was originally ﬁtted to electron cascade simulations in water, however it has been applied in Ref. [77, 132] to simulate impacts in various materials.

The TTM was originally developed for metals. Its application to materials with band gap is rather controversial as, unlike in metals, in these materials the density of the carriers has spatial and temporal dependence and they can diﬀuse. Nevertheless, up to date, TTM has been applied to model the formation of SHI-induced defects in insulators with reasonable success [78, 132]. The main inconvenience of this model is the great diﬃculty of estimating the electronic parameters in equation 5. Due to this limitation, oftentimes one or several of the electronic parameters are used as free parameters.

The heat capacity of electrons can be estimated accurately if the Density of States (DOS) g(r) is known as

Ce=

g(E)∂f(E, μ(Te), Te)

∂Te EdE. (7)

A rougher estimation can be used according to the classical limitCe= ³₂nkB, wherenis the carrier electron density in the material andkBthe Boltzmann constant. Alternatively,

(27)

in insulators Dufour et al. assume a linear heat capacity up to temperature corresponding to the optical bandgap of the material, and which the heat capacity takes the value of the free electron gas [133].

The electronic thermal conductivity and the electron-phonon coupling can be estimated as

ke = CeDe, (8)

G = Ce

τ , (9)

where De is the electronic diﬀusivity and τ is the relaxation time—the characteristic time of temperature decay in the electron subsystem—, in case these parameters are known. Both electron thermal diﬀusivity and the relaxation time can be estimated from sub-femtosecond pump probe experiments [134–136]. Whenever there is no experimental data available, in the literature the electron-phonon coupling is often estimated with the expression

G= CeDe

λ² (10)

, where λ is the so-called electron-phonon mean free path. Toulemonde fitted the G parameter to experimental data in different amorphizable insulators and proposed an empirical relation linking the mean free path and the band gap [137]. For metals, it is also possible to calculate the electron-phonon from first principles [138].

5.1.3 Combining Inelastic thermal spike with MD

The energy deposited by the ion in the electron subsystem is eventually transferred to the atoms, triggering material transformations around the ion path. Molecular Dynamics allows the study of material changes due to the fast and localized heating at the atomic level. The spatial distribution and the rate at which the energy is deposited to the atoms determine the structural defects. A simplistic approach to model the transfer, is to deposit homogeneously within a cylindrical volume the energy corresponding to the stopping power of the ion in the material [89]. In this case, the radius of energy deposition is a free parameter which can be adjusted to match the experimental track radius.

A more reﬁned approach, is to simulate the energy redistribution in the electronic and lattice systems with TTM and deposit the energy to the lattice in the MD cell. The deposition of energy is performed by adding the energy to the atoms as velocity in random directions. The energy can be deposited either instantaneously or in a time-dependent fashion. In the instantaneous scheme, often referred as two-temperature MD (TTMD), the time at which the energy is deposited is chosen such that the lattice temperature in the centre of the track reaches its maximum value. In the time dependent scheme, the energy is deposited instead at each MD time step; the spatial and temporal behavior of the deposition rate is extracted from the solution of the TTM diﬀerential equations. The

(28)

instantaneous TTMD model was employed in publicationIII,IVandVto simulate the SHI impacts in amorphous materials.

5.1.4 Two-temperature Molecular Dynamics model (2TMD)

In the TTM model, the dynamics of both the electron and lattice subsystems are described with two diﬀerential equations. In MD, the dynamics of the lattice are simulated explicitly by simulating the motion of each atom in the system. Therefore, the two-temperature model can be redeﬁned such that the electron dynamics are governed by equation 5, but the lattice dynamics are simulated explicitly with the MD algorithm. The 2TMD approach was introduced by Ivanov et al. [139] to model short-pulse laser irradiation in metals.

d²ri

dt² = F+σmiri (11)

σ = GVN(Te−Tl)

imi(vi)² , (12)

(13) wheremi, ri,vi,Fi are the mass, position, velocity and acting force of the atomi. The diﬀerential equation 5 is solved on a simulation grid divided into voxels;N indicates the number of atoms contained in a given voxel andVN its volume.

This model is used in publicationIto simulate the heat redistribution during the impact of a SHI in graphene. The grid employed for the solver had dimensions 21×21×1, as shown in Fig. 5a. The width of the simulation voxel was of 1 nm to ensure enough atom counts in each cell (∼50 atoms per cell in our case) and to avoid ﬂuctuations in the solution.

5.1.5 Simulation of SHI impacts

Before simulating the ion impacts, the atomic structures are relaxed in MD at the desired pressure and temperature (normally 0 GPa and 300 K) with a barostat and a thermostat.

This ensures that the structure is optimized before simulating the ion impact. The relaxation is normally performed in small simulation cells, and later replicated to create larger ones. The simulations are performed with periodic boundaries in all directions, except for 2D materials, where the out-of-plane direction is non-periodic to avoid self- interaction of the atomic layers. The cell dimensions are chosen large enough such that the waves induced in the material by the expansion of the molten track do not interact through the periodic boundaries.

The region near the ion impact point, at the center of the cell, is simulated with the NVE ensemble to not alter the dynamics of the atoms involved. We apply at the edges of the cell (schematics of typical setup shown in Fig. 5a-b) border cooling to 300K with

(29)

a b

Figure 5: Figures of irradiation setup in MD simulations. (a) Simulation setup of SHI irradiation of graphene with the 2TMD model. The grid used for solving the diﬀerential equation is shown as a projection on the bottom of the cell. (b) Simulation setup of SHI impacts on a-SiO₂. In both setups the ion impact takes place in the center of the cell, and the border regions (shown in a diﬀerent color) are cooled with a thermostat.

a thermostat to emulate the cooling eﬀect from the inﬁnite bulk material surrounding the track region. The simulations are performed for 100 or even 200 ps depending on the material, until the track region is cold and the dynamics of the atoms do not evolve further.

5.2 Non-equilibrium electron dynamics during SHI impacts

5.2.1 Monte Carlo-based model for electron cascades

Monte Carlo asymptotic trajectory event-by-event method approximates the propagation of individual particles as a sequence of scatterings. This approach was used already in the late 60s and in the 70s to study the ion-induced electron cascades in detectors [140, 141] as well as in biological tissues [142, 143]; later this techniques was applied also to other materials to investigate the excitation induced by SHIs [129, 144]. The Time- Resolved Electron Kinetics in SHI irradiated Solids (TREKIS) code [130] is based on the aforementioned model. In publicationII, we employed TREKIS to simulate the dynamics of the electron cascades initiated by ion impacts. Unlike other models, this code takes into consideration the collective modes of excitation in the solid by using the Complex Dielectric Function (CDF) formalism [145] to compute the scattering cross sections in the material. In the rest of this dissertation we refer to this model as MC-CDF.

The model describes the generation of primary electrons (δ-electrons) by the projectile as well as their propagation and subsequent scatterings which create additional electron-hole pairs (see schematics in Fig.6). The ion and the carriers scatter randomly in the solid with a mean free pathλ= 1/(nσ), wherenis the density of scattering centers andσ the scattering cross section. The projectiles—ion, electrons and holes—can interact either inelastically, i.e. exciting new electron-hole pairs in the material, or elastically, scattering

(30)

Figure 6: Schematics of generation and transport of carriers in TREKIS simulations. The carriers travel in the material until they reach the emission surface with enough energy and are emitted or until their kinetic energy drops below the energy cutoﬀ, after which they stop being simulated.

from the screened atomic nuclei of the target atoms. The inelastic cross section of the projectiles are obtained according to the expression

d²σ

d(ω)d(q) = 2(qeff(ν, Z)e)² nπ²ν²

1 qIm

−1 (ω, q)

, (14)

where qeff, Z and ν are the projectile eﬀective charge, atomic number and velocity, ω andqthe energy and momentum transfer in the scattering and Im ₍⁻¹_ω,q₎

the loss function in the material. The loss function can be decomposed as sum of oscillators [145]

as follows

Im

−1 ε(ω, q= 0)

=

N_osc

i=1

Aiγiω

(²ω²−E₀²_i)²+ (γiω)², (15)

whereNosc is the number of oscillators used to approximate the loss function, and Ai, E₀iandγi are the parameters of thei-th oscillator. These oscillator parameters can be obtained by ﬁtting them to experimental data. In this model, the holes interact only via elastic scatterings. The elastic scattering of carriers with atoms is described with the Mott’s cross section. Moreover, electrons can leave the material if they reach the surface with the kinetic energy suﬃcient to overcome the emission barrier [146] in the material. The model also includes deexcitation channels for the core-shell electrons, either via Auger or radiative decays.

(31)

5.2.2 Particle-In-Cell for simulation of charge dynamics in vacuum

Particle-In-Cell (PIC) is a simulation technique that solves self-consistently the motion of large number of charged particles under electromagnetic fields. It was first introduced in 1955 [147], and it became rapidly popular in the 60s in the field of plasma physics. In this method the particles are tracked in the continuous phase space, whereas the electromagnetic fields are solved self-consistently on a fixed grid. The algorithm is composed of the following steps:

1. Update particle trajectories

2. Interpolate density and currents on grid points 3. Solve self-consistently the ﬁeld on the grid 4. Interpolate the ﬁeld at the particle positions 5. Repeat the sequence.

The propagation of the individual particles can become computationally heavy when the number of particles becomes excessively large. This limitation can be avoided by grouping several particles forming what are called superparticles. In a similar manner, superparticles composed of a fractional number of physical particles can be used to solve electrostatic problems with low electron densities.

In publication II, we simulated with PIC the dynamics of the electrons emitted in graphene. The simulations were performed with the FEMOCS framework [148, 149], which uses Finite Element Methods (FEM) to solve the Poisson equation and obtain the electromagnetic ﬁeld for a given geometry.

5.2.3 Time-dependent Density Functional Theory

Time-dependent Density Functional Theory (TDDFT) is based on the Runge-Gross theorem [150, 151], the time-dependent equivalent of the Hohenberg-Kohn theorem [151]

in Density Functional Theory (DFT). This theorem states that in a quantum system, it is possible to map the time dependent density to the time dependent potential and vice versa. This relation becomes a unique mapping if the quantum system starts from the ground state [151]. Therefore, the manybody interacting quantum system can be transformed into a system of non-interacting particles under a time-dependent external potential vs(r, t). For a given system of interacting electrons with density n(r, t), the time-dependent external potentialvs(r, t)is postulated to be such that a system of non- interacting electrons following the Schrödinger equation present also the same density