Atomistic simulation of nearly defect-free models of amorphous silicon: An information-based approach

TL;DR

This paper introduces an information-based total-energy optimization method to generate nearly defect-free atomistic models of amorphous silicon, outperforming traditional simulation techniques and matching the quality of established models.

Contribution

The paper presents a novel hybrid data-driven approach that combines geometrical, structural, and topological information to produce superior amorphous silicon models.

Findings

Models are comparable to modified WWW models in structural and electronic properties.

The hybrid approach outperforms conventional reverse Monte Carlo and molecular dynamics methods.

Structural, electronic, and thermodynamic properties are validated against recent machine-learning-based models.

Abstract

We present an information-based total-energy optimization method to produce nearly defect-free structural models of amorphous silicon. Using geometrical, structural and topological information from disordered tetrahedral networks, we have shown that it is possible to generate structural configurations of amorphous silicon, which are superior than the models obtained from conventional reverse Monte Carlo and molecular-dynamics simulations. The new data-driven hybrid approach presented here is capable of producing atomistic models with structural and electronic properties which are on a par with those obtained from the modified Wooten-Winer-Weaire (WWW) models of amorphous silicon. Structural, electronic and thermodynamic properties of the hybrid models are compared with the best dynamical models obtained from using machine-intelligence-based potentials and efficient classical molecular-dynamics simulations, reported in the recent literature. We have shown that, together with the WWW models, our hybrid models represent one of the best structural models so far produced by total-energy-based Monte Carlo methods in conjunction with experimental diffraction data and a few structural constraints.

Tables

Table 1. Table 1: Structural properties of INDIA, SW-MD and WWW models: average bond length < r 𝑟 r >, average bond angle < θ 𝜃 \theta >, the root-mean-square (RMS) deviation of the bond angles Δ ​ θ Δ 𝜃 \Delta\theta , the percentage of 4-fold coordination number C 4 and the effective coordination number (ECN), respectively.

Model	N	<$r$>	$⟨\theta ⟩$	$\Delta \theta$	$\left(\Delta {\theta }_G\right)$111Values obtained from the Gaussian approximation of the bond-angle distribution.	C${}_4{}^{}(\%)$	ECN
INDIA	300	2.386	109.10	11.41	10.42	99.33	3.956
SW-MD		2.380	109.22	9.31	8.44	99.33	3.973
WWW		2.378	109.18	10.44	10.13	100.0	3.963
INDIA	512	2.387	109.10	11.48	10.64	99.60	3.955
SW-MD		2.379	109.27	9.12	8.58	99.22	3.976
ML-MD222From Ref. Deringer et al.,2018		2.371	109.19	9.69	9.36	98.44	3.940
WWW		2.365	109.11	10.69	10.47	100.0	3.974
INDIA	1024	2.390	109.01	11.96	10.78	98.34	3.958
SW-MD		2.381	109.27	8.94	8.46	99.22	3.968
WWW		2.371	109.14	10.63	10.30	100.0	3.969

Equations

$$V({\mathbf{R}}^N)=\frac{1}{2}\sum _{i=1}^N\sum _{\begin{array}{c}j=1\\ (j\ne i)\end{array}}^N{v}_2(r_{ij})+\sum _{i=1}^N\sum _{\begin{array}{c}j=1\\ (j\ne i)\end{array}}^N\sum _{\begin{array}{c}k=1\\ (k\ne i)\\ (k>j)\end{array}}^N{v}_3({\mathbf{r}}_{ij},{\mathbf{r}}_{ik}).$$

$$v_2(r_{ij})=ϵA\left[B{\left(\frac{r_{ij}}{\sigma }\right)}^{-p}-1\right]\exp \left(\frac{\sigma }{r_{ij}-a\sigma }\right)\Theta (a\sigma -r_{ij}),$$

$$\begin{array}{rl}{v}_3({\mathbf{r}}_{ij},{\mathbf{r}}_{ik})& =ϵ\lambda {\left(\cos {\theta }_{jik}+\frac{1}{3}\right)}^2\exp \left(\frac{\sigma \gamma }{r_{ij}-a\sigma }+\frac{\sigma \gamma }{r_{ik}-a\sigma }\right)\\ & \times \Theta (a\sigma -r_{ij})\Theta (a\sigma -r_{ik}),\end{array}$$

$$C_v(T)=k_B\sum _{\mathbf{k},j}\frac{(\frac{\hslash {\omega }_j(\mathbf{k})}{2k_{B}T}{)}^2}{{\sinh }^2(\frac{\hslash {\omega }_j(\mathbf{k})}{2k_{B}T})},$$

Full text

I Introduction

Amorphous silicon continues to play a major role in the application of silicon-based device technology. Street (2000) Recent developments of silicon-based heterojunction intrinsic technology (HIT) Taguchi et al. (2014) for photovoltaic cells and the two-qubit quantum logic gates Veldhorst et al. (2015) are indicative of the continuing importance of amorphous/crystalline silicon and silicon-based materials. While the structure of amorphous silicon (a-Si) can be readily described by using a continuous random network (CRN) model,Zachariasen (1932) which is characterized by the presence of high degree of tetrahedral ordering, the construction of CRN models with minimal strain and few coordination defects has been a vexing problem in structural modeling of amorphous silicon. Until recently, high-quality structural models of a-Si are best produced by the Monte Carlo based bond-switching algorithm of Wooten, Winer and Weaire (WWW). Wooten et al. (1985) The subsequent modification of the algorithm by Barkema and Mousseau Barkema and Mousseau (2000) further augments the capability of the WWW method in producing large structural models of a-Si, exhibiting experimentally compliant structural, electronic and vibrational properties. Other important approaches that are often employed to simulate a-Si are molecular-dynamics (MD) simulations Deringer et al. (2018); Atta-Fynn and Biswas (2018); Pedersen et al. (2017); Car and Parrinello (1988) (either ab initio/classical MD or an intermediate approach between the two), reverse Monte Carlo simulations (RMC), Biswas et al. (2004); McGreevy (2001); Gereben and Pusztai (1994) and the recently developed hybrid RMC methods in various flavors Cliffe et al. (2010); Pandey et al. (2015); Biswas et al. (2005) including the Force Enhanced Atomic Refinement (FEAR) approach. Pandey et al. (2015) In the following, we present briefly an information-driven inverse approach (INDIA), Limbu et al. (2018) where we have shown that high-quality structural models of amorphous silicon can be produced efficiently by inverting experimental scattering data, such as structure-factor or pair-correlation data, along with additional geometrical and topological information on the local bonding environment of silicon atoms in the amorphous state.

The rest of the paper is as follows. In Sect. II, we briefly outlined the computational method for structural modeling of amorphous silicon using classical and quantum-mechanical simulations. Section III addressed the results and discussion with an emphasis on the structural, electronic, vibrational and the thermodynamic properties of the models. This is followed by the conclusions of our work.

II Computational Method

Recently, we have implemented an information-driven inverse approach (INDIA) to model the atomistic structure of amorphous silicon. A detailed description of the method can be found in Ref. Limbu et al., 2018. Here, we primarily discussed the characteristic features of the resulting models with particular emphasis on the structural, topological, vibrational and thermodynamic properties of the models. Hereafter, we refer to these models as INDIA models. Following Ref. Limbu et al., 2018, it can be shown that high-quality structural models of a-Si can be generated by inverting experimental diffraction data and structural constraints in conjunction with an appropriate total-energy functional. For amorphous silicon, the method can be illustrated by employing the modified Stillinger-Weber (SW) potential, Stillinger and Weber (1985); Vink et al. (2001) which is given by,

[TABLE]

In Eq. 1, the two-body contribution, $v_{2}(r_{ij})$, is given by,

[TABLE]

and the three-body contribution follows from,

[TABLE]

with $\Theta$ being the Heaviside step function. The potential parameters, due to Vink *et al. *, Vink et al. (2001) were used in this work. The simulations were carried out by maintaining a minimum distance of 2.0 Å between silicon atoms and the mass density of the model was matched with the experimental density of 2.25 g/cm3 for amorphous silicon. The electronic and vibrational properties of models were studied using the local-basis density-functional theory (DFT) code Siesta. Soler et al. (2002) The latter employs pseudoatomic orbitals as basis functions (employing the double-zeta basis functions in the present work) and the norm-conserving Troullier-Martins pseudopotentials Troullier and Martins (1991) within the Perdew-Burke-Ernzerhof (PBE) formulation Perdew et al. (1996) of the generalized gradient approximation (GGA).

III Results and Discussion

In this section, we discussed the results from the INDIA models and compared with the corresponding WWW models. In addition, we also examined the results from molecular-dynamical models obtained using machine-intelligence-based algorithms (ML-MD) Deringer et al. (2018) and molecular-dynamics simulations (SW-MD). Atta-Fynn and Biswas (2018) We studied the structural properties of a-Si models consisting of up to 1024 atoms and presented the results in Table 1, which listed various structural properties of the relaxation-based INDIA models, WWW models and MD-based SW models showing the average bond length (<$r$>), average bond angle (<$\theta$>) and the percentage of 4-fold coordination number (C4), as well as the effective coordination number (ECN). The last two quantities were computed by using the relations described in Refs. Hoppe, 1979; Limbu et al., .

The static structure factor, S($k$), of the INDIA models was compared with the WWW models and the experimental structural-factor data from Laaziri *et al. * Laaziri et al. (1999) The structure factors presented in Fig. 1(a) show that the results from the 512-atom INDIA model agree with the WWW model of identical size, as far as the two-body correlations are concerned.

The results also match with the experimental data for as-deposited a-Si samples from Laaziri *et al. * Laaziri et al. (1999) To examine the interstitial volumes associated with the model networks, we computed the Voronoi volumes (of the atoms) of the 512-atom INDIA, WWW and SW-MD models and plotted the respective Voronoi-volume distributions in Fig. 1(b). The following observations are now in order. First, the Voronoi-volume distribution of a model can provide additional information on the degree of local ordering of a disordered network in a manner similar to the pair-correlation and bond-angle distributions. A narrow distribution with a sharp peak is indicative of more ordered networks, whereas networks with less order tend to show a relatively broad distribution. Second, since the Voronoi volume associated with an atom depends on the spatial positions of its nearest neighbors, the degree of inhomogeneity of a disordered network can be approximately gauged by analyzing the volume distribution, particularly the region away from the central peak. A large (or small) Voronoi volume is indicative of the presence of sparse (or dense) atomic environment in the network. not Thus, the presence of microvoids and the density-deficient regions in the network can be readily manifested in the Voronoi-volume distribution. It is apparent from Fig. 1(b) that, as far as the Voronoi volumes or the interstitial regions of the atoms are concerned, the SW-MD models exhibit more order than the corresponding WWW and INDIA models. This observation is quite consistent with the fact that MD models generally represent annealed samples of a-Si more closely, and hence more ordered, than those obtained from the relaxation-based Monte Carlo or similar approaches.

Since the bond angles between the nearest-neighbor atoms provide limited three-body correlations between atoms, we analyze the bond-angle distribution (BAD) of the atoms in the networks. The distribution is found to be essentially Gaussian in nature, which gives the average bond angle close to the value of ideally tetrahedrally bonded silicon atoms of 109.47∘ in crystalline silicon and the root-mean-square (RMS) deviation of the bond angles ($\Delta\theta$) in the experimentally observed range of 9-11∘. Beeman et al. (1985) To obtain $\Delta\theta$, we fitted the BAD with a Gaussian distribution. Although the Gaussian approximation to the bond-angle distribution somewhat underestimates the RMS value, due to the absence of a few large/small angles in the distribution, it describes the BAD fairly accurately for high-quality a-Si networks. This is reflected in Table 1.

Further characterization of the network is possible by examining the bond-orientational order parameter (BOP), $Q_{l}$, as defined by Steinhardt *et al. *, Steinhardt et al. (1983) which provides information on the orientation of a group of bonds. Figure 2 shows the $Q_{l}$ values, for $l$ = 1 to 8, for INDIA, WWW and SW-MD models consisting of 512 atoms. The BOP not only incorporates some aspects of structural information from higher-order correlation functions, but also provides a simple and effective measure for determining the presence of microcrystalline or paracrystalline structural units in the networks. Limbu et al. (2017) Here, we have used $Q_{l}$ ($l$=1,2,5) to determine the degree of crystallinity. While the magnitude of $Q_{1}$, $Q_{2}$ and $Q_{5}$ are exactly zero for ideal c-Si Atta-Fynn and Biswas (2018) networks, a high value of $Q_{5}$ for the INDIA and WWW models indicates more amorphous or disordered nature of these models in comparison to the SW-MD model. This observation is in agreement with the conclusion followed from the Voronoi-volume distribution as shown in Fig. 1(b).

We also examined the topological connectivity and the associated length scale(s) for various irreducible ring structures present in the amorphous networks by computing the radius of gyration ($R_{g}$) of the rings. The $R_{g}$ values and the numbers of the high-member rings ($\geq$ 7-member) can be indicative of the presence of intermediate range order (IRO) in amorphous networks. The distributions of $R_{g}$ of 6-member rings and 7-member rings are shown in Fig. 3(a) and Fig. 3(b) for the 512-atom models, respectively. The plot in Fig. 3(a) indicates the presence of more 6-member rings in the SW-MD model compared to the INDIA and WWW models of an identical size. This can be attributed to the somewhat more ordered nature of the SW-MD model, which has been obtained from the MD simulations and thus can be compared with annealed samples of a-Si. The average $R_{g}$ value of 6-member rings in c-Si is found to be 2.25 Å, which is somewhat lower that the average bond length (of 2.35 Å) in ideal c-Si, due to the non-planar nature of the hexagonal rings. Similarly, Fig. 3(b) shows the results for 7-member rings that are somewhat larger than the corresponding 6-member rings.

Finally, we addressed the electronic and thermodynamic properties of the INDIA models. Toward this end, we computed the electronic density of states and the variation of the specific heat at constant volume with temperature. Figure 4(a) shows the density of electronic states for the 300- and 512-atom INDIA models and compared the results with the same from the 512-atom WWW model. The INDIA models produced a remarkably clean electronic gap, almost identical to the WWW model. We studied the specific heat at constant volume ($C_{v}$) of the 300- and 512-atom INDIA models in the temperature range of 10 K–300 K. The specific heat was calculated from the (discrete) vibrational frequencies ($\omega$) obtained from direct diagonalization of the dynamical matrices in the harmonic approximation using the relation, Maradudin et al. (1971)

[TABLE]

and compared the resulting values with the experimental data obtained by Zink *et al. * Zink et al. (2006) at low temperature in the range from 10 K to 300 K. As mentioned in Sect. II, the dynamical matrices of the models were obtained from employing the density-functional code Siesta. Given the temperature range studied here, the harmonic approximation was found to be adequate in the present work. In Eq. 4, the sum over $\bm{k}$ corresponds to the $\Gamma$ point only. Figure 4(b) shows the low-temperature dependence of the specific heat ($C_{v}/T^{3}$) with temperature from 10 K to 300 K. The results suggest that the INDIA models are in good agreement with the experimental data for $T>$ 40 K. The inset shows the classical Dulong-Petit limit, $C_{v}\to 3R$, as the temperature approaches to 300 K.

IV Conclusions

In this paper, we have studied the structural, electronic, vibrational and thermodynamic properties of amorphous silicon obtained from using an information-driven inverse approach developed in Ref.Limbu et al., 2018, which simultaneously uses experimental structure factor and a few structural constraints along with a classical total-energy functional. By introducing a subspace optimization scheme, it was shown that structural constraints can be readily incorporated to produce high-quality model configurations of a-Si, which are otherwise very difficult to achieve using the conventional RMC and hybrid RMC schemes. The resulting data-driven optimization method produces high-quality (i.e., a few coordination defects and a narrow bond-angle distribution) CRN models of amorphous silicon of size up to 1024 atoms. The resulting relaxation-based models have the average bond angle of 109.2∘ and the root-mean-square (RMS) deviation of 9-11∘, which are within the range of experimental values obtained from Raman measurements. The data-driven computational approach can produce realistic models of a-Si, which exhibit accurate structural properties as produced by the WWW method and observed in experiments. The INDIA models show not only an excellent agreement with structural properties but also able to produce a pristine gap in the electronic density of states and the thermodynamic and vibrational properties, compatible with experimental measurements.

Acknowledgments

The work is partially supported by the US National Science Foundation (NSF) under Grants No. DMR 1507166 and No. DMR 1507118. PB thanks Prof. Stephen Elliott (Cambridge, UK) for discussions.