Wave optical simulation of the light trapping properties of black silicon surface textures

: Due to their low reflectivity and effective light trapping properties black silicon nanostructured surfaces are promising front side structures for thin crystalline silicon solar cells. For further optimization of the light trapping effect, particularly in combination with rear side structures, it is necessary to simulate the optical properties of black silicon. Especially, the angular distribution of light in the silicon bulk after passage through the front side structure is relevant. In this paper, a rigorous coupled wave analysis of black silicon is presented, where the black silicon needle shaped structure is approximated by a randomized cone structure. The simulated absorptance agrees well with measurement data. Furthermore, the simulated angular light distribution within the silicon bulk shows that about 70% of the light can be subjected to internal reflection, highlighting the good light trapping properties.


Introduction
Pyramidal textures are the industrial standard front side structure for monocrystalline silicon solar cells.In addition to antireflection properties, the pyramids redirect incoming light to different angles, thereby achieving a path length enhancement in the silicon and consequently increasing the absorption in the infrared part of the spectrum in comparison to a planar front side.For very thin solar cells light trapping is essential, whereas the application of a pyramidal front texture gets more and more difficult due to its size.Hence, other concepts are needed.
One promising approach is the use of so-called "black silicon".In general, this term describes a nanostructured, stochastic surface with very low reflectance over a wide range of the solar spectrum [1][2][3][4][5][6], thus making it look black.The height of the structure is typically in the range of only a few hundred nanometers to a few micrometers.There are different manufacturing methods available for black silicon leading to various structure geometries and optical properties [1][2][3][4][7][8][9][10]. In this paper, we investigate black silicon fabricated by inductively coupled plasma reactive ion etching (ICP-RIE) [4] as an example.A scanning electron microscope (SEM) picture of the structure is shown in Fig. 1(a).Figure 1(b) shows the measured reflectance.Even without an anti-reflection coating, the reflectance in the range of 280 nm to 1000 nm weighted by the photon flux given by the Air Mass 1.5 global spectrum (AM1.5g)[11] is only 1.2% compared to the 1.5% to 2% that can be reached with a combination of pyramidal front side texture and an antireflection coating [12].This highlights the excellent antireflection effect of black silicon.For the application of black silicon in solar cells, however, besides its optical properties, it is also important that the large surface can be electrically passivated.This has been achieved using aluminum oxide deposited by atomic layer deposition [2,4,7,13].The reflectance in the range of 280 nm to 1000 nm weighted by the Air Mass 1.5 global spectrum (AM1.5g)[11] is only 1.2%.
On top of the excellent anti-reflection properties, black silicon is also known to redirect light during its passage through the structured surface [2,3,7,14,15].Such redirection leads to path length enhancement and thus increased absorption, which is especially important in thin silicon solar cells.To further increase absorption, front-side textures and rear-side structures can be combined.Possible rear side structures could be hexagonal sphere gratings [16] nanoimprinted diffraction gratings [17] or random pyramids combined with a distributed Bragg reflector [15].
Clearly, the light redistribution induced by the black silicon affects functionality and optimal design of the rear-side structures.Therefore, to understand the optical behavior of black silicon, to calculate its impact on solar cell performance and for optimization purposes, especially in conjunction with a rear-side structure, a theoretical description of the structure is required.In particular, the angular distribution of light in the silicon bulk after passage through the surface texture is of interest [18,19].
For these purposes, Kroll et al. presented a formalism combining the Fourier Modal Method (FMM) and an incoherent propagation method to calculate the absorptance of a wafer with a black silicon front side texture [2].They also determined the angular distribution within the silicon bulk.To build the model, they measured the surface topography by taking cross section SEM micrographs between a series of milling steps using a focused ion beam.The angular distribution was then calculated using the finite differences time domain method (FDTD) [14].
In this paper, we introduce an alternative method for simulating the reflection, transmission and absorption of a wafer with black silicon as the front side texture and the angular distribution of light in the silicon bulk after passing through a black silicon surface.For this we approximate the black silicon with simple random cone structures and use the rigorous coupled wave analysis (RCWA) for computation.The simplified topography in our method allows for a more direct understanding of what creates the functionality of the black silicon and would also allow for a more targeted approach in future structure optimizations, for example by directly investigating the influence of feature height, density of "needles" (cones) etc.
In Section 2, the used methods are presented.The RCWA parameters, the implementation of the surface structure in the RCWA code and the characterization method are introduced.In Section 3 simulation results are discussed.First, reflectance, transmittance and absorptance are compared to measurement data and afterwards the angular distribution of light in the silicon bulk is presented.All results are compared to a Lambertian scatterer (diffusely reflecting surface with intensity of scattered light proportional to the cosine of the angle between the direction of the scattered light and the surface normal).Section 4 summarizes the main results.In the appendix, unit cell effects in the simulation results due to the periodic boundaries required for RCWA are discussed.

Rigorous coupled wave analysis (RCWA)
The simulation of light propagation in black silicon was conducted using RCWA [20][21][22][23][24], since it is a very efficient formalism for systems including thick homogeneous layers such as a silicon wafer.This wave optical simulation method provides theoretically an exact solution of Maxwell's equations.As it is an infinite summation series one has to truncate it at a certain point determined by convergence analysis.For this work the RCWA implementation RETICOLO CODE 2D written by J.P. Hugonin and P. Lalanne was used [25].The RCWA only works for periodic structures.A unit cell, which is continued periodically, has to be defined.The implementation of the stochastic black silicon structure in a periodically continued unit cell is described in section 2.2.
For the simulation, three regions have to be defined as sketched in Fig. 2(a).The first region is a homogenous incident medium with a constant refractive index.In this work we consider air as the incident medium with a refractive index of 1.The second region is the periodic structure.To consider the dependence of the refractive index in the z direction, several layers can be defined in the second region.Within each layer a periodic function of the complex refractive index in the x and y direction, respectively, is required.For the black silicon simulation we use 160 layers corresponding to a layer thickness of 6.25 nm at a total height of 1 µm.This value was determined by convergence analysis.The third region is the substrate, which has to be a homogenous medium again.
For the RCWA simulation of reflection and transmission, we use air as a substrate medium.The silicon bulk is considered as an additional 250 µm thick layer in the second region as shown in Fig. 2(b).Thus, the result of the simulation can be compared directly to the measurement data.As the RCWA simulation proceeds coherently, there are strong oscillations in reflectance and transmittance due to interference effects in the silicon bulk layer.Therefore, reflectance and transmittance are calculated with an increment in wavelength of 1 nm and then averaged over the 20 neighbouring points at each wavelength.
For the determination of the angular distribution of light in the silicon bulk, the substrate is simulated as silicon as shown in Fig. 2(c) without absorption (the extinction coefficient is 0).In principle, the difference in the extinction coefficient between the silicon in region II and III leads to additional reflection.However, the reflectance at a boundary between silicon without extinction coefficient and absorbing silicon is less than 10 7 in the wavelength range of 0.9 µm to 1.2 µm which is of interest for this work.
All values for the refractive index and extinction coefficient of silicon are taken from [26].In addition to the number of layers, the number of vertices with which circular elements within a layer are approximated has to be defined.After convergence analysis we fixed this value to 24.
Due to the periodicity, only a certain number of discrete diffraction modes can propagate in the considered medium, based on the grating equation for normally incident light Here, d is the periodicity, θ the polar angle of the diffracted wave, λ the wavelength and n the refractive index.p and q are the mode numbers in the x and y directions, respectively.As the sine term must not be greater than 1 the number of pairs (p,q) and thus the number of modes that can propagate is limited.The number of modes included in the computation was determined by convergence analysis for each unit cell size considered.For the simulation of the total reflection and transmission the sum of all reflected or transmitted modes is important.Based on the convergence analysis at least two evanescent modes are required.For the simulation of the angular distribution the diffraction efficiencies in each individual mode have to converge.According to the convergence analysis, in this case, at least three evanescent modes are required.This difference can be explained by the fact that a higher number of modes leads to a higher resolution of the structure representation in the model.These structure details are more important for the angular distribution.

Structure implementation
To implement black silicon in RCWA, the stochastic texture must be converted into a periodic structure.Therefore, a unit cell must be defined which can be repeated both in the x and y directions, respectively.Within this unit cell a random structure is generated.Quadratic unit cells are used.The black silicon needles are approximated with cones.The height of the cone and the diameter at the base are assigned randomly, such that all heights between 0.6 µm and 1 µm and base diameters between 0.2 µm and 0.4 µm, respectively, occur with the same probability.These values are extracted from scanning electron microscope (SEM) pictures (see Fig. 3).In addition, each cone is allocated to a random position (x and y) in the unit cell.
The SEM pictures also show a density of needles slightly above 20 needles in an area of 1 µm 2 .No significant planar areas can be seen in the SEM pictures.In order to match the actual needle density while avoiding planar areas, the structure generation starts with a cone density between 19 and 20 cones per µm 2 .For all planar areas with a diameter greater than 0.2 µm (resulting from the first randomized positioning of cones), an additional cone is added.Thus, in the end, the structure has between 20 and 23 cones per µm 2 , which corresponds to the value determined by the SEM pictures.If a cone overlaps the edge of the unit cell, the overlapping part is repeated on the opposite side, so that there are no discontinuities in the structure when repeating the unit cell.Figure 3 shows the SEM pictures, which the cone parameters are based on, and an example of a generated randomized cone structure.

Graded index layer
In order to assess the light trapping properties of black silicon, the absorption of a silicon wafer with a graded index layer at the front side has been simulated.A sufficiently thick graded index layer works as a perfect antireflection coating without any light path length enhancement.For the simulation the transfer matrix formalism [27] was used.1000 layers were considered, each 1 nm thick with linearly increasing refractive indices and extinction coefficients from the values of air to the values of silicon.A 250 µm thick silicon bulk layer was also incorporated.Accounting for the oscillations appearing in this thick silicon layer due to the coherent calculation, the curve was smoothed by averaging over an interval of 20 nm.

#251492
The absorptance of a graded index layer is similar to the absorptance of a regular cone structure with cone dimensions in the range of the black silicon needle sizes (diameter 0.2 µm to 0.4 µm).For a regular grating significant light trapping only occurs for periods larger than 0.5 µm [28].

Characterization
In order to compare the simulation result with measurement data, 250 µm thick wafers with black silicon on the front side and a planar rear side were characterized.The black silicon was fabricated by inductively coupled plasma reactive ion etching [4].Reflectance R and transmittance T were measured by photo spectrometry using a Varian Cary5000.Absorptance A is then determined from the measured values according to All experimentally determined spectra presented in this work were measured using a spectral resolution of 1 nm.

Reflectance, transmittance and absorptance
Reflectance, transmittance and absorptance were determined for unit cell sizes between 1.35 µm and 2.01 µm corresponding to between 69 and 161 propagating modes, e.g.pairs of mode numbers (p,q) at a wavelength of 1.0 µm.Since, in reality, black silicon is not a periodic structure, the unit cell size used in the simulation should be as large as possible.However, the size of the unit cell is limited due to high computing time.The larger the unit cell the more modes have to be considered in the simulation.This leads to a significant increase in the computation time and also in the size of temporary files generated during the calculation.For these reasons the simulations are limited to the spectral range between 1.0 µm and 1.1 µm.Within this range, light trapping effects are clearly visible but parasitic absorption is not yet as important as for longer wavelengths.Figure 4 shows reflectance, transmittance and absorptance of simulated structures for two different unit cell sizes.For the smaller cell with an edge length of 1.40 µm 15 random structures were generated as described in Section 2.2.For each structure, reflectance, transmittance and absorptance were calculated.The graph shows the mean value of these 15 structures.The standard deviation is represented by error bars.For the larger unit cell with an edge length of 2.01 µm, 10 random structures were considered.For comparison, the Yablonovitch limit [29], i.e., the absorptance of a silicon wafer with a perfect antireflection coating at the front side and a Lambertian scatterer at its rear, was analytically calculated.In addition, the absorptance of a silicon wafer with a graded index layer at its front is shown.
One can see that the absorptance is significantly higher than in the case of a graded index layer.Therefore, one can conclude that some light path enhancement in the silicon is achieved due to the black silicon front side.However, the absorptance is still lower than for a Lambertian scatterer, indicating that further absorption enhancement is possible by improved scattering properties of the black silicon or the application of light trapping structures to the rear side.For both unit cell sizes, the simulated spectra shown in Fig. 4 match the measured curves.This is, however, not necessarily the case for all unit cell sizes.The discrete angular distribution within the silicon bulk, composed of propagating modes that are induced by the periodicity, can lead to significant unit cell effects in the transmittance spectra.For a given unit cell size, depending on the wavelength, more or less modes can propagate.Some of the modes correspond to polar angels below the critical angle such that the light can leave the silicon at the rear side, while other modes are trapped by total internal reflection.Significant steps can be seen in the transmittance spectra at the wavelength where a mode overcomes the limit from transmission to total internal reflection.With increasing unit cell size these effects decrease because the number of propagating modes increases and thus a single mode is less important.A detailed analysis of these unit cell effects can be found in the Appendix.This analysis indicates that for periods larger than 2 µm the unit cell effects are negligible.

Angular distribution
Figure 5 shows the simulated angular distribution within the silicon bulk after passage through the modeled front side structure at a wavelength of 1.04 µm for the same unit cell sizes presented before.The graph shows the percentage of the total transmittance through the cone structure in the silicon bulk in dependence of the polar angle.The same structures as presented in Fig. 4 were used for computation and the results were averaged again.Due to the periodic structure only discrete modes can propagate.The corresponding polar angles are given by the grating equation.As we are interested in the deviation of normally incident light from the vertical direction, the irradiance is summed over all azimuth angles for every polar angle.The larger the unit cell the more modes can propagate and thus the more polar angles are possible.The two angular distributions presented in Fig. 5 show that different angular distributions can lead to similar reflectance, transmittance and absorptance curves (see Fig. 4).To determine the angular distribution in the silicon bulk after passage through a black silicon front side with higher accuracy, three different unit cell sizes with edge lengths larger than 2 µm were considered (2.01 µm, 2.20 µm and 2.40 µm, corresponding to 161, 193 and 233 propagating modes at a wavelength of 1.0 µm, respectively).According to the unit cell analysis presented in the appendix, undesirable effects should be negligible for these periods.At every size 20 random structures were created as described in Section 2.2.For each structure the angular distribution in the silicon bulk was first computed at wavelengths of 0.9 µm, 1.0 µm, 1.1 µm and 1.2 µm and the results then averaged for each unit cell size separately.Since the discrete polar angles are at different positions for the different unit cell sizes, a binning is necessary when averaging over the three angular distributions.Angular intervals of 5° each were chosen.The fraction of transmittance values was summed up in every interval for each unit cell size.Afterwards the obtained binned distributions of the three unit cell sizes were averaged.Figure 6 presents the results in comparison to the angular distribution of a Lambertian scatterer.To get comparable values, the cosine distributed intensity of the Lambertian scatterer has to be integrated over the solid angle associated with the appropriate polar angle interval [30].If θ is the polar angle and φ the azimuth angle, the proportion of the total intensity I in a polar angle interval [ 12 , ] is given according to The angular distributions qualitatively confirm the results of the simulated absorptance shown in Fig. 4.During transmission through the black silicon, light is redistributed to different polar angles.However, for high angles the fraction of transmitted light is lower than in the case of a Lambertian scatterer.Thus, the scattering effects due to black silicon are not as distinct as they would be with a Lambertian scatterer.
With increasing wavelength the proportion of light that is directly transmitted (polar angle 0 °) increases.The explanation is that for larger wavelengths the ratio between wavelength and structure size increases and the structure behaves more and more like a graded index layer without scattering effects.
The critical polar angle separating the escape cone of total internal reflection from the solid angle where total internal reflection occurs at the rear side is between 16.2° and 16.4° for the wavelength range considered.Depending on the wavelength, between 25% and 35% of the light is directed into the escape cone.These values are comparable to the results of Kroll et al., who investigated two different black silicon structures with between 20% and 40% of the light in the loss cone [14].

Summary
In this paper, we presented a way of implementing black silicon in rigorous coupled wave analysis (RCWA) code by describing the acicular structure with a rather simple structure of randomly sized cones.With this wave optical method, not only reflectance, transmittance and absorptance of a silicon wafer with a black silicon front side can be determined, but the angular distribution of light in the silicon bulk as well, which is crucial for further optical optimization of silicon solar cells with nanostructured surfaces.In particular, the combination of black silicon with diffractive rear side structures could lead to improved light path length enhancements.
Concerning reflectance, transmittance and absorptance we achieved very good agreement between simulation and measurement data.For the angular distribution of light within the silicon we found clear indication of scattering effects, however, less pronounced than for a Lambertian scatterer.The simulated angular distributions are also in good agreement with the data calculated by Kroll et al., who used the finite differences time domain method (FDTD) and a more laborious way to generate the modelled black silicon structure.

Appendix
Depending on the period, unit cell effects can be seen in the simulated spectra.Figure 7 shows differences between the simulated and measured transmittance for different unit cell sizes.15 random structures were considered for every period with the exception of the largest unit cell, where only 10 structures were taken into account due to the high computational effort.For periods of 1.35 µm and 1.40 µm the simulated spectra match the measured data within the statistical error.For the three subsequent unit cell sizes with periods of 1.46 µm, 1.49 µm and 1.54 µm, the simulated transmittance is significantly higher than the measured values at the beginning of the spectral range considered.With increasing wavelength, the differences vanish one by one.For the periods of 1.59 µm and 1.70 µm the simulated transmittance is too high compared to the measurement data, but with increasing unit cell size the difference decreases and for a period of 2.01 µm the simulated spectrum matches the measured curve well again.To explain the curve progressions in Fig. 7, one has to consider the polar angles of the propagating modes.Each propagating mode has a certain azimuth and polar angle.By integration over the azimuth angle, the intensity occurring under a certain polar angle can be calculated.The discrete set of modes results in a discrete set of polar angles.
Table 1 contains the number of different polar angles for which propagating modes exist for all periods considered in Fig. 7 and for three different wavelengths.The critical polar angle separating the loss cone from the solid angle of total internal reflection at the rear side is 16.27° at a wavelength of 1.00 µm, 16.33° at 1.04 µm and 16.38° at 1.08 µm.Within the loss cone (for polar angles below the critical angle), light can leave the silicon bulk at the rear side; for higher angles light is totally internally reflected.In Table 1, also the number of polar angles for which propagating modes exists that are below the critical angle is listed.
For the two smallest unit cells of periods 1.35 µm and 1.40 µm for two polar angles within the loss cone propagating modes exists, and modes associated with the third polar angle are already subjected to total internal reflection.For a period of 1.46 µm and a wavelength of 1.00 µm, three polar angles occur that are smaller than the critical angle, thus leading to outcoupling of the associated propagating modes.Hence, the transmittance is higher than for the smaller periods.At a wavelength of 1.04 µm only two polar angles occur that are smaller than the critical angle, thus the transmittance decreases.For a larger period of 1.49 µm, the transition between 3 and 2 occurring polar angles below the critical angles, takes place at a higher wavelength.Therefore, in Fig. 7 one can see a decrease in transmittance between 1.04 µm and 1.08 µm.The angular distributions for a period of 1.49 µm are shown in Fig. 8.With further increasing the period, the difference between the simulated and measured transmittance decreases, because the larger the unit cell, the more propagating modes exists altogether and the smaller the influence of one single mode.For the largest unit cell considered with a period of 2.01 µm the measurement data are matched well by the simulation.One can conclude that for simulation, unit cells of a period of at least 2 µm have to be used to get rid of unit cell effects.

Fig. 1 .
Fig. 1.(a) Scanning electron microscope (SEM) picture of black silicon fabricated by inductively coupled plasma reactive ion etching as considered in this work.(b) Measured reflectance of a 250 µm thick wafer with black silicon on the front side and a planar rear side.The reflectance in the range of 280 nm to 1000 nm weighted by the Air Mass 1.5 global spectrum (AM1.5g)[11] is only 1.2%.

Fig. 2 .
Fig. 2. (a) Sketch of the three regions required for the RCWA simulation.Incidence medium (I) and substrate (III) are homogenous.In between is a periodic structure (II).The z direction is discretized in a certain number of layers which are periodic in the x and y directions, respectively.(b) For the simulation of reflectance R, transmittance T and absorptance A, silicon is considered as a thick layer in region II.(c) For the simulation of the angular distribution in the silicon bulk, also the substrate is simulated as silicon.

Fig. 3 .
Fig. 3. SEM pictures of black silicon (left) compared to the implemented randomized cone structure (right).The cone parameters (height, diameter) were extracted from the cross section (a).The number of cones per area can be determined from the SEM picture in (b).In (c) the area covered by the cones is shown for four different sectional planes.z = 0 µm is the level of the silicon bulk.Blue represents silicon, white represents air.In (d) a complete unit cell is shown.The unit cell presented in (c) and (d) has an edge length of 2.01 µm filled with 89 cones corresponding to a cone density of 21.99 cones per µm 2 .

Fig. 4 .
Fig. 4. (a) The simulated reflectance and transmittance of a 250 µm thick silicon wafer with black silicon at the front side agrees well with the measured data.(b) Simulated and measured absorptance compared to the absorptance of a graded index layer and an ideal diffuse scatterer (Yablonovitch limit).

#251492 Received 7 Fig. 5 .
Fig. 5. Angular distribution at a wavelength of 1.04 µm for two unit cell sizes.In the case of the unit cell size of 1.40 µm, 15 structures were considered for averaging.For the unit cell size of 2.01 µm, 10 structures were averaged.The error bars represent the standard deviation.

Fig. 6 .
Fig. 6.Simulated angular distribution of light in the silicon bulk after passage through a black silicon structure for different wavelengths compared to the distribution of a Lambertian scatterer.Three different unit cell sizes were considered, each with 20 random structures.

Fig. 7 .
Fig. 7. Difference between simulated and measured transmittance for different periods.

Fig. 8 .
Fig. 8. Angular distributions for a period of 1.49 µm at wavelengths of 1.04 µm and 1.08 µm.At 1.04 µm three polar angles associated with propagating modes are smaller than the critical angle marked by the vertical line.At a wavelength of 1.08 µm only two polar angles corresponding to propagating modes are within the loss cone.