Nearly defect-free dynamical models of disordered solids: The case of amorphous silicon
aa r X i v : . [ c ond - m a t . d i s - nn ] M a r Nearly defect-free dynamical models of disordered solids: The case ofamorphous silicon
Raymond Atta-Fynn a) and Parthapratim Biswas b) Department of Physics, University of Texas at Arlington, Arlington, TX 76019 Department of Physics and Astronomy, The University of Southern Mississippi, Hattiesburg,MS 39406
It is widely accepted in the materials modeling community that defect-free realistic networks of amorphous siliconcannot be prepared by quenching from a molten state of silicon using classical or ab initio molecular-dynamics (MD)simulations. In this work, we address this long-standing problem by producing nearly defect-free ultra-large models ofamorphous silicon, consisting of up to half-a-million atoms, using classical molecular-dynamics simulations. The struc-tural, topological, electronic, and vibrational properties of the models are presented and compared with experimentaldata. A comparison of the models with those obtained from using the modified Wooten-Winer-Weaire bond-switchingalgorithm shows that the models are on par with the latter, which were generated via event-based total-energy relax-ations of atomistic networks in the configuration space. The MD models produced in this work represent the highestquality of amorphous-silicon networks so far reported in the literature using molecular-dynamics simulations.
I. INTRODUCTION
Silicon continues to play a major role in the technologicalrevolution of the 21st century. The recent developments ofhigh-efficiency silicon-based heterojunction solar cells andtwo-qubit quantum logic gates in silicon, with applicationsto quantum computation, are indicative of its lasting impor-tance for future technology. Unlike the crystalline state, theamorphous state of silicon is characterized by the presenceof disorder in the radial and bond-angle distributions, alongwith a distinctly different topological (dis)ordering from itscrystalline counterpart. Although a number of models basedon the presence of small clusters of ordered materials – var-iously known as the paracrystalline model, the significantstructure theory, and the crystalline hypothesis – have beenproposed from time to time, a CRN model appears to be con-ceptually simple and representative of most of the characteris-tic properties of a -Si and a -Si:H, observed in experiments. Itis now widely accepted that the structure of amorphous siliconcan be described fairly accurately by the continuous-random-network (CRN) model of Zachariasen. The chemistry of sil-icon demands that an ideal CRN model of a -Si should sat-isfy the following properties: (i) every Si atom in the networkmust be bonded to four neighboring Si atoms; (ii) the net-work must exhibit a well-defined short-range order and pos-sibly an intermediate-range order, the former being charac-terized by narrow bond-length and bond-angle distributions(e.g., the average bond angle and its variance should be closeto 109.47 ◦ and 9 ◦ –11 ◦ , respectively); (iii) any deviation fromthe ideal 4-fold coordination of the atoms must be as mini-mal as possible, preferably 1–100 in 10 atoms, so that the re-sulting electronic, optical, and vibrational properties, obtainedfrom such a model of a -Si, are in agreement with experimen-tal data. Following Weaire and Thorpe, and Heine, it canbe shown, by employing a simple tight-binding Hamiltonian,that an atomistic model with these properties exhibits a gap a) Electronic mail: [email protected] b) Electronic mail: [email protected] in the electronic spectrum. Despite these simple geometricalproperties, atomistic modeling of defect-free a -Si networks,using molecular-dynamics simulations, has been proved to beparticularly challenging and, at present, no defect-free largemolecular-dynamics models of a -Si exist in the literature toour knowledge.The current approaches to structural modeling of a -Si canbe broadly classified into three categories. The first and fore-most is based on Monte Carlo simulations, using the so-calledbond-switching algorithm of Wooten, Winer and Weaire (W3). Here, one starts from a disordered silicon crystal with100% 4-fold coordinated atoms and introduces a series ofbond switches in the network, which are followed by total-energy relaxations using the Keating potential. The bondswitching between a pair of atoms changes the network topol-ogy, but keeps the atomic-coordination number constant, byincorporating mostly 5- and 7-member rings in the network.A repeated application of bond switching, followed by total-energy relaxations, drives the system stochastically on thepotential-energy surface (PES) from one local minimum toanother over the energy barriers on the PES. In order for thesystem to be able to explore the relevant part of the config-uration space, which is associated with topologically-distinctconfigurations than a crystal, it is necessary to conduct a min-imal number of bond switches so that the system can escapefrom the initial (crystalline) state. An efficient implementationof the algorithm was presented by Barkema and Mousseau, where the method was modified to start from a random config-uration and introducing local relaxations, upon bond switch-ing, followed by intermittent total-energy relaxations of thenetwork. The method is capable of generating large 100%defect-free a -Si networks with a very narrow bond-angle dis-tribution. Molecular dynamics simulations, on the other hand, pro-vide an alternative route, where one attempts to generate a -Siconfigurations by quenching from a molten state of Si at hightemperature. Starting from a random configuration, with amass density close to the experimental density of 2.25 g/cm for a -Si, the temperature of the system is increased well abovethe melting point at 1687 K of c -Si. After thermalization, thetemperature of the system is gradually reduced, typically atthe rate of 5–10 K/ps, until the final temperature decreasesto 300 K. The approach is commonly known as ‘quench-from-the-melt’ and it is particularly effective for amorphoussolids having strong glass-forming ability. It has been ob-served that a -Si models obtained from the melt-quench ap-proach contain a high density ( ≥ Such a high density of coordination defects ad-versely affects the optoelectronic properties of the material,which renders MD models unsuitable for predictive studiesof electronic and optical properties of a -Si and a -Si:H. Theapparent failure of melt-quench approaches to produce satis-factory models of a -Si is not particularly surprising in viewof the fact that a -Si is not a glass and that it is generallyprepared in laboratories by vapor deposition on a cold sub-strate or by similar methods. The lack of glassy behavior in a -Si can be qualitatively understood from Phillips’ constraint-counting approach and Thorpe’s rigidity-percolation con-cept. In the latter, the glassy nature of a disordered solidwas shown to be connected with the mean coordination num-ber of the atoms in a network. Specifically, a disordered net-work with a mean coordination number r av < . behaves asa polymeric glass , whereas networks with r avg > . , such as a -Si with r avg ≈ , behave like rigid amorphous solids . Theseobservations, along with a high cooling rate and a rather shorttotal simulation time, appear to indicate that realistic modelsof a -Si is unlikely to be produced by molecular-dynamics sim-ulations. Thus, atomistic modeling of a -Si using direct MDsimulations has attracted little attention in recent years anda great majority of MD studies on a -Si were conducted in thepast decades. In the last decade, a number of hybrid approaches weredeveloped that ushered in a new direction for modelingcomplex materials using information paradigm.
Thesemethods employ prior knowledge of materials from exper-iments, often in conjunction with structural, chemical and electronic constraints, and combine the informationwith an appropriate total-energy functional for structural de-termination of complex solids. Motivated by the ReverseMonte Carlo method and its drawback, where a three-dimensional structural model is generated by inverting one-dimensional experimental diffraction data, hybrid approachesgo a step further by incorporating experimental data witha total-energy functional in search for structural solutionsin an augmented solution space, so that the final structureis in agreement with both experiments and theory. Exam-ples of such approaches include experimentally constrainedmolecular relaxation (ECMR), ab initio random struc-ture searching (AIRSS), first-principles assisted structuralsolutions (FPASS), and force-enhanced atomic relaxations(FEAR). While the methods have achieved a remark-able success in determining structures of a variety of com-plex solids, none of the methods has been successful so farin producing high-quality, defect-free configurations of a -Si. In this paper, we have addressed this problem and pre- sented a solution by producing nearly defect-free large modelsof amorphous silicon using molecular-dynamics simulations,supported by experimental evidence from X-ray diffractionand Raman spectroscopy.The layout of the paper is as follows. In section II, wediscuss a dynamical approach to modeling amorphous sili-con using classical molecular-dynamics simulations in canon-ical and microcanonical ensembles, followed by ab initio total-energy relaxation within the framework of the density-functional theory. We successfully demonstrate that high-quality structural models of amorphous silicon, with defectconcentrations as low as 0.7% and root-mean-square devia-tions (RMS) of bond angles below 10 ◦ , can be obtained di-rectly from ultra-long MD simulations that last for several tensof nanoseconds. Section III discusses the results of our sim-ulations with an emphasis on structural, electronic and vibra-tional properties of the models. We show that the results fromthe MD models are in excellent agreement with experimentaldata and that from the W3 models. This is followed by theconclusions of our work in section IV. II. COMPUTATIONAL METHODOLOGYA. Dynamical approach to a -Si using MD simulations Classical molecular-dynamics simulations were performedfor several systems ranging from 300 to 400,000 atoms. Thedetails of the simulations are as follows. (i)
Initial ran-dom configurations : Initial configurations of different sizes, N = of a -Si of 2.25 g/cm . (ii) MDsimulations : MD simulations, using the initial random con-figurations as input, were carried out in the canonical (NVT)and microcanonical (NVE) ensembles. The interatomic inter-action between Si atoms were described using the modifiedStillinger-Weber (SW) potential, which is given by, V ( R N ) = 12 N X i =1 N X j =1( j = i ) v ( r ij ) + N X i =1 N X j =1( j = i ) N X k =1( k = i )( k>j ) v ( r ij , r ik ) , (1)where R N indicates the atomic configuration and v ( r ij ) isthe two-body contribution to the potential energy given by, v ( r ij ) = ǫA (cid:20) B (cid:16) r ij σ (cid:17) − p − (cid:21) exp (cid:18) σr ij − aσ (cid:19) × Θ ( aσ − r ij ) , (2) TABLE I. Modified Stillinger-Weber potential energy parametersdue to Vink et al . ǫ (eV) λ σ ( ˚A) γ A B a p and v ( r ij , r ik ) is the three-body contribution to the potentialenergy, v ( r ij , r ik ) = ǫλ exp (cid:18) σγr ij − aσ + σγr ik − aσ (cid:19) × (cid:18) cos(ˆ r ij · ˆ r ik ) + 13 (cid:19) × Θ ( aσ − r ij ) Θ ( aσ − r ik ) , (3)with Θ being the Heaviside step function. Here, we have usedthe modified potential parameters, due to Vink et al , whichare listed in Table I.The equations of motion were integrated, using thevelocity-Verlet algorithm, with a time step of ∆ t = 1 fs, and achain of Nos´e-Hoover thermostats was used to control thesimulation temperature and the time evolution of the system incanonical ensembles. Each system was initially equilibratedat 1800 K for 50 ps and then it was cooled to 300 K, over a to-tal time period of 600 ps, by gradually decreasing the temper-ature in steps of 100 K with an average cooling rate of 2.5 Kper ps. At the end of the NVT-dynamics at 300 K, the systemwas subjected to evolve in a microcanonical ensemble (NVE)for an additional 50 ps so that the total simulation time for asingle NVT-NVE cycle from an initial temperature of 1800 Kto the final temperature of 300 K was 700 ps. The purpose ofthe NVE run was to eliminate any possible artifacts that couldoriginate from the thermostats and to collect structural con-figurations during the NVE dynamics, which were free fromany thermal fluctuations. Silicon configurations that satisfieda set of predefined criteria, such as those with a total potentialenergy and the density of coordination defects less than someprescribed values were collected during the 50-ps NVE run.On completion of the NVE run, the temperature of the sys-tem was increased to 1800 K to initiate the next NVT-NVEcycle. In order for the system to be able to explore a consider-able part of the phase space in searching for better amorphousconfigurations, the NVT-NVE cycle was repeated many timesuntil there were no significant changes in the total energy anddensity of coordination defects. In this study, we carried out20-60 such cycles – depending on the size of the systems –for a duration of 700 ps per cycle, resulting in a total simula-tion time of 14-42 nanoseconds. To demonstrate the efficacyof such NVT-NVE cycles, we have plotted in Fig. 1 the timeevolution of the potential energy per atom and the fraction ofthe 4-fold coordinated atoms at the end of each cycle versusthe number of cooling cycles. For the purpose of demonstra-tion only, the data shown in Fig. 1 were generated using a highvalue of the cooling rate, 5 K/ps, for a 1000-atom SW model.(iii) Optimized configurations : As mentioned before, a num-ber of low-energy configurations with desired structural prop-erties were collected during NVE evolution at the end of each -f o l d c oo r d i n a ti on ( % ) SW 1000 atoms-3.08-3.07-3.06-3.05-3.04-3.03-3.02 P E ( e V / a t o m ) SW 1000 atoms (a)(b)
FIG. 1. (a) Time evolution of the potential energy per atom versusnumber of NVT-NVE cycles. In this illustrative plot, each cycle cor-responds to a time period of 300 ps when the simulation temperaturedecreases from 1800 K to 300 K. (b) Evolution of the number oftotal atoms (in %) having 4-fold coordination versus the number ofNVT-NVE cycles. cycle in step (ii). The low-energy configurations were relaxedby minimizing the total energy using the modified SW po-tential with respect to the atomic positions. The energy min-imization was carried out using the limited-memory BFGSalgorithm.
For comparison, W3 configurations of sizes300, 1000, and 5000 atoms, each with a mass density of2.25 g/cm , were also generated following the prescription ofBarkema and Mousseau (see upper panel of Table II). B. Ab Initio relaxation using density-functional theory
The electronic structure of the SW-optimized models of a -Si with 300, 1000, and 5000 atoms was studied us-ing the local-basis density-functional code S IESTA . Norm-conserving Troullier-Martins pseudopotentials, which werefactorized in the Kleinman-Bylander form, were employed inthis work. The 300- and 1000-atom models were optimizedfully self-consistently using the Perdew-Burke-Ernzerhof(PBE) formulation of the generalized-gradient approxima-tion (GGA) of the exchange-correlation energy, and double- ζ with polarization (DZP) basis functions. For the modelwith 5000 atoms, we resorted to the Harris-functional ap-proach, which is based on the linearization of the Kohn-Sham equations in the density-functional theory. We also em-ployed single- ζ (SZ) basis functions and the local density ap- r (Å) -303691215 G ( r ) SW 1000 atomsW3 1000 atoms
FIG. 2. (Color online) The variation of the reduced pair-correlationfunctions ( G ( r ) ) with radial distances ( r ) for 1000-atom SW (red)and W3 (blue) models. proximation (LDA) to treat the exchange-correlation effectsby using the Perdew-Zunger formulation of the LDA. Sim-ilar S
IESTA -based ab initio total-energy optimizations of W3models of same sizes and mass density were also carried outfor comparison with the SW models (cf. lower panel of Ta-ble II). The vibrational density of states of 1000-atom SW andW3 models was computed, using the harmonic approxima-tion, by constructing the dynamical matrix of the system. Thelatter was obtained by computing the average atomic force oneach atom via numerical differentiation of the total energy,using an atomic displacement of × − ˚A. The resultingdynamical matrix was diagonalized to obtain the vibrationalfrequencies and eigenvectors. III. RESULTS AND DISCUSSION
We begin by addressing the structural properties of the SWand W3 models. As mentioned earlier, the quality of CRNmodels of a -Si is primarily determined by: a) the radial distri-bution function; b) the bond-angle distribution (BAD) and itswidth; and c) the number of 4-fold coordinated atoms presentin the network. A high-quality CRN model must satisfy allthese properties simultaneously. The bond-angle distributionmust have an width, ∆ θ , in the range 9–11 ◦ and it must nothave coordination defects more than 1 in 1000 atoms in orderto be experimentally compliant. Currently, no MD models ex-ist that can satisfy these requirements to our knowledge. Inthis section, we shall show that the MD models produced inour work, by combining NVT and NVE simulations discussedin Sec. IIA, closely satisfy all these requirements with a de-fect density of 0–1% and a root-mean-square (RMS) bond-angle width of 9–10 ◦ for models with size up to 1000 atoms.Figure 2 shows the reduced pair-correlation function (PCF), G ( r ) , versus radial distance ( r ) for S IESTA -optimized 1000-atom SW and W3 models of a -Si. The reduced PCF is definedas, G ( r ) = 4 πρ r ( g ( r ) − , where ρ is the number den-sity of Si atoms and g ( r ) is the conventional pair-correlationfunction. The reduced PCFs from the W3 and SW modelsmatched closely in Fig. 2, showing an almost identical nature Wave vector magnitude k (Å -1 ) S t r u c t u r e f ac t o r S ( k ) SW 1000 atomsExperiment
FIG. 3. (Color online) Comparison of the static structure factors of a -Si from molecular-dynamics simulations (SW) and high-energy X-ray diffraction. The experimental data (blue) correspond to annealedsamples of a -Si from Laaziri et al. of two-body correlations between the atoms in the respectivemodels. This observation is also reflected in Fig. 3, where wehave plotted the static structure factor, S ( k ) = * N N X i N X j exp [ ı k · ( r i − r j )] + , (4)of a 1000-atom SW model, along with the experimentalstructure-factor data for annealed samples of a -Si, due toLaaziri et al. The symbol hi in Eq. 4 indicates the rotationalaveraging of the wave-vector transfer k over a solid angle of π , along 20000 different directions. A slightly high valueof the first-sharp-diffraction peak (FSDP) for the SW modelcould be attributed partly to the use of the classical SW poten-tial, producing a pronounced structural ordering, and in part toinstrumental broadening associated with the collection of X-ray diffraction data. A relatively slow cooling rate could alsoplay a part in producing a somewhat more ordered structure.We shall return to this point later to discuss the topologicalorder and ring statistics of the SW models.Table II lists various structural properties of the SW andW3 models, from the average bond length and bond angle tothe percentage of n -fold coordinated atom with n = 3, 4, and5. We have used a cutoff value of 2.8 ˚A to calculate the co-ordination number of an atom, which corresponds to the firstminimum of the PCF and it is practically identical for all mod-els. The lower panel of Table II shows the results obtainedfrom ab initio relaxations of the corresponding classical mod-els. An examination of Table II shows that the SW modelscontain a small fraction of coordination defects. The percent-age of 4-fold coordinated atoms can be seen to decrease onlyby 2% as the system size increase from 300 atoms to 400,000atoms. The defect density in the 300- and 1000-atom modelswas found to be no more than 1% and, in each case, the RMSwidth ( ∆ θ ) of the bond-angle distribution was observed to beless than 10 ◦ . To our knowledge, these are the best configu-rations so far obtained from any molecular-dynamics simula-tions to date. Although the models with 25000 and 400,000Si atoms have a somewhat high defect concentration, in the TABLE II. Structural properties of SW and W3 models. N , h θ i , ∆ θ , h n i , and C k indicate the number of atoms, the average bond length( ˚A) and bond angle (degree), the RMS width of the bond angles (de-gree), the mean coordination number, and the percentage of k -foldcoordinated atoms, respectively. SW and W3 denote the results frommolecular-dynamics simulations and the bond-switching algorithmof Wooten, Winer and Weaire, respectively. Model
N d h θ i ∆ θ h n i C (%) C (%) C (%)W3 and SW-MD ModelsSW 300 2.38 109.22 9.06 4 0.3 99.3 0.3W3 300 2.36 109.26 9.74 4 0 100 0SW 1000 2.39 109.19 9.69 4 0.7 99 0.3W3 1000 2.37 109.25 9.06 4 0 100 0.0SW 5000 2.39 109.18 9.62 4 1.0 98 1.0W3 5000 2.36 109.22 9.84 4 0 100 0SW 25000 2.38 109.24 9.12 4 1.2 97.8 1.0SW 400000 2.39 109.22 9.46 4 1.2 97.6 1.2 Ab Initio -relaxed W3 and SW-MD ModelsSW 300 2.38 109.20 9.44 4 0.3 99.3 0.3W3 300 2.37 109.12 10.73 4 0 100 0SW 1000 2.39 109.23 9.40 4 0.4 99 0.6W3 1000 2.37 109.16 10.08 4 0 100 0.SW 5000 2.39 109.19 9.49 4 1.0 98 1.0W3 5000 2.37 109.11 10.71 4 0 100 0 range (2.2–2.4)%, the results can be readily appreciated inview of the large phase-space volume involved in the simu-lation of these ultra-large models. We should emphasize thefact that the MD models reported here were obtained fromunconstrained atomic dynamics without using any coordina-tion constraints unlike the W3 models, where a 4-fold nearest-neighbor list was always maintained during simulation.Earlier we have mentioned that the PCF and the concentra-tion of 4-fold coordinated atoms alone are not sufficient to de-termine the quality of CRN models of a -Si. For a high-qualityCRN model, the distribution of bond angles must be suffi-ciently narrow, with an RMS width of the bond-angle distri-bution in the range 9–11 ◦ , as observed in X-ray and Ramantransverse-optic (TO) peak measurements. Thus, in additionto the correct PCF, a CRN model of a -Si must have appropri-ate values of C , h θ i , and ∆ θ so that the computed propertiesfrom the model are consistent with electron spin resonance, X-ray diffraction, and Raman spectroscopy, respectively.It is notable that all the reported values of ∆ θ for the SWmodels, including the one with 400,000 Si atoms, in Table IIlie in the range of 9 ◦ –10 ◦ . The dihedral-angle distributionsfor the 1000-atom SW and W3 models are shown in Fig. 4.The distributions are quite similar with a characteristic peakat 60 ◦ .Further characterization of SW and W3 models of a -Si ispossible by examining the site-averaged orientational orderparameter (OOP), Q l , due to Steinhardt et al. In spherical φ [º]00.511.522.53 P ( φ ) SW 1000 atomsW3 1000 atoms
FIG. 4. (Color online) The distribution ( P ( φ ) ) of dihedral angles ( φ )in a -Si for SW (red) and W3 (blue) models of size 1000 atoms. polar coordinates, the site-projected OOP is defined as, Q il = vuuut π l + 1 l X m = − l (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n i X j ∈{ n i } Y ml ( θ ( r ij ) , φ ( r ij )) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (5)and the corresponding site-averaged value follows from, Q l = 1 N N X i =1 Q il . In Eq. (5), Q il is the site-projected OOP associated with site i , whose value depends on the orientation of n i bonds thatextend from site i to the neighboring sites j . Here, N is thesystem size, and θ and φ are the polar and azimuthal anglesof a bond r ij , respectively. Figure 5 show the magnitude of Q l , l = 1 to l = 8, for SW, W3, and crystalline networks, eachconsisting of 1000 atoms. A crystalline silicon network (indiamond structure) is characterized by a null or zero value of Q , Q , and Q . By contrast, SW and W3 models show astrong presence of Q indicating the amorphous character ofthe networks. It is notable that Q l values for SW and W3models practically match with each other, reflecting an iden-tical nature of local environment as far as the orientational or-dering of the atoms in the first coordination shell (of an atom)is concerned.Having addressed key structural properties of the SW mod-els, we now examine their topological character in relation tothe W3 models. Since the latter were generated by inducing,mostly, 5- and 7-member rings in a network, it is instructiveto examine to what extent the ring statistics from W3 mod-els, obtained from an event-based bond-switching algorithm,compare with the same from SW models. Toward that end,we have computed the number of n -member irreducible rings,from n =3 to n =10, and listed the values in Table III. Mathe-matically, an irreducible ring of size n is defined as the short-est, self-avoiding, irreversible path, which starts and ends atthe same atomic site in n steps. Here, irreducibility refers tothe fact that such a ring cannot be partitioned further into asmaller set of rings without changing the topology of the path l Q l SW 1000 atomsW3 1000 atoms c -Si 1000 atoms FIG. 5. (Color online) The site-averaged bond-orientational orderparameter, Q l , versus l for SW (red), W3 (blue), crystalline silicon(green) networks.TABLE III. Irreducible ring statistics (rings/per atom) for SW andW3 models.Model N Ab Initio -relaxed W3 and SW-MD ModelsSW 1000 0.002 0.247 0.662 0.372 0.062 0.012 0.003W3 1000 0.006 0.324 0.498 0.320 0.075 0.020 0.005SW 5000 0.003 0.297 0.782 0.453 0.078 0.017 0.001W3 5000 0.020 0.391 0.566 0.373 0.121 0.022 0.003 or circuit. Figure 6 presents a comparison of the ring-size dis-tribution between an SW model and a W3 model of size 1000atoms. In computing the irreducible ring-size distributions,we have not imposed periodic boundary conditions. Sincecrystalline silicon, without any defects, is characterized by 6-member rings only, the fraction of 6-member rings present in anetwork can be used as a measure of the degree of crystallinityof the network from a topological point of view. Thus, as faras the ring statistics of the W3 and SW models are concernedin Fig. 6, the latter can be described as having a somewhatpronounced topological ordering in comparison to the former.This observation is reflected in the FSDP of the SW models inFig. 3.A good atomistic configuration of a -Si must exhibit a gapin the electronic density of states (EDOS), irrespective of themethod of preparation or modeling. The size of the gap de-pends on a number of factors, such as the density of coordina-tion defects, the type of the defects (e.g., dangling and floatingbonds), and the degree of disorder in the network. In general,a high density of 3-fold coordinated atoms or dangling bondsresults in a large number of electronic states to appear nearthe Fermi level, which lead to a noisy or gap-less density ofstates. Figure 7 shows the EDOS for the 300- and 1000-atomSW and W3 models of a -Si, with the Fermi level at 0 eV. Forthe 300-atom SW model, we get a reasonably clean gap withtwo defects near the Fermi level, whereas the 1000-atom SW R i ng s p e r a t o m SW 1000 atomsW3 1000 atoms
FIG. 6. (Color online) The distribution of irreducible rings/atom inS
IESTA -relaxed 1000-atom SW (red) and W3 (blue) models. A highvalue of 6-member rings/atom in the SW model is indicative of amore topologically-ordered structure. See text for discussions. -1 -0.5 0 0.5 1 E (eV) DO S ( s t a t e s / e V ) SW 300 atomsW3 300 atoms -15 -10 -5 0 5 100510 (a) -1 -0.5 0 0.5 1 E (eV) DO S ( s t a t e s / e V ) SW 1000 atomsW3 1000 atoms -15 -10 -5 0 5 100100200300400 (b)
FIG. 7. (Color online) Depiction of the electronic densities of states(EDOS) of a -Si from SW (red) and W3 (blue) models in the regionsclose to their band-gaps. The full EDOS for the 300- and 1000-atommodels are shown as inset in (a) and (b), respectively. The Fermilevel is indicated by a dashed vertical line at 0 eV. model produces a somewhat small gap due to 1% coordina-tion defects. At any rate, a band gap of size 0.6-0.8 eV hasbeen realized in the SW models. A few defects states that arepresent in the vicinity of the band-gap region can be readilypassivated by hydrogenation to further improve the quality ofthe MD models.Finally, a true atomistic model must reproduce correctly the ω (cm -1 )0204060 VDO S Experiment SW 1000 atomsW3 1000 atoms
FIG. 8. (Color online)
Ab initio vibrational densities of states from1000-atom SW and W3 models. Inelastic neutron-scattering exper-imental data, due to Kamitakahara et al. , are also included in theplot for comparison. vibrational density of states (VDOS). Since vibrational exci-tations are more sensitive to local atomic environment thantheir electronic counterparts, any minute changes in the localatomic structure of a -Si can readily reflect on the VDOS atroom temperature. The VDOS was computed by construct-ing the force-constant (FC) matrix, using the harmonic ap-proximation, and diagonalizing the corresponding dynamicalmatrix. Figure 8 shows the VDOS of a -Si, obtained from1000-atom SW and W3 models. To compare the results withexperimental data, we have also plotted the experimental in-elastic neutron-scattering data from Kamitakahara et al. Thevibrational density of states from the SW and W3 modelsmatches closely with the experimental data.
IV. CONCLUSIONS
In this paper, we have addressed a long-standing problem ofstructural modeling of a -Si using molecular-dynamics (MD)simulations. It is generally believed that, since a -Si is not aglass, direct MD simulations of a -Si, by quenching from amolten state of Si at high temperature, are not possible. Areview of numerous MD studies on a -Si appears to supportthis conjecture by noting that MD simulations tend to producea high density of coordination defects, which is inconsistentwith the dangling-bond density observed in a -Si thin-films,by electron spin resonance measurements. This observation isconsistent with the fact that a -Si samples are prepared in lab-oratories not by melt-quenching but by vapor deposition on acold substrate or equivalent methods. Despite these observa-tions, in this paper we have shown conclusively that it is pos-sible to simulate high-quality CRN models of amorphous sili-con by employing ultra-long molecular-dynamics simulationsspanning several tens of nanoseconds. The size of the mod-els produced in this study ranges from 300 atoms to 400,000atoms and the corresponding defect density lies between 0.7%to 2.4%, respectively. All the MD models exhibit a very nar-row bond-angle distribution, characterized by an average bondangle of ∼ ∼ ◦ –9.7 ◦ . The latter is fully consistent with the valueextracted from Raman TO peak measurements. Our resultscan be fully appreciated by noting that the 300-atom modelhas only a pair of defects – one dangling bond and one float-ing bond – and a bond-angle width of 9.44 ◦ , which exhibits aclean gap in the electronic density of states. On the other hand,the largest model with 400,000 atoms has a defect density of2.4% only and an RMS width of 9.46 ◦ in the bond-angle dis-tribution.We conclude this section with the following observations.First, the availability of nearly defect-free MD models of a -Si brings considerable advantages in studying a -Si and a -Si-related materials. A molecular-dynamical approach providesa natural route to produce atomistic models of a -Si. The ap-proach is simple and intuitive and it can be implemented eas-ily and efficiently for large models, both in serial and parallelcomputing environments. Second, unlike the W3 method, amolecular-dynamical approach is free from any constraints.In the former, the system evolves from one configuration toanother, by employing a limited set of bond switches, suchthat the coordination number of the atoms remains constant.While such bond reassignments can efficiently induce the req-uisite 5- and 7-member rings to produce the characteristic net-work topology of a -Si, the system can only explore a coarse-grained configuration space due to the coordination constraintand a limited combination of bond switches. Thus, the topo-logical character of the resulting amorphous networks fromsuch an approach may not be necessarily identical to the onesobtained from unconstrained MD simulations. This observa-tion is reflected in the ring statistics obtained from the SWand W3 models. Since rings can play an important role indetermining the stability and topological ordering of amor-phous/disordered networks, it remains to be seen to what ex-tent MD models and W3 models differ from each other and af-fect the material properties, involving higher-order correlationfunctions. Finally, the availability of large MD models will beparticularly useful in studying silicon-heterojunction photo-voltaic devices based on amorphous/crystalline interfaces. ACKNOWLEDGMENTS
The work was supported by U.S. National Science Foun-dation under Grants No. DMR 1507166 and No. DMR1507118. The authors thank Profs. Gerard Barkema (Utrecht,The Netherlands) and Normand Mousseau (Quebec, Canada)for providing their modified WWW code. The authors greatlyappreciate the discussions with Profs. Stephen Elliott (Cam-bridge, UK) and David Drabold (Athens, Ohio). One of us(RAF) acknowledges the Texas Advanced Computing Center(TACC) at The University of Texas at Austin for providingHPC resources that have contributed to the results reported inthis work. M. Taguchi, A. Yano, S. Tohoda, K. Matsuyama, Y. Nakamura, T. Nishi-waki, K. Fujita, and E. Maruyama, IEEE J. Photovolt. , 96 (2014). T. Mishima, M. Taguchi, H. Sakata, and E. Maruyama, Sol. Energ. Mat.Sol. Cells , 18 (2011). M. Veldhorst, C. H. Yang, J. C. C. Hwang, W. Huang, J. P. Dehollain, J. T.Muhonen, S. Simmons, A. Laucht, F. E. Hudson, K. M. Itoh, A. Morello,and A. S. Dzurak, Nature , 410 EP (2015). R. Hosemann and S. Bagchi,
Direct Analysis of Diffraction by Matter , Se-ries in Physics (North-Holland Pub. Co., 1962). The readers may note that paracrystalline models, i.e., CRN models con-taining small clusters of ordered Si atoms, were employed in the past to ex-plain the origin of the intermediate-range order in a -Si, observed in Fluctu-ation Electron Microscopy. For an alternative explanation, based on voidsin CRN models, see Refs. 63 and 64. J. Walter and H. Eyring, The Journal of Chemical Physics , 393 (1941). G. Bartenev,
The structure and mechanical properties of inorganic glasses (Wolters-Noordhoff, 1970). W. H. Zachariasen, Journal of the American Chemical Society , 3841(1932). D. Weaire and M. F. Thorpe, Phys. Rev. B , 2508 (1971). V. Heine, Journal of Physics C: Solid State Physics , L221 (1971). F. Wooten, K. Winer, and D. Weaire, Phys. Rev. Lett. , 1392 (1985). B. R. Djordjevi´c, M. F. Thorpe, and F. Wooten,Phys. Rev. B , 5685 (1995). P. N. Keating, Phys. Rev. , 637 (1966). G. T. Barkema and N. Mousseau, Phys. Rev. B , 4985 (2000). R. L. C. Vink, G. T. Barkema, M. A. Stijnman, and R. H. Bisseling,Phys. Rev. B , 245214 (2001). J. S. Custer, M. O. Thompson, D. C. Jacobson, J. M. Poate, S. Roorda,W. C. Sinke, and F. Spaepen, Applied Physics Letters , 437 (1994). M. H. Brodsky and R. S. Title, Phys. Rev. Lett. , 581 (1969). J. Phillips, Journal of Non-Crystalline Solids , 153 (1979). M. Thorpe, Journal of Non-Crystalline Solids , 355 (1983). A. Pedersen, L. Pizzagalli, and H. Jnsson,New Journal of Physics , 063018 (2017). J. F. Justo, M. Z. Bazant, E. Kaxiras, V. V. Bulatov, and S. Yip,Phys. Rev. B , 2539 (1998). R. Car and M. Parrinello, Phys. Rev. Lett. , 204 (1988). D. A. Drabold, P. A. Fedders, O. F. Sankey, and J. D. Dow,Phys. Rev. B , 5135 (1990). N. Cooper, C. Goringe, and D. McKenzie,Computational Materials Science , 1 (2000). E. Kim and Y. H. Lee, Phys. Rev. B , 1743 (1994). I. ˇStich, R. Car, and M. Parrinello, Phys. Rev. B , 11092 (1991). P. Klein, H. M. Urbassek, and T. Frauenheim,Computational Materials Science , 252 (1999). M. D. Kluge, J. R. Ray, and A. Rahman, Phys. Rev. B , 4234 (1987). W. D. Luedtke and U. Landman, Phys. Rev. B , 4656 (1988). P. Biswas, R. Atta-Fynn, and D. A. Drabold,Phys. Rev. B , 125210 (2007). A. Pandey, P. Biswas, and D. A. Drabold,Scientific Reports , 33731 EP (2016). K. Prasai, P. Biswas, and D. A. Drabold,Scientific Reports , 15522 EP (2015). B. Meredig and C. Wolverton, Nature Materials , 123 EP (2012). C. J. Pickard and R. J. Needs, Journal of Physics: Condensed Matter , 053201 (2011). P. Biswas and S. R. Elliott, Journal of Physics: Condensed Matter , 435201 (2015). A. Pandey, P. Biswas, and D. A. Drabold, Phys. Rev. B , 155205 (2015). R. L. McGreevy and L. Pusztai, Molecular Simulation , 359 (1988). O. Gereben and L. Pusztai, Phys. Rev. B , 14136 (1994). J. K. Walters and R. J. Newport, Phys. Rev. B , 2405 (1996). P. Biswas, R. Atta-Fynn, and D. A. Drabold,Phys. Rev. B , 195207 (2004). M. J. Cliffe, A. P. Bart´ok, R. N. Kerber, C. P. Grey, G. Cs´anyi, and A. L.Goodwin, Phys. Rev. B , 224108 (2017). P. Biswas, D. N. Tafen, and D. A. Drabold,Phys. Rev. B , 054204 (2005). R. L. C. Vink, G. T. Barkema, W. F. van der Weg, and N. Mousseau,Journal of Non-Crystalline Solids , 248 (2001). F. H. Stillinger and T. A. Weber, Phys. Rev. B , 5262 (1985). G. J. Martyna, M. E. Tuckerman, D. J. Tobias, and M. L. Klein,Mol. Phys. , 1117 (1996). W. G. Hoover, Phys. Rev. A , 1695 (1985). S. Nos´e, J. Chem. Phys. , 511 (1984). J. Nocedal, Mathematics of Computation , 773 (1980). D. C. Liu and J. Nocedal, Mathematical Programming , 503 (1989). J. M. Soler, E. Artacho, J. D. Gale, A. Garci´a, J. Junquera, P. Ordej´on, andD. S´anchez-Portal, Journal of Physics: Condensed Matter , 2745 (2002). N. Troullier and J. L. Martins, Phys. Rev. B , 1993 (1991). L. Kleinman and D. M. Bylander, Phys. Rev. Lett. , 1425 (1982). J. P. Perdew, K. Burke, and M. Ernzerhof,Phys. Rev. Lett. , 3865 (1996). J. Harris, Phys. Rev. B , 1770 (1985). J. P. Perdew and A. Zunger, Phys. Rev. B , 5048 (1981). K. Laaziri, S. Kycia, S. Roorda, M. Chicoine, J. L. Robertson, J. Wang, andS. C. Moss, Phys. Rev. B , 13520 (1999). D. Beeman, R. Tsu, and M. F. Thorpe, Phys. Rev. B , 874 (1985). P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B , 784(1983). W. A. Kamitakahara, H. R. Shanks, J. F. McClelland, U. Buchenau,F. Gompf, and L. Pintschovius, Phys. Rev. Lett. , 644 (1984). Since the energy associated with vibrational excitations in solids, involvinga fraction of an electron-volt, is significantly smaller than their electroniccounterpart (with a few to several electron-volts), a small perturbation inthe local atomic environment of atoms can be readily manifested in thevibrational spectrum of a solid. At high temperature, such effects can beconsiderably weakened or washed out due to thermal broadening. S. V. Khare, S. M. Nakhmanson, P. M. Voyles, P. Keblinski, and J. R.Abelson, Applied Physics Letters , 745 (2004). J. Gibson, M. Treacy, and P. Voyles, Ultramicroscopy , 169 (2000). P. Biswas, R. Atta-Fynn, S. Chakraborty, and D. A. Drabold,Journal of Physics: Condensed Matter , 455202 (2007). P. Biswas, D. N. Tafen, R. Atta-Fynn, and D. Drabold,Journal of Physics: Condensed Matter16