Hydrogen-induced degradation dynamics in silicon heterojunction solar cells via machine learning

Among silicon-based solar cells, heterojunction cells hold the world efficiency record. However, their market acceptance is hindered by an initial 0.5\% per year degradation of their open circuit voltage which doubles the overall cell degradation rate. Here, we study the performance degradation of crystalline-Si/amorphous-Si:H heterojunction stacks. First, we experimentally measure the interface defect density over a year, the primary driver of the degradation. Second, we develop SolDeg, a multiscale, hierarchical simulator to analyze this degradation by combining Machine Learning, Molecular Dynamics, Density Functional Theory, and Nudged Elastic Band methods with analytical modeling. We discover that the chemical potential for mobile hydrogen develops a gradient, forcing the hydrogen to drift from the interface, leaving behind recombination-active defects. We find quantitative correspondence between the calculated and experimentally determined defect generation dynamics. Finally, we propose a reversed Si-density gradient architecture for the amorphous-Si:H layer that promises to reduce the initial open circuit voltage degradation from 0.5\% per year to 0.1\% per year.


INTRODUCTION
With record-breaking efficiencies approaching 27% [1], silicon heterojunction (Si HJ) solar cells are rapidly becoming one of the most promising next generation technologies. A key driver for their excellent performance is the electronic passivation of the crystalline silicon (c-Si) by a thin layer of hydrogenated amorphous silicon (a-Si:H). In spite of their great promise, Si HJ solar cells' market acceptance is lagging, primarily because their efficiency degradation was reported to be about 0.7% per yr [2], much higher than the usual rate of 0.2% per yr. [3] The excess 0.5% per yr rate was attributed to the degradation of the open circuit voltage V oc . * amdiggs@ucdavis.edu † ztzhao@ucdavis.edu ‡ Also at Center for Nanoscale Materials, Argonne National Laboratory, Lemont, IL 60439, USA [2,3] (It is noted that these are initial degradation rates, measured over limited time ranges. Nonlinear aspects over longer time ranges will be discussed below.) Therefore, understanding and subsequently minimizing the performance loss of Si HJ solar cells is key to accelerating their market acceptance. Towards this goal, this paper focuses on analyzing and potentially suppressing Si HJ solar cell degradation. For completeness, we add that their market is also limited by the need for disruptive changes in the manufacturing chain, increased cost of the n-type wafers, and the low-temperature Ag metallizations, among others. Although a-Si:H has been analyzed for more than 40 years, several questions of its physics still remain open. Unresolved issues include understanding the classes of structural and electronic defects, [4][5][6][7][8][9][10][11] their formation and statistical distributions, and the long-term structural defect dynamics of a-Si:H. These issues all play crucial roles in the degradation of the crystalline/amorphous (c-Si/a-Si:H) interface of Si HJs. [12][13][14] The primary form of degradation seen in Si HJ cells is an increase in surface recombination at the c-Si/a-Si:H interface. Here, dangling bonds (DBs) are known to be the primary recombination centers. [15] DBs are inherently present, created by disorder and thermodynamics. The recombination is suppressed by passivating these DBs with the introduction of hydrogen during fabrication. The degradation over the lifetime of the cell unfolds by the increase of the DB density. The DB density can increase by different mechanisms. In a previous paper we computed the increase of the DB density by the structural evolution of the Si matrix alone. [14] We found that the structural evolution of the Si matrix did increase the DB density, but quantitatively was not sufficient to account for the data. Therefore, in this paper we will focus on another mechanism that can drive the increase of the DB density: the loss of passivating hydrogen at the interface.
Over the past 30 years, many studies investigated hydrogen dynamics in bulk a-Si:H, focusing on its drift and diffusive motion and the formation of light induced defects (LID). [16][17][18][19][20][21][22][23][24][25][26][27][28][29] Also, there have been numerous papers on the dynamics of hydrogen in a-Si:H and its role in passivating the surface of c-Si. [12,[15][16][17][18][19][20][21][22][23][24][25][26][30][31][32] More recently, studies zoomed in on the dynamics at the c-Si/a-Si:H interface; specifically, on the Si-H bonding structure and restructuring during annealing and light soaking, [33][34][35] hydrogen motion at the interface at elevated temperatures, [36] the role of hydrogen in preventing epitaxial growth, [37,38] and interface defect types and formation. [15,[30][31][32] These works on the kinetics of hydrogen in a-Si:H have provided crucial understanding of key aspects of the electronic passivation of c-Si. However, several issues remain open. Theoretically, the structural variability of amorphous Si necessitates the study of large samples which are hard to access by numerical methods. Experimentally, it has proven challenging to track hydrogen motion in silicon structures. The above challenges demonstrate the need for fur-ther theoretical and numerical research to provide additional insight. The most effective numerical methods are first principles electronicstructure studies using density functional theory (DFT). [24,28,[39][40][41] However, the computational expense of DFT scales unfavorably with O(N 3 ). Thus, in order to reach larger sample sizes, many computational studies were performed with Molecular Dynamics (MD) methods using parameterized empirical interatomic potentials, whose cost scales linearly with O(N). [42] However, this gain in computational speed was achieved at the cost of precision.
In the last decade the adoption of machine learning (ML) methods made it possible to create new interatomic potentials that improved MD methods to attain DFT level accuracy. [43][44][45][46][47][48][49][50][51][52][53][54][55][56] Motivated by these successes, a subset of the present authors very recently developed the first ML-trained Gaussian Approximation Potential (GAP) for hydrogenated silicon, the Si-H GAP. [57] It was demonstrated that the Si-H GAP was able to simulate samples of sizes unattainable for DFT with DFT level accuracy, as described in detail below.
Motivated by the above challenges, needs, and considerations, in this work we combine experimental and computational efforts to perform a comprehensive analysis of the degradation of c-Si/a-Si:H heterojunction systems. Experimentally, we have tracked the degradation in Si HJ systems over a year. Theoretically, we upgraded our SolDeg degradation simulator to track hydrogen-induced degradation in the same Si HJ systems from femtosecond to gigasecond time scales. Remarkably, we have achieved a compelling quantitative correspondence of our simulation results with the experimental data. Finally, building on the understanding developed with this SolDeg analysis, we developed an actionable proposition on how to modify the structure of Si HJ cells to dramatically suppress the initial V oc degradation rate by as much as 80%. Such a substantial reduction may considerably help the market acceptance of Si HJ cells.

Experimental Degradation Study
The goal of the experimental part of the project was to investigate the degradation dynamics of stacks of hydrogenated amorphous silicon, a-Si:H, deposited on top of crystalline silicon, c-Si. We created a set of n-type c-Si/a-Si:H stacks and measured the carrier lifetime over the period of a year under dark and ambient temperature conditions. The sample preparation and experimental techniques are described in the Methods section.
The effective minority carrier lifetime (τ eff ) of a symmetrically passivated sample with low surface recombination velocity (S) is a combination of the individual lifetimes due to bulk, Auger, Shockley-Read-Hall (SRH), radiative, and surface recombination processes: where the surface recombination rate is the only one that depends of the wafer thickness W . We determined S from the slope of the 1/τ eff versus 1/W plot, constructed from the lifetime data taken on four wafers with varying thicknesses W . The S was determined for a range of injection levels ∆n and a range of temperatures T to construct S(∆n, T ). Finally, we determined the defect density N and the charge density Q at the interface by fitting S(∆n, T ) to the amphoteric model proposed by Olibet et al.. [32] We collected the S(∆n, T ) data every two weeks over the course of a year, and used these data to construct the time dependence of N (t) and Q(t). Fig. 1 shows that the effective carrier lifetime τ eff (t) decreased over the year. Fitting the surface component of τ eff to the amphoteric model of Olibet et al. [32] showed that the defect density at the interface N (t) increased over the year, while Q(t) remained constant. Motivated by these experiments, we developed a robust theoretical framework to simulate, to From analyzing these measurements at different temperatures and different injection levels, we determined the time-dependent defect density at the c-Si/a-Si:H interface N(t) (red diamonds).
model, and to analyze the time evolution of the interface defect density N (t).
Si HJ Stacks: Structure, Energies and Defects As mentioned earlier, in a previous paper we simulated the degradation dynamics in simplified, Si-only c-Si/a-Si stacks. [14] In that paper, we used the machine learning-based GAP for the Si-Si interatomic potential. However, a crucial ingredient in the experimental samples, hydrogen, was not included in that work. Naturally, a quantitative analysis must include hydrogen. In order to perform high precision modeling, we needed a Si-H interatomic potential with a precision comparable to the Si-Si GAP.
Since such a potential was not available, we used machine learning to train and construct the Si-H GAP de novo. Here we describe our methods only in bullet points, and refer to the Methods section for details.
First, we developed the Si-H GAP using ML, as reported very recently. [57] Most importantly, MD simulations performed with the Si-H GAP were able to reproduce DFT energies with a 4 meV per atom precision. The Si-H GAP enabled us to simulate a-Si:H samples with close to 5,000 atoms, or hundreds of average-sized samples in parallel. Such scales were necessary for the high quality determination of the distributions of key descriptors, needed to account for the measured experimental data.
Next, we used Si-H GAP-powered MD simulations to create separate c-Si and a-Si:H slabs, by validated quenching protocols.
Then, we merged the c-Si and a-Si:H slabs into stacks using the same Si-H GAP MD, wherein their densities were matched at the interface. We created 60 interfaces.
Finally, we determined µ H (z), the positiondependent energy of adding a hydrogen atom to the stacks at an interstitial location a distance z from the interface. The result of 25,000 such calculations is presented in Fig. 2, which shows our first central finding: We observed that µ H (z) exhibited a gradient across the interface. This energy gradient means that mobile hydrogen atoms at the interface experience a force that causes them to drift away from the interface into the a-Si:H layer. This hydrogen drift leaves behind unpassivated Si DB, and thus degrades the cell performance. In the remainder of this paper, we focus on this discovery as the most promising mechanism of the experimentally measured increase of the defect density in Fig. 1.
Next, we explored several different possible drivers of this energy gradient. Eventually, we established that the Si density, as characterized via the Radial Distribution Function (RDF), exhibits a very analogous gradient across the interface, see Fig. 2, and is therefore the most likely driver of the hydrogen chemical potential gradient. the Si density remarkably closely, making is plausible that the Si density gradient is the cause of the H chemical potential gradient.
We then investigated the electronic properties of the hydrogen-induced defects and established that in 94% of the cases, departing hydrogens create localized dangling bonds, and in 6% of the cases they leave behind delocalized electronic states. Moreover, we showed that the average Löwdin charge of the generated dangling bonds is -0.1 e, and thus they can be properly identified as the neutral defects seen in our experiment.

Energy Barrier Distributions Connecting Femtoseconds to Gigaseconds
We connected the femtosecond time scale of MD simulations to the gigasecond time scale of the experiments by adapting the Nudged Elastic Band Method (NEB) to determine the energy barriers that control the hydrogen dynamics.
We found that the key processes were hydrogen breaking free from Si-H bonds, and hydrogen hopping between interstitial sites. We computed the energy barriers for about 650 initialfinal state pairs for the breaking of Si-H bonds, and for more than 2000 initial-final state pairs for interstitial hopping, for both the forward and backward directions. Fig. 3 shows our second central finding: The energy barriers for the four processes (bond breaking and interstitial hopping, to and from the interface) all have wide distributions. The computed mean energy barriers of these processes, E bond =1.31±0.03 eV, E inter =0.49±0.08 eV, were consistent with previous experimental and computational values: [18,23,25,26,40] E inter being the average of E drift and E reverse .
In an important distinction, however, most previous studies reported only single values or narrow distributions (σ ≈ 0.025 eV) for these barriers. In contrast, we found that the width of these distributions were quite large, with σ in the range of 0.20 eV-0.30 eV. The substantial width of the energy barrier distributions is important because at room temperatures, activated processes over the mean bond-breaking barrier of 1.31 eV are frozen during the 1 week to 10 years time scale (using attempt frequencies in the 10 13 Hz range from FTIR vibrational data), and thus could not explain the defect density measurements in Fig. 1. Only barriers in the 1.1-1.2 eV range are active in the experimental time range of 1 week-1 year. Notably, such barriers are present in our simulation results.
We also computed the interstitial energy barriers: E drift controlling the drift away from the interface into the a-Si:H layer, and E reverse controlling the drift back form the a-Si:H layer toward the interface. Most previous papers did not recognize any difference between these processes and reported a single barrier value of E inter = 0.50 eV. [18,23,26,40] However, our first main result, the existence of the hydrogen chemical potential gradient, resulted in a bias between these barriers: E drift = 0.41 eV, and E reverse = 0.58 eV. While the average of these two barriers, E inter = 0.49 eV, remained in agreement with the previously published value of 0.50 eV, the bias between them presents a physical explanation for the drift of hydrogen away from the interface.
Furthermore, we found that E recap , the barrier that controls the recapture of a hydrogen atom from an interstitial position by a Si dangling bond, was centered at 0.28 eV. Published values identified this barrier also as 0.50 eV. [18] The difference in E recap from the other interstitial barriers is a result of the local low Si density region near a dangling bond.

Degradation Dynamics
Determining the energy barriers and their distributions created the foundation on which we next built the analysis of the dynamics of hydrogen drift. First, we simulated the hydrogen drift from the interface across the first 4-16 local energy minima along the z perpendicular direction, by both Accelerated Superbasin Kinetic Monte Carlo methods, [14] and as a series of thermally activated jumps, a representation of the energy landscape used in these simulations is shown in Supplementary Figure 2. The results of these simulations illuminated the dramatic effect of the hydrogen chemical potential gradient. The bias of the energy barriers relative to k B T , E drift −Ereverse k B T ≈ 6 made the probability of a hydrogen returning to the interface less than 0.1%. This very small return probability validated reducing the dynamics to the following three hydrogen states: H bonded to a Si at the interface, H at an interface-adjacent interstitial site, and H in the bulk of the a-Si layer. The three energy barriers controlling the dynamics between these three states are the Si-H bond-breaking E bond , the bond-recapture barrier E recap , and the barrier against forward drift into the a-Si:H layer E drift , for validation of the Three Barrier Model see Supplementary Figure 3.
We constructed the following set of rate equations for the three barrier model: where the H vector denotes the hydrogen densities in the three model states, and denotes the matrix of thermally activated rate constants across the above determined energy barriers. We solve this Eq. 3 for the hydrogen density at the interface H i (t) as: where where α ≡ k 1 + k 2 + k 3 , and β ≡ √ α 2 − 4k 1 k 3 . As we established before, when hydrogens leave the interface, they leave behind dangling bonds and thus increase the defect density as ∆N i (t) = −∆H i (t). We determined the total defect density N (t) by averaging N i (t) over the distributions of the three rate constants, which we determined from the barrier distributions computed using the NEB method with our Si-H GAP. From here on, we extend the scope of the SolDeg simulator of Ref. [14] to refer to the here-described multiscale, hierarchical combination of simulations and analytical methods. Fig. 4 shows our SolDegsimulated N (t) alongside the experimentally determined N (t) from Fig. 1.
The central result of our paper is the remarkable correspondence of the measured and the SolDeg-simulated defect density N (t), see Fig.  4. This compelling correspondence is validation that our SolDeg simulator provides a quan-titatively reliable description of the hydrogeninduced degradation dynamics of heterojunction systems. Fig. 4 further shows that the measured and simulated defect densities can also be described very well with a stretched exponential form, which is a hallmark of relaxation in disordered systems with broad barrier distributions.
The time-dependence of the defect density N (t) has direct experimentally measurable consequences. The increasing defect density makes an additive contribution to the inverse career lifetime 1/τ , causing it to decrease over time. Further, following the models of Green [58] and Olibet et al. [32], or Blank et al. [59], V oc is another quantity directly effected by the defect density N (t) via where C is a constant and q is the elementary charge, for a full derivation see Supplementary Equations [4][5][6][7][8][9][10][11][12][13] . This relation connects our SolDeg simulation results to the commercially relevant degradation of V oc (t) in Si HJ modules. Additional utility and relevance of our SolDeg platform can be appreciated, e.g., from the fact that the expression for V oc (t), constructed by inserting a stretched exponential N (t) into Eq. (5), fits the experimentally measured [2] degradation of V oc (t) in fielded modules over the entire measurement range of 8 years, see Supplementary Figure 1.

Reversed Gradient NoDeg Stacks
Having gathered the above insights into the physics of the degradation of existing Si HJ cells, we developed a proposition on how to slow down this degradation. Since the Si density gradient caused the hydrogen drift from the interface, we suspected that reversing this density gradient would do the opposite; it would pin the hydrogen atoms to the interface. To test this idea, we changed the c-Si to a-Si:H slabmerging protocol to create a Si density that was lower at the interface and increased into the a-Si:H layer.   5 shows that in these reversed-gradient stacks, the hydrogen atoms feel a force that pins and localizes them to the interface, thereby ensuring that the defect generation, and thus cell degradation by hydrogen drift is suppressed. From here on we refer to heterojunction cells which are fabricated with a Si density minimum at the interface as NoDeg cells.
We repeated our SolDeg simulation of the hydrogen-induced degradation in these NoDeg stacks. The key difference was that the E Drift barrier controlling the forward drift increased as the reversed Si density gradient increased. Fig. 6 shows the remarkable result of our Sol-Deg simulations of the NoDeg cells.
Most remarkably, in these NoDeg cells, after one year the defect density N (t) was reduced relative to the original stacks by about 80%. Here we recall that the original experiment that motivated our work reported that the open circuit voltage of Si HJ cells, V oc decreased at a rate of 0.5% per yr. The above simulation of the NoDeg cells suggests that by reversing the Si density gradient at the interface, the V oc degradation rate could be suppressed to a very encouraging 0.1% per yr. As a practical matter, the Si density and its gradient can be tuned in a wide range by suitably modifying deposition temperatures, partial pressures of the SiH4, H2, and possibly Si2H6 gases, and other deposition parameters.
After the completion of our work, we learned that some experimental groups fabricated and characterized an analog of the proposed reverse gradient heterojunction cells. These groups created a bilayer architecture of the a-Si:H layer: They deposited a thin underdense a-Si:H on the c-Si, followed by the deposition of a regular a-Si:H layer. All groups measuring heterojunction cells with such bilayer architecture demonstrated the improvement of surface passivation.
Lee et al.. reported V oc improvements of 10-20 mV when comparing bilayer cells with 5 nm of low density a-Si:H, plus 5 nm regular density a-Si:H layers with cells with a 10 nm a-Si:H layer. [60] Smets et al.. reported a 10 mV V oc increase for a 1nm+8nm bilayer design relative to a 9nm single layer. [61] Sai et al.. reported a 30 mV increase for a 2nm+8nm bilayer relative to a 10 nm single layer. [38] Other groups also reported improvements in passivation, fill factor, or V oc . [37,[62][63][64] The above groups did not perform systematic degradation studies. However, we speculate that the reported V oc improvements could have been the consequence of suppressed degradation that took place from the moment of fabrication to the time of the V oc measurements, days or weeks later.
The cited bilayer papers already established that fabricating HJ cells with underdense-dense a-Si:H bilayers improved the passivation. Beyond these, our NoDeg simulations propose that even greater passivation improvements can be achieved by depositing the a-Si:H layer with a continuous reversed Si density gradient. Of course, other physical processes, like the formation of higher hydrides and morphological variations could have played a role in the measured changes of the passivation, and their roles should be explored.

Conclusion
This paper reported a comprehensive study of the performance degradation of c-Si/a-Si:H heterojunction stacks. We carried out a yearlong experimental study of the performance degradation of these stacks and determined the time dependence of their interface defect density. Then we developed SolDeg, a multiscale, hierarchical simulator, to analyze this degradation. First, we used Machine Learning to construct the most accurate Si-H interatomic Gaussian Approximation Potential GAP that reproduced DFT energies with a 4 meV per atom precision. Then we used this Si-H GAP to carry out LAMMPS Molecular Dynamics simu-lations to construct c-Si/a-Si:H heterojunction stack interfaces. We discovered that the hydrogen chemical potential developed a gradient across the interface, which forced the hydrogen atoms to drift away from the interface. The departing hydrogens left behind dangling bond defects that acted as recombination centers for the charge carriers of the c-Si wafer. We implemented the Nudged Elastic Band method to determine the barriers that control the hydrogen dynamics, and their distributions. We identified the crucial processes as the breaking of the Si-H bonds and the drift between interstitial states. We determined that the barrier distributions were characterized with remarkably large widths of about 0.30 eV. We constructed and validated a simplified model of the dynamics of the hydrogen-induced defect generation across these barriers. We found a compelling, quantitative correspondence between the calculated and the experimentally measured defect generation dynamics, as well as a promising connection to experimentally measured V oc (t) degradation in similar fielded modules out to 8 years. Finally, we built on this analysis to propose a reversed Si-density gradient architecture for the a-Si:H layer of heterojunction cells to suppress the interface defect generation. In such NoDeg cells the initial V oc degradation rate may be reduced from 0.5% per yr to 0.1% per yr.

Experimental Degradation Study
For the lifetime measurements, c-Si wafers were bifacially coated with a-Si:H(i). The wafers were double-side polished float zone (FZ) quality n-type c-Si wafers with (100) crystal orientation, 2.5Ωcm resistivity, and thickness of The deposited a-Si:H(i) films were characterized with Fourier transform infrared spectroscopy (FTIR) using Nicolet 6700 spectrometer from Thermo Electron to determine the microstructure and hydrogen content of the films. A M2000 ellipsometer from JA Woollam was employed to characterize the thickness, bandgap (Eg), refractive index, and degree of crystallinity in the a-Si:H(i) layer. Ellipsometric spectra were collected at multiple incident angles of 65, 70 and 75 degrees in reflection mode. The resulting spectra were fitted with a single Tauc-Lorentz oscillator and yielded a gap of 1.68 eV.

Development of the Si-H GAP Using ML
The Si-H GAP was developed by simulating Si samples with 12% and 15% hydrogen concentrations. We performed 27 rounds of training of the Si-H GAP on a wide variety of relevant Si:H structures, including amorphous and liquid Si, and several defect and interface structures. The Si-H GAP was validated via several key descriptors. Specifically, our training successfully reduced the difference between the Si-H GAP energy-per-atom and the DFT energy-per-atom to below 4 meV per atom. Given that the differences of the energy-per-atom calculated by different DFT methods, as well as the differences between DFT vs. experiments repeatedly exceed this value, the precision of Si-H GAP we achieved is quite remarkable. Further, the differences of the force differentials between the Si-H GAP and DFT were reduced by about 40% in the course of the training. For further context, the energy-per-atom, force fields, and stress fields were also substantially reduced compared to the results of non-GAP interatomic potentials, such as the Tersoff potential. [57] In addition, the Radial Distribution Function (RDF), shown in Fig. 7, and the Bond Angle Distribution functions were also calculated by our Si-H GAP.

Creation of c-Si/a-Si:H Stacks
Armed with the freshly-developed Si-H GAP potential, we proceeded to create hydrogenated c-Si/a-Si:H stacks. For our work, we used Molecular Dynamics simulations, in particular the LAMMPS package. The c-Si slabs contained 162 Si atoms. The a-Si slabs were created by starting with a 216-atom c-Si slab, and inserting into it 28 hydrogen atoms at random locations, creating an a-Si:H layer with 13 at.% H, a typical value for the passivating layer of Si HJs. Next, these hydrogendoped Si slabs were heated to 1,800K by using the Si-H GAP to form liquid Si:H. The liquid Si:H was subsequently re-solidified by quenching to 1,500K at a rate of 10 13 Ks −1 , and then equilibrated at 1500K for 100 ps. This solid Si:H was then further quenched down to 500K at a rate of 10 12 Ks −1 , following protocols optimized in previous studies. [45,[65][66][67] The first quench was performed in the constantvolume and variable-pressure (NVT) ensemble, while the second quench was performed in the variable-volume and constant-pressure (NPT) ensemble with fixed x and y cell-dimensions to match the dimensions of the c-Si slab in the later steps. Both quenches used a Nose-Hoover thermostat and barostat. We minimized the structural energy using a Hessian-free truncated Newton (HFTN) algorithm to relax all atomic positions into their local minima. These relaxed a-Si:H structures were further optimized with DFT, for DFT methods and parameters see Supplementary Methods 2.
In the second stage, the a-Si:H slabs were fused with the c-Si slabs, thus creating an interface and forming the c-Si/a-Si:H heterojunction, see Fig. 9. The c-Si/a-Si:H structures will be also referred to as stacks. First, we created a set of stacks by positioning the slabs at a sequence of separations. Next, we determined the z-dependent density D(z) of each stack as the density averaged within a 2.5Å thick sliding cuboid, centered at z and integrated over the lateral r = (x, y) coordinate across the entire lateral extent of the stack. This density D(z) is shown for one of the stacks in Fig. 10. Then we selected the stack with the a-Si:H -c-Si slab separation that caused the density D(z) to have the smoothest transition from the c-Si slab into the a-Si:H slab. In effect, we matched the density of the top 2.5Å of the c-Si with that of the bottom 2.5Å of the a-Si:H, indicated by the blue and yellow rectangles in Fig. 10.
The zero of the z coordinate indicates the interface, i.e. location where the slabs were fused. This method was designed to simulate the commonly used transition zone deposition procedure of a-Si:H(i) layers. The fusing was completed by a final annealing step, where we annealed the fused stacks by heating them to 1500 K, and then quenched them at a rate of The blue atoms represent silicon and the red atoms represent hydrogen.
10 13 Ks −1 . This allowed the distortions and defects generated by the fusing at the interface to diffuse or relax into the bulk of the a-Si:H slab.

Hydrogen Chemical Potential Gradient
Next, we computed several characteristics of the stacks. We began by calculating the energy µ H (z) experienced by a probe hydrogen atom inserted in a stack at a depth z. We inserted these probe hydrogen atoms into each stack at 300-500 randomly selected locations (x, y, z) and determined the total energy µ H (x, y, z) of the resulting structure. All of the probe hydrogen atoms relaxed into interstitial sites at the center of Si-Si bonds as almost all dangling bonds in the stacks were passivated by existing hydrogens. The µ H (x, y, z) energies were smoothed by averaging over all insertion locations within a 2.5Å thick sliding cuboid, just like we did with the density D(z) earlier. Finally, we slid the center of this averaging cubiod in small increments across the entire stack to construct µ H (z). Carrying out this procedure for the 60 interfaces (2 interfaces per stack), our comprehensive µ H (z) is the result of about 25,000 calculations. Fig. 2 shows the so-calculated µ H (z) hydrogen chemical potential.

Origin of the Hydrogen Chemical Potential Gradient
Next, we explored several different possible drivers of this energy gradient. Eventually, we suspected that the Si density and microstructure played a crucial role. To test this hypothesis, we analyzed the radial distribution function RDF(r, z), which was again calculated using a 2.5Å thick cuboid centered on a z coordinate, as a function of the lateral radius r = (x, y). Specifically, we characterized the local microstructures with the height of the first peak of RDF(r, z) as a function of the lateral radius r, which we denoted by RDF1(z). A high RDF1(z) indicated a high degree of crystallinity and therefore high density, which translated to a low porosity, whereas a low RDF1(z) indicated the opposite. The red line in Fig. 2 shows RDF1(z) as a function of z. RDF1(z) visibly also exhibits a marked gradient across the interface at z = 0. The correlation between the gradient of the hydrogen chemical potential µ H (z) and the gradient of RDF1(z) is conspicuous. This makes us confident to conclude that the gradient of the hydrogen chemical potential was created by a corresponding gradient of the microstructural crystallinity, and thus density across the interface of the heterjunction, as the stack transitioned from the crystalline Si to the less crystalline, and thus more porous, amorphous Si.

The Nudged Elastic Band Method
The experimental study determined the degradation dynamics up to a year. Solar cell manufacturers make guarantees or warranties out to 20-30 years (the latter being equal to about 1 gigasecond). In contrast, the time step of our LAMMPS MD calculations is less than 1 femtosecond. As such, MD calculations have no chance of reaching anywhere near the experimentally relevant time scales. In order to bridge the 24 orders of magnitude chasm between these time scales, we recognize that thermally activated stochastic processes across energy barriers in the ≈ 1eV range are the ones that control the hydrogen dynamics on the slow 1-50 weeks experimental time scale. In order to describe the degradation dynamics on such a slow scale, we must determine the energy barriers the hy-drogen atoms face up to the eV scale. (Here we assumed reasonable attempt frequencies of about 10 13 Hz, corresponding to local vibration modes, as determined by Fourier Transform Infrared spectroscopy, or FTIR.) To compute the energy barriers, we turned to the Nudged Elastic Band (NEB) method. The NEB method is based on constructing multiple replicas of a given system, corresponding to a path from an initial to a final state over a high energy region. In an informative scenario, the initial state is a hydrogen in a local minimum (bound to a Si, or at an interstitial position), the final state is the same hydrogen at another local minimum, and the replicas are intermediate configurations where the hydrogen crosses a high energy region between these minima, tracking both the displacements of the hydrogen and surrounding Si atoms. The NEB method connects these replica configurations by imaginary springs which provide interreplica forces, F sp . The total force on the atoms in each replica is the sum of the interatomic forces and inter-replica forces. The NEB method minimizes the total energy of all replicas, and thereby determines the Minimum Energy Path (MEP) between the initial and final configurations. The activation energy of such a process is then determined as the energy difference between the saddle point/maximum energy configuration and the initial configuration. Specifically, we used the climbing image NEB method, [68,69] as implemented by LAMMPS, with the Si-H GAP [57] interatomic potential and the FIRE [70] damped minimizer, for a complete set of NEB parameters see Supplementary  Table 3.

NEB Initial and Final States
Previous studies have presented a description of hydrogen in a-Si:H, where hydrogen is located in either of two classes: a deep trap, a hydrogen atom bound to an silicon atom; or a shallow trap, a hydrogen at an interstitial Si-Si bond center. [18,29,36,71,72] Following this model, we organized our NEB calculations into two classes: the energy barriers for the breaking of a Si-H bond and the energy barriers for hopping between interstitial sites. The interstitial hopping barriers were calculated with initial and final states taken from the µ H (z) calculations, such that the H atom moves in the z direction by approximately one Si-Si bond length, away from the interface. The bond-breaking initial and final states were determined by an algorithm that scanned the heterojunction stacks to find all mono-hydride Si-H bonds. We identified the mono-hydride Si-H bonds by the coupled criteria that the hydrogen was separated from a single nearest Si by a bond length of about 1.50Å, in agreement with established bond length values [40,57,73,74], while the Si was not bonded to a second hydrogen. It was also confirmed that the hydrogen was not at a bond center by checking for multiple Si neighbors within the cutoff distance of 1.7Å. The mono-hydride Si-H were chosen as initial configurations for the NEB analysis. For the final configurations, the hydrogen was moved from this mono-hydride bond to a next nearest neighbor bond-centered position. To remain consistent with established Si-H-Si configurations, the two Si atoms were simultaneously displaced by 0.45Å along the direction of the bond. The energies of the initial and final states were then minimized using the FIRE minimizer to an energy tolerance of 10 −3 eV before commencing the NEB analysis. After developing and establishing the validity of this method we proceeded to compute the energy barriers for 657 initial-final state pairs for the breaking of Si-H bonds, and for more than 2000 initial-final state pairs for interstitial hopping. We determined the energy barriers for both the forward and backward direction of these processes.

Electronic Properties of Defects
We completed the characterization of these stacks by determining some of their electronic properties. We decided to measure the inverse participation ratio (IPR) and the partial charge density in each stack. In short, the IPR charac-terizes how localized an electron wavefunction is [14], for a details on the IPR see Supplementary Equations [16][17][18]. The ideal case of an IPR=1 indicates complete localization of the electron wavefunction on a single atom. Complete delocalization is indicated by the IPR value of the order of ( 1 N ), in our case, 1/406. To adopt a pragmatic convention, we took IPR values exceeding 0.1 as an indicator that the electronic wavefunction was substantially localized on one or just a few atoms. We have described our IPR method in great detail in a recent paper. [14] Using this IPR method, we determined the change in the number of localized defects when a hydrogen atom broke away from a Si-H bond and moved to an interstitial Si-Si bond center.
We computed 100 such processes and found that when a hydrogen broke away from a Si-H bond at the interface, on average it induced 0.94 localized dangling bonds. This justifies us identifying a hydrogen breaking away from a Si-H bond at the interface as a defect-generating event.
The statistics of these defect-generating events were found to be surprisingly broad numerically and spatially. Several bond-breaking events generated 2-4 localized dangling bonds. Also, the IPR value changed not only on the center Si, but also on atoms several interatomic spacings away. [75] Both of these facts are consistent with the softness and glassy interconnected dynamics of a-Si. Fig. 11(a) shows the IP R nkj s, calculated for all Kohn-Sham orbitals obtained by DFT as a function of their energy for a typical c-Si/a-Si:H stack.
For completeness, we also calculated the partial charge of these left-behind dangling bonds. We found that their partial charge was -0.1e ± 0.1e. As such, it is appropriate to approximately identify these dangling bonds as neutral defects. This brings the theoretical simulation into agreement with the experimental analysis that explained the degradation data in terms of neutral defects. Partial charges were obtained using the Lowdin population analysis, where plane wave wavefunctions are projected on to the atomic orbitals of atoms.

DATA AVAILABILITY
The data that support the findings of this study are available from the authors on reasonable request, see author contributions for specific data sets.
The GAP suite of programs is freely available for non-commercial use from http://github.com/libatoms/GAP. The Quantum Espresso software package is freely available from www.quantum-espresso.org.
The LAMMPS software package is freely available from lammps.sandia.gov.