E ﬀ ects of Incoherent Front Cover Glass on Current-Voltage Characteristics of Cu(In,Ga)Se 2 Solar Cells: Investigation into Calculation Accuracy for Cover Glass Modeled as Optically Coherent or Incoherent

: We numerically investigate the e ﬀ ects of an incoherent front cover glass on the current–voltage (J–V) characteristics of a Cu(In,Ga)Se 2 (CIGS) solar cell using an integrated optoelectronic model. A 3-mm cover glass—the thickness of which was larger than the coherence length of sunlight—was incoherently modeled based on the equispaced thickness averaging method, where coherent simulation results of the wave equation were averaged over a set of equispaced phase thicknesses. The changes in optical power dissipation, absorptivity and electron–hole pair generation rate were calculated depending on the variation of the equispaced phase thickness. The calculation results of the J–V curves were obtained through numerical solutions of the coupled Poisson and continuity equations. By comparing the J–V curves calculated between coherently and incoherently modeled cover glass, we obtained a maximum ± 0.54% deviation of the short-circuit current density. This demonstrates that the front cover glass should be modeled as optically incoherent to improve the calculation accuracy of the electrical J–V curves as well as the optical absorption characteristics in the optoelectronic modeling of CIGS solar cells.


Introduction
Thin-film solar cells based on a Cu(In,Ga)Se 2 (CIGS) absorption material have been intensively investigated as an auspicious alternative to traditional silicon-based solar cells owing to their low cost and high efficiency [1][2][3]. Because of their high optical absorption coefficient, CIGS solar cells with a wide absorption spectrum and high stability reached a power conversion efficiency of 19.9% in 2008 [1] and beyond 20% in more recent years [2,3]. In general, a CIGS solar cell consists of multiple sub-micrometer layers, which range from a metal back contact to transparent conducting oxide (TCO) top contact layers. In the module configuration, millimeter-thick front encapsulation and cover glass layers are added on the TCO top contact layer to implement glass-glass encapsulating technology, which can provide a higher mechanical durability and operational stability [4].
Because the coherence length of sunlight is approximately 0.6 µm [5] (preserving the phase information of the optical field), the multilayer structure of a CIGS solar cell module can be treated as mixed coherent and incoherent layers in optical modeling. Inside the coherent multilayer structure, the thickness of which is much less than the coherence of sunlight, the optical interference effect plays an important role in the optical characteristics of the reflectivity and absorptivity [6]. By contrast, light experiences a random phase change at each coherence length, such that the optical interference effect, on average, disappears in the incoherent layer [7,8]. Thus, this mixed coherent-incoherent multilayer structure of a CIGS solar cell greatly affects the absorption characteristics and has to be thoroughly considered in optoelectronic modeling.
There have been several studies conducted on the optoelectronic modeling of a CIGS solar cell, which have provided the current-voltage (J-V) characteristics obtained through the combination of optical and electrical simulations [9][10][11][12][13]. However, the previous optoelectronic modeling results of CIGS solar cells have not included the incoherent front encapsulation or cover glass layers in the optical modeling. Instead, they assumed sunlight to be directly incident from the air into the multilayer structure of the CIGS solar cell. The analytical formulation of a 2 × 2 transfer matrix method (TMM) based on the thin-film optics was used to calculate the optical interference effect within the coherent multilayer structure, which ranged from the TCO top contact layer to the metal back contact layer [11][12][13]. Thus, the previous optoelectronic modeling methods were insufficient when considering the realistic aspect of the optical and electrical characteristics of a CIGS solar cell because the optical effect of the incoherent front layer was not considered.
According to thin-film optics, the absorption spectrum of the mixed coherent-incoherent layer structure of organic solar cells was calculated based on the generalized transfer matrix method (GTMM). Here, the optical behavior of the electric field amplitude in the coherent multilayer was described by the 2 × 2 TMM, and the optical behavior of the electric field intensity in the incoherent layer was modeled by the 2 × 2 intensity matrix formulation [14,15]. In numerical methods, such as rigorous coupled-wave analysis (RCWA) and finite element method (FEM), coherent simulation results (obtained through the numerical solution of the wave equation in a small mesh or grid) were averaged out to provide the incoherent characteristics of light propagation in a mixed coherent-incoherent multilayer structure [16][17][18]. The effects of the incoherent front layer on the reflection spectrum of the CIGS solar cell were calculated based on the spectral averaging method in combination with the RCWA [16,17].
We recently used the equispaced thickness averaging method (ETAM) to efficiently calculate the effect of the incoherent front layer on the reflection spectra of a CIGS solar cell [19]. In the ETAM, coherent simulation results (obtained based on a set of equispaced phase thicknesses added to the initial thickness of the incoherent layer) were averaged to optically model the mixed coherent-incoherent multilayer structure. However, the above-mentioned optical modeling methods, which can only consider the effects of the incoherent layer on the absorption characteristics, were also insufficient to obtain the realistic J-V characteristics of a CIGS solar cell owing to their limitations in the electrical model. In addition, it has yet to be investigated how the optical effect of the front incoherent layer affects the J-V characteristics of a CIGS solar cell.
In this study, we investigated the effects of the incoherent front cover glass on the J-V characteristic of a CIGS solar cell. An incoherent cover glass was optically modeled as a 1-µm initial layer plus a set of equispaced phase layers (EPLs), the thicknesses of which are determined based on the ETAM. Both the optical and electrical characteristics of a CIGS solar cell were numerically calculated using an FEM-based integrated optoelectronic model. The changes in the optical power dissipation, absorptivity, and electron-hole pair generation rates were calculated depending on the variation of the equispaced phase thickness of the EPL. To obtain the ETAM-applied calculation, which includes the incoherent property of a thick cover glass layer, coherent optical simulation results were averaged over 10 equispaced phase thicknesses of the EPL. By comparing the J-V characteristics of the CIGS solar cell calculated between coherently and incoherently modeled cover glass layers, we calculated the percent of deviation of the short-circuit current density (J SC ) and open-circuit voltage (V OC ), which can be considered to be the inaccuracy of the calculation when a thick cover glass layer is inappropriately modeled as a coherent layer.  Figure 1 shows a schematic diagram of the multilayer structure of the CIGS solar cell used in this study [17,19]. It comprises a cover glass (3 mm) at the top, aluminum-doped zinc oxide (AZO) (50 nm) as a TCO, zinc sulfide (ZnS) (10 nm) as an n-doped layer, CIGS (2 µm) as a p-doped layer and molybdenum (Mo) (400 nm) as a bottom contact layer. The bottom glass substrate is omitted because it does not affect the optical and electrical characteristics of the CIGS solar cell.

Optical Modeling
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 16 Figure 1 shows a schematic diagram of the multilayer structure of the CIGS solar cell used in this study [17,19]. It comprises a cover glass (3 mm) at the top, aluminum-doped zinc oxide (AZO) (50 nm) as a TCO, zinc sulfide (ZnS) (10 nm) as an n-doped layer, CIGS (2 μm) as a p-doped layer and molybdenum (Mo) (400 nm) as a bottom contact layer. The bottom glass substrate is omitted because it does not affect the optical and electrical characteristics of the CIGS solar cell. Figure 1. Schematic diagram of multilayer structure of the Cu(In,Ga)Se2 (CIGS) solar cell used in this study. The 3-mm cover glass, the thickness of which is larger than the coherence length of sunlight (~0.6 μm), is considered as incoherent, and the remaining layers are modeled as coherent. Because the 3-mm cover glass has no absorption, the incoherent cover glass can be modeled as a 1-μm initial cover glass layer plus an equispaced phase layers (EPLs) whose thickness tb is determined based on Equation (17) when the equispaced thickness averaging method (ETAM) is applied to calculate the effect of the incoherent front glass layer in the simulation.
Sunlight is assumed to have a spectral irradiance of air mass (AM) 1.5 sunlight λ ( ) S with a total power of S0 = 100 mW/cm 2 . Sunlight is assumed to be normally incident from the air into the CIGS solar cell, which can be optically modeled as a mixed incoherent-coherent multilayer structure. The 3-mm cover glass, the thickness of which is larger than the coherence length of sunlight (~0.6 μm), is considered optically incoherent. In optical modeling based on the ETAM, the thickness of the transparent incoherent layer can be converted into the corresponding optical phase change, which has a period of 2π. Thus, the thickness of the transparent cover glass used in the simulation does not require an actual thickness of 3 mm [18,19]. To reduce the computation time, we use a 1-μm-thick initial cover glass layer in addition to the EPL, the role of which is to efficiently model a random phase change of light in the incoherent layer through a set of equispaced phase thicknesses in the EPL [18]. The remaining four layers were treated as coherent layers. Nearly all incident lights were absorbed around the interface between the ZnS and CIGS layers and there is no optical interference effect on the light absorption characteristics of the CIGS layer, which will be explained later. Thus, the CIGS layer can be treated as optically coherent although its 2-μm thickness is greater than the coherence length of sunlight.
In the optical modeling, we assume that all materials used in the CIGS solar cell were nonmagnetic, linear, isotropic and homogeneous media with an instantaneous response. If ˆx z r xa za = +  indicates the position on the x-z plane in Figure 1, the wave equation can be written as is the electric field amplitude at a wavelength of λ, and 0 2 k π λ = is the wave number of the free space. The complex refractive index of the materials is expressed as n n jκ = −  , where n and κ were the refractive index and extinction coefficient, respectively. If tb represents the Figure 1. Schematic diagram of multilayer structure of the Cu(In,Ga)Se 2 (CIGS) solar cell used in this study. The 3-mm cover glass, the thickness of which is larger than the coherence length of sunlight (~0.6 µm), is considered as incoherent, and the remaining layers are modeled as coherent. Because the 3-mm cover glass has no absorption, the incoherent cover glass can be modeled as a 1-µm initial cover glass layer plus an equispaced phase layers (EPLs) whose thickness t b is determined based on Equation (17) when the equispaced thickness averaging method (ETAM) is applied to calculate the effect of the incoherent front glass layer in the simulation.
Sunlight is assumed to have a spectral irradiance of air mass (AM) 1.5 sunlight S(λ) with a total power of S 0 = 100 mW/cm 2 . Sunlight is assumed to be normally incident from the air into the CIGS solar cell, which can be optically modeled as a mixed incoherent-coherent multilayer structure. The 3-mm cover glass, the thickness of which is larger than the coherence length of sunlight (~0.6 µm), is considered optically incoherent. In optical modeling based on the ETAM, the thickness of the transparent incoherent layer can be converted into the corresponding optical phase change, which has a period of 2π. Thus, the thickness of the transparent cover glass used in the simulation does not require an actual thickness of 3 mm [18,19]. To reduce the computation time, we use a 1-µm-thick initial cover glass layer in addition to the EPL, the role of which is to efficiently model a random phase change of light in the incoherent layer through a set of equispaced phase thicknesses in the EPL [18]. The remaining four layers were treated as coherent layers. Nearly all incident lights were absorbed around the interface between the ZnS and CIGS layers and there is no optical interference effect on the light absorption characteristics of the CIGS layer, which will be explained later. Thus, the CIGS layer can be treated as optically coherent although its 2-µm thickness is greater than the coherence length of sunlight.
In the optical modeling, we assume that all materials used in the CIGS solar cell were non-magnetic, linear, isotropic and homogeneous media with an instantaneous response. If → r = xâ x + zâ z indicates the position on the x-z plane in Figure 1, the wave equation can be written as where → E( → r , λ) is the electric field amplitude at a wavelength of λ, and k 0 = 2π/λ is the wave number of the free space. The complex refractive index of the materials is expressed as n = n − jκ, where n and κ were the refractive index and extinction coefficient, respectively. If t b represents the equispaced phase thickness of the EPL in Figure 1, the numerical solution of the wave equation in Equation (1) depends on the value of t b . The time-averaged Poynting vector → S ( → r , λ, t b ), indicating the amount of the optical power flow, is given by where is the complex conjugate of the magnetic field amplitude. The amount of optical absorption is described by the optical power dissipation, which is given by The absorptivity can then be obtained by taking the integration over the optical power dissipation within the entire region of the CIGS solar cell in reference to the wavelength-dependent incident power from the uppermost boundary, which corresponds to the solar irradiance of AM 1.5 sunlight. Thus, the absorptivity is written as According to the ETAM, the effect of the incoherent layer can be modeled by averaging the coherent simulation results over a set of equispaced phase thicknesses [18,19]. In the coherent simulations, the set of equispaced phase thicknesses of the EPL used in the cover glass layer is matched with the equally spaced optical phase changes ranging between 0 and 2π, which provide different optical calculation results of the CIGS solar cell. However, if the cover glass is appropriately modeled as optically incoherent, all calculation results with respect to the set of equispaced phase thicknesses should be equal to the ETAM-applied incoherent calculation result irrespective of these equispaced phase thicknesses. Thus, coherent simulation results over the set of equispaced phase thicknesses of the EPL can be interpreted as the variation of the calculation inaccuracy when an incoherent cover glass layer is inappropriately modeled as optically coherent. Then, the ETAM-applied optical power dissipation and absorptivity, which include the effect of the incoherent cover glass layer, can then be expressed as where N indicates the number of equispaced phase thicknesses. Finally, the optical power dissipation Q( → r , λ, t b ) is related with the generation rate of the electron-hole pair, which will be used as one of the input parameters in the continuity equation of the electrical modeling. The electron-hole-pair generation rate G( → r , t b ), depending on the thickness of the EPL, is given by Appl. Sci. 2020, 10, 3312

of 16
where h is the Plank constant and c is the speed of light in a free space. When the optically incoherent characteristics of the thick cover glass were considered, the ETAM-applied generation rate of the electron-hole pair G ETAM ( → r ) is written as

Electrical Modeling
Electrical modeling is applied based on the drift-diffusion model, which is widely used in the electrical modeling of CIGS solar cells [9]. The drift-diffusion model consists of coupled Poisson and continuity equations for the electrons and holes. First, the electron current density J n ( → r ) and hole current density J p ( → r ), composed of the drift and diffusion components, were where q is the basic electric charge and V( → r ) is the electric potential. The terms n( → r ) and p( → r ) indicate the electron and hole concentrations, respectively. When µ n and µ p represent the electron and hole mobilities, D n and D p indicate the diffusion coefficients for the electrons and holes, which were obtained based on the Einstein relation D n(p) = µ n(p) k B T/q. Here, k B is the Boltzmann constant and T is the absolute temperature. The Poisson equation is expressed as where ε 0 is the electric permittivity in a free space and ε r is the dielectric constant, whereas the acceptor and donor densities were N − A ( → r ) and N + D ( → r ). Finally, the continuity equation is given by where G( → r ) indicates the generation rate of the electron-hole pair, which is determined through Equations (7) or (8) depending on whether the thick cover glass is considered to be coherent (G( → r , t b )) or incoherent (G ETAM ( → r )), respectively. The term R( → r ) represents the recombination rate of the electron-hole pair, where we include the non-radiative Shockley-Read-Hall (SRH) and radiative recombination processes through R( The non-radiative SRH recombination rate is expressed as where n i ( → r ) is the intrinsic carrier concentration when τ p ( → r ) and τ n ( → r ) represent the carrier lifetimes of the electron and hole, respectively. The intrinsic carrier concentration is different in each material and can be obtained by Appl. Sci. 2020, 10, 3312 6 of 16 where N c ( → r ) and N v ( → r ) were the effective density of states in the conduction and valence bands, respectively, when E g ( → r ) represents the band gap energy of the material. The electron and hole concentrations at the trap energy level of E t (z) were expressed as n 1 ( the intrinsic Fermi energy level. In addition, the radiative recombination rate is written as where R B is the coefficient of the radiative recombination.

Simulation Results
The optical modeling of the CIGS solar cell in Figure 1 was conducted based on the two-dimensional (2D) FEM, which was implemented using the RF module of COMSOL Multiphysics [20]. The complex refractive index spectra of the materials, which were taken from [21][22][23], are shown in Figure 2. In the FEM simulation, wave equations shown in Equation (1) were numerically solved at very small meshes of the calculation domains with appropriate boundary conditions. We use the port boundary condition at the uppermost boundary of the air domain to generate a plane wave with a wavelength of λ, which was normally injected into the multilayer structure of the CIGS solar cell. In addition, the Floquet periodic boundary conditions were used at the leftmost and rightmost boundaries of the simulation domain.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 16 represents the intrinsic Fermi energy level. In addition, the radiative recombination rate is written as where B R is the coefficient of the radiative recombination.

Simulation Results
The optical modeling of the CIGS solar cell in Figure 1 was conducted based on the twodimensional (2D) FEM, which was implemented using the RF module of COMSOL Multiphysics [20]. The complex refractive index spectra of the materials, which were taken from [21][22][23], are shown in Figure 2. In the FEM simulation, wave equations shown in Equation (1) were numerically solved at very small meshes of the calculation domains with appropriate boundary conditions. We use the port boundary condition at the uppermost boundary of the air domain to generate a plane wave with a wavelength of λ, which was normally injected into the multilayer structure of the CIGS solar cell. In addition, the Floquet periodic boundary conditions were used at the leftmost and rightmost boundaries of the simulation domain. If the refractive index of the transparent incoherent cover glass was ng, the equispaced phase thickness of the EPL can be determined by When the ETAM is applied to the incoherent modeling of a planar organic solar cell, the accuracy of the calculation results converge if the number of the equispaced phase thicknesses is greater than five [18]. In this calculation, we set the number of equispaced phase thicknesses to 10. Hence, the corresponding equispaced phase thicknesses of the EPL were  Figure 3 shows the calculated 2D spatial profiles of the electric field amplitude, time-averaged Poynting vector, and optical power dissipation when the equispaced phase thickness of the EPL was set to t1 = 0 nm. All calculation results were obtained at a wavelength of λ = 400 nm, where the extinction coefficient of the CIGS material was the largest, as shown in Figure 2b. In Figure 3, all calculation results depend only on the longitudinal z-direction and have no spatial dependence of the transverse x-direction, which results from the Floquet periodic boundary condition applied to the If the refractive index of the transparent incoherent cover glass was n g , the equispaced phase thickness of the EPL can be determined by When the ETAM is applied to the incoherent modeling of a planar organic solar cell, the accuracy of the calculation results converge if the number of the equispaced phase thicknesses is greater than five [18]. In this calculation, we set the number of equispaced phase thicknesses to 10. Hence, the corresponding equispaced phase thicknesses of the EPL were Figure 3 shows the calculated 2D spatial profiles of the electric field amplitude, time-averaged Poynting vector, and optical power dissipation when the equispaced phase thickness of the EPL was set to t 1 = 0 nm. All calculation results were obtained at a wavelength of λ = 400 nm, where the extinction coefficient of the CIGS material was the largest, as shown in Figure 2b. In Figure 3, all calculation results depend only on the longitudinal z-direction and have no spatial dependence of the transverse x-direction, which results from the Floquet periodic boundary condition applied to the leftmost and rightmost boundaries in the optical simulation domain. In Figure 3b, the time-averaged Poynting vector remains unchanged at the transparent cover glass layer, but continues decreasing at the AZO, ZnS and CIGS layers owing to the gradual light absorption at these layers. Because of the very high absorption coefficient at the 2-µm thick CIGS layer, the optical power dissipation in Figure 3c was significant near the ZnS/CIGS interface and there was no light absorption near the CIGS/Mo interface. Hence, the CIGS layer can be optically modeled as a coherent layer although its 2-µm thickness was larger than the coherence length of sunlight.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 16 leftmost and rightmost boundaries in the optical simulation domain. In Figure 3b, the time-averaged Poynting vector remains unchanged at the transparent cover glass layer, but continues decreasing at the AZO, ZnS and CIGS layers owing to the gradual light absorption at these layers. Because of the very high absorption coefficient at the 2-μm thick CIGS layer, the optical power dissipation in Figure  3c was significant near the ZnS/CIGS interface and there was no light absorption near the CIGS/Mo interface. Hence, the CIGS layer can be optically modeled as a coherent layer although its 2-μm thickness was larger than the coherence length of sunlight. To investigate how the light absorption within the CIGS solar cell was changed depending on the equispaced phase thickness of the EPL, the spatial profiles of the optical power dissipations were calculated at 10 different equispaced phase thicknesses of t1, t2, ⋅⋅⋅ and t10 when the wavelengths of sunlight were 400 and 800 nm. Because the 2D calculation results obtained through the optical modeling, as shown in Figure 3, depend only on the longitudinal z-direction owing to the periodic boundary condition applied to the transverse x-direction, the calculated spatial profiles in Figure 4 were plotted in the longitudinal z-direction. In Figure 4a, the overall exponentially decaying shapes of the optical absorption calculated at different thicknesses of the EPL were comparable to each other except for their relative magnitudes. This was ascribed to the difference in the optical phase change induced by the respective thickness of the EPL in the cover glass layer. The ETAM-applied optical power dissipation, which includes the optically incoherent characteristics of the cover glass layer, lies in the middle of the 10 coherent calculation results. If the cover glass was considered optically incoherent, all the spatial profiles of the optical power dissipations calculated at 10 different equispaced phase thicknesses will be the same as the ETAM-applied calculation result. By comparing the aggregate sums of the optical power dissipations distributed at the CIGS layer between the coherent and ETAM-applied models, we obtain a deviation within the range of −9.65% to +10.22%, which can be considered as the inaccuracy of the calculation when a thick cover glass was inappropriately modeled as the coherent layer.
The variations in the optical power dissipations depending on the equispaced phase thickness of the EPL were also calculated in Figure 4b when the wavelength of sunlight was 800 nm. The overall shapes of the calculation results in Figure 4b are comparable to those in Figure 4a. When the aggregate sums of the optical power dissipations positioned at the CIGS layer were compared between the coherent and ETAM-applied models, a deviation of −9.60% to +10.17% was obtained. It was noticeable that the calculated optical power dissipations in Figure 4b decrease more slowly from the ZnS/CIGS interface toward the CIGS layer than the calculated optical power dissipations in Figure 4a. This can To investigate how the light absorption within the CIGS solar cell was changed depending on the equispaced phase thickness of the EPL, the spatial profiles of the optical power dissipations were calculated at 10 different equispaced phase thicknesses of t 1 , t 2 , ··· and t 10 when the wavelengths of sunlight were 400 and 800 nm. Because the 2D calculation results obtained through the optical modeling, as shown in Figure 3, depend only on the longitudinal z-direction owing to the periodic boundary condition applied to the transverse x-direction, the calculated spatial profiles in Figure 4 were plotted in the longitudinal z-direction. In Figure 4a, the overall exponentially decaying shapes of the optical absorption calculated at different thicknesses of the EPL were comparable to each other except for their relative magnitudes. This was ascribed to the difference in the optical phase change induced by the respective thickness of the EPL in the cover glass layer. The ETAM-applied optical power dissipation, which includes the optically incoherent characteristics of the cover glass layer, lies in the middle of the 10 coherent calculation results. If the cover glass was considered optically incoherent, all the spatial profiles of the optical power dissipations calculated at 10 different equispaced phase thicknesses will be the same as the ETAM-applied calculation result. By comparing the aggregate sums of the optical power dissipations distributed at the CIGS layer between the coherent and ETAM-applied models, we obtain a deviation within the range of −9.65% to +10.22%, which can be considered as the inaccuracy of the calculation when a thick cover glass was inappropriately modeled as the coherent layer. be explained by the fact that the extinction coefficient of the CIGS layer, which can be seen in Figure  2b, decreases as the wavelength increases from 400 to 800 nm.  Figure 5 shows the calculation results of the absorptivity spectra depending on 10 different equispaced phase thicknesses of t1, t2, ⋅⋅⋅ and t10 along with the ETAM-applied absorptivity. According to Equations (4) and (6), the absorptivity corresponds to the accumulated sum of the optical power dissipation in the longitudinal z-direction at a given wavelength. As indicated in Figure 5, the absorptivity spectrum calculated at each equispaced phase thickness shows a respective ripple. This results from the fact that the wave vector of 0 2 k π λ = and the complex refractive index in Figure 2 depend on the wavelength of light, which causes a wavelength-dependent optical phase change and an optical interference effect. Owing to the additional optical phase change induced by the different thicknesses of the EPL tb in Figure 1, the spectral positions of the peaks and valleys were shifted in the absorptivity spectra calculated at different equispaced phase thicknesses of t1, t2, ⋅⋅⋅ and t10. However, these ripples in the absorptivity spectra disappear in the calculation results obtained by the ETAM, where the calculation results of 10 coherently calculated absorptivity spectra were averaged to include the optically incoherent characteristics of the cover glass layer. For comparison, the absorptivity was also calculated based on the GTMM, which was considered another analytical model to include the incoherent effect of the cover glass layer. The numerical results based on the ETAM were in agreement with the analytical results based on the GTMM, which verifies the validity of the calculation results obtained using the ETAM. In addition, the ETAM-applied absorptivity spectrum has the average value among the 10 coherently modeled absorptivity spectra. At a wavelength of 400 nm, as shown in Figure 5, the variation of the coherently modeled absorptivity in reference to the ETAM-applied absorptivity was observed to be within the range of −9.52% to +10.44%, which was close to the range of the deviation of the optical power dissipation shown in Figure 4a.   Figure 4. Calculated spatial profiles of the optical power dissipations in the longitudinal z-direction Q(z, λ, t b ) with respect to 10 different equispaced phase thicknesses of t 1 , t 2 , ··· and t 10 . The calculation results were obtained at wavelengths of (a) 400 and (b) 800 nm. To consider the effect of the incoherent cover glass, the ETAM-applied optical power dissipation Q ETAM (z, λ), obtained by averaging the 10 coherent calculation results, was designated by the red dotted line.
The variations in the optical power dissipations depending on the equispaced phase thickness of the EPL were also calculated in Figure 4b when the wavelength of sunlight was 800 nm. The overall shapes of the calculation results in Figure 4b are comparable to those in Figure 4a. When the aggregate sums of the optical power dissipations positioned at the CIGS layer were compared between the coherent and ETAM-applied models, a deviation of −9.60% to +10.17% was obtained. It was noticeable that the calculated optical power dissipations in Figure 4b decrease more slowly from the ZnS/CIGS interface toward the CIGS layer than the calculated optical power dissipations in Figure 4a. This can be explained by the fact that the extinction coefficient of the CIGS layer, which can be seen in Figure 2b, decreases as the wavelength increases from 400 to 800 nm. Figure 5 shows the calculation results of the absorptivity spectra depending on 10 different equispaced phase thicknesses of t 1 , t 2 , ··· and t 10 along with the ETAM-applied absorptivity. According to Equations (4) and (6), the absorptivity corresponds to the accumulated sum of the optical power dissipation in the longitudinal z-direction at a given wavelength. As indicated in Figure 5, the absorptivity spectrum calculated at each equispaced phase thickness shows a respective ripple. This results from the fact that the wave vector of k 0 = 2π/λ and the complex refractive index in Figure 2 depend on the wavelength of light, which causes a wavelength-dependent optical phase change and an optical interference effect. Owing to the additional optical phase change induced by the different thicknesses of the EPL t b in Figure 1, the spectral positions of the peaks and valleys were shifted in the absorptivity spectra calculated at different equispaced phase thicknesses of t 1 , t 2 , ··· and t 10 . However, these ripples in the absorptivity spectra disappear in the calculation results obtained by the ETAM, where the calculation results of 10 coherently calculated absorptivity spectra were averaged to include the optically incoherent characteristics of the cover glass layer. For comparison, the absorptivity was also calculated based on the GTMM, which was considered another analytical model to include the incoherent effect of the cover glass layer. The numerical results based on the ETAM were in agreement with the analytical results based on the GTMM, which verifies the validity of the calculation results obtained using the ETAM. In addition, the ETAM-applied absorptivity spectrum has the average value among the 10 coherently modeled absorptivity spectra. At a wavelength of 400 nm, as shown in Figure 5, the variation of the coherently modeled absorptivity in reference to the ETAM-applied absorptivity was observed to be within the range of −9.52% to +10.44%, which was close to the range of the deviation of the optical power dissipation shown in Figure 4a. Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 16 The final step of the optical modeling was to calculate the spatial profile of the generation rate of the electron-hole pair, which was one of the input parameters of the continuity equation in the electrical modeling. According to Equations (7) and (8), the generation rate of the electron-hole pair was obtained by integrating the optical power dissipation weighted by the spectral irradiance of AM 1.5 sunlight, as shown in Figure 6a. Figure Figure 6b were comparable to those of the optical power dissipations in Figure 4. In comparison with the aggregate sums of the generation rates between the coherent and ETAM-applied models, the deviation was calculated to be between −0.67% and +0.9%. This value was much smaller than the maximal absolute deviations calculated under the optical power dissipations shown in Figure 4a (10.22%) and the absorptivity in Figure 5 (10.44%), which results from the significant dependence of both the wavelength of sunlight and the thickness of the EPL. By contrast, the generation rate corresponds to the accumulated sum of the optical power dissipation within the wavelength region ranging between 300 and 1100 nm. Thus, the generation rate calculated at each equispaced phase thickness was less sensitive to the change in wavelength, which causes a reduction of the percentage variation of the electron-hole pair generation rate. The final step of the optical modeling was to calculate the spatial profile of the generation rate of the electron-hole pair, which was one of the input parameters of the continuity equation in the electrical modeling. According to Equations (7) and (8), the generation rate of the electron-hole pair was obtained by integrating the optical power dissipation weighted by the spectral irradiance of AM 1.5 sunlight, as shown in Figure 6a. Figure 6b shows the calculated spatial profiles of the generation rate in the longitudinal z-direction G(z, t b ) with respect to 10 different equispaced phase thicknesses of t 1 , t 2 , ··· and t 10 . To include the incoherent effect of the thick cover glass, the ETAM-applied generation rate in the longitudinal z-direction G ETAM (z) was also plotted in the red dotted line. The overall exponentially decaying shapes of the generation rates in Figure 6b were comparable to those of the optical power dissipations in Figure 4. In comparison with the aggregate sums of the generation rates between the coherent and ETAM-applied models, the deviation was calculated to be between −0.67% and +0.9%. This value was much smaller than the maximal absolute deviations calculated under the optical power dissipations shown in Figure 4a (10.22%) and the absorptivity in Figure 5 (10.44%), which results from the significant dependence of both the wavelength of sunlight and the thickness of the EPL. By contrast, the generation rate corresponds to the accumulated sum of the optical power dissipation within the wavelength region ranging between 300 and 1100 nm. Thus, the generation rate calculated at each equispaced phase thickness was less sensitive to the change in wavelength, which causes a reduction of the percentage variation of the electron-hole pair generation rate.
To obtain the J-V characteristics of a CIGS solar cell, the coupled drift-diffusion and continuity equations, shown in Equations (9)-(16), were numerically solved through the use of the semiconductor module of COMSOL multiphysics [20]. In the 2D simulation domain of the CIGS solar cell in Figure 1, we only consider the regions containing the ZnS and CIGS layers for the electrical modeling. Both the bottom cathode terminal of the Mo/CIGS interface and the top anode terminal of the AZO/ZnS interface were assumed to be ideal Schottky contacts. In accordance with the Floquet periodic boundary condition in the optical modeling, the periodic conditions were used at the leftmost and rightmost boundaries of the simulation domain in the electrical modeling. The FEM-based electrical modeling was performed under the same mesh configuration as the FEM-based optical modeling to minimize the error of a mismatch. The simulation parameters used in the electrical modeling are shown in Table 1 [24]. Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 16 To obtain the J-V characteristics of a CIGS solar cell, the coupled drift-diffusion and continuity equations, shown in Equations (9)-(16), were numerically solved through the use of the semiconductor module of COMSOL multiphysics [20]. In the 2D simulation domain of the CIGS solar cell in Figure 1, we only consider the regions containing the ZnS and CIGS layers for the electrical modeling. Both the bottom cathode terminal of the Mo/CIGS interface and the top anode terminal of the AZO/ZnS interface were assumed to be ideal Schottky contacts. In accordance with the Floquet periodic boundary condition in the optical modeling, the periodic conditions were used at the leftmost and rightmost boundaries of the simulation domain in the electrical modeling. The FEMbased electrical modeling was performed under the same mesh configuration as the FEM-based optical modeling to minimize the error of a mismatch. The simulation parameters used in the electrical modeling are shown in Table 1 [24].    Figure 7a shows the calculated J-V characteristics of a CIGS solar cell depending on the equispaced phase thickness of the EPL. In addition, the ETAM-applied J-V curve, where the ETAM-applied generation rate was used in the electrical modeling, was also plotted along the red dotted line. None of the eleven calculation results of the J-V curves can be distinguished in Figure 7a due to their extremely small differences at this J-V scale. In the ETAM-applied J-V curve, the calculated values of J SC and V OC were −34.0 mA/cm 2 and 0.81 V, respectively, which were close to other calculation results of J SC and V OC in a similar CIGS solar cell structure [25]. A close view of the J-V characteristics at near the zero applied voltage is shown in Figure 7b. The deviation of J SC calculated between the coherent and ETAM-applied models lies within the range of −0.54% to +0.54%, which was close to the −0.67% to +0.9% deviation obtained by the aggregate sums of the electron-hole pair generation rates between the coherent and ETAM-applied models. In addition, the maximal deviation of V OC was as small as ±0.017%, which was ascribed to the fact that V OC has a logarithmic dependence on the carrier density [26]. characteristics at near the zero applied voltage is shown in Figure 7b. The deviation of JSC calculated between the coherent and ETAM-applied models lies within the range of −0.54% to +0.54%, which was close to the −0.67% to +0.9% deviation obtained by the aggregate sums of the electron-hole pair generation rates between the coherent and ETAM-applied models. In addition, the maximal deviation of VOC was as small as ± 0.017%, which was ascribed to the fact that VOC has a logarithmic dependence on the carrier density [26]. The maximal percent deviation of JSC calculated between the coherent and ETAM-applied models was ± 0.54%, the value of which can be understood to be the calculation inaccuracy when the thick cover glass layer was inappropriately modeled as the coherent layer.
It is notable that the ±0.54% maximal percent deviation of JSC can be understood to be the inaccuracy of the calculation when a thick cover glass layer was inappropriately modeled as a coherent layer. An arbitrarily determined thickness of the coherently modeled cover glass can provide a random optical phase change ranging between 0 and 2π. This will result in a different optical interference effect in the multilayer structure of a CIGS solar cell and affect all of the calculation results, such as the optical power dissipation in Figure 4, the generation rate in Figure 6b, and the J-V characteristics in Figure 7. Because 10 equispaced phase thicknesses of the EPL used in the ETAM correspond to equally spaced optical phase changes ranging between 0 and 2π, 10 coherent calculation results of the optical power dissipations, electron-hole pair generation rates and J-V characteristics with respect to the equispaced phase thicknesses of the EPL can be speculated to represent the range of variation of the calculation inaccuracy caused by the equispaced phase thicknesses of the coherently modeled cover glass. By contrast, the optical modeling of the ETAMapplied cover glass provides the exact solution to the optical power dissipation, generation rate and J-V characteristic, despite the thickness of the cover glass layer being arbitrarily determined. Therefore, the front cover glass must be modeled as optically incoherent to improve the calculation accuracy of the J-V and optical absorption characteristics in the optoelectronic modeling of a CIGS solar cell.

Discussion
To investigate which wavelength affects the variation of JSC in Figure 7 further, we calculate the EPL-thickness dependence of absorptivity in the CIGS absorber layer and external quantum efficiency (EQE), which are shown in Figure 8. The absorptivity in the CIGS absorber layer was The maximal percent deviation of J SC calculated between the coherent and ETAM-applied models was ± 0.54%, the value of which can be understood to be the calculation inaccuracy when the thick cover glass layer was inappropriately modeled as the coherent layer.
It is notable that the ±0.54% maximal percent deviation of J SC can be understood to be the inaccuracy of the calculation when a thick cover glass layer was inappropriately modeled as a coherent layer. An arbitrarily determined thickness of the coherently modeled cover glass can provide a random optical phase change ranging between 0 and 2π. This will result in a different optical interference effect in the multilayer structure of a CIGS solar cell and affect all of the calculation results, such as the optical power dissipation in Figure 4, the generation rate in Figure 6b, and the J-V characteristics in Figure 7. Because 10 equispaced phase thicknesses of the EPL used in the ETAM correspond to equally spaced optical phase changes ranging between 0 and 2π, 10 coherent calculation results of the optical power dissipations, electron-hole pair generation rates and J-V characteristics with respect to the equispaced phase thicknesses of the EPL can be speculated to represent the range of variation of the calculation inaccuracy caused by the equispaced phase thicknesses of the coherently modeled cover glass. By contrast, the optical modeling of the ETAM-applied cover glass provides the exact solution to the optical power dissipation, generation rate and J-V characteristic, despite the thickness of the cover glass layer being arbitrarily determined. Therefore, the front cover glass must be modeled as optically incoherent to improve the calculation accuracy of the J-V and optical absorption characteristics in the optoelectronic modeling of a CIGS solar cell.

Discussion
To investigate which wavelength affects the variation of J SC in Figure 7 further, we calculate the EPL-thickness dependence of absorptivity in the CIGS absorber layer and external quantum efficiency (EQE), which are shown in Figure 8. The absorptivity in the CIGS absorber layer was obtained by integrating the optical power dissipation in the longitudinal z-direction over the CIGS absorber layer and was written as The EQE was defined as the ratio of the number of charge carriers collected by the solar cell to the number of incident photons shining on the solar cell outside. In this calculation, the EQE spectra was expressed as where J SC (λ, t b ) indicates the EPL-thickness-dependent spectral short-circuit current density, which was obtained by the numerical solution of the coupled drift-diffusion and continuity equations when the incident optical power was set to the spectral irradiance of AM 1.5 sunlight S(λ). In Figure 8, the ETAM-applied calculation results were obtained by averaging 10 coherent calculation results over the equispaced phase thicknesses. According to the calculated EQE spectra in Figure 8b, the variation of the EQE with respect to the EPL thickness was noticeable around the wavelength range between 600 and 1000 nm.
the number of incident photons shining on the solar cell outside. In this calculation, the EQE spectra was expressed as where λ ( , ) indicates the EPL-thickness-dependent spectral short-circuit current density, which was obtained by the numerical solution of the coupled drift-diffusion and continuity equations when the incident optical power was set to the spectral irradiance of AM 1.5 sunlight λ ( ) S . In Figure 8, the ETAM-applied calculation results were obtained by averaging 10 coherent calculation results over the equispaced phase thicknesses. According to the calculated EQE spectra in Figure 8b, the variation of the EQE with respect to the EPL thickness was noticeable around the wavelength range between 600 and 1000 nm.  Figure 9a shows the ETAM-applied calculation results of absorptivity in the entire multilayer structure, the absorptivity in the CIGS absorber layer and the EQE, where a significant reduction of the EQE was observed within the relatively short (300-450 nm) and long (1000-1100 nm) wavelength ranges. First, the reduction of the EQE between 300 and 450 nm was ascribed to the optical absorption in the AZO-based TCO layer, which has a high extinction coefficient within the short wavelength range, as shown in Figure 2b. Hence, the difference in the absorptivity between the entire multilayer structure and the CIGS absorber layer becomes very high within the wavelength range of 300-450 nm. This fact was also confirmed by the ETAM-applied spatial profiles of the optical power dissipations shown in Figure 9b, where the optical power dissipation at the AZO layer was observed at wavelengths of 300, 350 and 400 nm. Second, the difference between the absorptivity and EQE in Figure 9a increases as the wavelength increases from 500 to 1100 nm. As shown in Figure 9b, this can be explained by the fact that the optical absorption, taking place farther away from the interface between the p-doped CIGS and n-doped ZnS layers at a longer wavelength, contributes less to the electric current at a longer wavelength due to the low diffusion length of the minority electron carriers [26].  Figure 9a shows the ETAM-applied calculation results of absorptivity in the entire multilayer structure, the absorptivity in the CIGS absorber layer and the EQE, where a significant reduction of the EQE was observed within the relatively short (300-450 nm) and long (1000-1100 nm) wavelength ranges. First, the reduction of the EQE between 300 and 450 nm was ascribed to the optical absorption in the AZO-based TCO layer, which has a high extinction coefficient within the short wavelength range, as shown in Figure 2b. Hence, the difference in the absorptivity between the entire multilayer structure and the CIGS absorber layer becomes very high within the wavelength range of 300-450 nm. This fact was also confirmed by the ETAM-applied spatial profiles of the optical power dissipations shown in Figure 9b, where the optical power dissipation at the AZO layer was observed at wavelengths of 300, 350 and 400 nm. Second, the difference between the absorptivity and EQE in Figure 9a increases as the wavelength increases from 500 to 1100 nm. As shown in Figure 9b, this can be explained by the fact that the optical absorption, taking place farther away from the interface between the p-doped CIGS and n-doped ZnS layers at a longer wavelength, contributes less to the electric current at a longer wavelength due to the low diffusion length of the minority electron carriers [26].
In the actual module configuration of the CIGS solar cell, the front encapsulation layer, made of ethylene vinyl acetate (EVA), lies between the front cover glass and the TCO top contact layer [4]. Figure 10 shows a modified schematic diagram of the CIGS solar cell, where a 0.5-mm EVA encapsulation layer was added into the multilayer structure of the CIGS solar cell shown in Figure 1 [17,19]. Because the thickness of the EVA layer (0.5 mm) was much larger than the coherence length of sunlight (~0.6 µm), the EVA layer should also be considered optically incoherent. The EVA layer was assumed to be transparent with a constant refractive index of 1.49 between 300 and 1100 nm [19]. According to the optical modeling based on the ETAM, the incoherent cover glass and EVA layers can be modeled as an initial 1-µm layer plus the EPL layer, respectively. The optically incoherent characteristics of the cover glass/EVA layer can then be considered by averaging the coherent simulation results over all combinations of the equispaced phase thicknesses in two incoherent layers [19]. For instance, the ETAM-applied optical power dissipation, which includes the effect of the incoherent cover glass/EVA layer, can be obtained by where M and N indicate the number of equispaced phase thicknesses in the EVA and cover glass layers, respectively. If the number of equispaced phase thicknesses was set to 10 for both incoherent layers, the ETAM-applied optical power dissipation can be calculated by averaging the coherent calculation results Q( In the actual module configuration of the CIGS solar cell, the front encapsulation layer, made of ethylene vinyl acetate (EVA), lies between the front cover glass and the TCO top contact layer [4]. Figure 10 shows a modified schematic diagram of the CIGS solar cell, where a 0.5-mm EVA encapsulation layer was added into the multilayer structure of the CIGS solar cell shown in Figure 1 [17,19]. Because the thickness of the EVA layer (0.5 mm) was much larger than the coherence length of sunlight (~0.6 μm), the EVA layer should also be considered optically incoherent. The EVA layer was assumed to be transparent with a constant refractive index of 1.49 between 300 and 1100 nm [19]. According to the optical modeling based on the ETAM, the incoherent cover glass and EVA layers can be modeled as an initial 1-μm layer plus the EPL layer, respectively. The optically incoherent characteristics of the cover glass/EVA layer can then be considered by averaging the coherent simulation results over all combinations of the equispaced phase thicknesses in two incoherent layers [19]. For instance, the ETAM-applied optical power dissipation, which includes the effect of the incoherent cover glass/EVA layer, can be obtained by where M and N indicate the number of equispaced phase thicknesses in the EVA and cover glass layers, respectively. If the number of equispaced phase thicknesses was set to 10 for both incoherent layers, the ETAM-applied optical power dissipation can be calculated by averaging the coherent calculation results  Figure 1. Because the thickness of the EVA layer (0.5 mm) was larger than the coherence length of sunlight (~0.6 μm), the EVA layer should also be considered optically incoherent. Because both a cover glass and EVA have no absorption, two incoherent layers can be modeled as a 1-μm initial layer plus the EPL, whose respective thicknesses of ta and tb can be determined based on Equation (17) to include the effect of the incoherent cover glass/EVA layers on the optical characteristics based on the ETAM. Figure 11 shows the calculation results of the ETAM-applied absorptivity spectra and J-V characteristics of a CIGS solar cell having only a cover glass (Figure 1) or a cover glass/EVA as an incoherent layer ( Figure 10). In Figure 11a, the absorptivity spectra calculated by the GTMM show a good agreement with those obtained by the ETAM, which verifies the validity of the calculation results obtained using the ETAM. There was a slight difference between the two calculation results Figure 10. Modified schematic diagram of CIGS solar cell, where a 0.5-mm EVA encapsulation layer was added into the multilayer structure of the CIGS solar cell in Figure 1. Because the thickness of the EVA layer (0.5 mm) was larger than the coherence length of sunlight (~0.6 µm), the EVA layer should also be considered optically incoherent. Because both a cover glass and EVA have no absorption, two incoherent layers can be modeled as a 1-µm initial layer plus the EPL, whose respective thicknesses of t a and t b can be determined based on Equation (17) to include the effect of the incoherent cover glass/EVA layers on the optical characteristics based on the ETAM. Figure 11 shows the calculation results of the ETAM-applied absorptivity spectra and J-V characteristics of a CIGS solar cell having only a cover glass (Figure 1) or a cover glass/EVA as an incoherent layer ( Figure 10). In Figure 11a, the absorptivity spectra calculated by the GTMM show a good agreement with those obtained by the ETAM, which verifies the validity of the calculation results obtained using the ETAM. There was a slight difference between the two calculation results of absorptivity because the difference in the refractive index between the cover glass (1.55) and EVA (1.49) was very small. Hence, the calculated values of J SC between the cover glass/AZO/ZnS/CIGS/Mo (−34 mA/cm 2 ) and cover glass/EVA/AZO/ZnS/CIGS/Mo (−33.91 mA/cm 2 ) structures were very close and deviate by 0.26%. This indicates that the addition of the incoherent EVA encapsulation layer into the multilayer structure of the CIGS solar cell in Figure 1 will provide little change in either the optical or electrical calculation results shown in Figures 3-7, where only a cover glass was considered as the incoherent layer.
was added into the multilayer structure of the CIGS solar cell in Figure 1. Because the thickness of the EVA layer (0.5 mm) was larger than the coherence length of sunlight (~0.6 μm), the EVA layer should also be considered optically incoherent. Because both a cover glass and EVA have no absorption, two incoherent layers can be modeled as a 1-μm initial layer plus the EPL, whose respective thicknesses of ta and tb can be determined based on Equation (17) to include the effect of the incoherent cover glass/EVA layers on the optical characteristics based on the ETAM. Figure 11 shows the calculation results of the ETAM-applied absorptivity spectra and J-V characteristics of a CIGS solar cell having only a cover glass (Figure 1) or a cover glass/EVA as an incoherent layer ( Figure 10). In Figure 11a, the absorptivity spectra calculated by the GTMM show a good agreement with those obtained by the ETAM, which verifies the validity of the calculation results obtained using the ETAM. There was a slight difference between the two calculation results of absorptivity because the difference in the refractive index between the cover glass (1.55) and EVA (1.49) was very small. Hence, the calculated values of JSC between the cover glass/AZO/ZnS/CIGS/Mo (−34 mA/cm 2 ) and cover glass/EVA/AZO/ZnS/CIGS/Mo (−33.91 mA/cm 2 ) structures were very close and deviate by 0.26%. This indicates that the addition of the incoherent EVA encapsulation layer into the multilayer structure of the CIGS solar cell in Figure 1 will provide little change in either the optical or electrical calculation results shown in Figures 3-7, where only a cover glass was considered as the incoherent layer.  Figure 10). The absorptivity spectra calculated using the GTMM show a good agreement with the results obtained using the ETAM. The inset shows a close view of the J-V characteristics at near a zero applied voltage.

Conclusions
We numerically investigated the effect of the incoherent front cover glass on the J-V characteristics of a CIGS solar cell using an FEM-based integrated optoelectronic model. The incoherent property of the 3-mm cover glass was optically modeled as a 1-µm initial layer plus the EPL, the thickness of which was determined based on the ETAM. Various absorption characteristics such as the optical power dissipation, absorptivity and electron-hole pair generation rate were calculated at 10 equispaced phase thicknesses when the cover glass was modeled as the coherent layer. The ETAM-applied calculation results were also obtained by averaging the coherent simulation results over 10 equispaced phase thicknesses. By comparing the J-V characteristics calculated between coherently modeled and ETAM-applied cover glasses, we obtained a maximum ±0.54% deviation of J SC , which can be considered a calculation inaccuracy when a thick cover glass layer was inappropriately modeled as the coherent layer. Therefore, it is important to treat the front cover glass as optically incoherent to improve the calculation accuracy of both the optical absorption and J-V characteristics in the optoelectronic modeling of the CIGS solar cell.