Monte Carlo simulation of solar cells

Solar Cells, I0 (1983) 211 - 222

211

MONTE CARLO SIMULATION OF SOLAR CELLS

C. M A N F R E D O T T I and M. MELIGA

Istituto di Fisica Superiore, Corso Massimo D'Azeglio 46, 10125 Turin (Italy) (Received December 8, 1982; accepted May 10, 1983)

Summary The Monte Carlo method was applied for the first time to calculate the performance of a solar cell. It was proved that the method is self-consistent even for a relatively low number of histories and, in the simple case of a Schottky barrier, the results are in good agreement with the existing simple theory. The possible application of this method is relatively straightforward and could give much more information than existing numerical methods.

1. Introduction Even though solar cells are semiconductor devices that are easy to produce, a complete analysis which takes into account the morphology of the cell and the influence of various parameters is rather complex. Generally it is necessary, for a complete understanding of the problem, to solve Poisson's equation and the continuity equation simultaneously for electrons and holes. Under illumination, these are non-linear partial differential equations. A problem of this kind may be solved essentially in two ways: by numerical analysis and by the Monte Carlo method (MCM). Numerical analysis, which has many examples in the literature [1 - 3], presents the advantage of short computation times, but it cannot avoid problems due to boundary conditions, discontinuities, complex models or structures. The MCM, which differs from numerical methods by its simplicity and versatility, can also overcome very easily the most complicated problems connected with structural complexity, electrical potential behaviour, discontinuities etc. without increasing the program complexity and, by taking the appropriate precautions, without even a marked increase in computation time. In fact this time can be kept to the same order of magnitude as for numerical analysis; moreover, it is not proportional to the complexity of the problem. Finally, the MCM, being generally based on a microscopic treatment, leads to an easy visualization of the physical processes occurring in solar cells and for this reason it can give various results without modification of the program itself. 0379-6787/83/$3.00

© Elsevier Seauoia/Print~rl i,, ~

~,^,L . . . . . . .

212

In our work, a Schottky barrier n-type Au-Si solar cell was chosen because of both structural simplicity and lack of literature of previous work in this field. However, the program is quite general and with slight modifications it may be applied to other kinds of solar cells, such as p - n junction cells, metal-insulator-semiconductor (MIS) cells or heterojunction cells, and it can easily be used for more complex structures, e . g . tandem cells.

2. The m e t h o d The MCM is very extensively used in various fields of physics, whenever processes or experiments have to be simulated. The best advantages result in cases where a few types of process are repeated several times (cooling down of neutrons, electromagnetic or nucleonic cascades etc.). In these cases each particle is followed until it has lost all its energy and the various processes t h a t occur are interpreted on a probability basis. Since the sum of the probabilities of the various processes at each step is of course equal to unity, a random number between zero and unity is consequently chosen in order to decide the kind o f process which should occur. In our case, for instance, in order to select a p h o t o n from the solar spectrum, the integral of its energy distribution is normalized to u n i t y and the random number chosen by the program is compared with the normalized partial integrals of the energy distribution. For the subsequent process, light absorption in the cell, the problem can be solved analytically since the probability of absorption between x and x + dx is given by

d/ p(x ) dx

= - - = o~

exp(--ax) dx

I0 which is obviously normalized to u n i t y , and thus the random number R is given by the following expression: d

f exp(--ax) dx = 1 -- exp(--~d) = 1 - - R 0

Therefore the penetration depth of the involved p h o t o n can be obtained from the equation In R d----

~(~,)

a(~) being the absorption coefficient at the selected p h o t o n wavelength.

3. The model In the present work, the program describes an n-type Au-Si Schottky barrier solar cell with a surface area of 1 cm 2. However, the various physical

213

TABLE 1 Symbols and parameters used in the Monte Carlo method Parameter

Symbol

Value

Metal thickness Solar cell thickness Barrier height {lowered by Schottky effect} Energy gap (as a function of temperature) Resistivity Electron lifetime Hole lifetime Hole diffusion length

SM SL qbB

100 A 300 p m 0.78 eV

Eg(T)

1.2

p Te

1 ~cm 3.95 x 10 -4 s 2.62 x 10-Ss 150 p m

Tp Lp

--

(3

X

10-4T) eV

0.780.

0.778

! antiref/~ coating mE

0. 774

uJ

r

barrier

base

0.770

tal 0.766

W $L

20

40 d(A}

60

8'0

Fig. 1. Geometry of the solar cell with the various zones considered by the program. The figure is not to scale. Fig. 2. The metal-semiconductor potential barrier used in the program as a function of depth. The Schottky effect is included and the lowest energy corresponds to a tunnelling transmission of 3%.

and geometrical parameters are given as external data and can be varied according to the characteristics o f the device under investigation. In Table 1 the parameters which were used in this work are given. From a geometrical point of view, the m o d e l structure may be divided into four different parts (Fig. 1). In the a n t i r e f l e c t i o n c o a t i n g a n d m e t a l e l e c t r o d e z o n e only the transmission of incident light is taken into account [4]. In fact, for typical metal thicknesses of a S c h o t t k y barrier solar cell the photoexcitation of electrons from the metal to the semiconductor is negligible. The barrier z o n e thickness is evaluated in the program [5]. The shape of the potential barrier is shown in Fig. 2: it was obtained b y considering

214 TABLE 2 Variation in Iasc and Voc as functions of the time interval used in the calculation of the current At

Number o f boxes

Iasc (mA)

Voc (V)

(1/30)Tp (1/100)Tp (1/150)Tp (1/200)rp (1/250)Tp

8 14 17 20 22

24.28 26.52 27.38 27.63 27.86

0.281 0.283 0.284 0.284 0.284

a barrier height o f 0.8 eV and determining the lowering due to the Schottky effect. In any case, it is given as an external vector and therefore it is possible to consider barriers of any shape and height. The barrier zone is divided into boxes in such a way as to obtain a constant electrical field in each box. In the base zone, carrier diffusion takes place and therefore the characteristic parameter will be the minority {hole) diffusion length Lp. The base zone is also divided into boxes and the box width depends on the time interval At chosen to calculate the instantaneous current. The diffusion length {given in terms of the number o f boxes) is determined by simulating a semi-infinite crystal with the same physical characteristics as the base and by considering hole injection at each time interval (see Table 2) and thus the n u m b e r of boxes necessary for the base zone is obtained. The back electrode is important because it gives rise to p h o t o n reflection and it is characterized by a reflectivity coefficient that depends on the metal of the electrode and varies with the p h o t o n wavelength.

4. The program The flow chart of the program is shown in Fig. 3. The MCM, if rigorously applied to our program, would take a lot o f central processing unit (CPU) time. In fact, under air mass {AM) 1 conditions about 3 × 101~ photons should impinge in 1 s on a cell of area 1 cm 2. Four of five carrier lifetimes are enough for the semiconductor to reach a stationary state. In fact it was proved t h a t the cell reaches stationary conditions in a time shorter than the lifetime and therefore the calculation is limited to an interval equal to one lifetime. A further reduction in the n u m b e r of photons can be obtained by considering a n u m b e r o f histories sufficient to give a good representation of the solar spectrum and o f the absorption coefficient. Usually, 10 4 histories are enough from this point o f view and therefore a multiplication coefficient is applied to the n u m b e r o f carriers reaching the barrier zone. This point will be discussed further below. The p h o t o n source considered exhibits an AM 1 solar spectrum which, for simplicity, is restricted to a wavelength interval

215

®

Fig. 3. Flow chart of the program (for nomenclature see Appendix A). TABLE 3 Comparison of the total energy and the average energy of photons in the real and in the simulated source

Data for the real source Data for the simulated source Total number ofphotons Totalpower density (mW cm-2) Average energy o f photons (eV)

2.49 × 1017 cm-2 s-1 72.62 1.8235

104 histories 72.6171 1.8230

between 3 0 0 0 and 12 000 A, corresponding t o a t ot al intensity o f 72.62 mW cm -2. Table 3 shows t ha t t h e simulated source (104 histories) gives a t ot al intensity and an average p h o t o n e n e r g y almost coincident with those o f t he real source. Th e dark c u r r e n t is simulated b y considering t h e pot ent i al barrier illust ra t ed in Fig. 2, taking f o r electrons in t he m et al a density o f states with a normal distribution and including t h e e f f e c t o f t h e F e r m i - D i r a c statistics. Of course, only thermionic emission o f electrons in t h e metal is taken into account. Th e energy interval f o r electrons is restricted in o r d e r t o neglect, on t h e high energy side (at a b o u t 1 eV above t he Fermi level), t o o low an

216

o.sj/ O 0.3

0.1 0.I

O.2 AE (eV)

0/3

Fig. 4. Q u a n t u m m e c h a n i c a l reflection coefficient Q as a f u n c t i o n o f t h e e n e r g y differe n c e between electrons a n d t h e t o p o f t h e barrier.

occupation probability and, on the low energy side, t o o low a tunnelling probability (approximately 3%). F o r each history, both the tunnelling probability and, for E > ~B, the quantum mechanical reflection coefficient [6] are taken into account (Fig. 4). The dark current is then calculated as the drift current of electrons injected into the semiconductor. The dark electric field is evaluated as a derivative of the electrical potential, and the carrier mobility is calculated as a function of the doping and of the electric field according to ref. 7. The simulation of the shortcircuit current Isc is obtained with 104 histories; for each history, the wavelength of an incident p h o t o n and its path length d = --(ln R)/o~(~) in the cell (where ~(~) is the absorption coefficient and R is a random number) are determined. If d ~< w t h e p h o t o n will be absorbed in the depletion layer (the q u a n t u m efficiency is taken as equal to unity); if w < d <~ SL the p h o t o n will be absorbed in the base zone; finally, if d > SL the p h o t o n will be reflected b y the back electrode, according to the relative reflectivity coefficient. After each event, the charge parameters C(I) and G(I), which represent the photogenerated carrier numbers in box I o f the depletion zone and o f the base zone respectively, are updated. In the calculation of the stationary short-circuit current, the diffusion process from each b o x of the base is followed at each time interval At, which is expressed as a fraction o f the minority carrier lifetime. This parameter is quite important since it implicitly determines, as previously shown, the number of boxes in the base zone. Subsequently, after a time lag equal to the lifetime rp the drift current and diffusion current are calculated, their sum being Isc. As shown in Fig. 5 (where At = (1/30)rp) I,c increases rapidly with time, reaching for t = Tp a value within 2% o f the stationary value. Moreover, it turns o u t that the average current Ia,~ calculated in the interval 0 - Tp is also within 3% o f the stationary value. The rapid increase in Is¢ can be under-

217

15: E

l{:t

I~J/30

20/30 t(Vp}

30/30

Fig. 5. Short-circuit current transient calculated at time intervals At = (I/30)~p.

stood if it is considered that the current transient is affected not only by the minority carrier lifetime but also by the drift time in the depletion layer, which is m u c h shorter than the lifetime. Even the diffusion time in the base zone is practically "weighted" according to the exponential distribution of the photogenerated carriers. As far as At is concerned, it is shown in Table 2 that it is not necessary to choose the shortest At since, for example by changing from At = (1/100)rp to At = (1/150)rp, only a 3 % increase in the average Isc value is obtained. Neverthless, in our calculation At = (1/250)rp was used. The present version of the program calculates the short-circuit current, the dark current, the open~ircuit voltage, the current-voltage (I-V) characteristics, the fillfactor FF and the conversion efficiency; the C P U time for 104 histories is about 5 rain. Moreover, the program also provides the possibility of evaluating all these data at various temperatures.

5. Results and discussion Since in the model the number o f histories is necessarily low, some checks are needed in order to confirm that this n u m b e r is sufficient to account for the various physical processes involved with good precision. Among the possible tests, the solar spectrum and the absorption coefficient at different wavelengths have been considered. The minority carrier diffusion (since it is simulated in a semi-infinite crystal) and the carrier transport in the depletion layer are clearly less d e p e n d e n t on the total n u m b e r o f events.

218 !

0.084 !

~' I

'I ii i I-iii Iii

I

i-i I

i

J

0.034

I

--If

- -

'I--

-] ........ I

~I

I-~ I

]

l

i

I I I

-T

E-

I-

I

-z

II

I.~ J'i-

~,

-i~

I

-r

I r

I

-,

(

I

II - i I

II I {

I

3' I 1

T[

0.3

0.7

I~

Fig. 6. Comparison between the real and simulated solar spectrum in the wavelength range 0.3 - 1.2 pro. O n the y axis the ratio of the total n u m b e r of photons in the real source to the total n u m b e r of photons in the simulated source is reported. The spectra have been shifted for comparison.

TABLE 4 Comparison between the simulated and the real absorption coefficient at three wavelengths and for two thicknesses

Thickness (~m)

1 10

X (A)

Numberof absorbed photons

~sinmla~d

~m~

4000 6000 8000 4000 6000 8000

9876 3429 952 10000 9850 6320

4.39 x 104 4.199 x 103 1.0004 x 103 -4.1997 × 103 0.99967 × 103

4.4 x 104 4.2 x 103 1.0 × 103 4.2 x 103 1.0 x 10 3

Figure 6 shows t hat the solar spectrum in t h e interval 0.3 - 1.2 pm is very well r e p r o d u c e d by simulation. It must be r e m e m b e r e d (see Table 3) t h a t t h e to tal solar energy simulated is, within 0.004%, coi nci dent with t h a t delivered to t h e cell b y taking into a c c o u n t t h e real n u m b e r o f phot ons. T h e a b s o r p t i o n coefficient ~(k) is also simulated at various wavelengths; it is derived f r o m t h e decrease in t h e n u m b e r o f p h o t o n s due t o t h e absorpt i o n in the various boxes and is c o m p a r e d with t he real absorption coefficient. As shown in Table 4, when t h e absorption in a crystal either 1 or 10 /~m thick is considered t h e m a x i m u m possible error in ~(k) is 0.2%. Two examples o f simulated p h o t o n a t t e n u a t i o n at 4 0 0 0 • for a thickness o f 1 ~m and at 8000 A f o r a thickness o f 10 #m are shown in Fig. 7

219

i+

|[-Jr-*

Nixie le

w

l-++ I~+.. [+. I-]...

N m

..~ 4, Z

"o. + lee

[--.

_.= ~i

J--]--

50+

10

30

50

10

30 50 ~x *box Fig. 7. Attenuation of 4000 A photons in a thickness of 1 p m (logarithmic scale). Fig. 8. Attenuation of 8000 A photons in a thickness of 10 pm (linear scale). TABLE 5 Comparison between the solar cell performances obtained by simulation and by theory Theoretical results a

MCM results With Schottky effect b

Without Schottky effect Average values Stationary values

Isc (A) I s (A) Voc (V) FF 7) (%)

27.86 x 10 -3 5.0 x 10 - 7 0.284 0.71 7.695

27.86 x 10 -3 2.2 x 10 - 7 0.304 0.72 8.39

28.67 x 10 -3 2.2 x 10 - 7 0.304 0.72 8.66

28.84 x 10 -3 2.253 × 10 - 7 0.304 0.72 8.69

aThe theory does not include the Schottky effect. The values given are stationary values. b Average values.

and Fig. 8 respectively: a truly exponential photon attenuation is clearly evident. The results obtained in the present work are summarized in Table 5, where I~, Is, Voc, FF and W are reported and compared with theoretical results obtained by using equations quoted in the literature [8]. In order to make the comparison more understandable, the results of the simulation are also reported for a barrier without the lowering due to the Schottky effect, since theoretical results can be obtained for this case only. For completeness, results are also presented for both the average and the stationary values of the parameters. In the theoretical calculation, an average value of the quantum mechanical reflection coefficient has been used [6]. The agreement between the Monte Carlo simulation and theoretical calculation is extremely good, the larger differences being 2% for I, (stationary values) and 3% for Isc and ~ (average values).

220 V (v} 0.1

0.2

0.3

J

10

/ 1

5 20

30l Fig. 9. I-V characteristics in A M 1 conditions as obtained by the M C M (o) (average values) and by the theory (4). The theory does not include the Schottky effect.

0.6

~ 0.4 0.2

0

O.5

0.7

O. 9

1.1

k('/~m)

Fig. 10. Responaivity R ( m A m W - 1 ) as a function of wavelength as obtained by the MCM ( e ) and by the theory (4) (reflection from the back electrode was not included). The theory does not take into account the Schottky effect.

The comparison was extended to the I - V characteristics (Fig. 9) and to the responsivity (Figs. 10 and 11). In Fig. 9 the differences can be attributed to the Schottky effect which lowers Voe. The differences which appear in Fig. 10 can be explained by looking at Fig. 11, which shows the drift current in the depletion layer and the diffusion current in the base zone. In the former case the disagreement is a m a x i m u m at 4 0 0 0 A, corresponding to a penetration depth comparable with the depletion layer thickness, while in the latter case the disagreement corresponds to a deep penetration of light in the base zone (about 30 pm). It is therefore probable that the former difference is due to the lower surface electric field in the MCM and the latter results because the diffusion current is calculated by the theory under the assumption o f complete depletion at the end o f the base zone.

221

~

0.3

0.5

i

0.7

o

n

0.9

1.1

Fig. 11. The drift current and the diffusion current as functions of wavelength as calculated by the MCM (o) and by the theory (A). Monochromatic light w i t h a p o w e r of 10 mW was incident on the cell.

6. Conclusions Our work m a y be considered as a first example of the possibility o f using the MCM in the simulation o f solar cells. It was proved that the MCM can work with good precision, even taking into account a reasonable n u m b e r of p h o t o n s and considering the diffusion process in microscopic detail. Since reliable or usable experimental results are not available, our results were compared with those o f the usual theory. The good agreement which was obtained can be considered as p r o o f o f the validity o f the MCM in a very simple physical case. However, as has been noted in the introduction, the MCM can easily be extended to more complex structures such as tandem cells, heterojunction cells, graded band gap cells and MIS cells and to transient effects in solar cells. This extension will be carried o u t in future work, particularly if more experimental data are available. From our point o f view, the MCM is probably the simplest m e t h o d by which to s t u d y the influence of the various parameters (lifetime, temperature, diffusion profiles, thickness of the various zones, absorption coefficient etc.) on the performance o f the cell, particularly in cases where these parameters are also d e p e n d e n t variables. Efforts will be made in future work to investigate all these possibilities o f the MCM. We are also confident that a two- or three
222 3 J. Michel and A. Mircea, Acta Electron., 18 (4) (1975) 331. 4 M. V. Schneider, Bell Syst. Tech. J., 45 (1966) 1611. 5 S. M. Sze, Physics o f Semiconductor Devices, Wiley-Interscience, New York, 1960, Chap. 8, p. 363. 6 C. R. Crowell and S. M. Sze, J. Appl. Phys., 37 (7) (1966) 2683. 7 S. Sharfetter and A. Gummel, IEEE Trans. Electron Devices, 16 (1969) 64. 8 H. J. Hovel, in R. K. Willardson and A. C. Beer (eds.), Semiconductors and Semimetals, Vol. 11, Solar Cells, Academic Press, New York, 1975, pp. 112 - 126.

Appendix A: nomenclature for Fig. 3

E(I) electric field FF Iasc Id /dr Iisc Is Isc I- V

fill factor average current diffusion current drift current instantaneous current dark current short-circuit current I - V characteristic Nbb number of boxes in the base Nbw number of boxes in the barrier Ndt number of time intervals Np number of time intervals NB number of histories V(I) potential Voc open,circuit voltage w barrier width At time interval efficiency Pe electron mobility #h hole mobility