AFORS-HET, A NUMERICAL PC-PROGRAM FOR SIMULATION OF HETEROJUNCTION SOLAR CELLS,

VERSION 1.1 (OPEN-SOURCE ON DEMAND), TO BE DISTRIBUTED FOR PUBLIC USE

R.Stangl* (1), A.Froitzheim (2), M.Kriegel (1),

T.Brammer (3), S.Kirste (1), L.Elstner (1), H.Stiebig (3), M.Schmidt(1), W.Fuhs (1)

(1) Hahn-Meitner-Institut Berlin (HMI), Abteilung Silizium Photovoltaik, Kekulé-Str. 5, D-12489 Berlin, Germany
(2) Shell-Solar GmbH, Otto-Hahn Ring 6, D-81739 München, Germany
(3) Photovoltaik-Institut, Forschungszentrum Jülich, D-52452 Jülich, Germany
*Corresponding author: e-mail: stangl@hmi.de, tel: +49/30/8062-1312, fax: +49/30/8062-1333

ABSTRACT: We offer the (new open-source on demand) Version 1.1 of AFORS-HET, a numerical computer simulation
program for modeling (thin film) heterojunction solar cells. It will be distributed free of charge on CD-ROM at the
conference site and can also be downloaded via internet: www.hmi.de/bereiche/SE/SE1/projects/aSicSi/AFORS-HET
An arbitrary sequence of semiconducting layers can be modeled, specifying the corresponding layer and interface
properties, i.e. the defect distribution of states (DOS). Using Shockley-Read-Hall recombination statistics, the one-
dimensional semiconductor equations are solved (1) for thermodynamic equilibrium, (2) for steady-state conditions
under an external bias and/or illumination, (3) for small additional sinusoidal perturbations of the applied
bias/illumination. Thus, the internal cell characteristics, such as band diagrams, local generation and recombination rates,
local cell currents, carrier densities and phase shifts can be calculated. Furthermore, a variety of characterization methods
can be simulated, such as current-voltage (I-V), internal and external quantum efficiency (IQE, EQE), intensity and
voltage dependant surface photovoltage (SPV), photo- and electroluminescence (PEL), impedance spectroscopy (IMP),
capacitance-voltage (C-V) and capacitance-temperature (C-T). A user-friendly graphical interface allows not only a
visualisation of the simulated data (all data can be imported and exported), but also arbitrary parameter variations can be
performed.
Version 1.1 has now been rewritten in a modulized form, such that new measurement methods and new numerical
modules can be implemented by external users (open-source on demand). Most measurement methods have been
improved and new ones have been added. Several numerical modules have been implemented. The optical generation
rate can now be calculated taking into account coherent/incoherent multiple reflections. If modelling crystalline silicon,
ionized impurity scattering and carrier-carrier scattering can be considered. The program will be described and
demonstrated showing selected results on the simulation of amorphous/crystalline silicon heterojunction solar cells.

Keywords: Simulation, Thin Film, Experimental Methods

1 INTRODUCTION hydrogenated silicon are deposited on top of a thick

(300 µm) crystalline c-Si(p) wafer.
In order to investigate (thin film) heterojunction solar
cells, a lot of different electrical measurement methods AFORS-HET is distributed free of charge on CD-
are used, ranging from standard solar cell ROM and can also be downloaded via internet:
characterisation techniques like current-voltage (I-V) or www.hmi.de/bereiche/SE/SE1/projects/aSicSi/
quantum efficiency (EQE, IQE) to more advanced AFORS-HET. The new version 1.1 of AFORS-HET,
characterization techniques like for example surface which will be launched, is now presented in a modulized
photovoltage (SPV), photo- or electroluminescence (PL, form, such that new measurements and new numerical
EL), capacitance-voltage (C-V), capacitance-temperature modules can be implemented also by external users
(C-T), impedance (IMP), or intensity modulated (open-source on demand).
photocurrent spectroscopy (IMPS).
We therefore developed a numerical simulation tool DOS at intefaces
(AFORS-HET, automat for simulation of hetero- (and bulk layers) Parameter Measurement
structures), which allows to simulate the output of Variations Methods
different measurement techniques for an arbitrary
sequence of semiconducting layers and interfaces, with
an arbitrary number of defects distributed within the
different band gaps. A user friendly interface allows to
perform multidimensional parameter variations and to
visualise and analyse the corresponding results see
(Fig.1). Different numerical modules allow to treat
different experimental situations, like the choice of a MS or MIS
metal/semiconductor or a metal/insulator/semiconductor front contact
front contact. steady-state &
We use the program mainly to simulate amorphous/ sinusoidal modulations
crystalline silicon heterojunction solar cells, of the
structure TCO/a-Si:H(n)/a-Si:H(i)/c-Si(p)/a-Si:H(p+)/Al, Figure 1: graphical interface of AFORS-HET.
where some ultrathin layers (5 nm) of amorphous
2 MODELLING CAPABILITIES The electron/hole currents jn / p are driven by the
gradient of the corresponding quasi Fermi energy EFn,p.
AFORS-HET numerically solves the one Within the bulk of a semiconducting layer, this is
dimensional semiconductor equations with the equivalent to the sum of a diffusion and a drift current
appropriate boundary conditions under steady-state with the corresponding mobility µn,p.
∂EFn (x, t ) µ kT ∂n(x, t ) ∂ϕ (x, t )
conditions and under additional small sinusoidal
perturbations. The set of coupled partial differential jn (x, t ) = qµ n n(x, t ) =− n + µ n n ( x, t )
∂x q ∂x ∂x
equations is transformed into a set of nonlinear algebraic
∂E Fp (x, t ) µ p kT ∂p(x, t ) ∂ϕ (x, t )
j p (x, t ) = qµ p p (x, t ) − µ p p (x, t )
equations by the method of finite differences. So far, the =−
grid on which the equations are solved is fixed at the ∂x q ∂x ∂x
beginning of the calculation (fixed x-discretization, non- The only time dependence which is treated in this model
adaptive meshing), but can be modified by the user, if is a small periodic sinusoidal perturbation of a steady-
needed. The free electron density ni, the free hole density state condition. That is, neglecting second order terms,
pi and the cell potential ϕi at each gridpoint are used as all time dependant quantities can be described by a small
independent variables. All other variables in the complex amplitude, i.e.
discretized differential equations / boundary conditions n(x, t ) = n(x ) + n~ (x ) eiω t p ( x, t ) = p ( x ) + ~
p ( x ) e iω t
are expressed in a way that they only depend on these
independent variables. The resulting nonlinear equations ϕ (x, t ) = ϕ (x ) + ϕ~ (x ) eiω t
are solved using the Newton-Raphson iteration scheme
thereby requiring a good starting solution. If equilibrium 2.1.1 Generation
conditions are chosen, the program supplies a starting If the heterostructure is illuminated (specifying the
solution from analytical approximations, otherwise the spectral distribution Φ(λ ) of the incoming photon flux),
last calculated solution serves as a starting solution for the super-bandgap optical generation rate Gn = G p (for
the new boundary conditions to be solved. Alternatively,
hc / λ ≥ E g ) from the valence band into the conduction
starting solutions can also be saved and loaded. For more
information on the numeric procedure see [1]. band of the semiconductor layers can be obtained in
With the ability to solve for the internal cell different ways: (1) Assuming simple Lambert-Beer
characteristics (band diagrams, local generation and absorption by specifying the spectral absorption
recombination rates, local cell currents, carrier densities coefficient α of each layer. (2) By taking
and phase shifts) under some specified boundary coherent/incoherent internal multiple reflections into
conditions, measurement methods can be defined by a account, specifying the dielectric properties (n,k) of each
specific variation of the external boundary conditions and layer. (3) Alternatively, the generation rate can also be
some additional post-processing data analysis. imported, using external programs for its calculation. An
In the following, the differential equations and optical sub-bandgap generation (for hc / λ < E g ) from a
corresponding boundary conditions, which are solved by defect to the conduction/valence band can be defined by
AFORS-HET under various circumstances, are stated.
specifying optical emission coefficients eno (E ), e op (E ) ≠ 0
New numerical modules, solving for modified
differential equations / boundary conditions may be for the defect state:
added by external users (open source on demand). eno (E , x ) = σ no N C Φ (λ , x ) ϑ (EC − E − hc λ )
e op (E , x ) = σ op NV Φ(λ , x ) ϑ (E − EV − hc λ )
2.1 Bulk
Poisson’s equation and the transport equation for ( σ no, p : optical capture cross sections, Φ ( λ ,x) spectral
electrons and holes are solved in one dimension. photon flux, with wavelength λ at the position x
ε 0ε r ∂ 2ϕ (x, t )
= p (x , t ) − n (x , t ) + N D − N A + ∑ ρ t (x , t ) N C ,V , EC ,V : conduction/valence band density/energy,
q ∂x 2
ϑ : ϑ (E ) = 1 for E ≤ 0 , ϑ (E ) = 0 for E > 0 )
defects

1 ∂jn (x, t ) ∂
− = Gn (x, t ) − Rn (x, t ) − n(x, t )
q ∂x ∂t 2.1.2 Recombination
1 ∂j p (x, t ) ∂ Recombination from the conduction band into the
= G p ( x, t ) − R p ( x, t ) − p ( x , t )
q ∂x ∂t valence band may occur directly (band to band
The electron density n, the hole density p, and the recombination, Auger recombination) and via trap states
(Shockley-Read-Hall recombination, SHR):
electric potential ϕ are the independent variables for
which the system of differential equations is solved. (q: Rn , p (x, t ) = Rndirect
, p (x, t ) + Rn , p (x, t )
SHR

electron charge, ε 0 , ε r the absolute/relative dielectric Direct recombination requires to specify the band to band
constant). N D / A are the concentrations of the rate constant r BB and the Auger rate constants rnA , rpA .

, p (x, t ) = Rn , p (x ) + Rn , p (x ) e
~ direct iω t
donors/acceptors, which are assumed to be completely Rndirect direct

ionized. The charge stored in the defects is described by

a distribution function f t , specifying the probability that , p (x ) = r
Rndirect [
BB
]{
+ rnA n(x ) + rpA p (x ) n(x ) p (x ) − N C NV e
− Eg kT
}
defects with defect density N t at the position E within
~
, p (x ) =
Rndirect [
r BB + rnA n(x ) + 2rpA n(x ) p(x ) ~
2
p (x ) ]
the bandgap are occupied with electrons: [
+ r + r p ( x ) + 2 r n (x ) p (x ) n (x )
BB A
p
2 A
n
~ ]
acceptor-type defect: ρ t (x, t ) = − ∫ dE f t (E , x, t ) N t (E ) An arbitrary number of defects distributed arbitrarily
donator-type defect: ρt (x, t ) = ∫ dE (1 − f t (E , x, t )) N t (E ) within the bandgap of the semiconductor layers can be
specified. SHR recombination requires to specify the
energetic distribution of the defect density N t (E ) of each
ε 0ε r−
∂ϕ
∂x
− ε 0ε r+
∂ϕ
∂x
=q ∑ ρit ( ) ( )
ϕ xit− = ϕ xit+
defect and its electron/hole capture coefficients x it- x it+ defects

cn , p = vn , pσ n , p ( vn, p : thermal velocity, σ n, p : capture cross ∆EC ∆EC

section). The emission rates en, p are then given by ( )

jn xit− = vn− n xit− e ( ) −
kT
ϑ (−∆EC )
− vn+ n xit+ e ( ) −
kT
ϑ (∆EC )
+ Rn xit− ( )
∆EC ∆EC
en (E , x ) = cn N c e −(Ec − E ) kT + eno (E , x ) ( )
jn xit+ = vn+ n xit+ e ( ) −
kT
ϑ (−∆EC )
− vn+ n xit+ e ( ) −
kT
ϑ (∆EC )
− Rn xit+ ( )
e p (E , x ) = c p NV e −( E − EV ) kT
+ e (E , x )
o
∆EV ∆EV

( ) ( ) − ϑ (−∆EV )
( ) − ϑ (∆EV )
( )
p
j p xit− = v −p n xit− e kT − v +p n xit+ e kT + R p xit−
and the SRH recombination rate due to the defects is
, p (x, t ) = Rn , p (x ) + Rn , p (x ) e
~ SRH iω t ∆EV ∆EV
RnSRH
( ) ( ) ϑ (−∆EV )
( ) ϑ (∆EV )
( )
SRH
− −
j p xit+ = v +p n xit− e kT
− v +p n xit+ e kT
− R p xit+
R SRH
n (x ) = ∫ dE {cn n(x ) Nt (E ) (1 − f t (E , x )) − en (E, x ) N t (E ) ft (E, x )} The recombination at the interface is treated
R p (x ) = ∫ dE {c p p (x ) N t (E ) f t (E , x ) − e p (E , x )N t (E ) (1 − f t (E , x ))}
SRH equivalently to the bulk recombination, with two
exceptions: (1) the interface states are given in cm-2
RnSRH (x ) = ∫ dE { cn N t (E ) (1 − f t (E , x )) n~ (x )
~
instead of cm-3, (2) the distribution function changes, as
− cn N t (E ) f t (E , x ) − en (E , x ) N t (E ) f t (E , x )}
~ ~ the interface states can interact with both adjacent
semiconductors:
R p (x ) = ∫ dE { c p N t (E ) f t (E , x ) ~ p (x )
~ SRH
f it (E ) =
+ c p N t (E ) f t (E , x ) + e p (E , x ) N t (E ) f t (E , x ) }
~ ~
( )
cn+ n xit+ + e +p (E ) + cn− n xit− + e −p (E ) ( )
The distribution function f t (E , x, t ) describing the cn+ n ( )
xit+ + en+ (E ) + c +p p(xit+ ) + e +p (E ) + cn− n(xit− ) + en− (E ) + c −p p(xit− ) + e −p (E )
f it (E ) =
occupied defect states is given by SHR theory: ~

f t (E , x, t ) = f t (E , x ) + f t (E , x ) eiωt ( ) ( ) ( ) ( )
~
cn+ {1 − f it (E )} n~ xit+ − c +p f it (E ) ~
p xit+ + cn− {1 − f it (E )} n~ xit− − c −p f it (E ) ~
p xit−

f t (E , x ) =
cn n(x ) + e p (E , x ) cn+ n ( )
xit+ + en+ (E ) + c +p p ( )
xit+ + e +p (E ) + cn− n ( ) (E ) + c p(x )+ e (E ) + iω
xit− + en− −
p
−
it
−
p
cn n(x ) + en (E , x ) + c p p(x ) + e p (E , x )
cn {1 − f t (E , x )} n~ (x ) − c p f t (E , x ) ~p (x ) 2.3 Boundaries
f t (E , x ) =
~
cn n(x ) + en (E , x ) + c p p (x ) + e p (E , x ) + iω The electric potential is fixed to zero at one contact
(for example the back contact). At the second contact a
boundary condition has to be specified, which relates the
2.2 Interfaces external cell voltage/current to internal quantities.
Furthermore, the electron and hole currents into the metal
2.2.1 Interface currents driven by drift diffusion contacts have to be modeled.
This models uses an interface layer of a certain
thickness (which can be specified), in which the material 2.3.1 Schottky contact, voltage controlled
properties change linearly from semiconductor I to ϕ (0 ) = φ front − φback + Vext ϕ (L ) = 0
semiconductor II. The specified interface defects are
distributed homogeneously within this layer. The jn (0) = qS n
front
(n(0) − n (0)) eq jn (L ) = −qS back
n (n(L) − n (L)) eq
transport across the heterojunction interface can then be
treated like in the bulk of a semiconducting layer
j p (0) = −qS pfront ( p(0) − p (0)) eq j p (L ) = qS back
p ( p(L) − p (L )) eq

(drift/diffusion driven). The cell currents have to be ( φ : metal work function, S : surface recombination
evaluated directly from the gradient of the corresponding velocity). The majority carrier density under equilibrium
quasi fermi energies, since the electron affinity χ , the ( neq or peq ) is given by the energy barrier of the
bandgap E g and the effective conduction/valence band metal/semiconductor contact, the corresponding minority
density of states N C , NV will depend on the position x carrier density ( peq or neq ) by the mass action law
qφ −qχ E g − qφ + qχ Eg
− − −
within the interface layer. neq = N c e kT
or peq = NV e kT
, neq peq = N C NV e 2 kT
n (x )
E Fn (x ) = −qχ (x ) + qϕ (x ) − kT ln
N C (x )
2.3.2 Schottky contact, current controlled
p (x )
E Fp (x ) = −qχ (x ) + qϕ (x ) − E g (x ) + kT ln The boundary condition for the external cell voltage
NV ( x ) is replaced by a condition for the external cell current:
jn (0 ) + j p (0 ) = jext
2.2.2 Interface currents driven by thermionic emission
The transport across a heterojunction interface can 2.3.3 MIS contact, voltage controlled
alternatively be modeled by thermionic emission over the In case of a metal/insulator/semiconductor front
energetic barrier of the interface, which is the
contact, the insulator capacity C I has to be specified. At
conduction/valence band offset ∆EC ,V (with χ −, + , E g−, + :
x=0 there are additional interface states, which are
electron affinity and bandgap of the semiconductor left or treated equivalent to the bulk, with the exception that the
right to the interface): interface states are given in cm-2 instead of cm-3.
∆EC = qχ + − qχ − , ∆EV = E g+ − E g− + qχ + − qχ − ∂ϕ
C I {Vext + ϕ (0 ) − ϕ (L )} − ε 0ε r+ =q ∑ ρit ϕ (L ) = 0
The interface is treated as a boundary condition. Thus, ∂x x =0 + defects

conditions for the potential and for the electron/hole jn (0 ) = 0 jn (L ) = 0

currents at both sides of the interface have to be stated
j p (0 ) = 0 j p (L ) = 0
( xit−, + : position infinitesimal left or right to the interface).
3 GRAPHICAL INTERFACE illumination), (2) a change of the external cell voltage
(variation of dc-voltage, ac-amplitude and frequency)
A user-friendly graphical interface allows an easy and (3) a change of the device temperature is possible.
visualization of all simulated quantities. All data (internal
cell results and measurement methods) can be imported 4.3 Simulation of measurement methods
and exported and thus graphically compared. The Measurements can be simulated, just by varying
heterostructures and some parts of it (layers, interfaces) external parameters like in a real experiment, and by
can be saved and loaded. The cell dicretization can be performing some post-processing data analysis. So far
modified by the user and also saved and loaded. The current-voltage (I-V), quantum efficiency (QE),
same holds for starting solutions used for solving the impedance (IMP), capacitance-voltage (C-V),
differential equations. The user can pre-select the capacitance-temperature (C-T), surface photovoltage
measurement methods he wants to use. Arbitrary (SPV), photo- and electroluminescence (PEL) and
parameter variations can be performed, selecting an electron beam induced current (EBIC) is implemented.
arbitrary number of external and internal input parameter New measurement methods may be added by external
of the program to be varied in a specified range. As an users (open source on demand).
output of a parameter variation, all possible internal cell As an example capacitance-temperature (C-T)
results as well as measurements can be selected. simulations are shown, see Fig.3. The capacitance/
Parameter variations can also be stored and reloaded. conductance at a given frequency ω and an external
~
applied ac-voltage Vext can be calculated from the
4 SELECTED RESULTS ~
resulting ac-current jext through the simulated device:

In the following section, we like to give a brief C (ω ) =

~
~
{ }
1 Im jext (ω )
, G (ω ) =
{ }
Re jext (ω )
~
~
ω Vext Vext
overview of the capabilities of the program and highlight
-4
its new features. So far, we used the program mainly to 1x10
D 0 -1 0
10
cm
-2
-5 in t
simulate amorphous/crystalline silicon heterojuction 1x10
D 10
12
cm
-2
in t
solar cells in order to study the influence of internal cell -6
Specific Conductance/ω (F/cm )

1x10
2

12 -2
D in t
4x10 cm
parameters on measurements and cell efficiency [2-9]. 1x10
-7
D in t
2x10
13
cm
-2

-8
1x10
4.1 Optical calculations 1x10
-9

In order to demonstrate the new optical capabilities 1x10

-1 0

of the program, Figure 2 shows a comparison with the 1x10

-1 1

commercial software window-coating-designer. The -1 2

1x10
reflection and transmission of the mulitlayer stack -1 3
1x10
air/ZnO(80nm)/c-Si(1µm)/air is plotted. Both values are -1 4
1x10
in perfect agreement with the results from window- 100 150 200 250 300
coating designer. The minimum of R between 500 nm T e m p e ra tu re (K )
and 600 nm is caused by the TCO with a thickness of Figure 3: C-T curves at 100 Hz for a-Si:H(n)/c-Si(p)
80 nm. The oscillating R and T for a wavelength above heterostructures simulated with different a-Si:H(n)/
500 nm are interference fringes. c-Si(p) interface state densities, Dint [7].
1.0

0.8
air TCO Si air 7 SUMMARY
80nm 1µm
The capabilities of the simulation program AFORS-
0.6
HET have been described. The new open-source on
T
R, T

demand Version 1.1 will allow external users to

0.4 R implement their own measurement methods and numeric
T
R modules. It will be distributed free of charge.
0.2

8 REFERENCES
0.0 [1] A.Froitzheim, HMI, Dissertation 2003.
400 600 800 1000
[2] Stangl, Froitzheim, Elstner, Fuhs, Proc. 17th Eur. PV Sol.
Wavelength (nm)
En. Conf., München, Germany 2001, 1387
Figure 2: Reflection (R) and transmission (T) of the [3] Froitzheim, Stangl, Elstner Schmidt, Fuhs,
Proc. 25th IEEE Conf., New Orleans, 2002
simple multilayerstack air/ZnO(80nm)/Si(1µm)/air.
[4] Stangl, Froitzheim, Fuhs.,
Comparison of the optical calulations: symbols: AFORS- Proc. PV in Europe Conf., Rome, Italy, 2002, 123-126
HET, lines: window-coating designer. [5] Froitzheim, Stangl, Kriegel, Elstner, Fuhs, Proc. WCPEC-3,
World Conf. PV En. Conv., Osaka, Japan, 2003, 1P-D3-34
4.2 Internal cell characteristics [6] Stangl, Froitzheim, Schmidt, Fuhs, Proc. WCPEC-3,
The internal cell characteristics, such as band World Conf. PV En. Conv., Osaka, Japan, 2003, 4P-A8-45
diagrams, local generation and recombination rates, local [7] Gudovskikh, Kleider, Stangl, Schmidt, Fuhs,
cell currents, carrier densities and phase shifts can be this conference, 2CV.1.21
[8] Stangl, Schaffarzik, Laades, Kliefoth, Schmidt, Fuhs,
calculated under various external boundary conditions:
this conference, 2CV.1.18
I.e., (1) a change of the illumination (variation of [9] Schmidt, Korte, Kliefoth, Schoepke, Stangl, Laades,
wavelength and intensity, choosing an arbitrary spectral Brendel, Scherff, this conference, 2DO.2.5
illumination file and an additional monochromatic