Ab initio study on localization and finite size effects in the structural, electronic, and optical properties of hydrogenated amorphous silicon

Computational Materials Science 155 (2018) 159–168

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Ab initio study on localization and finite size effects in the structural, electronic, and optical properties of hydrogenated amorphous silicon Philippe Czajaa, Massimo Celinob, Simone Giusepponib, Michele Gussoc, Urs Aeberharda,

T

⁎

a

IEK-5 Photovoltaik, Forschungszentrum Jülich, 52425 Jülich, Germany ENEA, C. R. Casaccia, Via Anguillarese 301, 00123 Rome, Italy c ENEA, C. R. Brindisi, 72100 Brindisi, Italy b

A R T I C LE I N FO

A B S T R A C T

Keywords: Hydrogenated amorphous silicon Molecular dynamics Electronic structure Optical properties

We present a first-principles study of the structural, electronic, and optical properties of hydrogenated amorphous silicon (a-Si:H). To this end, atomic configurations of a-Si:H with 72 and 576 atoms respectively are generated using ab initio molecular dynamics, where the larger structures are defect free, closely matching the experimental situation and enabling the comparison of the electronic and optical properties with experimental results. Density functional theory calculations are applied to both configurations in order to obtain the electronic wave functions. These are analyzed and characterized with respect to their localization and their contribution to the density of states, and are used for calculating ab initio absorption spectra of a-Si:H. The results show that both the size and the defect structure of the configurations modify the electronic and optical properties and in particular the value of the band gap. This value could be improved by calculating quasi-particle (QP) corrections to the single-particle spectra using the G0W0 method. We find that the QP corrections can be described by a set of scissors shift parameters, which can also be used in calculations of larger structures. The analysis of individual contributions to the absorption by evaluating the optical matrix elements indicates that strong localization enhances the optical coupling, but has little effect on the average transition probability 〈|vcv |2 〉, for which we find a dependence E 2 + const on the photon energy E, irrespective of the nature of the initial or final state.

1. Introduction Hydrogenated amorphous silicon (a-Si:H) has been used as a cheap and efficient absorber material in silicon thin-film solar cells for more than 40 years [1], and has lately found another application in photovoltaics as a passivation layer in silicon-heterojunction cells. Understanding its microscopic structure in order to optimize its macroscopic properties for the application in photovoltaics has motivated several ab initio studies of a-Si:H throughout the years [2–8]. Two principle challenges are thereby met. First, a model atomic structure has to be generated that correctly reproduces certain experimental features of aSi:H, such as the defect density, the radial pair correlation function, or the vibrational properties. Second, the electronic structure has to be calculated on a level that allows for the extraction of physically meaningful macroscopic properties. From the viewpoint of photovoltaics, special interest lies on the description of the optical properties and on the identification and characterization of localized defect states, which have a crucial impact on the device performance due to their role as recombination centers in non-radiative recombination [9].

⁎

While the generation of defective a-Si:H configurations, i.e., configurations containing dangling bonds, is instructive for studying the origin and the nature of localized defect states, these configurations are not well suited for obtaining realistic macroscopic properties, due to an overestimation of the defect density. In fact, structures containing one defect need to have a size of at least 106 atoms to yield realistic defect densities [10], which is out of the range of current studies dealing with structure sizes of the order of 1000 atoms. The generation of defect-free configurations is therefore an important step towards a full ab initio description of a-Si:H. However, for a long time defect-free configurations of a-Si and a-Si:H could be generated only with model approaches such as the Wooten-Winer-Weaire algorithm [11], the Bethe-lattice approach [2] or the Reverse Monte-Carlo approach [12]. Only recently, large-scale (∼500 atoms) atomistic simulations of a-Si:H using a quench-from-a-melt approach [7,13] combining both classical and ab initio molecular dynamics (MD) have been reported to yield configurations of low defect density [8]. The same approach was used to generate low-defect and even defect-free a-Si:H configurations of 72 atoms within ab initio MD [6]. Following this approach, we present

Corresponding author. E-mail address: [email protected] (U. Aeberhard).

https://doi.org/10.1016/j.commatsci.2018.08.027 Received 27 April 2018; Received in revised form 27 July 2018; Accepted 15 August 2018 0927-0256/ © 2018 Elsevier B.V. All rights reserved.

Computational Materials Science 155 (2018) 159–168

P. Czaja et al.

approximation (RPA) [32]. Electron-hole interaction is disregarded as it is generally assumed to have no significant effect on the absorption spectra of amorphous semiconductors [10].

atomistic ab initio MD simulations on a-Si:H, resulting in defective and defect-free configurations of 72 and 576 atoms, respectively. The calculation and analysis of the electronic structure and the optical properties of a-Si:H on the Density-Functional-Theory (DFT) level has been the subject of a number of recent works [4–8]. The focus of interest in these works has been mainly on the origin of mid-gap states and band tails, and on the effect of hydrogen concentration and structural features on the mobility gap and the optical gap respectively. Very little attention has however been paid to the effect of computational artifacts on the electronic and optical properties. In particular, two effects should be taken into account when trying to reproduce the experimental properties of a-Si:H, and are therefore investigated in this work: the effect of the super-cell size and the effect of many-body interactions. A recent work stated that finite size effects do not play any role for structures larger than 72 atoms [6], which however disagrees with our findings. The incomplete description of many-body effects on the other hand is a well-known problem of standard DFT [14], and is the reason why the optical and mobility gaps are severely underestimated in previous studies using the local density approximation (LDA) or the generalized gradient approximation (GGA). Good values for the gaps have however been achieved recently with hybrid functionals [8]. In this work we try to incorporate many-body interactions systematically by explicitly calculating the quasi-particle corrections [15] to the Kohn–Sham energies. These corrections are often described by a heuristic approach – termed scissors shift (SS) [16] – where the electron energies are simply shifted to fit the experimental band gap. Since a distinct experimental value of the band gap of a-Si:H does however not exist, a set of shifting parameters can only be determined from a GW calculation. Here we present the results of such a calculation.

3. Generation of the atomic structure The first step towards the ab initio description of a-Si:H is the generation of physically meaningful atomic configurations, i.e., configurations that reproduce the experimental properties of a-Si:H. Our method of choice to achieve this is a simulated-annealing quench-froma-melt protocol [7,13,33]. With this method two configurations are generated, one consisting of 64 Si + 8 H atoms, which we will refer to as the small system, and one consisting of 512 Si + 64 H atoms, which we will refer to as the large system. Both structures are cubic with a size of a = 11.06 Å and a = 22.12 Å respectively, resulting in a density of 2.214 g/cm3 that matches the experimental value [7]. A hydrogen concentration of about 11% is chosen as this is the nominal concentration set in experimental materials optimized for photovoltaic performance [34]. The small structure is generated by randomly inserting H atoms in a crystalline Si cell, avoiding placing atoms too near to each other to minimize the energy of the system. A long (≈ 50ps ) BOMD simulation is then performed on the system at constant volume and temperature T = 2000 K, controlled by an Andersen thermostat [35]. This simulation assures that any remaining crystalline symmetry in the atomic configuration is destroyed. During the high temperature simulation, atoms cover a distance of about 2 nm, which ensures that the final configuration retains no memory of the initial geometry. Afterwards the temperature is lowered to T = 1200 K within 5 ps and held constant for another 5 ps, as an intermediate step prior to the final quench to T = 300 K. This quench takes another 5 ps, which results in a high quench rate of ∼ 1014 K/s. The quenched configuration is then used as a starting point for a final BOMD simulation at room temperature to fully thermalize the amorphous system. After 30 ps of simulation another 5 ps are used to compute average physical quantities and to characterize the system with respect to structural and electronic properties. The final configuration at the end of the simulation is shown in Fig. 1. The starting configuration for the large system is produced by replicating the small system in all directions. Since a high temperature annealing was already performed for the small system, a low temperature annealing in order to minimize the spurious defects at the internal interfaces is sufficient. For that purpose BOMD simulations are performed for 80 ps in time steps of 20 a.u. at constant volume. The temperature is controlled by a Nosé thermostat [36]. Within the first 60 ps the temperature is modified from 300 K to 600 K and back to 300 K in steps of 100 K. Afterwards, an additional simulation run at T = 300 K is performed for 20 ps. Fig. 2 shows the final configuration at the end of the simulations.

2. Technical details Born–Oppenheimer molecular dynamics (BOMD) simulations and electronic structure calculations are performed on the DFT [17,18] level, using a PBE-GGA exchange–correlation functional [19] and periodic boundary conditions (PBC), meant to mimic an infinitely extended system. For the BOMD simulations of the small structure (72 atoms) the PWscf (Plane-Wave Self-Consistent Field) code of the Quantum ESPRESSO suite [20,21] is used with ultrasoft pseudopotentials, whereas for the large structure (576 atoms) the Quickstep code of the CP2K suite [22] is used with norm-conserving Goedecker-Tetter-Hutter pseudopotentials [23–25]. All MD simulations are restricted to the Γ -point, which is justified by the size of the super cell. The electronic structure is calculated with the PWscf code of the Quantum ESPRESSO package using norm-conserving pseudopotentials. k -point summations are carried out on a 4 × 4 × 4 grid for the small system, and on a 2 × 2 × 2 grid for the large system. The plane-wave cut-off energy is set to 52 Ry. These parameters were chosen by checking the convergence of the total energy of the system. Quasi-particle corrections to the Kohn–Sham energies for the small configuration are obtained by performing single-shot G0W0 calculations [26] with the BerkeleyGW code [27] within the generalized plasmonpole (GPP) approximation [28–30] on a 2 × 2 × 2 grid using 3000 bands and a kinetic energy cut-off of 10 Ry. These values were chosen by checking the convergence of the LUMO–HOMO gap with respect to all three parameters simultaneously. The Kohn–Sham wave functions are retained as they are assumed to differ very little from the quasiparticle wave functions [30]. The GPP approximation is chosen because it has the advantage of requiring the dielectric tensor ∊ (ω) only in the static limit ω → 0 , while still taking into account the effects of dynamical screening. This ensures an accuracy similar to a full-frequency calculation for many semiconductors, including c-Si [31,30]. The BerkeleyGW code is also used for calculating the optical properties within linear-response theory using the random phase

4. Structural properties In order to characterize the system and estimate its quality in terms

Fig. 1. Small a-Si:H configuration in the simulation box. The super cell consists of 64 Si atoms (yellow) and 8 H atoms (blue). Periodic boundary conditions are used in all three directions. The configuration is obtained by quenching from 2000 K to 300 K at a rate of ∼ 1014 K/s and subsequent thermalization at 300 K. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 160

Computational Materials Science 155 (2018) 159–168

P. Czaja et al.

Fig. 2. Large a-Si:H configuration in the simulation box. The super cell consists of 512 Si atoms (yellow) and 64 H atoms (blue). Periodic boundary conditions are used in all three directions. The configuration is obtained by replicating the small configuration, and then annealing the system while changing the temperature stepwise from 300 K to 600 K and back to 300 K. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

of how it compares to real a-Si:H we study the radial pair correlation function gαβ (r ) , which is a measure for the probability of finding an atom of species α at a distance r from an atom of species β . It is defined as

gαβ (r ) =

V

Nα

Nβ

∑ ∑

4πr 2Nα Nβ I = 1 J = 1

δ (r −|RI −RJ |),

(1)

where V is the super-cell volume, Nα / β are the numbers of atoms of the respective atomic species in the super cell, and R are the atomic positions. The characterization starts after the system has reached thermal equilibrium. Fig. 3 shows g (r ) for both systems averaged over the last 5 ps of the simulation. Both figures reveal sharp peaks at 2.37 Å for Si-Si pairs, and at 1.51 Å for Si-H pairs. The only difference is in the smoothness of the curves, whereas the position and height of the main peaks remain unchanged, indicating an equal structural quality of both configurations. The comparison of gSi − Si (r ) for the large system with another recent work by Khomyakov et al. [7] as well as experimental data for a-Si:H [37] (Fig. 4) shows very good agreement at all distances, indicating that our configurations reproduce the structural properties of real a-Si:H very well. Analyzing the data, we observe that the ab initio results of Khomyakov et al. [7] are characterized by a lower first peak. This is likely due to the presence of some structural defects in their system, since structural defects are able to break the tetrahedral symmetry leading to a lower first peak in the radial pair correlation function, and such defects are unavoidable with their method to generate atomic configurations. To gain a deeper understanding of the structural properties, a coordination analysis of the Si atoms is performed and reported in Table 1. On the basis of the calculated radial pair correlation functions, a geometrical criterion is used to identify the nearest neighbors in the coordination analysis, applying a distance cutoff of 2.85 Å and 1.7 Å for Si-Si pairs and for Si-H pairs, respectively. These values are used to characterize both systems and correspond to the distances where the first minima are located in the radial pair correlation functions. Concerning the small system, it is observed that the average number of neighbors of a Si atom is 3.96. In detail, 4 Si atoms have threefold coordination (6.3%), 58 Si atoms have fourfold coordination (90.6%) and the remaining 2 Si atoms have fivefold coordination (3.1%). In contrast, the large system has an average coordination number equal to 3.99 because 99.0% of the Si atoms have fourfold coordination. Only 4 Si atoms (0.8%) have threefold coordination and 1 Si atom is fivefold coordinated (0.2%). In conclusion, the large systems is more ordered with a higher percentage of fourfold coordinated Si atoms. This trend is

Fig. 3. Radial pair correlation functions g (r ) of the small system (a) and the large system (b) for Si-Si pairs (black), Si-H pairs (green), and H-H pairs (orange), averaged over the last 5 ps of the MD simulation. Both gSi − Si and gSi −H are almost identical for both systems, showing sharp peaks at 2.37 Å and 1.51 Å , respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. Comparison between experimental and computed (large system) Si-Si radial pair correlation functions gSi − Si (r ) . The experimental data is from a neutron scattering paper (Ref. [37]). The work of Khomyakov et al. (Ref. [7]) is reported for comparison but data are computed on an amorphous system with some structural defects, which could have broken the overall tetrahedral symmetry and affected the height of the first peak.

161

Computational Materials Science 155 (2018) 159–168

P. Czaja et al.

1

Table 1 Coordination analysis for Si atoms in both a-Si:H systems. The number of Si atoms with a certain coordination is computed in two ways: the geometric protocol and the ELF approach. The number in brackets gives the percentage. The values show that in the large system the structural order is significantly improved as compared to the small system. Coord

1 2 3 4 5

Small system

0.8

ELF

0.6

ELF

Geom.

ELF

0 0 4 (6.3) 58 (91) 2 (3.1)

0 1 (1.6) 0 63 (98) 0

0 0 4 (0.8) 507 (99) 1 (0.2)

0 0 0 512 (100) 0

0 0

Small system ELF

Geom.

ELF

Si1H0 Si0H1

0 0

0 0

0 0

0 0

2

Si2H0 Si1H1 Si0H2

0 0 0

1 0 0

0 0 0

0 0 0

3

Si3H0 Si2H1 Si1H2

4 0 0

0 0 0

4 0 0

0 0 0

4

Si4H0 Si3H1 Si2H2

51 7 0

55 8 0

443 64 0

448 64 0

5

Si5H0 Si4H1 Si3H2

1 1 0

0 0 0

1 0 0

0 0 0

1

1

1.5 distance [Å]

2

2.5

5. Electronic properties 5.1. Bonding In order to relate the electronic properties to the atomic structure, and in particular identify structural defects, we use the electron localization function [38], which enables us to determine the coordination of each atom by identifying covalent and dangling bonds. A covalent bond is formed by overlapping atomic orbitals, which results in an accumulation of charge between the bonded atoms. This accumulation shows as a broad maximum in the ELF along the bonding axis for nonpolar covalent semiconductors such as c-Si [39], for which it reaches a value of 0.95 (Fig. 5). As the bond breaks, charge is no longer localized between the atoms, but in atomic orbitals instead, which shows as peaks in the ELF near the atomic positions, forming a minimum in the center (Fig. 5). This difference in the behavior of the ELF along the axes between neighboring atoms can be used to distinguish whether a bond exists or not. As the different shape also goes along with a distinctly different value of the ELF in the center between the atoms, we can alternatively use this value as a simple criterion for identifying a dangling bond. Even though the values of the maxima and minima vary in amorphous semiconductors, a value of 0.8 has been found to separate the maxima from the minima for all analyzed Si-Si and Si-H bonds in aSi:H, and is therefore used as a threshold throughout this study. This threshold can also be applied to Si-H bonds, for which the ELF is also shown in Fig. 5. The reason why the ELF drops to zero near the Si but not near the H atoms is that the core electrons of Si are not taken into account. By systematically calculating the ELF along the bonding axis between all pairs of neighboring atoms and applying the criterion described above, we find that the small configuration contains exactly one Si atom that is only twofold coordinated, whereas all other Si atoms are perfectly fourfold coordinated, and all the H atoms are bonded to exactly one Si atom (see Table 1). This translates to a defect density of 1.5 × 1021 cm−3 , which is about five orders of magnitude higher than experimentally measured values [40,41]. The large configuration does not contain any dangling bonds, that is, all Si atoms are fourfold and all H atoms are onefold coordinated (see Tables 1 and 2). Thus the ELF approach, with its ability to ascertain the presence of covalent bonds, predicts that the large configuration is defect free and is hence promising from the viewpoint of reproducing the electronic properties of real a-Si:H. The discrepancy between the ELF and the geometrical approach can be explained by the presence of strain. The analysis of the bond lengths yields that the number of compressed and elongated bonds are on average the same. Even though the strain is small it is not negligible and

Large system

Geom.

0.5

Fig. 5. Analysis of bonding by means of the Electron localization function (ELF). This example shows the ELF along the bonding axes for a threefold bonded Si atom in a-Si:H. A maximum in the center between two atoms exceeding a threshold value of 0.8 indicates a bond. A minimum indicates a broken bond. One can distinguish two Si–Si bonds (dashed), one dangling bond (dotted), and one Si–H bond (dash-dotted). The ELF in c-Si is shown as a reference (solid).

Table 2 Environment analysis for Si atoms in both a-Si:H systems. The environment SixHy of an atom is defined by the number of Si atoms x and the number of H atoms y bonded to that atom. The number of Si atoms with a certain environment is computed in two ways: the geometric protocol and the ELF approach. Values different from zero are printed in bold for better readability. The results show that almost all H atoms saturate dangling bonds to form fourfold coordinated Si atoms, whereas there is no sign of H2 dimer formation. Envir

cSi Si−Si bond dangling bond Si−H bond

0.2

Geom.

confirmed by a more precise analysis of the bonds based on the computation of the Electron Localization Function (ELF) [38]. This method will be explained in the next section and the discrepancies with the geometric approach will be discussed. To understand which type of bonding is formed around every Si atom and which role the H atoms play, the environment of each Si atom is analyzed. To this end, we compute the number of SixHy configurations, where x and y are the number of neighboring Si and H atoms of the considered Si atom, respectively. The results are reported in Table 2, together with a comparison between the results obtained from the geometrical and the ELF method. In both systems it is clear that no H2 dimers are formed and no two H atoms are bonded to the same Si atoms, since all the SixH2 environments are equal to zero. On the contrary, almost all the H atoms are saturating dangling bonds of Si atoms to form fourfold coordinated Si atoms. Indeed, for example in the large system, all the 64 H atoms are involved in Si3H1 environments. This analysis confirms that the large system is more ordered than the small one and that our approach is effectively able to produce large and at the same time defect-free a-Si:H configurations.

Coord

0.4

Large system

162

Computational Materials Science 155 (2018) 159–168

P. Czaja et al. 2−fold bonded atom

can induce errors in the computation of the coordination.

4−fold bonded atom

We investigate the effect of the structure and the size of the configurations on the electronic properties in terms of the density and localization of states. As a quantitative measure for the localization of a wave function ψ we use the spread S, which is calculated as the square root of the variance of |ψ|2 with respect to one super cell:

= 2·

∫Ω dr|ψ (r)|2 r 2−(∫Ω dr|ψ (r)|2 r)

,

dr|ψ (r)|2 r,

(5)

11 10.5 10 9.5 9

Spread DOS −2

−1

0 E [eV]

1

Spread [Å]

DOS [1/eV]

which can be interpreted as the position where the wave function is localized. This definition allows us to identify localized states, and to locate them both in real and in energy space. In Fig. 6 the spread is shown as a function of the energy together with the DOS for the small configuration. The DOS reveals two distinct peaks inside the gap that originate from localized states, i.e., states with low spread. In order to understand the origin of these states we use the definition of the wave function center to determine the position of localized states with a spread up to 9.5 Å inside the super cell (Fig. 7). Each dot represents one wave function, where the color represents the spread (top) and the energy (bottom) of that wave function. Additionally, the position of the atoms are symbolized by circles, where the larger circle represents the single twofold bonded atom. The map shows that the strongly localized states inside the gap are localized directly at the under-coordinated atom, and can therefore be assigned to dangling bonds. Moreover, the bottom figure shows that the states below and above the Fermi energy are spatially separated and thus 0.88 eV

2

z

y

−1

originate from different dangling bonds. While the gap states can be easily identified, the distinction of the tail states in Fig. 6 is not so trivial. Even though a tail region with semilocalized states and a band region with mostly extended states can be identified, the relatively small difference in the spread makes the accurate determination of the mobility edges, i.e., the transitions between regions of different localization, very difficult. The way this is usually done is by setting a threshold for the localization that separates localized and extended states [8]. This gives an estimate for the mobility gap, which is defined as the gap between extended valence and conduction band states [42]. The threshold for the spread should be set such that band states far away from the gap exceed this value. Clearly this method becomes more accurate the more unambiguous the choice of the threshold is, and the steeper the localization changes at the mobility edges. By choosing a threshold value of S = 10.1 Å we obtain a mobility gap of roughly 0.9 eV for the small configuration. The difficulty in distinguishing band and tail states is to some extent an artifact of the finite super-cell size. This can be understood by considering that even the band states are not perfectly delocalized plane waves but have regions of higher and lower probability density, leading to a spread smaller than the extension of the super cell. On the other hand, an exponentially localized tail state close to the mobility edge can have a spread that is of the order of the super-cell size, causing an underestimation of the mobility gap. However, when the size of the super cell is increased, the spread of an extended state will increase accordingly, whereas the spread of an exponentially localized state will converge to a finite value. This separates tail and band states, and should lead to convergence of the mobility gap if the cell size is chosen large enough. Another effect of the finite size is a broadening of bands due to the interaction of semi-localized states with their periodic images, which could shift states out of the gap and therefore smear out the transition between band and tail states. The distribution of the spread in the large system (Fig. 8) shows indeed the effects of the larger super-cell size. The relative difference in the spread of the band states is much smaller, and the tails are more pronounced. Also the mobility gap is more distinct and slightly larger with a value of 1.0 eV (obtained with a threshold of S = 20.7 Å ), which is however still small as compared to the experimental value of around 1.9 eV [42]. The DOS does not show any mid-gap states, which is in agreement with the fact that the large configuration does not have any dangling bonds. There are, however, two small peaks at the border of the gap arising from weak bonds, i.e., bonds that show a comparably low value for the ELF between the atoms.

instead. Throughout the rest of the paper the term spread, or S respectively, will therefore refer to Smin . Using the integration volume Ωm that minimizes the spread we can also define a wave function center m

x

Fig. 7. Distribution of localized states with S ⩽ 9.5 Å in the super cell for the small system (projections onto the xy-, xz-, and yz-plane). Each dot represents one wave function, where the colors denote its spread (top) and its energy (bottom) respectively. The circles represent the Si atoms, where different sizes stand for different coordination. Gap states both below and above the Fermi level (0 eV) are strongly localized at the under-coordinated atom, but are spatially separated and can therefore be assigned to two different dangling bonds.

(4)

Ω

−3

E [eV]

y

(3)

Smin = min S

45 40 35 30 25 20 15 10 5 0

0 −0.5

x

∫Ω

8.5

0.5

where we assume that ψ is normalized with respect to the super-cell volume Ω. It can be easily seen that a maximally localized ψ (i.e., a Dirac delta function) gives S = 0 , whereas a wave function that is maximally delocalized over the super cell (i.e., a plane wave) will result in S = a , where a is the lattice constant of the super cell. This explains the factor 2 in the definition. As S is not uniquely defined but will in general depend on the choice of the integration volume Ω, we use the minimal spread

〈r〉 =

y

1

(2) 2

x

z

4(〈r 2〉−〈r〉2 )

S=

x

9

z

y

z

5.2. Energy and localization of wave functions

Spread [Å]

9.5

8.5 3

Fig. 6. Density of states (DOS) and wave-function spread around the Fermi level (0 eV) for the small system. Each dot represents one wave function. The peaks of the DOS inside the gap can be clearly associated with localized states (states with low spread). The horizontal dashed line indicates the threshold for the distinction of extended states from tail states, leading to a definition for the mobility band edges, marked by the vertical dashed lines. The value for the mobility gap obtained this way is 0.88 eV. 163

Computational Materials Science 155 (2018) 159–168

P. Czaja et al.

22 21 20 19 18 17 16 15 14 13 3

0.96 eV

DOS [1/eV]

300 250 200 150 100 50 0

Spread DOS −3

−2

−1

0 E [eV]

1

2

Table 3 Scissors-shift parameters for a-Si:H obtained from linearly fitting the results of the G0W0 calculation in different energy ranges around the Fermi level. The resulting absorption spectra are plotted in Fig. 10. The first parameter set is used in Fig. 9.

Spread [Å]

350

6.1. G0W0 calculations The G0W0 calculation for the small a-Si:H structure yields the quasiparticle corrected electron energies shown in Fig. 9 (inset). The results indicate that the effect of the quasi-particle corrections consists mainly in a spreading of valence and conduction band by approximately 0.27 eV, which suggests that the costly G0W0 calculation can be substituted by a simple scissors shift. A scissors shift is a linear shift QP ε v/c = a v/c ·ε v/c + b v/c of both the valence and the conduction band, where ε v/c are the uncorrected energies [16]. The respective shifting parameters are obtained by a linear fit of the G0W0 results. The choice IP SS GW

bc [eV]

[−1: 1] [−2: 2] [−3: 3]

1.088 1.044 1.025

−1.097 −0.858 −0.762

1.146 1.064 1.043

−1.228 −0.665 −0.509

3

105

2

104

1 ΔεQP = 0.56 eV 0 -1

Δε = 0.29 eV

-2 -3 -3

-2

-1

3

0

0.5

1

1.5

2 E [eV]

2.5

0 ε [eV]

3

1

2

3.5

α [1/cm]

εQP [eV]

α [1/cm]

ac

Fig. 9 shows the absorption spectrum of the small configuration calculated within the independent-particle (IP) approximation (i.e., with the uncorrected Kohn–Sham energies), the GW approximation, and the scissors-shift (SS) approximation. The IP spectrum shows two sub-gap absorption peaks at 0.34 eV and 0.60 eV. By comparison with the DOS the first peak can be related to absorption processes between two gap states, whereas the second peak arises from absorption processes between a gap state and a tail state. In order to estimate the optical gap we use a Tauc plot [43], where the linear regime of αE is extrapolated and Eg is determined as the intersection of the extrapolated line with the energy axis. This yields a value of Eg = 0.7 eV , which is small compared to the experimental value of approximately 1.7 eV [44]. The G0W0 correction modifies the absorption spectrum only in terms of a shift and a slight stretch, which results in a corrected optical gap of 1.0 eV . The figure shows that the G0W0 correction can be well approximated by a scissors shift, where the first parameter set in Table 3 was used. These parameters were chosen because they best approximate the G0W0 absorption spectrum, as can be seen in Fig. 10, where the spectra for all three parameter sets are compared to the G0W0 spectrum. After finding a suitable set of scissors-shift parameters for a-Si:H, we use these parameters to calculate a quasi-particle corrected absorption spectrum also for the large configuration, for which a G0W0 calculation would be too costly. The result is shown in Fig. 11, together with the uncorrected spectrum. As compared to the small configuration, the subgap absorption decreased significantly. Moreover the optical gap obtained from the Tauc plot increased to 1.0 eV in the IP approximation,

6. Optical properties

10

b v [eV]

6.2. Absorption spectrum

Both the results for the small and the large configuration support the use of the ELF for the detection of dangling bonds, as the density and localization of states in both cases matches very well the number of dangling bonds predicted with the ELF approach. In particular, the agreement is much better than with the geometric approach, which predicts a higher number of dangling bonds that is not reflected in the electronic structure.

6

av

of the right set of parameters thereby depends on the energy range of interest. By using different energy ranges for fitting we obtain different parameter sets, which are listed in Table 3.

Fig. 8. Density of states (DOS) and wave-function spread around the Fermi level (0 eV) for the large system. Each dot represents one wave function. The horizontal dashed line indicates the threshold for the distinction of extended states from tail states. The tails are more pronounced than in the small system, leading to more distinct mobility edges, marked by the vertical dashed lines. Unlike in the small system there are no mid-gap states (just some states near the valence band edge arising from weak bonds) and the mobility gap is slightly bigger with 0.96 eV.

10

Fit. range [eV]

3

10

5

4

Fig. 9. Absorption spectrum of the small configuration calculated in independent particle (IP, dotted), GW (dash-dotted), and scissors-shift (SS, solid) approximation. The GW and SS spectrum are based on the quasi-particle corrected electron energies ε QP , which are plotted against the uncorrected energies (inset). Δε refers to the energy difference between the lowest unoccupied and the highest occupied state. The main effect of the corrections is a spreading of the valence and conduction band by 0.27 eV, which appears in the spectrum as a shift of the onset. Qualitatively, all three spectra are similar, showing pronounced sub-gap absorption peaks.

10

GW SS1 SS2 SS3

4

0.5

1

1.5

2

E [eV]

Fig. 10. Comparison of the absorption spectrum of the small configuration calculated in G0W0 approximation (solid) and with the different sets of scissorsshift parameters from Table 3 (SS; dashed, dotted, and dash-dotted). The best agreement is found with the first parameter set. 164

Computational Materials Science 155 (2018) 159–168

P. Czaja et al.

10

0.08

IP SS

6

a)

0.07

JDOS 2 a1/[1+a2/(E-Eg) ] a1/a2(E-Eg)2

10

JDOS [1/eV]

105

6

105

10

10

10

4

10

3

3

0

0.5

1

1.5

0.05 0.04

Eg = 0.76 eV

0.03 0.02

4

0.01 0 2 Average coupling 〈|vcv| 〉 [a.u.]

α [1/cm]

0.06

0 0.5 1 1.5 2 2.5 3 3.5 4

2 E [eV]

2.5

3

3.5

4

Fig. 11. Absorption spectrum of the large configuration calculated in IP (dotted) and SS (solid) approximation, using the SS parameters obtained from the small system. Compared to the small configuration the optical gap is increased by 0.3 eV, whereas sub-gap absorption is reduced. Inset: Absorption spectra averaged over 10 large configurations. The optical gap remains unaffected but the sub-gap peaks vanish.

350

0

1

2

b)

3

300 250 200 150 2

100

〈|vcv| 〉 3 2 b1/[(E-Eg) +b2/(E +b3)] 2 b1/b2(E +b3)

50 0

0

1

2

3

∑ 〈|vcv |2 〉 (E ) =

On the microscopic level the optical properties of a material are determined by the probabilities of valence band electrons being excited to a conduction state under photon absorption. These probabilities depend on the optical coupling between the initial and the final states, given by the optical transition matrix element. The matrix element between a valence state |kv〉 and a conduction state |kc〉 is defined as vcv = 〈kc|e ·v|kv〉, where v = i [H , r] is the velocity operator and e the polarization vector. The sum over the transition probabilities |vcv |2 between all pairs of valence and conduction states defines the imaginary part of the dielectric function

|vcv |2 δ (Ec−Ev−E ),

cv

∊2 (E ) ∝

|vcv |2 δ (Ec−Ev−E )

cv

∑

δ (Ec−Ev−E )

,

J (E ) =

〈 |vcv |2 〉 (E ) J (E ) , E2

∑

δ (Ec−Ev−E ).

cv

(9)

(10)

Eq. (9) means that the energy dependence of both the DOS/JDOS and the coupling strength has to be understood in order to understand and model absorption in a-Si:H and to determine the optical band gap. To begin with, we investigate the JDOS of the large configuration, shown in Fig. 12(a). In the energy range up to 5 eV displayed here, J can be well described by a function of the form

(6)

which is related to the absorption via

E ∊2 (E ) n (E )

(8)

with the joint density of states (JDOS)

J (E ) ≈

α (E ) ∝

5

which describes the average transition probability for all transitions with a given energy E, Eq. (6) can be written as

6.3. Optical transition matrix elements

∑

4

Fig. 12. (a) Normalized joint density of states (JDOS) for the large configuration. Fitting a function of the form a1/[1 + a2 /(E −Eg )2] (red dashed line) yields a value of Eg = 0.76 eV for the optical gap. Close to Eg (in a range of about 1 eV) the JDOS can be approximated quadratically (purple dotted line), which is the common assumption being made when determining the optical gap. (b) Average optical coupling strength 〈|vcv |2 〉 as a function of the transition energy. Above ≈ 0.9 eV (i.e., at energies where transitions between localized states do not play a role) the average coupling can be described by a function of the form b1/[(E −Eg )3 + b2 /(E 2 + b3)] (red dashed line). This means that close to Eg (in a range of about 1 eV) the average coupling has the form 〈|vcv |2 〉 (E ) ∝ E 2 + b3 (purple dotted line), then assumes a maximum at E ≈ 2.9 eV , and falls off like 1/ E 3 at high energies. Below 0.9 eV the coupling is strongly increased due to transitions between localized states (cf. Fig. 14). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

cv

1 E2

5

E [eV]

and to 1.3 eV with scissors-shift corrections. While the relation between a larger super cell and a decreased sub-gap absorption can be easily understood in terms of the reduced defect density and the reduced spatial overlap of localized states, we do not have a straightforward explanation for the dependence of the optical gap on the super-cell size. The results suggest however that finite-size effects do have an influence on Eg , and that a closer agreement with the experimental value could be reached with even larger super cells. The spectra shown so far were obtained from single configurations and have therefore limited physical significance. In order to obtain a physically meaningful absorption spectrum of a-Si:H, i.e., a spectrum that can be compared to experimental data, the configurational average has to be taken. For that purpose we calculated spectra for 10 different large configurations and averaged over them. The result is shown in Fig. 11 (inset). While sub-gap absorption is still present in the averaged spectrum, the distinct peaks disappeared. This resembles more the experimental findings [44,45], even though the contribution of sub-gap absorption is still overestimated. The value of the optical gap is not affected by the averaging.

∊2 (E ) ∝

4

E [eV]

a1 , 1 + a2 /(E −Eg )2

(11)

where the ai are fitting parameters. Fitting this function to the data yields a value of Eg = 0.76 eV for the band gap. Close to Eg the function becomes approximately a1/ a2 (E −Eg )2 . This agrees with the commonly made assumption that the density of band states increases as E from the band edges [43], resulting in a quadratic energy dependence for the JDOS. The approximation however only holds in a small energy range

(7)

with n denoting the refractive index. By defining the average coupling strength 165

Computational Materials Science 155 (2018) 159–168

P. Czaja et al.

400

〈 |vcv |2 〉 (E ) ≈

(E −Eg

)3

b1 , + b2 /(E 2 + b3)

Average coupling 〈|vcv|2〉 [a.u.]

of about 1 eV above the gap. The average optical coupling strength is shown in Fig. 12(b). Above ≈ 0.9 eV , where transitions between localized states do not play a role, the data can be described by

(12)

where the bi are fitting parameters. In the vicinity of Eg this function becomes b1/ b2 (E 2 + b3) , where b3 is of the order of 3 eV 2 . This result clearly contradicts the common assumption that 〈|vcv |2 〉 is constant at low energies [43]. This was disproven experimentally already by Jackson et al. [46], who however stated that instead 〈|vcv |2 〉 ∝ E 2 , which also disagrees with our results. However, the energy range in which an E 2 behavior is displayed in the experiment is rather small, which, combined with the comparably large error estimates in the fitting range, would allow also alternative interpretations of the measurement results. It therefore remains unclear if the disagreement in E dependence is due to an actual discrepancy between measurement and simulation, or arises rather from a different interpretation of the results. Good agreement with the experimental findings on the other hand is reached at high energies, where we find that 〈|vcv |2 〉 falls off like 1/ E 3 . Also the position of the peak near the direct c-Si band gap, which is ≈ 2.6 eV in GGA-DFT, agrees with the experimental result. We will now discuss the determination of the optical band gap. The assumptions behind the Tauc fit are J (E ) ∝ (E −Eg )2 and 〈|vcv |2 〉 constant near the gap, i.e., ∊2 ∝ (E −Eg )2 / E 2 and therefore E ∊2 becomes linear and yields the optical gap. From Eqs. (11) and (12) we see however that E ∊2 ∝ (E −Eg ) E 2 + b3 near Eg , which increases super-linearly (unless b3 ≫ Eg2 , which is not the case here). As shown in Fig. 13, a linear fit ∼ to E ∊2 will therefore always yield a gap Eg that is larger than the value Eg obtained from the JDOS, where the difference increases with the fitting range. This is supported by the experimental results of Jackson et al. [46] who obtain 1.82 eV from fitting the DOS, and 1.86 eV from ∼ fitting E ∊2 . The discrepancy between Eg and our Tauc gap of 1.0 eV reported above can be explained with the fact that αE is only approximately proportional to E ∊2 because the refractive index n in Eq. (7) is also energy dependent. Also, a larger fitting range of about 3 eV was used in the Tauc fit, which is problematic because the quadratic approximation for the JDOS holds only up to ≈ Eg + 1 eV as stated before. A Tauc fit would in principle make sense therefore only within that range. The reason why it still gives sensible results is that the errors in the assumptions for J, 〈|vcv |2 〉, and n approximately compensate [46].

ext-ext loc-ext loc-loc

350 300 250 200 150 100 50 0

0

1

2

3

4

5

E [eV]

Fig. 14. Average optical coupling strength 〈|vcv |2 〉 resolved by type of state. A spread of S = 20.7 Å is used to distinguish localized from extended states. The error bars indicate the standard deviation of the average. The plot shows no difference in the qualitative behavior for the coupling of two extended states (solid purple line), and the coupling between an extended and a localized state (blue dashed line), respectively. The coupling between two localized states (orange dotted line) is significantly stronger only at E ≲ 0.9 eV and is the sole reason for the increase of the total coupling seen in Fig. 12b. At high energies its error becomes too large to make any statement about the energy dependence. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Another point of interest in the context of understanding and modeling absorption is the question, whether different types of transition exhibit different behavior, and if the transition probability also depends on the properties of the initial state. To address these issues, states are separated in localized and extended states using a threshold for the spread of S = 20.7 Å . 〈|vcv |2 〉 is then calculated separately for each kind of transition: extended to extended, localized to extended and vice versa, and localized to localized. The result is shown in Fig. 14. We see that there is no qualitative difference for extended-extended and extended-localized transitions, in agreement with the experimental findings [46]. Localized-localized transitions however show a significantly higher coupling below ≈ 0.9 eV . As these transitions dominate at low energies, this also leads to a large increase of the total coupling, as can be seen in Fig. 12(b). The anomalous energy dependence of the coupling between localized states, which is significantly enhanced, but only at low energies, suggests that the transition probability does not only depend on the energy of the transition, but also on the respective states. In order to investigate this dependence, we generalize the definition (8) for the average coupling such that, instead of restricting the average to all transitions with a given transition energy E, we average over all transitions with any given property. As an example, 〈|vcv |2 〉 is calculated as a function of the valence and the conduction state energy (Fig. 15). The plot shows that, even though brighter spots do exist at the highest valence state, there is no general dependence on Ev and Ec , but only on the energy difference E. On the other hand, a correlation with the localization can be found when plotting 〈|vcv |2 〉 as a function of the valencestate spread and the transition energy. For strongly localized valence states the coupling seems to be much stronger than for weakly localized and extended states, where 〈|vcv |2 〉 (Sv ) becomes approximately independent of Sv . Due to the small number of these states, this effect however becomes visible in the average 〈|vcv |2 〉 (E ) only at very low energies, where transitions involving strongly localized states dominate. The observation made above, that the probabilities for extendedextended and extended-localized transitions do not differ on average, therefore remains valid.

Fig. 13. Determination of the optical gap from the dielectric function ∊2 . According to the fits of J and 〈|vcv |2 〉 (Fig. 12), E ∊2 (blue solid line) assumes a form ∝ (E −E g) E 2 + b3 close to Eg (red dashed line). The determination of the optical gap from a linear fit, which implies E ∊2 ∝ E −Eg (purple dotted line), ∼ therefore leads to a slightly larger value Eg = 0.83 eV . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 166

Computational Materials Science 155 (2018) 159–168

P. Czaja et al.

dependence suggested in the literature. The authors gratefully acknowledge funding from the European Commission Horizon 2020 research and innovation program under grant agreement No. 676629, support through the COST action MP1406 MultiscaleSolar, as well as the computing time granted on the supercomputers JURECA [47] at Jülich Supercomputing Centre and CRESCO [48] on the ENEA-GRID infrastructure. The processed data required to reproduce these findings are available to download from https://data.mendeley.com. 8. CRediT authorship contribution statement Simone Giusepponi and Michele Gusso performed the ab initio molecular dynamics simulations to generate the atomic structure of the a-Si:H configurations and analyzed their quality. Philippe Czaja performed the DFT and GW computations to characterize the electronic structure and optical properties of the material. Urs Aeberhard and Massimo Celino conceived the idea for the work and supervised its execution. All authors contributed to the writing of the manuscript.

Fig. 15. Left.: 〈|vcv |2 〉 as a function of the conduction state energy Ec and the valence state energy E v (with respect to the Fermi level). The white diagonal indicates the line of constant transition energy E where 〈|vcv |2 〉 (E ) has a maximum. Along the lines of constant E the coupling does not vary much, showing that the coupling between two states does not depend significantly on the energies E v and Ec , but mostly on the difference E = Ec−E v . Right: 〈|vcv |2 〉 as a function of the valence state spread Sv and the transition energy E. Coupling is particularly strong at all energies for strongly localized states, whereas for weakly localized and extended states it becomes almost independent from the localization.

References [1] [2] [3] [4] [5]

7. Conclusions and outlook We succeeded in obtaining large-scale defect free a-Si:H configurations from ab initio molecular dynamics, which are a good approximation to real a-Si:H structures and are therefore an ideal starting point for the electronic and optical characterization of a-Si:H from first principles. The structural, electronic, and optical properties of the generated configurations were calculated and subsequently analyzed. We found that the finite size of the super cell affects both the localization of the electronic states and the absorption spectrum. In particular, a larger super cell improved the values we calculated for the mobility gap and the optical gap, even though the experimental values are still underestimated. Furthermore, we successfully performed G0W0 calculations for an a-Si:H configuration with 72 atoms and found that the quasi-particle corrections can be approximated by a scissors shift. This approximation makes calculations for larger – and thus physically more representative – configurations possible. The extracted set of scissors-shift parameters was used for calculating quasi-particle corrected absorption spectra for all configurations, which improved the optical gap by roughly 0.3 eV, but could not completely bridge the difference to the experimental value. This discrepancy is in contradiction to the good agreement of the structural properties, which suggests that it could be a computational artifact. One candidate for causing this are the finite size effects, which were demonstrated and discussed qualitatively, and which might be eliminated by making the super cell large enough. The other potential source of error is the negligence of the quasi-particle effects on the wave function in the G0W0 calculations, where only the energies are corrected. This error might be eliminated by performing fully self-consistent GW calculations instead. Finally, we analyzed the optical transition matrix elements for the large configuration. Our results agree with experiment in that the average transition probability decreases with 1/ E 3 at high energies, and that it is the same for transitions between two extended states and transitions between a localized and an extended state. Nevertheless, it seems that strong localization actually does increase the transition probability, but the effect is not visible on average due to the small number of strongly localized states. Concerning the often discussed energy dependence of the matrix elements near the gap, we find a dependence E 2 + c with some constant c. This disagrees with the common assumption of constant matrix elements, but also with the E 2

[6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20]

[21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44]

167

D.E. Carlson, C.R. Wronski, Appl. Phys. Lett. 28 (1976) 671–673. D.C. Allan, J.D. Joannopoulos, W.B. Pollard, Phys. Rev. B 25 (1982) 1065–1080. P.A. Fedders, D.A. Drabold, Phys. Rev. B 47 (1993) 13277–13282. B. Tuttle, J.B. Adams, Phys. Rev. B 57 (1998) 12859–12868. A.A. Valladares, F. Alvarez, Z. Liu, J. Sticht, J. Harris, Eur. Phys. J.B 22 (2001) 443–453. K. Jarolimek, R.A. de Groot, G.A. de Wijs, M. Zeman, Phys. Rev. B 79 (2009) 155206. P.A. Khomyakov, W. Andreoni, N.D. Afify, A. Curioni, Phys. Rev. Lett. 107 (2011) 255502. M. Legesse, M. Nolan, G. Fagas, J. Phys. Chem. C 117 (2013) 23956–23963. W. Shockley, W.T. Read, Phys. Rev. 87 (1952) 835–842. R.A. Street, Hydrogenated Amorphous Silicon, Cambridge University Pr., Cambridge, 1991. M. Durandurdu, D.A. Drabold, N. Mousseau, Phys. Rev. B 62 (2000) 15307–15310. V. Rosato, M. Celino, J. Appl. Phys. 86 (1999) 6826–6834. M. Nolan, M. Legesse, G. Fagas, Phys. Chem. Chem. Phys. 14 (2012) 15173–15179. J.P. Perdew, Int. J. Quantum Chem. 28 (1985) 497–523. M.S. Hybertsen, S.G. Louie, Phys. Rev. Lett. 55 (1985) 1418–1421. R.W. Godby, M. Schlüter, L.J. Sham, Phys. Rev. B 37 (1988) 10159–10175. P. Hohenberg, W. Kohn, Phys. Rev. 136 (1964) B864–B871. W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133–A1138. J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G.L. Chiarotti, M. Cococcioni, I. Dabo, A.D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A.P. Seitsonen, A. Smogunov, P. Umari, R.M. Wentzcovitch, J. Phys.: Condens. Matter 21 (2009) 395502. Quantum ESPRESSO, http://www.quantum-espresso.org. CP2K, http://www.cp2k.org. S. Goedecker, M. Teter, J. Hutter, Phys. Rev. B 54 (1996) 1703–1710. C. Hartwigsen, S. Goedecker, J. Hutter, Phys. Rev. B 58 (1998) 3641–3662. M. Krack, Theoret. Chem. Acc. 114 (2005) 145–152. L. Hedin, Phys. Rev. 139 (1965) A796–A823. J. Deslippe, G. Samsonidze, D.A. Strubbe, M. Jain, M.L. Cohen, S.G. Louie, Comput. Phys. Commun. 183 (2012) 1269–1289. B.I. Lundqvist, Physik der kondensierten Materie 6 (1967) 193–205. A.W. Overhauser, Phys. Rev. B 3 (1971) 1888–1898. M.S. Hybertsen, S.G. Louie, Phys. Rev. B 34 (1986) 5390–5413. P. Larson, M. Dvorak, Z. Wu, Phys. Rev. B 88 (2013) 125205. H. Ehrenreich, The Optical Properties of Solids, Academic, New York, 1965. R. Car, M. Parrinello, Phys. Rev. Lett. 60 (1988) 204–207. E. Johlin, L.K. Wagner, T. Buonassisi, J.C. Grossman, Phys. Rev. Lett. 110 (2013) 146805. H.C. Andersen, J. Chem. Phys. 72 (1980) 2384–2393. S. Nosé, J. Chem. Phys. 81 (1984) 511–519. R. Bellisent, A. Chenevas-Paule, P. Chieux, A. Menelle, J. of Non-Cryst. Solids 77–78 (1985) 213–216. A.D. Becke, K.E. Edgecombe, J. Chem. Phys. 92 (1990) 5397–5403. A. Savin, O. Jepsen, J. Flad, O.K. Andersen, H. Preuss, H.G. von Schnering, Angew. Chem. Int. Ed. English 31 (1992) 187–188. M. Favre, H. Curtins, A. Shah, J. Non-Cryst. Solids 97 (1987) 731–734. R.A. Street, K. Winer, Phys. Rev. B 40 (1989) 6236–6249. K. Winer, Annu. Rev. Mater. Sci. 21 (1991) 1–21. J. Tauc, R. Grigorovici, A. Vancu, Phys. Status Solidi (b) 15 (1966) 627–637. G.D. Cody, T. Tiedje, B. Abeles, B. Brooks, Y. Goldstein, Phys. Rev. Lett. 47 (1981)

Computational Materials Science 155 (2018) 159–168

P. Czaja et al.

Cucurullo, P. Dangelo, M.D. Rosa, P.D. Michele, A. Funel, G. Furini, D. Giammattei, S. Giusepponi, R. Guadagni, G. Guarnieri, A. Italiano, S. Magagnino, A. Mariano, G. Mencuccini, C. Mercuri, S. Migliori, P. Ornelli, S. Pecoraro, A. Perozziello, S. Pierattini, S. Podda, F. Poggi, A. Quintiliani, A. Rocchi, C. Sciò, F. Simoni, A. Vita, in: High Performance Computing Simulation (HPCS), 2014 International Conference on, 2014, pp. 1030–1033.

1480–1483. [45] C. Tsai, H. Fritzsche, Solar Energy Mater. 1 (1979) 29–42. [46] W.B. Jackson, S.M. Kelso, C.C. Tsai, J.W. Allen, S.-J. Oh, Phys. Rev. B 31 (1985) 5187–5198. [47] Jülich Supercomputing Centre, J. Large-scale Res. Facil. 2 (2016) A62. [48] G. Ponti, F. Palombi, D. Abate, F. Ambrosino, G. Aprea, T. Bastianelli, F. Beone, R. Bertini, G. Bracco, M. Caporicci, B. Calosso, M. Chinnici, A. Colavincenzo, A.

168