TWO-DIMENSIONAL OPTICAL MODEL FOR SIMULATING PERIODIC OPTICAL STRUCTURES IN THIN-FILM
SOLAR CELLS
Andrej Čampa, Janez Krč and Marko Topič University of Ljubljana, Faculty of Electrical Engineering, Ljubljana, Slovenia
Key words: optical modelling, thin-film solar cells, periodic structures, light management
Abstract: Novel optical concepts based on periodic structures such as diffractive gratings are important for thin-film solar cells since they exhibit high potential of light management. In order to simulate and optimize the periodic structures different approach is needed in comparison to widely used one-dimensional approaches. In this paper a two-dimensional finite element numerical approach is described to solve the Maxwell's differential equations of the propagating light in thin-film structures. The different boundary conditions that can be applied at the borders of the simulation domain are presented. The concepts how to efficiently solve the system and how to obtain the final optical parameters of the solar cells - reflectance, absorptances in the layers, charge carriers generation rates - are described. The model is verified with the measurements of the total reflectance at realistic diffractive gratings. Simulation results of a complete amorphous silicon solar cell with the diffractive gratings are presented.
Dvo-dimenzionalni optični model za simulacijo periodičnih optičnih struktur v tankoplastnih sončnih celicah
Kjučne besede: optično modeliranje, tankoplastne sončne celice, periodične strukture, vodenje svetlobe
Izvleček: Novi optični pristopi osnovani na periodičnih strukturah, kot so uklonske mrežice, so v tankoplastnih sončnih celicah pomembni, saj izkazujejo velik potencial pri vodenju svetlobe v strukturi. Za simulacijo periodičnih struktur ne zadostuje samo uporaba klasičnih eno-dimenzionalnih pristopov pač pa se moramo poslužiti dvo-dimenzionalnih simulacij. V prispevku je predstavljen dvo-dimenzionalen simulator, ki temelji na reševanju valovne Maxwellove enačbe širjenja svetlobe znotraj sončne celice na podlagi metode končnih elementov. Opisani so različni robni pogoji, ki jih uporabimo v simulacijah za zaključitev računskega območja. Podani so načini, kako določiti izhodne optične parametre simuliranih sončnih celic (odbojnost, absorpcija svetlobe v plasteh, profili generacij prostih nosilcev naboja). Rezultati simulacij so ovrednoteni s primerjavo izmerjenih odbojnosti izdelane uklonske mrežice. Na koncu so prikazani rezultati optične simulacije celotne tankoplastne amorfno-silicijev sončne celice.
1. Introduction
In thin-film (TF) solar cells light trapping is of great importance, since absorber layers are very thin (from one hundred nanometers to few micrometers). To capture light inside the absorber layer efficiently the scattering of light at internal interfaces is needed. In this way light paths throughout the layers are prolonged and the number of light passes increased due to enhanced back reflectances at internal interfaces. In current state-of-the-art silicon TF solar cells randomly textured substrates are used to introduce interface texture in the cell structures /1-3/. However, to further improve and use the potential of light scattering process, new approaches of light management are becoming important. One of them are periodic structures, such as diffractive gratings, which can be applied to the interfaces to efficiently scatter light only into specific, but very large angles /4/. Large angles can consequently lead to a significant increase of the light paths and to the total internal reflectance of the scattered light rays at the front interfaces of the solar cell structures. The integration of diffractive gratings in thin-film solar cell structures have already been investigated experimentally /5, 6/ and by numerical simulations /7-10/. The results of current structures of solar cells, including the gratings, indicate that the potential of the diffractive gratings has not been exploited fully. Thus,
further investigations of diffractive gratings and especially of their integration and optimization in the TF solar cell structures is of great importance. Numerical modelling and simulation present here an important tool.
To analyse optical situation at the grating, implemented in the structure of thin-film solar cells, at least a two dimensional simulator is needed. Some two- or three- dimensional (2-D, 3-D) simulators for the analysis of the diffractive gratings in the solar cells have already been reported e.g. /7-10/. However, the simulators have certain drawbacks, since they are mostly developed for general usage in the field of electromagnetic wave propagation and are not adapted and optimized to the specific structures, such as TF solar cells. One of the drawbacks is, for example, that is difficulty to simulate not-rectangular shapes of the grating features (triangles, sine shape etc.). The slopes of the different shapes are usually roughly approximated by rectangular steps, which can not give good representative results especially if only a few steps are used for approximation. Sharp edges of even small rectangles can affect the simulation results noticeably. Another drawback is that the output results are not the quantities that are important for TF solar cells (like optical reflectance, R, absorptance, A, charge carrier generation rate, Gl).
In this paper we present our 2-D optical simulator which has been developed for the analysis of TF solar cell structures including diffractive gratings and other periodic structures. The model is based on the finite element numerical method (FEM) /11/. One of the specialities of the model is that the simulation domain is represented by triangular elements rather than rectangular. These way different shapes can be more easily and effectively described. The details of the model and other advantages of the simulator, based on the model are presented. The boundary conditions that can be applied at the borders of the simulation domain are described. The determination of the electric and magnetic field distribution in the structure will be explained. The methods to calculate standard optical quantities used in the field of photovoltaics, such as R, A, Gl will be presented. Mathematical methods and techniques how to manipulate large data matrix and how to calculate out the desired results are briefly described. The verification results of the simulator and the simulations of the complete amorphous silicon solar cell are presented.
2. Optical model
2.1 Electromagnetic background
The optical model is based on Maxwell's wave equations in the frequency domain (Eqs. 1).
V X (—V X E) - CO ^eE =-yco j
Vx(-VxH)-o)Vh = Vx(-J) Eqs 1. e	e
where E and H are the complex vectors of the electric and magnetic field of propagating waves, ji = ßrßa = (ßRe -and £ = Er eo = (ers- jeim) £0 are the complex permeability and the complex permittivity of the material (jo = 4nx10"7 H/m, £o = 8.854x 10"12 F/m), (o is the angular frequency and J is the complex vector of the current density. In a 2-D space problem it is assumed that the field components (E and H) and media (jr, £r) do not change in the third dimension (axis z in our case). Thus, the Eqs. 1 can be written in the form of differential equations for two dimensions as given in Eqs. 2.
E =
dx ox dy Hr qy	.
= jkgZ^J^ - TE wave
dx dx dy Er qy
=	+	) - TM wave
dx E, dy E,
H =
Eqs 2
where x and y are the spatial directions (x - lateral, y -vertical), k0 is the wavenumber in vacuum and is defined
as ro-^eoM-o, ^0 is the impedance of vacuum ^ ^iijl^). The subscripts of the field and current components indi-
cate the directional component of the vectors. In Eqs. 2 it is considered that the electromagnetic wave is transversal (electric and magnetic field vectors perpendicular to each other) /11/. The wave is divided into the transverse Electric (TE) and transverse magnetic (TM) wave. In the case of TE wave, the wave is entirely represented by electric field component in z axis, Ez, whereas in the case of TM wave, the wave is entirely represented by the magnetic field in the z axis, Hz. Since the incident light consists both components of polarization (TE and TM), both equations included in Eq. 2 need to be considered.
The geometry of the periodic structure was described by choosing the triangular elements as basic building elements of mesh (Fig. 1a and 1b). Eqs. 2 have to be solved for each element in the mesh-grid, where we assume that the properties of material (j and £r) are constant within the element. The alignment of the triangular elements in the grid is shown in Fig. 1c. All the elements are of the same size.
a)
Fig. 1:
Basic elements used in our optical model: a) triangular element for linear approximation of unknown function (field), b) triangular elements for quadratic approximation of unknown function and c) discretization of simulation domain with triangular elements (with either linear or quadratic approximation).
For each triangular element the unknown function ^ (representing Ez or Hz) at the position of the element has to be determined along the three borders of the element. In our simulator we have implemented two options: a) linear approximation (Fig. 1a, Eq. 3) and b) quadratic approximation (Fig. 1b, Eq. 4) of the unknown function along the borders of the element.
if'{_x,y) = a' +b'x + c'y <lf'{x,y) = a' +b'x + c'y + d'x^ + e'xy + fy^
Eq. 3 Eq. 4
In Eq. 3 and 4, x and y are the 2-D spatial directions (see Fig. 1), where the superscript e indicates that the approximation is related to one-element level. The symbols a-f are the constants that have to be determined in the calculation process. In the case of linear approximation the function 0 is determined for the three nodes presenting the corners of the triangle (Fig. 1a). In the case of the quadratic approximation at each border line of the element an additional node is added, resulting in six calculation nodes
Fig. 2: Different boundary conditions for different configurations of simulation domain: a) periodic symmetric configuration, b) general periodic configuration and c) open region configuration - bounded system.
(Fig. 1b). The use of quadratic approximation has found to be useful when the abrupt changes in the field or spikes in the field are expected.
2.2 Boundary conditions
Determination of the size of the simulation domain and application of the mesh-grid to our periodic structure are crucial steps, which are both related to the boundary conditions applied to the borders of simulation domain.
Typically the structures are periodically repeated (infinitely) in lateral direction in our case. Two different types of boundary conditions applicable to the problem will be explained here. One type is related to the periodicity and lateral symmetry (left and right border according to examples in Fig. 2), whereas the other one is related to the incident (radiated) and outgoing field (top and bottom border of the simulation domain).
In the case of first type of boundary conditions (left, right border), application of three different boundary conditions is illustrated in Fig. 2a, 2b. In Fig. 2a the example of the application of the homogeneous Neumann condition /11/ is shown. In this case the first derivative of the field at the left and right border is zero (the divergence of field is zero). According to the figure this condition enables that only half of the period (P/2) is included in the simulation domain. Next, taking the whole period of the structure in the simulation domain (Fig. 2b) the field at the left and the right borders should be set to the equal values (virtually connected boundary system). The reason is the periodicity of the structure acquired in the simulation domain. Therefore, this condition is assigned to the periodic boundary condition /11/. The third condition which can be applied to the left and right border is absorbing boundary condition (ABC, Fig. 1c) /11/. This case is assigned to an open-region configuration in the lateral direction (structure with limited lateral dimensions, the whole structure is included in the simulation domain). An ideal boundary condition here would be zero reflectivity of waves at the borders. To approach this ideal case different orders of the ABC conditions can be used /11/.
For the second type of the boundary conditions (top, bottom border of the simulation domain) the above mentioned ABC condition is applied at the top and bottom border in all the cases.
In Fig. 2a and Fig. 2b for the illumination source the perpendicular plane wave is used, in Fig. 2c, a laterally limited illumination is applied (e.g. Gaussian beam). This enables us to simulate also the laterally limited non-periodic structures with limited area of light illumination.
In Figs. 2 the structures consists of only one layer (the diffractive grating below) and the incident medium above. However, in the model more layers (e.g. complete solar cell structures) can be applied with different type of the interface morphology. The vertical and lateral geometry are described with optical properties (e and ß) for each node in the mesh-grid.
2.3 Solving the system
To solve the differential equations (Eqs. 2) for each node of the element the Ritz method is applied /11/. Boundary conditions, linear or quadratic approximation of O, illumination source (plane wave or Gaussian beam) are considered. In the Ritz method the differential equations are transformed in the mathematical functional. Finding solutions for 0 at each node is based on minimising the functional / 11/. As a final result of the Ritz method the following matrix description is obtained (Eq. 5)
K0=b	Eq. 5
where K is the matrix of the system coefficients (the size of N x N; N - number of all nodes in the simulation domain), 0 is the vector (the size of N) of the system unknowns (E or H) and b is the vector (the size of N) of the light sources (different than zero only at the nodes at the top border, where the incident field is applied). The matrix K is sparse, symmetric, with complex numbers and it is not positive-definite. Its elements consist information of the material properties (e and ß) at the nodes.
It has to be noted that N in our case can even be more than a million, depending on the problem. Due to the large
size of the matrix K (Nx N) the memory consumption in the computer program can become a problem, especially if the direct methods for solving (e.g. Gauss elimination, LU decomposition /12/) are used, due to generation of new non-zero elements. To solve such a system, iterative methods are recommended /12, 13/. We found out that in our case the solution can be efficiently obtained by using non-stationary gradient iterative methods: the special form of Conjugate Gradient (CG) method /11/, Bi-Conjugate Gradient (Bi-CG) and Quasi-minimal residual (QMR) method / 12/. The CG method was found to be slow, but relatively stable, however the Bi-CG and QMR were found to be fast but less stable. In order to improve the condition of matrix and also to stabilize the method different preconditions were implemented into CG, QMR and Bi-CG algorithms. The easiest precondition used was the Jacobian precondition /13/, leading to a good convergence. However, in some special cases with the Bi-CG method we could still obtain divergence of the method, especially when the sharp metal structures were simulated, where the high spikes in the field might occur. In order to stabilize the method we have implemented the Symmetric Successive Over-Relaxation (SSOR) precondition /13/. The simulations with the SSOR precondition needed less iterations to obtain the result compared to Jacobian precondition.
The very large sparse system needed to be efficiently solved in the fastest time as possible and also with using very low memory consumption. Special attention has to be paid on the description of sparse matrixes in order to fasten the computation time when calculating the product of the sparse matrices with vectors in the calculation procedure, since these products have been found as the most time consuming. Another simplification considering the solving of system is to use mesh-grid consisted of the elements of the same size and of same orientation in the grid (pre-defined regular grid, Fig 1c). This way we do not need additional large matrix to describe the mesh.
2.4 Determination of the output parameters of the simulation
After the field Ez or Hz has been obtained at each node, the calculation of the final output parameters is performed. In our optical analysis, optoelectronic structures like thin-film solar cells are investigated. The following parameters are defined as the main output parameters: the total reflectance from the structure, Rtot, the absorptance inside individual layers of the structure, flayer, and the 2-D generation rate across the structure, Gl, of the photo generated electrons and holes in the active layers.
The Rtot of the structure is calculated at the top border, where also the incident wave is generated. The basis for the Rtot calculation presents the Poynting vector S = (E x H*)/2 ("*" presents the conjugated value). By considering only its normal direction (y axis) the Rtot can be determined as Rtot = Sy_refl / Sy_inc, where the subscripts "refl" and "inc" correspond to the reflected and the incident component
of the Poynting vector. These two components can be obtained from the calculated total and known incident field values at the top border (E refl = Etot - Einc , H refl = Htot -Hinc). This results in the final equation for Rtot (Eq. 6):



=
E dy

Eq. 6
Absorption inside the single element in the layer and inside entire layer A, can be calculated as (Eq. 7), which is derived from Poynting's theorem.
j^^^^HH'+j^EE'	Eq. 7
where ßi and e are imaginary parts of permeability and permittivity corresponding to the element. To obtain absorption of a layer one has to sum up all the absorptions of the elements which are composing specific layer.
Generation rate profile is calculated from the local absorption (absorption of the elements) and is given by Eq. 8 for one element e.

Eq. 8
where h is Planck constant, c is speed of light and /in illumination power density.
is
3. Verification results and simulations of the solar cells
Based on the presented optical model a computer simulator FEMOS-2D was developed. A user friendly interface enables simple simulation of the multilayer structures, including diffractive gratings, in the entire solar spectrum. The results of simulations were verified by comparing them to the measurements obtained on realistic samples. In Fig. 3 the structure of one of the samples as well as the Atomic Force Microscopy (AFM) scan of its surface is shown. On the polycarbonate substrate, which is typically used for CDs, DVDs or BDs, the periodic grating structure was embossed with the periodicity of P = 700 nm. On the top of the polycarbonate substrate 100 nm thick layer of an aluminium alloy was deposited. By means of AFM the shape of the grating surface was determined. In our case the sine shape was used in simulator to describe the grating shape with the height of h = 40 nm as determined from the AFM measurements. The sine shape agrees with the measured AFM profile very well, however, realistic shape of the profile can be imported in our simulator. The other input pa-
rameters for simulation were realistic wavelength dependent refractive indices of the layers and polarization of light. In simulations we used un-polarized light, 50 % of TE and 50 % of TM polarization, approaching to the realistic illumination in our measurement. In Fig. 4 the measurements and simulation of the Rtot of the grating as a function of light wavelength is shown. All measurements of Rtot were done with Lambda950 spectrophotometer from which unpolarized monochromatic light in the range from 400 to 1000 nm was obtained.
1.0
100 nm of Aluminium
Fig. 3: Analysed grating structure with thin film
aluminium layer (100 nm) on the top of the polycarbonate substrate. On the right the AFM measurement of the left structure is presented.
Good agreement is observed between the measured (dashed curve) and simulated (full curve) Rtot of the grating structure. For comparison also the simulation of the flat structure is shown (dash-dot curve). In this relatively simple structure the effect of the grating is related to the decreasing spike at the wavelength of 700 nm. This spike is due to additional absorption in aluminium alloy at the grating structure in the mentioned wavelength region. In this wavelength region the anti-reflective effect and light scattering in the first diffraction order /4/ occur. In the simulation and measurement (not shown in the figure) of the flat structure no spike is observed.
In the next step we simulated the whole thin-film amorphous silicon solar cell deposited on the 2-D grating (substrate configuration). The structure of the cell is as follows: Al/n-a-Si:H(20 nm) / i-a-Si:H(200 nm) / p-a-Si:H(10 nm) / ZnO:Al(500 nm) on the top (see insert in Fig. 5). At all interfaces the sine shape of gratings was used, with the height of 150 nm and with two different periods (300 and 400nm). The wavelength dependent refractive indices, N(k)
0.7
0.6
- Measurement - grating structure
---Simulation 2D - grating structure
- - Simulation 1D - flat structure
400 500 600 700 800 wavelength, X [nm]
900
1000
Fig. 4: Measured and simulated total reflectance of the aluminium grating structure with P=700 nm and h = 40 nm showing in Fig. 3.
of the realistic layers were used to determine the e(k) of the materials needed in the simulation.
In Fig. 5 the simulated absorptance in the i-a-Si:H layer of the solar cell is plotted for selected grating parameters. The grating should act as an efficient scatterer in the cell, leading to enhanced absorptance in the i-a-Si:H layer. Higher absorptance in the mentioned active layer leads to a higher short-circuit current and quantum efficiency of the solar cell. The 2-D simulations reveal the increase in absorptance is partly due to antireflective effect of the textured front interfaces and partly due to scattering effect of gratings at the interfaces. The simulations showed that by changing the period of the grating (in Fig. 5 shown for the case of P = 300 nm and 400 nm) in this case it is affecting the position of the interferences in the absorptance curve.
The presented simulator enables the direct study of the relation between the (periodic) surface morphology and light scattering. This is a very important point in the simulations of photovoltaic devices with textured interfaces (most of them), where the optical situation inside the structure cannot be measured, thus, simulations are needed to optimize the structures. This way the direction towards the optimal texturing can be indicated. Extending the simulator to the randomly textured interfaces it can be used to evaluate and select different TCO substrates from the light scattering point of view. However, analysis of the regularly textured interfaces (such as gratings) can already give useful information about improvements of randomly textured interfaces. Further on, special optical effects in the solar cell structure can be investigated with the simulator, such as plasmon absorption at textured metal back contacts. Minimizing the optical losses in metal, related to this absorption, by possibly optimizing the texture shapes is one of the important issues, other advantage of our simulator is also that it is highly customized for a solar cell application, but can be also used for other EM problem.
Fig. 5:
Simulated absorptance of the i-a-Si:H layer in the complete amorphous silicon solar cell structure with the sine grating applied to all interfaces.
4. Conclusions
Two-dimensional optical model for solving electromagnetic wave equations at periodical structures - diffractive gratings - was presented. The model is based on robust FEM method and solves differential equation for both (TE and TM) polarizations. Special attention was paid on boundary conditions and how to effectively solve the system discre-tized differential equations. This way we were able to obtain accurate results of simulation in the shortest time.
One of the main advantages of the optical model is simulation of arbitrary (periodic) interfaces shapes, where a good approximation of the interface texture can be achieved by using triangular elements instead of rectangular. The simulator based on the developed model is especially dedicated to simulation of optoelectronic structures such as thin-film solar cells. It automatically calculates the main optical output parameters from the field, such as absorption in individual layer of the structure, total reflectance and 2-D generation rate profile at each discrete element. The simulator is optimised for speed of the calculation and for low memory consumption to allow large number of discretization points.
Simulation results of Al based grating structure with the period of 700 nm is compared with the measured total reflectance of the sample. Good agreement is observed, indicating the validity of the simulations. The result of optical simulation of a complete amorphous silicon solar cell with different periods of the incorporated grating is presented. The developed optical simulator presents a powerful tool for further investigation of light management in
thin-film solar cell with the diffractive gratings and other textures.
5. Acknowledgement
The authors thank Miro Zeman from DIMES institute, Delft University of Technology, Netherlands for providing the grating sample and for useful discussions on optical modeling and OM&T B.V., Netherlands for providing polycarbonate substrates.
References:
/1/ Mueller et al., Solar Energy 77, 2004, p. 917
/2/ Fa' et al., Proceedings of EU PVSEC, Glasgow, 2000, p. 361
/3/ Kambe et al., Proceedings of WCPEC-3, Osaka, 2003, p. 1812
/4/ C. Heine, R. H. Morf, Applied Optics 34, 1995, p. 2476
/5/ S. H. Zaidi, J. M. Gee, D. S. Ruby, Diffraction grating structures
in solar cells, IEEE, 2000, pp. 395-398 /6/ N. Senoussaoui, M. Krause, J. Müller, E. Bunte, T. Brammer, H.
Stiebig, Thin Solid Film 451-452, 2004, p. 397 /7/ M. Niggemann, B. Bläsi, A. Gombert, A. Hinsch. H. Hoppe, P. Lalanne, D. Meisnner, V. Wittwer, 17th European Photovoltaic Solar Energy Conference Proceedings, 2002, p. 284 /8/ C. Haase, H. Stiebig, Prog. Photovolt: Res. Appl. 14, 2006, p. 629
/9/ R. H. Morf, J. Opt. Soc. Am. A 12, 1995, p. 1043
/10/ C. Haase, D. Knipp, H. Stiebig, Proceedings of the 22nd EU
PVSEC, Milano, 2007, p. 2186 /11/ J. Jin, The Finite Element Method in Electromagnetics, 2nd edition, John Wiley & Sons, New York, 2002 /12/ R. Barrett, et al., Templates for the Solution of Linear Systems,
2nd Edition, SIAM, Philadelphia, 1994 /13/ Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd Edition, SIAM, Philadelphia, 2003
Andrej Čampa, univ. dipl. inž. el. Asst. Prof.. Dr. Janez Krč, univ. dipl. inž. el. Prof. Dr. Marko Topič, univ. dipl. inž. el.
University of Ljubljana, Faculty of Electrical Engineering, Laboratory of Photovoltaics and Optoelectronics, Tržaška cesta 25, SI-1000 Ljubljana, Slovenia
E-mail: andrej.campa@fe.uni-lj.si