Temperature-induced nanostructural evolution of hydrogen-rich voids in amorphous silicon: a first-principles study

Parthapratim Biswas *a, Durga Paudel a, Raymond Atta-Fynn b and Stephen R. Elliott c
aDepartment of Physics and Astronomy, The University of Southern Mississippi, Hattiesburg, MS 39406, USA. E-mail: partha.biswas@usm.edu
bDepartment of Physics, University of Texas at Arlington, Arlington, TX 76019, USA
cDepartment of Chemistry, University of Cambridge, Cambridge, CB2 1EW, UK

Received 24th September 2019 , Accepted 14th November 2019

First published on 14th November 2019

The paper presents an ab initio study of temperature-induced nanostructural evolution of hydrogen-rich voids in amorphous silicon. By using large a-Si models, obtained from classical molecular-dynamics simulations, with a realistic void-volume density of 0.2%, the dynamics of Si and H atoms on the surface of the nanometer-size cavities were studied and their effects on the shape and size of the voids were examined using first-principles density-functional simulations. The results from ab initio calculations were compared with those obtained from using the modified Stillinger–Weber potential. The temperature-induced nanostructural evolution of the voids was examined by analyzing the three-dimensional distribution of Si and H atoms on/near void surfaces using the convex-hull approximation, and computing the radius of gyration of the corresponding convex hulls. A comparison of the results with those from the simulated values of the intensity in small-angle X-ray scattering of a-Si/a-Si:H in the Guinier approximation is also provided, along with a discussion on the dynamics of bonded and non-bonded hydrogen in the vicinity of voids.

I. Introduction

Amorphous silicon and its hydrogenated counterpart have a wide range of applications, from photovoltaics to thin-film technology.1,2 Thin films of hydrogenated amorphous silicon (a-Si:H) are extensively used for passivation of crystalline silicon (c-Si) surfaces, which is essential to produce a high open-circuit voltage in silicon-based heterojunction solar cells.3,4 The passivation of Si dangling bonds by H atoms in the vicinity of a-Si:H/c-Si interfaces plays a key role in improving the properties of silicon heterostructures.5 Post-deposition annealing is routinely used for structural relaxation of a-Si:H samples and incorporation of hydrogen at the a-Si:H/c-Si interfaces. The presence of hydrogen to an extent thus is often preferred6 for the preparation of a-Si:H samples with good interface properties. Experimental studies, such as spectroscopic ellipsometry7,8 and Raman spectroscopy,9 indicate that, although the presence of nanometer-size voids can degrade the electronic quality of the samples, a small concentration of voids can improve the degree of surface passivation. This is generally believed to be due to the presence of mobile H atoms in a void-rich environment, where hydrogen can reach a-Si:H/c-Si interfaces, via diffusion or other mechanisms, in comparison to dense a-Si:H samples.

Heat treatment or annealing is an effective tool for nanostructural relaxation of laboratory-grown samples. Annealing at low to medium temperature (400–700 K) can considerably improve the quality of as-deposited samples, by removing unwanted impurities and reducing imperfections (e.g., defects) in the samples.10 For amorphous materials, annealing also increases the degree of local ordering in the amorphous environment. Earlier studies by small-angle X-ray scattering (SAXS) reported the presence of columnar-like geometric structure11 and blister formation in annealed samples of a-Si:H.12 Likewise, experiments on a-Si/a-Si:H, using positron-annihilation lifetime spectroscopy, have indicated that annealing above 673 K can lead to the formation of voids via vacancy clustering, which affects the nanostructural properties of the network.13,14 It has been suggested that at higher temperature, near 800 K, a considerable restructuring can take place in the vicinity of the void boundary that can modify the shape and size of the voids.15 Although tilting SAXS15,16 can provide some information on the geometrical shape of the voids, a direct experimental determination of annealing effects on the morphology of voids in the amorphous environment is highly nontrivial, as scattering experiments generally provide integrated or scalar information on the size and shape of the voids. By contrast, computational modeling can yield reliable structural information on the morphology of voids,17,18 provided that high-quality large models of a-Si/a-Si:H – with a linear dimension of several nanometers – are readily available and accurate total-energy functionals to describe the atomic dynamics of Si and H atoms are in place. However, despite numerous studies on computational modeling of a-Si and a-Si:H over the past several decades, there exist only a few computational studies17–21 that truly attempted to address the structure and morphology of voids on the nanometer length scale in recent years.

The purpose of the present study is to conduct ab initio simulations of hydrogen-rich voids in large models of a-Si at low and high annealing temperatures. We are particularly interested in the structural reconstruction of the void surfaces, induced by annealing, and the resulting effects on the intensity in small-angle X-ray scattering (SAXS). Further, we intend to study the evolution of the shape of voids with temperature in the presence of hydrogen inside the voids. To this end, we develop significantly large models of a-Si, which are characterized by realistic distributions of voids as far as the results from SAXS,16,22 nuclear magnetic resonance (NMR),23 infrared (IR),24,25 and positron-annihilation26 measurements are concerned. The term realistic here refers to the fact that the models must be able to accommodate nanometer-size voids with a void-volume fraction of 0.2–0.3%, as observed in SAXS, IR, and NMR measurements. The three-dimensional shapes of the voids are studied by using the convex-hull approximation17 of the boundary region of the voids, upon annealing at 300 K and 800 K, and the effects of the shape and size of the voids on the simulated SAXS intensity from the models are examined. The presence of molecular hydrogen and the statistics of various silicon-hydrogen bonding configurations near the void surfaces are also examined by introducing H atoms within the voids, with a number density of hydrogen (inside cavities) as observed in infrared spectroscopy and NMR measurements.24,27

The rest of the paper is as follows. In section II, we briefly describe the simulation of large models of a-Si with a linear dimension of a few tens of nanometers and the generation of nanometer-size voids in the resulting models in order to incorporate a realistic distribution of voids, with a void-volume fraction as observed in SAXS and IR measurements. Section III discusses the results of our work, with particular emphasis on the three-dimensional shape and size of the voids upon annealing at 300 K and 800 K. The correlation between the three-dimensional shape and size of the voids and the intensity of SAXS from simulations is discussed here. A comparison of the results with those from positron-annihilation lifetime (PAL) spectroscopy13 and Doppler-broadening positron-annihilation spectroscopy (DB-PAS)28 are presented. This is followed by a discussion on the atomic dynamics of hydrogen and the time evolution of monohydride silicon-hydrogen (Si–H) bond formation in the vicinity of the voids, using ab initio molecular-dynamics (AIMD) simulations over several picoseconds. Section IV presents the conclusions of our work.

II. Computational methods

A. Generation of large models of a-Si

We began by generating a large model of amorphous silicon in order to be able to study the intensity due to SAXS in the small-angle region of scattering wavevector k ≤ 1.0 Å−1. Since the key purpose of the study is to examine the structural evolution of the voids in a-Si on the nanometer length scale, by joint analyses of simulated SAXS intensities and three-dimensional distributions of atoms near void surfaces, it is necessary to simulate a sufficiently large model that can accommodate a realistic distribution of nanometer-size voids with the correct void-volume fraction, in the range from 0.01% to 0.3%, as observed from SAXS and IR measurements.15,16 To this end, a model a-Si network, consisting of 262[thin space (1/6-em)]400 atoms, was generated using classical molecular-dynamics (MD) simulations, followed by geometry relaxation, by employing the modified Stillinger–Weber potential.29,30 The initial configuration was obtained by randomly generating atomic positions in a cubic box of length 176.123 Å, which corresponds to the experimental mass density, 2.24 g cm−3, of a-Si.31 Thereafter, constant-temperature MD simulations were performed using a chain of Nosé–Hoover thermostats32,33 and a time step of Δt = 1 fs. During the MD runs, the system was heated to a temperature of 1800 K for 10 ps, and then it was cooled from 1800 K to 300 K at the rate of 10 K ps−1. The heating–cooling cycle was carried out for at least 20 iterations. For each cycle, structural properties of the system, such as the coordination-number statistics, the bond-angle distribution, and the structure factor of the evolving configurations were examined. The procedure continues until it satisfies the convergence criteria that entail imposing a restriction on the upper limit of the root-mean-square (RMS) width of the bond-angle distribution (Δθmax ≤ 10°) and the coordination-defect density (of 3% or less) in the network. It may be noted that these variables are in conflict with each other, meaning that one can be reduced at the expense of increasing the other, and vice versa. Thus, care was taken to ensure that these criteria were simultaneously satisfied. Once the structural properties satisfied the aforementioned convergence criteria, the MD run was terminated and the system was relaxed by minimizing the total energy of the system with respect to the atomic positions. The resulting configuration was then used as a structural basis for the generation of additional models of a-Si by introducing voids in the network. The details of the simulation procedure and an analysis of the structural and electronic properties of these large a-Si models can be found in ref. 34.

B. Void geometry and distribution in a-Si

Having obtained a large model of a-Si, we proceed to generate a realistic distribution of voids. Since computational modeling of void formation in a-Si, by mimicking the actual deposition and growth processes, is highly nontrivial, we generated a void distribution using experimental information on the shape, size, and the number density of voids. Experimental data from an array of measurements, such as small-angle scattering (SAS),16,22,35,36 infrared (IR) spectroscopy,24 nuclear magnetic resonance (NMR),23 positron-annihilation lifetime (PAL) spectroscopy,26 and calorimetry measurements,27,37,38 suggest that the linear size of the voids in a-Si typically lies between 10 Å and 40 Å and that the void-volume fraction can range from 0.01% to 0.3% of the total volume, depending upon the growth process, deposition rate, and the consequent electronic quality of the samples. Following these observations, we constructed a void distribution in the large model with the following properties: (a) the voids are spherical in shape and have a diameter of 12 Å; (b) the number density of the voids corresponds to a void-volume fraction of 0.2%; (c) the voids are sparsely and randomly distributed, so that the distance between any two void centers is significantly larger than the dimension of the voids.

Fig. 1 shows an example of a void embedded in a two-dimensional amorphous network. The void region is indicated by the gray circle of radius Rv, which is surrounded by a layer of Si atoms (yellow) with a thickness d. Silicon atoms in the layer constitute the surface/wall of the void, which we shall refer to as the void-surface atoms. The remaining Si atoms, within the cubic box (of length l), will be referred to as bulk Si atoms and are indicated as black and blue circles. To study the effect of hydrogen on the evolution of void surfaces, H atoms were introduced inside the voids so that no two H atoms were at a distance less than 1 Å from each other. The presence of H atoms in the network requires that the problem must be treated at the quantum-mechanical level using, for example, first-principles density-functional theory. However, given the size of the large model, such a task is hopelessly difficult and computationally infeasible, and an approximation of some sort is necessary. Here, we employed the local approach to electronic-structure calculations by invoking the principle of nearsightedness of an equilibrium system, as proposed by Kohn,39 and address the problem by studying an appropriately large sub-system, defined by a cubical box of length l surrounding each void (see Fig. 1). Following an (approximate) invariance theorem on electronic structure of disordered alloys, due to Heine and coworkers,40,41 the quantum local density of states near the voids is essentially independent of the boundary conditions at a sufficiently large distance,42 for example, the presence of a fixed layer of Si atoms (blue) on the edge of the box. Referring to Fig. 1, we chose Rv = 6 Å, d = 3 Å, and l = 36 Å. To ensure that the void-volume fraction is consistent with the experimentally observed value of 0.01%–0.3%, we created 12 voids, reflecting a void-volume fraction of 0.2%, which were randomly distributed in the large a-Si model of linear dimension 176.123 Å. An outer layer of Si atoms (blue) of thickness 3 Å was held fixed at the edge of the box of length l for each void during first-principles simulations and total-energy relaxation of the sub-system so that the latter can be put back into the original large model for the calculation of the SAXS intensity, upon completion of annealing and geometry relaxation without boundary perturbations. To ensure that the SAXS intensity from the large a-Si sample with voids is not affected by the atoms on the fixed boundary layer (shown in blue color in Fig. 1) of the sub-systems or supercells, the (undercoordinated) atoms on the boundary layer were left undisturbed during simulations. All other atoms were allowed to move during classical and quantum-mechanical simulations. Each sub-system, consisting of a cubical box of length 36 Å with a central void region of radius 6 Å, contained about 2200 atoms. The size of the cubical box was chosen carefully so that the local structural and electronic properties of the atoms in the vicinity of void surfaces were not affected by the boundary layer on the edge of the box, at a distance of 12 Å from the void surfaces. Care was also taken to ensure that the centers of the voids were sparsely distributed in the original large model so that the individual supercells did not overlap with each other. Specifically, the positions of the center of the voids were generated randomly in such a way that the center-to-center distance between any two neighboring voids in the large model is at least 36 Å or more. To study the effect of hydrogen on void surfaces, each void was loaded with 30 H atoms, as observed in IR and NMR measurements.43 The initial distributions of the H atoms were generated to satisfy the following two criteria: (a) the H atoms were randomly distributed so that no two H atoms were separated by a distance less than 1 Å; (b) the distance between any H atoms inside a void and a Si atom on the void surface was greater than or equal to 2 Å. The latter was achieved by confining the H atoms within a sphere of radius 4 Å from the center of the void. These criteria prevent the system from developing a sudden large force on light H atoms and an instantaneous formation of silicon-hydrogen bonds on void surfaces or shooting off H atoms from the void region at the beginning of AIMD simulations.

image file: c9nr08209c-f1.tif
Fig. 1 A schematic representation of a void (gray region) in two dimensions. The annular region with Si atoms (yellow circles) indicates the void boundary, with a few H atoms (red circles) inside the void. The bulk Si atoms are shown as black and blue circles – the latter indicate a layer of fixed Si atoms for the purpose of local AIMD simulations, as discussed in the text.

C. Ab initio simulations of voids in a-Si

In the preceding section, we have generated a number of models (supercells/subsystems), containing as many as 2200 atoms and a central void of radius 6 Å, using the modified Stillinger–Weber (SW) potential. However, since the SW potential may not accurately describe the chemistry of silicon surfaces, it is necessary to treat these models using a first-principles total-energy functional from density functional theory (DFT) in order to reduce artifacts, such as the number of excessive dangling bonds or defects, associated with the SW-generated void surfaces. The standard DFT protocol here is to compute the self-consistent-field (SCF) solution of the Kohn–Sham (KS) equation, and the presence of voids or inhomogeneities in the system suggests that the generalized gradient approximation (GGA) should be employed in order to deal with the atomic density fluctuations near void surfaces. However, it is simply not feasible to solve the KS equations self-consistently for a problem of this size (i.e., about 2200 atoms) and complexity using the existing plane-wave-based/local-basis DFT methodology. In particular, the computational cost for conducting SCF AIMD runs for several picoseconds is prohibitively large, even for a single configuration. To reduce the computational complexity of the problem, it is thus necessary to employ a suite of approximations for solving the KS equation. Here, we invoke the non-self-consistent Harris-functional approach44,45 in the local density approximation (LDA), using a minimal set of basis functions for Si atoms. While these approximations may affect the accurate determination of some aspects of the electronic structure of void surfaces, the results from test calculations on 1000-atom models with a void of radius 5 Å, obtained via the full SCF procedure for solving the KS equation in the LDA and GGA, indicate that the structural properties of the void surfaces in a-Si are minimally affected by these approximations. In particular, the shape and size of the voids, as well as the values of the average bond angle and bond length of the silicon atoms on the reconstructed void surface, exhibit a very little to no changes with the varying level of DFT approximations (see ESI for results from test calculations). A discussion on the accuracy and computational efficiency of the Harris-functional approach to treat bulk a-Si and a-Si:H can be found in ref. 46.

Ab initio calculations and total-energy optimization were conducted using the first-principles density-functional code SIESTA.47 The latter employs a localized basis set and norm-conserving Troullier–Martins pseudopotentials,48 which are factorized in the Kleinman–Bylander49 form to remove the effect of core electrons. Electronic and exchange correlations between electrons are handled using the LDA, by employing the Perdew–Zunger parameterization of the exchange–correlation functional.50 As stated earlier, we employed the Harris-functional approach44,45 that involves the linearization of the Kohn–Sham equations for structural relaxation and AIMD simulations. The latter were conducted over a period of 10 ps, with a time step of 1 fs. Single-zeta (SZ) and double-zeta-polarized (DZP) basis functions were employed for Si and H, respectively. Here, SZ functions correspond to one s and three p orbitals for Si atoms, whereas DZP functions refer to two s and three p orbitals for H atoms. The size of the simulation box or supercell is found to be large enough in this work so that the Brillouin zone collapses onto the Γ point, which was used for the Brillouin-zone sampling.

D. Simulation of SAXS intensity and void-size determination

The SAXS intensity for a binary system can be expressed with the aid of the partial structure factors and the atomic form factors of the constituent atoms. For a binary system, the scattering intensity, I(k), can be expressed as,51
image file: c9nr08209c-t1.tif(1)
image file: c9nr08209c-t2.tif(2)
is the partial structure factor associated with the atoms of type i and j such that Sij = Sji, and gij(r) is the partial pair-distribution function. Here, image file: c9nr08209c-t3.tif and fi are the atomic fraction and the atomic form factor of the atoms of type i, respectively. In eqn (2), the partial pair-distribution function, gij(r), is defined as gij(r) = nij/nj, where nij is the average number density of atoms of type j with an atom of type i at the center, nj = xjn0, and n0 = N/V. Once the scattering intensity is available, either numerically or experimentally, the size of the voids, or any extended inhomogeneities, can be estimated by invoking the Guinier approximation52 in the small-angle limit, kRG ≤ 1,
image file: c9nr08209c-t4.tif(3)
where RG is the characteristic size of the inhomogeneities in the dilute concentration limit. While eqn (3) provides a reasonable estimate of void sizes from experimental data, difficulties arise in extracting the RG value from a ln[thin space (1/6-em)]Ik2 plot in simulations, owing to the finite size of the models. Since the effective small-angle region of k is approximately given by 4π/L < k < 1/RG, where L is the linear dimension of the sample, a sufficiently large model is needed in order to obtain the mean value of RG, via averaging over several independent configurations. However, it has been observed53 that, even for very large models, the lower limit of k cannot be reduced arbitrarily by increasing L due to the presence of noise in rg(r) beyond r ≈ 25–30 Å. The term r[g(r) − 1] in eqn (2) for r > 25 Å introduces artifacts in S(k) for small values of k, which make it difficult to accurately compute RG from the simulated intensity data. This necessitates the construction of alternative measures of void sizes and shapes directly from the distribution of atoms in the vicinity of voids. A simple but effective measure is to compute the gyrational radius from the distribution of the atoms on void surfaces between radii Rv and Rv + d (=Rv), where Rv is the void radius and d is the width of the void surface, as discussed in section IIB. This is schematically illustrated in two dimensions in Fig. 2, where the positions of the void-surface atoms before and after annealing are shown in yellow and red colors, respectively. Assuming that the void-surface atoms (yellow) can displace, in the mean-square sense, from their initial position by x during annealing, the radius of gyration, Rg, for the assembly of void atoms (red) can be expressed as,
image file: c9nr08209c-t5.tif(4)

image file: c9nr08209c-f2.tif
Fig. 2 An illustration of the (re)construction of the shape and size of a void formed by an assembly of void-surface atoms in two dimensions. The distributions of the atoms before and after annealing are shown in yellow (image file: c9nr08209c-u1.tif) and red (image file: c9nr08209c-u2.tif) colors, respectively. The convex polygon (red line) approximates the void region in two dimensions after annealing.

Here, ri and rcm are the position of the ith atom and the center of mass of ns void-surface atoms, respectively. The atoms on the reconstructed surface are distributed within a spherical region of radius Rv + x and H[y] is the Heaviside step function that asserts that only the atoms with y ≥ 0 contribute to Rg. Here, we chose a value of x = 1.5 Å for Si atoms in order to define the boundary of a reconstructed void surface.53 Since annealing at high temperature can introduce considerable restructuring of void surfaces, it is often convenient to invoke a suitable convex approximation to estimate the size and shape of the voids. To this end, we employ the convex-hull approximation that entails constructing the minimal convex polyhedron that includes all the void-surface atoms on the boundary. The size of the convex region can be expressed as the radius of gyration of the convex polyhedron, or convex hull,

image file: c9nr08209c-t6.tif(5)

In eqn (5), nH is the number of atoms (or vertices) that defines the convex polyhedral surface, and [r with combining macron]h is the center of mass of the polyhedron. To determine the degree of deviation of the shape of a (convex) void from an ideal sphere, we employ the convex-hull volume (VH) and surface area (AH), and obtain the sphericity parameter,54

image file: c9nr08209c-t7.tif(6)

The definition above is frequently used to measure the shape and roundness of sedimentary quartz particles, and it expresses the ratio of the surface area of a sphere to that of a non-spherical particle having an identical volume.54 It may be noted that a highly distorted void may not be accurately represented by a convex polyhedron and that the convex-hull volume provides an upper bound of the actual void volume, leading to RHRg. Likewise, as observed by Porod,55 the RG value obtained from the Guinier approximation may not be accurate for very anisometric voids.

III. Results and discussion

Before addressing the temperature-induced nanostructural changes in voids, we briefly examine the large model of a-Si, developed in section IIA, in order to validate its structural properties. Fig. 3 shows the simulated values of the X-ray scattering intensity, along with the experimental values of the intensity for as-implanted (green) and annealed (red) samples of pure a-Si, reported by Laaziri et al.56 The simulated values closely match with the experimental data in Fig. 3, especially for the annealed samples, in the wide-angle region from 1 Å−1 to 10 Å−1. Some characteristic structural properties of the large a-Si model are also listed in Table 1. Together with the X-ray intensity, the average bond angle and its deviation, 109.23° ± 9.26°, and the number of four-fold-coordinated atoms, 97.6%, suggest that the structural properties of the model obtained from classical molecular-dynamics simulations in section IIA are consistent with experimental results. It may be reasonably assumed that the presence of a small amount of coordination defects (≈2.4%), which involve a length scale of 2–3 Å and are sparsely distributed in the amorphous environment, would not affect the scattering intensity appreciably in the small-angle region of the scattering wavevector.57 A minor deviation in the height of the first peak in the simulated data in Fig. 3 from experimental values can be attributed in part to the use of the classical SW potential and partly to the large size of the model. A discussion on the validation of the structural and electronic properties of these large models can be found in ref. 34.
image file: c9nr08209c-f3.tif
Fig. 3 The X-ray scattering intensity for a-Si from experiments and classical molecular-dynamics simulations. The experimental data56 correspond to the results from as-implanted (image file: c9nr08209c-u3.tif) and annealed (image file: c9nr08209c-u4.tif) samples, while the simulated values (image file: c9nr08209c-u5.tif) of the intensity correspond to the modified Stillinger–Weber model of a-Si.
Table 1 Structural properties of the large model of a-Si. N, L, ρ, C4, r, θ, and Δθ are the number of atoms, simulation box length, mass density, four-fold coordination (in %), average bond length, average bond angle (in degree), and the root-mean-square width of the bond angles, respectively
N L (Å) ρ (g cm−3) C 4 (%) r (Å) θ Δθ
262[thin space (1/6-em)]400 176.123 2.24 97.6 2.39 109.23° 9.26°

A. Structural evolution of voids in a-Si at 300 K and 800 K

In this section, we discuss the three-dimensional shape and size of the voids in pure a-Si obtained from classical and quantum-mechanical annealing of the sub-systems at 300 K and 800 K, followed by total-energy optimization. Starting with an examination of the classical- and ab initio-generated void surfaces, we analyze the computed SAXS spectrum obtained from the large model, which is embedded with the supercells containing voids annealed at 300 K and 800 K.

Fig. 4 shows the structure of two representative voids (V3 and V16) after thermalization at 300 K for 10 ps, using classical and ab initio MD simulations. The void surfaces shown in Fig. 4 corresponds to Gaussian-convoluted isosurfaces from the Xcrysden package.58 The surface has been generated by placing a three-dimensional Gaussian function at the center of each atom and ensuring that the function has a value of 2.0 and 1.0 at the center and on the surface of the atom, respectively. By choosing an appropriate isovalue, between 0 to 2, and a suitable radius for Si atom, from standard crystallographic databases (e.g., Cambridge Structural Database), it is possible to examine the morphological changes of voids, associated with the structural relaxation of the void-surface atoms at different temperature. A close examination of Fig. 4(a) and (b) suggests that the resultant structures are very similar to each other as far as the void surfaces are concerned, except for small changes as indicated by white rings. A similar observation can be made for V16 in Fig. 4(c) and (d). This observation was found to be true for the rest of the voids as well, which showed very little changes on the void surfaces. These results are consistent with the variation of the simulated intensity values with wavevectors shown in Fig. 5. A somewhat more pronounced scattering from the classically-treated voids (CMD-R) in comparison to the ab initio-treated voids (AIMD-R) can be attributed to a minor expansion of the voids during classical simulations. This is clearly evident from Fig. 6(a) and (b), where the radius of gyration and the volume of the voids, respectively, are obtained from the convex-hull approximation of the void region. The polyhedral radii (RH) for the voids from classical simulations were found to be consistently larger than the corresponding radii from ab initio simulations in Fig. 6(a). A similar conclusion applies to polyhedral volumes (VH), which are plotted in Fig. 6(b). Once the linear size, area, and volume of the voids are available in the convex approximation, the sphericity of voids, ΦS, can be determined from eqn (6). A further measure of the size of the voids follows from the variation of the SAXS intensity, I(k), with the wavevector, k, in the small-angle region. By invoking the Guinier approximation (see eqn (3)), the Guinier radius, RG, can be obtained from ln I(k) vs. k2 plots. Fig. 6(c) shows the Guinier fits in the wavevector range, 0.15 Å ≤ k ≤ 0.3 Å, which resulted in an RG value of 6.38 Å for the CMD-R model and 8.07 Å for the AIMD-R model. It may be noted that a somewhat smaller value of RG was obtained for classically-treated voids compared to their AIMD counterparts (cf.Table 2), despite increased scattering by the former in the small-angle region (see Fig. 5). This minor discrepancy between the two RG values originated from the difficulty in computing the Guinier radius from the SAXS data in the wavevector region of 0.15–0.3 Å−1. For computational consistency, we employed an identical wavevector range, namely 0.15–0.3 Å−1, for the evaluation of RG. However, this range may not necessarily be the same for CMD-R and AIMD-R cases and can affect the RG value obtained from the Guinier's fit.

image file: c9nr08209c-f4.tif
Fig. 4 An isosurface representation of two voids (V3 and V16), obtained from the Gaussian-convolution of atomic positions, after thermalization at 300 K from classical and AIMD simulations: (a) classical V3; (b) AIMD V3; (c) classical V16; and (d) AIMD V16. The white rings in (a) and (b) indicate minor structural differences between the two surfaces in V3. See section IIIA for details.

image file: c9nr08209c-f5.tif
Fig. 5 The SAXS intensity of the large model of a-Si from classical (CMD-R) and ab initio (AIMD-R) MD simulations of voids at 300 K, followed by geometry relaxation.

image file: c9nr08209c-f6.tif
Fig. 6 (a) The radius of gyration of the voids from the convex-hull approximation at 300 K. The results for the CMD-R (image file: c9nr08209c-u6.tif) and AIMD-R (image file: c9nr08209c-u7.tif) models are presented here. (b) Estimated void volumes in the convex approximation. (c) Guinier plots for the CMD-R and AIMD-R models at 300 K. The Guinier radii from the plots correspond to a value of 6.38 Å (CMD-R) and 8.07 Å (AIMD-R).
Table 2 Estimated volumes (Å3), linear sizes (Å), and the average sphericities (ΦS) of the voids in a-Si with and without hydrogen inside voids. RX indicates the radius of the void obtained from the convex-hull (X = H) and Guinier approximations (X = G), and the gyrational radius of the void-surface atoms (X = g)
Models R g R G R H V H Φ S
a-Si (without hydrogen in voids)
CMD-R300 7.60 6.38 8.34 2100 0.66
AIMD-R300 7.55 8.07 8.28 2040 0.67
CMD-R800 7.56 6.57 8.31 2080 0.63
AIMD-R800 6.84 9.48 8.67 2250 0.68
a-Si (with hydrogen in voids)
AIMD-R300 6.66 8.66 8.36 2120 0.56
AIMD-R800 6.61 10.01 8.76 2290 0.71

Having studied the microstructure of voids at 300 K, we now examine the structure of voids at the high temperature of 800 K. The results for the SAXS intensity obtained from the large model at 800 K are presented in Fig. 7. It is evident from Fig. 7 that the high temperature annealing at 800 K via AIMD simulations produces notable changes in the SAXS spectrum by reducing intensity values in the small-angle region. This reduction in the SAXS intensity can be understood by examining the shape of the individual voids, obtained from AIMD-R models at 800 K. Fig. 8 shows the shape of a representative void (V3) in a-Si after annealing at 800 K from classical and ab initio MD simulations. A comparison between the two structures in Fig. 8(a) and (b) shows the regions on the two surfaces with a varying degree of reconstruction. The shape of the respective void surfaces obtained from using the convex-hull approximation is also shown in Fig. 8(c) and (d). Table 2 lists some of the parameters that characterize the void geometry in terms of the convex-hull radius (RH), the radius of gyration of the void-surface atoms (Rg), the Guinier radius (RG), the hull volume (VH), and the average sphericity of the voids (ΦS) after annealing at 300 K and 800 K. Although these scalar parameters do not vary much in the convex approximation, it is evident from the Gaussian isosurface representation of the void in Fig. 8(a) and (b), and the respective convex hulls in Fig. 8(c) and (d), that classical results noticeably differ from ab initio results, indicating the limitation of the classical force field in describing the restructuring of void surfaces at high temperature, namely at 800 K. This observation is also reflected in Fig. 9a–c, where the convex-hull radii, the hull volumes, and the intensity of SAXS at 800 K, respectively, are presented.

image file: c9nr08209c-f7.tif
Fig. 7 The SAXS intensity of the large model of a-Si from classical (CMD-R) and ab initio (AIMD-R) MD simulations of voids at 800 K, followed by geometry relaxation.

image file: c9nr08209c-f8.tif
Fig. 8 An isosurface representation of a void (V3) in a-Si obtained from: (a) classical; (b) ab initio simulations upon annealing at 800 K, followed by geometry relaxation. The corresponding shapes of the voids in the convex-hull approximation are shown in (c) and (d), respectively.

image file: c9nr08209c-f9.tif
Fig. 9 (a) The convex-hull radius for the voids from classical CMD-R (image file: c9nr08209c-u8.tif) and AIMD-R (image file: c9nr08209c-u9.tif) simulations at 800 K, followed by geometry relaxation. (b) The corresponding volume of voids from the convex-hull approximation. (c) Guinier plots for the models at 800 K. The Guinier radii from the plots correspond to a value of about 6.57 Å (CMD-R) and 9.48 Å (AIMD-R).

B. Structural evolution of hydrogen-rich voids in a-Si

Small-angle X-ray scattering measurements on hot-wire chemically vapor deposited (HWCVD) films indicate that the presence of preexisting H clusters can affect the nanostructure of voids due to H2 bubble pressure in the cavities, which, upon annealing, can lead to geometric changes, as observed via tilting SAXS and surface transmission electron microscopy (STEM).15 Although the degree of such nano-structural changes may depend on the growth mechanism and the rate of film growth, it is instructive to address the problem from a computational viewpoint using CRN models of a-Si with nanometer-size voids, which are characterized by an experimentally consistent void-volume fraction and clusters of H atoms inside the voids.

Fig. 10 shows the structure of a hydrogen-rich nanovoid surface (of V3) obtained after annealing at 300 K and 800 K, followed by total-energy minimization. The evolution of Si and H atoms during annealing was described using ab initio density-functional forces and energies, from the local-basis DFT package SIESTA.47 The procedure to conduct such AIMD simulations of voids, which are embedded in a large model of a-Si, has been described in section IIC. A comparison of Fig. 10(a) (with H atoms inside the void at 300 K) with Fig. 4(b) (with no H atoms at 300 K) reveals that the thermalization of the void (V3) and its neighboring region at 300 K does not introduce notable changes on the void surface, with the exception of the distribution of H atoms. The initial random distribution of H atoms, which was confined within a radius of 4 Å from the void center, evolved to produce several monohydrides (SiH) and a few dihydrides (SiH2), along with a few H2 molecules inside the cavity. The majority of H atoms continue to stay near the void at 300 K. By contrast, a considerable number of Si and H atoms have been found to be thermally driven out of the void region V3 at 800 K. An analysis of the distribution of Si atoms on the void surface, shown as yellow blob in Fig. 10, and H atoms (red) inside the void reveals that approximately 12.4% of the original Si atoms and 23.3% of the H atoms left the void region V3 of radius 10.5 Å after annealing at 800 K. This is apparent from the isosurface plot in Fig. 10(b), where a considerable number of Si atoms can be seen to move away from the void-surface region, creating a somewhat diffused or scattered isosurface. This observation is found to be true for most of the other voids. On average, over 12 independent voids, approximately 8.3% of the Si atoms and 16.1% of the H atoms were found to leave the void regions at 800 K. By contrast, the corresponding average values for the same at 300 K for Si and H atoms are 0% and 2.6%, respectively. It thus follows that at 800 K the restructuring of a void surface can be considerably affected via thermal and H-induced motion of Si atoms, as well as through the formation of various silicon-hydrogen bonding configurations on void surfaces. These changes can be also observed in the convex-hull approximation of the void surface (V3) in the presence and absence of H atoms at 800 K. The convex polyhedral region in Fig. 11b, with H atoms present inside the cavity, consists of 78 triangular faces and 41 vertices. These values are noticeably higher than the corresponding values of 62 faces and 32 vertices for the same cavity with no H atoms in Fig. 11a. The structural changes on the void surface have been reflected in the number of triangular faces, which are needed to construct the surface in the convex-hull approximation. Table 3 presents the statistics of various silicon-hydrogen bonding configurations in the vicinity of each void, which is defined by a radius of 10.5 Å.

image file: c9nr08209c-f10.tif
Fig. 10 The structure of a representative void (V3), in isosurface representation, from AIMD simulations in the presence of H atoms inside the cavity at: (a) 300 K; (b) 800 K. The yellow and red blobs represent Si and H atoms, respectively. For visualization and comparison on the same footing, an identical set of isosurface parameters was used to generate the surfaces.

image file: c9nr08209c-f11.tif
Fig. 11 The shape of a void (V3) obtained from the convex-hull approximation in AIMD-R simulations at 800 K in the presence and absence of H atoms in the void. (a) The convex surface with no H atoms. (b) The convex surface with 30 H atoms inside the void.
Table 3 Statistics of various silicon-hydrogen bonding configurations near hydrogen-rich voids at 300 K and 800 K from AIMD simulations. RH and ΦS indicate the convex-hull radius and sphericity of the voids, whereas asterisked entries in column 1 indicate at least one H atom left the void region of radial size 10.5 Å
Void Hiso H2 SiH SiH2 SiH3 R H Φ S
300 K
V1* 2 1 20 2 0 8.34 0.51
V2 0 3 22 1 0 8.33 0.56
V3 0 3 17 3 0 8.5 0.59
V4 0 0 26 1 0 8.4 0.58
V5 0 1 25 1 0 8.2 0.62
V6* 1 2 22 1 0 8.34 0.54
V7 0 1 26 1 0 8.31 0.54
V8 1 1 23 2 0 8.5 0.57
V9 0 3 20 0 1 8.35 0.51
V10* 1 3 18 2 0 8.28 0.60
V11 2 4 19 0 0 8.34 0.59
V12 2 2 23 0 0 8.52 0.56
800 K
V1 0 0 20 2 0 8.36 0.90
V2 0 2 20 2 0 8.84 0.66
V3* 1 2 16 1 0 9.07 0.79
V4 0 1 23 1 0 8.64 0.62
V5 0 0 20 2 0 8.95 0.70
V6 0 2 14 1 1 8.67 0.83
V7 0 3 14 3 0 8.94 0.64
V8* 1 1 14 3 0 8.95 0.73
V9 0 1 22 1 0 9.07 0.62
V10 0 1 21 2 0 8.71 0.62
V11 2 3 15 1 0 8.38 0.81
V12 1 0 21 2 0 8.56 0.56

To obtain the linear size of the voids, we computed the radius of gyration, Rg, using eqn (4). Likewise, the convex-hull approximation of a void region, defined by a set of void-surface atoms, provides the radius of the convex polyhedron, RH, from eqn (5). Fig. 12 shows the values of RH and the corresponding hull volumes, VH, of 12 voids, each of which has been computed upon annealing and geometry relaxation. Similarly, the Rg values obtained from the real-space distribution of the void atoms are listed in Table 2. The apparent deviation of RH from Rg can be attributed to the fact that the computation of Rg involves all the atoms on the surface and interior regions of the voids, whereas the convex-hull approximation includes only those atoms on the void-surface region that define the minimal convex polyhedral volume. Thus, RH provides an upper bound of the radial size of the voids. In order to compare RH and Rg values with the linear size of the voids from the simulated intensity plots in SAXS, we have invoked the Guinier approximation to estimate the Guinier radius, RG, from ln[thin space (1/6-em)]I(k)–k2 plots, as shown in Fig. 12(c). Table 2 suggests that the RG values obtained from the simulated intensity plots reasonably match with the corresponding average values of RH at 300 K and 800 K. It may be noted that, while the intensity plot in Fig. 12(c) is not particularly sensitive to the shape of the voids, the slightly pronounced scattering in the region k2 ≤ 0.05 Å−2 is possibly indicative of the expansion of the void volume at 800 K, as observed in the convex-hull volume in Fig. 12(b). The exact cause of this increased scattering intensity is difficult to determine but it appears that a combination of the pressure due to H2 bubbles, surface roughening of the voids, and the displacement of the atoms on the void surfaces at high temperature could lead to this additional scattering. This observation is corroborated by the experimental results from a small-angle X-ray scattering study of nanovoids in HWCVD amorphous Si films that indicated an increase of individual void volumes upon annealing at 813 K.15

image file: c9nr08209c-f12.tif
Fig. 12 (a) The radius of gyration of hydrogen-rich voids from AIMD simulations at 300 K (image file: c9nr08209c-u10.tif) and 800 K (image file: c9nr08209c-u11.tif) in the convex-hull approximation. (b) Estimated convex-hull volumes at 300 K and 800 K. (c) Guinier fits for the SAXS intensities in AIMD-R models of a-Si:H at 300 K and 800 K.

C. Hydrogen dynamics near voids in a-Si

In this section, we shall make a few observations on the dynamics of hydrogen atoms in the vicinity of void surfaces from ab initio molecular-dynamics simulations. Since the mass of H atoms is significantly smaller than that of Si atoms, the motion of H atoms is more susceptible to the annealing temperature and, more importantly, to the resulting temperature-induced structural changes on the void surfaces. This is particularly noticeable at high temperature. Having established earlier that the structural changes on the void surfaces are more pronounced at 800 K, we shall now examine the character of hydrogen dynamics inside (and near the surface of) the voids at 800 K. Despite the limited timescale of the present AIMD simulations, spanning 10 ps only, it is reasonable to assume that, at high temperature, the dynamical character of hydrogen motion inside the voids would be reflected in the first several picoseconds. It may be noted that, unlike structural properties, the dynamics of hydrogen in a-Si may depend on the accuracy of total energies and forces from DFT calculations. The use of the Harris-functional approximation in the present work, in an effort to address the structure of voids on the nanometer length scale, suggests that the results on hydrogen dynamics presented here are somewhat approximate in nature. A full self-consistent-field DFT calculation, over a period of a few tens of picoseconds, using plane-wave or even local basis sets for 2200-atom models is computationally overkill and outside scope of the present work on nanostructural evolution.

Fig. 13(a) depicts the time evolution of hydrogen atoms at 800 K within two voids, V12 and V7. Initially, all hydrogen atoms/molecules were uniformly distributed inside the voids within a radius of 4 Å from the center of the voids (see Fig. 1). The mean-square displacements (MSD), 〈R2(t)〉, of hydrogen atoms – averaged over all H atoms in the respective void – are shown in Fig. 13(a) from 0 to 10 ps. The corresponding real-space distributions of Si and H atoms on the void surfaces and their immediate vicinity are also shown in Fig. 13(b) and (c) for V12 and V7, respectively. The early-time behavior of the dynamics, from 0 to 0.5 ps, can be roughly understood from elementary considerations. Given that the majority of H atoms are away from the void surface at the beginning of simulation, by at least 2 Å or more, the early-time behavior of the dynamics would be primarily determined by the temperature of the system for a reasonably well-defined initial H distribution. Since the root-mean-square (RMS) speed of H atoms at 800 K for the Maxwell-Boltzmann distribution is about 0.045 Å fs−1, it would be reasonable to assume that the initial behavior can last for about a fraction of a picosecond, after which H atoms become affected by the presence of the void surface. It is the marked difference in the value of 〈R2(t)〉 in Fig. 13(a), from 1 ps to 8 ps, that we find most interesting. This notable difference between the two sets of MSD values beyond 0.5 ps can be understood by taking into account the real-space structure of the void surfaces. In addition to the trivial thermal and configurational fluctuations – caused by temperature and the disorder in H distributions, respectively – the motion of H atoms is also affected by the varying degree of re-structuring of Si atoms and the disorder of the void surfaces after 0.5 ps. The surface structure of V7 can be seen to be more open or scattered compared to its counterpart in V12, which is more compact or well-defined, as shown in Fig. 13(c) and (b), respectively. This has been also confirmed by analyzing the distributions of void-surface atoms in the interior of V12 and V7. The presence of Si atoms inside V7 may reduce the effective diffusion rate of H atoms, due to additional scattering from interior Si atoms, during intermediate-time evolution from 1 ps onward. By contrast, the more well-defined void surface of V12, with few Si atoms inside, provides less resistance to H atoms during their intermediate-time evolution. Thus, hydrogen diffusion in V12 proceeds relatively uninterruptedly, in the presence of few Si atoms from void surfaces, leading to a notable difference in their MSD values during the first several picoseconds of simulations. This observation has been found to be true for other pairs of voids as well, for which a notable difference in the MSD values exists. Additional results on the time evolution of H atoms in voids V5, V9 and V11, which exhibit a similar MSD behavior, are provided as ESI. In view of our observation that the intermediate-time behavior (i.e., from 1 ps to 7–8 ps) of the H dynamics can be influenced by the roughness of void surfaces, it might be necessary to take into account the appropriate size and the accurate surface structure of the voids in the calculations. However, a full SCF study of hydrogen dynamics inside nanometer-size voids in a-Si for several tens of picoseconds is computationally prohibitive and outside the scope of the present study.

image file: c9nr08209c-f13.tif
Fig. 13 (a) The time evolution of the mean-square displacements (MSD) of H inside two voids, V12 (image file: c9nr08209c-u12.tif) and V7 (image file: c9nr08209c-u13.tif), at 800 K. The difference between the two sets of MSD values, from 1 ps to 8 ps, can be attributed to the degree of void-surface restructuring, as discussed in section IIIC. (b) The compact or smooth structure of the V12 surface, obtained from an isosurface representation of the void, compared to a relatively diffused or scattered distribution of Si atoms in (c) V7 at 800 K.

D. Silicon-hydrogen bonding configurations on void surfaces

Experimental data from positron-annihilation lifetime (PAL) spectroscopy,13 Rutherford backscattering spectrometry (RBS),59 and Fourier-transform infrared spectroscopy-attenuated total reflection (FTIR-ATR)59 indicate that bonded and non-bonded hydrogens play an important role in characterizing the structural and optical properties of a-Si:H. Sekimoto et al.13,59 recently demonstrated that the presence of non-bonded hydrogens (NBH) at concentrations beyond 2.8 at% in amorphous silicon networks led to changes in the vacancy-size distribution and induced formation of nanometer-size voids in the network via relaxation of internal stress. The formation of H2 molecules and their evolution can be probed by an analysis of the low-temperature (LT) peak of the hydrogen-effusion profile13,60,61 near 750 K, whereas the number density of H2 molecules can be estimated from the collision-induced weak IR absorption in the frequency region of 4100 cm−1 to 5500 cm−1. Chabal and Patel24 reported a value of 1021 H2 per cm3 by analyzing data obtained from IR measurements at room temperature. Although the present study, based on pure a-Si networks, does not permit us to address the crucial role of silicon-hydrogen bonding configurations in a-Si:H in general, it does provide new insights into the role of bonded and non-bonded hydrogens in the vicinity of the walls of H-rich voids. Tables 3 and 4 list the statistics and the average concentration of various bonded and non-bonded hydride configurations inside the voids, respectively. Non-bonded hydrogens chiefly comprise H2 molecules, along with one or two isolated H atoms near the voids. The latter possibly originate from the short annealing time of 10 ps of AIMD simulations. The number of H2 molecules per void ranges from 0 to 4 and 0 to 3 at 300 K and 800 K, respectively, which translates into an average density of 1.5–2.2 × 1021 H2 per cm3 for nanometer-size voids. This value matches very well with the experimental value of 1021 H2 per cm3 from IR measurements by Chabal and Patel,24 as mentioned earlier. These results are also consistent with the reported data from RBS measurements.13 On the other hand, the great majority of bonded hydrogens appeared in monohydride (SiH) configurations. An examination of the local bonding environment of silicon atoms near the void surfaces showed that approximately 72.5 at% and 61.1 at% of total H atoms formed Si–H bonds at 300 K and 800 K, respectively (see Table 4). Of the remaining H atoms, about 8.3% (11.7%) were found to be bonded with Si as SiH2 configurations at 300 K (800 K). Apart from a few isolated H atoms, the rest of the hydrogen has left the spherical void region of radius 10.5 Å. A graphic distribution of the presence of H2 molecules and the bonded SiH/SiH2 configurations on the wall of V3 at 300 K is depicted in Fig. 14. Several SiH and a few SiH2 configurations can be seen to have formed on the surface of the void (V3). This observation is consistent with experimental results from IR measurements24 and the recent computational studies based on information-driven atomistic modeling of a-Si:H.17,21Fig. 15 shows the concentration of hydrogen associated with various silicon-hydrogen bonding configurations. The results clearly establish that a considerable number of H2 molecules, of about 9–13% of the total H, can form in the vicinity of voids, depending upon the annealing temperature, concentration of H atoms, the source of hydrogen, and the preparation methods of a-Si samples.
image file: c9nr08209c-f14.tif
Fig. 14 Bonded and non-bonded hydrogens inside a void (V3) from AIMD simulations at 300 K. (a) The interior wall of the void can be seen to be decorated with several SiH (yellow-red) and three SiH2 (blue-red) configurations, as well as three H2 molecules (white). (b) A magnified view, showing protruding monohydrides (yellow-red) and dihydrides (blue-red) on the walls, and a hydrogen molecule (white) inside the void.

image file: c9nr08209c-f15.tif
Fig. 15 Average statistics of silicon-hydrogen bonding configurations near void surfaces in a-Si:H after annealing at 300 K and 800 K. The results were obtained by averaging over 12 independent voids, distributed over a distance of several nanometers in the large a-Si model.
Table 4 The average concentration of hydrogen (in % of total H atoms) in silicon-hydrogen bonding configurations, H2 molecules and isolated H near voids (within a radius of 10.5 Å). Bonded and non-bonded hydrogen that left the voids during annealing are listed under the category Ex-void
Temp. Hiso H2 SiH SiH2 SiH3 Ex-void
300 K 2.5 13.3 72.5 8.3 0.8 2.6
800 K 1.4 8.9 61.1 11.7 0.8 16.1

We conclude this section with the following observation on silicon-hydrogen bond formation during AIMD simulations. The results from the preceding paragraph suggest that the internal surfaces of hydrogen-rich voids in a-Si are principally decorated with monohydride Si–H bonds. It is thus instructive to examine the time evolution of the number of Si–H bonds during annealing via MD simulations at low and high temperatures. Since a statistically significant number of H atoms participate in forming Si–H bonds on void surfaces, we restrict ourselves to address the formation of monohydride Si–H configurations only. Fig. 16 shows the evolution of the average hydrogen content in SiH bonds in the vicinity of voids at 300 K (blue) and 800 K (red). The results were obtained by sampling atomic configurations (of voids) at an interval of every 10 fs during the course of AIMD simulations, for a period of 7–8 ps, which were then averaged over 12 individual voids. It is apparent from Fig. 16 that the formation of monohydride Si–H bonds takes place very rapidly at 800 K, which can saturate the void surfaces within the first two picoseconds. By contrast, it takes several picoseconds for hydrogen atoms inside the cavities to passivate the void surfaces at 300 K. A comparison of the saturated values of the SiH bonds in Fig. 16, obtained from the MD-annealed models at 300 K and 800 K, with those from the corresponding final AIMD-R models in Table 4 suggests that the results are consistent with each other. A small difference between the saturated values of the number of SiH bonds at 300 K and 800 K (in Fig. 16) with those in Table 4 can be attributed to the total-energy relaxation of the final MD-annealed models, which were used in the calculation of silicon-hydrogen bonding statistics in Table 4.

image file: c9nr08209c-f16.tif
Fig. 16 The time evolution of monohydride SiH bonds during AIMD simulations at 300 K (blue) and 800 K (red). The vertical axis indicates the amount of hydrogen (in percent of total H atoms) that contribute to SiH bond formation. The results were obtained via configurational averaging over 12 voids.

IV. Conclusions

In this paper, we studied the temperature-induced nano-structural changes of the voids in a-Si at low and high annealing temperature in the range 300–800 K using classical and quantum-mechanical simulations, in the absence and presence of hydrogen near voids. This was achieved by generating a large model of a-Si, consisting of more than 262[thin space (1/6-em)]000 atoms, in order to be able to produce a realistic distribution of nanometer-size voids in the amorphous environment of Si atoms, as observed in SAXS, PAL, IR, and hydrogen-effusion measurements. An examination of the distribution of the atoms near the voids reveals that the reconstruction of void surfaces in pure a-Si at 300 K led to minimal changes in the shape and size of the voids, which can be readily understood from a classical treatment of the problem. By contrast, the high-temperature annealing at 800 K caused significant changes in the shape and size of the voids, with noticeable structural differences between classical and quantum-mechanical results. This observation appears to indicate that classical potentials, such as the modified SW potential, might not be particularly suitable for an accurate description of the dynamics of Si atoms near voids and the resulting reconstruction of the void surfaces at high temperatures.

An important outcome of the present study is that the dynamics of bonded and non-bonded hydrogens near the voids and the degree of void-surface reconstruction are found to be intrinsically related to each other. Ab initio annealing of a-Si networks with a void-volume fraction of 0.2% at 300 K and 800 K suggests that the presence of hydrogen within voids can facilitate surface reconstruction through the formation of silicon monohydride and dihydride configurations, as well as the displacement of Si atoms on the void surfaces. The rearrangement of atoms on void surfaces affects the diffusion of non-bonded hydrogens, which in turn produce surface bumps and changes the shape and size of the original void. This observation is reflected in the time evolution of hydrogen during annealing. The presence of Si atoms inside a heavily reconstructed rough void surface can reduce the effective diffusion rate of hydrogen, due to additional scattering from the interior Si atoms, in their intermediate-time evolution from 1 ps onward. In contrast, a compact void surface provides less resistance to hydrogen during their evolution. Thus, hydrogen diffusion inside a smooth or well-defined void proceeds with little or no hindrance, leading to a notable difference in their mean-square-displacement (MSD) values during their intermediate-time evolution. Finally, the results from the study show the presence of bonded hydrogens (BH), mostly SiH, SiH2, and non-bonded hydrogens (NBH) in the form of H2 molecules. The concentration of BH and NBH found in our study is consistent with the experimental values obtained from infrared (IR), Rutherford back scattering (RBS), and hydrogen-forward scattering (HFS) measurements. An examination of the silicon-hydrogen bonding configurations near the voids suggests that the interior walls of the voids are decorated with SiH2, which is supported by experimental results from infrared and ellipsometric studies. Likewise, the concentration of H2 molecules obtained from the first-principles density-functional simulations in the present study is found to be consistent with the experimental value estimated from the collision-induced weak IR absorption by H2 molecules in the frequency region of 4100–5500 cm−1. A preliminary study using more accurate DFT calculations on 1000-atom models, via the SCF solutions of the Kohn–Sham equation in the generalized-gradient approximation, with 20–40 H atoms inside the voids indicates that more SiH2 configurations tend to form at high concentration of hydrogen in the cavities. However, further studies are needed for an accurate determination of SiH, SiH2 and H2 molecules with varying concentrations of H atoms inside the voids.

Conflicts of interest

There are no conflicts of interest to declare.


This work was partially supported by the U.S. National Science Foundation under Grants No. DMR 1507166 and No. DMR 1507118. We acknowledge the Texas Advanced Computing Center at the University of Texas at Austin for providing HPC resources that have contributed to the results reported in this work. One of us (P. B.) thanks Prof. David Drabold (Ohio University) and Mr Dil Limbu (USM) for discussions.


  1. M. J. Thompson and H. C. Tuan, Amorphous Silicon Electronic Devices and their Applications, IEEE Aero. Electron. Sys. Mag., 1987, 2, 25–29 Search PubMed.
  2. R. A. Street, Technology and Applications of Amorphous Silicon, Springer-Verlag, 2000, vol. 37 Search PubMed.
  3. S. De Wolf, A. Descoeudres, Z. Holman and C. Ballif, High-efficiency silicon heterojunction solar cells: A review, Green, 2012, 2, 7–24 CAS.
  4. A. Kanevce and W. K. Metzger, The role of amorphous silicon and tunneling in heterojunction with intrinsic thin layer (HIT) solar cells, J. Appl. Phys., 2009, 105, 094507 CrossRef.
  5. M. Zeman and D. Zhang, Heterojunction Silicon Based Solar Cells, Springer Berlin Heidelberg, 2012, pp. 13–43 Search PubMed.
  6. J. Ge, Z. P. Ling, J. Wong, T. Mueller and A. Aberle, Optimisation of Intrinsic a-Si:H Passivation Layers in Crystalline-amorphous Silicon Heterojunction Solar Cells, Energy Procedia, 2012, 15, 107–117 CrossRef CAS.
  7. H. Fujiwara, M. Kondo and A. Matsuda, Real-time spectroscopic ellipsometry studies of the nucleation and grain growth processes in microcrystalline silicon thin films, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 115306 CrossRef.
  8. R. Collins, A. Ferlauto, G. Ferreira, C. Chen, J. Koh, R. Koval, Y. Lee, J. Pearce and C. Wronski, Evolution of microstructure and phase in amorphous, protocrystalline, and microcrystalline silicon studied by real time spectroscopic ellipsometry, Sol. Energy Mater. Sol. Cells, 2003, 78, 143–180 CrossRef CAS.
  9. H. Daxing, J. D. Lorentzen, J. Weinberg-Wolf, L. E. McNeil and Q. Wang, Raman study of thin films of amorphous-to-microcrystalline silicon prepared by hot-wire chemical vapor deposition, J. Appl. Phys., 2003, 94, 2930–2936 CrossRef.
  10. B. Macco, J. Melskens, N. Podraza, K. Arts, C. Pugh, O. Thomas and W. Kessels, Correlating the silicon surface passivation to the nanostructure of low-temperature a-Si:H after rapid thermal annealing, J. Appl. Phys., 2017, 122, 035302 CrossRef.
  11. D. L. Williamson, S. J. Jones and Y. Chen, Small-Angle x-ray Scattering Studies of Microvoids in Amorphous-Silicon-Based Semiconductors, 1994, p. 6395, NREL/TP-451-6395.UC Category 271. DE94006931 Search PubMed.
  12. M. Serényi, C. Frigeri, Z. Szekrényes, K. Kamarás, L. Nasi, A. Csik and N. Q. Khánh, On the formation of blisters in annealed hydrogenated a-Si layers, Nanoscale Res. Lett., 2013, 8, 84 CrossRef.
  13. T. Sekimoto, M. Matsumoto, A. Sagara, M. Hishida and A. Terakawa, Changes in the vacancy size distribution induced by non-bonded hydrogens in hydrogenated amorphous silicon, J. Non-Cryst. Solids, 2016, 447, 207–211 CrossRef CAS.
  14. D. T. Britton, A. Hempel, M. Härting, G. Kögel, P. Sperr, W. Triftshäuser, C. Arendse and D. Knoesen, Annealing and recrystallization of hydrogenated amorphous silicon, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 075403 CrossRef.
  15. D. L. Young, P. Stradins, Y. Xu, L. M. Gedvilas, E. Iwaniczko, Y. Yan, H. M. Branz, Q. Wang and D. L. Williamson, Nanostructure evolution in hydrogenated amorphous silicon during hydrogen effusion and crystallization, Appl. Phys. Lett., 2007, 90, 081923 CrossRef.
  16. D. L. Williamson, A. H. Mahan, B. P. Nelson and R. S. Crandall, The observation of microvoids in device quality hydrogenated amorphous silicon, J. Non-Cryst. Solids, 1989, 114, 226–228 CrossRef CAS.
  17. P. Biswas and S. R. Elliott, Nanoscale structure of microvoids in a-Si:H: a first-principles study, J. Phys.: Condens. Matter, 2015, 27, 435201 CrossRef.
  18. P. Biswas, D. Paudel, R. Atta-Fynn, D. A. Drabold and S. R. Elliott, Morphology and Number Density of Voids in Hydrogenated Amorphous Silicon: An Ab Initio Study, Phys. Rev. Appl., 2017, 7, 024013 CrossRef.
  19. R. Biswas, I. Kwon, A. M. Bouchard, C. M. Soukoulis and G. S. Grest, Intense small wave-vector scattering from voids in amorphous silicon: A theoretical simulation, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 39, 5101–5106 CrossRef CAS.
  20. R. Brahim and A. Chehaidar, Small-Angle X-Ray Scattering of Amorphous Germanium: Numerical Modeling, Adv. Mater. Phys. Chem., 2013, 3, 19–30 CrossRef.
  21. P. Biswas, D. A. Drabold and R. Atta-Fynn, Microstructure from joint analysis of experimental data and ab initio interactions: Hydrogenated amorphous silicon, J. Appl. Phys., 2014, 116, 244305 CrossRef.
  22. S. Muramatsu, S. Matsubara, T. Watanabe, T. Shimada, T. Kamiyama and K. Suzuki, Small-angle X-ray scattering studies of a-Si1−xGex: H alloys, J. Non-Cryst. Solids, 1992, 150, 163–166 CrossRef CAS.
  23. J. B. Boyce and M. Stutzmann, Orientational Ordering and Melting of Molecular H2 in an a-Si Matrix: NMR Studies, Phys. Rev. Lett., 1985, 54, 562–565 CrossRef CAS PubMed.
  24. Y. J. Chabal and C. K. N. Patel, Infrared Absorption in $a$-Si: H: First Observation of Gaseous Molecular ${\mathrm{H}}_{2}$ and Si-H Overtone, Phys. Rev. Lett., 1984, 53, 210–213 CrossRef CAS.
  25. A. von Keudell and J. R. Abelson, Evidence for atomic H insertion into strained Si-Si bonds in the amorphous hydrogenated silicon subsurface from in situ infrared spectroscopy, Appl. Phys. Lett., 1997, 71, 3832–3834 CrossRef CAS.
  26. Y. J. He, M. Hasegawa, R. Lee, S. Berko, D. Adler and A.-L. Jung, Positron-annihilation study of voids in a-Si and a-Si:H, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 5924–5927 CrossRef CAS PubMed.
  27. J. E. Graebner, B. Golding, L. C. Allen, D. K. Biegelsen and M. Stutzmann, Solid Hydrogen in Hydrogenated Amorphous Silicon, Phys. Rev. Lett., 1984, 52, 553–556 CrossRef CAS.
  28. J. Melskens, S. W. H. Eijt, M. Schouten, A. S. Vullers, A. Mannheim, H. Schut, B. Macco, M. Zeman and A. H. M. Smets, Migration of Open Volume Deficiencies in Hydrogenated Amorphous Silicon During Annealing, IEEE J. Photovoltaics, 2017, 7, 421–429 Search PubMed.
  29. R. L. C. Vink, G. T. Barkema, W. F. van der Weg and N. Mousseau, Fitting the Stillinger-Weber potential to amorphous silicon, J. Non-Cryst. Solids, 2001, 282, 248–255 CrossRef CAS.
  30. F. H. Stillinger and T. A. Weber, Computer simulation of local order in condensed phases of silicon, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 5262–5271 CrossRef CAS.
  31. J. S. Custer, M. O. Thompson, D. C. Jacobson, J. M. Poate, S. Roorda, W. C. Sinke and F. Spaepen, Density of amorphous Si, Appl. Phys. Lett., 1994, 64, 437–439 CrossRef CAS.
  32. S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys., 1984, 81, 511–519 CrossRef.
  33. W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A, 1985, 31, 1695–1697 CrossRef.
  34. R. Atta-Fynn and P. Biswas, Nearly defect-free dynamical models of disordered solids: The case of amorphous silicon, J. Chem. Phys., 2018, 148, 204503 CrossRef.
  35. P. D'Antonio and J. H. Konnert, Small-Angle-Scattering Evidence of Voids in Hydrogenated Amorphous Silicon, Phys. Rev. Lett., 1979, 43, 1161–1163 CrossRef.
  36. A. C. Wright, A. C. Hannon, R. N. Sinclair, T. M. Brunier, C. A. Guy, R. J. Stewart, M. B. Strobel and F. Jansen, Neutron scattering studies of hydrogenated, deuterated and fluorinated amorphous silicon, J. Phys.: Condens. Matter, 2007, 19, 415109 CrossRef.
  37. M. S. Conradi and R. E. Norberg, Molecular h2: Nuclear-spin-relaxation centers for protons in a-Si: H, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 24, 2285–2288 CrossRef CAS.
  38. H. V. Löhneysen, H. J. Schink and W. Beyer, Direct Experimental Evidence for Molecular Hydrogen in Amorphous Si: H, Phys. Rev. Lett., 1984, 52, 549–552 CrossRef.
  39. W. Kohn, Density Functional and Density Matrix Method Scaling Linearly with the Number of Atoms, Phys. Rev. Lett., 1996, 76, 3168–3171 CrossRef CAS.
  40. H. Ehrenreich, F. Seitz and D. Turnbull, Solid State Physics, Elsevier Science, 1980 Search PubMed.
  41. R. Haydock, V. Heine and M. Kelly, Electronic structure based on the local atomic environment for tight-binding bands. II, J. Phys. C: Solid State Phys., 1975, 8, 2591–2605 CrossRef.
  42. Strictly speaking, the boundary must be located at a distance greater than the characteristic decay length of the density matrix of the system. The presence of an electronic gap in a-Si asserts that its density matrix 62 can be characterized by a well-defined decay length.
  43. The rationale behind adding 30 H atoms inside the voids is based on calorimetry and infrared measurements, 24 which suggest the presence of about 12 to 15 “bulk” H2 molecules inside nanometer-size voids.
  44. J. Harris, Simplified method for calculating the energy of weakly interacting fragments, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 1770–1779 CrossRef CAS.
  45. E. Zaremba, Extremal properties of the Harris energy functional, J. Phys.: Condens. Matter, 1990, 2, 2479–2486 CrossRef.
  46. R. Atta-Fynn, P. Biswas, P. Ordejón and D. A. Drabold, Systematic study of electron localization in an amorphous semiconductor, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 085207 CrossRef.
  47. J. M. Soler, E. Artacho, J. D. Gale, A. Garciá, J. Junquera, P. Ordejón and D. Sánchez-Portal, The SIESTA method for ab initio order- N materials simulation, J. Phys.: Condens. Matter, 2002, 14, 2745 CrossRef CAS.
  48. N. Troullier and J. L. Martins, Efficient pseudopotentials for plane-wave calculations, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 43, 1993–2006 CrossRef CAS.
  49. L. Kleinman and D. M. Bylander, Efficacious Form for Model Pseudopotentials, Phys. Rev. Lett., 1982, 48, 1425–1428 CrossRef CAS.
  50. J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 5048–5079 CrossRef CAS.
  51. N. E. Cusack, The Physics of Structurally Disordered Matter, Adam Hilger, 1987 Search PubMed.
  52. A. Guinier, G. Fournet and K. L. Yudowitch, Small-angle scattering of X-rays, Wiley, New York, 1995 Search PubMed.
  53. D. Paudel, R. Atta-Fynn, D. A. Drabold, S. R. Elliott and P. Biswas, Small-angle X-ray scattering in amorphous silicon: A computational study, Phys. Rev. B, 2018, 97, 184202 CrossRef CAS.
  54. H. Wadell, Volume, Shape, and Roundness of Quartz Particles, J. Geol., 1935, 43, 250–280 CrossRef.
  55. O. G. Porod, Small Angle X-ray Scattering, ed. O. Glatter and O. Kratky, Academic Press, New York, 1982, pp. 17–51 Search PubMed.
  56. K. Laaziri, S. Kycia, S. Roorda, M. Chicoine, J. L. Robertson, J. Wang and S. C. Moss, High-energy x-ray diffraction study of pure amorphous silicon, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 60, 13520–13533 CrossRef CAS.
  57. The observation rests on the assumptions that the variation of the local mass density near the defect sites and the corresponding length scale of 2–3 Å, associated with a sparse distribution of a few isolated defects, would not be too large to induce additional scattering in the small-angle region of the scattering wavevector. Suffice it to mention that, for clustered defects, this may not be necessarily true. These assumptions have been verified computationally by creating artificial defects in 100% 4-fold-coordinated a-Si networks here.
  58. A. Kokalj, XCrySDen¿a new program for displaying crystalline structures and electron densities, J. Mol. Graphics Modell., 1999, 17, 176–179 CrossRef CAS.
  59. T. Sekimoto, M. Matsumoto and A. Terakawa, Dense restructuring of amorphous silicon network induced by non-bonded hydrogen, Jpn. J. Appl. Phys., 2018, 57, 08RB07 CrossRef.
  60. W. Beyer, Hydrogen effusion: a probe for surface desorption and diffusion, Phys. B, 1991, 170, 105–114 CrossRef CAS.
  61. W. Beyer, Diffusion and evolution of hydrogen in hydrogenated amorphous and microcrystalline silicon, Sol. Energy Mater. Sol. Cells, 2003, 78, 235–267 CrossRef CAS.
  62. S. Goedecker, Linear scaling electronic structure methods, Rev. Mod. Phys., 1999, 71, 1085–1123 CrossRef CAS.


Electronic supplementary information (ESI) available. See DOI: 10.1039/C9NR08209C

This journal is © The Royal Society of Chemistry 2020