## Abstract

Finite element method is coupled with Huygens’ expansion to determine light intensity distribution of scattered light in solar cells and other optoelectronic devices. The rigorous foundation of the modelling enables calculation of the light intensity distribution at a chosen distance and surface of observation in chosen material in reflection or in transmission for given wavelength of the incident light. The calculation of scattering or anti-reflection properties is not limited to a single textured interface, but can be done above more complex structures with several scattering interfaces or even with particles involved. Both scattering at periodic and at random textures can be efficiently handled with the modelling approach. A procedure for minimisation of the effect of small-area sample, which is considered in the finite element method calculation, is proposed and implemented into the modelling. Angular distribution function, total transmission and total reflection of the scattering interface or structure can be determined using the model. Furthermore, a method for calculation of the haze parameter of reflected or transmitted light is proposed. The modelling approach is applied to periodic and random nano-textured samples for photovoltaic applications, showing good agreement with measured data.

© 2015 Optical Society of America

## 1. Introduction

Determination of the angular distribution function (*ADF*, a.k.a. angular intensity distribution – *AID*), the ratio of diffused over the total light (*Haze*) and total transmission and reflection of a single textured interface or complex optical structures is essential for indication about how such a structure will perform inside a solar cell or other optoelectronic device. While the mentioned parameters can be measured in air using angular resolved spectroscopy (ARS) systems [1–3], for analysis and optimisation of devices, determination of the scattering and antireflection parameters directly from the optical properties of the materials, and the geometry (morphology and position in the structure) of the scattering centres is of great interest. Using optical modelling light, scattering parameters, such as *ADF* and *Haze* inside or outside the device, in different materials and either in near-, transition- or far-field region can be determined. Although the approaches are applied to (thin-film) photovoltaic devices here, they are also applicable to other optoelectronic devices such as organic light emitting diodes, where the optical extraction efficiency and light intensity distribution outside the device is of great importance [4].

In the recent years, a few models for calculation of *ADF* (and haze parameter) of random nano-textured surfaces used in thin-film silicon solar cells have already been developed. The models have been based on different approximations of optical situation at a single nano-textured interface. Jäger and Zeman [5] presented a model for calculation of scattered light distribution at an interface based on the first order Born approximation. Dominé et al. [6] formulated a model based on the Rayleigh-Sommerfeld diffraction integral. Bittkau et al. [7] applied discrete Fourier transform (FT) to the scattering surface and then superimposed partial *ADF* functions of individual FT components. Finally, Jäger et al. [8] presented a further evaluated model based on the scalar scattering theory, achieving a good agreement with measurements for *ADF* of not only transmitted but also reflected light. While these models show relatively good agreement in angular intensity distribution for given randomly nano-textured interface, they can only handle a single (isolated) interface at a time. Furthermore, secondary (and subsequent) reflections from the same textured interface are not considered in these calculations. In addition, anti-reflection effect of nano-textured interfaces was not included in the presented models, thus total reflectance and transmittance of the interface cannot be determined. While, due to usage of Fourier transform (FT) for conversion of initial function to *ADF*, these methods are fast, FT also introduces a limited range of validity [9]. Since expansion using FT only considers first scattering (diffraction) order (for each sinusoidal component of the texture), higher deviations from real values are expected for random textures of higher root-mean-square (RMS) values [8,10] and even more so for periodic textures, which have become of interest in optoelectronic devices in the last years [11–13]. Furthermore, these models assume calculation of *ADF* at infinite distance from scattering structure (far-field), while in the case of thin-film and thinner wafer-based solar cells, distribution of the light close to the cell might be of more interest.

Besides approximative approaches, rigorous methods for determination of the scattered field have been proposed. They are based on combining rigorous calculation of Maxwell’s equations with either finite difference time domain (FDTD) [14], finite element method (FEM) [15] or others, with a method of near-field to far-field expansion, such as employing Huygens’ principle [16]. Observation at infinite far-field distance is often assumed in these cases while the effect of small simulation domain in FDTD or FEM simulations on angular distribution of light is not minimised (see section 3.3).

In this paper, a rigorous method was implemented coupling FEM with Huygens’ expansion in its base formulation. In contrast to typical implementations of such coupling, this form enables calculation of the *ADF* and other scattering parameters at a chosen distance in either near-, transition- or far-field region, making it much more applicable to the situation inside layers of the devices, with thicknesses starting from a few nm to several µm. We present an approach to minimise the effects of small calculation area limited by FEM, which appears to be the predominate artefact in the expanded *ADF*. Furthermore, the observation surface is not necessarily fixed to a hemisphere, but can be modified to best suit the analysed case (e.g. observation on a plane or at morphology of a neighbouring interface). In comparison to the approximate models, rigorous methods also enable calculation of the (anti-)reflective properties of the interfaces and do not limit the analysis to a single interface, but can handle entire optical structures (stack of layers, volume scatterers etc.) at once. In addition, the presented method for determination of the power of specular and diffuse light enables straightforward calculation of the haze parameter from results of any rigorous near-field simulation.

The modelling approach is applied to periodic nano-textures first as well as verified on selected samples with random textures comonly used in thin-film silicon solar cells.

## 2. Experimental

With the purpose of verification of the presented modelling approach, nano-textured samples were used in our analysis. In particular, two transparent conductive oxide (TCO) substrates used in thin-film silicon solar cell technology were analysed: (i) Low-pressure-chemical-vapour-deposited (LP-CVD) plasma-etched ZnO:B TCO of thickness around 2.3 µm, root-mean-square roughness *σ*_{RMS} ≈100 nm and autocorrelation length *ACL* ≈250 nm [17] and (ii) magnetron-sputtered and acid-etched (SP-E) ZnO:Al of *σ*_{RMS} ≈100 nm and *ACL* ≈500 nm [18]. Both TCO layers were deposited on a 1 mm thick glass substrate. Details of sample fabrications can be found in [16, 17]. The AFM images of the nano-textured TCO surfaces are presented in Fig. 1.

For optical characterisation of the samples, integrating sphere (IS) and angular resolved scattering (ARS) systems were used. For IS measurements, where total and diffused reflection (*R*_{TOT} and *R*_{DIFF}) and transmission (*T*_{TOT} and *T*_{DIFF}) were determined, the Perkin Elmer spectrometer Lambda 950 was utilised. For the ARS measurements, we used a goniometric system, where the *ADF* is measured in a horizontal plane around the sample. A 10 mW polarized He-Ne laser source with the wavelength of 633 nm was used for illumination. The distance from the sample to the Si detector rotating around the sample was 40 cm in all presented measurements. The samples were illuminated from the glass side for transmission and from the TCO side for reflection measurements.

## 3. Modelling

#### 3.1 The coupled rigorous modelling approach

In this approach rigorous FEM simulations [19] are employed to calculate the electric – *E*, and magnetic – *H*, field (hereinafter electromagnetic field) on a given surface – in our case plane *S* in the vicinity of the analysed structure. The plane *S* can be located in the half space of the transmitted light, as presented in Fig. 2, or of the reflected light, depending on calculation of scattered light in transmission or in reflection, respectively. The FEM part of the calculation is depicted in the left side of Fig. 2, whereas on the right side the Huygens’ expansion [20] is indicated.

Optical simulations of the electromagnetic field strength near the structure were done in COMSOL Multiphysics simulation software [15] in our case. Direct import of surface topology by atomic force microscopy (AFM) scans or otherwise generated surfaces (e.g. by means of non-conformal growth modelling [21]) is available in the program. In this work, full 3-D modelling with realistic optical properties (complex refractive indices) of the materials was used.

The structures in FEM simulations were enclosed on transmissive and reflective side using standard perfectly matched layer (PML) absorbing boundaries to model the required infinity of the incident and outgoing medium, assuring minimal back reflection into the simulation domain. On the remaining (lateral) sides of the simulation domain, symmetrical boundary conditions were used (perfect electric and perfect magnetic conductor). While for a given structure this selection can introduce a small error (by assuming a symmetric mirroring of the texture), it is still far better choice than, for example, applying the PML at the sides where a part of the light, which is propagating at high angles would be lost/absorbed resulting in erroneous *ADF*. Both transverse electric (TE) and transverse magnetic (TM) polarisation of the light were simulated for each case with equally weighted superposition of the electromagnetic field strength.

In the case of the calculation of the ADF in reflection, electromagnetic field strength has to be determined behind the excitation plane representing the illumination source. Since typical excitation used in COMSOL (surface current condition) generates waves in both, forward and backward direction, the electromagnetic field strength behind excitation plane would (besides reflected waves) also include the contribution of the backward plane wave. To avoid this, an additional excitation plane was added to the model in these cases, shifted by *λ*/4, both spatially and by phase (*λ* being the wavelength of the light in the material of the excitation). This results in a constructive interference of the plane waves in forward direction and in a destructive interference in the backward direction [20]. Thus, the plane wave appears only in the forward direction, while only reflected light is present behind the second excitation plane, as required.

In all presented cases, the interface or the structure lateral dimensions (simulation domain) and the surface *S* were both chosen to be of the fixed size for all cases, in particular 10 µm$\times $10 µm (with the only exception in section 3.2). While this size of the domain is still acceptable for simulation of a 3-D thin-film structure with a few interfaces, for more complex structures (larger simulation domains), compromises in the size become necessary. Average mesh density was about 5 elements per relative wavelength. All simulations were made considering light of a wavelength of 633 nm. Typical FEM simulation time was between 2 and 5 hours.

Huygens’ principle states, that the electromagnetic field outside a volume region *V* is completely determined by tangential fields specified over the surface *S* enclosing *V*, given that all contributing radiating sources are enclosed by *S* [16]. The electric and magnetic field strengths, represented by a complex vectors $\overline{E}$ and $\overline{H}$ outside the closed surface *S* are given as:

*S*is situated,

*µ*and

*ϵ*are permeability and permittivity of the said material, respectively, $\widehat{n}$ is normal vector, pointing out of

*V*. $\overline{r}$ and ${\overline{r}}^{\text{'}}$ denote the positions of the point of observation and of the contributing point on the

*S*, respectively. $\overline{\overline{G}}$ is dyadic Green’s function, given by:

*k*being wavenumber $(k=\omega \cdot \sqrt{\mu \u03f5})$ in outer media and

*g*being scalar Green’s function:

If electromagnetic field strength is known on a surface enclosing the analysed structure, resulting field at a chosen point outside *V* can be therefore calculated. However, one must consider that sources on (in general) enclosed surface *S* only contribute to the outward direction [22]. Typically, *S* is chosen in a manner, that the parts on *S* contributing to the field in a certain observation point are easily distinguished from the non-contributing ones. Nevertheless, solution for cases where this distinction is not trivial, was also proposed [22]. The surface *S* can be chosen as a sphere of infinite radius, with a surface point fixed in vicinity (> λ/2) of the structure – an infinite plane above the structure. In our case the infinite plane is furthermore reduced to a finite scan of electromagnetic field strength with the following argumentation: In general, infinite plane *S* could be approximated by numerous mirroring of the single scan of the electromagnetic field. However, such replication of the single scan does not only fall short of improving accuracy of the *ADF* calculation, but even leads to erroneous increase of the intensity in the angles linked to the period of the simulation domain. Therefore, only a single scan in the size of FEM simulation domain is used as the input for Huygens’ expansion – plane *S*.

From the FEM simulation, tangential components of *E* and *H* across the surface *S* are used as the input for Huygens’ expansion (see Eqs. (1)–(4)). For *S* situated in air and for chosen wavelength of light (633 nm), 100 nm sampling was found to be a good choice (no significant improvements in accuracy of ADF calculation were found for denser grid). However, the sampling should be adjusted accordingly to any changes of the relative wavelength of the light in the media where *S* is situated. Although, due to the symmetry boundary condition used, the distance between the structure and the surface *S* should not play a significant role, some differences were observed in our calculations for distances below *λ*/2, measured from the *S* to the closest feature point of the structure. Hence the distance was chosen to surpass this interval in our simulations (the average distance from the closest interface to the surface *S* was approx. 500 nm).

The calculation of Huygens’ expansion was implemented in MATLAB software [23]. To obtain the *ADF* in transmission or reflection, points of observation were chosen on a hemispherical grid for the analysed cases, to compare the results directly with the ARS measurements, which were done at a fixed distance around the sample. Once the field information on a grid of observation was obtained, the Poynting vector was used to calculate the power flow [16] – the intensity of light from electromagnetic field strength.

#### 3.2 Evaluation of the modelling approach

In order to evaluate the modelling approach, the *ADF* of light transmitted through a laterally limited 1-D grating with four slits (insert in Fig. 3) was calculated first. Ideal conductor is considered as the blocking material in this case. This example was chosen, since it can be compared to analytical solution for the finite-width slit grating derived from the Fraunhofer diffraction formula [20] for angular distribution function giving the intensity as:

*θ*is angle of the light propagation (0° being the normal direction to the slit-line),

*I*

_{0}is intensity of light at

*θ*= 0°,

*W*and

*P*are the width and the period of the slits, respectively.

*λ*is the wavelength of the light in exiting media, in this case air.

*N*is the number of the slits. This equation is, however only valid for the far-field region. As the grating is 1-D, the diffraction of the light is defined by a single angle in the space. The chosen four-slit grating was of

*W*= 500 nm,

*P*= 2 µm and the light of

*λ*

_{0}= 633 nm was applied in our calculation and simulation. In Fig. 3 the results of the far-field

*ADF*calculation are presented for the coupled modelling approach in comparison with the one calculated with the Fraunhofer diffraction formula. Distance of observation points from the sample was 1 m in our model, which is considered to be in the far-field region, as the calculated

*ADF*has no longer changed when increasing the distance. The size of

*S*for the 1-D grating was limited to 8 µm, presenting also the lateral size of the 2-D simulation domain. Since the structure is finite in this particular case, PML boundary conditions were also used at the sides. One can observe a good agreement between the calculations with the formula and simulations obtained with the coupled model. Some deviations are observed in the dips, which we mainly attribute to the finite discretisation in the FEM simulation and the Huygens’ expansion.

#### 3.3 Minimisation of the small-area effect

In the analysed case in Fig. 3 a small area of the structure was optically active (the region with the four slits) and thus it was entirely included in the simulation domain. In realistic cases, the structures and illumination spots are usually laterally much larger than the sizes we can include in simulation domains of rigorous FEM or other type (e.g. FDTD) of simulations. Limiting the size of illuminated area to a few µm or tens of µm in simulations induces an undesired effect. When light is travelling through an aperture it is subjected to diffraction [20]. The effect is only apparent for relatively small apertures (below a couple hundred wavelengths of light) and is negligible e.g. for the relatively large laser spot of the ARS measurements (the laser spot represents the aperture in this case as the sample itself is usually larger). In the case of Huygens’ expansion from a small-size surface *S* (limited by either the AFM scan size or FEM simulation domain), this effect appears to be significant and needs to be suppressed in order to obtain realistic *ADF*.

The diffraction effect due to the aperture and the scattering due to realistic features in the structure can be mathematically expressed as a convolution of the (sought) distribution function (*f*) with an *aperture function* (*g*) [20]:

*x*and

*y*denote the spatial coordinates. Function

*f*is related to the scattering structure (the sought-after

*ADF*), whereas the aperture function

*g*is related to the illuminated area. If the aperture function is known, however, deconvolution can be applied and the

*ADF*corresponding to the scattering features extracted. In the next paragraphs, a method for determination of the aperture function is described for the case of a normal incidence of plane wave.

If light is perpendicularly incident on a modelled structure, which does not include any scattering or refraction features (e.g. one or more flat interfaces), no component of electromagnetic field contributing to the diffracted light is present in the scan of the electromagnetic field on the *S.* If *S* is parallel to the interfaces, electromagnetic field strength is constant on the whole surface in this case. Therefore, spatially constant electromagnetic field is contributing solely to the specular component of light. This spatially constant electromagnetic field component, corresponding to the specular light, we call *static* component hereinafter. If the structure comprises scattering features, besides the static component, spatially changing – *dynamic* component is also present in the resulting electromagnetic field. Therefore, dynamic component contributes solely to the scattered light. Thus, the field on *S* can in general be written as:

If only static component of any electromagnetic field is considered and used in the calculation of the *ADF* with the coupled approach, the waves should propagate only in a single direction. Therefore, the sought function *f* is an impulse function (Dirac delta function) in this case. Since the impulse function is an identity element for convolution, the convoluted function *f* * *g* is equal to the aperture function (*g*) itself (see Eq. (6)). After *g* is determined, the deconvolution can be applied to the entire simulated *ADF* (calculated using the entire electromagnetic field), extracting the *ADF* component corresponding to the scattering features only and suppressing the artefact due to small aperture. Since the normal incidence of the light is prerequisite for Eqs. (7) and (8), for oblique incidence, analytical methods for the aperture function determination should be used [24].

To demonstrate the importance of the suppression of the small aperture effect on *ADF* we analyse an example of a 2-D sinusoidal texture in hexagonal lattice applied to a single interface (top view shown in the insert in Fig. 4). The interface geometry was generated as a superposition of three 1-D sine surfaces, with 120° angular shift between them. The interface was formed of a material with refractive index *n*_{1} = 2 (close to the refractive index of TCO’s) and *n*_{2} = 1 (air) in simulations. Light was incident as a plane wave perpendicularly to the surface from the *n*_{1} side. In Fig. 4 the results corresponding to the following cases are shown: (i) the *ADF* calculated from the entire field (including static and dynamic field), (ii) for the static component only (i.e. the 2-D aperture function) and (iii) for the deconvoluted *ADF* component, approaching the realistic case (of laterally large grating and large illumination spot). The deconvoluted curve corresponds to a line in 3-D *ADF* in Fig. 5(d) on the left-hand side. Blue dotted lines represent the analytically determined scattering orders (only angle, not the intensity). The dips in the *ADF* of the entire field (and the aperture function) are correlated to the size of the surface *S* (10 µm$\times $10 µm).

While in theory deconvolution leads to complete elimination of the unwanted aperture diffraction, this is not true for practical implementations (finite sampling and finite number of iterations). In cases shown here, the deconvolution was carried out in MATLAB, using the implemented Lucy-Richardson iterative method [24–26]. While the presented results of deconvoluted *ADF* are still not exact solutions, the effect of the aperture is significantly suppressed (intensity outside the diffraction modes is significantly lowered, note the logarithmic scale).

Although deconvolution can effectively minimise aperture diffraction, the studied section still has to be a good statistical representation of the texture. Recent publication showed that this condition is satisfied when using sample size of 5 times *ACL* or larger [27].

#### 3.4 Haze and reflective properties

Besides the angular distribution, anti-reflective properties of the textured interfaces and structures are also of interest. Since the calculation of an electromagnetic field above a structure is calculated by rigorous simulation, total transmittance (*T*_{TOT}) or total reflectance (*R*_{TOT}) can also be extracted from the relative power flow in given direction directly from the FEM simulations. Furthermore, besides *T*_{TOT} and *R*_{TOT} often the *Haze* parameter, giving the ratio between diffuse and total power flow in given direction (reflective/transmissive) is of interest for nano-textures. Since electromagnetic field strength contributing to the specular and to the diffuse part of the light, can be separated (see Eqs. (7) and (8)), associated power flow correlated to the specular $(\underset{\xaf}{P})$ and diffuse part $(\tilde{P})$ can be calculated for a given wavelength of light directly from the near-field:

*λ*as:

## 4. Application of the modelling approach

#### 4.1 Application to periodically textured interfaces

In the previous section, the coupled rigorous approach has already been applied to selected periodic textures for the purpose of its explanation. In next sections, it is furthermore applied to different textures and structures to demonstrate its usability and to validate it on selected samples (section 4.2). In this section, the hexagonal periodic texture is used again to simulate the 3-D *ADF*. Applicability of the modelling on the periodic textures is considered important since these types of textures show the potential to compete or even to outperform the improvements in cells brought by best random textures, if properly optimised [11]. As in section 3.3, the media on the side of the light incidence is chosen to have a refractive index of *n*_{1} = 2 in all cases. To demonstrate the usability of the model, different variations were done here. First, two different materials on the transmission side were selected: *n*_{2} = 1 (air – as before) and *n*_{2} = 3.5 (in the range of silicon layers). Next, the observation distance from the sample (radius of the hemisphere grid of the observation points) was varied: 10 µm ((A) – near-field region), 30 µm, 50 µm ((B) and (C) – transition region) and 1 m ((D) – far-field region). In our case we consider far-field region to be where the distribution pattern no longer changes significantly with increasing the distance from the sample (radius of the sphere) and near-field region where the pattern is strongly correlated to the structure morphology [28]. While here the materials are chosen not to absorb any light, complex refractive indices of absorptive materials can be used in general (see section 4.2). Simulation results of the *ADF* (only deconvoluted ones) in transmission in a 3-D half-space are presented in Fig. 5. *ϕ* is the *azimuthal angle* and *ϴ* is the angle of deviation from the normal direction – *deviation angle* (hemisphere grid is stretched to a disk, keeping *ϴ* spacing equidistant).

The results show that at 10 µm distance, the near-field intensity pattern is strongly correlated to the geometry of the interface nano-texture (see insert in Fig. 4), while evolving towards far-field pattern when increasing the distance towards 1 m. For distances above 1 m changes in *ADF* are no longer observed (not shown here). Considering these results, one should point out that for the case of thinner wafer-based and thin-film solar cells, near-field and transition region *ADF* might be more appropriate for evaluation of distribution of the light in cell than the far-field pattern.

Furthermore, the influence of the refractive index on the *ADF* is studied. This is significant, as the media inside a device is not air but other material (such as silicon). The results in Fig. 5 (right halves of (A)-(D)) show that the refractive index changes the near-, transition and far-field region *ADF* considerably. As seen for *n*_{2} = 3.5 the *ADF* pattern is narrowed towards smaller angles *ϴ* (towards the centre of the polar plot). Some intensities at higher angles are observed at distances 30 µm and 50 µm at right hand side, which are in fact numerical artefacts. With the presented approach however, an actual optical situation inside or outside the device can be determined; at the right distance, at a properly chosen surface and in actual material.

#### 4.2 Application to randomly textured interfaces

Randomly nano-textured substrates are commonly used in thin-film solar cells. To validate the model for random textures, the two TCO substrates (see section 2.) were used in calculation and results compared with measurements. Realistic complex refractive indices of the magnetron-sputtered ZnO:Al and the LP-CVD ZnO:B were used in the FEM simulation [29], while for Huygens’ expansion, absorption in transmission material was discarded, since it would only lead to proportional attenuation of the *ADF*. The corresponding refractive indices for the materials used are provided in the Table 1.

Besides the textured interface of air/ZnO (or p-*a*:SiC:H/ZnO), flat interfaces of ZnO/glass and glass/air were also taken into account in the modelling in this case – the complete structure of the actual measured sample. Furthermore, thick glass (1 mm) was considered in FEM simulations according to [30]. Light was incident as a plane wave perpendicularly to the sample from the specified side of the sample.

Simulated hemispherical *ADF*s for both TCO samples are shown in Fig. 6 (middle and the right hand side graphs). They were calculated for a distance of 1 m. Corresponding averaged cross-sections (azimuthal averaging around normal, for all angles *ϕ*) of the *ADF*s are shown on the left hand side and compared to the measurements of the samples (single horizontal scan). Roman numerals are used to pair the averaged cross-sectional to the corresponding hemispherical simulated *ADF*s. Since the same power of the incident plane wave was used in FEM simulations and the same observation distance was chosen for the *ADF* calculations, Figs. 6(a)–6(c) are directly comparable.

In Fig. 6(a) the *ADF* in transmission is shown for the LP-CVD ZnO:B sample (illumination applied from the glass side of the sample, transmission observed in air on the TCO side of the sample). In this case the *ADF* is shown after (graphs I) and also before the deconvolution (graphs II), to show the effect of the small-area sample correction (see section 3.3). Due to the additional diffraction caused by a relatively small-area surface *S,* overestimation of the intensity on the convoluted *ADF* in low and in high angles can be observed (graphs II). Since single direction propagation of each light wave is diffracted/spread across the neighbouring angles, the convolution is effectively expressed as a type of weighted angular (*ϴ*) averaging. This results in spreading of the specular component in low-angle area and prevents a fast drop of the intensity at high angles. Additionally it causes a smoother *ADF* in general. Although for larger part (except for the low and high angles), this mimics a better agreement to the measurement in this case, this is not the more accurate solution for the statistical properties given by the chosen size of the sample. Since the effective area of the calculation (*S*) is much smaller than the spot size of the laser, less dense frequency space (less components) of the sample is considered than in measurement. Therefore, the resulting calculated *ADF* exhibits some oscillation indicating the discrete modes in frequency spectrum (better statistical representation – larger sample – would smoothen these out). To reduce the oscillation of the deconvoluted *ADF*, one could enlarge the lateral dimensions of the simulation domain (and get a better statistical assessment of the sample). The second option is to analyse multiple scans of the texture in FEM separately and then average the obtained electromagnetic field scans.

In Fig. 6(b), simulated and measured *ADF*s are presented for LP-CVD ZnO:B texture in reflection. The sample was illuminated from the TCO side and reflected light was detected in air above the TCO. Additionally to this reflection (graphs III, comparable to the measurement) the *ADF* of reflected light is also shown for the material with *n*_{2} = 3.3 (IV), as it may be the case of p-*a*-SiC:H layer (without absorption) in thin-film silicon solar cell (TCO/p interface). Although one is typically interested in scattering parameters occurring in the absorbing material (e.g. i-*a*-Si:H), the main scattering is introduced at the interface between TCO and p layer, rather than p and i, since the optical properties of p and i layer are typically more evenly matched; hence p-*a*-SiC:H material was used. Good agreement of simulation to measurement is achieved also in this case. Furthermore, these results show similar trends that were observed by Jäger et. al [8], showing shift of intensity towards lower angles for the case where the exiting material is silicon.

In Fig. 6(c), *ADF*s in transmission are presented for the case of SP-E ZnO:Al. Illumination is applied from the glass side. *ADF*s are calculated in air (graphs V) and in p-*a*-SiC:H (without absorption) (graphs VI) above the TCO. Simulation and measurement are once again closely matched. In addition, for the case of the p-*a*-SiC:H for transmission, the intensity is pushed towards lower angles, which is in agreement with the measurement shown by Schulte et al. [31].

Comparison of the graphs I and V in Fig. 6 indicates that a large portion of the intensity shifts towards lower angles for SP-E ZnO:Al, which is consistent with its significantly larger *ACL* (approx. 500 nm) in comparison to the analysed LP-CVD ZnO:B (approx. 250 nm), since larger features tend to scatter light into lower angles. The higher propagation angle in general leads to higher generated photocurrent of solar cells due to the higher average optical path of the photons passing through the absorber layer. This, in agreement with previous studies [32,33], indicates that LP-CVD ZnO:B is typically a better choice of the two for state-of-the art thin-film silicon solar cells regarding the optical properties.

Finally, for all analysed random structures, total transmittances/reflectances and the haze parameter were also calculated and compared to the measurements where available. The results are shown in Table 2. Good agreement between simulated quantities and the measured-ones is observed. Comparison of the transmission of the LP-CVD ZnO:B and SP-E ZnO:Al sample in air (*n*_{2} = 1) shows only minor difference for reflective and scattering properties of both samples, therefore *ADF* carries important information regarding optical properties of given structures. Severe changes in cases where observation was changed from air to p-*a*-SiC:H indicate the importance of determination of the internal optical properties. This gives such modelling even greater significance.

## 5. Conclusions

A rigorous calculation of light scattering in optical structures was implemented; coupling FEM based optical simulations with the Huygens’ expansion. Such approach enables calculation of light intensity distribution in reflection and transmission at a chosen distance and chosen surface of observation in an arbitrary material. Total reflection/transmission of the structures can be determined as well. The optical analysis is possible above complex structures comprising many textured or flat interfaces or even volume scatterers. Furthermore, analysis of the textures of higher roughnesses and, besides random, also of periodic textures can be performed. A solution for minimisation of the effect of small-area sample is proposed. The coupled approach with small-sample effect minimisation presents an essential tool either for semi-coherent 1-D modelling approaches [34] as well as for novel coupled modelling approaches under development [35]. Good agreement is observed for the calculated and measured *ADF* and haze for all given cases, making such modelling a powerful, universal and versatile tool for optical analysis and optimisation of different optical devices.

## Acknowledgments

The authors thank EPFL (PV-LAB Neuchatel) and TU Delft (PVMD group) for providing TCO samples. M. Jošt is acknowledged for his help in performing *ADF* and haze measurements. The authors acknowledge the financial support from the Slovenian Research Agency (P2-0197 and J2-5466). M. Sever personally acknowledges the Slovenian Research Agency for providing PhD funding.

## References and links

**1. **K. Jäger, O. Isabella, L. Zhao, and M. Zeman, “Light scattering properties of surface-textured substrates,” Phys. Status Solidi, C Conf. Crit. Rev. **7**, 945–948 (2010).

**2. **S. Schröder, T. Herffurth, H. Blaschke, and A. Duparré, “Angle-resolved scattering: an effective method for characterizing thin-film coatings,” Appl. Opt. **50**(9), C164–C171 (2011). [CrossRef] [PubMed]

**3. **M. Jošt, J. Krč, and M. Topič, “Camera-based angular resolved spectroscopy system for spatial measurements of scattered light,” Appl. Opt. **53**(21), 4795–4803 (2014). [CrossRef] [PubMed]

**4. **W. Gaynor, S. Hofmann, M. G. Christoforo, C. Sachse, S. Mehra, A. Salleo, M. D. McGehee, M. C. Gather, B. Lüssem, L. Müller-Meskamp, P. Peumans, and K. Leo, “Color in the corners: ITO-free white OLEDs with angular color stability,” Adv. Mater. **25**(29), 4006–4013 (2013). [CrossRef] [PubMed]

**5. **K. Jäger and M. Zeman, “A scattering model for surface-textured thin films,” Appl. Phys. Lett. **95**(17), 171108 (2009). [CrossRef]

**6. **D. Dominé, F.-J. Haug, C. Battaglia, and C. Ballif, “Modeling of light scattering from micro- and nanotextured surfaces,” J. Appl. Phys. **107**(4), 044504 (2010). [CrossRef]

**7. **K. Bittkau, M. Schulte, M. Klein, T. Beckers, and R. Carius, “Modeling of light scattering properties from surface profile in thin-film solar cells by Fourier transform techniques,” Thin Solid Films **519**(19), 6538–6543 (2011). [CrossRef]

**8. **K. Jäger, M. Fischer, R. A. C. M. M. van Swaaij, and M. Zeman, “A scattering model for nano-textured interfaces and its application in opto-electrical simulations of thin-film silicon solar cells,” J. Appl. Phys. **111**(8), 083108 (2012). [CrossRef]

**9. **J. A. C. Veerman, J. J. Rusch, and H. P. Urbach, “Calculation of the Rayleigh-Sommerfeld diffraction integral by exact integration of the fast oscillating factor,” J. Opt. Soc. Am. A **22**(4), 636–646 (2005). [CrossRef] [PubMed]

**10. **D. C. O’Shea and S. of P. I. Engineers, *Diffractive Optics: Design, Fabrication, and Test* (SPIE, 2004).

**11. **C. Battaglia, C.-M. Hsu, K. Söderström, J. Escarré, F.-J. Haug, M. Charrière, M. Boccard, M. Despeisse, D. T. L. Alexander, M. Cantoni, Y. Cui, and C. Ballif, “Light trapping in solar cells: Can periodic beat random?” ACS Nano **6**(3), 2790–2797 (2012). [CrossRef] [PubMed]

**12. **L. Müller-Meskamp, Y. H. Kim, T. Roch, S. Hofmann, R. Scholz, S. Eckardt, K. Leo, and A. F. Lasagni, “Efficiency enhancement of organic solar cells by fabricating periodic surface textures using direct laser interference patterning,” Adv. Mater. **24**(7), 906–910 (2012). [CrossRef] [PubMed]

**13. **A. Abass, C. Trompoukis, S. Leyre, M. Burgelman, and B. Maes, “Modeling combined coherent and incoherent scattering structures for light trapping in solar cells,” J. Appl. Phys. **114**(3), 033101 (2013). [CrossRef]

**14. **A. Taflove, A. Oskooi, and S. G. Johnson, *Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology* (Artech House, 2013).

**15. **“COMSOL Multiphysics Modeling and Simulation Software,” http://www.comsol.com/.

**16. **J. A. Kong, *Electromagnetic Wave Theory*, 2 Sub (Wiley-Interscience, 1990).

**17. **O. Kluth, B. Rech, L. Houben, S. Wieder, G. Schope, C. Beneking, H. Wagner, A. Loffl, and H. W. Schock, “Texture etched ZnO : Al coated glass substrates for silicon based thin film solar cells,” Thin Solid Films **351**(1-2), 247–253 (1999). [CrossRef]

**18. **S. Faÿ, J. Steinhauser, N. Oliveira, E. Vallat-Sauvain, and C. Ballif, “Opto-electronic properties of rough LP-CVD ZnO:B for use as TCO in thin-film silicon solar cells,” Thin Solid Films **515**(24), 8558–8561 (2007). [CrossRef]

**19. **J.-M. Jin, *The Finite Element Method in Electromagnetics*, 2nd ed. (Wiley-IEEE, 2002).

**20. **M. Born, E. Wolf, A. B. Bhatia, P. C. Clemmow, D. Gabor, A. R. Stokes, A. M. Taylor, P. A. Wayman, and W. L. Wilcock, *Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light*, 7th edition (Cambridge University, 1999).

**21. **M. Sever, B. Lipovšek, J. Krč, A. Čampa, G. Sánchez Plaza, F.-J. Haug, M. Duchamp, W. Soppe, and M. Topič, “Combined model of non-conformal layer growth for accurate optical simulation of thin-film silicon solar cells,” Sol. Energy Mater. Sol. Cells **119**, 59–66 (2013). [CrossRef]

**22. **D. A. B. Miller, “Huygens’s wave propagation principle corrected,” Opt. Lett. **16**(18), 1370–1372 (1991). [CrossRef] [PubMed]

**23. **“MATLAB - The Language of Technical Computing,” http://www.mathworks.com/products/matlab/.

**24. **P. A. Jansson, *Deconvolution of Images and Spectra: Second Edition* (Courier Corporation, 2014).

**25. **W. H. Richardson, “Bayesian-based iterative method of image restoration,” J. Opt. Soc. Am. **62**(1), 55–59 (1972). [CrossRef]

**26. **L. B. Lucy, “An iterative technique for the rectification of observed distributions,” Astron. J. **79**, 745 (1974). [CrossRef]

**27. **K. Jäger, D. N. P. Linssen, O. Isabella, and M. Zeman, “Ambiguities in optical simulations of nanotextured thin-film solar cells using the finite-element method,” Opt. Express **23**(19), A1060–A1071 (2015). [CrossRef] [PubMed]

**28. **R. C. Johnson, H. A. Ecker, and J. S. Hollis, “Determination of Far-Field Antenna Patterns from Near-Field Measurements,” Proc. IEEE **61**(12), 1668–1694 (1973). [CrossRef]

**29. **J. Müller, B. Rech, J. Springer, and M. Vanecek, “TCO and light trapping in silicon thin film solar cells,” Sol. Energy **77**(6), 917–930 (2004). [CrossRef]

**30. **A. Čampa, J. Krč, and M. Topič, “Two approaches for incoherent propagation of light in rigorous numerical simulations,” Prog. Electromagnetics Res. **137**, 187–202 (2013). [CrossRef]

**31. **M. Schulte, K. Bittkau, K. Jäger, M. Ermes, M. Zeman, and B. E. Pieters, “Angular resolved scattering by a nano-textured ZnO/silicon interface,” Appl. Phys. Lett. **99**(11), 111107 (2011). [CrossRef]

**32. **A. Čampa, M. Meier, M. Boccard, L. V. Mercaldo, M. Ghosh, C. Zhang, T. Merdzhanova, J. Krč, F.-J. Haug, and M. Topič, “Micromorph silicon solar cell optical performance: Influence of intermediate reflector and front electrode surface texture,” Sol. Energy Mater. Sol. Cells **130**, 401–409 (2014). [CrossRef]

**33. **Y. Kim, R. Santbergen, K. Jäger, M. Sever, J. Krč, M. Topič, S. Hänni, C. Zhang, A. Heidt, M. Meier, R. A. C. M. M. van Swaaij, and M. Zeman, “Effect of substrate morphology slope distributions on light scattering, nc-Si:H film growth, and solar cell performance,” ACS Appl. Mater. Interfaces **6**(24), 22061–22068 (2014). [CrossRef] [PubMed]

**34. **J. Krč, F. Smole, and M. Topič, “Analysis of light scattering in amorphous Si:H solar cells by a one‐dimensional semi‐coherent optical model,” Prog. Photovolt. Res. Appl. **11**, 15–26 (2003). [CrossRef]

**35. **M. Topič, M. Sever, B. Lipovšek, A. Čampa, and J. Krč, “Approaches and challenges in optical modelling and simulation of thin-film solar cells,” Sol. Energy Mater. Sol. Cells **135**, 57–66 (2015). [CrossRef]