Quokka Manual
Quokka is a free software tool specialized for the fast simulation of silicon solar cell devices in one to three dimensions. It employs simplifications to the general semiconductor carrier transport model, which results in much less computational effort compared to alternative simulation software. Those simplifications, namely quasineutrality and conductive boundaries, are well validated to not impose notable loss of generality or accuracy for waferbased silicon solar cell devices. Thus Quokka enables to simulate even moderately complex 3D silicon solar cell device geometries in short computation times on standard computers, while providing a similar level of accuracy and generality as stateoftheart commercial device simulation software.
The main modelling difference to most other device simulation software comes with the conductive boundary simplification. Here nearsurface regions, for example an emitter diffusion, are not modelled in detail by defining the doping profile and surface recombination etc. The inputs required are rather lumped properties of such regions, most importantly the sheet resistance and the effective recombination characteristic like for example the emitter saturation current density J_{0e}. This is well suitable if those inputs are the ones derived e.g. experimentally, but makes Quokka not (directly) applicable if for example an optimization of the doping profile is of interest.
Quokka essentially solves for the steadystate electrical characteristics of the device, and is capable to derive various typical solar cell characteristics: fixed terminal voltage, fixed terminal current, opencircuit (OC) conditions, maximumpowerpoint (MPP) conditions, shortcircuitcurrent (Jsc) conditions, light and darkIV curve, quantumefficiency (QE) curve, sunsVoc curve and seriesresistance (Rs) curve. A notable extension is the inclusion of luminescence modelling, which allows Quokka to simulate for example spatially and spectrally resolved photoluminescence characteristics.
While being powerful enough to simulate many silicon solar cell designs and conditions of practical interest, Quokka’s scope is to be highly accessible for both newcomers and simulation experts. This is realized by the meshing and the solver numerics being largely automated and by using a predefined but flexible geometry layout, so that the required user inputs are reduced to a minimum and specifically designed for typical silicon solar cell conditions. Furthermore, Quokka automatically finds typical cell conditions of interest like maximum power point, and features an optimizer routine to fit unknown device properties to output characteristics. The software is hosted on pvlihgthouse.com.au, where Quokka can be downloaded for free and support resources as well as a webbased settings file generator can be found.
Although these measures enable to perform extensive simulation tasks with little knowledge about fundamentals of numerical simulation, it is reminded that Quokka still is a full multidimensional solver with a large space of possible input combinations. As with any such complex tool the principle of “garbage in, garbage out” is still highly applicable, and limits of solvable conditions exist. In particular users intending to produce sensitive simulation results, for example for scientific publications, are strongly encouraged to follow the considerations in Computational Performance and Accuracy.
Quokka is a compiled MATLAB program and requires the Matlab Comiler Runtime (MCR) 2013b to be installed. The MCR is freely available from the Mathworks website, be sure to download the correct (NOT the latest) version, 2013b:
The latest Quokka version is available for download as Windows 32bit / 64bit, Linux and OSX version:
The zip archive contains the Quokka executable as well as some example settings and input data files.
All settings required to define the simulation setup are given in the settings file. This is a Matlab script file (.m), which is however in ascii format and can therefore be modified with any text editor. While it is possible to directly edit the settings file, it is recommended to use the settings file generator on PVLighthouse, a webbased GUI which creates the settings file in a more userfriendly and less errorprone way:
Note that any input data files required for the simulation must be in the same directory as the settings file. Also all solution data will be stored in this directory.
The easiest way to run the simulation is using the little GUI, which opens by default when starting Quokka.
GUI Element  Description and comments 

Select different settings file 
Browse and select one or multiple settings files.
No reload is required when the settings file is modified. 
Display Geometry 
This is good practice to do before running the simulation as a quick sanity check for the geometry definition and resulting number of elements
For multiple settings files only the first one will be displayed Sweep and optimizer settings are ignored In the Figure which opens up (see Figure 2) one can rotate and zoom, and toggle the visibility of selected features: xyz axis, F front side, R rear side, M metal coverage and S symmetry sides For multiple conductive boundaries of the same type a higher sheet resistance will be of lighter colour than the lower one, which can be useful for validating the arrangement of the respective boundaries The latter must not be mistaken with shaded colours, which indicates a contacted region 
Start Simulation 
Starts the simulation, plots and saves results if finished successfully
Multiple settings files will be run sequentially 
Abort 
Aborts the simulation at selected points
Only limited functionality for parallel simulations 
Overwrite results files  When checked, Quokka will not append date and time to the result files 
Sweep: save spatial data  When checked, Quokka will save all spatial data for each parameter during a sweep 
Messages  Displays status and error messages 
When starting the Quokka executable from the command prompt, additional messages are displayed, and some additional functionality can be enabled by command line options, see Table 2.
For example, to run settings1.m and settings2.m parallelized on 5 CPUs in Windows:
Starting those settings file in Linux as a background job without graphical output and redirecting command line output to nohup.out:
Command line option  Description and comments 

*.m 
Settings file to be simulated, *.m must be replaced by an actual unique filename
Multiple settings files can be assigned and will be simulated sequentially; this saves the overhead to startup Quokka for each simulation 
sweep_store_spatial_data  Quokka will save all spatial data for each parameter during a sweep 
coresX 
Enables parallel simulation mode, where X needs to be replaced by the desired number of parallel threads (max. 12 in Matlab 2013b)
Only QEcurves, nonautomatic IVcurves, sweeps and sequential optimization are executed in parallel, and the parallel mode will consequently not increase the speed of other simulations Limited messaging and abort functionality 
no_display  Prevents the GUI and any graphic figures to open 
overwrite_results  Quokka will not append date and time to the result files 
sweep_store_spatial_data  Quokka will save all spatial data for each parameter during a sweep 
The input parameters are organized in several groups, which is loosely also reflected in the PVL settings file generator. These comprise geom (overall solution domain and contact geometry), bulk (bulk properties), bound (boundary conditions, i.e. diffusion and surface properties), generation (optics / generation), circuit (external circuit), lumi (luminescence modelling), sweep (parameter sweep) and optim (iterative optimization / curve fitting).
A standout input parameter is circuit.terminal, which controls at what operating point(s) the device will be simulated, and consequently what output characteristic is achieved. There is a qualitative difference between single operating point conditions, like a fixed terminal voltage, and “curves” consisting of many different operating conditions, like an IVcurve.
The sweep functionality sits “on top” and must be distinguished from the “curve” circuit.terminal options. If a curve terminal option is set together with a parameter sweep, the respective curve will be simulated for each sweep parameter.
The optimizer sits further “on top” on the sweep functionality, and the sweep can therefore not be used to run multiple optimizations with different input parameters. Such a “sequential optimization” must rather be defined using vector assignments to the optim.override inputs. The sweep can be used to create xy value pairs for curve fitting. However, most common targets are readily produced by the curve terminal options, for example an IV curve by “IV_curve” and effective lifetime curve by “sunsVoc_curve”.
Parameter name  Units  Typical value range  Description  Applicability / dependencies 

Geometry  
version.design  ‘FRC’ ‘IBC’  predefined geometry layout to use, see Predefined Layouts.  
geom.dimensions  1, 2 or 3  1D, 2D or 3D simulation  
geom.Wz  µm  50 … 500  device thickness (zdirection)  
geom.Wx  µm  1 … 10000  unit cell width in xdirection  2D and 3D only 
geom.Wxfront  µm  1 … 10000  front unit cell width in xdirection  if different, values must be integer and must have a reasonably low lcm, see Predefined Layouts. 
geom.Wxrear  µm  1 … 10000  rear unit cell width in xdirection  
geom.Wy  µm  1 … 10000  unit cell width in ydirection  3D only 
geom.frontcont.shape  ‘circle’ ‘rectangle’ ‘line’ ‘full’  shape of front (emitter) contact(s) in 2D, ‘circle’, ‘rectangle’ and ‘line’ produce identical geometry in 3D, ‘line’ will be orientated in ydirection  2D and 3D only  
geom.rearcont.shape  ‘circle’ ‘rectangle’ ‘line’ ‘full’  shape of rear contact(s) in 2D, ‘circle’, ‘rectangle’ and ‘line’ produce identical geometry in 3D, ‘line’ will be orientated depending on the defined contact position, see Predefined Layouts.  2D and 3D only  
geom.frontcont.wx  µm  1 … 10000  halfwidth of front contact(s) in xdirection radius for ‘circle’  2D and 3D only 
geom.frontcont.wy  µm  1 … 10000  halfwidth of front contact(s) in ydirection  3D only only for ‘rectangle’ contact shape 
geom.rearcont.wx  µm  1 … 10000  halfwidth of rear contact(s) radius for ‘circle’  2D and 3D only 
geom.rearcont.wy  µm  1 … 10000  halfwidth of front contact(s) in ydirection  3D only only for ‘rectangle’ contact shape 
geom.rearcont.position  [1] … [1 3] … [1 2 3 4]  vector with any combination of positions 1 to 4; defines the corner(s) where a rear contact exists, see Predefined Layouts; all contacts have same shape and dimensions.  2D and 3D only positions 3 and 4 in 3D only  
geom.leftcont.shape  ‘circle’ ‘rectangle’ ‘line’  shape of left (emitter) contact(s) in 2D, ‘circle’, ‘rectangle’ and ‘line’ produce identical geometry in 3D, ‘line’ will be orientated in ydirection  
geom.rightcont.shape  ‘circle’ ‘rectangle’ ‘line’  shape of right contact(s) in 2D, ‘circle’, ‘rectangle’ and ‘line’ produce identical geometry in 3D, ‘line’ will be orientated depending on the defined contact position, see Predefined Layouts.  
geom.lefttcont.wx  µm  1 … 10000  halfwidth of left contact(s) in xdirection radius for ‘circle’  
geom.leftcont.wy  µm  1 … 10000  halfwidth of left contact(s) in ydirection  3D only only for ‘rectangle’ contact shape 
geom.rightcont.wx  µm  1 … 10000  halfwidth of right contact(s) radius for ‘circle’  2D and 3D only 
geom.rightcont.wy  µm  1 … 10000  halfwidth of right contact(s) in ydirection  3D only only for ‘rectangle’ contact shape 
geom.leftcont.numberx  0.5, 1, 1.5 …  multiple of 0.5 defining the number of “half” left (emitter) contacts in xdirection within the unit cell, see Predefined Layouts.  
geom.leftcont.pitchx  µm  1 … 1000  full pitch of multiple left (emitter) contacts within the unit cell  only for numberx > 0.5 
geom.rightcont.numberx  0.5, 1, 1.5 …  multiple of 0.5 defining the number of “half”right contacts in xdirection within the unit cell, see Predefined Layouts.  
geom.rightcont.pitchx  µm  1 … 1000  full pitch of multiple right contacts within the unit cell  only for numberx > 0.5 
geom.leftcont.w_metal  µm  0 … 10000  left (emitter) metal half width, influences optics only  
geom.rightcont.w_metal  µm  0 … 10000  right metal half width, influences optics only  
geom.leftcont.y_position  ‘aligned’ ‘opposite’ ‘double’ ‘half’  position of left (emitter) contact(s) relative to right contact(s) in ydirection, see Predefined Layouts.  3D only  
geom.meshquality  1,2,3 or ‘user’  1: coarse (sufficient for most simulations), 2: medium, 3: fine, ‘user’: use expert settings below  
geom.dzminfront  µm  0.2 … 2  element size in zdirection at the front surface  for ‘user’ mesh quality only 
geom.dzminfront  µm  0.2 … 2  element size in zdirection at the rear surface  for ‘user’ mesh quality only 
geom.dxmin  µm  0.2 … 20  minimum element size in xdirection  2D and 3D only for ‘user’ mesh quality only 
geom.dymin  µm  0.2 … 20  minimum element size in ydirection  3D only for ‘user’ mesh quality only 
geom.dxmax  µm  1 … 500  maximum element size in xdirection  2D and 3D only for ‘user’ mesh quality only 
geom.dymax  µm  1 … 500  maximum element size in ydirection  3D only for ‘user’ mesh quality only 
geom.scale  3 … 20  determines minimum element sizes in x and ydirection via dividing the smallest respective features (contact dimensions, gap etc.) by this number  2D and 3D only for ‘user’ mesh quality only  
geom.inflation  1 … 5  maximum allowable ratio of neighbouring element sizes; effectively controls how fast mesh sizes are increased away from feature edges  for ‘user’ mesh quality only  
Bulk properties
Note that all recombination mechanisms are always active and additive, see Quokka Physics 

bulk.type  ‘ptype’ ‘ntype’  doping type  
bulk.rho  Ohm.cm  0.1 … 1000  bulk resistivity  
bulk.T  K  300K*  temperature; carefully check validity of applied models and inputs if changing temperature  
bulk.nk  ‘Green08_Nguyen14’* ‘Green08_300K’  silicon optical properties, used for 1D_model generation and luminescence modeling  
bulk.mobility  ‘Klaassen’* ‘Arora’ 'Fixed'  mobility model, see Quokka Physics  
bulk.mun bulk.mup 
cm²/Vs  constant electron / hole mobility  ‘fixed’ mobility only  
bulk.nieff  ‘default’* ‘fixed’  ‘default’: see Quokka Physics ‘fixed’: user value  
bulk.nieffvalue  cm^{3}  ~1e10  fixed value of nieff  for ‘fixed’ bulk.nieff only 
bulk.taub  µs  1 … 1e5  fixed “background” bulk lifetime set to very high value (e.g. 1e10) to disable  
bulk.Auger  ‘Richter2012’* ‘Altermatt2011’ ‘Kerr2002’ ‘Sinton1987’ ‘off’  Auger recombination model, see Quokka Physics  
bulk.Brad  cm^{3}/s  4.73e15*  radiative recombination coefficient set to 0 to disable radiative recombination  
bulk.SRH.BO.Nt  cm^{3}  0* … 1e18  oxygen concentration for BOcomplex recombination, see Quokka Physics. delete or set to zero to disable  for ptype bulk only 
bulk.SRH.BO.m  2 … 4  processing dependent parameter for BOcomplex recombination  for ptype bulk only  
bulk.SRH.midgap.taun0  µs  1 … 1e5  fundamental electron lifetime for midgap (defect level = Fermi level) SRH recombination  
bulk.SRH.midgap.taup0  µs  1 … 1e5  fundamental hole lifetime for midgap (defect level = Fermi level) SRH recombination  
bulk.SRH.custom{i}.Nt  cm^{3}  defect density of ith SRH defect  
bulk.SRH.custom{i}.Et_Ei  eV  0.5 … 0.5  defect level of ith SRH defect, relative to intrinsic energy  
bulk.SRH.custom{i}.sigman  cm²  electron capture cross section of ith SRH defect  
bulk.SRH.custom{i}.sigmap  cm²  hole capture cross section of ith SRH defect  
boundary conditions
Index {i} stands for the ith conductive / nonconductive boundary condition, see Defining Multiple Boundary Regions. Each boundary can have a contacted and a noncontacted region with different recombination properties assigned to .cont. and .noncont, respectively The IBC version does not support contacted nonconductive boundaries 

bound.conduct{i}.location  ‘front’ ‘rear’ ‘left’ ‘right’  location of the ith conductive boundary front and left located conductive boundaries are the “emitter” for FRC and IBC version, respectively  
bound.nonconducti{i}.location  ‘front’ ‘rear’  location of nonconductive boundary  
bound.conduct{i}.cont.rec bound.nonconduct{i}.cont.rec bound.conduct{i}.noncont.rec bound.nonconduct{i}.noncont.rec  ‘J0’ ‘S’ ‘expr’  How to model boundary recombination; by constant saturation current density J0, constant effective surface recombination velocity S or analytical expression  
bound.conduct{i}.cont.J0 bound.nonconduct{i}.cont.J0 bound.conduct{i}.noncont.J0 bound.nonconduct{i}.noncont.J0  A/cm²  0 … 1e10  saturation current density note the unit is A/cm², NOT fA/cm²  only for ‘J0’ recombination model 
bound.conduct{i}.cont.S bound.nonconduct{i}.cont.S bound.conduct{i}.noncont.S bound.nonconduct{i}.noncont.S  cm/s  0 … 1e6  effective surface recombination velocity 1e6 represents “infinite” recombination limited by carrier transport to the boundary; higher values therefore do not make any difference other than potentially causing numerical problems  only for ‘S’ recombination model 
bound.conduct{i}.cont.expr bound.nonconduct{i}.cont.expr bound.conduct{i}.noncont.expr bound.nonconduct{i}.noncont.expr  A/cm²  ‘expression’ 
analytical expression resulting in a boundary recombination current in A/cm²
Variables allowed to be used:
dn
: excess carrier density at the boundary in cm^{3}
const.N
: bulk doping density in cm^{3}
const.nieff
: bulk intrinsic carrier density in cm^{3 }
const.n0
: equilibrium minority carrier density in cm^{3 }
const.T
: temperature in K^{}
const.Vt
: thermal voltage in V^{ }
const.q
: elementary charge in As
e.g.: ‘10*dn*const.q‘
(corresponds to Seff = 10 cm/s) 
only for ‘expr’ recombination model 
bound.conduct{i}.cont.J02 bound.nonconduct{i}.cont.J02 bound.conduct{i}.noncont.J02 bound.nonconduct{i}.noncont.J02  A/cm²  0* … 1e6  n=2 saturation current density  only for ‘J0’ recombination model 
bound.conduct{i}.type  ‘ptype‘ ‘ntype‘  conduction type of front conductive boundary in IBC version  
bound.conduct{i}.Rsheet  Ω  1 … 10000  Sheet resistance  
bound.conduct{i}.shape  ‘full‘ ‘line‘ ‘rectangle‘ ‘circle‘ ‘contact‘  shape of conductive boundary is always aligned to any contact(s) of the same type ‘contact’ sets the shape identical to the same type contact defined in the geometry group  2D and 3D only not applicable for IBC front conductive boundary (always full area)  
bound.conduct{i}.wx  µm  1 … 10000  width of conductive boundary in xdirection; radius for ‘circle’  2D and 3D only 
bound.conduct{i}.wy  µm  1 … 10000  width of conductive boundary in ydirection  3D only ‘rectangle’ shape only 
bound.conduct{i}.jctdepth  µm  0* … 10  junction depth  
bound.conduct{i}.colleff  ‘fixed’ ‘ext_file’  how to model collection efficiency, see Quokka Physics  junction depth >0  
bound.conduct{i}.colleff_value  0 … 1*  fixed value for collection efficiency  ‘fixed’ collection efficiency only  
bound.conduct{i}.colleff_filename  ‘filename’  load collection efficiency from external file, see External File Format  ‘ext_file’ collection efficiency only  
generation  
generation.type  ‘1D_model’ ‘Jgen_surface’ ‘Jgen_uniform’ ‘ext_file’ ‘customdata’ ‘off’  how to derive generation profile ‘1D_model’: models generation profile by optical properties, see Quokka Physics ‘Jgen_surface’, ‘Jgen_uniform’: applies userdefined generation current to the illuminated surface or uniformly distributed through the device thickness  
generation.Jgen  mA/cm²  ~40  user defined total generation current  ‘Jgen_surface’, ‘Jgen_uniform’ generation or with Jgen_correction=1 
generation.Jgen_correction  0* 1  1: scales the external or ‘1D_model’ generation profile to match a desired total generation current defined by generation.Jgen  ‘1D_model’, ‘ext_file’ and ‘customdata’ generation only  
generation.ext_file  ‘filename’  load generation profile from external file, see External File Format  ‘ext_file’ generation only  
generation.customdata  µm cm^{3}s^{1}  [z1, z2, … zn; G1, G2, … Gn]  generation profile with n value pairs; z: distance to illuminated surface G: generation rate  ‘customdata’ generation only 
generation.spectrum  ‘AM1.5g‘ ‘monochromatic‘ ‘custom‘  defines the illumination spectrum  ‘1D_model’ generation only  
generation.monochromatic.wavelength  nm  250 … 1450  monochromatic illumination wavelength  ‘monochromatic’ spectrum only 
generation.monochromatic.flux  cm^{2}s^{1}  ~2e17  photon flux of monochromatic illumination  ‘monochromatic’ spectrum only 
generation.facet_angle  °  0 or 54.7  facet angle of illuminated surface texture; set to 0 for planar surface  ‘1D_model’ generation only 
generation.spectrum_custom  nm Wm^{2}nm^{1}  [λ1, λ2, … λn; I1, I2, … In]  custom spectrum with n value pairs: λ: wavelength I: spectral intensity  ‘custom’ spectrum only 
generation.transmission  ‘fixed’ ‘ext_file’ ‘custom’  how (wavelength dependent) transmission at the illuminated surface is defined ‘fixed’: fixed value for all wavelengths  ‘1D_model’ generation only  
generation.transmission_value  0 … 1  fixed value for front transmission  ‘fixed’ transmission only  
generation.transmission_filename  ‘filename’  load front transmission data from external file, see External File Format  ‘ext_file’ transmission only  
generation.transmission_custom  nm  [λ1, λ2, … λn; T1, T2, … Tn]  transmission curve with n value pairs: λ: wavelength T: transmission  ‘custom’ transmission only 
generation.Z  ‘fixed‘ ‘ext_file‘ ‘custom‘ ‘param‘ ‘limit_4n2‘ ‘limit_Green02‘  how (wavelength dependent) pathlength enhancement is defined  ‘1D_model’ generation only  
generation.Z_value  0 … 1  fixed value for pathlength enhancement  ‘fixed’ Z only  
generation.Z_filename  ‘filename’  load pathlength enhancement data from external file, see External File Format  ‘ext_file’ Z only  
generation.Z_custom  nm  [λ1, λ2, … λn; Z1, Z2, … Zn]  pathlength enhancement curve with n value pairs: λ: wavelength Z: pathlength enhancement  ‘custom’ Z only 
generation.Z0  1 … 1000  Z0 for Zparameterization  ‘param’ Z only  
generation.Zinf  1 … 3  Zinf for Zparameterization  ‘param’ Z only  
generation.Zp  1 … 10  Zp for Zparameterization  ‘param’ Z only  
generation.suns  0 ... 1 ... 10  scales the generation (AND intensity for efficiency calculation)  
generation.illum_side  ‘front’ ‘rear’  defines the illuminated side of the device  
generation.shading_width  µm  0 … 10000  half shading width by metal fingers, independent of the contact width; Is applied as a line shading in ydirection centred to each contact on the illuminated side  2D and 3D only 
generation.profile_type  ‘standard’ ‘cumulative’  ‘cumulative’ expects PC1Dtype profile  ‘ext_file’ and ‘customdata’ only  
external circuit  
circuit.terminal  ‘Vuc’ ‘Vterm’ ‘Jterm’ ‘OC’ ‘MPP’ ‘Jsc’ ‘light_IV_auto’ ‘IV_curve’ ‘QE_curve’ ‘sunsVoc_curve’ ‘Rs_curve’  what “terminal” condition to solve single operating point conditions: ‘Vuc’: fixed unit cell voltage (fastest because no iterations required) ‘Vterm’: fixed terminal voltage ‘Jterm’: fixed terminal current ‘OC’: opencircuit (do NOT set high series / contact resistance as in PC1D) ‘MPP’: maximum power point ‘Jsc’: shortcircuitcurrent (NOT shortcircuit, voltage will be >0, but extracted current will be Jsc) “curve” terminal condition: ‘light_IV_auto’: automated algorithm to quickly but accurately derive light IVcurve and parameters (Voc, Jsc, FF, eta) ‘IV_curve’: solve at multiple user defined voltages ‘QE_curve’: quantum efficiency at multiple user defined wavelengths ‘sunsVoc_curve’: solve at multiple user defined sunsvalues ‘Rs_curve’: extract voltage dependent series resistance via doublelight method; solves two automated lightIV curves with slightly different generation; advice: set circuit.IV_accuracy to at least 5 for smooth and accurate results  
circuit.Vuc.value  V  0.3 … 1  unit cell voltage (BEFORE external circuit elements); avoid values close to zero, see Shortcircuit Limitations  ‘Vuc’ terminal only 
circuit.Vterm.value  V  0.3 … 1  terminal voltage (AFTER external circuit elements); avoid values close to zero, see Shortcircuit Limitations  ‘Vterm’ terminal only 
circuit.Jterm.value  mA/cm²  40 … 1000  terminal current  ‘Jterm’ terminal only 
circuit.IV.mode  ‘Vterm’ ‘Vuc’  whether to define unit cell (faster) or terminal voltages (slower) for IV_curve solution  ‘IV_curve’ terminal only  
circuit.IV.V_values  V  [V1, V2, … Vn]  n voltage values for IV_curve  ‘IV_curve’ terminal only 
circuit.IV.init_previous  0* 1  1: use solution from previous results as starting guess; faster and more stable for closely spaced IV points; disables parallel simulation of IVcurve 0: initialize start values for each IV point; more robust for largely spaced IV points  ‘IV_curve’ terminal only  
circuit.QE.wavelength_values  nm  [λ1, λ2, … λn]  n wavelength values for QE_curve  ‘QE_curve’ terminal only 
circuit.sunsVoc.suns_values  [s1, s2, … sn]  n suns values for sunsVoc_curve  ‘sunsVoc_curve’ terminal only  
circuit.DJ0  A/cm²  0*… 1e7  saturation current density of external diode, see Quokka Physics  requires circuit.Dn 
circuit.Dn  1 … 4  ideality factor of external diode, see Quokka Physics  required if circuit.DJ0 is defined  
circuit.Voc_guess  V  0.67*  guess of Voc for iteration starting point; better guess can speed up the simulation  ‘OC’, ‘MPP’ and ‘light_IV_auto’ terminal only 
circuit.IV_accuracy  1* … 10  higher value increases the number of points simulated on the light IVcurve does NOT influence the accuracy of IVcurve parameters (Voc, Jsc, FF, eta) >= 5 is recommended for ‘Rs_curve’ terminal option  ‘light_IV_auto’ and ‘Rs_curve’ terminal only  
luminescence modelling
For luminescence modelling “front” always refers to the detection side which could actually be the rear side 

lumi.enable  1 0  switch luminescence on (1) or off (0) if off, no further inputs are required  
lumi.scale  a.u.  simply scales the simulated luminescence signal can be used to convert simulated detected photon flux into measurement units, e.g. counts/s  
lumi.normalize_signal  1 0  whether to normalize (1) the simulated signal or not (0); will normalize the spectrum and / or intensity map to its respective peak value  
lumi.detection_side  ‘illuminated’ ‘opposite’  whether the sensor detects luminescence from the illuminated or nonilluminated side  
lumi.filter  ‘none’ ‘ext_file’  how to define wavelength dependent transmission from the sample surface to the sensor (i.e. lens and filter transmission) ‘none’: 100% transmission  
lumi.filter_filename  ‘filename’  load optics transmission data from external file, see External File Format  ‘ext_file’ filter only  
lumi.sensor_EQE  ‘unity’ ‘silicon’ ‘ext_file’  how to define wavelength dependent external quantum efficiency of sensor ‘unity’: 100% EQE ‘silicon’: calculates EQE assuming LambertBeerLaw for absorption with silicon absorption coefficient and sensor thickness  
lumi.sensor_thickness  µm  thickness of sensor for calculation of EQE; more correctly, this is the total pathlength which is directly inserted into the LamberBeerLaw  ‘silicon’ sensor EQE only  
lumi.sensor_EQE_filename  ‘filename’  load sensor EQE from external file, see External File Format  ‘ext_file’ sensor EQE only  
lumi.wavelength_start  nm  800* … 1400  lower bound of wavelength to include into modeling  
lumi.wavelength_end  nm  800 … 1400*  upper bound of wavelength to include into modeling  
lumi.emission_function  1 2 3 4  how to model internal optics, see Quokka Physics 1: both sides planar, quasi1D 2: both sides textured, quasi1D 3: general, quasi1D; further inputs required, see below; set of default values approximately applicable for typical textured front and semiplanar rear solar cell devices 4: both sides lambertian, 2D and 3D geometrical blurring; userdefined mesh settings with relatively low dxmax (and dymax) recommended for smooth and accurate results; can be highly computationally intensive  4: 2D and 3D only  
lumi.internal_refl  ‘default’ ‘ext_file’  how to define wavelength dependent internal surface reflectivity ‘default’: assumes planar silicon – air interface on both sides  not applicable for emission function 2  
lumi.internal_refl_filename  ‘filename’  load internal reflectivity for front and rear surface from external file, see External File Format  ‘ext_file’ internal reflectivity only  
lumi.general.facet_angle  °  0 … 54.7  facet angle of front surface texture set to zero for planar surface  emission function 3 only 
lumi.general.lambertian_factor  0 … 1  lambertian factor of rear side 0: completely specular 1: completely diffuse  emission function 3 only  
lumi.general.refRear  0 … 0.6* … 1  internal rear reflectivity  emission function 3 only  
lumi.general.refFrontS  0 … 0.62* … 1  first internal front reflectivity  emission function 3 only  
lumi.general.refFrontN  0 … 0.93* … 1  nth internal front reflectivity  emission function 3 only  
lumi.geometrical.number_reflections  0 … 20  number of internal reflections to trace high number might add significant computational effort  emission function 4 only  
parameter sweep
Index {i} stands for the ith dependent sweep parameter All dependent sweep parameters must have the same number of values To have only one independent (group of) sweep parameter(s), assign values only to one (group) (…_1 or …_2) 

sweep.enable  1 0  switch parameter sweep on (1) or off (0) if off, no further inputs are required  
sweep.param_1{i} sweep.param_2{i}  ‘input parameter‘  input parameter to sweep within the first / second independent group of sweep parameters, e.g.: ‘bulk.rho‘ must exactly match an input parameter existing within the settings file  
sweep.values_1{i} sweep.values_2{i}  same as ith sweep parameter  [v1, v2, … vn] {‘S1‘, ‘S2‘, … ‘Sn‘} 
sweep values for ith sweep parameter of the first / second independent group of sweep parameters
can be numeric (v) or strings (S)
sweeping of some string values including external file names NOT possible yet
to disable the second (group of) independent sweep parameter(s):
sweep.value_2{1}=[]; 

optimizer
Index {i} stands for the ith parameter to optimize / the ith goal to achieve / the ith override rule to apply For sequential optimization, each vector input must have the same length n 

optim.enable  2 1 0  switch off / change mode of optimizer 2: curve fitting mode 1: goal mode 0: off; not further input required  
optim.maxiter  5 … 50  maximum number of iterations; aborts once this number is reached  
optim.accuracy  0.1 … 1* … 10  influence termination tolerances higher values potentially increases accuracy but at the cost of increased number of iterations and increased likelihood of “stuck” optimizations  
optim.param{i}  ‘input parameter‘  input parameter to optimize, e.g.: ‘bulk.rho‘ must exactly match an input parameter existing within the settings file  
optim.start{i}  same as ith optim. parameter  same as ith optim. parameter  start value of parameter to optimize should be carefully chosen not far from the final result for good convergence  
optim.lb{i}  same as ith optim. parameter  same as ith optim. parameter  lower bound of allowed values for parameter to optimize <optim.start{i}  
optim.ub{i}  same as ith optim. parameter  same as ith optim. parameter  upper bound of allowed values for parameter to optimize >optim.start{i}  
optim.goal{i}  ‘scalar result‘ e.g. ‘results.eta’  name of ith goal to achieve; must be a scalar within the results structure produced by the simulation setup, see Table 5  mode 1 only  
optim.value{i}  same as scalar result defined in ith goal  ‘max’ ‘min’ v [v1, v2, … vn]  value of ith goal to achieve; ‘max’: maximises the result ‘min’: minimises the result v: (vector of) numerical value(s)  mode 1 only ‘min’ and ‘max’ only for a single goal i=1 
optim.goalx 
‘vector result’
e.g.
‘results.curves
.navg’ 
name of xvalues to achieve must be a vector within the results structure produced by the simulation setup, see Table 5  mode 2 only  
optim.goaly 
‘vector result’
e.g.
‘results.curves
.taueff’ 
name of yvalues to achieve must be a vector within the results structure produced by the simulation setup with the same length as the one defined in optim.goalx, see Table 5  mode 2 only  
optim.xyfile  ‘filename’  external file containing xyvalue pairs which define the target curve to fit to  mode 2 only  
optim.override.param{i}  ‘input parameter’  name of ith parameter to override during the optimization  
optim.override.value{i}  ‘expression’ v [v1, v2, … vn]  override value to assign during the optimization ‘expression’: analytical expression containing the parameter(s) to optimize, e.g.: ‘bulk.rho*5’ v: (vector of) numerical value(s) 
Several inputs like for example a generation profile can be defined via external files. Those files must be located in the same directory as the settings file and the containing data must be correctly organized to be usable by Quokka.
Accepted file formats are either Excel files (.xls or .xlsx) or ASCII files (.txt or .csv). Note that only “proper” Excel files can be read, while those produced with alternative tools, like for example by the PVLighthouse website, can cause reading errors. This can be fixed by opening and saving them in Excel. For the ASCII files, most common delimiters are supported. Furthermore, a comma decimal separator is supported, which takes slightly longer to import due to internal string replacement.
The data must be organized in columns, where the first two (three) columns contain the input data and the first row can optionally be a header row or the first data row. Table 4 summarizes the data organization for the different external file inputs.
Assigned to input parameter  First column  Second column  Third column 

generation.ext_file  distance to illuminated surface in µm  generation rate in cm3s1  
generation.transmission_filename  wavelength in nm  front transmission 0 … 1  
generation.Z_filename  wavelength in nm  pathlength enhancement Z  
bound.conduct{i}.colleff_filename  wavelength in nm  collection efficiency 0 … 1  
lumi.filter_filename  wavelength in nm  transmission 0 … 1  
lumi.sensor_EQE_filename  wavelength in nm  EQE 0 … 1  
lumi.internal_refl_filename  wavelength in nm  front reflectivity 0 … 1  rear reflectivity 0 … 1 
optim.xyfile  same quantity and unit as defined in optim.goalx  same quantity and unit as defined in optim.goaly 
Naturally, what output data is available depends on what characteristics the device is solved for. A basic distinction is between single operating point and multiple operating points characteristics. Only for a single operating point a unique electrical condition with unique spatial results, including luminescence results, is present and will be stored. This is not the case for the curve terminal conditions ‘IV_curve’, ‘sunsVoc_curve’ and ‘Rs_curve’, and consequently no spatial and luminescence output is available. Outstanding operating points exist for ‘light_IV_auto’ (maximum power point) and ‘QE_curve’ (short circuit current conditions) and therefore spatial data is available for those, as well as for all single operating point terminal conditions ‘Vuc’, ‘Vterm’, ‘Jterm’, ‘OC’, ‘MPP’ and ‘Jsc’.
When a sweep is performed, no spatial data will be available by default, however can be stored for each parameter by activating the corresponding switch, see Running A Simulation.
The kind of graphical output is predefined and dependent on what terminal condition was solved, and whether a sweep or optimization task was performed. Spatial data will only be plotted when a single point terminal condition was set, respective curves will be shown for curve terminal conditions, and the variation of few major scalar outputs for a parameter sweep. Beyond that, there is no influence by the user on the graphical output. All important results are stored in a csv and Matlab file, the latter containing also spatial output if applicable, from which users must create custom graphics themselves.
Quokka saves two results files, a .csv file containing most scalar and curve outputs, and a .mat file containing additional results including spatial output data (if available). The .mat file is a Matlab file format, which can be read by several software tools, but is naturally best processed in Matlab. If the optimizer is used, the results will be saved in files which include ‘optim’ in the file name (with some additional outputs not listed below).
Output quantity  Units  Description 

Scalar outputs  
Vterm  mV  terminal voltage (including external circuit elements) 
Jterm  mA/cm²  terminal current density (including external circuit elements) 
Vuc  mV  unit cell voltage (excluding external circuit elements) 
Juc  mA/cm²  unit cell current density (excluding external circuit elements) 
Voc  mV  open circuit voltage (@ 1 sun for sunsVoc curve) 
Jsc  mA/cm²  shortcircuit current density (total generation current density @ 1 sun for sunsVoc curve) 
FF  %  fill factor of lightIV curve (pseudo fill factor for sunsVoc curve) 
eta  %  efficiency (pseudo efficiency for sunsVoc curve) 
Vmpp  mV  terminal voltage at maximum power point 
Jmpp  mA/cm²  terminal current density at maximum power point 
navg  cm3  average minority carrier density 
taueff  µs  effective lifetime, definition see Quokka Physics 
conduct  S  excess conductivity, definition see Quokka Physics 
JL_two_diode  mA/cm²  current density source from twodiode model fit 
J01_two_diode  A/cm²  ideal saturation current density from twodiode model fit 
J02_two_diode  A/cm²  nonideal (n=2) saturation current density from twodiode model fit 
Rseries_two_diode  Ωcm²  series resistance from twodiode model fit 
Rshunt_two_diode  Ωcm²  shunt resistance from twodiode model fit 
Rseries_mpp_pow  Ωcm²  series resistance @ MPP from power loss, definition see Quokka Physics 
Rseries_mpp_DLM  Ωcm²  series resistance @ MPP from double light method 
lumi_mean  a.u.  average luminescence signal (detected photon flux times scaling factor) 
computing_time  s  total computing time (not correct when optimizer is used) 
JGbulk  mA/cm²  generation current density of the bulk 
JGfront  mA/cm²  generation current density of the front conductive boundaries 
JGfront  mA/cm²  generation current density of the rear conductive boundaries 
N  cm3  bulk doping density 
nieff  cm3  effective intrinsic carrier density 
n0  cm3  equilibrium minority carrier density 
FELA outputs (see Quokka Physics for definitions)  
FELA.gen  mW/cm²  generated power 
FELA.res_bulk_n  mW/cm²  resistive loss from electron transport in the bulk 
FELA.res_bulk_p  mW/cm²  resistive loss from hole transport in the bulk 
FELA.res_front_conduct FELA.res_left_conduct  mW/cm²  resistive loss in the front / left conductive boundary(ies) 
FELA.res_front_cont FELA.res_left_cont  mW/cm²  resistive loss from front / left contact resistance 
FELA.res_rear_conduct FELA.res_right_conduct  mW/cm²  resistive loss in the rear / right conductive boundary(ies) 
FELA.res_rear_cont FELA.res_right_cont  mW/cm²  resistive loss from rear / right contact resistance 
FELA.res_front  mW/cm²  IBC only: resistive loss in the front conductive boundary 
FELA.rec_Auger  mW/cm²  intrinsic bulk recombination loss (Auger and radiative) 
FELA.rec_SRH  mW/cm²  other bulk recombination loss (SRH and fixed bulk lifetime contributions) 
FELA.rec_front_conduct FELA.rec_left_conduct  mW/cm²  recombination in the front /left conductive boundary(ies), includes loss due to nonoptimal collection efficiency 
FELA.rec_front_cont FELA.rec_left_cont  mW/cm²  recombination at the front / left contacts 
FELA.rec_rear_conduct FELA.rec_right_conduct  mW/cm²  recombination in the rear / right conductive boundary(ies), includes loss due to nonoptimal collection efficiency 
FELA.rec_rear_cont FELA.rec_right_cont  mW/cm²  recombination at the rear / right contacts 
FELA.rec_front_nonconduct  mW/cm²  recombination in the front nonconductive boundary 
FELA.rec_rear_nonconduct  mW/cm²  recombination in the rear nonconductive boundary 
FELA.ext_series  mW/cm²  loss at the external series resistance 
FELA.ext_shunt  mW/cm²  loss at the external shunt resistance 
FELA.ext_diode  mW/cm²  loss at the external diode 
FELA.error  %  FELA balance error, see Verify Accuracy 
curve outputs  
curves.Vterm  mV  terminal voltage (including external circuit elements) 
curves.Jterm  mA/cm²  terminal current density (including external circuit elements) 
curves.Vuc  mV  unit cell voltage (excluding external circuit elements) 
curves.Juc  mA/cm²  unit cell current density (excluding external circuit elements) 
curves.QE_lambda  nm  wavelength values for QE curve 
curves.CE  %  collection efficiency, definition see Quokka Physics 
curves.IQE  %  internal quantum efficiency 
curves.EQE  %  external quantum efficiency 
curves.suns  suns  values for sunsVoc curve 
curves.Voc  mV  open circuit voltage 
curves.navg  cm3  average minority carrier density 
curves.taueff  µs  effective lifetime, definition see Quokka Physics 
curves.Vi  mV  terminal voltage, finely spaced values for Rscurve 
curves.Rseries_DLM  Ωcm²  series resistance from double light method 
main spatial outputs  
grid.XX / YY / ZZ  m  vector with coordinate values of grid 
grid.X / Y / Z  m  3D array with coordinate values of grid 
grid.dx / dy / dz  m  3D array with elements sizes 
spatial.dn  cm3  3D array of excess carrier density 
spatial.nfront / nrear  cm3  2D array of excess carrier density at the front / rear boundary 
spatial.QFn / QFp  V  3D array of electron / hole quasiFermi potential 
spatial.Pel  V  3D array of electric potential 
spatial.Ffront / Frear  V  2D array of potential in front / rear conductive boundary(ies) 
spatial.G  cm3s1  3D array of bulk generation rate 
spatial.R  cm3s1  3D array of bulk recombination rate (sum of all contributions) 
spatial.Jnx / Jny / Jnz mA/cm²  3D array of electron current density component  These are defined at element faces and therefore contain one more value in each coordinate direction than other 3D arrays 
spatial.Jpx / Jpy / Jpz mA/cm²  3D array of hole current density component  These are defined at element faces and therefore contain one more value in each coordinate direction than other 3D arrays 
spatial.sigman / sigmap  S/cm  3D array of electron / hole conductivity 
spatial.Jrec3  mA/cm²  3D array of recombination current density at solution domain boundaries (contains nonzero values at outer elements only) 
Quokka’s (relative) ease of setup, speed and capability for extensive sweep and optimizing tasks was found to raise expectations on stability and convergence compared to more complex tools. It is therefore reminded that this is still a full multidimensional simulation tool allowing for a wide range of conditions, so convergence naturally is a challenging topic. Yet Quokka, due to its simplifications and specialization, actually is more robust and faster for the conditions it is designed for than full device simulators.
Quokka’s numeric are designed and optimized for typical solar cell conditions, and bad or no convergence can be expected for very badly performing cell structures, exotic test structures or extreme conditions. This should especially be considered when using large parameter sweeps which may cover such extreme conditions. Quokka’s philosophy is to keep settings simple and don’t allow users “messing around” too much (potentially producing wrong results), therefore there is limited possibility for the user tweak numerics.
In general, finer mesh quality gives higher accuracy at the cost of computational time. It must be considered however, that a finer mesh in general also worsens convergence. Usually this results in a crash, but potentially also leads to a solution with low accuracy.
While meshing in Quokka is largely automated and even the coarse mesh setting will yield acceptable accuracy in many cases, verifying numerical accuracy is ultimately the responsibility of the user. The comments below give hints on how to judge and reach numerical accuracy:
Conditions with illumination and close to short circuit are challenging because of the physically existing large gradients. in particular the relatively coarse mesh close to the surface compared to other device simulators, as a result of the conductive boundary simplification, may lead to significant inaccuracies and />or convergence issues in Quokka. A finer mesh increases accuracy but can also worsen the convergence. It is strongly advised to avoid voltages close to zero, i.e. Vuc < ~200 mV. For representing shortcircuit conditions it is usually sufficient to use higher voltages which still extract Jsc. This is what Quokka does for automatic lightIV and Jsc terminal conditions, and is an important consideration when performing a manual lightIV curve.
A typical numerical artefact related to this shortcircuit inaccuracy is an apparent minimum of the current density at around 100 mV – 300 mV of the lightIV curve, with the minimum value being closest to the actual Jsc value.
The statistical escape function models 1 – 3 come with negligible computational effort relative to the device simulation. For an accurate spectrum, a finer zresolution of the mesh may be required, which is up to the user to check and verify.
The geometrical blurring model 4 can be computationally more expensive than the electrical simulation, in particular in 3D. Significant inaccuracy in spatial results can be present if large lateral element sizes are set, observable by unphysical peaks and nonsmooth results. To overcome this, the use of expert mesh settings with small geom.dxmax and geom.dymax is recommended. Another source of inaccuracy is present when the typical “blurring length” is in the order of the unit cell width. More than one repetition of the unit cell would be necessary, which would however be computationally expensive and is so far not implemented. Also note that this blurring model has not yet been extensively validated and results should be interpreted and used with care.
Using the optimizer functionality in Quokka means iteratively performing full numerical device modelling within a large parameter space. As with any such complex and potentially strongly nonlinear function to optimize, it is of high importance to carefully setup the optimization task. This comprises choosing sensible start values for the parameter(s) to optimize, tight upper and lower bounds to it, and having the confidence that the optimal value is within the bounds and that the target is sufficiently sensitive to the parameter(s) to optimize. Furthermore, the number of parameters and goals, though in principle unlimited, should be limited to a minimum, to prevent bad convergence and / or finding of local rather than global minima. For example, while choosing more than two or three input parameters as unknown’s to maximise efficiency is easily set up, it will hardly converge to the true global maximum, unless the start values are already very close to the solution. For curve fitting, care needs to be taken by the user to ensure overlapping of the userdefined and simulated xvalues within the given bounds of the unknown parameter, as leastsquares calculation can only be performed in the overlapping region.
As mentioned, Quokkas scope is to only provide a minimum number of input parameters to define and solve typical wafer based silicon solar cell devices. Being the main task for most multidimensional device simulations, Quokka is restricted to solve an element of symmetry or 'unit cell', which means that the sides (2 sides in 2D, 4 in 3D) are always set to symmetry boundary conditions. Furthermore, coarsely predefined geometry layouts are used, with the flexibility to change relevant dimensions and contact patterns. However, this means that Quokka's input capabilities do not allow the definition of arbitrary device geometries. The two supported layouts are called 'front and rear contact' (FRC) and 'interdigitated back contact' (IBC) version.
The FRC version fixes the front to only have contacts and conductive boundary types of opposite polarity to the bulk, and at the rear of the same polarity, respectively. In other words, all emitter diffusions and emitter contacts must be located at the front, and all base contacts on the rear. Rear emitter cells can still be simulated by setting illumination to the rear. The FRC version features multiple rear contacts located at any of the unit cell corners / edges, which allows for example to define a hexagonal rear contact pattern. Furthermore it supports a different unit cell width in xdirection for the front and the rear side. This effectively means that one can define a different front and rear contact pitch. This results not in a valid unit cell geometry, which Quokka handles by repeatedly stitching the front and rear geometry until they exactly match up. This larger valid unit cell has an xdimension which is equal to the least common multiplier (LCM) of the defined front and rear width. This means that no arbitrary numbers should be chosen which could result in extremely large solution domains, but care has to be taken by the user to ensure a reasonably low LCM.
The IBC versions has both polarity contacts on the rear, and no contacts on the front. It is organized in such way that the opposite polarity conductive boundaries and contacts ('emitter') are always positioned at the left side, and the respective base features on the right side. Within the IBC unit cell one can define an arbitrary number of in xdirection equally spaced and shaped contacts for the left and right side independently. For 3D no such feature is available in ydirection, here it can only be defined whether the contacts are positioned at one or both of the edges.
Quokka allows to define an arbitrary number of conductive boundaries with different shapes and properties by indexing. In areas where conductive boundaries overlap, the one with the highest index is applied. Their centre position is always aligned to every contact at the respective side (front / rear for FRC, left / right for IBC). An exception exists in the IBC version: if the width of a conductive boundary is larger than the xpitch of the respective contacts, it is centred at the respective edge of the solution domain (left / right). Adjacent conductive boundaries of the same type are electrically connected. In regions where no conductive boundary is applied, properties for the nonconductive boundary are applied. Figure 5 illustrates this functionality.
Quokka uses an orthogonal, nonequidistant mesh for discretization of the solution domain, which is well suited for the typically cuboidal shapes of silicon solar cell devices. The automated meshing algorithm works as follows: first, a minimum element size in each coordinate direction is determined from the minimum of either the user input or the minimum feature / gap size divided by a userdefined scale factor. This minimum element size is applied to all feature edges, including front and rear side, but NOT the symmetry sides. From there the mesh is inflated where the size of adjacent elements is increased by the inflation factor. Typical settings for scale factor, inflation factor and minimum element sizes are predefined in the meshquality settings 1 (coarse), 2 (medium) and 3 (fine), and can further be defined by the user. Note that for many purposes a 'coarse' mesh will yield acceptably accurate results, and is the recommended initial setting due to numerical speed and robustness. However, mesh independency should be verified for sensitive results, see Verify Accuracy.
A complication arises when different junction depths for conductive boundaries are defined. The physical solution domain in Quokka covers the quasineutral region only, i.e. excludes the junction depths. At the same time the cuboidal solution domain can only accept a unique value for its size in zdirection. Quokka handles this by setting the solution domain size to the defined thickness less the areaaverage junction depth, see Figure 6. Care is taken to use the originally defined junction depths when calculating the generated and collected current in the boundary, and compensating for the spatial gap / overlap by accordingly adjusting the generation in the bulk next to the boundary. This approach is accurate if the junction depths are much smaller than the device thickness and the size of the conductive boundaries.
The solver algorithm in Quokka requires a voltage to be applied, which makes the definition of both polarity contacts a fundamental requirement . Often it is desirable to simulate test structures without contacts, i.e. lifetime or photoluminescence test structures. This can still be done in Quokka by defining 'virtual contacts' and solving for open circuit conditions. If the properties of the virtual contacts are defined in a way that they represent noncontacted regions, their influence at open circuit eventually vanishes and the simulation represents a sample without any contacts. This is achieved by setting recombination of the contacted region equal to the noncontacted region it represents, and choosing a small contact size to not significantly supress lateral gradients by pinning the respective quasiFermipotential to the constant voltage of the boundary. However, the virtual contact should not be too small to prevent extremely small element sizes and consequently high computational effort. A typically good value for the size of the virtual contact is a few percent of the solution domain size.
Quokka has the builtin functionality to perform parameter sweeps. It can sweep one or two independent groups of dependent parameters indexed. Here, dependent parameters mean that both must have the same number ofvalues, which is the number of simulations to be performed using the defined value pairs, and are indexed in the settings file by curly brackets {i}. For the two independent groups of parameters in contrast, each value of the first parameter group will be combined with each value of the second parameter group. For example, if each of the parameters in the first independent group has 5 values assigned sweep.values_1{1}=[v1 v2 v3 v4 v5], and 4 in the second independent group sweep.values_2{1}=[v6 v7 v8 v9], a total of 5x4=20 simulations will be performed. Either of the two independent groups can be disabled by assigning an empty matrix to its values, e.g. for disabling the second group: sweep.values_2{1}=[].
Sweep parameter can be any parameter assigned in the settings file, with the exception of some parameters with string values, in particular external filenames.
If the parallel mode is in, a sweep will be simulated in parallel, see Command Line Options.
Inbuilt in Quokka is a generic optimizer algorithm, essentially using Matlab′s bounded nonlinear optimizer algorithm′s, which allows to fit unknown input parameter(s) to known output characteristics. For a successful optimization, a careful setup is required as with any other such complex function to optimize, see Optimizer.
There are two different optimizer modes. Mode 1 allows to define multiple goal values to achieve, or one goal to maximize / minimize. For a minimize / maximize goal the number of unknown parameters can be more than one. Otherwise the number of unknown parameters have to be equal or less than the number of goals, where in the latter case leastsquare fitting is performed rather than attempting exact goal achievement. Mode 2 performs curve fitting, where the simulated curve will be fitted by least squares to a userdefined curve. This is most commonly done in conjunction with curve terminal conditions to produce the desired curve (IV, effective lifetime curve by sunsVoc, to), while the sweep functionality can also be used to produce the desired curve.
Another functionality within the optimizer is given by the 'override parameters'. This is needed when other input parameters are dependent on the unknown parameter, and need therefore be adjusted during the iterative optimization. For example, if one wanted to fit front and rear J0 of a symmetric lifetime structure, one would define only the front J0 as the unknown parameter, and override the rear J0 with the front J0 to ensure a symmetric structure for each simulation.
A second application for those overrides is to envoke 'sequential optimization' to perform multiple identical optimization tasks with different conditions. In the above example, one might want to investigate the influence of the wafer thicknesses on the derived J0. For this, the override can be used to assign multiple values to any input parameter other than the unknown one, and Quokka will perform the defined optimization for each of those values. Alternatively / additionally start, bound and goal inputs can be assigned multiple values (with the same number) to perform sequential optimization. Note that sequential optimization can NOT be envoked by the sweep functionality, as the optimizer is implemented 'on top'. Such a sequential optimization is performed in parallel, when the parallel mode is on, see Command Line Options.
Table 6: Common problems with suggested solutions and / or explanations
Problem / question  Solution / description 

Inconsistent n_{i,eff} and J_{0} 
While (bulk) n_{i,eff} does actually have little direct impact on most important results like IVcharacteristics, it can indirectly show large influence in combination with “J_{0}” boundary recombination. The user must take care that a consistent assumption for n_{i,eff} for both deriving J_{0} and for the subsequent device simulation is used. J_{0} values can often be adjusted to a different n_{i,eff} by the assumption of J_{0}⁄n_{i,eff}^{2} = const, well valid for diffusions.
As the most prominent example, the Excel spreadsheets on Sinton lifetime testers typically assume n_{i,eff} = 8.6e^9 cm^{3} when calculating J_{0} values, which results in ca. 30% underestimation of recombination when the same J_{0} values and default n_{i,eff} settings in Quokka are used at 300K. 
Quokka stuck / bad or no convergence 
In estimated 95% of cases this is caused by a mistake in the settings or by setting very untypical conditions for a solar cell, not being reported by Quokka as such. Below is a list of common cases leading to convergence problems. If you still think Quokka should be able to solve your setup, please contact support@pvlighthouse.com.au.

Cryptic error message  Not being the fun part and lacking commitment within a free tool, the implementation of error handling is far from comprehensive and a “cryptic error message” is the dominant behaviour when something goes wrong. Mostly this is due to errors in the settings file. If you created / changed it manually, a convenient way for a sanity check is uploading it to the settings file generator, which might find and fix wrong or missing assignments. If you can’t find anything wrong with your settings file, please send it to support@pvlighthouse.com.au. 
Should I use ‘S’ or ‘J_{0}’ model for surface recombination?  Which model is more appropriate is up to the user to decide and depends on the detailed physics within the nearsurface region. As a simple guide, most highly doped regions or highly charged surfaces follow a J_{0} recombination, while some passivated surfaces are well approximated by an effective surface recombination velocity. In low injection, both models yield equivalent results. For low and moderately charged surfaces, neither of the two models may be applicable, and an analytic expression for the actual injection dependence needs to be used. 
What exactly do the inputs suns and intensity do? 
Once a generation profile is defined by either of the generation types, it is scaled by the user input for suns, which therefore effectively states a scaling factor. However, further influence is given on the postcalculated efficiency as it also scales the intensity.
The user input for intensity does have no influence on the generation. It is solely required to postcalculate the efficiency. Furthermore, the intensity user input will be disregarded in the ‘1D_model’, as it is then defined by the chosen spectrum. The formula below highlights the influence of suns and intensity on the efficiency: efficiency = terminal power / (intensity × suns) A complication arises when one wants to scale an external generation profile to fit a desired current. Using the sunsvalue does this job, but also scales the intensity for calculating the efficiency, which is consequently not the desired (but physical correct) value. Therefore suns should not be used for this purpose, but rather generation.Jgen_correction should be enabled with the desired generation current defined as generation.Jgen. 
© 2014 Andreas Fell
Quokka numerically solves the 1D/2D/3D steadystate charge carrier transport in a quasineutral silicon device in an efficient and fast manner. It uses so called 'conductive boundaries' to account for increased lateral conductivity in near surface regions (e.g. diffusion or inversion layer) and thus to simulate most silicon solar cell devices without major loss of generality. This page summarizes the implemented models and material properties.
The equations here are shown for a ptype bulk as an example, and can be analogously derived and applied for an ntype bulk.
The basic equations solved in Quokka are a simplified set of two differential equations describing the steadystate charge carrier transport in a quasineutral semiconductor. The variables to solve for are the quasiFermi potentials of electrons and holes φ_{Fn} and φ_{Fp} respectively.
(1)The conductivities and the recombination rate in general is a function of the carrier density. They can be derived from the quasiFermi potentials by the quasineutrality condition p ≈ n + N_{A}  n_{0} and equation (18) to [1]:
(3)For both conductive and nonconductive boundaries, the recombination current at or into the boundary is evaluated by either a J_{01}/J_{02} model, S model, or a userdefined analytic expression.
(4)At the symmetry boundaries, i.e. all boundaries except the front and the rear side, the component of the gradient of the quasiFermi potentials in direction of the surface normal, defined by the normal vector , is zero:
(7)The electron current equals the negative recombination current:
(8)At noncontacted regions the total current into the boundary needs to be zero:
(9)An electrical contact with the voltage V at the nonconductive boundaries is included by setting the current density through the contact resistance r_{c} equal to the total current density at the boundary:
(10)At regions without external contacts the contact resistance becomes infinite and sets the total current density into the boundary consistently to zero. Care should be taken when simulating a lowly recombining contact, as the quasineutrality approximation is not valid when significant currents are extracted and can thus result in significant errors.
Note that this implementation of the contact resistance, also in case of a conductive boundary, accounts correctly for current transfer effects into the contacts.
The conduction of carriers in a conductive layer at the surface (e.g. diffusion or inversion layer) is accounted for by a 2D continuity equation of the majority quasiFermi potential φ_{Fdiff} with the sheet resistance R_{sheet}.
(11)For a ntype conductive boundary (emitter) this gives the boundary condition for the bulk electron quasiFermi potential:
(12)The source term J_{diff} equals the total current density from the base into the emitter less the current density into the contact:
(13)The hole current density into the boundary equals the recombination current density less the current density generated within the boundary J_{G} times its collection efficiency η_{coll}, which gives the boundary condition for the hole quasiFermi potential:
(14)For a ptype diffusion (back surface field) the expressions are equivalent with altered carrier types and an inverted sign of the recombination current density.
The effective band gap energy is calculated by a temperature dependent E_{g,0} and a doping dependent band gap narrowing (BGN) ΔE_{g}.
(15)ΔE_{g} is determined by a lookup table with values derived from [2] which considers a doping dependence at 300 K only, i.e. representing the curve in Fig. 9 of [2].
To calculate carrier densities Boltzman statistics are assumed. This is valid as Quokka solves only in the relatively lowly doped bulk, where the difference to FermiDirac statistics is negligible.
(16)Here φ_{n} and φ_{p} are the electron and hole quasiFermi potentials, N_{c} and N_{v} the density of states (DOS) in the conduction and valence band respectively, and V_{T} = kT/q the thermal voltage. Note that what is denoted as the electric potential φ_{e} is more accurately the conduction band edge potential, which differs to the electric potential by the electron affinity and the reference electric potential.
It follows the wellknown relation between pn product and Fermilevel splitting which is applied to solve the bulk carrier transport, with the effective intrinsic carrier density n_{i,eff}^{2} being the only decisive value for accurate modelling of carrier transport.
(18)As default settings Quokka uses the values for DOS published by Green [3], and the temperature dependent band gap expression given in the Sentaurus user manual with slightly changed parameters to match the commonly accepted best value for the intrinsic carrier density at 300 K, n_{i,0}=9.65 × 10^{9} cm^{3}.
(20)The conductivities σ of the carriers are functions of their mobilities μ.
(23)Implemented in Quokka are the mobility models by Klaassen [4] and Arora [5], where in the latter case the excess carrier density is added to the doping density to account for injection dependence, which is the way it is implemented in PC1D.
Quokka accounts for radiative, Auger and ShockleyReadHall (SRH) recombination in the bulk to calculate the recombination rate R. The sum of radiative and Auger recombination is called intrinsic recombination. The Auger models implemented in Quokka comprise the parameterizations by Kerr and Cuevas [6], Altermatt [7] and Richter et al. [8] (eq. (18) without B_{rel}), as well as the constant recombination coefficient model by Sinton and Swanson [9], which is the one used in the spreadsheet of the Sinton testers.
Implemented in Quokka is the general SRH expression:
(25)Multiple custom defects can be defined by their defect level E_{t} relative to the intrinsic energy E_{i} (approx. defined as the energy level in the middle of the effective bandgap), the capture cross sections c_{n/p} and the defect density N_{t}, with the thermal velocity v_{th} = 1.1 × 10^{7} cm s^{1}.
(26)Additionally a simplified midgap defect energy recombination is implemented where the user inputs τ_{p0}/τ_{n0}, and n_{1}/p_{1} are equal to n_{i,eff}.
As a further contribution is given by the boronoxygen recombination in ptype material as given by Altermatt [7], where the user needs to set the oxygen concentration N_{t,0} and a processing dependent parameter m while the boron concentration is given by the doping density N_{A} and the energy level is 0.41 eV below the conduction band edge.
(29)The last contribution to bulk recombination is by a fixed, i.e. injection independent, minority carrier lifetime τ_{b,fixed}. The total bulk recombination is then calculated as the sum of all contributions, where each contribution can be either switched off, or needs to be assigned very high lifetime values to effectively disable:
(31)Wavelength dependent refractive index / absorption data is required for determining the generation profile when the ‘1D_model’ is used, as well as for reabsorption within luminescence modelling. By default, Quokka uses the optical properties of silicon from Nguyen [10], and combines it with the data from Green [11] to cover the full wavelength range. Quokka incorporates the published temperature dependence for deriving optical properties at arbitrary temperatures. As a further option, which is prominently used in literature and other software tools, the data from Green [11] fixed at 300K can be set.
Quokka’s focus in on the simulation of carrier transport and excludes the capability to perform optical modelling from physical optical properties of the cell and consequently derive the required generation profile. There are several options to define generation, comprising the direct definition of the generation profile by the user, a userdefined generation current, and a onedimensional model accounting for wavelengthdependent effective properties, which was published in [12].
The onedimensional model ‘1D_model’ derives a generation profile from a defined spectrum, front transmission T_{front}(λ) and pathlength enhancement Z(λ). Note that those inputs already define the total generation current (for a given wafer thickness W_{z} and refractive index data of silicon). Consequently total generation can be considered a userinput rather than a quantity modelled by Quokka, which ensures consistency with for example raytracing programs when deriving the required inputs. Furthermore, the facet angle of the front surface texture must be given, which is zero for a planar surface.
For each wavelength bin (5nm in Quokka), the incoming photon flux N_{ph}(λ) is derived from the incident spectrum. Each monochromatic ray is then traced through the wafer accounting for its angle to enhance its travel length and absorption. For a textured surface, a different angle is applied inside and outside of the junction depth as suggested in [13], see figure below. After traveling to the rear of the wafer, the remaining absorption from the difference of total pathlength and firstpass pathlength is equally distributed throughout the thickness of the device. It has been found that this model provides high accuracy for waferbased silicon solar cells.
The free energy loss analysis (FELA) as comprehensively described in [14] is implemented in Quokka. It enables to compare resistive losses and recombination losses in the common power loss unit mW/cm². Incidentally, for the typical 100 mW/cm² incident light intensity, the power loss values are equivalent to efficiency loss in %abs.
A drawback of the FELA is that recombination losses are in general underestimated, as they are a function of quasiFermi level splitting, which in turn is a result of recombination. The losses add up to the generated power density, which is also a function of quasiFermi level splitting and therefore overall recombination. Effectively, the part of the entropy loss which is caused by avoidable recombination is missing in FELA when aiming for a loss breakdown towards the AugerLimit [15]. However, FELA is very suitable to identify the main resistive and recombination loss contributions and the potential for incremental efficiency improvements. Table 1 summarizes the expressions used for the different loss contributions considered in Quokka.
A notable quantity Quokka calculates is the FELA balance ε_{FELA}. Mathematically the output power must equal the generated power less the sum of the FELA contributions. The FELA balance gives a metric how close this is numerically achieved, and is consequently a useful indicator for the consistency and accuracy of the numerical solution to the carrier transport equations. See the user manual for details on how this should be interpreted. It is defined as
(32)FELA name / comment  Equation 

generated power  
bulk recombination loss  
boundary recombination loss including suboptimal collection efficiency loss for conductive boundaries 

bulk electron resistive loss  
bulk hole resistive loss  
conductive boundary resistive loss note that boundary current density has units of A/m, opposed to A/m² in the bulk 

contact resistance loss takes current transfer effects into account 

external shunt resistance loss  V_{uc}^{2} / R_{shunt} 
external series resistance loss  J_{term}^{2}R_{series} 
external diode loss  V_{uc}J_{D} 
Quokka 2.0 simulates the detected luminescence map by integrating over the z axis of the simulated device:
(33)where r_{sp} is the rate of spontaneous emission via bandtoband transitions in silicon, f_{esc} the escape function defining the probability that a photon will be emitted from the device and Φ the detector sensitivity. As Quokka separates the electrical and optical simulations, luminescence maps can be produced under either photoluminescence or electroluminescence conditions. Each component of Equation ( 31 ) is described in detail in the subsequent sections.
The rate of spontaneous emission via bandtoband transitions is calculated using the generalized form of Planck's Law as described by Wurfel [16]. The rate of spontaneous emission (in units of photons nm^{1} s^{1} cm^{2}) emitted isotropically into 4π steradians is expressed as:
(34)where α_{Si}(λ) is the wavelengthdependent absorption coefficient of silicon, BB(λ,T) the blackbody photon flux and Δη(x,y,z) = φ_{Fn}  φ_{Fp} the local separation of the quasiFermi energies. Quokka employs the optical data of Daub over the wavelength range of the emission spectrum of silicon as it has demonstrated the best fit to experiment [17]. The electronic simulation of Quokka calculates the local splitting of the quasiFermi potential for every element of the device and from Equation ( 32 ), an accompanying radiative emission rate, an example of which is produced in the figure below.
This spectrum, however, is not what is ultimately emitted from the silicon wafer or solar cell, which is determined by the optical losses within the structure.
Numerous socalled emission functions exist in literature. Several [1820] have been implemented in Quokka. These functions describe the probability with which a spontaneously emitted photon will escape from the emitting medium without being 'lost', either to reabsorption or total internal reflection.
In the case of perfectly planar devices, only a relatively small proportion of emitted photons fall within the narrow escapecone that is determined by the relative refractive indices of the emitting and ambient media. In Quokka, it is assumed that the device is in an air ambient (n_{air}(λ) = 1), and the fraction 1/n_{Si}^{2} of total emitted photons falls within the escape cone.
If the device surfaces are perfectly planar and parallel, the photons incident on the siliconair interface beyond the escape cone are totally internally reflected and assumed to be unavailable to an external detector. However, not all photons within the escape cone will be transmitted on their first interaction with the wafer surface and are instead reflected internally. Owing to subsequent internal reflections these photons will have a further opportunity to be emitted from the front surface of the device, but on each pass will be attenuated by reabsorption in the silicon bulk modeled after BeerLambert absorption. This brings us to the escape function proposed by Schick [18], which is implemented in Quokka to describe the emitted luminescence flux for planar devices.
(35)It has been noted by Schinke et al. [20] that the escape function of Equation (33) is functionally equivalent to those described by Green [21], Daub [22] and Trupke [23], and thus is representative of a function broadly used in the PV literature.
Quokka also includes escape functions that account for diffuse and/or specularly reflecting surfaces. The model of Rudiger et al. [19] assumes that the front and rear surfaces of the device are fully diffuse reflectors, as might be encountered with front and rear textured devices. The escape function is described by:
(36)The figure below plots the emitted luminescence spectrum simulated by Quokka for the escape functions of Schick and Rudiger.
A more general escape function has been developed by Schinke et al. [20], which describes the probability of photon emission for any combination of planar (specular) and textured (diffuse) reflecting surfaces. It contains the functions of Schick and Rudiger as special cases, determined by appropriate input variables. It has been implemented in Quokka, with the functions of Schick and Rudiger retained for the relative ease of use. For details on the general emission function, refer to [20], which also provides a valuable overview of escape functions employed in the PV literature.
It is important to note that the above escape functions are statistical predictions of emission, and do not account for the spatial distribution of the emitted photons. Thus, luminescence maps, particularly in the case of diffusely reflecting samples, need to be interpreted with caution. Quokka assumes that the emitted photons are constrained within the zaxis of the element from which they are emitted.
To address this last mentioned limitation, the geometrical blurring model by Padilla et al. [24] is implemented in Quokka and can be chosen instead of the escape functions. It assumes both the front and the rear internal surfaces to be Lambertian reflectors. Note that it requires a relatively fine lateral mesh to give accurate and smooth results, which, in particular in 3D, can be highly computationally expensive.
The detection function describes the wavelengthdependent external quantum efficiency (EQE) of the imaging sensor, including the presence of any short or longpass filters. The response is modelled after BeerLambert absorption assuming perfect collection efficiency for absorbed photons. Quokka will determine the QE for a siliconbased detector of userdefined thickness, or will accept a definition of the sensor's EQE via an external excel spreadsheet.
Below is a sketch of the external circuit elements and current directions. The convention is that a positive unit cell voltage V_{uc} in dark condition results in a positive unit cell current J_{uc}, regardless of whether it is a ntype or ptype solar cell, and V_{uc} is applied to the contact of opposite polarity to the bulk. The external diode is in forward bias and its current is given by the saturation current density J_{0D} and the ideality factor n_{D}.
(37)Quokka has an inbuilt algorithm to derive a biased quantum efficiency curve. First, using the bias illumination only, short circuit current conditions are searched for iteratively. Note that this does not mean actual short circuit conditions, but searching for a higher unit cell voltage which still results in a negligible difference of extracted current to short circuit current J_{sc,bias}, which is usually the case up to a couple of 100 mV. The reason for this is to prevent numerical inaccuracies, which can be significant at actual short circuit. This so derived unit cell voltage is used as a fast to solve for fixed terminal condition for the following simulations. A monochromatic illumination is then added with a photon flux representing a small fraction of the bias generation current J_{ph}=0.005 J_{gen,bias} and a sweep of the wavelength results in a wavelength dependent short circuit current J_{sc}(λ).
The fundamental quantum efficiency characteristic Quokka simulates is the collection efficiency CE(λ), which is largely independent of the optical properties of the device. The other more common characteristics internal and external quantum efficiency IQE(λ) and EQE(λ) are essentially derived from the collection efficiency by the optical input parameters front surface transmission T_{front} and pathlength enhancement Z(λ). The respective formulas are given below, with the effective transmission T_{eff}(λ) being defined as the (unshaded) front transmission times the unshaded area fraction T_{eff}(λ) = T_{front}(λ) A_{unshaded}⁄A_{total}:
(38)[1]  A. Fell, "A free and fast threedimensional/twodeimensional solar cell simulator featuring conductive boundary and quasineutrality approximations," IEEE Transactions on Electron Devices, 60, pp. 733–738, 2013. 
[2]  A. Schenk, "Finitetemperature full randomphase approximation of band gap narrowing for silicon device simulation," Journal of Applied Physics, 84, pp. 36843695, 1998. 
[3]  M.A. Green, "Intrinsic concentration, effective densities of states, and effective mass in silicon," Journal of Applied Physics, 67, pp. 29442954, 1990. 
[4]  D.B.M. Klaassen, "A unified mobility model for device simulation  I. Model equations and concentration dependence," SolidState Electronics, 35, pp. 953959, 1992. 
[5]  N.D. Arora, J.R. Hauser and D.J. Roulston, "Electron and hole mobilities in silicon as a function of concentration and temperature," IEEE Transactions on Electron Devices, 29, pp. 292295, 1982. 
[6]  M.J. Kerr and A. Cuevas, "General parameterization of Auger recombination in crystalline silicon," Journal of Applied Physics, 91, pp. 24732480, 2002. 
[7]  P.P. Altermatt, "Models for numerical device simulations of crystalline silicon solar cells  a review," Journal of Computational Electronics, 10, pp. 314330, 2011. 
[8]  A. Richter, S.W. Glunz, F. Werner, J.Schmidt and A. Cuevas, "Improved quantitative description of Auger recombination in crystalline silicon," Physical review B, 86, 165202, 2012. 
[9]  R.A. Sinton and R.M. Swanson, "Recombination in highly injected silicon," IEEE Transations on Electron Devices, 34, pp. 13801389, 1987. 
[10]  H. T. Nguyen, F. E. Rougieux, B. Mitchell, and D. Macdonald, "Temperature dependence of the bandband absorption coefficient in crystalline silicon from photoluminescence," Journal of Applied Physics115, pp. , 2014. 
[11]  M.A. Green, "Selfconsistent optical parameters of intrinsic silicon at 300 K including temperature coefficients," Solar Energy Materials and Solar Cells, 92, pp. 13051310, 2008. 
[12]  A. Fell and K. R. McIntosh, "Deriving the generation profile for silicon solar cells from lumped optical parameters," in 42rd IEEE Photovoltaic Specialists Conference, New Orleans, USA, 2015. 
[13]  P. A. Basore, "Numerical modeling of textured silicon solar cells using PC1D," Electron Devices, IEEE Transactions on, vol. 37, pp. 337343, 1990. 
[14]  R. Brendel, S. Dreissigacker, N. P. Harder, and P. P. Altermatt, "Theory of analyzing free energy losses in solar cells," Applied Physics Letters, vol. 93, pp. 1735031735033, 2008. 
[15]  J. Greulich, H. Höffler, U. Würfel, and S. Rein, "Numerical power balance and free energy loss analysis for solar cells including optical, thermodynamic, and electrical aspects," Journal of Applied Physics, vol. 114, pp. , 2013. 
[16]  P. Wurfel, S.Finkbeiner and E. Daub, "Generalized Planck's radiation law for luminescence via indirect transitions," Applied Physics A, 60, pp. 6770, 1995. 
[17]  B. Mitchell, M.K. Juhl, M.A. Green and T. Trupke, "Full spectrim photoluminescence lifetime analysis on silicon bricks," IEEE Journal of Photovoltaics, 3, pp. 962969, 2013. 
[18]  K. Schick, E. Daub, S. Finkbeiner and P. Wurfel, "Verification of a generalized Planck law for luminescence radiation from silicon solar cells," Applied Physics A, 54, pp. 962969, 2013. 
[19]  M. Rudiger, T. Trupke, P. Wurfel, T. Roth and S.W. Glunz, "Influence of photon reabsorption on temperature dependent quasisteadystate photoluminescence lifetime measurements on crystalline silicon," Applied Physics Letters, 92, 222112, 2008. 
[20]  C. Schinke, D. Hinken, J. Schmidt, K. Bothe and R. Brendel, "Modeling the spectral luminescence emission of silicon solar cells and wafers," IEEE Journal of Photovoltaics, 3, pp. 10381052, 2013. 
[21]  M.A. Green, "Analytical expressions for spectral composition of band photoluminescence from silicon wafers and bricks," Applied Physics Letters, 99, 131112, 2011. 
[22]  E. Daub and P. Wurfel, "Ultralow values of the absorption coefficient for bandband transitions in moderately doped Si obtained from luminescence," Journal of Applied Physics, 80, pp. 53255331, 1996. 
[23]  T. Trupke, E. Daub and P. Wurfel, "Absorptivity of silicon solar cells obtained from luminescence," Solar Energy Materials and Solar Cells, 53, pp. 103114, 1998. 
[24]  M. Padilla, H. Höffler, C. Reichel, H. Chu, J. Greulich, S. Rein, W. Warta, M. Hermle, and M. C. Schubert, "Surface recombination parameters of interdigitatedbackcontact silicon solar cells obtained by modeling luminescence images," Solar Energy Materials and Solar Cells, 120, pp. 363–375, 2014. 