Simulated annealing electro-photonic optimization of organic solar cells

5 downloads 265 Views 3MB Size Report
Sep 6, 2012 - Simulated annealing electro-photonic optimization of organic solar cells ...... cuit voltage Voc, as well as golden section search45 to deter-.
JOURNAL OF APPLIED PHYSICS 112, 054502 (2012)

Simulated annealing electro-photonic optimization of organic solar cells Christoph Kirscha) and Sorin Mitranb) Department of Mathematics, University of North Carolina, Chapel Hill, North Carolina 27599-3250, USA

(Received 21 March 2012; accepted 25 July 2012; published online 6 September 2012) Micro-patterned organic solar cells can exhibit enhanced light absorption properties due to a photonic crystal effect [Tumbleston et al., Appl. Phys. Lett. 94, 043305 (2009)]. Here, a three-dimensional model of light absorption and charge carrier transport in micro-patterned materials is presented and applied to the design of organic bulk heterojunction (BHJ) solar cells. Rigorous coupled wave analysis is used to simulate the multiple scattering and absorption of light in a layered solar cell device. The non-linearly coupled steady-state electric field and charge transport equations are solved iteratively by a sequence of linear partial differential equations (PDEs). Each linear PDE is discretized by an exponential upwind finite difference scheme, and the preconditioned conjugate gradient method is applied to the resulting algebraic system. The electro-photonic solver is coupled with the simulated annealing optimization algorithm to investigate the effect of micro-patterning upon performance of BHJ solar cells. Starting from the baseline configuration of a cell formed from flat layers of optimal thickness, the optimization algorithm leads to improvements of up to 15% in energy conversion efficiency by patterning both the front and back electrodes with a periodic C 2012 American Institute of Physics. ridge shape, with conformally coated layers in-between. V [http://dx.doi.org/10.1063/1.4748314] I. INTRODUCTION

Photovoltaic cells are typically made of layers of different materials. One-dimensional mathematical models1 may capture phenomena in cells consisting of flat layers, but the device simulation of solar cells with a spatial variation of material properties in the directions orthogonal to the direction of the incident light requires higher-dimensional models. Such variation of material parameters is interesting, because a periodic pattern with a characteristic length of 10!7 m (i.e., on the order of the wavelength of light) may lead to enhanced light absorption in the photoactive layer, compared to a flat configuration.2 As a result, the electron-hole pair generation rate in the solar cell will be higher and a larger current—and thus higher efficiency—is expected. Such geometric performance enhancement complements material design improvements, and needs to be investigated from both the optics and electronics point-of-views. We shall present three-dimensional mathematical models for the light scattering and absorption, as well as for the charge carrier transport in a patterned organic bulk heterojunction (BHJ) solar cell. These models are coupled via the exciton generation rate density, which is computed from the optics simulation and used as a source term in the charge transport simulation. For a given shape of the material layers, our multiphysics simulation will yield the current-voltage characteristics of the corresponding solar cell. In particular, key figures such as the short circuit current, open circuit voltage, and maximum power point can be determined, from which the fill factor of the solar cell and also its energy conversion efficiency may be computed.

The simulated annealing algorithm3,4 is used to maximize the energy conversion efficiency by varying the shape of the solar cell device. This allows us to identify optimal layer thicknesses for flat cells with given material properties, as well as optimal patterns. We also analyze the potential for further improvement and identify the main loss mechanisms remaining after shape optimization. The research code (written in C) developed to do the simulations is publicly available.5 This paper is organized as follows: In Sec. II, we review the solution of the time-harmonic Maxwell’s equations to compute the electric field by rigorous coupled wave analysis,6,7 in a multi-layered solar cell device. Post-processing of the numerical solution yields the exciton generation rate density everywhere in the photoactive layer. In Sec. III, we describe the charge transport in the organic bulk heterojunction material based on the classical semiconductor equations,8 adopting an effective medium approach. We also discuss the transformation to dimensionless variables9 and derive an expression for the net charge carrier generation rate, which includes electric field-dependent exciton dissociation.10,11 From the numerical solution, current-voltage characteristics of the solar cell device may be determined, and in particular its energy conversion efficiency can be computed. In Sec. IV, we combine our multiphysics simulation with the simulated annealing metaheuristic to identify optimal shapes. While we find no discernible performance enhancement from micro-patterning of the bulk-heterojunction layer only, combined patterning of multiple layers can lead to a 15% improvement in power per unit area. Numerical implementation details are presented in the appendices. II. LIGHT SCATTERING, ABSORPTION, AND EXCITON GENERATION

a)

[email protected]. [email protected].

b)

0021-8979/2012/112(5)/054502/13/$30.00

When light enters a solar cell device, it is scattered at the various interfaces and partially absorbed in the material 112, 054502-1

C 2012 American Institute of Physics V

054502-2

C. Kirsch and S. Mitran

J. Appl. Phys. 112, 054502 (2012)

FIG. 1. A conventional organic solar cell device with four flat layers of finite thickness (ITO, Poly(3,4-ethylenedioxythiophene) poly(styrenesulfonate) (PEDOT:PSS), P3HT:PCBM bulk heterojunction, nc-ZnO) and two semiinfinite layers (glass, Al). Incoming light is scattered at the material interfaces and partially absorbed in the photoactive layer (3). Excitons (7) are generated there which dissociate into free electrons (–) and holes (þ) that are transported to the electrodes (anode: 1, 2, cathode: 4, 5) by an electric field, !rw.

layers (Fig. 1). The exciton generation rate density G ½m!3 s!1 # is computed as an integral involving the timeaveraged power density of absorbed light in the photoactive layer, hQix ½Wm!3 #, over the angular frequency x :¼ 2pc=k, where k ½m# denotes the wavelength, and where c denotes the speed of light in vacuum G :¼

1 ð 0

hQix dx; hx !

hQix :¼

x !0 Imð!r;x ÞjEx j2 : 2

(1)

The quantities G, hQix , the complex relative permittivity !r;x , and the time-harmonic electric field Ex ½Vm!1 # vary in space; ! h denotes the reduced Planck constant and !0 denotes the vacuum permittivity. Equation (1) means that each absorbed photon generates one exciton, i.e., we have not included any multiple exciton generation from high-energy photons, an effect which has been observed in quantum dots.12 The electric field Ex in the device satisfies the secondorder time-harmonic Maxwell equations curl curl Ex ¼

x2 !r;x Ex ; c2

x > 0:

(2)

The complex relative permittivity !r;x ¼ ðnx þ ijx Þ2 , with i2 ¼ !1, has been measured13 in the range of wavelengths from 400 nm to 700 nm for all materials present in the solar cell device (Fig. 2). Equations of the form (2) must be solved for many different wavelengths in order to approximate the integral in Eq. (1). The structure of the solar cell device—layered in the direction of the incident light and periodic in the orthogonal directions—motivates the use of rigorous coupled wave analysis6,7 to solve (2), instead of a full-wave simulation. We review this method in the remainder of this section; details on the numerical implementation are given in Appendix A. A. Series expansion in terms of plane waves, Fourier transform, generalized eigenvalue problem

Throughout this section, we consider a single fixed angular frequency x > 0, and we omit all indices x to simplify notation. We restrict our considerations to light propagating in the z-direction (Fig. 1), and we assume discrete

FIG. 2. Measured optical constants13 for the six materials used in the solar cell device. Refractive index n (a) and extinction coefficient j (b) vs. wavelength k ½nm#. The extinction coefficients of glass and nc-ZnO are zero in the range of wavelengths considered here, and therefore not shown.

translational invariance of Eq. (2) on a lattice K ( R2 in the xy-plane orthogonal to the direction of propagation. This also requires all interfaces present in the solar cell device to be either parallel or perpendicular to the xy-plane; if multiple directions of propagation are to be considered, an integral over the first Brillouin zone is also involved.14 The coefficients in a series expansion of E, in terms of plane waves in (x, y)-directions, are given by ð 1 e!iG)r Eðr; zÞ dr; (3) EG ðzÞ ¼ jX0 j X0

with r :¼ ðx; yÞ> 2 X0 , where X0 denotes the primitive cell of the lattice K, containing the origin, and where G 2 K0 is any vector in the reciprocal lattice: eiG)R ¼ 1; 8 R 2 K. This reduction to X0 * R is due to Bloch’s theorem.15 With Maxwell’s equations (2), we obtain a system of linear ordinary differential equations (ODEs) for EG ¼ ðExG ; EyG ; EzG Þ> k > ¼: ðE? G ; EG Þ . The coefficients in these ODEs, ð 0 1 e!iðG!G Þ)r !r ðr; zÞ dr; (4) !r;G;G0 ðzÞ ¼ jX0 j X0

054502-3

C. Kirsch and S. Mitran

J. Appl. Phys. 112, 054502 (2012)

G; G0 2 K0 are piecewise constant due to the assumed layered structure of the solar cell device in the z-direction. Thus in each layer, the linear ODEs may be solved using the Fourier transform. For the series expansion coefficients of the orthogonal field components, E ? :¼ fE? G gG2K0 , we then obtain E ? ðzÞ ¼

1 X ikj z !ikj z ^ ? ðAþ þ A! ÞE j ; j e j e

(5)

j¼1

^ ? Þ; j 2 N are solutions of a generalized where the pairs ðkj ; E j eigenvalue problem of the form ^ ? ¼ k2 N E ^ ?: ME

(6)

The complex matrices M ¼ fM G;G0 gG;G0 2K0 and N ¼ fN G;G0 gG;G0 2K0 consist of 2 * 2 blocks with entries depending on the angular frequency x and on the layer values of !r;G;G0 (4). The square root of the generalized eigenvalues kj2 is taken such that arg kj 2 ð!p=2; p=2#; with this choice of kj , !Imðkj Þz i Reðkj Þz ^ ? the “forward” partial waves Aþ e E j in Eq. (5) j e travel in the positive z-direction, whereas the “backward” Imðkj Þz !i Reðkj Þz ^ ? e partial waves A! E j travel in the negative j e z-direction—they are evanescent if Imðkj Þ > 0. ^ ? Þ of Eq. (6) have been Once the solution pairs ðkj2 ; E j identified in each layer, it remains to match the solutions (5) across the interfaces and with the boundary conditions, in order to find the values of the amplitudes A6 j in each layer. For that purpose, the orthogonal components of the series expansion coefficients of the magnetic field B ¼ l0 H ½T# (assuming non-magnetized materials with vacuum permeability l0 ) are also used. The coefficients H? :¼ fH? G gG2K0 are related to E ? via the time-harmonic Maxwell-Faraday equation: curl E ¼ i xl0 H. Consider M > 0 layers of finite thickness with two additional, semi-infinite layers on the front and back, so that the interfaces are located at the vertical positions z0 ; …; zM . These interfaces coincide with the material interfaces in the flat case (Fig. 1), but not in general (Fig. 3). The orthogonal

components of both the electric and magnetic fields must be continuous across the interfaces between layers. Defining Y :¼ ðE ? ; H? Þ> , we have in the mth layer ~ m; ~ m eiK~ m z A Y m ðzÞ ¼ F

m ¼ 0; …; M þ 1;

(7)

^? , ~ containing the eigenvectors E with a complex matrix F m;j m ~ a diagonal complex matrix Km containing the eigenvalues ~ m containing the scatter6km;j , and with a complex vector A 6 ing amplitudes Am;j . The tilde above these quantities indicates permutation, i.e., a change in the order of summation in Eq. (5). The choice of this permutation is crucial for the numerical stability of the transfer matrix method, as discussed in Sec. II B. B. Transfer matrix method

The matching conditions at the interfaces, Y m ðzm Þ ¼ Y mþ1 ðzm Þ;

m ¼ 0; …; M;

(8)

yield a system of M þ 1 linear equations for the M þ 2 scattering amplitudes A6 m ; m ¼ 0; …; M þ 1. With Eq. (7) we also have ~

eiK m dm Y m ðzm!1 Þ ¼ Y m ðzm Þ;

m ¼ 1; …; M;

(9)

where dm :¼ zm ! zm!1 > 0 denotes the thickness of the mth layer. 6 The forward scattering matrices S0;m relate A6 0 and Am . 0;0 They can be computed recursively starting from S ¼ I (identity matrix), by alternate application of transfer and propagation steps (8) and (9). A similar recursion may be used to compute the backward scattering matrices Sm;Mþ1 6 relating A6 m and AMþ1 . The permutation mentioned in Sec. II A is used to prevent numerical instability during repeated application of propagation steps (9), when M is large.16 The appropriate permutation exchanges the positions ~ if Imðkm;j Þ < 0. of km;j and !km;j in K m With the scattering matrices, the system of linear equations (8) can be transformed into block-diagonal form,17 so that the unknown amplitudes in each finite layer, A6 m; , and the transmism ¼ 1; …; M, as well as the reflection, A! 0 sion, Aþ Mþ1 , are finally related to the boundary values ! Aþ 0 ; AMþ1 (Sec. II C) via M þ 1 equations. C. Incoming sunlight and radiation condition

FIG. 3. A solar cell device with patterned front and back electrodes and conformally coated layers. Also indicated are the interfaces used in the optics simulation (Sec. II), the computational domain X and boundaries R6 ; Rp used in the charge transport simulation (Sec. III C), as well as the geometry parameters used in the shape optimization (Sec. IV).

The incoming light is modeled as a plane wave coming from !1 and propagating in the positive z-direction. Onesun illumination is simulated by use of the AM1.5 reference solar spectral irradiance.18 We assume that the electric permittivity in the front semi-infinite layer is a positive constant: !r ðr; zÞ ¼ !r;0 > 0; z < z0 (Fig. 1). Then, the matrices M; N in the generalized eigenvalue problem (6) are block diagonal, and we obtain Einc G + 0; G 6¼ 0, and 0 þ 1 A0;j1 x pffiffiffiffiffiffi B þ C ikz inc (10) !r;0 > 0 E0 ðzÞ ¼ @ A0;j2 Ae ; k ¼ c 0

054502-4

C. Kirsch and S. Mitran

J. Appl. Phys. 112, 054502 (2012)

with Aþ 0;j1;2 2 C. Here, j1 ; j2 are the indices of eigenbasis vec> ^? tors (solutions of Eq. (6)) with E and j1 ;0 ¼ ð1; 0Þ > þ ? ^ j ;0 ¼ ð0; 1Þ . The two non-zero components A 2 C of E 0;j1;2 2 characterize the polarization of the incoming the vector Aþ 0 light, and they are only constrained by the prescribed irradiance18 I ½Wm!2 # for the angular frequency x 2 þ 2 jAþ 0;j1 j þ jA0;j2 j ¼

2xl0 I: k

(11)

In order to incorporate some of the unpolarized nature of actual sunlight into the simulation, we take an average over field computations for several linear polarizations rffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffi 2xl0 2xl0 þ þ A0;j1 ¼ I cos h; A0;j2 ¼ I sin h; (12) k k where h denotes the angle of the polarization plane with respect to the x-axis. For the other boundary value, we use A! Mþ1 ¼ 0. This corresponds to imposing a radiation condition at infinity, i.e., there should be no waves coming from þ1 and propagating in the negative z-direction. D. Exciton generation rate density, light absorption efficiency

When the scattering amplitudes have been computed, the exciton generation rate density G is evaluated at every point (r, z) in the photoactive material with (1) and (3) !0 X iðG!G0 Þ)r BHJ e IG;G0 ðzÞ; (13) Gðr; zÞ ¼ 2!h 0 0

FIG. 4. Relative reflection and absorption in each material, obtained from the optics simulation (Sec. II). (a) Flat benchmark cell (gabs ¼ 0:288) and (b) double ridge pattern with maximum energy conversion efficiency (gabs ¼ 0:329) (Sec. IV). Percentage values in legends state the cumulative relative reflection/absorption in the window 400 nm - k - 659 nm, where P3HT:PCBM is absorbing (Fig. 2).

Pabs

G;G 2K

where the z-dependence involves an integral over the angular frequency x BHJ IG;G 0 ðzÞ



1 ð

Imð!BHJ r;x ÞEx;G ðzÞ ) Ex;G0 ðzÞ dx;

0

!BHJ r;x

ð1 ð X 0

hQix ðxÞ dx dx

denotes the complex relative for G; G 2 K . In Eq. (14), permittivity of the photoactive (bulk heterojunction) mateBHJ BHJ rial. Because the imaginary part, Imð!BHJ r;x Þ ¼ 2nx jx , has compact support (Fig. 2), the integral in Eq. (14) involves only a finite interval in a practical computation. The parallel k k > electric field components Ex;G in Ex;G ¼ ðE? x;G ; Ex;G Þ are computed from the orthogonal magnetic field components H? x;G using Ampere’s circuital law, curl H x ¼ !i x!0 !r;x Ex . Thus, no further derivatives are required, and all field components necessary to evaluate the exciton generation rate (13) and (14) in the photoactive part of the mth layer fzm!1 < z < zm g (Fig. 3) are contained in Y m (7), m ¼ 1; …; M. Other quantities of interest, which are not required for the coupling with the charge transport simulation (Sec. III), can also be computed from the optics simulation results, such as the reflection, transmission, and absorption in each material (Fig. 4). We denote the domain occupied by the photoactive material by X ( X0 * R. The power per unit area of light absorbed in the photoactive material is given by

½Wm!2 #:

(15)

With the reference solar spectral irradiance18 Ik (for wavelengths k , 400 nm), we also compute the light absorption efficiency of the photoactive region

(14)

0

0

1 :¼ jX0 j

gabs

Pabs :¼ ; Pin

Pin :¼

1 ð

Ik dk ’ 954 Wm!2 :

(16)

400 nm

III. CHARGE TRANSPORT, EXCITON DISSOCIATION, AND CHARGE CARRIER RECOMBINATION

We use the exciton generation rate density G computed in Sec. II as a source term in the charge transport simulation. The photoactive material in an organic solar cell device is regarded as an effective insulator, with a band gap Egap ½J# given by the difference between the lowest unoccupied molecular orbital (LUMO) level of the electron acceptor material (such as [6,6]-phenyl-C61-butyric acid methyl ester (PCBM)) and the highest occupied molecular orbital (HOMO) level of the electron donor material (such as P3HT).19 Electric current is produced as photogenerated excitons dissociate into free charge carriers and are transported to the electrodes by an electric field. For the charge transport in patterned solar cell devices such as the one illustrated in Fig. 3, we extend a one-dimensional metal-insulator-metal model1 to three space dimensions.

054502-5

C. Kirsch and S. Mitran

J. Appl. Phys. 112, 054502 (2012)

We note that the model described below includes charge transport in the photoactive material only. An approach to incorporate the non-active materials into the charge transport model as effective resistors has been proposed recently.20 A. Coupled charge carrier transport equations

where U~ :¼ U=Nint ½s!1 #, and where the characteristic Debye length kD ½m# is given by sffiffiffiffiffiffiffiffiffi !Vt kD :¼ : (26) qNint

The semiconductor equations8 involve the electric potential w ½V# and the charge carrier number densities n, p ½m!3 # of electrons and holes, respectively. These coupled partial differential equations (PDEs) describe the drift and diffusion of each charge carrier species in the material at the continuum scale, and they also relate the total charge density to the electric field. The steady-state semiconductor equations in nondoped material are given by

A number of monographs cover semiconductor device physics and modeling;21,22 for a mathematical analysis of Eqs. (17)–(19), we refer the reader to Jerome.23,24 For the time-dependent version of Eqs. (17)–(19), exponential convergence to the unique steady state has been shown analytically.25,26

div ð!rwÞ ¼ qðn ! pÞ;

Photogenerated excitons may dissociate into free charge carriers (electron-hole pair) if a donor/acceptor interface is located within a distance of the order of the exciton diffusion length (.10!8 m) from the point of generation. In a polymerfullerene blend such as P3HT:PCBM, the characteristic length scales of both phases are comparable to the exciton diffusion length,27 and therefore these blends are interesting candidates for the photoactive layer in a solar cell device. Thus, we neglect the exciton diffusion in the photoactive layer, and we consider the following local generation and loss mechanisms

div Jn ¼ Uðjrwj; n; pÞ; div Jp ¼ Uðjrwj; n; pÞ;

(17)

Jn :¼ nln rw ! Dn rn;

(18)

Jp :¼ !plp rw ! Dp rp:

(19)

In Eqs. (17)–(19), the (static) electric permittivity and the elementary charge are denoted by ! ½Fm!1 # and q ½C#, respectively. The charge carrier mobilities ln;p ½m2 ðVsÞ!1 # and diffusion coefficients Dn;p ½m2 s!1 # satisfy the Einstein relation20 Dn;p ¼ ln;p Vt ;

Vt :¼

kB T ½V#; q

(20)

where the thermal voltage Vt involves the Boltzmann constant kB ½JK!1 # and the absolute temperature, T [K]. The net charge carrier generation rate density U ½m!3 s!1 # is discussed in Sec. III B, and the boundary conditions are presented in Sec. III C. At thermal equilibrium, the electron and hole fluxes Jn ; Jp ½m!2 s!1 # vanish, and we find exponential relationships between the charge carrier number densities n, p and the electric potential w. This motivates a transformation to dimensionless variables (u, v, w), which are related to the natural variables ðw; n; pÞ via9 w ¼ uVt ;

n ¼ Nint veu ;

p ¼ Nint we!u :

(21)

In Eq. (21), the intrinsic charge carrier number density Nint ½m!3 # depends on the effective densities of states Nc ; Nv ½m!3 # in the conduction and valence band, respectively, as well as on the band gap; it is defined by # $ pffiffiffiffiffiffiffiffiffiffi 1 Egap : (22) Nint :¼ Nc Nv exp ! 2 qVt After division by Nint , the steady-state semiconductor equations (17)–(19) take the form FðuÞ ¼ 0, where u :¼ ðu; v; wÞ> and F :¼ ðFu ; Fv ; Fw Þ> with Fu ðuÞ :¼ div ðk2D ruÞ ! veu þ we!u ;

(23)

~ v; wÞ; Fv ðuÞ :¼ div ðDn eu rvÞ þ Uðjruj;

(24)

~ v; wÞ; Fw ðuÞ :¼ div ðDp e!u rwÞ þ Uðjruj;

(25)

B. Net charge carrier generation rate

ground state

G

!! kf X

X free charge excitons kdiss !! X carriers n; p:

(27)

R

In Eq. (27), X ½m!3 # denotes the exciton number density, G denotes the exciton generation rate density (Sec. II), R ½m!3 s!1 # denotes the charge carrier recombination rate density, and kdiss ; kf ½s!1 # denote the exciton dissociation and decay rates, respectively. At steady state, these rate densities must be balanced ðkdiss þ kf ÞX ¼ G þ R:

(28)

Therefore, we may obtain the net charge carrier generation rate density U from a convex combination of exciton generation and charge carrier recombination ð27Þ

ð28Þ

U ¼ kdiss X ! R ¼ Pdiss G ! ð1 ! Pdiss ÞR

(29)

with the dissociation probability (the fraction of excitons which dissociates into free charge carriers) Pdiss :¼

kdiss : kdiss þ kf

(30)

The bimolecular model28 is commonly used for the recombination rate density R 2 Þ; R :¼ cðnp ! Nint

% & q c ¼ min ln ; lp : !

(31)

Several different expressions have been proposed for the Langevin recombination constant c ½m3 s!1 # in the case of

054502-6

C. Kirsch and S. Mitran

J. Appl. Phys. 112, 054502 (2012)

polymer-fullerene blends.29,30 The exciton dissociation rate kdiss ½s!1 # depends on the electric field strength and on the electron-hole pair separation distance x > 0 ½m#, and is given by10,11 # $ 3c Eb ðxÞ I1 ð2bÞ : (32) kdiss ðjrwj; xÞ ¼ exp ! 4px3 qVt b Equation (32) involves the first-order modified Bessel function of the first kind, I1 sffiffiffiffiffiffiffiffiffiffiffiffi 1 I1 ð2bÞ X b2k qjrwj ; b :¼ ¼ ; (33) k!ðk þ 1Þ! b 4p!Vt2 k¼0

odic boundary conditions at Rp ; the boundary conditions imposed at R6 are specified below. We still use the metal-insulator-metal picture, so that the following five positive parameters describe the energy levels in the materials (Fig. 5(a)): the electron affinity v, the work function /, and the band gap Egap of the insulator, as well as the work functions /a and /c of the anode and cathode, respectively, with /c < /a . When the thermal equilibrium is established after contact, the Fermi levels are aligned, and Ec ; Ev vary across the insulator accordingly (Fig. 5(b)). Upon applying an external voltage Va ½V# (bias), the thermal equilibrium is disturbed and we consider quasiFermi levels EFn ; EFp for the electrons and holes, respectively (Fig. 5(c)). They are assumed to be pinned to the metal

and the coulombic electron-hole pair binding energy is given by Eb ðxÞ :¼ q2 =ð4p!xÞ ½J#. The disorder in the polymer-fullerene blend is taken into account by assuming a random electron-hole pair separation distance with a Maxwell-Boltzmann distribution31 # 2$ 4x2 x (34) f ðx; aÞ :¼ pffiffiffi 3 exp ! 2 ; x > 0: a pa The parameter a > 0 ½m# is the location of the maximum of the probability density function f ð ) ; aÞ; the expected value of the pffiffiffi electron-hole pair separation distance is given by 2a= p. Then, the exciton dissociation rate kdiss and thus the dissociation probability (30) also become random variables. The value used for Pdiss in Eq. (29) is the expected value Pdiss ðjrwj; a; kf Þ :¼

1 ð 0

kdiss ðjrwj; xÞ f ðx; aÞ dx: kdiss ðjrwj; xÞ þ kf

(35)

All material properties of P3HT:PCBM which we have used in our simulations are listed in Table I. C. Boundary conditions

The semiconductor equations (17)–(19) have to be solved in the domain X occupied by the polymer-fullerene blend. This domain is bounded by the electrodes and by the faces of the right prism X0 * R, where X0 denotes the primitive cell of the lattice used in the optics simulation (Sec. II A). The boundary of X is decomposed into @X ¼ Rþ [ R! [ Rp , where Rþ and R! denote the cathode and anode boundaries, respectively (Fig. 3). We impose periTABLE I. Material parameters used for the charge transport simulation in P3HT:PCBM, with references to the equations where they first appear. !

3:4 !0 Fm!1

(17)

ln

2 * 10!7 m2 ðVsÞ!1

(18)

10!8 m2 ðVsÞ!1 300 K

(19) (20)

2:5 * 1025 m!3 0.95 eV 2 * 104 s!1 1:8 * 10!9 m

(22) (22) (27) (34)

lp T pffiffiffiffiffiffiffiffiffiffi Nc Nv Egap kf a

FIG. 5. Illustration of the energy levels in the metal-insulator-metal model used for the charge transport simulation: the five parameters v; Egap ; /; /a ; /c describing the energetic properties of the metals and of the insulator (a), the thermal equilibrium after contact at R6 (b), and the situation under forward bias (c). The band bending near the contacts is not shown in these illustrations.

054502-7

C. Kirsch and S. Mitran

J. Appl. Phys. 112, 054502 (2012)

Fermi levels at each boundary and to be separated across the insulator by the applied voltage EFn jR! ¼ EFp jR! ;

EFn jRþ ¼ EFp jRþ ;

(36)

ðEi ! EFp ÞjR! ¼ /a ! /;

(37)

ðEFn ! Ei ÞjRþ ¼ / ! /c ;

(38)

EFn jRþ ! EFp jR! ¼ qVa :

(39)

Iþ increases with the applied voltage Va , and it changes sign at the open circuit voltage Voc ½V#. For the simulation of the current-voltage characteristic of a solar cell device, we divide the current through the cathode, Iþ ðVa Þ, by the area of the primitive cell X0 used in the optics simulation (Sec. II A), i.e., we plot the current per unit area, jðVa Þ :¼

Iþ ðVa Þ jX0 j

½Am!2 #;

(45)

Equations (36)–(39) translate to Dirichlet boundary conditions for the unknowns ðw; n; pÞ if we assume a Boltzmann distribution for the occupation probability of electron states at the boundaries32

vs. the applied voltage Va (Fig. 6). The electric power per unit area is then given by

1 wjRþ ! wjR! ¼ ð/a ! /c Þ ! Va ; q # $ p '' Nint '' /a ! / ; ' ¼ ' ¼ exp Nint R! n R! qVt # $ n '' Nint '' / ! /c : ' þ¼ ' þ ¼ exp p R qVt Nint R

which is positive for 0 < Va < Voc , under illumination. Together with the short circuit current Jsc :¼ !jð0Þ ½Am!2 #, the fill factor FF and the energy conversion efficiency g of the solar cell device can be computed

(40) (41)

PðVa Þ :¼ !Va jðVa Þ

FF :¼ (42)

The unspecified additive constant in the electric potential w may be neglected because the semiconductor equations (17)–(19) involve only the gradient, rw; therefore, we may add the equation wjRþ þ wjR! ¼ 0. The boundary conditions (40)–(42) are further simplified by assuming that /a ! / ¼ / ! /c ¼ 12 Egap . Because the boundary conditions for the charge carrier number densities (41) and (42) are independent of the applied voltage, they have been called locked density boundary conditions, and their validity in polymer-fullerene solar cell device models has recently been questioned.33 We would like to emphasize that the three-dimensional charge transport/recombination model described in Secs. III A and III B may be completed with any type of boundary condition for the metal-insulator contacts.34 D. Post-processing and efficiency factors

Pout ; Jsc Voc

g :¼

Pout ; Pin

½Wm!2 #;

Pout :¼ maxVa PðVa Þ;

(46)

(47)

where the input power per unit area, Pin , was defined in Eq. (16). The energy conversion efficiency g is to be maximized via shape optimization (Sec. IV). Other interesting quantities obtained from postprocessing include the exciton dissociation and charge carrier collection efficiencies, which are defined by gdiss :¼

hPdiss Gi ; hGi

gcoll :¼

hUi : hPdiss Gi

(48)

The angled brackets h ) i in Eq. (48) denote spatial averages over the domain X occupied by the polymer-fullerene blend. We validate our charge transport simulation by reproducing the results obtained by Koster et al.1 for a flat PPV:PCBM layer of 120 nm thickness, with a constant exciton generation rate density. Our simulation results are stated in Table II— they are in good agreement with these earlier results, which were also verified by experimental measurements.1

From the solution ðw; n; pÞ of the semiconductor equations (17)–(19) with the boundary conditions (40)–(42), several quantities of practical interest may be computed. From the electric current density, j :¼ !qðJn ! Jp Þ

½Am!2 #;

(43)

the electric current through the electrode boundaries may be computed by integration ð 6 (44) I :¼ j ) n ds ½A#; R6

where n denotes the outward unit normal vector on @X. The currents through the cathode and anode are of equal magnitude with opposite signs. We have Iþ < 0 under illumination at short circuit (Va ¼ 0), indicating that more electrons than holes are flowing outwards through the cathode. The current

FIG. 6. Current-voltage characteristics of the flat benchmark cell (dashed line) and of the double ridge pattern with maximum energy convergence efficiency (solid line) (Sec. IV), obtained from the electro-photonic simulation (Secs. II and III). Current density j vs. applied voltage Va . The short circuit current Jsc , open circuit voltage Voc , and fill factor FF (Sec. III D) for the two solar cell devices are given in the legend.

054502-8

C. Kirsch and S. Mitran

J. Appl. Phys. 112, 054502 (2012)

TABLE II. Simulation results for a flat PPV:PCBM layer of 120 nm thickness, with constant exciton generation rate density.1 Applied voltage Va and current per unit area j (45) at short circuit, maximum power, and open circuit (Sec. III D). Also reported are the values of the dissociation efficiency gdiss and of the recombination loss, which is computed from the collection efficiency gcoll (48).

Short circuit Maximum power Open circuit

Va ½V#

!jðVa Þ½Am!2 #

gdiss ½%#

1 ! gcoll [%]

0 0.658 0.844

28.8 19.5 0

60.2 50.4 46.5

6.37 23.7 97.5

IV. SHAPE OPTIMIZATION

With the electro-photonic simulation described in Secs. II and III, we are now able to compute the energy conversion efficiency of solar cell devices under realistic illumination (47), given the periodic structure of the material layers and certain optical and electrical material properties (Fig. 2, Table I). We combine this multiphysics simulation with the simulated annealing metaheuristic (Appendix C) to find the cell geometry which maximizes the energy conversion efficiency for given material properties. This optimization aims at identifying photonic crystal structures which result in solar cell devices with a higher energy conversion efficiency than flat cells. A. Results for flat cells and flat benchmark

We begin with a one-dimensional parameter space, where only the P3HT:PCBM layer thickness h3 is variable (Fig. 1). The other layer thicknesses are fixed to ðh1 ; h2 ; h4 Þ ¼ ð178; 50; 0Þ nm. In Fig. 7, we plot the light absorption and charge carrier collection efficiencies gabs (16), gdiss (48) at the maximum power point, as well as the energy conversion efficiency g (47) vs. h3 . While more light is absorbed in the photoactive material as its volume increases, the loss of charge carriers due to recombination increases, too. The optimal active layer thickness is found at h3 ¼ 77 nm ðg ¼ 00393Þ. Similar

FIG. 7. Simulation results for a flat layered solar cell device with variable thickness of the P3HT:PCBM layer, h3 (Fig. 1). The other layer thicknesses are fixed to ðh1 ; h2 ; h4 Þ ¼ ð178; 50; 0Þ nm. Light absorption and charge carrier collection efficiencies gabs (16), gdiss (48) at the maximum power point, as well as energy conversion efficiency g (47) vs. h3 . Optimal thickness found at h3 ¼ 77 nm ðg ¼ 00393Þ.

values for the optimal P3HT:PCBM layer thickness in flat cells have been reported before.35,36 Next, we consider a two-dimensional parameter space, where both the P3HT:PCBM and the nc-ZnO layer thicknesses h3 ; h4 are variable. The nc-ZnO material is non-absorbing in the range of wavelengths considered here (Fig. 2) and thus acts as an optical spacer. Fig. 8 illustrates how the simulated annealing metaheuristic automatically focuses in the “most interesting” region of the parameter space. The optimal parameter pair is found at ðh3 ; h4 Þ ¼ ð53; 27Þ nm ðg ¼ 0:0400Þ. We also note from Fig. 8 that the energy conversion efficiency is mainly governed by the sum h3 þ h4 , an observation previously made also for the short circuit current.37 Patterned cells shall be compared with the best possible flat cell for this solar cell device model (Secs. II and III) and given material properties (Fig. 2, Table I). Therefore, we now vary all four layer thicknesses, h1 ; …; h4 (Fig. 1) and show the results in Fig. 9. The gray dots correspond to points in the parameter space generated by the simulated annealing metaheuristic, shown together with the values of the energy conversion efficiency g. We find the maximum g ¼ 0:0433 at ðh1 ; h2 ; h3 ; h4 Þ ¼ ð78; 39; 51; 32Þ nm. This is our flat benchmark cell. B. Ridge pattern in the back electrode

We consider patterned solar cell devices with a planar front electrode (flat glass, ITO, and Poly(3,4-ethylenedioxythiophene) poly(styrenesulfonate) (PEDOT:PSS) layers), and a periodic ridge pattern imprinted into the P3HT:PCBM bulk heterojunction material. The space between the ridges is filled with nc-ZnO, and a flat aluminum layer closes the device. For one particular pattern of this kind, the electric field lines and the net charge carrier generation rate at maximum power, obtained from the electro-photonic simulation, are

FIG. 8. Points in the two-dimensional parameter space (variable P3HT:PCBM and nc-ZnO layer thicknesses, h3 ; h4 (Fig. 1)) generated by the simulated annealing metaheuristic. Other layer thicknesses fixed to ðh1 ; h2 Þ ¼ ð178; 50Þ nm. Gray values indicate energy conversion efficiency g (47) of the corresponding solar cell device. Optimal layer thicknesses found at ðh3 ; h4 Þ ¼ ð53; 27Þ nm ðg ¼ 00400Þ.

054502-9

C. Kirsch and S. Mitran

J. Appl. Phys. 112, 054502 (2012)

FIG. 9. Gray dots: energy conversion efficiency g (47) vs. flat layer thicknesses h1 ; …; h4 (Fig. 1). Optimal values (flat benchmark): ðh1 ; h2 ; h3 ; h4 Þ ¼ ð78; 39; 51; 32Þ nm ðg ¼ 00433Þ (Sec. IV A). Black pluses: g vs. volume per unit area of each material in a solar cell device with a double ridge pattern (Fig. 3). Maximum g ¼ 0:0497 found at ðw1 ; w2 ; h1 ; h2 ; h3 ; h4 ; h5 Þ ¼ ð0; 0; 168; 4; 0; 28; 4Þ nm (Sec. IV C).

shown in Fig. 10. While these patterned solar cell devices may absorb more light than the flat cell with the same amount of P3HT:PCBM bulk heterojunction material,2 they also suffer from a high recombination loss due to a charge carrier sink (U < 0) in the photoactive layer. In our shape optimization efforts, we have not found any pattern of this kind with a higher energy conversion efficiency than the flat benchmark cell. C. Ridge pattern in both front and back electrodes

Instead, we propose the pattern outlined in Fig. 3, where ridges of height h1 are etched into both the glass front layer

and the aluminum back layer. The remaining layers with thicknesses h2 ; …; h5 are coated conformally in-between. The widths of the two ridges, w1 ; w2 , are also variable. Such a pattern should offer the possibility of increased light absorption in the photoactive material due to a photonic crystal effect,14 while at the same time providing a fairly uniform travel distance for the charge carriers. We use the simulated annealing metaheuristic on the 7-dimensional parameter space, and we plot the material volumes per unit area as black pluses in Fig. 9, together with the corresponding value of the energy conversion efficiency g (47). It turns out that the optimal pattern allows for about three times more P3HT:PCBM material than the flat benchmark cell. Together with a much smaller volume of the front layers (ITO and PEDOT:PSS) the pattern increases the light absorption efficiency in the photoactive material by 14% compared to the flat benchmark cell (Fig. 4). The energy conversion efficiency g is improved by 15%. Optimal parameter values are given by ðw1 ; w2 ; h1 ; h2 ; h3 ; h4 ; h5 Þ ¼ ð0; 0; 168; 4; 0; 28; 4Þ nm. This means pure ITO and ncZnO ridges of height 168 nm and width 8 nm on the front and back, respectively, with a P3HT:PCBM layer of 28 nm thickness in-between (no PEDOT:PSS). The electric field lines in P3HT:PCBM for this pattern are shown in Fig. 11. Comb structures of this kind have also been found in blend morphology optimization.38

D. Potential for further improvement

FIG. 10. Electric field lines (solid) and labeled contour lines (dotted) of the net charge carrier generation rate density U ½1027 m!3 s!1 # at maximum power (obtained from the electro-photonic simulation), for a solar cell device with a ridge pattern imprinted into the back side of the P3HT:PCBM bulk heterojunction material and filled with nc-ZnO (Sec. IV B). The unit of length is ½nm#. This solar cell device suffers from a high recombination loss due to a charge carrier sink (U < 0) in the photoactive layer.

The energy conversion efficiency g of an organic bulk heterojunction solar cell device may be written as a product of five factors g ¼ gabs ggen gdiss gcoll gout ;

(49)

where the exciton generation and output efficiencies are defined by

054502-10

C. Kirsch and S. Mitran

J. Appl. Phys. 112, 054502 (2012) TABLE III. Efficiency factors computed for the flat benchmark cell (Sec. IV A) and for the best pattern (Sec. IV C). The main energy loss mechanisms remaining after shape optimization are the light absorption and the generation of excitons from absorbed photons. Flat benchmark

Best pattern

gabs (16)

0.288

0.329

ggen (50) gdiss (48) gcoll (48) gout (50)

0.206 0.839 0.929 0.935

0.200 0.880 0.963 0.891

g (47)

0.0433

0.0497

(Fig. 4); a substantial increase of the absorption efficiency gabs can therefore be achieved by using materials which absorb at higher wavelengths. Materials with long and thin absorption tails seem particularly promising, since the photonic crystal structure may lead to large gains in absorption especially in the tail region. Furthermore, it is necessary to make use of effects which generate multiple excitons from a single photon,12 in order to increase the energy conversion efficiency of exciton generation from absorbed photons, ggen . V. CONCLUSIONS

FIG. 11. Electric field lines in P3HT:PCBM for the optimal double ridge pattern (Sec. IV C). The unit of length is ½nm#.

ggen :¼

jXj qVmp jX hGi 0j

Pabs

;

gout :¼

Pout jXj qVmp jX hUi 0j

:

(50)

In Eq. (50), Vmp ½V# denotes the applied voltage at the maximum power point. ggen and gout relate the exciton generation rate to the power of absorbed light and the maximum output power to the charge carrier collection rate, respectively. Computed efficiency factors for both the benchmark flat cell (Sec. IV A) and the best patterned cell (Sec. IV C) are shown in Table III. The main energy loss mechanisms remaining after shape optimization are the light absorption (loss due to reflection, transmission, and absorption in non-active materials) and the exciton generation from absorbed photons (loss due to excess photon energy beyond the threshold required to produce an exciton): these two loss mechanisms alone limit the efficiency of our solar cell device at gabs ggen ’ 6 %. The absorption in the photoactive material within the window of wavelengths 400 nm - k - 659 nm is already high

Using a combined optics and charge transport simulation we have computed the energy conversion efficiency for given solar cell geometries and material parameters, under one-sun illumination. We have combined this simulation with the simulated annealing metaheuristic to identify optimal geometry parameters. A periodic pattern consisting of ridges of ITO and nc-ZnO, with a thin layer of P3HT:PCBM coated conformally in-between, improves the energy conversion efficiency by 15% compared to the best possible flat cell for the given material parameters (also found by simulated annealing). Further investigation showed that after shape optimization, the main energy loss occurs in the light absorption and in the exciton generation from absorbed photons. These two loss mechanisms already limit the energy conversion efficiency of our solar cell devices at .6 %. ACKNOWLEDGMENTS

The authors would like to thank John R. Tumbleston (now at North Carolina State University) and Rene Lopez of the Department of Physics and Astronomy at the University of North Carolina at Chapel Hill for optical measurements13 and fruitful discussions throughout the project. This research was funded by the National Science Foundation within the Solar Energy Initiative (DMR 0934433). APPENDIX A: OPTICS SIMULATION

The rigorous coupled wave analysis (Sec. II) allows for a mesh-free simulation of the light scattering and absorption in the solar cell device. We restrict the reciprocal lattice K0 to the N > 0 shortest lattice vectors; then the generalized eigenvalue problem (6) is of size 2N, and thus the series (5) becomes a finite sum. If the complex relative permittivity !r

054502-11

C. Kirsch and S. Mitran

J. Appl. Phys. 112, 054502 (2012)

does not vary in the horizontal directions, such as in flat layered devices (Fig. 1), then N ¼ 1 (K0 ¼ f0g) is sufficient, and the problem becomes one-dimensional. For N > 1, we compute the layer values of !r;G;G0 (4) analytically, since for a pattern such as the one illustrated in Fig. 3, !r is piecewise constant in (x, y)-directions, in each layer. We have used LAPACK39 to do the matrix computations described in Secs. II A and II B. The computational cost of the QZ algorithm used for solving the generalized eigenvalue problem (6) and of the LU decompositions used to solve the linear systems in the transfer matrix method (for 10 different / / / linear polarization angles h 2 f0 ; 10 ; …; 90 g (12)) is OðN 3 Þ floating point operations, as N ! 1. These problems need to be solved for a large number of wavelengths in order to approximate the integral in Eq. (13). We have used N ¼ 7 for the ridge patterns illustrated in Fig. 3. The exciton generation rate (13) needs to be evaluated at each point in the computational mesh used for the charge transport simulation (Appendix B); the trapezoid rule was used for the numerical approximation of the integral . APPENDIX B: CHARGE TRANSPORT SIMULATION 1. Iterative scheme

Iteration is applied to transform the nonlinear system of three PDEs F(u) ¼ 0 (23)–(25) into a sequence of linear partial differential equations. A Gauss-Seidel-type iteration leads to a sequence of scalar nonlinear PDEs. The two charge transport equations (24) and (25) become linear upon decoupling, and a single Newton step is used to solve (23) in each iteration.40 Thus, we solve in each step the following sequence of three linear, second-order, elliptic PDEs for u0 , given the current iterate u divðk2D ru0 Þ þ Aðu; v; wÞu0 ¼ f ðu; v; wÞ; 0

divðDn eu rv0 Þ þ Bðu0 Þv0 w ¼ gðu0 Þ; 0

divðDp e!u rw0 Þ þ Bðu0 Þv0 w0 ¼ gðu0 Þ;

(B1) (B2) (B3)

we have opted for the more robust Gauss-Seidel-type method. Iteration (B1)–(B3) is run until the 2-norm of the nonlinear residual, kFðuÞk2 , falls below a prescribed tolerance. 2. Discretization

As we have mentioned in Sec. III C, the computational domain X is bounded by the electrode boundaries R6 and by the faces of the right prism X0 * R, where X0 ( R2 denotes the primitive cell of the lattice K used in the optics simulation (Sec. II A). We have found it convenient to discretize a part of the full right prism first, and to mask out the mesh cells which do not contain polymer-fullerene blend material (Fig. 12). These cells are used to store the boundary values. The second-order differential operators in Eqs. (B1)– (B3) are of the form divðDr/Þ, with the unknown function / and with a variable coefficient D > 0. An approximation of the divergence at the cell centers using centered finite differences requires the values of the flux Dr/ at the cell faces. Assuming a uniform mesh of cubic cells with edge length h > 0, we use the following approximation of the flux in one space dimension 0

ðD/ Þiþ1=2

0

1 ’@ h

xð iþ1 xi

1!1 1 /iþ1 ! /i ; dxA h DðxÞ

(B8)

where /i denote the values of the unknown function / at the cell centers xi . In higher space dimensions, we may simply discretize in the same way in each direction, because the differential operator separates in cartesian coordinates. The discretization (B8) applied to Eqs. (B1)–(B3) leads to an exponential upwind scheme, which is commonly used to solve the semiconductor equations.41 The discrete versions of (B1)–(B3) each give rise to a sparse linear system (no more than 7 non-zero entries in each matrix row in 3D), which we solve iteratively using the

where the coefficients and source terms are given by Aðu; v; wÞ ¼ !veu ! we!u ; f ðu; v; wÞ ¼ veu ! we!u þ Aðu; v; wÞu; ( ) Bðu0 Þ ¼ ! 1 ! Pdiss ðVt jru0 j; a; kf Þ cNint ;

gðu0 Þ ¼ Bðu0 Þ ! Pdiss ðVt jru0 j; a; kf ÞG=Nint ;

(B4) (B5) (B6) (B7)

with the dissociation probability Pdiss (35). Alternatively, a Newton-type iteration can be applied to F(u) ¼ 0, which linearizes F (23)–(25) about the current iterate u so that a system of linear PDEs needs to be solved in each step. While the Newton-type iteration is expected to converge faster than the Gauss-Seidel-type iteration (B1)–(B3), it may also fail due to a bad starting point or inexact linearization of F, for example. Since we need to compute solutions for several different boundary values (to find the maximum power point) and also for various domain geometries (shape optimization),

FIG. 12. Illustration of the computational mesh for a patterned solar cell device (cf., Fig. 3). Mesh cells which do not contain polymer-fullerene blend material are masked out (*); they are used to store the electrode boundary values (Sec. III C).

054502-12

C. Kirsch and S. Mitran

conjugate gradient method.42 We also employ a simple preconditioner given by the diagonal approximate inverse,43 to accelerate convergence. The evaluation of the residual Fðu0 Þ at the end of each iteration step is straightforward, because the operators appearing in Eqs. (23)–(25) are discretized anyway during the iteration. 3. Dissociation probability

J. Appl. Phys. 112, 054502 (2012)

density is chosen such that its location is the current state x, and there may be additional parameters involved, such as the variance. These additional parameters may also change during the course of the MH algorithm, to reduce the mixing time.48 The algorithm constructs a sample sequence of a new Markov chain fX k gk2N0 with transition probability * + pðx0 Þqðxjx0 Þ 0 0 ; (C1) pðx jxÞ ¼ qðx jxÞ min 1; pðxÞqðx0 jxÞ

The dissociation probability Pdiss ðVt jru0 j; a; kf Þ needs to be evaluated before step (B2) in the iteration, from the solution u0 of Eq. (B1). We use centered finite differences to compute the local electric field strength from the cell-center values of u0 . The integral (35) is transformed into an integral over the whole real line by writing x ¼ expðyÞ; y 2 R. Then, the integrand exhibits a doubly-exponential decay as y ! 61. Thus, we may use the trapezoid rule for the approximation, which converges exponentially in this case.44 The dissociation rate (32) is evaluated using the fastconverging series in Eq. (33), by truncation.

by accepting proposal states x0 sampled from qð ) jxÞ with a probability that depends on the ratio of probabilities of the states x and x0 . It is easily verified that p is indeed a stationary distribution of the Markov chain fX k gk2N0 . Therefore, after a sufficiently long “burn-in” period, the MH algorithm generates samples from p. Introducing the temperature T > 0 again, the typical choice for the probability density function fT of Y in simulated annealing is proportional to the Boltzmann factor,

4. Post-processing

with the thermal beta bT :¼ ðkB TÞ!1 , where kB denotes the Boltzmann constant. Notice that the value of the normalization constant in Eq. (C2) is not important, because only ratios of values of pT are used in Eq. (C1). This is a benefit of the MH algorithm, because the partition function is difficult to compute in general. pT is thus interpreted as the probability density function of states of a system in thermodynamic equilibrium at temperature T, where the thermodynamic potential is given by the objective function g. In simulated annealing, the temperature T is gradually reduced according to a “cooling schedule” and the MH algorithm is used to sample pT for each value of T. Notice that for T > 0, proposal states x0 with gðx0 Þ > gðxÞ are occasionally accepted, which should prevent the algorithm from getting stuck in a local minimum of g. As T decreases, so does the expectation EpT ½gðXÞ#, and therefore we expect to reach a global minimum of g as T ! 0.

We vary the applied voltage Va in a one-dimensional search and use the false position method to find the open circuit voltage Voc , as well as golden section search45 to determine the maximum power Pout (47) without requiring derivatives of j (45). APPENDIX C: SIMULATED ANNEALING

The simulated annealing metaheuristic3,4 is used in global optimization, such as for maximizing the power output per unit area of a solar cell device in our application, with a variable cell geometry. Simulated annealing is based on the Metropolis-Hastings (MH) algorithm,46,47 which is a Markov chain Monte Carlo method for obtaining random samples from a (typically complicated) probability distribution. For N > 0 geometry parameters (such as layer thicknesses and ridge width (Fig. 3)), we consider a real-valued, N-dimensional random vector X, as well as an objective function g : RN ! R which is to be minimized (here, the negative of the maximum power per unit area, !Pout (47), for a solar cell geometry (“state”) described by realizations x 2 RN of X). In our application, evaluating g(x) corresponds to running the electro-photonic simulation to find the maximum power point for a cell geometry described by x. We define the random variable Y: ¼ g(X) for which we prescribe a probability density function fT . The function fT depends on the “temperature” T > 0 in the simulated annealing metaheuristic. The probability density function of X is then given by pT :¼ fT g, and the MH algorithm is used to construct a sample sequence of a Markov chain with equilibrium distribution pT . We assume a fixed T > 0 for the moment and omit the index T. The MH algorithm employs an arbitrary Markov chain with transition probability qðx0 jxÞ; x; x0 2 RN . For every fixed x 2 RN , the function qð ) jxÞ is a probability density function in the state space, which is called the proposal density. Typically, the proposal !

fT ðyÞ / expð!bT yÞ ) pT ðxÞ / expð!bT gðxÞÞ;

1

(C2)

L. J. A. Koster, E. C. P. Smits, V. D. Mihailetchi, and P. W. M. Blom, “Device model for the operation of polymer/fullerene bulk heterojunction solar cells,” Phys. Rev. B 72, 085205 (2005). 2 J. R. Tumbleston, D.-H. Ko, E. T. Samulski, and R. Lopez, “Electrophotonic enhancement of bulk heterojunction organic solar cells through photonic crystal photoactive layer,” Appl. Phys. Lett. 94, 043305 (2009). 3 S. Kirkpatrick, C. D. Gelatt, Jr., and M. P. Vecchi, “Optimization by simulated annealing,” Science 220, 671–680 (1983). 4 " y, “Thermodynamical approach to the traveling salesman problem: V. Cern# An efficient simulation algorithm,” J. Optim. Theory Appl. 45, 41–51 (1985). 5 C. Kirsch, “Solar cell research code,” http://mitran-lab.amath.unc.edu:8081/ subversion/PolymerPhotovoltaic/C-ChargeTransport/code (2012), written in C. 6 M. G. Moharam and T. K. Gaylord, “Rigorous coupled-wave analysis of planar-grating diffraction,” J. Opt. Soc. Am. 71, 811–818 (1981). 7 S. G. Tikhodeev, A. L. Yablonskii, E. A. Muljarov, N. A. Gippius, and T. Ishihara, “Quasiguided modes and optical properties of photonic crystal slabs,” Phys. Rev. B 66, 045102 (2002). 8 W. van Roosbroeck, “Theory of the flow of electrons and holes in germanium and other semiconductors,” Bell Syst. Tech. J. 29, 560–607 (1950). 9 J. W. Slotboom, “Iterative scheme for 1- and 2-dimensional d.c.-transistor simulation,” Electron. Lett. 5, 677–678 (1969).

054502-13 10

C. Kirsch and S. Mitran

L. Onsager, “Deviations from ohm’s law in weak electrolytes,” J. Chem. Phys. 2, 599–615 (1934). 11 C. L. Braun, “Electric field assisted dissociation of charge transfer states as a mechanism of photocarrier production,” J. Chem. Phys. 80, 4157– 4161 (1984). 12 M. C. Beard, K. P. Knutsen, P. Yu, J. M. Luther, Q. Song, W. K. Metzger, R. J. Ellingson, and A. J. Nozik, “Multiple exciton generation in colloidal silicon nanocrystals,” Nano Lett. 7, 2506–2512 (2007). 13 J. R. Tumbleston, personal communication (2011). 14 J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light (Princeton University Press, 2008). 15 € F. Bloch, “Uber die quantenmechanik der elektronen in kristallgittern,” Z. Phys. A 52, 555–600 (1929). 16 D. Y. K. Ko and J. C. Inkson, “Matrix method for tunneling in heterostructures: Resonant tunneling in multilayer systems,” Phys. Rev. B 38, 9945– 9951 (1988). 17 F. Abele`s, “La th#eorie g#en#erale des couches minces,” J. Phys. Radium 11, 307–309 (1950). 18 American Society for Testing and Materials, “Reference solar spectral irradiance: Air mass 1.5,” (2003), http://rredc.nrel.gov/solar/spectra/am1.5. 19 C. J. Brabec, N. S. Sariciftci, and J. C. Hummelen, “Plastic solar cells,” Adv. Funct. Mater. 11, 15–26 (2001). 20 € A. Einstein, “Uber die von der molekularkinetischen theorie der w€arme geforderte bewegung von in ruhenden fl€ ussigkeiten suspendierten teilchen,” Ann. Phys. 322, 549–560 (1905). 21 K. Hess, Advanced Theory of Semiconductor Devices (Prentice-Hall, 1988). 22 S. Selberherr, Analysis and Simulation of Semiconductor Devices (Springer, 1984). 23 J. W. Jerome, “Consistency of semiconductor modeling: An existence/stability analysis for the stationary van roosbroeck system,” SIAM J. Appl. Math. 45, 565–590 (1985). 24 J. W. Jerome, Analysis of Charge Transport (Springer, 1996). 25 N. B. Abdallah, F. M#ehats, and N. Vauchelet, “A note on the long time behavior for the drift-diffusion-poisson system,” C. R. Acad. Sci., Paris Ser. I 339, 683–688 (2004). 26 A. Arnold, P. Markowich, and G. Toscani, “On large time asymptotics for drift-diffusion-poisson systems,” Transp. Theory Stat. Phys. 29, 571–581 (2000). 27 J. K. J. van Duren, X. Yang, J. Loos, C. W. T. Bulle-Lieuwma, A. B. Sieval, J. C. Hummelen, and R. A. J. Janssen, “Relating the morphology of poly(p-phenylene vinylene)/methanofullerene blends to solar-cell performance,” Adv. Funct. Mater. 14, 425–434 (2004). 28 P. Langevin, “Recombinaison et mobilit#es des ions dans les gaz,” Ann. Chim. Phys. 28, 433–530 (1903). 29 L. J. A. Koster, V. D. Mihailetchi, and P. W. M. Blom, “Bimolecular recombination in polymer/fullerene bulk heterojunction solar cells,” Appl. Phys. Lett. 88, 052104 (2006).

J. Appl. Phys. 112, 054502 (2012) 30

M. M. Mandoc, L. J. A. Koster, and P. W. M. Blom, “Optimum charge carrier mobility in organic solar cells,” Appl. Phys. Lett. 90, 133504 (2007). 31 V. D. Mihailetchi, L. J. A. Koster, J. Hummelen, and P. W. M. Blom, “Photocurrent generation in polymer-fullerene bulk heterojunctions,” Phys. Rev. Lett. 93, 216601 (2004). 32 J. Nelson, The Physics of Solar Cells (Imperial College Press, 2003). 33 J. Bisquert and G. Garcia-Belmonte, “On voltage, photovoltage, and photocurrent in bulk heterojunction organic solar cells,” J. Phys. Chem. Lett. 2, 1950–1964 (2011). 34 E. H. Rhoderick, Metal-Semiconductor Contacts (Oxford University Press, 1978). 35 Y. M. Nam, J. Huh, and W. H. Jo, “Optimization of thickness and morphology of active layer for high performance of bulk-heterojunction organic solar cells,” Sol. Energy Mater. Sol. Cells 94, 1118–1124 (2010). 36 H. Jin, J. Olkkonen, M. Tuomikoski, P. Kopola, A. Maaninen, and J. Hast, “Thickness dependence and solution-degradation effect in poly (3-hexylthiophene):phenyl-c61-butyric acid methyl ester based solar cells,” Sol. Energy Mater. Sol. Cells 94, 465–470 (2010). 37 B. V. Andersson, D. M. Huang, A. J. Moul#e, and O. Ingan€as, “An optical spacer is no panacea for light collection in organic solar cells,” Appl. Phys. Lett. 94, 043302 (2009). 38 G. A. Buxton and N. Clarke, “Predicting structure and property relations in polymeric photovoltaic devices,” Phys. Rev. B 74, 085207 (2006). 39 E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen, LAPACK Users’ Guide, 3rd ed. (Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999). 40 H. K. Gummel, “A self-consistent iterative scheme for one-dimensional steady state transistor calculations,” IEEE Trans. Electron Devices 11, 455–465 (1964). 41 D. L. Scharfetter and H. K. Gummel, “Large-signal analysis of a silicon read diode oscillator,” IEEE Trans. Electron Devices 16, 64–77 (1969). 42 M. R. Hestenes and E. Stiefel, “Methods of conjugate gradients for solving linear systems,” J. Res. Natl. Bur. Stand. 49, 409–436 (1952). 43 M. J. Grote and T. Huckle, “Parallel preconditioning with sparse approximate inverses,” SIAM J. Sci. Comput. 18, 838–853 (1997). 44 J. Waldvogel, “Towards a general error theory of the trapezoidal rule,” in Approximation and Computation, Springer Optimization and Its Applications Vol. 42, edited by W. Gautschi, G. Mastroianni, and T. M. Rassias (Springer, 2011), Part 3, pp. 267–282. 45 J. Kiefer, “Sequential minimax search for a maximum,” Proc. Am. Math. Soc. 4, 502–506 (1953). 46 N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” J. Chem. Phys. 21, 1087–1092 (1953). 47 W. K. Hastings, “Monte carlo sampling methods using markov chains and their applications,” Biometrika 57, 97–109 (1970). 48 L. Ingber, “Very fast simulated re-annealing,” Math. Comput. Modell. 12, 967–973 (1989).