Also at ]Center for Nanoscale Materials, Argonne National Laboratory, Lemont,
IL 60439, USA

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

Andrew Diggs<EMAIL_ADDRESS>Zitong Zhao<EMAIL_ADDRESS>Davis Unruh
[  Gergely T. Zimányi  Physics Department, University of California
Davis, CA 95616, USA  Reza Vatan Meidanshahi  Salman Manzoor  Mariana Bertoni
Stephen M. Goodnick  School of Electrical, Computer and Energy Engineering,
Arizona State University
Tempe, AZ 85287, USA

###### Abstract

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.

silicon, heterojunction, degradation, passivation, hydrogen, bilayer, machine
learning, microstructure

## 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 $\rm{V_{oc}}$. [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. [16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 15, 30, 31, 32, 12]

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 further theoretical
and numerical research to provide additional insight. The most effective
numerical methods are first principles electronic-structure studies using
density functional theory (DFT). [39, 40, 24, 41, 28] However, the
computational expense of DFT scales unfavorably with O(N3). 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 $\rm{V_{oc}}$ degradation rate by
as much as 80%. Such a substantial reduction may considerably help the market
acceptance of Si HJ cells.

## Results and Discussion

### 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 $(\tau_{\rm{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:

$\frac{1}{\tau_{\rm{eff}}}=\frac{1}{\tau_{aug}}+\frac{1}{\tau_{rad}}+\frac{1}{\tau_{SHR}}+\frac{2S}{W}$
(1)

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/\tau_{\rm{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 $\Delta n$ and a range of temperatures $T$ to construct S$(\Delta
n,T)$. Finally, we determined the defect density $N$ and the charge density
$Q$ at the interface by fitting S$(\Delta n,T)$ to the amphoteric model
proposed by Olibet et al.. [32] We collected the S$(\Delta 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)}$.

Figure 1: Degradation of the surface passivation of silicon heterojunction
solar cells. The quality of the surface passivation of four silicon
heterojunction cells with varying thickness was tracked over the course of a
year by monitoring the effective minority carrier lifetime $\tau_{\rm{eff}}$
(blue circles). 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 $\rm{N}(t)$ (red diamonds).

Fig. 1 shows that the effective carrier lifetime $\tau_{\rm{eff}}(t)$
decreased over the year. Fitting the surface component of $\tau_{\rm{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 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.

Figure 2: Silicon density and hydrogen chemical potential at the c-Si/a-Si:H
interface. The hydrogen chemical potential $\mu_{H}(z)$ (blue), and the Si
density as indicated by the first peak height of the radial distribution
function $\rm{RDF}1(z)$ (red) across the interface of the Si HJ are shown.
$\mu_{H}(z)$ tracks the Si density remarkably closely, making is plausible
that the Si density gradient is the cause of the H chemical potential
gradient.

Finally, we determined $\mu_{H}(z)$, the position-dependent 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 $\mu_{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.

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 initial-final 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, $\rm{E_{bond}}$=1.31$\pm$0.03 eV, $\rm{E_{inter}}$=0.49$\pm$0.08
eV, were consistent with previous experimental and computational values:
$\rm{E_{bond}}$=1.2-1.4 eV, $\rm{E_{inter}}$= 0.48-0.50 eV,[18, 23, 40, 25,
26] $\rm{E_{inter}}$ being the average of $\rm{E_{drift}}$ and
$\rm{E_{reverse}}$.

Figure 3: Energy barrier distributions that control hydrogen dynamics. The
distribution, mean value and standard deviation of the energy barriers are
shown for the four key processes that control the activated thermal motion of
hydrogen at the c-Si/a-Si:H interface: a, Freeing a hydrogen atom by breaking
a Si-H bond at interface; b, Recapture of the freed H atom back to the Si-H
bond; c, Activated drift of H atom away from the interface; d, Reverse drift
toward interface.

In an important distinction, however, most previous studies reported only
single values or narrow distributions ($\sigma\approx 0.025~{}\rm{eV}$) for
these barriers. In contrast, we found that the width of these distributions
were quite large, with $\sigma$ 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: $\rm{E_{drift}}$
controlling the drift away from the interface into the a-Si:H layer, and
$\rm{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
$\rm{E_{inter}}=0.50~{}\rm{eV}$. [40, 26, 18, 23] However, our first main
result, the existence of the hydrogen chemical potential gradient, resulted in
a bias between these barriers: $\rm{E_{drift}}=0.41~{}\rm{eV}$, and
$\rm{E_{reverse}}=0.58~{}\rm{eV}$. While the average of these two barriers,
$\rm{E_{inter}}=0.49~{}\rm{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 $\rm{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 $\rm{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$,
$\frac{\rm{E_{drift}}-\rm{E_{reverse}}}{k_{B}T}\approx 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 $\rm{E_{bond}}$, the bond-recapture barrier $\rm{E_{recap}}$,
and the barrier against forward drift into the a-Si:H layer $\rm{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:

$\frac{\partial}{\partial t}\mathbf{H}=\mathbf{M}\mathbf{H}$ (2)

where the $\mathbf{H}$ vector denotes the hydrogen densities in the three
model states, and

$\mathbf{M}=\begin{pmatrix}-k_{1}&k_{2}&0\\\ k_{1}&-(k_{2}+k_{3})&0\\\
0&k_{3}&0\end{pmatrix}$ (3)

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:

$H_{i}(t)=H_{i}(0)e^{-\frac{\alpha t}{2}}\bigg{(}\rm{cosh}\big{(}\frac{\beta
t}{2}\big{)}+\frac{\alpha-2k_{1}}{\beta}\rm{sinh}\big{(}\frac{\beta
t}{2}\big{)}\bigg{)}$ (4)

where where $\alpha\equiv k_{1}+k_{2}+k_{3}$, and
$\beta\equiv\sqrt{\alpha^{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 $\Delta N_{i}(t)=-\Delta 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 SolDeg-simulated $N(t)$ alongside the
experimentally determined $N(t)$ from Fig. 1.

Figure 4: Defect density $N(t)$ determined by SolDeg. The blue line shows our
SolDeg-simulated defect density as determined by the expectation value of the
solution of the three barrier model using the computed barrier distributions
in Fig. 3. The experimentally measured $N(t)$ is shown with the turquoise
triangles. The SolDeg-simulated blue line is not a fit: only the overall scale
factor was selected to match the measured data. Also shown: stretched
exponential fit to the experimental data, with red dashed line.

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
quantitatively reliable description of the hydrogen-induced 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/\tau$, causing it to decrease
over time. Further, following the models of Green [58] and Olibet et al. [32],
or Blank et al. [59], $\rm{V_{oc}}$ is another quantity directly effected by
the defect density $N(t)$ via

$\rm{V_{oc}}(t)=C-\frac{k_{\rm{B}}T}{q}\rm{ln}(N(t)),$ (5)

where $C$ is a constant and $q$ is the elementary charge, for a full
derivation see Supplementary Equations [4-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 slab-merging protocol to create a Si density
that was lower at the interface and increased into the a-Si:H layer.

Figure 5: Hydrogen chemical potential for a reverse gradient Si HJ interface.
The reversed-gradient stacks were constructed so that the Si density (red), as
described by the RDF1(z), exhibited a minimum at the interface. The Si density
minimum induced a corresponding minimum in $\mu_{H}(z)$ (blue). The reversal
of the $\mu_{H}(z)$ gradient induced a force on the mobile H atoms toward the
interface.

Fig. 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 $\rm{E_{Drift}}$ barrier
controlling the forward drift increased as the reversed Si density gradient
increased. Fig. 6 shows the remarkable result of our SolDeg simulations of the
NoDeg cells.

Figure 6: Defect growth over time with different $\mu_{H}(z)$ gradients.
SolDeg-simulated results of the interface defect density $N(t)$ with mean
values of $\rm{E_{Drift}}$ representing a $\mu_{H}(z)$ gradient pointing away
from the interface - blue; or a reversed gradient pointing towards the
interface: moderate gradient - red, stronger gradient - green, comparable to
the magnitude of the original gradient (blue). The defect density $N(t)$ at 1
year was suppressed by 80% in reversed gradient 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, $\rm{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 $\rm{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 $\rm{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 $\rm{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 $\rm{V_{oc}}$. [62, 63,
64, 37]

The above groups did not perform systematic degradation studies. However, we
speculate that the reported $\rm{V_{oc}}$ improvements could have been the
consequence of suppressed degradation that took place from the moment of
fabrication to the time of the $\rm{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
simulations 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 $\rm{V_{oc}}$
degradation rate may be reduced from 0.5% per yr to 0.1% per yr.

## Methods

### 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\Omega cm$ resistivity, and
thickness of $\approx 275\mu m$. The wafers were rigorously cleaned in Piranha
$(\rm{H_{2}SO_{4}}\colon\rm{H_{2}O_{2}}|4\colon 1)$ and RCA-b
$(\rm{H_{2}O}\colon\rm{HCl}\colon\rm{H_{2}O_{2}}|6\colon 1\colon 1)$
solutions. An important note on crystal orientation: deWolf et al.. showed
that on long time scales the time dependence of recombination lifetime is
independent of the orientation of the crystal axes. Therefore, we believe our
results are relevant for HJ cells which use c-Si wafers with (111) orientation
at the interface.

Wafers were formed with four different thicknesses W in the range of
$160-260\mu$ by etching them over different times in an HNA mixture
$(\rm{HF}\colon\rm{HNO_{3}}\colon\rm{CH_{3}COOH}|10\colon 73\colon 17)$,
ending with a dip in a buffered oxide $(\rm{HF}\colon\rm{H_{2}O}|10\colon 1)$
etch solution, for specific wafer cleaning and etching protocols see
Supplementary Methods 1. Subsequently, 50 nm of a-Si:(i) was deposited on both
sides of the wafers using plasma-enhanced chemical vapor deposition (PECVD)
technique in Octopus I tool from INDEOtec S.A. The effective minority carrier
lifetime was measured with the WCT-120TS tool from Sinton Instruments,
lowering the temperature from $230^{o}$ C to $30^{o}$ C.

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.

Figure 7: Si-H pair correlation function in liquid Si:H. The liquid Si-H RDF
computed with the Si-H GAP and with DFT are barely distinguishable, while both
differ from the Tersoff-computed RDF markedly. These results further
demonstrate the quality of ML based potentials such as the Si-H GAP.

Finally, we demonstrated the remarkable benefit of having constructed the Si-H
GAP by creating a hydrogenated amorphous slab of 4,096 Si atoms and 558 H
atoms with our Si-H GAP, see Fig. 8. Such simulations are utterly inaccessible
by DFT calculations. The Si-H GAP makes it possible to create exceptionally
large samples, and thus to determine quantities inaccessible for smaller DFT
simulations, an example of which is the statistical distribution of void
sizes.

Figure 8: a-Si:H supercell, with 4096 Si atoms and 558 H atoms. Simulations of
this size are utterly inaccessible by DFT calculations. The Si-H GAP makes it
possible to create exceptionally large samples, and thus to determine
quantities that would be missed by smaller DFT simulations, an example of
which is the statistical distribution of void sizes.

### 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 hydrogen-doped 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 constant-volume 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.

Figure 9: c-Si/a-Si:H interface structure. Shown is a stick and ball model of
one of the final interface structures used in our study. The blue atoms
represent silicon and the red atoms represent hydrogen.

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.

Figure 10: Density matching at the left interface. The z-dependent density of
a stack 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. The shaded regions show the 2.5 Å cuboid representing
(blue) the top of c-Si, and (yellow) the bottom of the a-Si:H. The peak height
in both regions shows the smooth transition of the Si density across the
interface, similar to the desired structure of the transition zone passivating
layer.

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 $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 $\mu_{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 $\mu_{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 $\mu_{H}(x,y,z)$ energies were smoothed by averaging
over all insertion locations within a 2.5 $\rm{\AA}$ 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
$\mu_{H}(z)$. Carrying out this procedure for the 60 interfaces (2 interfaces
per stack), our comprehensive $\mu_{H}(z)$ is the result of about 25,000
calculations. Fig. 2 shows the so-calculated $\mu_{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 $\rm{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 $\rm{RDF}(r,z)$ as a function of the lateral
radius r, which we denoted by $\rm{RDF}1(z)$. A high $\rm{RDF}1(z)$ indicated
a high degree of crystallinity and therefore high density, which translated to
a low porosity, whereas a low $\rm{RDF}1(z)$ indicated the opposite. The red
line in Fig. 2 shows $\rm{RDF}1(z)$ as a function of z. $\rm{RDF}1(z)$ visibly
also exhibits a marked gradient across the interface at $z=0$. The correlation
between the gradient of the hydrogen chemical potential $\mu_{H}(z)$ and the
gradient of $\rm{RDF}1(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 $\approx 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 hydrogen 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 inter-replica forces, $\mathbf{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.[29, 36, 18, 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 $\mu_{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\rm{\AA}$, in agreement with established bond length values [40, 73, 74,
57], 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\rm{\AA}$. 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\rm{\AA}$ 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}\rm{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
characterizes how localized an electron wavefunction is [14], for a details on
the IPR see Supplementary Equations [16-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
($\rm{\frac{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 $IPR_{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.

Figure 11: Defect characterization at the c-Si/a-Si:H interface. IPR of atoms
at the c-Si/a-Si:H interface and in aSi:H regions a, before and b, after a
hydrogen atom is taken out of a dangling bond and placed at a Si-Si bond
center.

For completeness, we also calculated the partial charge of these left-behind
dangling bonds. We found that their partial charge was -0.1e $\pm$ 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.

## Code availability

The Si-H GAP can be found at https://github.com/dgunruh/Si-H-GAP and
https://libatoms.github.io/GAP/data.html.

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.

## References

  * Yoshikawa  _et al._ [2017] K. Yoshikawa, H. Kawasaki, W. Yoshida, T. Irie, K. Konishi, K. Nakano, T. Uto, D. Adachi, M. Kanematsu, H. Uzu, and K. Yamamoto, Silicon heterojunction solar cell with interdigitated back contacts for a photoconversion efficiency over 26%, Nature Energy 2, 10.1038/nenergy.2017.32 (2017).
  * Jordan  _et al._ [2018] D. C. Jordan, C. Deline, S. Johnston, S. R. Rummel, B. Sekulic, P. Hacke, S. R. Kurtz, K. O. Davis, E. J. Schneller, X. Sun, M. Alam, and R. Sinton, Silicon heterojunction system field performance, IEEE J. Photovolt 8, 177 (2018).
  * Gao  _et al._ [2022] B. Gao, Y. Shao, W. Liu, H. Xiang, Y. Yu, and Z. Liu, Out-door reliability and degradation of hit, cigs, n-type multi-busbar, perc, and cdte modules in shanghai, china, Solar Energy Materials and Solar Cells 236, 111490 (2022).
  * Staebler and Wronski [1977] D. L. Staebler and C. R. Wronski, Reversible conductivity changes in discharge‐produced amorphous si, Appl. Phys. Lett. 31, 292 (1977).
  * Johnson  _et al._ [1991] N. Johnson, C. Doland, F. Ponce, J. Walker, and G. Anderson, Hydrogen in crystalline semiconductors: A review of experimental results, Physica B: Condensed Matter 170, 3 (1991).
  * Shimizu [2004] T. Shimizu, Staebler-wronski effect in hydrogenated amorphous silicon and related alloy films, Japanese Journal of Applied Physics 43, 3257 (2004).
  * Wronski and Niu [2014] C. R. Wronski and X. Niu, The limited relevance of swe dangling bonds to degradation in high-quality a-si:h solar cells, IEEE Journal of Photovoltaics 4, 778 (2014).
  * Melskens  _et al._ [2014] J. Melskens, M. Schouten, A. Mannheim, A. S. Vullers, Y. Mohammadian, S. W. H. Eijt, H. Schut, T. Matsui, M. Zeman, and A. H. M. Smets, The nature and the kinetics of light-induced defect creation in hydrogenated amorphous silicon films and solar cells, IEEE Journal of Photovoltaics 4, 1331 (2014).
  * Guha  _et al._ [1992] S. Guha, J. Yang, S. J. Jones, Y. Chen, and D. L. Williamson, Effect of microvoids on initial and light‐degraded efficiencies of hydrogenated amorphous silicon alloy solar cells, Applied Physics Letters 61, 1444 (1992), https://doi.org/10.1063/1.107564 .
  * Stutzmann [1985] M. Stutzmann, Role of mechanical stress in the light‐induced degradation of hydrogenated amorphous silicon, Applied Physics Letters 47, 21 (1985), https://doi.org/10.1063/1.96415 .
  * Redfield and Bube [1989] D. Redfield and R. H. Bube, Reinterpretation of degradation kinetics of amorphous silicon, Applied Physics Letters 54, 1037 (1989), https://doi.org/10.1063/1.100790 .
  * Bernardini and Bertoni [2019] S. Bernardini and M. I. Bertoni, Insights into the degradation of amorphous silicon passivation layer for heterojunction solar cells, phys. status solidi (a) 216, 1800705 (2019).
  * Holovský  _et al._ [2020] J. Holovský, S. Martín De Nicolás, S. De Wolf, and C. Ballif, Amorphous/crystalline silicon interface stability: Correlation between infrared spectroscopy and electronic passivation properties, Adv. Mater. Interfaces 7, 2000957 (2020).
  * Unruh  _et al._ [2021] D. Unruh, R. V. Meidanshahi, C. Hansen, S. Manzoor, S. M. Goodnick, M. I. Bertoni, and G. T. Zimanyi, From femtoseconds to gigaseconds: The soldeg platform for the performance degradation analysis of silicon heterojunction solar cells, ACS Appl. Mater. Interfaces 13, 32424–32434 (2021).
  * De Wolf  _et al._ [2012] S. De Wolf, C. Ballif, and M. Kondo, Kinetics of $a$-si:h bulk defect and $a$-si:h/$c$-si interface-state reduction, Phys. Rev. B 85, 113302 (2012).
  * Beyer [1991] W. Beyer, Hydrogen effusion: a probe for surface desorption and diffusion, Physica B: Condensed Matter 170, 105 (1991).
  * Beyer [1996] W. Beyer, Physica B (1996).
  * Santos  _et al._ [1993] P. V. Santos, N. M. Johnson, R. A. Street, M. Hack, R. Thompson, and C. C. Tsai, Hydrogen migration and electronic carriers in a-si:h, Phys. Rev. B 47, 10244 (1993).
  * Santos  _et al._ [1992] P. Santos, N. Johnson, and R. Street, Enhanced hydrogen diffusion under illumination in hydrogenated amorphous silicon, MRS Proceedings 258, 353 (1992).
  * Street  _et al._ [1987] R. A. Street, C. C. Tsai, J. Kakalios, and W. B. Jackson, Hydrogen diffusion in amorphous silicon, Philosophical Magazine B 56, 305 (1987), https://doi.org/10.1080/13642818708221319 .
  * Van de Walle [1996] C. G. Van de Walle, Stretched-exponential relaxation modeled without invoking statistical distributions, Phys. Rev. B 53, 11292 (1996).
  * Kakalios  _et al._ [1987a] J. Kakalios, R. A. Street, and W. B. Jackson, Stretched-exponential relaxation arising from dispersive diffusion of hydrogen in amorphous silicon, Phys. Rev. Lett. 59, 1037 (1987a).
  * Biswas  _et al._ [1998] R. Biswas, Q. Li, B. C. Pan, and Y. Yoon, Mechanism for hydrogen diffusion in amorphous silicon, Phys. Rev. B 57, 2253 (1998).
  * Van de Walle and Street [1995] C. G. Van de Walle and R. A. Street, Silicon-hydrogen bonding and hydrogen diffusion in amorphous silicon, Phys. Rev. B 51, 10615 (1995).
  * Van De Walle and Street [1995] C. G. Van De Walle and R. A. Street, Silicon-hydrogen bonding and hydrogen diffusion in amorphous silicon, MRS Proceedings 377, 389 (1995).
  * van Wieringen and Warmholtz [1960] A. van Wieringen and N. Warmholtz, Physics, 22, 849 (1956); rc frank and je thomas, jr, J. Phys. Chem. Solids 16, 144 (1960).
  * Carlson and Wronski [1976] D. E. Carlson and C. R. Wronski, Amorphous silicon solar cell, Appl. Phys. Lett 28, 671 (1976), https://doi.org/10.1063/1.88617 .
  * Wronski [1984] C. Wronski, Chapter 10 the staebler-wronski effect, in  _Hydrogenated Amorphous Silicon_, Semiconductors and Semimetals, Vol. 21, edited by J. I. Pankove (Elsevier, 1984) pp. 347–374.
  * Herring  _et al._ [2001a] C. Herring, N. M. Johnson, and C. G. Van de Walle, Energy levels of isolated interstitial hydrogen in silicon, Phys. Rev. B 64, 125209 (2001a).
  * Burrows  _et al._ [2008] M. Z. Burrows, U. K. Das, R. L. Opila, S. De Wolf, and R. W. Birkmire, Role of hydrogen bonding environment in a-si:h films for c-si surface passivation, Journal of Vacuum Science & Technology A 26, 683 (2008), https://doi.org/10.1116/1.2897929 .
  * De Wolf and Kondo [2007] S. De Wolf and M. Kondo, Abruptness of a-si:h/c-si interface revealed by carrier lifetime measurements, Applied Physics Letters 90, 042111 (2007), https://doi.org/10.1063/1.2432297 .
  * Olibet  _et al._ [2007] S. Olibet, E. Vallat-Sauvain, and C. Ballif, Model for a-si:h/c-si interface recombination based on the amphoteric nature of silicon dangling bonds, Phys. Rev. B 76, 035326 (2007).
  * Cattin  _et al._ [2021] J. Cattin, L.-L. Senaud, J. Haschke, B. Paviet-Salomon, M. Despeisse, C. Ballif, and M. Boccard, Influence of light soaking on silicon heterojunction solar cells with various architectures, IEEE Journal of Photovoltaics 11, 575 (2021).
  * Veirman  _et al._ [2022] J. Veirman, A. J. K. Leoga, L. Basset, W. Favre, O. Bonino, A. L. Priol, N. Rochat, and D. Rouchon, Understanding the improvement of silicon heterojunction solar cells under light soaking, AIP Conference Proceedings 2487, 020017 (2022), https://aip.scitation.org/doi/pdf/10.1063/5.0089284 .
  * Olibet  _et al._ [2006] S. Olibet, E. Vallat, and C. Ballif, Effect of light induced degradation on passivating properties of a-si:h layers deposited on crystalline si, (2006).
  * De Wolf  _et al._ [2008] S. De Wolf, S. Olibet, and C. Ballif, Stretched-exponential a-si:h/c-si interface recombination decay, Applied Physics Letters 93, 032101 (2008), https://doi.org/10.1063/1.2956668 .
  * Liu  _et al._ [2016] W. Liu, L. Zhang, R. Chen, F. Meng, W. Guo, J. Bao, and Z. Liu, Underdense a-si:h film capped by a dense film as the passivation layer of a silicon heterojunction solar cell, Journal of Applied Physics 120, 175301 (2016), https://doi.org/10.1063/1.4966941 .
  * Sai  _et al._ [2018] H. Sai, P.-W. Chen, H.-J. Hsu, T. Matsui, S. Nunomura, and K. Matsubara, Impact of intrinsic amorphous silicon bilayers in silicon heterojunction solar cells, Journal of Applied Physics 124, 103102 (2018), https://doi.org/10.1063/1.5045155 .
  * Geerlings  _et al._ [2003] P. Geerlings, F. De Proft, and W. Langenaeker, Conceptual density functional theory, Chem. Rev. 103, 1793 (2003).
  * Herring  _et al._ [2001b] C. Herring, N. M. Johnson, and C. G. Van de Walle, Energy levels of isolated interstitial hydrogen in silicon, Phys. Rev. B 64, 125209 (2001b).
  * Wagner and Grossman [2008] L. K. Wagner and J. C. Grossman, Microscopic description of light induced defects in amorphous silicon solar cells, Phys. Rev. Lett. 101, 265501 (2008).
  * Finnis [2003] M. Finnis,  _Interatomic Forces in Condensed Matter_ (Oxford University Press, 2003).
  * Bartók  _et al._ [2010] A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons, Phys. Rev. Lett. 104, 136403 (2010).
  * Bartók  _et al._ [2013] A. P. Bartók, R. Kondor, and G. Csányi, On representing chemical environments, Phys. Rev. B 87, 184115 (2013).
  * Bartók  _et al._ [2018] A. P. Bartók, J. Kermode, N. Bernstein, and G. Csányi, Machine learning a general-purpose interatomic potential for silicon, Phys. Rev. X 8, 041048 (2018).
  * Thompson  _et al._ [2015] A. Thompson, L. Swiler, C. Trott, S. Foiles, and G. Tucker, Spectral neighbor analysis method for automated generation of quantum-accurate interatomic potentials, J. Comput. Phys 285, 316 (2015).
  * Chen  _et al._ [2017] C. Chen, Z. Deng, R. Tran, H. Tang, I.-H. Chu, and S. P. Ong, Accurate force field for molybdenum by machine learning large materials data, Phys. Rev. Mater. 1, 043603 (2017).
  * Behler and Parrinello [2007] J. Behler and M. Parrinello, Generalized neural-network representation of high-dimensional potential-energy surfaces, Phys. Rev. Lett. 98, 146401 (2007).
  * Behler [2011] J. Behler, Neural network potential-energy surfaces in chemistry: a tool for large-scale simulations, Phys. Chem. Chem. Phys. 13, 17930 (2011).
  * Behler [2016] J. Behler, Perspective: Machine learning potentials for atomistic simulations, J. Chem. Phys. 145, 170901 (2016), https://doi.org/10.1063/1.4966192 .
  * Kocer  _et al._ [2019] E. Kocer, J. K. Mason, and H. Erturk, A novel approach to describe chemical environments in high-dimensional neural network potentials, J. Chem. Phys 150, 154102 (2019), https://doi.org/10.1063/1.5086167 .
  * Shapeev [2016] A. V. Shapeev, Moment tensor potentials: A class of systematically improvable interatomic potentials, Multiscale Model. Simul. 14, 1153 (2016).
  * Jinnouchi  _et al._ [2019] R. Jinnouchi, F. Karsai, and G. Kresse, On-the-fly machine learning force field generation: Application to melting points, Phys. Rev. B 100, 014105 (2019).
  * Chmiela  _et al._ [2017] S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Schütt, and K.-R. Müller, Machine learning of accurate energy-conserving molecular force fields, Sci. Adv. 3, e1603015 (2017).
  * Smith  _et al._ [2017] J. S. Smith, O. Isayev, and A. E. Roitberg, Ani-1: an extensible neural network potential with dft accuracy at force field computational cost, Chem. Sci. 8, 3192 (2017).
  * Drautz [2019] R. Drautz, Atomic cluster expansion for accurate and transferable interatomic potentials, Phys. Rev. B 99, 014104 (2019).
  * Unruh  _et al._ [2022] D. Unruh, R. V. Meidanshahi, S. M. Goodnick, G. Csányi, and G. T. Zimányi, Gaussian approximation potential for amorphous Si:H, Phys. Rev. Materials 6, 065603 (2022).
  * Green [2006] M. A. Green,  _Third Generation Photovoltaics_ (Springer Berlin, Heidelberg, 2006).
  * Blank  _et al._ [2017] B. Blank, T. Kirchartz, S. Lany, and U. Rau, Selection metric for photovoltaic materials screening based on detailed-balance analysis, Phys. Rev. Applied 8, 024032 (2017).
  * Lee  _et al._ [2014] K.-S. Lee, C. B. Yeon, S. J. Yun, K. H. Jung, and J. W. Lim, Improved surface passivation using dual-layered a-si:h for silicon heterojunction solar cells, ECS Solid State Letters 3, P33 (2014).
  * Zhao  _et al._ [2022] Y. Zhao, P. Procel, A. Smets, L. Mazzarella, C. Han, L. Cao, G. Yang, Z. Yao, A. Weeber, M. Zeman, and O. Isabella, Effects of (i)a-si:h deposition temperature on passivation quality and performance of high-efficiency silicon heterojunction solar cells, in  _Proc. 49th IEEE PVSC_ , Photovoltaic Materials and Devices Group, Delft University of Technology (IEEE, Philadelphia, PA, 2022).
  * Duan  _et al._ [2022] W. Duan, A. Lambertz, K. Bittkau, D. Qiu, K. Qiu, U. Rau, and K. Ding, A route towards high-efficiency silicon heterojunction solar cells, Progress in Photovoltaics: Research and Applications 30, 384 (2022).
  * Ru  _et al._ [2020] X. Ru, M. Qu, J. Wang, T. Ruan, M. Yang, F. Peng, W. Long, K. Zheng, H. Yan, and X. Xu, 25.11% efficiency silicon heterojunction solar cell with low deposition rate intrinsic amorphous silicon buffer layers, Solar Energy Materials and Solar Cells 215, 110643 (2020).
  * Zhang  _et al._ [2017] Y. Zhang, C. Yu, M. Yang, L.-R. Zhang, Y.-C. He, J.-Y. Zhang, X.-X. Xu, Y.-Z. Zhang, X.-M. Song, and H. Yan, Significant improvement of passivation performance by two-step preparation of amorphous silicon passivation layers in silicon heterojunction solar cells, Chinese Physics Letters 34, 038101 (2017).
  * Ishimaru  _et al._ [1997] M. Ishimaru, S. Munetoh, and T. Motooka, Generation of amorphous silicon structures by rapid quenching: A molecular-dynamics study, Phys. Rev. B 56, 15133 (1997).
  * Štich  _et al._ [1991] I. Štich, R. Car, and M. Parrinello, Amorphous silicon studied by ab initio molecular dynamics: Preparation, structure, and properties, Phys. Rev. B 44, 11092 (1991).
  * Jarolimek  _et al._ [2009] K. Jarolimek, R. A. de Groot, G. A. de Wijs, and M. Zeman, First-principles study of hydrogenated amorphous silicon, Phys. Rev. B 79, 155206 (2009).
  * Henkelman and JÃ³nsson [2000] G. Henkelman and H. JÃ³nsson, Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points, The Journal of Chemical Physics 113, 9978 (2000).
  * Henkelman and Jónsson [2000] G. Henkelman and H. Jónsson, Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points, The Journal of Chemical Physics 113, 9978 (2000).
  * Bitzek  _et al._ [2006] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Structural relaxation made simple, Phys. Rev. Lett. 97, 170201 (2006).
  * Kakalios  _et al._ [1987b] J. Kakalios, R. A. Street, and W. B. Jackson, Stretched-exponential relaxation arising from dispersive diffusion of hydrogen in amorphous silicon, Phys. Rev. Lett. 59, 1037 (1987b).
  * Shlesinger and Montroll [1984] M. F. Shlesinger and E. W. Montroll, On the williams & watts function of dielectric relaxation, Proceedings of the National Academy of Sciences 81, 1280 (1984).
  * Panzarini and Colombo [1994] G. Panzarini and L. Colombo, Hydrogen diffusion in silicon from tight-binding molecular dynamics, Phys. Rev. Lett. 73, 1636 (1994).
  * Tersoff [1989] J. Tersoff, Modeling solid-state chemistry: Interatomic potentials for multicomponent systems, Phys. Rev. B 39, 5566 (1989).
  * Meidanshahi  _et al._ [2021] R. V. Meidanshahi, D. Vasileska, and S. M. Goodnick, Role of hydrogen in the electronic properties of a-si:h/c-si heterostructures, The Journal of Physical Chemistry C 125, 13050 (2021), https://doi.org/10.1021/acs.jpcc.1c03173 .

## Acknowledgements

We acknowledge useful discussions with Christophe Balliff and Chase Hansen.
This material is based upon work supported by the U.S. Department of Energy’s
Office of Energy Efficiency and Renewable Energy (EERE) under the Solar Energy
Technologies Office Award Numbers DE-EE0008979 and DE-EE0009835. The views
expressed herein do not necessarily represent the views of the U.S. Department
of Energy or the United States Government. Work performed at the Center for
Nanoscale Materials, a U.S. Department of Energy Office of Science User
Facility, was supported by the U.S. DOE, Office of Basic Energy Sciences,
under Contract No. DE-AC02-06CH11357. This research used resources of the
National Energy Research Scientific Computing Center (NERSC), a U.S.
Department of Energy Office of Science User Facility located at Lawrence
Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231
using NERSC award BES-ERCAP m3634.

## Author contributions

AD developed the three barrier model, analyzed and fitted the experimental
data, and played an energetic role in writing the paper. ZZ created the slabs,
fused them into stacks, invented the hydrogen chemical potential analysis, and
analyzed the number and charge of the generated defects; his contribution is
entirely comparable to AD, and thus both are equally deserving to be first
authors. RVM performed the stack-structure optimizations with DFT methods and
the IPR analysis. DU developed the sample creation protocols and the Kinetic
Monte Carlo simulations. SM performed the degradation experiments. MB guided
and oversaw the degradation experiments. SMG contributed to the development of
the physical picture and the analysis and provided useful guidance throughout.
GTZ conceived the SolDeg method and the physical analysis, directed the entire
project, and played a primary role in writing the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

The Supplementary Information is not included in the arXiv version of this
article.