Solcore: a multi-scale, Python-based library for modelling solar cells and semiconductor materials

Computational models can provide significant insight into the operation mechanisms and deficiencies of photovoltaic solar cells. Solcore is a modular set of computational tools, written in Python 3, for the design and simulation of photovoltaic solar cells. Calculations can be performed on ideal, thermodynamic limiting behaviour, through to fitting experimentally accessible parameters such as dark and light IV curves and luminescence. Uniquely, it combines a complete semiconductor solver capable of modelling the optical and electrical properties of a wide range of solar cells, from quantum well devices to multi-junction solar cells. The model is a multi-scale simulation accounting for nanoscale phenomena such as the quantum confinement effects of semiconductor nanostructures, to micron level propagation of light through to the overall performance of solar arrays, including the modelling of the spectral irradiance based on atmospheric conditions. In this article, we summarize the capabilities in addition to providing the physical insight and mathematical formulation behind the software with the purpose of serving as both a research and teaching tool.


Introduction
Computer aided design and device models are valuable tools when developing and evaluating photovoltaic solar cells.Laboratory scale tests can be usefully compared against detailed models that account for all relevant processes or with ideal, thermodynamically limited behaviour.Over the years, and with different degrees of sophistication, many pieces of software have been developed and published to tackle different aspects of solar energy research.For example, to calculate the solar spectrum as a function of the atmospheric conditions a traditional solution is to use SMARTS [1]; the light absorption profile in the solar cell or even at module level could be addressed by OPTOS [2] or OPAL2 [3]; while to solve the transport equations of a solar cell one could use PC1D [4], SCAPS [5] or Quokka [6].Several free and commercial programs, not specifically designed for solar energy research, have also been used historically, including AFORS-HET [7], Nextnano [8], ATLAS [9] or SENTAURUS [10], with the first two focused on the device and semiconductor properties and the latter two also solving the optics of the solar cells, among other properties.An extensive list of software for solar energy research -both online calculators and downloadable programs -has been compiled by PV Lighthouse [11].In general, programs like ATLAS and SENTAURUS provide a general-purpose, easy to use interface -often solving multi-physics problems, such as electrical transport coupled with thermal transport -to the detriment of performance.On the contrary, specific programs like AFORS-HET or PC1D are extremely fast and efficient, but limited in the problems they can solve, in this case 1D heterostructures and solar cells.
Apart from a few exceptions, such as PVlib [12], all these solvers are high-level, self-contained applications.While users can provide their own inputs and, in some cases, access the source code of the programs and customize some aspects of them, they are not designed with that purpose in mind.
Solcore is a multiscale, modular simulation framework for solar energy research, written mostly in Python.Solcore evolved from SOL, a Fortran-based, quantum well solar cell solver developed by Nelson and Connelly [13], with the explicit purpose of simplifying its integration in other programs, its expansion with custom routines and algorithms, and being didactic and informative.It is a teaching and learning tool as much as a rigorous research tool.Solcore is also extremely flexible.It integrates several algorithms for the multiscale modelling of semiconductors and solar cells, allowing a user without access to other methods, but with some experience with Python, to solve many different problems out of the box.On the other hand, most of its functionality can be interfaced with external tools, optimized to solve a specific problem, which are more advanced and accurate or that use an approach not considered in Solcore.The most recent version has been released under the GNU Lesser General Public License (GNU-LGPL) and can be found on GitHub [14].
Solcore's capabilities can be grouped into four categories: materials science (Section 2), light sources (Section 3), solar cells (Sections 4 to 6) and large-scale calculators (Section 7), each of them tackling a different area and scale relevant for research in solar energy.Figure 1 shows how these parts relate to each other and summarizes some of their content.

Materials science
The materials science modules in Solcore deal with the retrieval and calculation of material properties as well as those of quantum nanostructures, in particular quantum wells.They form the building blocks necessary to create the structures and calculate the performance of full solar cell devices.While focused on its application for solar cells, this part of Solcore is widely applicable to any research area related to semiconductor materials, as a way of managing the material properties, customising them and using them in other calculations.
It should be noted that, in reality, electronic and optical properties are not independent but related to each other via the material band structure.As the case in most programs, Solcore uses a non-consistent approach, with independent electronic and optical parameters obtained from different sources.The reader should consider full band structure methods, like pseudopotential or tight-binding, for a consistent set of electronic and optical properties.

Parameters database
The parameters database contains the basic properties of many semiconductor materials, including silicon, germanium and many III-V semiconductor binary and ternary alloys.Among other parameters, it includes the energy bandgap, the electron and hole effective masses, the lattice constants and the elastic constants.
The main sources of data are the article by I. Vurgaftman on Band parameters for III-V semiconductors [15] and the Handbook Series on Semiconductor Parameters by Levinshtein et al. [16].The carrier mobility calculator is based on the empirical low-field mobility model by Sotoodeh et al. [17] and it is available only for some materials where the inputs for the model are available.
There are two methods for retrieving parameters from the database.The first one consists simply of getting the data using the get parameter function with the required inputs.For example: get_parameter ( " GaAsP " , " band_gap " , P =0.45 , T =300) will return the bandgap of GaAsP for a phosphorus concentration of 45% at a temperature of 300 K, equal to 1.988 eV.This method only uses the existing data.
Another method is to create a material object which will contain all the properties existing in the database for that material, as well as those included as input, which will override the value of the database parameter, if it exists.The following example creates a GaAs object Fig. 1: General structure and workflow of Solcore.
and an AlGaAs object, using a custom electron effective mass in the latter: GaAs = material ( " GaAs " ) ( T =300 , Na =1 e24 ) AlGaAs = material ( " AlGaAs " ) ( T =300 , Al =0.3 , Nd =1 e23 , e f f _ m a s s _ e l e c t r o n =0.1) Now, any parameter -including the custom ones -are attributes that can be easily accessed and used anywhere in the program.For example GaAs.band gap is the GaAs bandgap and AlGaAs.latticeconstant is the AlGaAs lattice constant, both at the composition and temperature chosen when creating the objects.
Figure 2 shows the well-known bandgap vs. lattice constant map of all semiconductor materials and alloys (only ternary compounds) currently implemented into Solcore.However, any other material can be used in all of the Solcore functions, as long as the necessary input parameters are provided.This can be done by overriding all the properties of an existing material during the creation as above, or adding it as an external material in the configuration files.

Optical Properties Database
In order to calculate and model the optical response of potential solar cell architectures and material systems, access to a library of accurate optical constant data is essential.Therefore, Solcore incorporates a resource of freely available optical constant data measured by Sopra S. A. and provided by Software Spectra Inc. [18].The refractive index n and extinction coefficient k are provided for over 200 materials, including many III-V, II-VI and group IV compounds in addition to a range of common metals, glasses and dielectrics.
Any material within the Sopra S. A. optical constant database can be used with the "material" function described above, but they will have only the optical parameters n and k.In the case of materials that are in Fig. 2: Materials with most parameters included in Solcore's database (excluding the optical properties).both databases, the keyword "sopra" will need to be set to "True" when creating the material.Once a material is loaded its n, k and absorption coefficient data is returned by calling the appropriate method, for example SiO2.n(wavelength) and SiO2.k(wavelength).For certain materials in the database, the optical constants are provided for a range of alloy compositions.In these cases, any desired composition within the range can be specified and the interpolated n and k data is returned.Several examples of materials created from the Sopra database are shown in the Listing 1.
Figure 3 highlights example output from the Sopra materials library with 3a showing GaAs and Ge optical constant data and 3b showing interpolated AlGaAs extinction coefficients for a range of aluminium fractions.

Quantum solver
The electronic band structure of semiconductor materials is responsible for their light absorption and emission properties as well as for many of their transport properties, ultimately depending on the carriers' effective masses.These properties are not intrinsic to the material, but depend on external factors, too, most notably the strain and the quantum confinement.
Given the crystalline nature of most semiconductor materials, there will be strain whenever two materials with different crystal lattice constants are grown on top of each other pseudomorphically.Even those with the same lattice constant might be under stress due to other effects such as the presence of impurities or if used at different temperatures having dissimilar thermal expansion coefficients.Quantum confinement, in turn, takes place when the size of the semiconductor material in one or more dimensions is reduced to a few nanometres.In that situation, the energy levels available to the carriers become quantized in the direction of confinement, also changing the density of states.Both conditions take place simultaneously when dealing with strain-balanced quantum wells (QW).
Quantum wells -and more recently quantum wireshave been employed to tune the absorption properties of high efficiency solar cells for the past two decades.The need for appropriate tools to study them in the context of photovoltaics led to the development of the simulation models that were the seed of Solcore [19,20,21,22].As strained materials with quantum confinement, special care must be taken to obtain a sensible set of parameters for the QW structures, including the band edges with confined energy levels, the effective masses and the absorption coefficient.
Solcore's approach for evaluating the properties of QWs involves calculating first the effect of strain using a 8-band Pikus-Bir Hamiltonian (Section 2.3.1),treating each material in the structure as bulk, and then using the shifted bands and effective masses to solve the 1D Schödinger equation, after a proper alignment between layers (Section 2.3.2) [23].Finally, the absorption coefficient is calculated based on the 2D density of states, including the effect of excitons (Section 2.5).

Bulk 8-band k•p calculator
There are many numerical methods to calculate the band structure of a material with a varied degree of sophistication and complexity, such as the tight binding, pseudopotential, Green's function or k• p methods.Solcore includes a modified 8-band Pikus-Bir Hamiltonian to calculate the band structure of bulk materials under biaxial strain [24], considering the coupling between the conduction, heavy hole, light hole and splitoff bands.The eigenfunctions Ψ and eigenstates E are the solutions of the following equation, where H is the Pikus-Bir hamiltonian: Here, the sub-diagonal elements are just the complex conjugate of the corresponding upper elements.Di-Fig.3: Example output from the Sopra S. A. optical constant database accessed from Solcore.(a) Refractive index and extinction coefficient data for GaAs (solid lines) and Ge (dashed lines).(b) Interpolated AlGaAs extinction coefficient data with aluminium fractions ranging from 10 to 100%.agonal terms have three components: the information about the unstrained band edges, a kinematic term, and a strain term.As an example, the term E cb is given by: where E c0 is the position of the unstrained conduction band edge, γ c a modified Luttinger parameter, the k i the momenta in the different directions of the reciprocal space, a c the conduction band hydrostatic deformation potential and the ij s the strain tensor components.Offdiagonal terms have similar expressions, not involving the unstrained band edges.A detailed description of all these terms and their origin can be found in Tomic [24].
This system is readily solved for the given k i using Numpy's linalg.eigfunction, which provides as output the eigenfunctions and the corresponding eigenvalues.Typically, we are interested in the new band edges due to the effect of strain and the resulting effective masses, given by the curvature of the bands near k i = 0. Fig. 4 shows an example of the bands calculated in this way for the case of a strained InGaAs layer grown pseudomorphically on GaAs, and the resulting dependence of the effective mass with the indium content of the layer.Notice that, due to the effect of strain, the heavy and light hole bands are no longer degenerate at the gamma point Γ (k = 0).

1D Schrödinger equation
Once the new band edges and effective masses for each of the materials forming the quantum well structure are known, the quantum properties can be calculated by solving the 1-dimensional Schrödinger's equation.Solcore uses the method described by Frensley [26], which allows calculation of the eigenvalues and eigenvectors of an arbitrary potential.However, only closed and periodic boundary conditions are implemented.Solcore does not incorporate the Quantum Transmitting Boundary Method (QTBM) described by Frensley, meaning that unbounded states will be, in general, not correct.
A tridiagonal matrix is constructed by writing the variable effective mass Schrödinger's equation over a series of mesh points.The eigenvalues of the matrix correspond to the allowed energy levels of the system.Thus, the system of equations to solve over the mesh points is given by: where s i and d i depend on the mesh spacing ∆ and the position dependent potential V i and effective masses m i as: This system is solved using the tools available in the Scipy package for solving sparse linear systems of equations, in particular sparse.linalg.eigs.
Fig. 5 shows two examples of the band profile and wavefunctions calculated by Solcore.The first one (Fig. 5a and b) is a single InGaAs QW with GaAs interlayers and GaAsP barriers.The strain and quantum confinement shift the light hole valence band (dashed line) with respect to the heavy hole valence band (continuous line).In the GaAsP barriers, under tensile strain, this shift is in the opposite direction to the shift inside the QW, which is under compressive strain and experiences the effects of confinement.Fig. 5c and d compares the position of the energy levels predicted by Solcore in this structure with the more rigorous treatment of the 1-dimensional 8-band kp solver implemented by Nextnano++.As shown, the electron energy levels are barely modified, but hole levels are shifted due to the coupling between the bands [8].The in-plane dispersion when using a 8-band solver will no longer be parabolic and we would expect this to have an impact into the absorption coefficient and especially the selection rules for the transitions due to the band mixing effects (see Section 2.5).
Finally, Fig. 5e and f shows 1D the local density of states (LDOS) of a a multi-QW structure with thin barriers, including a Lorentzian broadening of 5 meV.For the electrons, there is strong coupling between the QWs, resulting in a range of energies for the ground state.The heavy hole ground states are too deep, resulting in lower coupling between wells.This figure also shows the artefacts in the LDOS of unbound states due to the use of closed boundary conditions rather than QTBM.

Critical-Point Parabolic-Band Optical Constant Model
Understanding the optical response of both established and novel materials is crucial to effective solar cell design.To efficiently model the complex dielectric function of a material Solcore incorporates an optical constant calculator based on the well-known Critical-Point Parabolic-Band (CPPB) formalism popularised by Adachi [27,28,29].In this model, contributions to 2 (ω) from critical points in the Brillouin Zone at which the probability for optical transitions is large (van Hove singularities) are considered.The transition probability for such transitions is proportional to the joint density of states (JDOS) J cv (ω), which describes the number of available electronic states between the valence and conduction bands at given photon energy.The imaginary part of the complex dielectric function is related to the JDOS by: Where | c|p|v | is the momentum matrix element for transitions from the valence band (v) to the conduction band (c).Critical point transitions are considered at the following points of symmetry in the band structure: E 0 corresponds to the optical transition at the Γ point and E 0 + ∆ 0 to the transition from the spin-orbit split off band to the conduction band at the Γ point.E 1 and E 1 + ∆ 1 denote the transitions from the valence heavyhole (HH) band and the valence light-hole (LH) band respectively to the conduction band at the L point.The E 0 triplet and E 2 transitions occur at higher energies, between the HH band and the split conduction bands at the Γ point as well as across the wide gap X valley.The model also includes contributions from the lowest energy indirect band-gap transition and the exciton absorption at the E 0 critical point.The contributions listed above are summed to compute the overall value of 2 (ω).The real and imaginary components of the overall complex dielectric function (ω) = 1 (ω)−i 2 (ω) are then related via the Kramers-Kronig relations; The CPPB model included with Solcore also incorporates a modification to the critical point broadening present in Adachi's description, which is shown to produce a poor fit to experimental data in the vicinity of the E 0 and E 1 critical points [30].To give a more accurate description of the broadening of the optical dielectric function, Kim et al. proposed that a frequencydependent damping parameter be used to replace the damping constant given by Adachi at each critical point [31,32]; Where Γ is the damping constant used by Adachi and α describes the shape of the lineshape broadening with α = 0 producing purely Lorentzian character and α = 0.3 producing a good approximation to Gaussian broadening.
The Solcore module absorption_calculator contains the CPPB model within the Custom_CPPB class.The class offers a flexible way to build up the optical constant model by adding indivdual critical point contributions through the Oscillator structure type within Solcore.In addition to the oscillator functions described by Adachi the Custom_CPPB class also provides additional oscilator models and the Sellmeier equation for describing the real part of the dielectric function for non-absorbing materials [33].For example, the code in Listing 2 calculates the complex dielectric function of GaAs. Figure 6 shows the real and imaginary components of the complex dielectric function of GaAs as calculated by the Custom_CPPB class.The model, using a set of parameters for GaAs similar to those specified in [28], shows excellent agreement with the experimental data taken from Palik [34].For a recent demonstration of Solcore's CPPB model, please refer to the discussion on Wilson et al. [35].

QW absorption calculator
For modelling the optical properties of QWs we use the method described by S. Chuang [36].The absorption coefficient at thermal equilibrium in a QW is given by: where |I en hm | 2 is the overlap integral between the holes in level m and the electrons in level n; H is a step function, H(x) = 1 for x > 0, 0 and 0 for x < 0, ρ 2D rmn is the 2D joint density of states, C 0 a proportionality constant dependent on the energy, and F the excitonic contribution, which will be discussed later.
Here, n r is the refractive index of the material, m rmn = m en m hm /(m en + m hm ) the reduced, in-plane, effective mass and L an effective period of the quantum wells.The in-plane effective mass of each type of carriers is calculated for each level, accounting for the spread of the wavefunction into the barriers as [37]: This in-plane effective mass is also used to calculate the local density of states shown in Figure 5b.In Eq. 12, |ê•p| 2 is the momentum matrix element, which depends on the polarization of the light and on the Kane's energy E p , specific to each material and determined experimentally.For band edge absorption, where k = 0, the matrix elements for the absorption of TE and TM polarized light for the transitions involving the conduction band and the heavy and light holes bands are given in Table 1.As can be deduced from this table, transitions involving heavy holes cannot absorb TM polarised light.
In addition to the band-to-band transitions, QWs usually have strong excitonic absorption, included in Eq. 13 in the term F nm .This term is a Lorenzian (or Gaussian) defined by an energy E nmx,j and oscillator strength f ex,j .It is zero except for m = n ≡ j where it is given by Klipstein et al. [38]: Fig. 7: Absorption coefficient of a single 7.2 nm-thick In-GaAs QW with 3 nm GaAs interlayers and GaAs 0.9 P 0.1 barriers as a function of the indium content.
Here, ν is a constant with a value between 0 and 0.5 and σ is the width of the Lorentzian, both often adjusted to fit some experimental data.In Solcore, they have default values of ν = 0.15 and σ = 6 meV.R is the exciton Rydberg energy [36].Fig. 7 shows the absorption coefficient of a range of InGaAs/GaAsP QWs with a GaAs interlayer and different In content.Higher indium content increases the depth of the well, allowing the absorption of less energetic light and more transitions.

Light sources
Transforming sunlight into electricity is the final goal of any solar cell and it is therefore necessary to have a convenient way of creating, manipulating and modifying the properties of the spectrum of the light.Ideally, solar cells will be designed and evaluated under a standard solar spectrum -e.g. the air mass 1.5 direct solar spectrum, AM1.5D -but practical light sources are not standard.More often than not, we are interested in modelling the performance of a solar cell under the experimental spectrum of a solar simulator or lamp in a laboratory, simulated data calculated from atmospheric conditions (temperature, humidity, aerosol content, etc.) or even under real irradiance data measured at different locations worldwide.This can then be compared with the experimental performance or tailored to work best under certain conditions.
The Solcore module light source is designed to deal easily with different light sources.It has direct support for: -Gaussian emission, typical of lasers and light emitting diodes.-Black-body radiation, characteristic of halogen lamps defined by a temperature, but also used very often to simulate the spectrum of the Sun, very close to a black body source at 5800 K. -Standard solar spectra: the extraterrestial spectrum AM0 and the two terrestial ones, AM1.5D and AM1.5G as defined by the ASTM G173 -03(2008) standard.-Irradiance models, using location, time and atmospheric parameters to calculate a synthetic solar spectrum.Solcore includes two models: SPECTRAL2, fully implemented in Python, and an interface to SMARTS binaries (which need to be installed separately), which greatly simplifies its use in batch mode.-User-defined irradiances, provided externally from a database or any other source, allowing for maximum flexibility.
The syntax in all cases is simple and intuitive considering the type of source that needs to be created.In the case of the irradiance models, which often have a large number of inputs, Solcore defines a set of default values, so only those that are different need to be provided.The code in Listing 3 illustrates the creation of several light sources using the minimum required input in each case.A plot of those light sources is shown in Figure 8.  Once created, specific parameters of the light sources can be easily modified without the need for creating the source from scratch.That is particularly useful for the irradiance models, where we might be interested in getting the spectrum as a function of a certain parameter (e.g. the hour of the day, or the humidity) without changing the others.For example, smarts.spectrum(HOUR=11)and smarts.spectrum(HOUR=17)will provide the spectrum of the SMART light source defined above calculated at 11h and at 17h, respectively; all additional parameters have the default values.This method has been used to model experimental solar irradiances measured by different spectroradiometers based on the local atmospheric conditions [39].
A final, very convenient feature of the LightSource class is the ability to request the spectrum in a range of different units.The default is power density per nanometer, but other common units are power density per eV or photon flux per nanometer, among others.While these unit conversions are straightforward, it is often an initial source of errors due to missing constants or incompatible magnitudes.
The light source module has been described in the context of the solar spectrum, but it can be applied broadly where there is spectral data involved, such as the fitting of photoluminescence, electroluminescence or Raman spectra.

Optical solvers
The purpose of the optical solvers is to obtain the fraction of incoming light reflected, absorbed and transmitted in a solar cell as a function of the wavelength of the light and the position inside the structure.Solcore includes three models to tackle this problem: Beer-Lambert law (BL), transfer matrix method (TMM) and rigorous coupled wave analysis (RCWA).At the moment, Solcore does not have explicit support for light trapping effects using general textured surfaces, which are usually present in silicon solar cells.However, this can be implemented to a large extent using the RCWA method, although not very efficiently.Additionally, the reflected, absorbed and transmitted light can be calculated externally and then provided as input to Solcore to obtain the electrical properties of a solar cell structure, giving it full flexibility.
All the optical solvers apply to the solar cell structure as a whole, providing as output the fraction of light reflected (R(λ)), transmitted (T (λ)) and absorbed per unit length at a depth z from the front surface (A(λ, z)).
Figure 9 shows a comparison of the quantum efficiency of a thin GaAs solar cell with a distributed Bragg reflector (DBR) and an array of TiO 2 nanoparticles (NP) on top calculated using the three optical solvers described in the following sections.Supplementary Information shows the full code needed to produce these curves [40].

Beer-Lambert law (BL)
This is the simplest model to calculate the absorption in a multi-layer structure.It ignores all reflection at the interfaces -the front surface reflection can be provided externally, and is zero otherwise at all wavelengthsand the absorption per unit length as a function of the wavelength λ and the position z in layer n is given by: where α n is the absorption coefficient of layer n, w n its thickness and z n the position of the beginning of the Fig. 9: Comparison of the quantum efficiency of a thin GaAs solar cell with a distributed Bragg reflector (DBR) and an array of TiO2 nanoparticles (NP) on top calculated using the three optical solvers.The BL model ignores the NP and the DBR and does not include the front surface reflection, overestimating the result at all wavelengths.TMM correctly accounts for the reflection and the DBR but cannot model the effect of diffraction due to the NP layer.RCWA takes into account scattering from the NP, although it is significantly more time consuming.layer.Due to its simplicity, the BL law is used widely in photovoltaics but in reality it is only applicable when the contrast in the refractive index between layers can be ignored and when there is strong absorption, reducing the effects of light reflection at the interfaces.

Transfer matrix method (TMM)
In order to evaluate the realistic optical behaviour of a solar cell design it is important to consider the interaction of incident electromagnetic (EM) radiation with a succession of both absorbing and non-absorbing planar layers.The combined optical response of such a layered structure is crucial when considering the minimisation of extrinsic front surface reflection losses [41], the emissivity in the mid-IR of low emissivity coatings for hybrid PV-thermal applications [42] and also when studying the optical constants and layer thicknesses of material using the experimental technique of spectroscopic ellipsometry.Therefore, Solcore evaluates the interaction of incident EM radiation through a layered structure using the TMM.The incident light radiation takes the form of homogeneous, electromagnetic plane-polarised waves and is represented by components describing the electric, E and magnetic field strengths, H : Where E and H denote the electric and magnetic field amplitudes respectively, N the complex refractive index, z the distance in the direction of propagation, ω the angular frequency and λ the wavelength of radiation.ϕ and ϕ represent arbitrary phase angles for both the electric and magnetic travelling wave components and are not independent of each other.The characteristic transfer matrix for evaluating the interaction between planar electric and magnetic waves at the interface of n thin films on a semi-infinite substrate is derived in detail elsewhere [43] and given by the equation: Where δ q is defined as the phase factor of the q th planar layer, η q the optical admittance of the q th layer and η m the optical admittance of the substrate.The layer closest to the incident medium is evaluated first before working through the n layer structure in order.The term δ q describes the phase shift required to translate the z coordinate of the E and H interactions by the thickness of each layer, q.The spectrally varying Fresnel coefficients describing reflection, transmission and absorption of the multi-layer structure can be calculated from the solutions to equation 23 at discrete wavelengths: The implementation of the TMM in Solcore uses the freely available tmm module developed by Byrnes [44].The multi-layer optical stack is built up using Solcore's Structure object.Some example code evaluating the TMM for a triple layer anti-reflection coating (ARC) on top of conventional multi-junction solar cell materials AlInP and GaInP is included in Listing 4.
Figure 10a depicts calculated reflection and transmission from the solutions to the characteristic TMM equation over a range of incident angles for an optimised triple layer ARC design reported in [41].The solid lines indicate the optical reflection at the front surface whilst the dashed lines correspond to the transmitted light into the substrate, in this case taken to be the optically thick top sub-cell material of GaInP.
In addition to using the TMM for calculating the reflection, transmission and absorption in a multi-layer optical stack it can also be applied to the popular spectroscopic technique of ellipsometry.This measures a change in the polarisation of incident light reflected at the surface of a sample.A more detailed description of ellipsometry and its uses can be found elsewhere [33].The measured values are expressed as the angles Psi (Ψ ) and Delta (∆), which are related to R s and R p (the Fresnel reflection coefficients for s and p-polarized light, respectively) by: As the ratio between R s and R p is a complex quantity, phase change information is contained within ∆ in Eq. 27.The change in phase of the reflected electric 1 The definition of the phase ∆ in the definition of ρ (Eq.27), which affects the sign of 2 in Eq. 29, is the same here as in [33].However, other conventions are in use; for instance, ∆ is defined with the opposite sign in [44].
and magnetic plane waves can be evaluated from the solution of Eq. 23 and is given by: The complex dielectric function of a sample can be calculated from the experimental ellipsometry results with the solutions to Eq. 27: Some example code, calculating the complex dielectric function from the ellipsometric response of a sample of Ge substrate is shown in Listing 5.The output of the model is shown in Figure 10b and is compared with experimentally obtained data at an incident angle of 79 • .Good agreement with the experimental data is observed when a thin 4.5 nm layer of germanium oxide (GeO 2 ) is included in the layer model.

Rigorous coupled-wave analysis (RCWA)
Finally, Solcore includes an interface to the S 4 solver (which must be installed separately), developed at Stanford University, in order to model solar cells with nanophotonic designs [45].S 4 is an implementation of RCWA, also sometimes referred to as the Fourier Model Method (FMM), which solves the linear Maxwell's equations in structures containing 2D periodicity.Structures with 2D periodicity can be found in advanced solar cell designs aiming, for example, to reduce the solar cell thickness by scattering the incoming light using a periodic diffraction grating at the front or rear of the absorbing layer [46,47].S 4 defines structures by creating a layer stack of the desired materials using Solcore's Layer and Junction classes, in which each layer's composition can be modified by adding circles, rectangles, ellipses or a generalised polygon made of a specified Solcore material.Each layer is assumed to be infinitely periodic in the x and y direction, and uniform in the z direction.A unit cell must be defined for the whole structure (using the size attribute in the user options), and each shape is placed at a specified location (and, where relevant, angular orientation) in the unit cell; as many shapes as necessary can be added, supplied as a list of dictionaries with the relevant parameters for each shape.An example for each type of shape supported by S 4 is shown in Listing 6.The size of the unit cell in the x and y direc-tions and the base size of the unit cell must be given in nm.The number of Fourier components used in the calculation [45] must also be specified; this is done using the orders attribute in the user options passed to the solver of choice.Listing 6: Creating geometry objects for use with the S 4 RCWA solver, which are added to Solcore Layer objects using the geometry attribute.

Single-junction solar cells
Solcore includes four solvers to calculate the electrical properties of a single-junction device.In order of increasing accuracy, these are: detailed balance, 2-diode equation, depletion approximation and Poisson-driftdiffusion.

Detailed balance (DB)
This solver calculates the electrical properties of the junction by balancing the elementary processes taking place in the solar cell, carrier generation and radiative recombination, using the formalism described by Araújo and Martí [48].The method is widely used by the photovoltaic community to calculate the limiting conversion efficiencies of the different solar cell architectures or materials.The simplest DB formulation only needs an absorption edge energy and an absorptivity value above that edge.Out of this, the carrier generation and radiative recombination are calculated for different internal chemical potentials, equal to the external electrical bias, in the ideal case.Solcore includes this basic model, but also allows the user to provide a more complex absorption profile.
The radiative recombination or thermal generation current J rad from the solar cell is calculated following the formalism described by Nelson et al. [20], considering all the possible paths of the light absorbed by the cell and the reciprocity relation between emission and absorption.The total radiative current, using the generalized Planck equation, is given by: where A(E, θ, s) is the probability that a photon of energy E will be emitted (absorbed) from the point s on the surface at an internal angle θ, and A f ront (E) and A back (E) are the combined probability of the photon to be emitted (absorbed) by the front and the back of the cell, respectively.
The different paths of the absorbed light are depicted in Fig. 11.Path A represents light that reaches the front surface within the escape cone (θ < θ c ) and that crosses the structure.Path B is the light that reaches the back surface outside the escape cone of the front surface (θ > θ c ), being totally internally reflected and crossing the structure twice.Light reaching the back surface within the escape cone (θ < θ c ) can either escape through the front (path C) or be reflected (path D).With these considerations, the contribution to the surface integral of the four terms will be [20]: These equations can be written in a more compact form by noting that r = 1 for θ > θ c .In this situation, B, C and D can be combined and the integral extended from 0 to 1, resulting simply in: The factor S back representing the area of the back of the cell has been omitted here as we are interested in the current density, and therefore independent of the area.Likewise, A f ront (E) will be given simply by the component A in Eq. 31: J rad (V, T cell ) will represent the radiative recombination of the cell at a bias V and temperature T cell while J rad (0, T a ) will be the carrier generation due to thermal radiation from the ambient at a temperature T a .Typically, T cell = T a .
Carrier generation in the solar cell due to the absorption of the solar irradiance H(E) can be written simply as: where R = r(E, 0) is the normal incidence reflection and A n is the normal incidence absorptivity of the cell, given by A(E) = 1−exp(−α(E)w) where w is the thickness of the junction and α is the absorption coefficient.Combining all these equations, the total current as a function of the bias calculated with the DB model will be given by: If T a = T cell and E >> k b T , Eq. 35 simplifies, resulting in: with J 01 the reverse saturation current, given by: While in these equations the term exp(−αw) is used, it should be noted that none of the two parameters alpha and w are needed as the product is calculated internally by Solcore from the normal incidence absorptivity A n (E), which is the value given as input: This is the simplest method for simulating the behaviour of a solar cell, using electrical components to model the different transport and recombination mechanisms of the device.The 2D model is widely applied when modelling solar cells at the most engineering end of the topic, when a detailed knowledge of the solar cell structure (layers, absorption coefficient, etc.) are not known or sought.It is often used to fit experimental IV curves and find approximate, general information on the solar cell quality without entering on the fundamental processes.It can provide valuable information to engineers, when designing solar modules for example, or for diagnostic purposes The complete form of the equation is: Generally, the photocurrent is modelled as a current source (J sc ), with radiative and non-radiative recombination modelled as two diodes with reverse saturation currents J 01 and J 02 , and ideality factors n 1 ≈ 1 and n 2 ≈ 2, respectively.The shunt resistance R sh accounts for alternative current paths between the contacts of the solar cell, being infinite in the ideal case, and the series resistance R s accounts for the other transport losses.The values of the saturation currents and ideality factors can, ultimately, be calculated from the material properties and device structure, as is done in the depletion approximation model (Section 5.4), but the 2D model allows them to be provided directly as input, obtained from a fit to experimental data, for example.They can also be calculated internally, using the DB solver to obtain J 01 and J sc , and then using a radiative efficiency coefficient to obtain J 02 .The radiative efficiency η is defined as the fraction of radiative current J rad at a given reference total current J ref : The reference voltage V ref can be written as a function of η and J ref as: On the other hand, the radiative coefficient can also be written as: Combining Eq. 41 and 42 and using the expression for the diode with ideality factor n 2 , J 02 can be written as: In the common situation of very large shunt resistance and V ref >> k b T /q, this equation further simplifies to: This process can, of course, be reversed to use knowledge of J 01 and J 02 at a given reference current to calculate the radiative efficiency of a solar cell, which is useful to compare different materials, technologies or processing methods.This was done by Chan et al. using J ref = 30 mA/cm 2 , obtaining η values of 20% for InGaP, 22% for GaAs, and 27% for InGaAs devices [49].It should be pointed out that this method is only valid under the assumption that J 01 corresponds only to radiative recombination and J 02 only to non-radiative recombination, which is generally true for QW solar cells and some III-V solar cells, like those made of GaAs or InGaP, but not for Si or Ge, for example.Other definitions of the radiative efficiency are based on the external quantum efficiency, the I sc and V oc of the cell, as described by Green [50].
Despite the simplicity of the 2-diode model, it is very useful to guide the design of new solar cells and explore the performance of new materials, such as dilute bismuth alloys [51], or to asses the performance of large arrays of solar cells, as will be shown in Section 7 [52].

Poisson-drift-diffusion (PDD)
This method solves the Poisson equation for the electrostatic potential coupled with the transport equations for electrons and holes and suitable boundary conditions in the steady state.It is the standard method for calculating the electrical properties of most semiconductor devices, including solar cells, transistors or light emitting diodes.It is also the only method included in most software packages for simulating semiconductor devices, such as PC-1D, Nextnano, SCAPS and AFORS-HET.
Figure 12a shows the flow chart of Solcore's PDD solver, which currently only solves the time-independent PDD equations (steady state).Any simulation starts by calculating the band structure under equilibrium conditions (no illumination or bias).If the simulation includes illumination, the photogeneration as a function of the position in the structure is calculated externally to the PDD solver using any of the models described in Section 4. To aid convergence, the solution at short circuit conditions is calculated by increasing the light intensity from zero to the nominal value in small steps.Similarly, the solution at any bias is obtained by solving the problem first at zero bias and then increasing it in small steps, using the previous solution as the initial condition for the next one.Re-meshing is performed several times during the simulation of the current-voltage characteristics (see section 5.3.3).
To calculate the internal quantum efficiency (IQE), a small differential increase is included in the photogeneration profile as a function of wavelength.The IQE is then calculated as the ratio between the resulting increase in the photocurrent and the increase in the photogeneration at that wavelength.This procedure is comparable to the actual experimental measurement of the quantum efficiency.

Solver assumptions and formulation
The Poisson's and drift-diffusion equations relate the electrostatic potential created by the free and fixed charges with the carrier densities and their variation across the structure due to generation, recombination and externally applied bias.The reader is referred to any semiconductor textbook for a detailed description of their derivation (see for example [53] or [54]).The solver uses the Boltzmann approximation for the carrier distribution with the following assumptions: -All carrier populations are in quasi-thermal equilibrium.-The mobility of carriers is independent of the electric field.
-Temperature is uniform.
-There are no magnetic fields.
As a consequence of the field-independent mobility, Solcore's PDD solver will be valid only in situations where electric field is not very strong.Poisson's equation relates the electrostatic potential φ and the electrical charges in the structure.In one dimension, it can be written as: where N A , N D , n and p are the density of ionized acceptors, donors, the density of free electrons and holes, respectively, and the dielectric constant.The current density equations account for the movement of carriers due to the electric field (drift component) and the carrier concentration gradient (diffusion component): where µ is the carrier mobility and F = −dφ/dx the electric field.Finally, the continuity equations ensure particle conservation, balancing the carriers that enter and leave any point of the structure.Under steady state, this means that the variation of the current must be equal to the generation G and recombination R processes, which are equal for electrons and holes since they are created and annihilated in pairs: Combining Eq. 46, 47, 48 and 49 gives: Poisson's equation (Eq.45) and the continuity equations (Eq.50 and Eq.51), together with the definitions for n, p and R, represent the complete system that needs to be solved in order to obtain the performance of the solar cell.The PDD solver included in Solcore uses the same discretization scheme used by PC-1D [4,55] taking as independent variables the electrostatic potential φ and the quasi-Fermi potentials for electrons and holes, φ n and φ p , respectively.These three variables are continuous across the whole structure and have comparable magnitudes in the voltage range.Details of the discretization process are included in [4,55], but they are based on minimising the total electrostatic energy in the case of the Poisson's equation, and in the Scharfetter-Gummel discretization scheme for the driftdiffusion equations [56].
The bulk recombination models included in Solcore are Shockley-Read-Hall recombination, radiative recombination and Auger recombination, as well as a surface recombination velocity model for the recombination at the contacts.The three recombination models are given by the following equations, respectively: with τ n and τ p the non-radiative lifetimes of electrons and holes, B the radiative recombination coefficient and C n and C p the Auger recombination coefficient for electrons and holes.By default, the radiative recombination coefficient is calculated internally by Solcore based on the absorption coefficient, as described by Nelson [54], and given by: where n i is the intrinsic carrier concentration, n s the refractive index and α the absorption coefficient.At the moment, Solcore's PDD solver cannot include interface charges or bandgap narrowing due to heavy doping and only implements Ohmic contacts.Additionally, it only includes local carrier recombination processes and therefore cannot deal with tunnel transport, which relates remote nodes.This is the only part of Solcore implemented in Fortran using quadruple precision variables in order to increase the numerical accuracy and improve convergence.

QWs in the PDD solver
Quantum wells have been developed in the context of solar cells mainly to tailor the absorption edge of the sub-cells in multi-junction devices to their optimum values [57].Typically, achieving the proper performance requires a delicate trade-off between carrier collection and light absorption [58,59] .Solcore includes a simplified QW structure in the PDD solver in order to calculate the performance of solar cells containing them.Contrary to other programs like Nextnano, Solcore does not solve the Schrödinger equation and the PDD equations self-consistently: first, the energy levels of the quantum wells are solved using a flat-band condition, considering also the strain in the materials, and then an effective band structure is used to solve the transport equations in a bulk-like fashion.This is illustrated in Figure 12b.
From the perspective of the PDD solver, the actual bandgap and electron affinity of each layer in a quantum well depend on the energy levels, i.e. the minimum energy for electrons is not the band edge of the conduction band, but the ground confined level.The same applies to holes, with the actual band edge being the maximum between the ground states of light holes and heavy holes.The resulting band profiles used in the PDD solver are shown in the right picture of Figure 12b.
To use QWs in the PDD solver, we create an effective electron affinity and bandgaps for all layers in the QW.For the barriers, the electron affinity and band gap are the same as they are in bulk, modified by the strain, if necessary.For interlayers, if present, it depends on what is higher, the band edges of the interlayer or the confined carrier levels.
The density of states and the absorption profile need to be modified in a similar way.For the density of states: -Barriers have the bulk density of states and absorption profile.-Interlayers only have the bulk density of states above the barrier and the bulk absorption from the barrier energy and zero below that.-Wells have all the density of states associated with the confined states and the bulk density of states above the barrier, while they have the absorption of the confined levels below the barrier energy and of the bulk above it.
These simplifications are similar to those in Nelson et al. [54] and in Cabrera et al. [60] and allow us to keep the bulk-like form of the carrier densities in the drift diffusion equations under the Boltzmann approximation.A more rigorous treatment will be necessary in the presence of tunnel transport across a supperlattice, tunnel escape from the QWs to the barriers -possible in the presence of high electric fields -and in the case of very deep QWs, when carrier escape from the less confined levels might be possible but not from the deeper ones.In these situations, a set of rate equations linking the different levels, as well as a self-consistent solution of the transport and Schrödinger equations would be required, besides using more advanced methods such as a non-equilibrium Green's functions (NEGF) formalism [61].

Mesh creation and dynamic meshing
The PDD solver discretizes the device into a finite number of mesh points at which to calculate the band structure, carrier densities, the generation and recombina-tion.The mesh can be static and homogeneous, but the default configuration and the one that results in the least number of mesh points, most accurate result and best convergence uses an inhomogeneous mesh with dynamic re-meshing.
There are two types of nodes in the mesh: masternodes and normal nodes.There are two masternodes at the ends of the device and two more at each side of an interface separated by 0.1 nm.These nodes are static and not affected by re-meshing.The rest of the nodes are automatically distributed depending on the distance to the closest masternode, as illustrated in Figure 12c.During re-meshing, nodes can be added or removed according to the following rules: An element will be divided into smaller elements by adding new nodes if any of the following statements is true: -The variation of the potentials or the carrier densities across the element is large.-The element is too close to the masternodes limiting the layer.-The element is too big for the region.A node will be removed if it fulfils all the following conditions: -It is not a masternode.
-The variation of the potentials or the carrier densities with respect the previous and next nodes is small.-It is not too close to the masternodes limiting the layer.-Removing it does not create an element too big for the region.
After the initial meshing and every time there is a re-meshing, the position of the nodes (except that of the masternodes) is smoothed to avoid having adjacent elements too different in size.This re-meshing process is controlled by a growth parameter, which can be adjusted by the user.
Using the inhomogeneous mesh in addition to the dynamic re-meshing ensures that those regions where material properties change abruptly are modelled with more detail, aiding the convergence.It also allows the modelling of devices which have layers with very different thickness, such as QWs a few nanometers thick and bulk absorbers of several microns, without increasing the number of nodes significantly.

Depletion approximation
The depletion approximation provides an analytical -or semi-analytical -solution to the Poisson-drift-diffusion equations described in the previous section applied to simple PN homojunction solar cells.Historically, it has been used extensively to model solar cells and it is still valid, to a large extent, for traditional PN junctions.More importantly, it requires less input parameters than the PDD solver and these can be easily related to macroscopic measurable quantities, like mobility or diffusion lengths.The DA model is based on the assumption that around the junction between the P and N regions, there are no free carriers and therefore all the electric field is due to the fixed, ionized dopants.This "depletion" of free carriers reaches a certain depth towards the N and P sides; beyond this region, free and fixed carriers of opposite charges balance and the regions are neutral.Under these conditions, Poisson's equation decouples from the drift and diffusion equations and it can be solved analytically for each region.For example, for a PN junction with the interface between the two regions at z = 0, the solution to Eq. 45 will be: where w n and w p are the extensions of the depletion region towards the N and P sides, respectively, and can be found by the requirement that the electric field F and the potential φ need to be continuous at z = 0. V bi is the built-in voltage, which can be expressed in terms of the doping concentration on each side, N d and N a , and the intrinsic carrier concentration in the material, n 2 i : Another consequence of the depletion approximation is that the quasi-Fermi level energies are constant throughout the corresponding neutral regions and also constant in the depletion region, where their separation is equal to the external bias qV .Based on these assumptions, Eq. 46 and 47 simplify and an analytical expression can be found for the dependence of the recombination and generation currents on the applied voltage.A full derivation of these expressions is included in Nelson [54].
Solcore's implementation of the depletion approximation includes two modifications to the basic equations.The first one is allowing for an intrinsic region to be included between the P and N regions to form a PIN junction.For low injection conditions (low illumination or low bias) this situation can be treated as described before, simply considering that the depletion region is now widened by the thickness of the intrinsic region.Corrently, no low doping level is allowed for this region.
The second modification is related to the generation profile, which in the equations provided by Nelson is given by the BL law (Eq.20) which has an explicit dependence on z and results in analytic expressions for the current densities.In Solcore, we integrate the expressions for the drift-diffusion equations under the depletion approximation numerically to allow for an arbitrary generation profile calculated with any of the methods described in Section 4. It should be noted that although the equations are integrated numerically this will not be a self-consistent solution of the Poisson-driftdiffusion equations, as is achieved by the PDD solver (Section 5.3).

Multi-junction solar cells
A complete photovoltaic solar cell can include one or more junctions, metal contacts, optical layers (including anti-reflective coatings and nano-photonic structures) and tunnel junctions.The junctions, in turn, might range from simple PN homojunctions to complex heterojunctions, including multi-quantum well structures.The solvers described in Section 5 only calculate the properties of single junction devices.To combine them into a multi-junction device, it is necessary to consider that the individual junctions are electrically connected in series and the potential coupling of light emitted by the wider bandgap junctions into those with smaller bandgap.The Suplementary Information includes a full step-by-step example of the modelling of a dual junction solar cell with QWs, anti-reflecting coating and a tunnel junction, calculating the external quantum efficiency, the IV characteristics under illumination and and the performance of the solar cell as a function of light concentration.

No radiative coupling
We first consider the case of no radiative coupling between junctions.This is a good approximation for materials which do not radiate efficiently or radiative materials working at low concentration, when the fraction of radiative recombination compared to non-radiative recombination is low.In this case, the IV curve of each junction can be calculated independently of each other and the current flowing through the MJ structure is limited by the junction with the lowest current at any given voltage.Series resistances defined for each junction are now added together and included as a single term.The operating voltage of each of the junctions is finally back-calculated and added together to get the voltage of the MJ device.
The pseudocode for this solver is: 1. Calculate the I j (V ) of each junction j in the structure.2. Find the current flowing through the MJ device as Calculate the voltage of each junction by interpolating its IV curve at the new current values, V j (I M J ), and the voltage dropped due to the series resistances, V Rs = R s I M J . 4. Calculate the total voltage at a given current as V M J = V Rs + j V j . 5. Interpolate the I M J (V M J ) and the I M J (V j ) to the desired output voltage values.
Fig. 13 shows the simulated IV curve of a 3J solar cell made of Ge/InGaAs/GaAsP.The electrical properties of the three junctions were calculated using the depletion approximation solver.In the dark (Fig. 13a) the voltages of each of the junctions at a given current add together, resulting in the total voltage of the MJ structure.The R s contribution to the voltage goes in the same direction as those of the junctions.Under illumination (Fig. 13b) the junction producing the lower current (the top junction in this case) limits the overall current of the MJ device.At zero bias, or even at some negative bias, the non-limiting junctions are positively biased, recombining all the photocurrent that cannot be extracted because of the limiting top cell.The contribution of the R s to the voltage of the MJ device is negative, resulting in a reduction of the fill factor and the overall efficiency of the solar cell.

With radiative coupling
Radiative coupling takes place when the light produced by a high bandgap junction due to radiative recombination is absorbed by a lower bandgap junction, contributing to its photocurrent and changing the operating point.It has been identified in numerous highly radiative materials such as quantum well solar cells and III-V MJ solar cells [62,63,64].It appears as an artefact during the QE measurements of MJ solar cells [65], but it is also an effect that can be exploited to increase the performance of MJ devices [57] and their tolerance to spectral changes, resulting in superior annual energy yield [66].  [20,66].It is implemented only for the DB junction model and for the 2D model when it is defined in terms of a radiative efficiency and the parameters calculated form the DB model.The current of a junction j including radiative coupling from the junction immediately above it j − 1 is given by: This current depends on two factors: the amount of radiation effectively emitted downwards, towards the lower junction, and the fraction of it that is absorbed and converted into electricity.If we ignore the possible reflection of light at the interface between both junctions, this current can be written by using a modified version of Eq. 30 that considers only the radiation emitted towards the back: with A j−1→j (E) given by: As discussed previously, any information related to total internal reflection will be contained in the r(E, θ) term, and therefore the integral over cos θ can be done from 0 to 1.In the case of thin junctions, some radiation could reach the next junction j +1.The coupled current in that case can be easily calculated by modifying Eq. 60 to account for the fraction of light absorbed by the junctions between the emitting junction and the junction of interest.In the general case, the current coupled into junction j will be given by: Radiative coupling might change the junction that is current limiting the MJ device, so the process to obtain the IV curves in this case proceeds in two steps.First, the IV of the junctions and the total IV are calculated without coupling.The resulting IV curves are then used as the initial conditions for the numerical solver that will calculate the correct voltage of each junction including the radiative coupling.
Fig. 14 shows the IV curve under the AM1.5G solar spectrum of a three junction solar cell (a) without and (b) with radiative coupling.Without coupling, the middle junction severely limits the current of the MJ solar cell.When coupling is enabled, the middle junction is still the limiting one but part of the excess current of the top junction is transferred to it, increasing its photocurrent by around 20 A/m 2 .Part of the radiative recombination is also transferred to the bottom cell, increasing slightly its photocurrent.In this case, given that the junction was overproducing current already, such coupling is only visible as an increase in the voltage.Altogether, the radiative coupling results in an enhancement of the V oc of 30 mV and of the efficiency η of 5.3%.This example uses junctions with 100% radiative efficiency to illustrate the effect, but this phenomenon is always present to some extent, becoming especially important under concentration [57,66].

Restrictions in the Junction definitions
Having multiple methods for modelling the junctions gives a lot of freedom and flexibility but it also imposes some restrictions in how and when they can be combined in order to create a MJ solar cell.The following compatibility rules apply: -When there is no radiative coupling and we are interested only in the dark IV characteristics, all junction models can be combined with each other.This tion models that can be used are DB and 2D, as long as the latter includes an absorptivity value.

Tunnel junctions
Solcore includes partial support for tunnel junctions.They represent an optical loss due to parasitic absorption in the layers, but also an electrical loss.At the moment, there are two models for tunnel junctions.The first one is a simple resistive model, where the tunnel junction is simply modelled as a series resistance.This approximation should be valid in most cases, but will break down if the current is close to or higher than the peak current density of the junction.
The second model is a parametric model, based on the simple formalism described by Sze [53].In this model, the total current of the tunnel junction will have three components: the tunnel current J T accounting for bandto-band transport, the excess current J ex related to transport across states inside the forbidden gap, and the diffusion current J D , which is the usual minority- As illustrated in Fig. 15, J P and V P are the peak current and voltage, J V and V V are the valley current and voltages, C is a prefactor of the exponent and J 0 the reverse saturation current.In this simple implementation, these 6 parameters need to be provided as inputs, and can be used as fitting parameters to reproduce experimental data.This allows to correctly account for the break down of the tunnel junction in situations when the current is above the peak current.
Solcore can also accept external IV data for the tunnel junctions and the implementation of the more rigorous, but still analytic model, described by Louarn et al. is currently under way in order to relate the tunnel currents with the actual materials and layer structure used in the tunnel junction definition [67].

Large circuit solver
When the two diode model is used to define the junctions in a MJ solar cell, then larger scale circuits can be constructed.Solcore includes two levels of large scale equivalent circuits: quasi-3D solar cell modelling and solar array modelling.Both solvers are based on the interface between Solcore and SPICE, allowing for a fast calculation of complex structures with many elements.

Quasi-3D solar cell model
The quasi-3D solar cell model included in Solcore uses a SPICE-based electrical network to model the flow of injected current through the solar cell, as depicted in Fig. 16.The plane of the cell is discretized into many elements, each of them representing a small portion of the cell.Depending on the location of the element -exposed to the sunlight or underneath a metal finger -the IV curve of the cell will be the light IV or the dark IV.Each element is linked to their neighbours with resistors, representing the lateral current flow and dependent on the sheet resistance of the cells.This method can be applied to any number of junctions.
This type of formalism is widely used to simulate the performance of solar cells when the effect of a spatial variable needs to be incorporated in the model.This variable can be the design of the front metal grid, in order to minimise the effect of series resistances [68]; the inhomogeneous illumination profile in concentrator devices; the impact of such inhomogeneity on the transport through the tunnel junctions [69,70]; or the distribution of defects or inhomogeneities [71,72].Recently, this formalism was used to model the photoluminescence and the electroluminescence based IV curves of MJ devices, accounting for the limited lateral carrier transport [73].
Specifically for the modelling and optimization of the front grid of solar cells in order to minimise shading losses and series resistance, there are two packages already available: PVMOS, developed by B. E. Pieters in C and released as open source [74,75], and Griddler, developed by J. Wong using Matlab and available at PV Lighthouse [11,76].

In-plane discretization
As shown in Fig. 16, there are two regions in the plane: the metal and the aperture.These two are provided to Solcore as grey scale images that will work as masks.The resolution of the images, in pixels, will define the in-plane discretization.By default, the aspect ratio of the pixels in the image will be 1:1, but this can be Fig.16: Schema of the quasi-3D solar cell modelling included in Solcore.The solar cell is discretized (a) in the plane and (b) in the vertical direction.Illuminated and dark regions are then modelled using electrical components that, when combined, form a 3D electrical mesh giving the voltages and currents at any point of the structure.(c) An example of the vertical discretization of a N-junction solar cell.set to a different value in order to reduce the number of elements and improve speed.For example, the in-plane discretization of Fig. 16a has an aspect ratio A r = L y /L x = 4, with L x and L y the pixel size in each direction.The values of the pixels in the metal mask are 0 where there is no metal (the aperture area), 255 where there is metal and the external electrical contacts (the boundaries with fixed, externally set voltage values) and any other value in between to represent regions with metal but not fixed voltage.The pixels of the illumination mask -which become the aperture mask after removing the areas shadowed by the metal -can have any value between 0 and 255.These values divided by 255 will indicate the intensity of the sunlight at that pixel relative to the maximum intensity.Fig. 17 illustrates two examples of metal masks (a and b) and an illumination mask (c) with 120×120 pixels.As it can be seen, while rectangular metal fingers are well reproduced, diagonal fingers are less accurate and could require a finer discretization.The illumination mask is mostly homogeneous except around the edges and in the corners, where intensity is much lower.This pattern could be produced, for example, by the secondary optics of a concentration system.
The minimum total number of nodes where SPICE will need to calculate the voltages will be N×M×2×Q, with N and M the number of pixels in both in-plane directions and Q the number of junctions, which require  and (b) are grey, indicating that there is metal in those pixels but that their bias is not set to be equal to the external bias.
2 nodes each.To this, the front and back metal contacts could add a maximum of 2M×M nodes.Exploiting symmetries of the problem as well as choosing an appropriate pixel aspect ratio will significantly reduce the number of nodes and therefore the time required for the computation of the problem.

Vertical discretization
First, the solar cell is solved as described in Section 5 and 6 in order to obtain the parameters for the 2-diode model at the given illumination conditions.These parameters are then used to replicate the 2-diode model in SPICE.The I SC is scaled in each pixel by the intensity of the illumination given by the illumination mask.Sheet resistances above and below each junction, R sh (top) and R sh (bot), account for the lateral transport.Beneath the metal, there is no current source, as the region is in the dark, and there are extra resistances accounting for the contact between the metal and the semiconductor R c and the transport along the metal finger R s [68].Given that the pixels can be asymmetric, these resistances need to be defined in both in-plane directions, x and y: R y sh = A r R sh (67) where h is the height of the metal, ρ m their linear resistivity and ρ c the contact resistivity between metal and semiconductor.The sheet resistance of a stack of semiconductor layers R sh is equal to the combination in parallel of the individual sheet resistances.Using the single junction example of Fig. 16, R sh (top) will be given by: 1 R sh (top) = 1 R sh (window) + 1 R sh (emitter) (71) Each of these can be estimated from the thickness of the layer d, the majority carrier mobility µ and the doping N as [70]: If the solar cell has been defined using only the DA and PDD junction models, this information is already available for all the layers of the structure.For junctions using the DB and two diode models, R sh will need to be provided for the top and bottom regions of each junction.Intrinsic layers will be ignored as they do not contribute to the lateral current transport.

Solar array model
The ability to use Solcore to build a SPICE equivalent circuit allows entire PV systems to be simulated from the bottom up [77].Each photovoltaic solar cell is described using an equivalent circuit which can then be arranged in strings of series and parallel cells to represent the entire system.An example for a triple junction solar cell, complete with a bypass diode is shown in figure 18; this unit is the basic building block for a concentrator PV module [52].The diode and resistance values for the equivalent circuit are determined from solar cell testing, while the current source is evaluated by integrating the product of the spectral irradiance (estimated using an appropriate radiative transfer code e.g.SPCTRL2 or SMARTS) and the quantum efficiency which in turn can be calculated dynamically as a function of temperature by Solcore [78].
Since the entire module (and subsequently the system) is assembled from individual solar cell components, it is possible (and indeed, necessary) to distribute the component values to accommodate for manufacturing tolerances.This enables a close match between the modelled output power and that measured experimentally and has been used to determine how both aerosols and precipitable water affect the electricity yield from concentrator PV systems [79,80].Where system IV data is available, the emergence of electrical faults, (e.g.shunts or shading) can also be accounted for [81].

Closing remarks
In this article we have described the capabilities of Solcore, a multi-scale, Python-based, modular simulation framework for semiconductor materials and solar cells.Its main strengths are: -Flexibility: Provides a variety of tools, rather than a single solution, for the study of traditional and novel semiconductor materials and devices.-Modularity: Can be expanded with new capabilities, innovative solvers and tools.
-Accessibility: Not only is it open source, but it is also designed to be easy to learn and to use, serving as a teaching tool as much as a research tool.-Rigour: The physics behind every functionality are well understood and supported by numerous references, as are the approximations made in order to simplify the implementation of the problem or the interpretation of the results.-Integrated: All of Solcore's features are designed to be compatible with one another to allow for truly multi-scale modelling in an integrated way.

Fig. 4 :
Fig. 4: Band structure of In 0.2 Ga 0.8 As calculated with the bulk k• p solver.The inset shows the effective mass determined for a range of indium fractions and a comparison with the experimental data from Volk et al.[25]

Fig. 8 :
Fig. 8: Plot of the spectra of different light sources.

Fig. 10 :
Fig.10: Example solutions from the TMM module in Solcore calculating Reflection, Transmission and complex dielectric function.(a) Spectral reflection (solid lines) and transmission (dashed lines) over a range of angles for the optimised triple layer ARC reported in[41].(b) The calculated real and imaginary parts of the complex dielectric function obtained from spectroscopic ellipsometry of a Ge substrate at 79 • (points), and the dielectric constant calculated using TMM for a model consisting of a semi-infinite Ge substrate and 4.5 nm GeO 2 (lines).

Fig. 11 :
Fig. 11: All the paths that radiation can follow withing the cell, used to calculate the absorptivity≡emissivity as a function of the angle.

Fig. 12 :
Fig. 12: (a) Work flow of Solcore's PDD solver.(b) Process of obtaining the effective band structure of QWs to use in the PDD solver.From left to right: simple sequence of layers; band profile and energy levels after considering the strain and quantum confinement; effective band structure.(c) Description of the inhomogeneous mesh scheme used in Solcore.

Fig. 13 :
Fig. 13: (a) Dark IV curve of a MJ solar cell, including the IV of the individual junctions separately (continuous lines) and the junctions as part of the MJ structure.(b) Light IV curve of the same MJ solar cell.

Fig. 14 :
Fig. 14: Light IV curve of a 3J solar cell (a) without and (b) with radiative coupling.Continuous lines represent the individual IV curves of the junctions isolated and the dash lines when they are inside the MJ device, illustrating the effect of the coupling.

Fig. 15 :
Fig. 15: IV curve of a tunnel junction defined according to the parametric model.

Fig. 17
Fig. 17: (a) and (b) two examples of metal masks and (c) an illumination mask.The thin metal fingers in (a)and (b) are grey, indicating that there is metal in those pixels but that their bias is not set to be equal to the external bias.
Fig. 17: (a) and (b) two examples of metal masks and (c) an illumination mask.The thin metal fingers in (a)and (b) are grey, indicating that there is metal in those pixels but that their bias is not set to be equal to the external bias.

Fig. 18 :
Fig. 18: Equivalent circuit for a triple junction solar cell.