1. Introduction

The rise and the decline of human civilizations has been intimately linked to the availability of large and reliable sources of energy. With our Sun being the most abundant source of energy in our solar system, life on Earth has most likely harnessed solar energy for billions of years already. Humanity itself has developed increasingly sophisticated methods to make use of solar energy over several millennia [1].

History also shows that the more human technologies have advanced, the more power hungry these technologies became. The simplest solution was to consume unsustainable but vast sources of energy provided by fossil fuels, which lead to severe environmental problems. This is not a new issue, and photovoltaic devices, which directly turn sunlight into electric power, have already been suggested as an alternative and environmentally benign energy source since the late 19^th century [1].

At the time of writing this article, the humble beginnings of photovoltaics in the 19^th century have developed into an established field at the forefront of materials science, which combines state-of-the-art device fabrication methods with the most sophisticated material analysis techniques, as well as with detailed device modelling approaches. As a result, the production costs for photovoltaic devices are dropping quite rapidly, while their efficiencies have shown substantial increases over the years, many of which are based on a much better understanding of the underlying light-matter interactions.

Due to the importance of fundamental light-matter interactions, the simulation of photovoltaic devices turns out to be a typical multi-scale problem, which must combine numerical studies of light-matter interactions at the atomic scale with a phenomenological description of huge numbers of photons impinging on a solar cell, which are generating huge numbers of free charge carriers that will ultimately give rise to a macroscopic electromotive force (emf). This multi-scale approach requires the knowledge of a variety of materials properties [2], which can either be taken from experiment, or they may as well be determined using first principles materials simulations. The latter will be the topic of this study.

This paper will focus on the internal processes of solar cells and are their simulation. We know that integration of these cells into external circuitry, panels and other overall experimental/industrial set-ups are equally important topics, but independent from the topic of this paper. Another important aspect, which is the overall life-cycle and degradation of solar cell devices, cannot be covered here as well. To our knowledge, no adequate simulation method exists for this important issue.

Before we embark on a more detailed description of ab initio based solar cell simulations, we will first discuss the generic features that are common to the different types of photovoltaic devices currently in use, and the various technological challenges that these devices must overcome.

A generic layout of a typical photovoltaic device is shown in Fig. 1(a), which refers to a p-n junction made of n-doped and p-doped bulk semiconductor materials. Region I is where the impinging photons create electron-hole pairs by a fundamental process that is shown in Fig. 1(b). If the incoming photons have an energy smaller than the band gap ${E_g}$ of the bulk semiconductor materials, then these photons will not be able to generate free charge carriers and simply heat up the solar cell, which leads to considerable losses.

 

Fig. 1. (a) Structure of a typical solar cell device: Region I is characterized by the generation of free charge carriers due to fundamental light-matter interactions; region II describes the migration of free charge carriers towards the contacts, and region III describes the extraction of free charge carriers at the contacts, which subsequently drive a current through the attached wires. (b) Major losses in a solar cell: photons with energies lower than the band gap ${E_{gap}}$ are not able to generate electron-hole pairs, while higher energetic photons generate electron-hole pairs and free charge carriers, which lose a lot of their energy due to thermalization.

Download Full Size | PDF

But even in the case where an electron-hole pair is created by an incoming photon, this process will always be accompanied by the time-reverse process of recombination and the re-emission of a photon. Thus, in order to generate a current, the layout of the solar cell and the basic properties of the bulk materials must provide a mechanism to separate the electrons from the holes, and to generate free charge carriers that are migrating towards the contacts.

In the case of a p-n junction this separation of free charge carriers is facilitated by a large electric field that is provided by fixed ions in the so-called depletion region, whose formation is caused by the diffusion of mobile majority charge carriers from the doped bulk materials. The existence of a depletion region is independent of the solar illumination of the device, and its size may be changed by an applied voltage [2,3]. Furthermore, within the doped bulk materials of the p-n junction, there are drastically different mobilities for electrons and holes. Electrons have high mobilities on the n-side and low mobilities on the p-side, and the opposite applies to holes. In fact, some authors even picture the resulting separation of charge carriers as some sort of membrane process [3].

After the charge carriers have been separated, they will pass through region II on their way to the contacts, which is marked by all kinds of fundamental losses. One problem is the trapping of free charge carriers, f.e. by impurities and dislocations [2]. But there are also losses due to basic inelastic interactions between the free charge carriers and the lattice, which lead to a thermalization of electrons and holes as illustrated in Fig. 1(b). Due to these thermalization processes, the charge carriers will effectively lose a good fraction of their energies to the ionic lattice, which reduces the current and simply heats up the solar cell.

In region III the free charge carriers have finally reached the contacts and must be extracted in the most efficient fashion. This is accompanied by various technological challenges. First of all, the charge carriers must overcome a potential barrier that is always formed at the interface between two different materials. Therefore, the probability for a charge carrier to be extracted at the contacts is intimately linked to the nature of the potential barrier. Furthermore, the interfaces between different materials are never perfect, but marked by lattice mismatches and all sorts of defects. This leads to surface trapping of free charge carriers [2].

After developing a general feeling for the various physical processes in a typical solar cell device, the rest of this article will be devoted to a much more detailed study of these fundamental processes and of the technological challenges involved. In Section 2 we will discuss simple model approaches to describe the fundamental physics that takes place in regions I-III. We will point out how these models can be made more realistic using model parameters based on ab initio data. In Section 3 we will describe the physics and the numerics behind a typical state-of-the-art solar cell device simulation, where we also discuss the standard set of model parameters that are required to carry out such simulations. And we indicate how they may be obtained using ab initio simulations.

Section 4 will be devoted to ab initio materials simulations. We will present some of the concepts and numerical approaches currently used to study fundamental light-matter interactions, and we will also describe in some detail how to determine various material parameters that feed into a typical solar cell device simulation. In Section 5 we will summarize our findings and briefly discuss the implementation of advanced light management features like frequency-conversion layers and plasmonic nanoparticles. This section will also contain a brief discussion of future developments in the field of ab initio methods and the new opportunities that these developments may offer for solar cell simulations.

However, before we embark on a lot of technical discussions, there are also some cautionary remarks to be made here. Although we will show that with the help of ab initio methods, one could in principle run a solar cell device simulation entirely on a computer and never see a lab, the reader should not get the false impression that the methods involved are some sort of magical wand to study solar cells. First, the numerical accuracy of ab initio methods is limited. Second, we will see that there are complex materials parameters and processes, which we will never be able to simulate in sufficient detail using ab initio methods.

On the other hand, it must also be pointed out that ab initio methods have the distinct advantage that they provide us with many idealized fundamental scenarios, against which one must critically gauge any “spectacular” experimental data, and hopefully get a clue about the underlying physics or about the possible sources of error that might have led to those “spectacular” results.

2. Model approaches

In this section, we will describe three simple models that will cover the essential physics in regions I-III of a typical working solar cell. A similar discussion can be found in [4]. Each of these models will impose certain physical limits for the practical efficiency of a given solar cell.

Let us start our discussion by assuming that a voltage V has been applied to a model photovoltaic device, which is generating a current density J. Then the power density P of that device will be $P = J \cdot V$. Given a typical solar irradiation ${I_{sol}} \approx \frac{{1000\; W}}{{{m^2}}}$ the efficiency $\eta $ (in %) of such a device is basically given by:

Let us now focus on region I, where free charge carriers are generated by a flux of incoming photons. A suitable phenomenological approach for region I may be based on the diode equation, which describes the current density-voltage (J-V) characteristics of a typical p-n junction (see Fig. 1(a)).

However, we will augment the diode equation by the addition of a light current density ${J_l}$, which is generated by the solar cell under solar irradiation:

Here e is the elementary charge, ${k_B}$ is the Boltzmann constant and T is the temperature of the device. The so-called saturation current density ${J_s}$ is about twelve orders of magnitude smaller than the light current density ${J_l}$, and it can be approximated by an educated guess within the right order of magnitude, like ${J_s} \approx {10^{ - 8}}\; A/{m^2}$. As illustrated in [5] the efficiencies of a model solar cell based on the diode equation will depend very little on the precise value of ${J_s}$.

The light current ${J_l}$ may be obtained from

Here d is the size of the depletion region, where most of the electron-hole pairs are generated. The precise value of d depends strongly on the doping of the bulk materials in the p-n junction, but it is usually of the order of ${10^{ - 4}}\; m$. The symbol $G({\omega ,z} )$ stands for the generation rate, which depends on the position z measured from the top of the p-n junction, which is shown in Fig. 1(a). The generation rate also depends on the angular frequency $\omega $ of the incoming photons. Let us assume that we dispose of a proper model for the incoming photon flux $\mathrm{\Phi }(\omega )$ (see [2] for photon flux models referring to solar irradiation). Then the generation rate follows from:

We see from this equation that the initial photon flux $\mathrm{\Phi }(\omega )$ decreases according to Beer’s law [2], where $\alpha (\omega )$ is the attenuation coefficient of the bulk material, which can be determined using ab initio methods. In Eq. (4) and for the rest of Sections 2 and 3, we will emphasize in boldface all model parameters that can be determined using ab initio methods. Most of these model parameters are discussed in more detail in Section 4 below. Based on a suitable (average) $\alpha (\omega )$ (f.e. as taken from ab initio simulations) we can then embark on the simulation of J-V characteristics using the diode equation, and we will obtain results that look like the curve shown in Fig. 2(a). Note that the maximum power density provided by the device is located at the maximum power point ${P_m}$ [2].

 

Fig. 2. (a) J-V characteristics of a solar cell: The blue curve describes the “dark” characteristics of a solar cell, whereas the lower red curve describes its characteristics under illumination. ${P_m}$ marks the maximum working point. (b) Model for the trapping of free negative charge carriers by a fixed positive charge. The parameter ${a_\ast }$ characterizes the corresponding cross section. (c) Model contact between ideal metal and an n-type semiconductor (Schottky barrier): We observe band bending and the occurrence of a barrier of height ${\phi _B}$. $\Delta \phi $ denotes the difference in work function between the metal and the n-type semiconductor.

Download Full Size | PDF

The generation rates predicted by Eq. (4), and thus the efficiencies of the corresponding model solar cells, are usually very high. To describe more realistic systems, we must correct the ideal generation rate G, because some charge carriers are always lost due to recombination processes. Corrections of the ideal generation rates are usually described by a recombination rate R, and then the latter must be subtracted from the generation rate G. One of the simplest models is based on the assumption that the recombination rate is related to the excess charge carrier density $\Delta n$ and to the carrier lifetimes $\tau $ via [2]:

Note that we have shown the carrier lifetimes in boldface, because they can at least be partially determined using ab initio methods, which was shown in [6].

Now we come to region II. To model the physics in this region, we invoke Beer’s law again, and assume that the current densities J will decay exponentially away from the depletion region, where all these free charge carriers were originally created. The decay of the current density is caused by the presence of traps, which could be impurities or dislocations. Let ${z_0}$ mark the location of the depletion region, then the current densities will decay away from the depletion region according to:

The absorption coefficient for these processes is given by $\alpha^{\prime} = {N_{trap}} \cdot {\boldsymbol{\sigma }_{\boldsymbol{trap}}}$, where ${N_{trap}}$ is the density of traps, and ${\sigma _{trap}}$ is the cross section for the trapping process. A simple model for such a trap is shown in Fig. 2(b), where a fixed positive charge catches an electron. As discussed in [4], with a typical density of traps ${N_{trap}}\; \approx {10^{20}}\; {m^{ - 3}}$ for dirty Si and assuming a simple exciton-like model for the process in Fig. 2(b), we should expect cross sections of the order of ${\sigma _{trap}} \approx {10^{ - 18}}\; {m^2}$. Such exciton-like models and the corresponding cross sections for carrier trapping can in principle be parameterized using ab initio methods, based on the knowledge of the ab initio dielectric function and of the effective masses of the charge carriers extracted from ab initio band structures, see [6]. Further details concerning the modelling of excitons will be given in Section 4.

Finally, to model free carrier extraction at the contacts in region III, we will entirely focus on the transmission probability T, which relates the incoming current density ${J_{inc}}$ to the transmitted current density ${J_{trans}}$. The transmission probability is roughly given by the following expression:

Thus, the transmission probably through the barrier is obviously determined by yet another absorption coefficient $\alpha^{\prime\prime}$ and by the barrier width w. For the barrier shown in Fig. 2(c), the corresponding barrier width may be approximated by:

We see that both, the absorption coefficient $\alpha^{\prime\prime}$ and the barrier width w depend on the barrier height ${\phi _B}$, which is illustrated in Fig. 2(c). We also see from that picture, that a more accurate way to estimate the barrier height involves the difference in work functions $\Delta \phi $ between the metal and the semiconductor, and we will study this scenario in Section 4 using ab initio methods.

According to Eq. (8), the barrier width will depend on the donor concentration ${N_D}$, the vacuum permittivity ${\varepsilon _0}$ and on the dielectric function ${\boldsymbol{\varepsilon }_{\boldsymbol{r}}}$. In Section 4 we will show that the latter can be determined quite accurately using ab initio methods. In a recent study [4] we have also estimated some of the typical data that characterizes the Schottky barrier shown in Fig. 2(c). Assuming n-doped Si with a donor concentration ${N_D} \approx {10^{26}}\; {m^{ - 3}}$ and a typical dielectric function of ${\varepsilon _r} \approx 10$, we saw that ${\phi _B} \approx 0.7\; eV,\; \; w \approx {10^{- 9}}\; m$, and $\alpha^{\prime\prime} \approx {10^9}\; {m^{ - 1}}$. According to these numbers and by referring to Eq. (7), we may guess that a high transmission probability requires a very fine balance between a huge absorption coefficient $\alpha^{\prime\prime}$ and a tiny barrier width w. Any mismatch will lead to considerable losses, which makes carrier extraction at the contacts a very delicate process.

3. Device simulations

The idea of state-of-the-art solar cell device simulations is to model all the optical and electronic processes, which take place in a photovoltaic device. These simulations will produce J-V characteristics that look like the schematic curves shown in Fig. 2(a). But beyond that, a typical device simulation will also provide us with detailed information about the flow of photons and charge carriers through a model solar cell.

Device simulations may also serve as a basis for systematic studies that explore the importance of various design parameters [2]. Among the most important of these design parameters are the sequence and the size of the different layers of a solar cell and their dielectric and electronic properties, doping levels and trapping densities. Some of these design parameters must be set or estimated by the designer or by the experimentalist fabricating such a device. Other key materials parameters may as well be taken from ab initio simulations of the underlying bulk materials.

Before we describe the main features of a typical solar cell device simulation, there are some general issues that must be pointed out right from the beginning. First, we are now in a situation where we are no longer dealing with single electrons or few electron systems, but with macroscopic carrier densities, which are best described by some sort of electron and hole gases. Interactions of these electron gases with incoming light is generating macroscopic currents, which are flowing through the solar cell and drive the system out of thermal equilibrium. We will model this situation in the following way: despite the fact that the system is driven out of equilibrium, the electrons themselves are still at equilibrium with each other. We describe this quasi-equilibrium by a quasi-Fermi level that is traditionally denoted by the symbol ${E_{Fn}}$. The holes are also at equilibrium with each other, which is denoted by a quasi-Fermi level ${E_{Fp}}$. However, the corresponding charge carrier gases are no longer at equilibrium with each other, and therefore we find that ${E_{Fn}} \ne {E_{Fp}}$.

Furthermore, for a typical model device like the one shown in Fig. 1(a), we will make another drastic simplification and assume that all physical properties vary only in the vertical $z$-direction, which makes these types of device simulations quasi-one dimensional. This implies that the quasi-Fermi levels vary as ${E_{Fn}}(z ),\; {E_{Fp}}(z )$, and we will take both quasi-Fermi levels will as basic variables for our device simulations. Local values of the quasi-Fermi levels will determine the local free carrier densities $n(z ),\; p(z )$.

Another basic variable for our simulations will be the electrostatic potential $\varphi (z )$, from which we can determine the overall cell potential V. This means that we are now dealing with a set of three basic variables ${E_{Fn}}(z ),\; {E_{Fp}}(z )$ and $\varphi (z )$. The potential $\varphi (z )$ may be determined from the following Poisson equation:

The donor and acceptor densities ${N_D}(z ),\; {N_A}(z )$, as well as the trapped carrier densities ${p_t}(z ),\; {n_t}(z )$ are dependent on macroscopic features of the bulk materials, and they also depend on the basic design of a given device. The dielectric functions $\varepsilon (z )$ can either be taken from experiment, or it may be calculated using ab initio methods, as discussed in Section 4.

The electron and hole current densities ${J_n}(z ),\; {J_p}(z )$ may be obtained by solving two related continuity equations:

The symbols $G(z )$ and $R(z )$ denote the $z$-dependent generation and recombination rates as discussed in the previous section. In a typical device simulation, the generation rates $G(z )$ used in Eq. (10) involve a suitable averaging over the absorbed frequencies, the details of which are described in [2] and other textbooks. Furthermore, note that the generation rates are shown in boldface. They are obviously related to the frequency dependent absorption coefficient according to Eq. (4), and we have already emphasized that the absorption coefficient can be determined using ab initio methods.

In a standard type of device simulation, the implementation of recombination rates $R(z )$ is based on the Shockley-Read mechanism [7], whose linearized version is given by Eq. (5). Alternative approaches using ab initio methods would be extremely helpful in this regard, and therefore we have tentatively shown $R(z )$ in boldface as well.

The current densities ${J_n}$ and ${J_p}$ in Eq. (10) are directly related to the quasi-Fermi levels ${E_{Fn}}(z ),\; {E_{Fp}}(z )$ by the following expressions:

Here, ${\mu _n},\; {\mu _p}$ denote electron and hole mobilities, which can also be determined using ab initio calculations.

From a numerical point of view, the system of Eqs. (9)–(11) must be solved for each point on a grid in the vertical $z$-direction. As we are dealing with ordinary differential equations, it will be necessary to specify some proper boundary conditions. But according to Fig. 1(a), the boundaries of the electronic part of a typical device are the contacts, which means that any specification of boundary conditions will also imply the modelling of carrier extraction at these contacts.

The specification of these boundary conditions turns out to be a very challenging task, and naïve modelling often leads to the prediction of efficiencies, which are drastically reduced in comparison to experimental results or predictions made using the simple model approaches discussed in Section 2. Further details about the suite of models used to simulate carrier extraction can be found in the corresponding literature [2]. Many of these models typically depend on the electron affinity of the bulk materials, the latter being the energetic difference between the vacuum energy level of escaping free electrons and the bottom of the conduction band. Electron affinities can be determined from first principles, but we guess that alternative modelling of carrier extraction based on systematic ab initio simulations of bulk-contact interfaces might be a viable alternative in this regard. An example will be shown in Section 4.

Let us finally have a look at a typical device simulation, which is shown in Fig. 3. Here we have simulated the current density vs. voltage curve for a perovskite solar cell using the GPVDM simulation code [8], based on the provided default parameters. To exemplify a typical modelling application, we have varied the perovskite band gap from its default value of 1.6 eV. We see from Fig. 3(a) that the band gap is near optimal, and only very small improvements can be achieved by slightly lowering the band gap.

 

Fig. 3. Current density vs. voltage curve for a perovskite solar cell simulated with GPVDM a) using default values with the band gap being varied around its default value of E_g = 1.6 eV b) using ab initio values for the optical parameters (see text).

Download Full Size | PDF

To show the applicability of ab initio data, we have also replaced some default simulation parameters with values computed from ab initio simulations. Among the many device simulation codes that we have tested over the years, the GPVDM software allows for one of the easiest and most flexible input options concerning material-specific simulation parameters that are taken from ab initio simulations. Other programs often require different and less flexible data formats, and they may also use slightly different models for charge carrier trapping or carrier extractions at the contacts, which makes quantitative comparisons between different simulation programs often very difficult. Therefore, the main focus of device simulations should always be on systematic studies using the same simulation package.

For our present study, we took ab initio data into account by changing the static permittivity from 20 to 27.8 to include the phonon contributions and we have also changed the band gap from 1.6 eV to 1.44 eV. Finally, we have substituted the build-in refractive index and absorption coefficient with calculated values using the methods described in Section 4. We see from Fig. 3(b) that changes in the static dielectric constant have virtually no influence on the device performance, which makes it unnecessary to obtain very precise values here. However, the exact band gap and the absorption spectrum are of great importance, and they must be determined with the highest possible precision.

4. Ab initio data sets

After the discussion of two different modelling approaches and their most important parameters in Sections 2 and 3, we now describe how some of these modelling parameters can be taken from ab initio simulations. This should not come as a big surprise. In fact, decades of research in computational materials science have shown that quantum mechanical simulations can provide us with a plethora of valuable insight into material properties. In many cases ab initio modelling can already alleviate the necessity for synthesis and characterization in the Lab, and these methods yield a complementary set of materials data, thanks to the direct access to the electron wave function, which allows for a fine-grained control of the investigated system.

The usefulness of ab initio simulations is generally limited by simplifying assumptions to deal with the complexity of the many-body problem. For example, ab initio simulations most commonly use ideal crystalline structures and ignore temperature effects. Defects, impurities and grain boundaries are difficult to investigate, as they require fairly large simulation cells, which contain many atoms and even more electrons. Improvements in computing technologies are constantly pushing the limits of what is still doable, but clear numerical limitations remain.

Most ab initio simulations of solids are based on Density Functional Theory (DFT). DFT was developed in the late 1960’s by Hohenberg, Kohn and Sham [9,10]. It simplifies the many-body problem by replacing the 3N electron coordinates in the Schrödinger equation by an electron density $n(\boldsymbol{r} )$, and mapping a system of real interacting electrons onto an equivalent system of non-interacting pseudo-electrons [11], which is described by the following energy functional:

The electron density $n(\boldsymbol{r} )= \mathop \sum \nolimits_\sigma \mathop \sum \nolimits_i {|{\psi_i^\sigma (\boldsymbol{r} )} |^2}$ is given by the sum over all (occupied) orbitals and over all spins of the non-interacting reference system. This is the essence of the Kohn-Sham approach (KS-DFT). The independent particle kinetic energy ${T_S}[n ]={-} \frac{1}{2}\mathop \sum \nolimits_\sigma \mathop \sum \nolimits_i \psi _i^\sigma \textrm{|}{\nabla ^2}\textrm{|}\psi _i^\sigma $ is similarly given by the expectation value of a Laplacian related to these orbitals. The external potential ${V_{ext}}(\boldsymbol{r} )$ includes contributions from the ionic potentials and other external fields. The Hartree energy ${E_H}[n ]= \frac{1}{2}\smallint \frac{{n(\boldsymbol{r} )n({\boldsymbol{r^{\prime}}} )}}{{|{\boldsymbol{r} - \boldsymbol{r^{\prime}}} |}}{d^3}r\; {d^3}r^{\prime}$ contains the mean-field electron-electron Coulomb repulsions marked by an unphysical self-interaction. The necessary self-interaction corrections and all many-body effects are conveniently pushed into the exchange-correlation term ${E_{xc}}[n ]$. Finally, ${E_{II}}$ describes the repulsions between the various ions in the system.

The setup of DFT is ideally suited for computers. In principle it is even exact, provided one would know the exact exchange-correlation functional ${E_{xc}}[n ]$, which yields the exact total energy and electron density $n(\boldsymbol{r} )$ for a given system. A large variety of approximate exchange-correlation functionals have been developed over the decades, but an exact expression has not been found yet. While it is easy to show that there exists an exact exchange-correlation functional for every electronic system, to our knowledge it has never been demonstrated in practice that there is indeed something like a universal functional, which would provide the same accuracy for all many-electron systems. Despite these conceptional problems, modern DFT calculations are fairly reliable, and a lot is known about their practical capabilities and limitations. For a more detailed insight into DFT we refer to the books of Martin [11] and Kohanoff [12].

By construction, KS-DFT and other DFT approaches are only guaranteed to yield good ground-state electron densities and total energies. Unoccupied bands do not contribute to the electron density and are therefore not properly optimized. As a result, KS-DFT is generally not very reliable for any excited states properties, which leads to a tendency of KS-DFT to vastly underestimate the band gap of materials [13].

Amongst the ground-state properties that are most often investigated by KS-DFT are elastic properties (bulk modulus, hardness), phonons, and surface properties like work functions and adsorption energies. In the context of solar cell materials, the most important surfaces are the interfaces between bulk materials and the contacts, where many losses can occur. The work function $\phi $ for such a surface is defined (numerically) as the difference between the vacuum level ${V_{vac}}$ and the Fermi energy ${E_F}$ of the bulk material [14]:

The Fermi energy depends on the doping level of the bulk materials, a fact that is nicely illustrated in [15,16]. Differences $\Delta \phi $ in the work functions across an interface may serve as an approximation for the contact barrier height between the two adjacent materials [17]. As an example, we have studied the barrier heights for an Au-Si interface, using fcc $<111>$ surfaces for both bulk materials. Figure 4 shows the average electrostatic potential of thin Si (left) and gold (right) slabs as calculated with the GPAW code [18,19] and using the PBE exchange-correlation functional [20]. As with all other simulations in this article, relevant input files and results are available in Dataset 1 (see Ref. [21]). The valleys in the potential occur at the position of the atoms. On the left and on the right-hand side of each slab there is vacuum. The corresponding Fermi levels are represented by horizontal yellow lines. Work functions of about 4.6 eV for Si and 5.2 eV for Au do agree well with the experimental values, but great care must be taken to check that calculations are converged, especially with regards to the number of atomic layers in the slab [14]. Based on this information, we would estimate the Au-Si contact barrier height to be around 0.6 eV, which is a little higher than the experimental value of about 0.5 eV [17].

 

Fig. 4. Work function $\phi $ of a Si $<111>$ slab (a) and Au $<111>$ slab (b) as calculated with DFT and the PBE xc-correlation functional using the GPAW code in the finite difference mode (h=0.18, k=(9,9,1)). The orange line in the Fermi level and the blue curve is the plane-averaged electrostatic potential of the valence electrons.

Download Full Size | PDF

When moving away from pure ground-state properties of materials, there is no good reason to put too much faith into KS-DFT results. This clearly affects fundamental materials properties like the electronic band-structure, which comprise both ground state (valence bands) and excited state (conduction bands) information. But despite these formal limitations, KS-DFT has proven to be a surprisingly valuable tool for determining the electronic band structure and the density of states of materials.

The electronic band structure contains information about the position and size of direct and indirect band gaps, and about the symmetry of the bands. Total and projected density of states can give further inside into the character of the bands with respect to atomic reference states. For solar cell materials it is also interesting to extract the curvature of the valence band maximum and of the conduction band minimum, which determine the effective electron and hole masses. They can be also used to determine the electron and hole mobilities (see [5] for details).

Figure 5 shows the band structures of crystalline silicon and a CH₃NH₃PbI₃ perovskite. Our calculations were performed using Quantum Espresso [22] with the PBEsol xc-functional [23] and the SG15 norm-conserving pseudo-potentials [24]. The band structure of fcc Si shows four valence bands as well as the lowest four conduction bands. For the sake of comparison, results from the more accurate GW method are also included, which will be discussed in more detail below. The main message here is that KS-DFT does predict the band dispersions astonishingly well, but the biggest shortcoming is the underestimation of the band gap. Our calculations yield an indirect band gap of 0.34 eV, which is considerably smaller than the experimental value of about 1.1 eV, and which illustrates quite nicely the systematic underestimation of band gaps by KS-DFT.

 

Fig. 5. Band structures of Si (a) and the CH₃NH₃PbI₃ perovskite (b) as calculated with KS-DFT and the PBEsol xc-correlation functional using the Quantum Espresso code (blue) and as calculated with the G₀W₀ method using the Yambo code. Bands are aligned to the valence band maximum.

Download Full Size | PDF

The perovskite band structure in Fig. 5 comprises the six highest valence bands and the four lowest conduction bands. Here, spin-orbit coupling must be taken into account, due to the presence of heavy ions in such materials, and the corresponding modifications of the SG15 potentials were used [25]. One can see quite nicely the band splitting effect of the spin-orbit coupling, which moves the conduction band minimum slightly away from the A point. The CH₃NH₃ molecule at the center of the unit cell introduces a slight anisotropy, which changes the shape of the cell from cubic to orthorhombic. This minute anisotropy is hardly visible when comparing the X and Z directions in the band structure. Quite obviously, the CH₃NH₃ molecule does not really partake in the electronic or optical properties of the perovskite. Its sole purpose consists of providing strain to the unit cell.

Now that we have seen the limitations of KS-DFT calculations, the question is how one may overcome some these difficulties. One elegant solution is to use Hedin’s GW method [26]. In the technical framework of this method, the electronic many-body problem can be written in terms of a single particle Green’s function G [12]:

The self-energy operator Σ is non-local, non-Hermitian and energy-dependent, and it contains the static and dynamic exchange-correlation effects. This operator should correct for the short-comings of KS-DFT by including self-interaction corrections and by correcting the lack of a discontinuity in the exchange-correlation potential when adding electrons [13]. The latter is the main cause for the underestimated band gaps.

The Green’s function of Eq. (14) can be computed using a Dyson equation, which may be derived from many-body perturbation theory. In practice however, the exact perturbative formalism is far too cumbersome, and instead one is falling back on the Green’s function for the non-interacting electron system G₀ [12]:

Here we use the eigenstates ${\psi _n}(\boldsymbol{r} )$ and eigenvalues ${\epsilon _n}$ from KS-DFT calculations as input. The self-energy Σ may be calculated systematically with the help of some higher order two-particle Green’s functions, which is a very elegant approach, but again not feasible in practice. Instead, the Hedin approximation is used, which defines the self-energy as the product of the single particle Green’s function G and the dynamically screened Coulomb interaction W [12]:

This combination gives the GW method its name. The screened Coulomb interaction W is given by a Dyson equation that involves the polarization, which in turn is computed using the Green’s function. Thus, the whole technical procedure becomes quite involved, and practical numerical solutions are comprehensively described in a review article by Onida et al. [27].

The GW method itself is very resource intensive from a computational point of view, and therefore most of the standard numerical codes limit themselves to the (single-shot) G₀W₀ approximation, which uses the KS-DFT directly and simply skips the self-consistent solution of the Dyson equations. This approach usually works quite well in practice, but unfortunately, as with KS-DFT, there is no known analytical way to estimate the error for a given solution of a given many-electron system. Concerning the Si example from above, G₀W₀ as implemented in the Yambo code [28] predicts an indirect band gap of 0.85 eV, which is much better than the KS-DFT value, but still below the experimental value. For the CH₃NH₃PbI₃ perovskite the overall picture is even more drastic, where the KS-DFT band gap of 0.43 eV is vastly smaller than the experimental gap of about 1.55 eV [29]. G₀W₀ predicts a gap of 1.44 eV, which is fairly close to the experimental value.

In the framework of the GW scheme, we must consider electrons and holes to be quasi-particles (QP), and thus account for perturbations due to other quasi-particles being added in their vicinity. The quasi-particle picture also accounts for the dressing of quasi-particles by a polarization cloud of all surrounding particles. Thus, the QP bands energies will involve systematic corrections of the KS-DFT band energies by a term that involves the expectation value of the self- energy:

In general, the self-energy will be complex, and the imaginary part of the resulting complex quasiparticle band energies can in principle be used to calculate electron and phonon mediated quasi-particle lifetimes [28].

All in all, density functional theory in combination with a many-body approach based on Green’s functions is a very powerful combination of methods, which can be used to compute most relevant materials properties to a sufficiently high level of accuracy. The main problem for inexperienced newcomers is therefore no longer the availability such methods, but a proper knowledge when the combination of these different methods will fail in practice.

Let us now discuss the last group of materials properties to be covered in this article, which is the linear dielectric response. We have seen in previous sections that the frequency-dependent dielectric constant $\varepsilon (\omega )$ and the related absorption coefficient $\alpha (\omega )$ are important inputs for conventional solar cell simulations.

The dielectric response and thus the dielectric constant can be computed in a straightforward manner using linear response time-dependent density functional theory (LR-TDDFT) [30]. This approach determines the transition probability between DFT bands using the transition dipole matrix elements. The corresponding transition energies are the energy differences between the KS-DFT bands.

In general, this is not a very good approximation. The KS-DFT band energies are typically far off as we saw before, and they should be systematically replaced with the corresponding GW quasiparticle energies. But even after using GW there is no guarantee that the transition probabilities based on KS-DFT wave functions are meaningful. Most importantly LR-TDDFT misses out on a proper description of the interaction between a negatively charged excited electron and a positively charged hole, which is created when the excited electron has left the valence band. This interaction strongly influences the dielectric properties and the lifetime of the excited states. In materials with strong particle-hole interactions, the two types of quasi-particles can form bound excitonic states, which binding energies that can be considerably lower than the corresponding transition energy. Silicon f.e. is known to exhibit strong exciton effects in its absorption spectrum, which makes it notoriously hard to model.

The best way to describe excitons and the corresponding proper dielectric response of a material is provided by the Bethe-Salpeter equation (BSE). The latter is the two-particle equivalent of the one-particle Dyson equation. As such it describes the electron-hole attraction and the formation of an exciton, which is missing from other schemes. The general numerical framework of the Bethe-Salpeter equation can be stated in the following form [31]:

Equation (18) provides an eigenvalue problem to determines the coupling coefficients ${A_n}({{\boldsymbol{k}_c},{\boldsymbol{k}_v}} )$ of the electron-hole pairs, and to determine the corresponding eigenvalues ${\Omega _n}$. In general, the solution of this eigenvalue problem involves the diagonalization of a non-Hermitian matrix, which means that the eigenvalues ${\Omega _n}$ may be complex. However, many standard implementations of the Bethe-Salpeter use approximations which yield real eigenvalues ${\Omega _n}$, only.

The kernel $K({\boldsymbol{k}_c^{\prime},\boldsymbol{k}_v^{\prime};{\boldsymbol{k}_c},{\boldsymbol{k}_v}} )$ in Eq. (18) describes the interactions between the various electrons and holes that contribute to each two-particle eigenstate ${\psi _n}({{\boldsymbol{r}_e},{\boldsymbol{r}_h}} )$. As indicated in Eq. (19), these two-particle eigenstates can be expanded into a one-particle basis of electron states ${\varphi _{{\boldsymbol{k}_c}}}({{\boldsymbol{r}_e}} )$ and hole states ${\varphi _{{\boldsymbol{k}_v}}}({{\boldsymbol{r}_h}} )$ with eigenvalues $\tilde{E}({{\boldsymbol{k}_c}} )$ and $\tilde{E}({{\boldsymbol{k}_v}} )$, and they may be determined using density functional theory (DFT) [27]. The interaction kernel can also be written as the difference $K = V - W$ between the Coulomb kernel $V({\boldsymbol{k}_c^{\prime},\boldsymbol{k}_v^{\prime};{\boldsymbol{k}_c},{\boldsymbol{k}_v}} )$ and the screened kernel

The frequency dependence of this screened potential is usually neglected, assuming an instantaneous type of screening.

Equation (18) is best solved in reciprocal space, where further approximations can be made for the kernel K to simplify the numerical procedures. The macroscopic dielectric constant itself can be computed directly using [27]

Figure 6 shows the corresponding absorption spectrum for a CH₃NH₃PbI₃ perovskite, which was computed with Yambo [28] using a 8 × 8 × 8 k-point grid. The linear response and BSE calculations both use DFT bands, but with shifted conduction band energies to reproduce the G₀W₀ band gap. The BSE results are visibly red-shifted due to the electron-hole interaction, and they are in much better agreement with the experimental results [32]. The overestimated absorption intensity as well as the lack of smoothness are a result of insufficient k-point sampling. This turns out to be a major problem with all BSE calculations, because these methods scale quite badly with the number of electrons in the system.

 

Fig. 6. Absorption spectrum of the CH₃NH₃PbI₃ perovskite as computed with DFT and BSE using a 8 × 8 × 8 k-point grid. Experimental data [32] is shown for comparison.

Download Full Size | PDF

The selected ab initio simulations shown in this section are by far not exhaustive. Static dielectric constants, with or without phonon contributions, can also be computed quite efficiently using Density Functional Perturbation Theory (DFPT) [22]. Electron-phonon relaxation times [28] can be computed at the DFPT or GW level, but both types of calculations involve a huge numerical effort. But together with the effective electron/hole masses taken from DFT or GW band structures [6], we can use these relaxation times to compute the electron/hole mobilities.

5. Summary and outlook

Solar cell device simulations are an interesting, but also quite challenging topic of applied physics, which contribute to a rapidly growing field that provides us with a whole variety of efficient and environmentally benign energy technologies.

Although we made it clear that a conventional solar cell simulation can entirely be run based on experimental data, we have also pointed out several distinct advantages of using numerically determined ab initio data sets as part of simple model approaches, or as part of full-fledged device simulations. The best news might be that this whole field is not overcrowded yet, and therefore there are still plenty of opportunities to make important contributions, ranging from the study of fundamental light-matter interactions over the study of the flow of charge carriers and their extraction at the contacts, to the systematic improvement of ab initio methods in terms of accuracy and efficiency.

Other opportunities for the use of ab initio methods may arise in the context of functional augmentations of solar cells. As pointed out in [33], the current trend in solar cell design points towards more advanced photon management schemes, such that future types of solar cells might even become as complex as microchips. A general description of the implementation of the necessary light guiding features in the framework of conventional solar cell device simulations was given in [34].

In recent years, we have explored two of the most promising augmentation approaches and implemented them into solar cell simulations. The first strategy implies the use of plasmonic nanoparticles to concentrate incoming light, and to scatter that light forward into the active layer of the solar cell, rather than reflecting it. Review [35] summarizes the most recent trends in this field. The description of embedded plasmonic nanoparticles using a combination of classical electromagnetic simulation methods and ab initio data were also the topic of two extended book chapters [36]. The easiest way to estimate the influence of plasmonic nanoparticles on the flow of incoming photons goes over modifications of the light current in the diode equation, where one must work out a modified absorption coefficient to cater for the presence of plasmonic nanoparticles, see Eqs. (3–4). Such a modified absorption coefficient may best be obtained from standard electromagnetic simulation techniques, see [36]. A full implementation of plasmonic nanoparticles in a conventional type of solar cell simulation would require techniques that are described in [34], and which go well beyond the simplest 1D type of device simulations described in Section 3.

Another augmentation strategy is the use of frequency-conversion layers containing rare-earth ions, which has first been suggested in [37]. The idea is to avoid the spectral losses shown in Fig. 1(b) by using two-photon processes to combine lower energetic photons and generate a photon with energy close to the band gap energy, or to split a higher energetic photon into two lower-energetic photons with energies close to the band gap energy. This leads to corrections of the photon flux $\mathrm{\Phi }(\omega )$ and thus of the corresponding light current, as can be seen again from Eqs. (3–4). The additional photon flux generated by frequency-conversion layers may be calculated from the fundamental rate equations underlying such processes, which was described in [5] and [38]. In general, the use of ab initio methods to describe highly correlated systems like rare-earth ions is certainly quite challenging, but it could become a real blessing for the simulation of frequency-conversion processes.

That said, one should always keep in mind that ab initio methods are not, and they will never be, a magic tool that can do everything in photovoltaics. Many processes in a solar cell are of macroscopic nature, and these processes must be simulated on the level of a device simulation rather than on the level of a microscopic quantum mechanical simulation. On the other hand, we have pointed out that some of the most important processes like carrier generation rates, carrier recombination rates or the trapping of charge carriers will strongly depend on atomistic data. Here the use of ab initio methods in combination with some of the simplest model approaches may already provide experimentalists with important reference data to check and verify the validity of their measurements on new types of photovoltaic devices.

Furthermore, advances in hardware and software development will continuously push the limits of ab initio materials simulations. Electron-phonon interaction, non-linear processes, impurities and real-time optical excitation are already within reach for those, who have generous access to larger supercomputers. One such example are BSE calculations outside the optical limit ($\boldsymbol{q} \ne 0)$ [28].

Still, when reliable experimental data is available for a device simulation, then the corresponding parameter sets are usually sufficient to carry out device simulations. However, ab initio data can be always used to provide parameters that might have been measured with very low precision. Based on the accuracy of these materials parameters, one can carry out extensive numerical studies to identify technical bottlenecks for the overall efficiency of photovoltaic devices, which may become the basis for a virtual prototyping of novel types of solar cells that is entirely run on a computer.

But we think that the real strength of ab initio simulations consists of the fact, that they provide parameter sets for ideal materials and specific ideal processes. In combination with device modeling approaches, we can use this information to easily obtain very useful estimates of overall device efficiencies and check, whether a given solar cell really works as expected. If measured values deviate quite largely from the predictions made using ab initio data (in combination with simplified modeling approaches), then we might either talk about a real breakthrough in the field, or it might be high time to check the lab equipment!

As always, one should adhere to the KISS-principle. Decide on a simulation strategy and about the information you need from such a simulation. Then choose the simplest simulation approach that will yield the desired information. More complex simulations will always require more free parameters, some of which may not always be very well known, and such a simulation will provide you with less reliable information compared to simpler, but well understood simulations.