Coupled Optoelectronic Simulation and Optimization of Thin-Film Photovoltaic Solar Cells

TL;DR

This paper presents a comprehensive simulation and optimization framework for thin-film photovoltaic solar cells, integrating optical and electrical modeling to enhance efficiency through advanced computational techniques.

Contribution

It introduces a coupled optoelectronic simulation tool that combines rigorous optical modeling with charge transport analysis, optimizing layer parameters for improved solar cell performance.

Findings

Optimized layer dimensions and bandgaps increase efficiency.

Hybridizable discontinuous Galerkin method accurately models charge transport.

Differential evolution effectively finds optimal design parameters.

Abstract

A design tool was formulated for optimizing the efficiency of inorganic, thin-film, photovoltaic solar cells. The solar cell can have multiple semiconductor layers in addition to antireflection coatings, passivation layers, and buffer layers. The solar cell is backed by a metallic grating which is periodic along a fixed direction. The rigorous coupled-wave approach is used to calculate the electron-hole-pair generation rate. The hybridizable discontinuous Galerkin method is used to solve the drift-diffusion equations that govern charge-carrier transport in the semiconductor layers. The chief output is the solar-cell efficiency which is maximized using the differential evolution algorithm to determine the optimal dimensions and bandgaps of the semiconductor layers.

Tables

Table 1. Table 1: Table of scaling parameters and applicable variables.

Scaling	Symbol	Representative Value	Variable(s)
Length	$L_s$	$𝒪\left(L_z\right)\approx 2\times {10}^{-4}$ cm	$z$
Mobility	${\mu }_s$	$1.5$ cm2 V-1 s-1	${\mu }_{n,p}$
Voltage	${\varphi }_s=V_{\mathrm{th}}$	$0.0259$ V	$\varphi$, ${\varphi }_{n,p}$, $\chi$, $V_{\mathrm{ext}}$,
			${𝖤}_{i,c,F_n,F_p,g,v}$
Density	$N_s$	$\max \left(\left|N_f\right|\right)\approx {10}^{16}$ cm-3	$n$, $p$, $\overline{n}$,
			$n_{0,1}$, $p_{0,1}$,$N_{c,f,v}$
Poisson const.	${\lambda }^2$	$4.1717\times {10}^{-4}$	–
Density rate	$G_s=\frac{{\mu }_s{N}_s{\varphi }_s}{L_s^2}$	$9.6945\times {10}^{21}$ cm-3 s-1	$G$, $R$
Current	$J_s=\frac{q_e{\mu }_s{N}_s{\varphi }_s}{L_s}$	$0.3106$ mA cm-2	$J$, $J_{n,p}$

Table 2. Table 2: Correspondence between the generalized transport system of Eqs. ( 102 )–( 104 ) and Eqs. ( 94 )–( 101 ) of the dimensionless DDE model.

	Eqs. (94) and (96)	Eqs. (95) and (97)	Eqs. (100) and (101)
$𝖩$	$J_n$	$J_p$	$E$
$𝗆$	$n$	$p$	$\varphi$
$𝖹$	$\varphi$	$\varphi$	0
${𝗏}_{𝗆}$	${\varphi }_n$	${\varphi }_p$	0
${𝖼}_1$	$1/{\mu }_n$	$-1/{\mu }_p$	$-1/{\lambda }^2$
${𝖼}_2$	$-1$	$1$	$0$
${𝖼}_3$	0	0	1

Table 3. Table 3: Electrical parameters of CIGS.

Parameter	Symbol	Value
Bandgap	$E_g$	1.3 V
Electron affinity	$\chi$	4.5 V
Conduction-band density of states	$N_c$	$2.22\times {10}^{18}$ cm-3
Valence-band density of states	$N_v$	$1.78\times {10}^{19}$ cm-3
DC relative permittivity	${\epsilon }_{dc}^0$	13.6
Electron mobility	${\mu }_n$	$100$ cm2 V-1s-1
Hole mobility	${\mu }_p$	25 cm2 V-1s-1
SRH lifetime parameter (electrons)	${\tau }_n\left(z\right)$	$1\times {10}^{-9}$s-1
SRH lifetime parameter (holes)	${\tau }_p\left(z\right)$	$1\times {10}^{-9}$s-1

Table 4. Table 4: Baseline parameter settings in Sec.  5 for the HDG simulation of the model CIGS p-i-n solar cell in Fig.  1 (a). Note that N z = N zp + N zi + N zn subscript 𝑁 z subscript 𝑁 zp subscript 𝑁 zi subscript 𝑁 zn N_{\rm z}=N_{\rm zp}+N_{\rm zi}+N_{\rm zn} . The state vector is the vector of all solution values in the electronic step.

Symbol	Value	Description
$P_{\mathrm{deg}}$	variable	Degree of Lobatto polynomials
$N_{\mathrm{zp}}$	variable	Number of points in $p$-type layer
$N_{\mathrm{zi}}$	variable	Number of points in $i$-type layer
$N_{\mathrm{zn}}$	variable	Number of points in $n$-type layer
$I_{\mathrm{deg}}$	variable	Nonlinear term integration degree
${\overline{\tau }}_1$	${10}^{-3}$	First parameter for hybridization
${\overline{\tau }}_2$	${10}^3$	Second parameter for hybridization
$\delta V$	0.01	Maximum step for $V_{\mathrm{ext}}$
${\delta }_{\min }$	${10}^{-2}$	Initial homotopy damping
${\epsilon }_{\mathrm{homo}}$	$1.2$	Homotopy damping relaxation
${\delta }_{\mathrm{num}}$	$1$	No. of homotopy attempts
${\delta }_{\mathrm{inc}}$	10	Homotopy increase on new attempt
$n_{\mathrm{loop}}$	10	Maximum number of iterations for Newton’s method
$E_{\mathrm{rel},\mathrm{tol}}$	${10}^{-4}$	Relative change in state vector allowed
$E_{\mathrm{abs},\mathrm{tol}}$	${10}^{-6}$	Absolute change in state vector allowed
$J_{\mathrm{tol}}$	$0.1$	Relative noise allowed in $J$
$P_{\mathrm{ref}}$	100	Max iteration to find $P_{\max }$
$P_{\mathrm{tol}}$	${10}^{-6}$	Allowed variation around $P_{\max }$

Table 5. Table 5: Baseline parameter settings for DEA optimization in Sec.  7 . These are laxer than those used in Sec.  5 as rapid computation is necessary. Other settings are the same as in Table  4 .

Symbol	Value	Description
$N_t$	10	Fourier order for RCWA
$P_{\mathrm{deg}}$	3	Degree of Lobatto polynomials
$N_{\mathrm{zp}}$	10	Number of points in $p$-type layer
$N_{\mathrm{zi}}$	50	Number of points in $i$-type layer
$N_{\mathrm{zn}}$	10	Number of points in $n$-type layer
$I_{\mathrm{deg}}$	4	Nonlinear term integration degree
${\overline{\tau }}_1$	${10}^{-6}$	First parameter for hybridization
${\overline{\tau }}_2$	${10}^6$	Second parameter for hybridization
${\delta }_{\min }$	${10}^{-2}$	Initial homotopy damping
${\epsilon }_{\mathrm{homo}}$	1.1	Homotopy damping relaxation
${\delta }_{\mathrm{num}}$	40	No. of homotopy attempts
${\delta }_{\mathrm{inc}}$	5	Homotopy increase on new attempt
$n_{\mathrm{loop}}$	40	Maximum number of iterations for Newton’s method
$E_{\mathrm{rel},\mathrm{tol}}$	${10}^{-2}$	Relative change in state vector allowed
$E_{\mathrm{abs},\mathrm{tol}}$	${10}^{-4}$	Absolute change in state vector allowed
$P_{\mathrm{tol}}$	${10}^{-3}$	Allowed variation around $P_{\max }$
$t_{\max }$	300 s	Maximum computer time

Table 6. Table 6: Parameters varied in the optimal design study described in Sec.  6 . The optimal value refers to the value at the end of the DEA computation.

Variable Parameter	Optimal Value
$L_x$	632 nm
$L_g$	167 nm
${\zeta }_g$	0.3
$L_i$	596 nm
${𝖤}_{g,p}$	1.62 V
${𝖤}_{g,i}$	1.467  V
${𝖤}_{g,n}$	1.618 V

Equations

$$\mathcal{X}:\{\underline{r}\equiv (x,y,z)|-\infty <x<\infty ,-\infty <y<\infty ,-L_{\mathrm{m}}<z<L_{\mathrm{t}}\}$$

$${\underline{E}}_{\mathrm{inc}}(x,z,{\lambda }_0)=E_0\frac{a_{\mathrm{p}}\hat{\underline{x}}+a_{\mathrm{s}}\hat{\underline{y}}}{\sqrt{a_{\mathrm{p}}^2+a_{\mathrm{s}}^2}}\exp (-i{k}_{0}z).$$

$${\epsilon }_{\mathrm{rel}}(x,z,{\lambda }_0)$$

$${\beta }_{\mathrm{rel}}(x,z,{\lambda }_0)$$

$$\underline{E}(x,z,{\lambda }_0)$$

$$\underline{H}(x,z,{\lambda }_0)$$

$${\breve{\epsilon }}_{\mathrm{rel}}(z,{\lambda }_0)={\left[{\epsilon }_{\mathrm{rel}}^{(-N_{\mathrm{t}})}(z,{\lambda }_0),{\epsilon }_{\mathrm{rel}}^{(-N_{\mathrm{t}}+1)}(z,{\lambda }_0),...,{\epsilon }_{\mathrm{rel}}^{(N_{\mathrm{t}}-1)}(z,{\lambda }_0),{\epsilon }_{\mathrm{rel}}^{(N_{\mathrm{t}})}(z,{\lambda }_0)\right]}^T,$$

$${\breve{\beta }}_{\mathrm{rel}}(z,{\lambda }_0)={\left[{\beta }_{\mathrm{rel}}^{(-N_{\mathrm{t}})}(z,{\lambda }_0),{\beta }_{\mathrm{rel}}^{(-N_{\mathrm{t}}+1)}(z,{\lambda }_0),...,{\beta }_{\mathrm{rel}}^{(N_{\mathrm{t}}-1)}(z,{\lambda }_0),{\beta }_{\mathrm{rel}}^{(N_{\mathrm{t}})}(z,{\lambda }_0)\right]}^T,$$

$${\breve{e}}_{\sigma }(z,{\lambda }_0)={\left[e_{\sigma }^{(-N_{\mathrm{t}})}(z,{\lambda }_0),e_{\sigma }^{(-N_{\mathrm{t}}+1)}(z,{\lambda }_0),...,e_{\sigma }^{(N_{\mathrm{t}}-1)}(z,{\lambda }_0),e_{\sigma }^{(N_{\mathrm{t}})}(z,{\lambda }_0)\right]}^T,$$

$$\sigma \in \{x,y,z\},$$

$${\breve{h}}_{\sigma }(z,{\lambda }_0)={\left[h_{\sigma }^{(-N_{\mathrm{t}})}(z,{\lambda }_0),h_{\sigma }^{(-N_{\mathrm{t}}+1)}(z,{\lambda }_0),...,h_{\sigma }^{(N_{\mathrm{t}}-1)}(z,{\lambda }_0),h_{\sigma }^{(N_{\mathrm{t}})}(z,{\lambda }_0)\right]}^T$$

$$\sigma \in \{x,y,z\},$$

$$\breve{K}=\mathrm{diag}\left[{\kappa }^{(-N_{\mathrm{t}})},{\kappa }^{(-N_{\mathrm{t}}+1)},...,{\kappa }^{(N_{\mathrm{t}}-1)},{\kappa }^{(N_{\mathrm{t}})}\right]$$

$$\displaystyle{\cal{}T}_{N_{\rm t}}(\xi)=\left[\begin{array}[]{ccccc}\xi^{(0)}&\xi^{(-1)}&...&\xi^{(-2N_{\rm t}+1)}&\xi^{(-2N_{\rm t})}\\[5.0pt] \xi^{(1)}&\xi^{(0)}&...&\xi^{(-2N_{\rm t}+2)}&\xi^{(-2N_{\rm t}+1)}\\[5.0pt] \vdots&\vdots&\ddots&\vdots&\vdots\\[5.0pt] \xi^{(2N_{\rm t}-1)}&\xi^{(2N_{\rm t}-2,\lambda_{\scriptscriptstyle 0})}&...&\xi^{(0)}&\xi^{(-1)}\\[5.0pt] \xi^{(2N_{\rm t})}&\xi^{(2N_{\rm t}-1)}&...&\xi^{(1)}&\xi^{(0)}\end{array}\right]\,.$$

$$\displaystyle\frac{d}{dz}\left[\begin{array}[]{c}\breve{{e}}_{x}(z,\lambda_{\scriptscriptstyle 0})\\[3.0pt] \breve{{h}}_{y}(z,\lambda_{\scriptscriptstyle 0})\end{array}\right]=$$

$$\displaystyle i\left[\begin{array}[]{cccc}\breve{{0}}&\omega\mu_{\scriptscriptstyle 0}\breve{{I}}-\left(\omega\varepsilon_{\scriptscriptstyle 0}\right)^{-1}\breve{{K}}\breve{{N}}(z,\lambda_{\scriptscriptstyle 0})\breve{{K}}\\[3.0pt] \omega\varepsilon_{\scriptscriptstyle 0}\breve{{M}}(z,\lambda_{\scriptscriptstyle 0})&\breve{{0}}\end{array}\right]\left[\begin{array}[]{c}\breve{{e}}_{x}(z,\lambda_{\scriptscriptstyle 0})\\[3.0pt] \breve{{h}}_{y}(z,\lambda_{\scriptscriptstyle 0})\end{array}\right]$$

$$\displaystyle\frac{d}{dz}\left[\begin{array}[]{c}\breve{{e}}_{y}(z,\lambda_{\scriptscriptstyle 0})\\[3.0pt] \breve{{h}}_{x}(z,\lambda_{\scriptscriptstyle 0})\end{array}\right]=$$

$$\displaystyle-i\left[\begin{array}[]{cccc}\breve{{0}}&\omega\mu_{\scriptscriptstyle 0}\breve{{I}}\\[3.0pt] \omega\varepsilon_{\scriptscriptstyle 0}\breve{{M}}(z,\lambda_{\scriptscriptstyle 0})-(\omega\mu_{\scriptscriptstyle 0})^{-1}\breve{{K}}^{2}&\breve{{0}}\end{array}\right]\left[\begin{array}[]{c}\breve{{e}}_{y}(z,\lambda_{\scriptscriptstyle 0})\\[3.0pt] \breve{{h}}_{x}(z,\lambda_{\scriptscriptstyle 0})\end{array}\right]$$

$$\displaystyle\left.\begin{array}[]{l}\breve{{e}}_{z}(z,\lambda_{\scriptscriptstyle 0})=-(\omega\varepsilon_{\scriptscriptstyle 0})^{-1}\breve{{N}}(z,\lambda_{\scriptscriptstyle 0})\textbf{$\breve{{K}}$}\breve{{h}}_{y}(z,\lambda_{\scriptscriptstyle 0})\\[5.0pt] \breve{{h}}_{z}(z,\lambda_{\scriptscriptstyle 0})=(\omega\mu_{\scriptscriptstyle 0})^{-1}\textbf{$\breve{{K}}$}\breve{{e}}_{y}(z,\lambda_{\scriptscriptstyle 0})\end{array}\right\}.$$

$$\breve{E}(z,{\lambda }_0)={\mathcal{T}}_{N_{\mathrm{t}}}[{\breve{\epsilon }}_{\mathrm{rel}}(z,{\lambda }_0)]$$

$$\breve{B}(z,{\lambda }_0)={\mathcal{T}}_{N_{\mathrm{t}}}[{\breve{\beta }}_{\mathrm{rel}}(z,{\lambda }_0)].$$

$$G(x,z)=\frac{1}{\hslash c_0}{\int }_{{\lambda }_{0_{min}}}^{{\lambda }_{0_{max}}}\text{Im}\left\{{\epsilon }_{\mathrm{rel}}(x,z,{\lambda }_0)\right\}{\left|\frac{\underline{E}(x,z,{\lambda }_0)}{E_o}\right|}^{2}S({\lambda }_0)d{\lambda }_0,$$

$$J_{\mathrm{SC}}^{\mathrm{Opt}}=q_{\mathrm{e}}\frac{1}{L_{\mathrm{x}}}{\int }_0^{L_{\mathrm{z}}}{\int }_{-L_{\mathrm{x}}/2}^{L_{\mathrm{x}}/2}G(x,z)dxdz,$$

$$e_{E_s}(N_{\mathrm{t}})=\frac{\Vert E_{y,RCWA}(N_{\mathrm{t}})-E_{y,FE}{\Vert }_2}{\Vert E_{y,FEM}{\Vert }_2}$$

$$\Vert f{\Vert }_2^2={\int }_{{\Omega }_{\mathrm{ph}}}|f(x,y){|}^2\mathrm{d}z\mathrm{d}x$$

$$e_{E_p}(N_{\mathrm{t}})=\frac{\sqrt{\Vert E_{x,RCWA}(N_{\mathrm{t}})-E_{x,FEM}{\Vert }_2^2+\Vert E_{z,RCWA}(N_{\mathrm{t}})-E_{z,FEM}{\Vert }_2^2}}{\sqrt{\Vert E_{x,FEM}{\Vert }_2^2+\Vert E_{z,FEM}{\Vert }_2^2}},$$

$$\Vert f{\Vert }_2^2={\int }_{x=-L_{\mathrm{x}}/2}^{L_{\mathrm{x}}/2}{\int }_{z=0}^{L_{\mathrm{z}}}|f(x,y){|}^2\mathrm{d}z\mathrm{d}x$$

$$G(z)=\frac{1}{L_{\mathrm{x}}}{\int }_{-L_{\mathrm{x}}/2}^{L_{\mathrm{x}}/2}G(x,z)\mathrm{d}x.$$

$$\frac{d}{dz}{J}_{\mathrm{n}}(z)=-q_{\mathrm{e}}\left[G(z)-R(n,p;z)\right],$$

$$\frac{d}{dz}{J}_{\mathrm{p}}(z)=q_{\mathrm{e}}\left[G(z)-R(n,p;z)\right],$$

$${\epsilon }_0\frac{d}{dz}\left[{\epsilon }_{dc}(z)\frac{d}{dz}\varphi (z)\right]=-q_{\mathrm{e}}\left[N_{\mathrm{f}}(z)+p(z)-n(z)\right].$$

$$J_{\mathrm{n}}(z)=q_{\mathrm{e}}{\mu }_{\mathrm{n}}(z)n(z)\frac{d}{dz}{\mathsf{E}}_{{\mathrm{F}}_{\mathrm{n}}}(z)$$

$$J_{\mathrm{p}}(z)=q_{\mathrm{e}}{\mu }_{\mathrm{p}}(z)p(z)\frac{d}{dz}{\mathsf{E}}_{{\mathrm{F}}_{\mathrm{p}}}(z),$$

$${\mathsf{E}}_{{\mathrm{F}}_{\mathrm{n}}}(z)$$

$${\mathsf{E}}_{{\mathrm{F}}_{\mathrm{p}}}(z)$$

Full text

1 Introduction

The simulation of thin-film photovoltaic solar cells requires the coupling of (i) an optical model capable of capturing the absorption of photons with (ii) an electrical model capable of simulating the transport of charge carriers throughout the solar cell [1, 2]. To optimize solar cell designs, both models needed to be tailored to be rapidly computable, accurate, and robust across the relevant parameter space. To this end, we have developed a coupled optoelectronic simulation technique for thin-film photovoltaic solar cells containing periodic structures, and we have used this simulation technique together with the differential evolution algorithm (DEA) [3, 4] to optimize device performance.

The solar cell is assumed to be infinitely extended in the $xy$-plane, periodic in the $x$-direction, translation invariant in the $y$-direction, and of finite extent in the $z$-direction. These assumptions allow the simulation to be limited to two-dimensional (2D) space. However, we emphasize that our approach will work in three-dimensional (3D) space at the cost of greater computational complexity and cost.

The first step of the coupled optoelectronic simulation involves modeling the photonic characteristics of the solar cell. In the photonic step, the rigorous coupled-wave approach (RCWA) [5, 6, 7] is used to determine the electromagnetic fields in the solar cell due to incident solar radiation. The RCWA is an efficient computational technique to solve the frequency-domain Maxwell equations and determine the electric field in a representative subsection of the periodic (2D or 3D) domain. The frequency-domain electric field is used to determine the photon absorption rate, and therefore the generation rate of electron–hole pairs, throughout the solar cell. Rapid calculation across the solar spectrum [8, 9] is possible due to RCWA being a pseudo-spectral method, i.e., it avoids spatial meshing in the periodic dimensions.

The second step of the coupled optoelectronic simulation involves modeling the electronic characteristics of the solar cell. In the electronic step, the electron-hole-pair generation rate is used as the input to a drift-diffusion electronic (DDE) model in order to calculate the current-voltage curve of the solar cell and, hence, the efficiency of the solar cell. Because of the cost of solving the DDE model, we average the generation rate across one period in the $xy$-plane and implement the DDE model for variations of the electric potential and charge-carrier transport along the $z$-axis. The resulting one-dimensional (1D) DDE model consists of six coupled differential equations. The nonlinear Shockley–Read–Hall, Auger, and radiative terms are included to model electron-hole recombination [1, 2]. A hybridizable discontinuous Galerkin (HDG) method [10, 11, 12, 13] is developed, and the Newton–Raphson method [14] is used to find a solution of the resulting nonlinear system.

A version of the HDG method has already been used to simulate organic solar cells [10], based on the work of Lehrenfeld [11]. We use a similar scheme based on the classical HDG method [12], following Fu et al. [13] who have analyzed the method for a linear convection-diffusion system. Our study points to the advantages of this method in simulating inorganic solar cells having a spatially variable electron-hole-pair generation rate in the presence of semiconductor heterojunctions.

Anderson et al. [15] recently reported a coupled optoelectronic simulation of a Schottky-barrier thin-film solar cell, with the finite-element method (FEM) [16] used for both the photonic and the electronic steps of the simulation. Later, Anderson et al. [17] used the RCWA for the photonic step and the standard FEM for the electronic step for simulating the same Schottky-barrier thin-film solar cell (which is not considered in this paper). Other researchers have also developed techniques for optoelectronic simulation [18, 19, 20].

Our emphasis in this paper is different in two ways: first, our simulations lead to optimization of thin-film photovoltaic solar cells containing periodic structures; second, we concentrate on understanding the numerical performance of each of the two steps in the coupled optoelectronic simulation. For example, we investigate the convergence rate of RCWA for two model problems. We also propose, test, and use a high-order HDG scheme using fifth-degree polynomials for the DDE model, quite possibly the first time that such a high-order scheme has been used. Since the use of DEA for optimization of solar cells has been investigated elsewhere by us [8, 17], we devote little space to it in this paper.

The photonic and electronic steps described in this paper are closely related to components of the Solcore software package described recently by Alonso-Álvarez et al. [20]. That software package includes an RCWA solver as well as a 1D DD solver based on quasi-Fermi levels. In contrast, our DD solver uses the density of holes and electrons and is based on a high-order HDG scheme.

The ultimate goal is to maximize the efficiency of the solar cell, and this goal distinguishes our work from that of Alonso-Álvarez et al. [20]. Because of the presence of local minima, we use the DEA [3, 4] for optimization. This requires the evaluation of the efficiency of solar-cell designs from all regions of the multi-dimensional space encompassing appropriate parameters, which thus necessitates the development of a robust and efficient solver. See https://www.pvlighthouse.com.au for a list of other available photovoltaic-design software. Our work here should be seen as a complement to those efforts by extending the use of coupled optoelectronic simulation to the optimal design of solar cells as well as testing new simulation techniques.

In this paper, we assume that the solar cell occupies the region

[TABLE]

in $\mathbb{R}^{3}$, with air occupying the region $\mathbb{R}^{3}\setminus{\cal X}$. As stated previously, the device is assumed to be periodic along the $x$-direction with period $L_{\rm x}$ and translation invariant along the $y$-direction. Solely for illustrating various simulation issues, we performed calculations for the device shown Fig. 1(a). The current-generating region of this device comprises a thin layer of $p$-type semiconductor, a thicker layer of $i$-type semiconductor, and a thin layer of $n$-type semiconductor [1, 2]. The plane $z=0$ is taken to be the bottom face of the $n$-type semiconductor and the plane $z=L_{\rm z}$ is taken to be the top face of the $p$-type semiconductor. Below the plane $z=0$ is a grating formed by infinitely long strips of a metal and a dielectric material, the grating period being $L_{\rm x}$. Below the grating is a metal layer that is sufficiently thick to prevent light from escaping into air below it. Above the $p$-type semiconductor is an antireflection coating. The plane $z=L_{\rm t}$ is the top face of the antireflection coating and the plane $z=-L_{\rm m}$ is the bottom face of the metal layer below the grating. The grating shape in Fig. 1(a) is an academic example. A triangular grating used to test RCWA is shown in Fig. 1(b).

The layout of the paper is as follows. In Section 2, we start by describing the optical model solved in the photonic step, summarize the 2D RCWA, and define the charge-carrier generation rate that is used as the input to the DDE model solved in the electronic step. In Section 2.2 we show two examples of RCWA convergence to illustrate the advantages and drawbacks of this approach for solar cells. Section 3 is devoted to the electronic step. The DDE model is described in Section 3.1. The boundary conditions derived from local quasi-thermal equilibrium (LQTE) are discussed in Section 3.2, while the electron-hole recombination terms are presented in Section 3.3. We discuss heterojunctions in Section 3.4, and summarize a non-dimensionalized DDE model in Section 3.5. Then we proceed to describe the HDG method in Section 4, starting with the formulation of a generalized transport system in Section 4.1, followed by the discretization of 1D space (along the $z$ axis) in Section 4.2 and implementational aspects of the convection term in Section 4.3. Next, in Section 4.4, we present the upwinding strategy needed to stabilize the DDE model, followed in Section 4.5 by the interpolation method for implementing the non-linear recombination terms. The implementation of jump conditions is presented in Section 4.6. The numerical model is completed by discussing the homotopy (or continuation) method for solving the nonlinear system in Section 4.7. A numerical test is presented in Section 5. Together, the photonic and electronic steps provide information used by the DEA optimization scheme, which is described briefly in Section 6. The paper culminates in a presentation of results for the model solar cell presented in Fig. 1(a). We end with some conclusions in Section 8.

A note about notation: Underlined quantities represent vectors, with $\hat{\underline{x}}$, $\hat{\underline{y}}$ and $\hat{\underline{z}}$ being the unit vectors along the $x$-, $y$-, and ${z}$-axes. The $\;\breve{{}}\;$ mark is used to denote quantities emerging from the implementation of the spatial Fourier transform [21] with respect to $x$. For the frequency-domain calculations carried out in the photonic step (Section 2), electromagnetic fields are taken to vary with time $t$ harmonically as $\exp(-i\omega t)$, where $\omega$ is the angular frequency and $i=\sqrt{-1}$. Air is assumed to have the same electromagnetic properties as free space or vacuum. The wavenumber in air is denoted by $k_{\scriptscriptstyle 0}=\omega/c_{\scriptscriptstyle 0}$ where $c_{\scriptscriptstyle 0}$ is the speed of light in air, and $\lambda_{\scriptscriptstyle 0}=2\pi/k_{\scriptscriptstyle 0}$ is the wavelength of light in air. The permittivity and permeability of free space are denoted, respectively, by $\varepsilon_{\scriptscriptstyle 0}=8.854188\times 10^{-12}$ F m*-1* and $\mu_{\scriptscriptstyle 0}=4\pi\times 10^{-7}$ H m*-1*, and $\eta_{\scriptscriptstyle 0}=\sqrt{\mu_{\scriptscriptstyle 0}/\varepsilon_{\scriptscriptstyle 0}}$ is the intrinsic impedance of free space.

2 Photonic Step

In this section, first we summarize the RCWA for solving the frequency-domain Maxwell equations, and then we present some numerical results to quantify the performance of RCWA, which is the workhorse for computational investigations of optical gratings [22, 23]. The adequacy of this approach for solar cells has been established by comparison to other approaches [24, 25, 26]. Oddly, the numerical convergence of this method has not been investigated before; hence, a part of this section is devoted to a numerical investigation of RCWA convergence.

2.1 RCWA Formulation

The solar cell is taken to be illuminated from the half space $z>L_{\rm t}$ by a normally incident plane wave whose electric field phasor is denoted by

[TABLE]

The parameters $a_{\rm s}$ and $a_{\rm p}$ determine the polarization state of the incident light; thus, $a_{\rm p}=1$ and $a_{\rm s}=0$ for $p$-polarized light, but $a_{\rm s}=1$ and $a_{\rm p}=0$ for $s$-polarized light. For solar-cell calculations, we set $a_{\rm p}=a_{\rm s}=1$ because direct sunlight can be assumed to be unpolarized.

Due to the periodicity of the solar cell in the $x$-direction, the optical relative permittivity $\varepsilon_{\rm rel}(x,z,\lambda_{\scriptscriptstyle 0})$ is represented everywhere by the Fourier series

[TABLE]

where ${\varepsilon}_{\rm rel}^{(\ell)}(z,\lambda_{\scriptscriptstyle 0})$ are Fourier coefficients and $\kappa^{(\ell)}=2\pi\ell/L_{\rm x}$.In a similar way, the optical relative impermittivity $\beta_{\rm rel}(x,z,\lambda_{\scriptscriptstyle 0})=1/\varepsilon_{\rm rel}(x,z,\lambda_{\scriptscriptstyle 0})$ is repesented everywhere by the Fourier series

[TABLE]

where ${\beta}_{\rm rel}^{(\ell)}(z,\lambda_{\scriptscriptstyle 0})$ are Fourier coefficients. We can also express the $x$-dependences of the electric and magnetic field phasors by their Fourier series as

[TABLE]

where ${\underline{e}}^{(\ell)}\equiv{e}_{x}^{(\ell)}\hat{\underline{x}}+{e}_{y}^{(\ell)}\hat{\underline{y}}+{e}_{z}^{(\ell)}\hat{\underline{z}}$ and ${\underline{h}}^{(\ell)}\equiv{h}_{x}^{(\ell)}\hat{\underline{x}}+{h}_{y}^{(\ell)}\hat{\underline{y}}+{h}_{z}^{(\ell)}\hat{\underline{z}}$ are Fourier coefficients and normal incidence is implicit.

For computational tractability, all four of the foregoing series are truncated to include only $\ell\in\left\{-N_{\rm t},...,N_{\rm t}\right\}$, $N_{\rm t}\geq 0$. The Fourier coefficients are collected in the $(2N_{\rm t}+1)$-column vectors

[TABLE]

where the superscript $T$ denotes the transpose. Furthermore, the $(2N_{\rm t}+1)\times(2N_{\rm t}+1)$ matrix

[TABLE]

is defined for convenience. Finally, the Toeplitz matrix of the Fourier coefficients $\{\xi^{(\ell)}\}_{\ell=-\infty}^{\infty}$ of a periodic function $\xi(x)=\xi(x\pm L_{\rm x})$ is defined as

[TABLE]

Substitution of the four Fourier series into the frequency-domain Maxwell curl equations yields two matrix ordinary differential equations [7, Sec. 2.3.4]. The equation

[TABLE]

must be solved if $a_{\rm p}\neq 0$, and the equation

[TABLE]

if $a_{\rm s}\neq 0$. In these equations, $\breve{{0}}$ is the $(2N_{\rm t}+1)\times(2N_{\rm t}+1)$ null matrix, and $\breve{{I}}$ is the $(2N_{\rm t}+1)\times(2N_{\rm t}+1)$ identity matrix. As a corollary to the mutual independence of Eqs. (25) and (33), $p$-polarized light does not interact with $s$-polarized light and vice versa.

We solved both Eqs. (25) and (33) because direct sunlight can be assumed to be unpolarized. Thereafter, the Fourier coefficients of the $z$-components of the electric and magnetic field phasors were obtained as

[TABLE]

The matrices $\breve{{M}}$ and $\breve{{N}}$ in Eqs. (25) and (33) encode the Fourier series provided in Eqs. (3) and (4) for $\varepsilon_{\rm rel}$ and $\beta_{\rm rel}$. There are four choices possible for the combination of $\breve{{M}}$ and $\breve{{N}}$, the best choice depending on the polarization state of the incident plane wave [27]. Following Lalanne and Morris [5], we first define the Toeplitz matrices

[TABLE]

and

[TABLE]

Then the four possible choices are identified as follows:

Choice 1:

The most direct choice is for $\breve{{M}}$ to encode $\varepsilon_{\rm rel}$ and $\breve{{N}}$ to encode $\varepsilon_{\rm rel}^{-1}$ so that $\breve{{M}}=\breve{{E}}$ and $\breve{{N}}=\breve{{B}}$.

Choice 2:

We can write $\varepsilon_{\rm rel}=1/\beta_{\rm rel}$ and arrive at the choice $\breve{{M}}=\breve{{B}}^{-1}$ and $\breve{{N}}=\breve{{B}}$.

Choice 3:

In contrast to Choice 2, we can write $\beta_{\rm rel}=1/\varepsilon_{\rm rel}$ to choose $\breve{{M}}=\breve{{E}}$ and $\breve{{N}}=\breve{{E}}^{-1}$.

Choice 4:

The final choice is $\breve{{M}}=\breve{{B}}^{-1}$ and $\breve{{N}}=\breve{{E}}^{-1}$, in direct contrast to Choice 1.

Choices 1 and 3 are identical on the one hand, and Choice 2 is the same as Choice 4 on the other, for $s$-polarized illumination. All four choices are distinct for $p$-polarized illumination.

In order to solve Eqs. (25) and (33), it is usual to discretize in the $z$-direction by subdividing the solar cell into thin slices of infinite extent in the $xy$-plane [7]. In each slice, $\varepsilon_{\rm rel}$ is uniform in the $z$-direction, but it is either uniform or piecewise uniform in the $x$-direction. Depending on the profile of the grating, this may require piecewise-uniform approximations of $\varepsilon_{\rm rel}$ and $\beta_{\rm rel}$. Equations (25) and (33) can then be solved exactly in each slice using a stabilized stepping algorithm  [7, 23] akin to Gaussian elimination without pivoting [14], so that $\breve{{e}}_{x}$, $\breve{{e}}_{y}$ , $\breve{{h}}_{x}$, and $\breve{{h}}_{y}$ are determined everywhere in the solar cell. Equations (36) then let us determine $\breve{{e}}_{z}$ and $\breve{{h}}_{z}$ throughout the solar cell.

The convergence rates for the Fourier series of $\varepsilon_{\rm rel}$ and $\beta_{\rm rel}$ are limited by the discontinuities in the relative permittivity in the slices in the grating region. Additionally, approximating the interfaces of general grating profiles using stairstepping can also give low-order convergence with respect to the slice thickness. However, if the grating is in fact a superposition of a finite number of square waves, the only approximation is due to the truncation of the Fourier series.

2.1.1 Electron-Hole-Pair Generation Rate

With the assumption that the absorption of every photon in any semiconductor layer excites an electron-hole pair, the electron-hole-pair generation rate in the semiconductor region $\left\{|{x}|<L_{\rm x}/2,0<z<L_{\rm z}\right\}$ can be calculated as

[TABLE]

where $\hbar$ is the reduced Planck constant and $S(\lambda_{\scriptscriptstyle 0})$ is the incident power spectrum. We chose $S(\lambda_{\scriptscriptstyle 0})$ to be the standard AM1.5G spectrum [28]. As $S(\lambda_{\scriptscriptstyle 0})$ falls off rapidly at shorter free-space wavelengths, ${\lambda_{\scriptscriptstyle 0_{min}}}=350$ nm was fixed. Once the free-space wavelength becomes too large, the energy in a photon is not capable of exciting an electron sufficiently to cross the bandgap. Consequently, the upper cut off was taken to be $\lambda_{\scriptscriptstyle 0_{max}}=1240\leavevmode\nobreak\ \text{nm V}/{\sf E}_{\rm g,min}$, where ${\sf E}_{\rm g,min}$ (in V) is the minimum value of the bandgap ${\sf E}_{\rm g}$ in the solar cell. These choices are typical values for studies on solar cells.

After assuming full quantum efficiency, i.e. each excited electron or hole is collected at the appropriate terminal, the optical short-circuit current density can be calculated as

[TABLE]

where $q_{\rm e}=1.60218\times 10^{-19}$ C is the elementary charge. Let us note that $J_{\rm SC}^{\rm Opt}$ provides only a rough benchmark for the device efficiency, although it is often used to assess the efficiency of solar cells [29, 30, 31]. However, as recombination is neglected, $J_{\rm SC}^{\rm Opt}$ is necessarily larger than the attainable short-circuit current density $J_{\rm SC}$, which is the density of current that flows when the solar cell is illuminated and no external bias is applied (i.e., when $V_{\rm ext}=0$). We have already shown that, when the bandgap is one of the simulation parameters, optimizing Eq. (40) can lead to designs with high values of $J_{\rm SC}^{\rm Opt}$ but low values of extractable electrical power [8]. A better approach is to complement the photonic step by the electronic step, as accomplished in Section 3.

2.2 RCWA Convergence

There are two sources of approximation error in the RCWA. The first is the truncation of the Fourier series to include only $\ell\in\{-N_{\rm t},...,N_{\rm t}\}$, while the second is through the discretization of space in the $z$-direction whenever there is a need to approximate the relative permittivity slice by slice. As it is not possible to solve Eqs. (25) and (33) analytically, in order to investigate the convergence of RCWA, the electric field phasor ${\underline{E}}(x,z,\lambda_{\scriptscriptstyle 0})$ was calculated using an FEM code [24, 25] and compared with the RCWA results for various values of $N_{\rm t}$, the slice thickness being fixed at a sufficiently low value.

The FEM data were obtained using an adaptive method implemented in NGSolve [32]. The adaptive algorithm uses mesh bisection and the Zienkiewicz–Zhu a-posteriori error estimator [33]. The FEM calculations were carried out using 5th-order continuous finite elements. Mesh adaptivity was terminated whenever the algorithm reached 100,000 degrees of freedom. The simulated domain was sandwiched between two perfectly matched layers (PMLs). Both of the PMLs were taken to be one-wavelength thick and with a constant PML parameter of $(1.5+2.5i)$ [34], the reflection coefficient of the PMLs being then \sim$$3\times 10^{-12}.

2.2.1 Global Convergence of Electric Field Phasor

In order to investigate the global convergence of the RCWA, a simple but adequate problem was chosen: reflection by the grating shown in Fig. 1(b). The region above the grating was taken to be occupied by air with relative permittivity $\varepsilon_{\rm rel}=1+10^{-9}i$, the minuscule imaginary part needed for the stability of the RCWA algorithm. The grating, made from a fictitious metallic material with relative permittivity $\varepsilon_{\rm rel}=-22+0.4i$, comprises triangular protrusions on a $50$-nm-thick film. Each protrusion is of height 50 nm and base 250 nm. Calculations were made with $\lambda_{\scriptscriptstyle 0}=600$ nm and $L_{\rm x}=500$ nm.

The triangular profile was chosen for the grating, because it is approximated by a stairstep when the piecewise-uniform approximation is used to solve Eqs. (25) and (33). In addition, the true solution has a singularity in the field at the tip of the triangle which tests the algorithm’s ability to handle such behavior.

As the only non-zero Cartesian component of $\underline{E}(x,z,\lambda_{\scriptscriptstyle 0})$ for $s$-polarized illumination is the $y$-component, a relative global error was defined as

[TABLE]

as a function of $N_{\rm t}$, with $E_{y,RCWA}$ denoting the value of $E_{y}$ yielded by RCWA, and $E_{y,FEM}$ the value yielded by FEM. In addition, for any function $f(x,y)$,

[TABLE]

and where ${\Omega_{\rm ph}}$ is the entire simulation domain in the photonic step except the two PMLs. For $p$-polarized illumination, the relative global error was analogously defined as

[TABLE]

because only the $y$-component of $\underline{E}(x,z,\lambda_{\scriptscriptstyle 0})$ is null valued.

Four different choices of $\left\{\breve{{M}},\breve{{N}}\right\}$, as described after Eq. (38), are possible when implementing the RCWA. Accordingly, four different studies were conducted to investigate the relative global error in the electric field phasor for $s$- and $p$-polarized illuminations.

Plots of $e_{E_{s}}$ and $e_{E_{p}}$ vs. $N_{\rm t}$ are presented in Fig. 2 for each of the four choices of $\left\{\breve{{M}},\breve{{N}}\right\}$. For $s$-polarized illumination, the relative global error reduces as $\sim{N_{\rm t}^{-2.1}}$ until $N_{\rm t}\approx 8$ and saturates thereafter in Fig. 2(a), when $\breve{{M}}=\breve{{E}}$ (i.e., for Choices 1 and 3). The error saturation may be due to the limited accuracy of the FEM solution and/or due to stairstepping error. In contrast, $e_{E_{s}}$ reduces as $\sim{N_{\rm t}^{-1}}$ when $\breve{{M}}=\breve{{B}}^{-1}$ (i.e., for Choices 2 and 4). For $p$-polarized illumination, the relative global error decays slower than the reciprocal of $N_{\rm t}$ in Fig. 2(b), when $\breve{{M}}=\breve{{E}}$. The decay is somewhat faster when $\breve{{M}}=\breve{{B}}^{-1}$.

These results suggest that for calculating the electric field phasor due to $s$-polarized illumination, it is advantageous to choose $\breve{{M}}=\breve{{E}}$ because that should give second-order convergence to the true solution. For $p$-polarized illumination, it is preferable to choose $\breve{{M}}=\breve{{B}}^{-1}$ and $\breve{{N}}=\breve{{E}}^{-1}$, as that should give almost first-order convergence to the true solution even with stairstepping. The choice of Fourier representations for $\left\{\breve{{M}},\breve{{N}}\right\}$ thus has a profound effect on the accuracy of the RCWA [5, 27].

2.2.2 Convergence of Electric Field Phasor in the Semiconductor Region

Correct calculation of the electric field phasor in the semiconductor region is paramount for photovoltaic solar cells, as is clear from Eq. (39). This motivated the second convergence study using the device in Fig. 1(a).

The chosen device contains a grating made from a representative metal with relative permittivity $\varepsilon_{\rm rel}=-22+0.4i$ and has period $L_{\rm x}=500$ nm. The grating comprises rectangular protrusions on a $100$-nm-thick film. Each protrusion is of height 50 nm and base 250 nm. The grooves of the grating are entirely filled with a dielectric material with relative permittivity $\varepsilon_{\rm rel}=3.33+0.016i$. Above this structure lies a $240$-nm-thick semiconductor region with relative permittivity $\varepsilon_{\rm rel}=9.5+1.25i$, which is a representative value for a semiconductor at $\lambda_{\scriptscriptstyle 0}=680$ nm. Above this region is a 75-nm-thick layer of a dielectric material with relative permittivity $\varepsilon_{\rm rel}=3.33+0.016i$, which acts as an anti-reflection coating. The top layer is again air-like, with relative permittivity $\varepsilon_{\rm rel}=1+10^{-9}i$.

For this study, we used the definition

[TABLE]

in Eqs. (41) and (43). Four different choices of $\left\{\breve{{M}},\breve{{N}}\right\}$, as described after Eq. (38), were used in four different studies, whose results are presented in Fig. 3.

Convergence rates for the second study are much slower than in the first study. The best observed rates are $O(N_{\rm t}^{-1.6})$ for $s$-polarized illumination and $O(N_{\rm t}^{-0.35})$ for $p$-polarized illumination. The slower convergence is due to the strong singularities in the solution near the corners of the grating. Those singularities lead to slower convergence of the Fourier series. It should be noted that these same singularities are also observed in the adaptive FEM solution.

Figure 3(a) contains a split in $e_{E_{s}}$ between successive odd and even values of $N_{\rm t}$, when $\breve{{M}}=\breve{{E}}$ and $\breve{{N}}=\breve{{E}}^{-1}$ (Choice 3). This split is also present, although less extreme, in an earlier work [5]. It is also present for $p$-polarized illumination but dramatically reduced, especially when $\breve{{M}}=\breve{{E}}$ and $\breve{{N}}=\breve{{B}}$ (Choice 1).

Lalanne [6] has already discussed in detail the complex issues behind RCWA convergence. Our contribution is to suggest that stairstepping is not a major cause of error in RCWA, but it is necessary to examine convergence for different grating morphologies. Also, RCWA is reliably convergent for $s$-polarized illumination of solar cells, and provides slower convergence for $p$-polarized illumination. We agree with previous studies [24, 25, 26] that RCWA can provide the accuracy needed for photonics simulations, but our studies also show that the convergence rates can be quite low. With an appropriate choice of parameters we can achieve almost 10% error with $N_{\rm t}=10$, which became our choice for the optimization study reported in Section 6.

3 Electronic Step

Taking into account the size of the features in the solar cell, and under standard assumptions on the electron transport, such as the device being in quasi-thermal equilibrium (QTE), the transport of electrons and holes in semiconductors can be modeled using the drift-diffusion equations [1, 2, 35]. To speed up the calculations, we neglected variations with respect to $x$ and therefore used the 1D electron-hole-pair generation rate

[TABLE]

3.1 Drift-Diffusion Electronic Model

The DDE model comprises the following three differential equations holding in the semiconductor region \Omega_{\rm el}=\left\{z\;\big{|}\;0<z<L_{\rm z}\right\}:

[TABLE]

In these equations, the gradients of the electron current density $J_{n}(z)$, hole current density $J_{\rm p}(z)$, and electric potential $\phi(z)$ are related to the electron density $n(z)$, hole density $p(z)$, fixed charge (charged traps and doping) density $N_{\rm f}(z)$, electron-hole-pair generation rate $G(z)$, and electron-hole-pair recombination rate $R(n,p;z)$. The zero-frequency (dc) relative permittivity is denoted by $\varepsilon_{dc}(z)$.

The two current densities are defined as

[TABLE]

and

[TABLE]

where $\mu_{\rm n}(z)$ and $\mu_{\rm p}(z)$ are the electron mobility and hole mobility, respectively;

[TABLE]

is the electron quasi-Fermi level; and

[TABLE]

is the hole quasi-Fermi level. In the expressions for the quasi-Fermi levels, $V_{\rm th}=k_{\rm B}T/q_{\rm e}$ is the thermal voltage of the electrons and holes with $T$ as the temperature and $k_{\rm B}=1.380649\times 10^{-23}$ J K*-1* as the Boltzmann constant, $N_{\rm c}(z)$ is the conduction-band density of states, $N_{\rm v}(z)$ is the valence-band density of states, ${\sf E}_{\rm c}(z)$ is the conduction-band energy, and ${\sf E}_{\rm v}(z)$ is the valence-band energy. The bandgap is defined as

[TABLE]

and the intrinsic energy as

[TABLE]

In the Boltzmann approximation [1] the quasi-Fermi levels are assumed to lie far from the edges of the conduction and valence bands. Accordingly,

[TABLE]

and

[TABLE]

Here, the intrinsic carrier density is given by

[TABLE]

and the conduction-band energy is given by

[TABLE]

where $\chi(z)$ is the electron affinity of the semiconductor, and the (arbitrary) reference energy level ${\sf E}_{0}$ is often chosen as [1]

[TABLE]

Equations (49) and (50) may now be simplified to

[TABLE]

and

[TABLE]

respectively. Here,

[TABLE]

and

[TABLE]

are the built-in potentials for the electrons and holes (due to variations in the material properties), respectively. Both built-in potentials as well as the electron affinity may be discontinuous with respect to position at heterojunctions in the semiconductor region. The baseline number density $N_{0}$ is arbitrary because potentials are defined up to a constant.

Equations (46), (47), (48), (61), and (62), have to be solved concurrently for $z\in(0,L_{\rm z})$, in conjunction with a set of boundary conditions for $n(z)$, $p(z)$, and $\phi(z)$. We opted for the Dirichlet choice

[TABLE]

because it models an ideal ohmic contact [1, 2]. Herein, the functions $n_{0}(z)$, $p_{0}(z)$ and $\phi_{0}(z)$ are the electron density, hole density and potential at LQTE, as discussed in Section 3.2, while $V_{\rm ext}$ is the externally applied voltage across the terminals at the boundaries of the semiconductor region. Solution of the system of five equations enables the calculation of the current density

[TABLE]

flowing uniformly through the solar cell.

The current density $J$ depends on the choice of $V_{\rm ext}$, i.e., $J\equiv J(V_{\rm ext})$. Repeating calculations for various values of $V_{\rm ext}$ produces the $JV$-curve, with the maximum value of the power density

[TABLE]

indicating the maximum power density

[TABLE]

obtainable from the solar cell. This in turn gives the efficiency of the solar cell as

[TABLE]

where $P_{\rm in}=1000$ W m*-2* is the incident solar power density. It is the efficiency of the solar cell that we wish to optimize.

3.2 Local Quasi-Thermal Equilibrium

If a region containing a homogeneous semiconductor is charge-free and isolated from any external influences (e.g., $G\equiv 0$ and $V_{\rm ext}=0$), the distributions of electrons and holes with respect to energy tend to the Fermi distribution, with the Fermi-level lying at the intrinsic energy level, i.e., ${\sf E}_{\rm F_{\rm n}}={\sf E}_{\rm F_{\rm p}}={\sf E}_{F}$. Then, $J_{\rm n}=J_{\rm p}$ are identically zero and the semiconductor is said to be in LQTE. This concept is used to model the ideal ohmic boundary conditions mentioned in Section 3.1.

With $n_{0}(z)$ and $p_{0}(z)$ denoting the electron and hole densities, respectively, at LQTE, Eqs. (56)–(58) can be combined to give the mass-action equation

[TABLE]

furthermore,

[TABLE]

as the semiconductor is charge-free. Combining Eqs. (72) and (73), we get

[TABLE]

At LQTE, the electron and hole quasi-Fermi levels coincide and are uniform. Thus,

[TABLE]

must be equal to each other for all $\upsilon\in(0,L_{\rm z})$ and $\zeta\in(0,L_{\rm z})$. Substitution of Eqs. (54) and (59) into ${\sf E}_{\rm F_{\rm n}}(\upsilon)={\sf E}_{\rm F_{\rm p}}(\zeta)$ yields

[TABLE]

Equation (78) simplifies to

[TABLE]

with the choice $\upsilon=0$ and $\zeta=z$, and to

[TABLE]

with the choice $\upsilon=z$ and $\zeta=0$, where we have chosen $\phi_{0}(0)=0$ without loss of generality. Subtracting Eq. (79) from Eq. (80), we get

[TABLE]

Equations (74), (75), and (81) give the initial functions needed to specify the boundary data identified in Eqs. (65)–(67). However, in view of the LQTE assumption (73), these initial functions do not satisfy the system of ODEs (46)–(48).

3.3 Recombination

Recombination of an electron and a hole can take place through several different mechanisms [1, 2]. The electronic step in our program accommodates three different contributions to $R(n,p;z)$.

3.3.1 Radiative Recombination

Radiative recombination occurs when an electrons and a hole recombine across the full bandgap, releasing the energy as a photon with energy equal to the bandgap. At QTE, radiative recombination is identically balanced by electrons being thermally excited across the bandgap. The (net) radiative recombination rate is given by

[TABLE]

where $\alpha_{rad}(z)$ depends on the semiconductor material. It should be noted that $R_{rad}(z)\equiv 0\,\forall z\in\left[0,L_{\rm z}\right]$ at LQTE.

3.3.2 Shockley–Read–Hall (SRH) Recombination

The SRH recombination contribution is due to electrons and holes recombining via an intermediate gap state. It is modeled by

[TABLE]

where $n_{1}(z)$ and $p_{1}(z)$ are, respectively, the electron and hole densities at the trap energy level. If this level is the intrinsic energy level ${\sf E}_{\rm i}(z)$, then $n_{1}(z)=p_{1}(z)=\bar{n}(z)$ from Eqs. (56) and (57). The functions $\tau_{\rm n}(z)$ and $\tau_{\rm p}(z)$ are material dependent.

3.3.3 Auger Recombination

The Auger recombination contribution arises from a three-particle recombination pathway, occurring when an electron and a hole recombine across the bandgap, with the released energy transferred to a third charge carrier which is excited away from both band edges. This recombination rate is given by

[TABLE]

where the functions $C_{\rm n}(z)$ and $C_{\rm p}(z)$ are material dependent.

All three contributions add, so that

[TABLE]

It can be seen that the total recombination $R=0$ at LQTE irrespective of the choice of the material parameters. For all results reported in this paper, we set $\alpha_{rad}(z)\equiv 0$ and $C_{\rm n}(z)=C_{\rm p}(z)\equiv 0\,\forall z\in\Omega_{\rm el}$ and only the SRH recombination was activated, in order to present illustrative results without the expenditure of significant computational time.

3.4 Heterojunctions: Continuous Quasi-Fermi Levels

When there is a jump in either $\chi(z)$ or ${\sf E}_{\rm g}(z)$ or both, and thus in $\phi_{\rm n}(z)$ or $\phi_{\rm p}(z)$, then $n(z)$ and $p(z)$ may also have jumps. A heterojunction is formed at the discontinuity [36]. A commonly used way to model this discontinuity requires the assumption of continuous quasi-Fermi levels, whereas another commonly used way requires consideration of thermionic emission at the discontinuity [37]. We chose to implement the continuous quasi-Fermi level (CQFL) model, although the thermionic-emission model is also compatible with the HDG method.

The CQFL model uses the limit of a continuum model to quantify the jump at a heterojunction. To understand the resulting jump conditions, let us examine a small region $0<z\leq\delta$ in which the electron affinity $\chi(z)$ changes linearly by $\Delta\chi$, while $N_{\rm c}(z)$ is uniform. Then, differentiating Eq. (63) with respect to $z$ yields

[TABLE]

in this region. Furthermore, we assume that $G\delta\approx R\delta\approx 0$ in this region so that Eq. (46) simplifies to

[TABLE]

hence, $J_{\rm n}(z)=J_{\rm n}^{\circ}$. Finally, we also assume that both $\mu_{\rm n}(z)\approx\mu_{\rm n}^{\circ}$ and $\phi(z)\approx\phi^{\circ}$ are uniform in this region, so that Eq. (61) simplifies to

[TABLE]

after using Eq. (86).

The solution of Eq. (88) is

[TABLE]

where $C$ is the constant of integration. As $n(z)=n_{\rm L}$ at $z=0$, we get

[TABLE]

from Eq. (89). As $n(z)=n_{\rm R}$ at $z=\delta$, Eq. (90) yields

[TABLE]

In the limit of $\delta\to 0$, we get the jump condition

[TABLE]

from Eq. (91). The same analysis can be performed for the hole density $p(z)$, which results in the jump condition

[TABLE]

where $\Delta{\sf E}_{\rm g}$ is the change in the bandgap over the region with thickness $\delta$. Equations (92) and (93) need to be enforced across any junction over which $\chi$ or ${\sf E}_{\rm g}$ jump. This is accomplished using the HDG method described in Section 4.

3.5 Dimensionless Formulation

The equations of the DDE model are scaled in order to ensure that the internal variables are dimensionless [35], the scaling parameters being provided in Table 1. Equations (46) and (47) are thus non-dimensionalized to111At the risk of confusion, we have not used new notation for the dimensionless variables in order to keep the text clean.

[TABLE]

respectively, and Eqs. (61) and (62) to

[TABLE]

respectively. Electing to keep $\varepsilon_{dc}(z)\equiv\varepsilon_{dc}^{0}\,\forall z\in\left[0,L_{\rm z}\right]$ uniform in the semiconductor region, we obtain the non-dimensionalized form of Eq. (48) as

[TABLE]

where the Poisson constant

[TABLE]

is dimensionless. We define the dimensionless d.c. electric field

[TABLE]

and transform Eq. (98) to

[TABLE]

which has the same form as Eqs. (94) and (95). These three ODEs are discretized using the HDG method, as described in Section 4.

4 Hybridizable Discontinuous Galerkin Formulation

The hybridizable discontinuous Galerkin (HDG) method [13, 34] possesses several features which are advantageous for the implementation of the DDE model for solar cells [10]. The relaxation of the requirement that the solution is continuous at the boundaries of elements permits solutions where the variables have strong gradients and second derivatives. It also naturally allows discontinuities in the solution due to discontinuous material parameters, such as those which occur at heterojunctions discussed in Section 3.4.

4.1 Generalized Transport System

To simplify the description of the numerical scheme used to solve Eqs. (94)–(97), (100), and (101), it is useful to partition them into three sets, each with the form

[TABLE]

along the associated Dirichlet boundary conditions

[TABLE]

Here, ${\sf Z}$ is either given by the solution of another set or vanishes, ${\sf v}_{\sf m}$ is a given function of position, and ${\sf f}(n,p)$ is also given. To get each of the three sets of equations in the dimensionless DDE model provided in Section 3.5, we choose the parameters and functions as in Table 2.

Henceforth, we refer to Eqs. (102)–(104) as the generalized transport system. We now proceed to design a numerical scheme to discretize this system of equations.

4.2 Discretization of Electronic Domain

The electronic domain $\Omega_{\rm el}$ is covered by a mesh $\mathcal{S}=\{\mathcal{N},\mathcal{T}\}$ of $N_{\rm z}+1$ nodes $\mathcal{N}$ and $N_{\rm z}$ elements $\mathcal{T}$. The element $s_{\gamma}=({z}_{\gamma-1},{z}_{\gamma})$ lies between the nodes ${z}_{\gamma-1}$ and ${z}_{\gamma}$,

[TABLE]

is the set of all nodes, and

[TABLE]

is the set of all elements. The generalized transport system is discretized using the space

[TABLE]

where $\mathbb{P}_{P_{\rm deg}}$ is the set of polynomials of degree $P_{\rm deg}$ in $z$ (and so that $\mathbb{V}_{h}$ is a space of discontinuous piecewise polynomials),

[TABLE]

and $\mathbb{R}^{+}$ is the set of positive real numbers. Although a different value of $P_{\rm deg}$ may be chosen for each element $s_{\gamma}$, for simplicity of notation we choose the same value of $P_{\rm deg}$ for all elements.

We seek the numerical solutions $({\sf J}_{h},\hat{{\sf J}}_{h},{\sf m}_{h},\hat{{\sf m}}_{h})\in\mathbb{V}_{h}\times\mathbb{W}_{h}\times\mathbb{V}_{h}\times\mathbb{W}_{h}$ where $\hat{{\sf J}}_{h}$ and $\hat{{\sf m}}_{h}$ are the vectors of numerical approximations to ${\sf J}$ and $m$ at the nodes in $\cal N$. Multiplying Eq. (102) by a test function $\psi\in\mathbb{V}_{h}$ and Eq. (103) by a test function $\eta\in\mathbb{V}_{h}$, and then integrating over ${\Omega_{\rm el}}$, we obtain

[TABLE]

where the notation

[TABLE]

For functions $\eta\in\mathbb{V}_{h}$, $\psi\in\mathbb{V}_{h}$, and $\hat{\eta}\in\mathbb{W}_{h}$, we define

[TABLE]

via integration by parts, with

[TABLE]

and

[TABLE]

Then, Eq. (109) can be written as

[TABLE]

To complete the discretization, we need to enforce approximate continuity of the discrete flux. Therefore, at each node we define the discrete flux as

[TABLE]

where $\zeta({z}_{\gamma})=\lim\limits_{z\to{z}^{\pm}_{\gamma}}\zeta(z)$ and $\tau_{m}({z}_{\gamma})\geq 0$ is a hybridization (or penalty) function to be chosen. Continuity of the flux requires

[TABLE]

at each interior node. Finally, the boundary conditions are enforced by requiring

[TABLE]

at the ends of $\Omega$.

By defining the jump operator as

[TABLE]

the problem becomes that of finding a discrete solution ${({\sf J}_{h},{\sf m}_{h},\hat{{\sf m}})}\in\mathbb{V}_{h}\times\mathbb{V}_{h}\times\mathbb{W}_{h}$ such that

[TABLE]

for all $(\psi,\eta)\in\mathbb{V}_{h}\times\mathbb{V}_{h}$, and with $\hat{{\sf m}}={\sf m}_{0}\mbox{ for }{z\in\left\{0,L_{\rm z}\right\}}$. Obviously, there are $2N_{\rm z}(P_{\rm deg}+1)$ degrees of freedom for the functions ${\sf J}$ and ${\sf m}$, and $N_{\rm z}-1$ degrees for $\hat{{\sf m}}$. In turn, the test functions $\psi$ and $\eta$ give $2N_{\rm z}(P_{\rm deg}+1)$ equations from Eqs. (121) and (122, while Eqs. (123) give a further $N_{\rm z}-1$ equations. For a linear problem of this type, existence and uniqueness of the solution have been proven and the error estimates have been derived, for example, by Fu et al. [13]; see also Cockburn et al. [12]. The nonlinear problem defined here remains to be analyzed.

To help with conditioning and also to provide a useful choice of discretization points for the nonlinear problem, we choose the interpolation points for $\mathbb{V}_{h}$ on each element to be the $P_{\rm deg}+1$ Gauss–Lobatto points [38]. We now replace functions ${\sf J}_{h}$, ${\sf m}_{h}$, ${\sf Z}_{h}$, $n_{h}$, and $p_{h}$ by their finite-element approximations with the form

[TABLE]

where $\psi_{\ell}\in\mathbb{V}_{h}$ are piecewise polynomial finite-element basis functions that are each compactly supported on one of the elements $s_{\gamma}$, and $\zeta_{\ell}\in\mathbb{R}$ are constants to be determined. Upon substitution of these into Eqs. (121)–(123), and choosing the test functions $\eta$ and $\psi$ to also be the finite-element basis functions $\psi_{j}\in\mathbb{V}_{h}$ to be piecewise polynomials that are each compactly supported on element $s_{\gamma}$, we obtain the HDG equations

[TABLE]

for $1\leq j\leq N_{\rm z}(P_{\rm deg}+1)$, $1\leq k\leq N_{\rm z}(P_{\rm deg}+1)$, and $1\leq\ell\leq N_{\rm z}(P_{\rm deg}+1)$. In Eqs. (125)–(128) and henceforth, we use the Einstein summation convention (i.e., summation is implied over any duplicated index), column vectors denoted as ${\underline{v}}$ comprise elements $v_{k}$, and matrices denoted as $\underline{\underline{M}}$ comprise elements $M_{k,\ell}$.

The full set of drift-diffusion equations (94)–(101) along with the boundary conditions (65)–(67) is now discretized by three copies of the above equations with function and parameter choices from Table 2. There are two usual ways to solve the resulting nonlinear system: Gummel iteration [35] and Newton’s method [39]. We found that straightforward Gummel iteration did not converge reliably for the parameter values in our simulations. So we elected to use Newton’s method with a scheme of the homotopy type [40] to help provide a good initial guess.

4.3 Nonlinear Convection

A practical difficulty with implementing the scheme arises when dealing with the nonlinear term

[TABLE]

The integral

[TABLE]

on the right side of Eq. (129) is an element of a $[N_{\rm z}(P_{\rm deg}+1)]^{3}$ tensor denoted by $\bm{{{\bar{\Psi}}}}^{0,0,1}$, where the $0,0,1$ superscript specifies that $\psi_{k}$ is differentiated once. The non-linear term can then be seen as the product

[TABLE]

where the indices $\ell$, $j$, and $k$ select the elements of the tensor. Although manipulation of tensors with rank greater than $2$ is not natively supported in Matlab, the tensor $\bm{{{\bar{\Psi}}}}^{0,0,1}$ is very sparse as it is block diagonal; hence, the calculation speed can be increased by storing it as a non-square matrix.

As each of the finite-element basis functions $\psi_{j}$ is compactly supported on just one interval, a local tensor

[TABLE]

with dimensions $[(P_{\rm deg}+1)]^{3}$ can be formed. This local tensor, which is not necessarily sparse, can be rewritten as a $[P_{\rm deg}+1]\times[(P_{\rm deg}+1)^{2}]$ rectangular matrix $\bm{{{\bar{\Psi}}}}^{0,0,1}_{\gamma}$. Finally, a global matrix $\bm{{{\bar{\Psi}}}}^{0,0,1}$ with dimension $N_{\rm z}(P_{\rm deg}+1)\times(N_{\rm z}(P_{\rm deg}+1))^{2}$ to represent $\bm{{{\bar{\Psi}}}}^{0,0,1}$ is formed by repeating the local matrix as a block diagonal.

In order to calculate the term $\underline{\psi}_{\sf m}$, the vector of coefficients $\underline{{\sf m}}_{\ell}$ is reshaped into a $\left[N_{\rm z}(P_{\rm deg}+1)\right]^{2}\times N_{\rm z}(P_{\rm deg}+1)$ matrix ${\underline{\underline{M}}}$: the block diagonal is formed from $P_{\rm deg}+1$ blocks of ${\sf m}_{\ell}$, with each block repeated $N_{\rm z}$ times. The required solution is then given by

[TABLE]

While this is more convoluted than simply performing the tensor multiplication, this matrix multiplication method is comparatively very fast because $\bm{{{\bar{\Psi}}}}^{0,0,1}$ is independent of the solution (i.e., constant) and so can be precomputed.

4.4 Upwinding

One of the reasons for choosing the HDG method is to include upwinding in a natural way through the choice of hybridization functions $\tau_{m}$ [13]. Upwinding was included in our simulations as follows.

First, the dimensionless effective average electron speed $u_{n}(z)=-J_{\rm n}(z)/n(z)$ and the effective average hole speed $u_{p}(z)=J_{\rm p}(z)/p(z)$ were formulated. Next, Eqs. (96) and (97) were used to obtain

[TABLE]

respectively. The terms containing $\ln[...]$ are identifiable as the diffusion terms, with the remaining terms modeling drift in the effective electric field. In drift-dominated regions of the solar cell, which often constitute most of $\Omega_{\rm el}$, we therefore get

[TABLE]

where

[TABLE]

is the dimensionless effective electric field acting on electrons and

[TABLE]

is the dimensionless effective electric field acting on holes. These two fields arise from material nonhomogeneity.

Information in this drift-dominated system travels in the directions of the electron and hole velocities. At each node, we therefore want to take values of $n(z)$ and $p(z)$ from the inflow side and pass them to the outflow side. As $\mu_{\rm n}(z)>0$ and $\mu_{\rm p}(z)>0$, the respective directions entirely depend on the signs of $E_{\rm dc}(z)+E_{n}(z)$ and $E_{\rm dc}(z)+E_{p}(z)$. While the hybridization function $\tau_{m}(z)$ is defined as a function of position, only the values at the ends of the elements contribute to the model. Consequently, we create $\bar{\tau}(z)$ as a piecewise linear function, defined by the limiting values at each side of each node. Two possible limiting values $\bar{\tau}_{1}$ and $\bar{\tau}_{2}$ are chosen such that $\bar{\tau}_{1}\ll\bar{\tau}_{2}$. Then, to ensure that the correct information flow is achieved,

• •

$\tau_{m}|_{\gamma,R}=\bar{\tau}_{1}$ and $\tau_{m}|_{\gamma,L}=\bar{\tau}_{2}$ are used if $E_{\rm dc}(z)+E_{n}(z)>0$, or

• •

$\tau_{m}|_{\gamma,L}=\bar{\tau}_{1}$ and $\tau_{m}|_{\gamma,R}=\bar{\tau}_{2}$ are used $E_{\rm dc}(z)+E_{n}(z)<0$.

Our choices of $\bar{\tau}_{1}$ and $\bar{\tau}_{2}$ are given in Table 5.

4.5 Recombination

The recombination term $R(n,p;z)$ is a nonlinear function of $n(z)$ and $p(z)$. Consequently, care needs to be taken on how to incorporate this term into the HDG method. For example, if radiative recombination given by Eq. (82) is incorporated in the HDG method, we get

[TABLE]

The first term on the right side is nonlinear in the basis function, giving a third-rank tensor with indices $j$, $k$ and $\ell$. Note that $\alpha_{rad}$ and $\bar{n}$ are also projected onto $\mathbb{V}_{h}$ but, because both are material properties (and are therefore independent of the solution), they do not increase the rank of the tensor produced. The term can thus be implemented in a similar manner to the nonlinear drift term discussed in Section 4.3.

The SRH recombination is given by Eq. (83). If we test this term against the standard polynomial basis $\psi_{j}$ and integrate over $\Omega_{\rm el}$, we cannot write the result as a tensor, because the SRH term is not a polynomial in the basis functions. Consequently, we integrate using quadrature. For this, we chose Gauss–Lobatto quadrature 222The interpolation points for the finite-element basis functions are the Gauss–Lobatto points, as discussed in Sec. 4.2. to maximize the order of the polynomials that can be exactly integrated, while maintaining use of the values at the nodes. In particular, for a general function $f$, the use of Gauss–Lobatto quadrature to integrate over the element $s_{\gamma}$ gives

[TABLE]

where $\{w_{j}\}_{j=1}^{P_{\rm deg}+1}$ are the Gauss–Lobatto quadrature weights for the element $s_{\gamma}$ and $I_{\rm deg}$ is the degree of integration. A computationally quick way is to use $I_{\rm deg}=P_{\rm deg}+1$ quadrature points per element [41, 42] as this does not require interpolation of the solution.

The Auger recombination is given by Eq. (84). As with the SRH term, we cannot write this as a tensor. Either mass-lumping  [41] or quadrature [38] must then be used to calculate the integrals formed when the Auger term is incorporated into the HDG method.

4.6 Heterojunctions

To implement the jump conditions (92) or (93), we ensure that a node falls at each discontinuity in the material parameters. Every jump in $n(z)$ can be taken into account by redefining the jump operator as

[TABLE]

where $\Delta\chi_{\gamma}$ is the jump in electron affinity $\chi$ across node $\gamma$ and $\zeta$ is a place-holder function. This is equivalent to defining to values of $\hat{n}_{\gamma}$ at each node, separated in value by the necessary jump induced by the discontinuous electron affinity. Jumps in $p(z)$ can be handled similarly.

4.7 Homotopy

In order to aid convergence of the highly nonlinear discrete system, homotopy is employed in our simulation. The fixed charge density $N_{\rm f}$, recombination rate $R(n,p;z)$, and bandgap non-homogeneity, i.e. $E_{g}(z)-E_{g,av}$ where $E_{g,av}$ is the mean bandgap, are all multiplied by a constant $\delta_{\rm homo}$. The simulation is started with $\delta_{\rm homo}^{(0)}=\delta_{\rm min}$. Once the simulation is deemed to have converged, a larger value

[TABLE]

is chosen for $\delta_{\rm homo}$, with $\varepsilon_{\rm homo}>1$ as the homotopy damping relaxation rate. Thus, the iteration

[TABLE]

is employed, full convergence being deemed to occur in the iteration in which $\delta_{\rm homo}$ becomes equal to unity.

5 Numerical Test of HDG Method

As a model test problem to show the behavior of the HDG method (but not to present a solar-cell design), the solar cell was taken to comprise a p-i-n junction made from copper-indium-gallium-(di)selenide (CIGS), as shown in Fig. 1(a). The electrical parameters of CIGS [43] are given in Table 3. For this test problem, the generation rate was taken to be uniform, with $J_{\rm SC}^{\rm Opt}=10$ mA cm*-2*; thus RCWA results were not used. The $p$-type layer of thickness $L_{\rm p}=20$ nm was taken to be doped with acceptor atoms with concentration specified via $N_{\rm f}=-10^{17}$ cm*-3*, and the $n$-type layer of thickness $L_{\rm n}=20$ nm with donor atoms with concentration specified via $N_{\rm f}=10^{17}$ cm*-3*. The undoped $i$-type layer was taken to be $L_{\rm i}=200$ nm thick; thus, $L_{\rm z}=L_{\rm p}+L_{\rm i}+L_{\rm n}=240$ nm. All three layers were taken to be homogeneous. Neither the radiative nor the Auger recombination mechanism was activated for the test problem.

Chiefly, three parameters affect the accuracy of the HDG method: the degree of the interpolating polynomials $P_{\rm deg}$, the length of each element $d_{z}$, and the degree of quadrature integration $I_{\rm deg}$. The convergence of the method with respect to each of these parameters was investigated. The parameters used for the convergence study are given in Table 4.

The convergences of the short-circuit current density $J_{\rm SC}$, the open-circuit voltage $V_{\rm OC}$, the maximum power $P_{\rm max}$, and the fill factor $FF$ with respect to $P_{\rm deg}$, $d_{\rm z}$, and $I_{\rm deg}$, were investigated. Figure 4 shows the error in each of four electrical characteristics relative to the value for $P_{\rm deg}=9$ when $d_{\rm z}=2$ nm and $I_{\rm deg}=10$. As $P_{\rm deg}$ increases from $2$, each of the four errors is seen to decrease exponentially, with a rate approximately proportional to $\exp(-3P_{\rm deg})$. This reduction in error is seen to saturate at around $P_{\rm deg}=5$.

Figure 5 shows the error in each of four electrical characteristics relative to the value for $d_{\rm z}=1$ nm when $P_{\rm deg}=5$ and $I_{\rm deg}=10$. As $d_{\rm z}$ decreases from $20$ nm, all four errors decrease as $\mathcal{O}(d_{\rm z}^{4})$ which is suboptimal.

Finally, Fig. 6 shows the error in each of four electrical characteristics relative to the value for $I_{\rm deg}=12$ when $P_{\rm deg}=5$ and $d_{\rm z}=2$ nm. As $I_{\rm deg}$ increases, each error is seen to decrease exponentially, with a rate approximately proportional to $\exp(-4I_{\rm deg})$. This reduction in error is seen to saturate at around $I_{\rm deg}=5$.

Thus, Figs. 4–6 show that the short-circuit current density $J_{\rm SC}$, the open-circuit voltage $V_{\rm OC}$, the maximum power $P_{\rm max}$, and the fill factor $FF$, all converge at approximately the same rate with respect to $P_{\rm deg}$, $d_{\rm z}$, and $I_{\rm deg}$. The order of convergence is not optimal, as might be expected since we used polynomials of the same degree in all elements.

6 Differential Evolution Algorithm for Optimization

The differential evolution algorithm [44] is a gradient-free optimization algorithm suited for maximizing an objective function over a high-dimensional parameter space. Given the objective function $C:{\mathbb{S}}\to\mathbb{R}$, where ${\mathbb{S}}\subset\mathbb{R}^{\tilde{N}}$ is the set of all possible choices of the $\tilde{N}$ input parameters, the DEA starts by selecting $N_{\rm P}$ points to form an initial population $\textbf{P}_{0}\subset{\mathbb{S}}$. The objective function $C$ is evaluated at each of these points, with the results used by DEA in a mutation-recombination-selection process to build a new population $\textbf{P}_{1}\subset{\mathbb{S}}$. The objective function is then evaluated at each point of this new population to develop a new population $\textbf{P}_{2}\subset{\mathbb{S}}$ in the second mutation-recombination-selection step. This process continues iteratively until the objective function appears to stabilize.

The efficiency $\eta$ defined in Eq. (71) is the appropriate objective function for designing solar cells. Following suggestions in the DEA documentation [44], we fix the crossover fraction as $C_{\rm R}=0.6$. The step size $F$ used in the mutation steps is set to be randomly distributed in $\left[0.5,1\right]$ uniformly. Allowing $F$ to vary randomly with each iteration has been termed dither, and its use has been shown to improve convergence for many problems [4].

7 Numerical Test of DEA

As a model test problem, solar cell was taken to comprise a p-i-n junction made from CIGS, as shown in Fig. 1(a). The electrical parameters of CIGS [43] are given in Table 3. The $p$-type layer of thickness $L_{\rm p}=20$ nm was taken to be doped with acceptor atoms with concentration specified via $N_{\rm f}=-10^{17}$ cm*-3*, and the $n$-type layer of thickness $L_{\rm n}=20$ nm with donor atoms with concentration specified via $N_{\rm f}=10^{17}$ cm*-3*. The thickness $L_{\rm i}$ of the undoped $i$-type layer was taken as a variable for optimization. All three layers were taken to be homogeneous, with bandgaps ${\sf E}_{\rm g,p}$, ${\sf E}_{\rm g,n}$, and ${\sf E}_{\rm g,i}$ taken to be variables for optimization. Neither the radiative nor the Auger recombination mechanism was activated for the test problem.

The antireflection coating in Fig. 1(a) was ignored. We fixed $L_{\rm m}=150$ nm. The period $L_{\rm x}$ was kept variable. The height $L_{\rm g}$ and base $\zeta_{\rm g}L_{\rm x}$, $\zeta_{\rm g}\in(0,1)$, of the metallic protuberance in each period were also kept variable. The $\tilde{N}=7$ parameters chosen to be optimized in our test problem thus are: $L_{\rm x}$, $L_{\rm g}$, $\zeta_{\rm g}$, $L_{\rm i}$, ${\sf E}_{\rm g,p}$, ${\sf E}_{\rm g,n}$, and ${\sf E}_{\rm g,i}$.

Both the photonic and the electronic steps were implemented. The population number was set as $N_{\rm P}=100$. The remaining parameters used for the optimization algorithm were as described in the previous sections, together with values given in Table 5. These values were chosen to maintain reasonable accuracy, aid convergence for a wide range of solar cell designs, and to allow rapid computation of the efficiency at each step of DEA. These choices resulted in very quick evaluation of $\eta$ at a particular choice of the parameter values in roughly 6 min, with four evaluations running concurrently in MATLAB on a 20-processor (Intel Xeon Gold 61382GHz) Linux (Ubuntu 17.10) computer. The actual time per parameter set depends on the HDG solver (being longer if more homotopy steps are needed).

Figure 7 summarizes the progression of the DEA towards an optimal result. Each of panels (a)-(g) in this figure is a graph of the efficiency $\eta$ as a function of one of the parameters in the optimization. Thus, each point corresponds to a 7-dimensional vector of parameter values used by DEA. The optimal efficiency found is marked by a large red disk. The remaining panel, Fig. 7(h), shows the progress of optimization. In particular it shows the efficiency as a function of DEA step. The optimal values of the seven parameters in this study are given in Table 6 which result in $\eta=15.7\%$.

8 Concluding Remarks

We have presented details of a design tool for optimizing the efficiency of thin-film photovoltaic solar cells. The solar cell can have multiple semiconductor layers in addition to dielectric layers serving as antireflection coatings, passivation layers, and buffer layers. The solar cell is backed by a metallic grating which is periodic along a fixed direction.

The heart of the design tool is a coupled optoelectronic simulation. The photonic step of the simulation delivers the 2D variation of the electron-hole-pair generation rate inside the semiconductor layers of the solar cell. After averaging along the direction of the periodicity of the grating, the electron-hole-pair generation rate becomes the input to the electronic step of the simulation. The chief output of the electronic step is the efficiency of the solar cell. The design tool uses the differential evolution algorithm to determine the dimensions and bandgaps of the semiconductor layers to maximize the efficiency of the solar cell.

The design tool can be augmented to incorporate a biperiodic metallic grating at the cost of increased computation time [8]. Whether biperiodicity should be incorporated will depend on the consequent augmentation of the efficiency of the solar cell, and is therefore a matter of further research. Given that the period of the grating is in the 400-to-1000 nm range in order to invoke guided-wave phenomena for increased photon trapping [8, 9, 45, 46, 47, 48] but the electronic step is electrostatic in character, the averaging of the electron-hole-pair generation rate delivered by the photonic step along the periodicity directions of the grating is appropriate; elimination of that averaging will significantly increase computational time without significant gain in the electronic step.

We have begun to apply the design tool for diverse types of thin-film photovoltaic solar cells and will report our results in due course of time.

Acknowledgements

A. Lakhtakia thanks the Charles Godfrey Binder Endowment at the Pennsylvania State University for ongoing support of his research. The research of T. H. Anderson, B. J. Civiletti, and P. B. Monk was partially supported by the US National Science Foundation under grant number DMS-1619904. The research of A. Lakhtakia was partially supported by US National Science Foundation under grant number DMS-1619901.