A Drift-Diffusion-Reaction Model for Excitonic Photovoltaic Bilayers: Asymptotic Analysis and a 2-D Hdg Finite-Element Scheme


We present and discuss a mathematical model for the operation of bilayer organic photovoltaic devices. Our model couples drift-diffusion-recombination equations for the charge carriers (specifically, electrons and holes) with a reaction-diffusion equation for the excitons/polaron pairs and Poisson’s equation for the self-consistent electrostatic potential. The material difference (i.e. the HOMO/LUMO gap) of the two organic substrates forming the bilayer device are included as a work-function potential.

Firstly, we perform an asymptotic analysis of the scaled one-dimensional stationary state system i) with focus on the dynamics on the interface and ii) with the goal of simplifying the bulk dynamics away for the interface. Secondly, we present a two-dimensional Hybrid Discontinuous Galerkin Finite Element numerical scheme which is very well suited to resolve i) the material changes ii) the resulting strong variation over the interface and iii) the necessary upwinding in the discretization of drift-diffusion equations. Finally, we compare the numerical results with the approximating asymptotics.

Photovoltaics; drift-diffusion-reaction equations; finite element methods; asymptotic methods.

AMS Subject Classification: 35K57, 35A35, 65N30

1 Introduction

The search for cheap, environmentally sustainable energy solutions has lead to intense interest and research in the area of organic photovoltaic (OPV) devices. Currently, these devices have solar efficiencies which are significantly below those of modern inorganic semiconductor devices (i.e. at present 8.3% compared to 27.6% for “one-sun” Gallium Arsenide and 42.3% for concentrator cells[9]), currently limiting commercial implementation.

A main difference of OPV devices (in contrast to inorganic PV) is the generation mechanism of free charge carriers, which typically occurs through the generation and dissociation of so called excitons, excited energy states created by incoming light. (We shall discuss excitons and the generation of free charge carriers in more detail below.)

Many mathematical models have been proposed to study the behavior of bilayer OPV cells as simple implementable organic devices: Refs. \refciteBC,Cetal,BRG,K. These models involve various approximations of the generation of free charges. A mathematical basis for the these models comes from extensive literature on inorganic semiconductor models. For these devices, the standard macroscopic models couple two drift-diffusion-reaction equations for the behavior of the free charge carrier densities (i.e. electrons and holes) in a system with the Poisson equation governing the self-consistent electrical potential.[15] Much is known about such systems, including existence and uniqueness results for the steady-state (see e.g. Refs. \refciteM,Mss and the references therein). The steady state is of particular interest for photovoltaic devices which are expected to generate power on time scales much longer than the duration of the transient (i.e. the switching-on) dynamics.

Drift-diffusion-reaction type models for OPV devices must crucially take into account the specific material properties of inorganic and organic semiconductor materials, which involves in particular electric field dependent mobilities and specific recombination and dissociation rates (see Sec. 2 below). A second crucial amendment has to model the exciton dynamics and couple it appropriately to the original system.

We shall thus consider the following evolution equations for the four main components, i.e. the concentrations of electrons (), holes (), and excitons () and the electric potential in the device (), which must be calculated self-consistently.

The evolution of the free charge carrier densities and shall be described by typical drift-diffusion-recombination equations and an additional source term due to the excitons:


Here denotes the positive fundamental charge, , and , are the diffusion coefficients and mobilities of and , respectively, is the recombination rate of free electrons and free holes, and is the rate of dissociation of excitons into a pair consisting of a free electron and hole. Functional expressions for diffusion coefficients, mobilities and recombination/dissociation rates in organic semiconductor materials will be specified in Sec. 2.

For a bilayer OPV device formed of two different organic semiconductor materials (see Fig. 1), the potential (which is not an electric potential and doesn’t contribute to the electric field) models an additional convection term that comes from the differences in the energy levels of the two materials. In an organic device, the energy levels involved in charge transfer are the Highest Occupied Molecular Orbital (HOMO) and the Lowest Unoccupied Molecular Orbital (LUMO). The HOMO roughly corresponds to the valence band for classical semiconductors, and the LUMO resembles the conduction band. In general, the HOMO-LUMO gap corresponds to the band-gap.

For most organic materials, the HOMO-LUMO gap is too large for a photon to create free electrons and holes. Instead, it creates a bound electron/hole pair, a so-called exciton. In the bulk material, these excitons usually recombine without producing free carriers (with rate , see below). However, at an interface between two materials with suitable differences in the HOMO/LUMO properties, the excitons tend to align and split over the interface so that the electrons and holes are in separate (energetically favorable) materials. This effect is included in the model equations (1) and (2) by the term , or more precisely, by the change in over the interface region.

Excitons are typically called polaron pairs (also referred to as charge-transfer states or coulomb bound pairs) when aligned across the interface and are much more stable with respect to recombination. On the other hand, the aligned polaron pairs will either be pushed together or pulled apart by an approximately parallel electric field. Accordingly the dissociation rate () will be heavily field dependent. The combination of these two effects (reduced polaron pair recombination and field-driven dissociation) can lead to very high quantum efficiency (proportion of light creating free charges) for polaron pairs at an interface under an appropriately applied external potential.[1] Note that this is an internal quantum efficiency for polaron pair dissociation, not an external quantum efficiency, which further includes how many of the incoming photons actually create polaron pairs.

Following the above discussion we shall assume that every exciton becomes immediately a polaron pair at the interface (see further explanation below). We therefore consider the following diffusion equation for the evolution of a combined density of excitons and polaron pairs (to be identified by their position), which lacks a convection term since the excitons are electrically neutral:


The coefficient is the photo-generation rate for excitons, is the rate of geminate recombination of the excitons, and represents the proportion of recombining free charges that form excitonic states (as opposed to recombination leading to emission of light, etc).

Finally, we have the Poisson equation for the self-consistent electric potential:


where denotes the vacuum permittivity and is the relative permittivity of the material. Furthermore, represents the derivative normal to the direction of the material interface and denotes the indicator function of the interface region while is the separation length of the charge pairs of polaron pairs in the interface region. In fact, the extra term proportional to is the field contribution due to the alignment of the polaron pairs at the material interface. The gradient results from taking the continuum limit of a sum of Dirac delta primes (corresponding to point dipoles). The exciton distribution includes both polaron pairs and excitons, with their identity based on their location in the system. Polaron-pairs are split over the interface such that the electron and hole (although still bound) each lie in an energetically favorable material. This has the dual effect of aligning a sheet of dipoles (resulting in the aforementioned term in Eq. (4)) and of making polaron pair states extremely energetically favorable to excitons.[1] We thus assume that any in an interface region are polaron pairs and that any (the complement of ) are excitons. In actual devices, organic interfaces are not sharp, and can be blended over a variety of length scales. For our model, we take to be the region within of the theoretical sharp interface. Note that if we take to be a straight line parallel to the contacts of the device and assume homogeneity in the parallel direction, we exactly recover the model in Ref. \refciteBRG.

Figure 1: Schematic of the 2-D device giving labels to the quantities used in the asymptotics and a graphical representation of the function . For the numerics we take to be a piecewise quadratic connecting the values at either side so that . Note that for the asymptotics the vertical direction of the first diagram is suppressed and the problem becomes 1-D.

We supplement the system (1) to (4) with the following boundary conditions. Dividing the boundary into disjoint Dirichlet and Neumann parts , we shall consider


The boundary conditions for in Eq. (5) correspond to the conducting ends of the device (the left and right of Fig. 1). The difference in the values for the potential at the two Dirichlet contacts correspond to the potential difference in the device (in Volts) whose offset does not affect the dynamics of the device because only the electric field enters the other equations. The usual notation for a photovoltaic device follows the convention scheme of a diode, but in the working case current flows in the direction opposite to that of a usual diode. For this reason, a consistent applied potential is given at the left boundary of the device. We thus take the potential at the right side of the device to be zero, and take the potential at the left side to be the potential difference. Note that this notation has the convenient consequence that the average field in the device has the same sign as the potential difference.

For the boundary conditions of electrons and holes one has to take into account that the HOMO/LUMO levels of the device determine different energies for the electrons and holes implying that the concentrations of and at the metallic contacts can be very different. Thus, for majority carriers in their energetically favorable material ( to the left of the interface and to the right as indicated in Fig. 1), we take the boundary conditions given by Scott and Malliaras in Ref. \refciteSM. For the minority carriers ( to the left and to the right), homogeneous Dirichlet boundary conditions constitute a very good approximation as the Scott-Malliaras boundary conditions give values approximately six orders of magnitude smaller than typical concentrations (see Table 1).

On the insulated boundary we take homogeneous Neumann conditions for all of the variables because no particles should be able to pass into the insulator and it should be electrically neutral. For the excitons, we take homogeneous Neumann conditions on the whole boundary because the excitons are not be able to transmit from the organic material into the inorganic contacts.

The results of this paper is as follows: In Sec. 2 we will scale the physical parameters of the equations and detail the specific expressions for the mobilities, dissociation rates and recombination rates in an organic semiconductor material. Next, in Sec. 4 we perform an asymptotic analysis of the 1-D equations for some typical parameters. The asymptotic analysis quantifies approximately how the exciton production rate yields an electric current and gives an expression for the shunt resistance of the device. Moreover, we present the numerical system used for the 2-D simulations in Sec. 5. In Sec. 6 we plot the calculated concentrations in the device at three important points in the device operation - short circuit, optimal power point, and open circuit. We interpret these results in terms of their relation to the asymptotic results and furthermore show the current-voltage characteristic curves generated by both the numerical and asymptotic methods.

2 Scaling and Models for Mobilities, Recombination, Dissociation

We will scale (1) - (4). The basic scaling introduces the following dimensionless quantities:


where is the thermal voltage, is a characteristic length (usually on the order of the device length), is a reference concentration, and is a reference mobility. We assume the Einstein’s relation . Although this might not be justified for some organic materials, it greatly simplifies the equations and as of yet a generally accepted alternative model has not been found. We further introduce the dimensionless quantities:


where denotes a reference time scale. Note that the Debye length is typically not a small parameter - see Table 1.

For an OPV device under light we follow Ref. \refciteBRG and choose the Langevin recombination rate , which rescales with the reference time as follows:

Setting this gives the following system:


For the charge carrier mobilities and we shall apply the Poole-Frenkel model[8] postulating an exponential dependence on the square-root of the electric field :


Moreover, the exciton dissociation rate shall be given by Ref. \refciteBRG as a function of , which is a scaled square-root of electric field :

where the , superscripts represent positive and negative fields (with respect to the interface normal).

Device Parameter values

The physical scaling parameters we take are:

The values for the various constants are taken from Ref. \refciteBRG (except which is converted to a spatial density instead of an area density). Table 1 collects all the used dimensionless parameters.


Values for physical parameters. \toprule \colrule \colrule \colrule \botrule {tabnote} , corresponding to the short circuit potential difference of . The expressions and refer to the rate-values in- and outside of the interface region.

The values of and are not easy to calculate physically and no consensus seems to exist in the literature for their values, but we take the given values based on reasonable estimates given in the sources listed above.

3 Steady-State Equations

For usual device operation, a solar cell will be producing current steadily for time-scales on the order of hours. Any transient behavior occurs over such short time-scales (i.e. microseconds[11]) that we neglect them for modeling the long-term behavior and efficiency of the device.

Therefore, for the remainder of the paper, we will consider only the steady-state equations where we drop all the tildes and absorb the time-scaling into the rates:


We expect that this system of equations will have a unique steady-state solution. If we insist on the physically reasonable requirement that and are smooth and bounded from above (corresponding to physical device saturation), then the first three equations fit very nearly into the standard semiconductor framework (see, for example, Ref. \refciteM). The only remaining difficulty is the exciton equation, but given as in the usual case, we see that for smooth we have and the exciton terms in the first three equations should not pose any difficulty. Proofs of the existence and uniqueness and the corresponding results for the time-dependent are under current investigation and subject to an upcoming paper.

Concerning the steady-state solutions, we emphasize that the light generation term implies a non-zero right-hand-side in the n and p equations of (13) and thus non-constant drift-diffusion fluxes of and on the left-hand-side. Thus, we expect a steady flux of electrons and holes away from the interface with nearly negligible bimolecular recombination (due to the work-function considerations discussed in the introduction). Even for , we will only recover constant stationary drift-diffusion fluxes in such particular cases such as and , which is not realistic given the considered parameters. As a consequence, the model system (13) is not designed to accurately predict the behavior of the device in the dark .

In the following section, we shall investigate further the steady state solutions of (13) for in a 1-D setting.

4 1-D Stationary State Asymptotics

4.1 Large Applied Field Approximation

We begin by considering the steady state equations in 1-D:


where we have assumed that the and are spatially homogeneous.

In reverse bias, which constitutes the operating state of an OPV bilayer device, a negative voltage is given at the left boundary of the device, i.e. , while zero voltage is given at the right boundary of the device, . Note that with this notation, is the potential difference in the device. For working photovoltaic cells, this potential difference is the sum of the potential applied to the device and a built-in potential from the metallic contacts.

With a constant relative permittivity we introduce the Debye length , and rescale the potential according to this potential difference, i.e. which rewrites the Poisson equation from (14) as


Since in a working device a typical short circuit voltage is about Volts,[3] which is in our units, we find that .

Hence (assuming for the moment that , where is the fraction of interface width to device width) we expect that the electric potential is in zeroth -order dominated by the potential difference and, thus that the electric field is constant in the zeroth order of , i.e.


We shall see in the following that the approximation (16) remains consistent with the assumption after being carried over to the equations for the charge carriers and the excitons .

However, first we quote from Ref. \refciteSM the typical order of magnitudes of the boundary values of electron- and hole- densities in the working state of the device:

Next, we remark that with the zeroth order approximation also the mobilities and are constant in zeroth. For the short circuit value , we calculate that and and thus . Thus, by rescaling


we obtain with the following rescaled charge carrier equations


Moreover, the rescaled equation for the excitons is


where is small, typically of order .

Neglecting the quadratic term, we can solve this equation explicitly in terms of hyperbolic sines and cosines. This allows us to check the consistency the approximation (16) and the rescaling (17) with the underlying assumption : With and outside the interface and with and (for typical ) inside the interface we have:


and thus if and .

Zeroth Order Charge Carriers

Now, we proceed to study the zeroth order solutions and of (18) and (19). Introducing the zeroth order fluxes

and neglecting the recombination term, the equations (18) and (19) integrate immediately as


Utilizing the Slotboom variables and , we can explicitly solve for and :


where and we define

Next, we can determine the parameters by evaluating the boundary conditions and . Upon rearrangement, we have with :


and obtain explicit formulas for and :


which satisfy both boundary conditions since .

These equations, although explicit, do not yet intuitively present the leading order contributions. However, using the fact that is constant outside of the interface, one can explicitly calculate the integrals and (see A).

In the following, we use these formulas to identify the leading contributions to and and thus and . More precisely, we introduce the parameters


and observe that for the considered device geometry since we have (see Fig 1). With and as in the working device, these parameters are small: and for with the numerical values given in Table 1.

However, using these parameters directly in (26) and (27) does not directly yield insights into the behavior of and . Because can change over many orders of magnitude, the behavior of and depends highly on competing exponential terms, none of which can be easily eliminated. One possible way to proceed splits the domain into three areas: left of the interface, the interface, and right of the interface. For each of these areas, we would obtain formulas for and of the form of (23) and (24) in terms of the values and and the boundary terms , and , . Then, imposing continuity of and and continuity of the fluxes and at and would allow us to determine the values and .

However, a simpler way to proceed considers the zeroth order currents, which are sufficient to understand the produced electric current. Using the explicit formulas for and in terms of and given in A we have from (25):


where denotes a mean-value of within the interface.

Because , the sum can be used to calculate the total current in the device (as predicted by the asymptotics). We shall plot and discuss the predicted relationship between the current and the applied field as asymptotic IV-curve in Sec. 6 (see Fig. 9). Note that for the given plots we take to be piecewise linear so that is explicitly defined.

Short-Circuit Current

As mentioned earlier, the sign of determines if the parameters are large or small. For the short circuit current, we have and . In order to determine the lowest order terms of these equations, we must also calculate a more precise form of and thus its order. At first we remark that (since on the interface with width and outside the interface with device length ). Thus, one can verify with calculations similar to those below and in A that the exponentials of will determine the leading contributions to the integrals . More precisely, since is strictly increasing in the regime (recall that ) the leading order contributions derive from

Applying the definitions of Eq. (28), integration by parts and the mean value theorem yields

for a mean value . Together with similar calculations for and , we obtain

Taking only the dominating terms in (29) and in (30) yields:

The term gives a current which is proportional to the field, indicating that it represents a resistor. For a usual working device this term will be small ( in Sec. 2). Because the term does not depend on the length of the device, we interpret it as a shunt resistance (parallel to the device). We observe from Eqs. (29) and (30) that we have another current term with the form of a shunt resistance: . This factor is negligible in the short circuit case, but becomes large as we move to positive fields, where . We investigate the total shunt resistance further in Sec. 6.

The second two terms represent the interaction of the excitons on the system. For usual device parameters, (for the parameters given in Sec. 2, we have ). Since on the interface and outside the interface, we can neglect the contributions from outside the interface (recall that and ) compared to the contribution from within the interface.

Thus the most important contribution to the current for negative fields is:


This approximation works very well in the short circuit case (as examined in Sec. 6, in the discussion preceding Fig. 4), where it replicates within 2% and is actually closer to the numerically calculated .

Note that comes entirely from the term. Because we have evaluated and at the point , we expect that the hole current will be dominant. However, using Eq. (22) we can calculate and , and that the only change from the currents at is that the term appears in the electron current. Thus, as expected, in the region where the electrons are favored, the primary contribution to the current comes from the electrons and the total current is conserved.

First Order Terms

Since we have explicit forms for , , and , we can integrate twice to calculate (with an additional linear term to account for the two boundary values of ). Normally we don’t plot since it is generally dominated by its boundary terms. However, can also be calculated explicitly (by a single integration) and we include the first-order field in the plots in Sec. 6. The explicit form is not very illuminating and thus not written here. Note that we need to add a constant value to the field to insure that its integral gives the correct potential difference in the device (corresponding to the slope of the aforementioned linear term).

In theory we could use the second order form for the potential to calculate the next order and solutions. In fact, this is more or less the essence of the iteration scheme outlined in Sec. 5. However, without a constant field, it is no longer possible to explicitly integrate the continuity equations (especially given the sub-exponential forms for the electron and hole mobilities). It would be possible to do these calculations numerically, but this would simply become a simplified version of our numerical iteration scheme, and thus we do not pursue this further in the asymptotic case. Instead, we will present another method for asymptotic calculations.

4.2 Unipolar approximation

The discussion of the previous section, in particular Eq. (31), can be summarized by the statement that the electron- and hole- fluxes and are approximately constant outside the interface but feature a strong variation over the interface with a magnitude of .

As a consequence, using the Eqs. (26) and (27) and similar calculations as in Sec. 4.1.2 and A, it follows that the zero-order electron- and hole-densities and vary also strongly over the interface. In fact, we obtain for the hole density :