A versatile open-source analysis of the limiting efficiency of photo electrochemical water-splitting

Understanding the fundamental thermodynamic limits of photo-electrochemical (PEC) water splitting is of great scientific and practical importance. In this work, a ‘detailed balance’ type model of solar quantum energy converters and non-linear circuit analysis is used to calculate the thermodynamic limiting efficiency of various configurations of PEC design. This model is released as freely accessible open-source (GNU GPL v3) code written in MATLAB with a graphical user interface (GUI). The capabilities of the model are demonstrated by simulating selected permutations of PEC design and results are validated against previous literature. This tool will enable solar fuel researchers to easily compare experimental results to theoretical limits to assess its realised performance using the GUI. Furthermore, the code itself is intended to be extendable and so can be modified to include non-ideal losses such as the over-potential required or complex optical phenomena.

Efficiency of series connected tandem PV system Figure S1. Maximum efficiency of series connected (optically and electrically) tandem photovoltaic system 2/9 Description of the model Current voltage relationship of a single photo-absorber There are many approaches to determining the limiting efficiency of a solar fuel process as described by Bolton et al. 1 . In this work, we approach the analysis from a 'detailed balance' type methodology pioneered by Shockely and Queiser 2 . Following the derivation from various sources [3][4][5] , it is possible to derive the theoretical current-voltage behaviour of a single bandgap absorber cell. Net current density is given by the detailed balance of the current from absorbed photons minus the current lost due to radiative recombination (eq. (1)). Photons absorbed are from the ambient surrounding (radiating blackbody radiation at a temperature of T a ) and from the sun (at a temperature of T s ). Radiative recombination is unavoidable and means such a conversion system cannot exceed the ideal Carnot efficiency limit (η carnot = 1 − T a /T s ). Analyses that omit this term will therefore be inaccurate as explained by Bolton et al. 1 .
The generalised form of Planck's radiation law, given in eq. (2), can be used to calculate the ideal blackbody photon emission flux b i , normal to the surface per energy interval dE. This generalised form includes the effect of the chemical potential ∆µ due to excitation separation on the radiated photon flux 5 . For an ideal solar cell where transport is lossless 6 , the chemical potential difference generated can be expressed as qV . For the sun and ambient surrounding, ∆µ = 0. The subscript i denotes the blackbody such as the sun (s), emission from the solar cell (e), the ambient surroundings (a). The geometric factor F i is equal to π for hemispherical acceptance or emission.
Assuming the ambient and the cell is in thermal equilibrium with its surroundings and that one photon excites only one electron (i.e no multiple carrier generation).
where b s (E) and b a (E) is the solar photon flux and ambient photon flux normal to the surface per energy interval dE. If the cell has an abrupt absorption threshold, and that the angular range of the sun is small in comparison to the ambient (F s /F e ≈ 0) we can write this as follows. At thermal equilibrium in the dark, it can be shown that b a (E, 0, T a ) must equal b e (E, 0, T a ) and henceforth the photon flux from ambient is substituted for the photon flux emitted by the solar cell in the dark (i.e ∆µ = 0). And therefore eq. (3) can be written as eq. (4), assuming that all excitations (from absorption of above bandgap photons) that survive radiative recombination are collected.
The radiative recombination current under illumination (so ∆µ = qV ) can be similarly written if we assume the cell is a selective blackbody, with an abrupt absorption threshold: And therefore, eq. (1) can be expressed as: Equation (7) shows an approximation that can be made for large enough bandgaps (> 0.2 eV) for the radiative recombination of the solar cell under illumination 7,8 . Also, here f g is introduced where it is the geometric factor that specifies whether the planar plate can radiate from both sides ( f g = 2) or just the front side ( f g = 1).
Therefore we can approximate the ideal current voltage relationship of the solar cell with ideal diode equation: where J o is given by: The solar spectrum of light can be approximated by a blackbody at a temperature of ∼6000 K. In this case, J L = q ∞ E g b s (E, 0, 6000K)dE and Planck's law, eq. (2), can be used to calculate b s . Alternatively, a standardised dataset could be used which is a more accurate representation of the terrestrial solar spectrum, such as ASTM G173-03 reference spectra -37 • tilted hemispherical (AM 1.5G) 9 and therefore: The photovoltaic efficiency can be defined as η PV = P PV P solar where P PV = J(V )V . The maximum efficiency will be found at the so called maximum power point at a particular voltage that satisfies d(P PV ) dV = 0.

Solar driven photolysis
If the photovoltaic device is coupled to an electrochemical process the free energy of the electrons can be stored in the chemical potential of a species produced: i.e a fuel such as H 2 .
The efficiency of this solar to fuel process can be expressed as: Where J(V ) is the current provide at some voltage V from the photo-absorber(s), E rxn is the thermodynamic potential (gibbs rather than thermoneutral), V = E rxn +V o . For a ideal cell, the overpotential V o will be equal to zero. In reality an overpotential is required due to energy lost to entropy in transport and kinetic processes. For solar water splitting, E rxn = E H 2 = 1.229 V at 300 K. The faradaic efficiency η F is the efficiency by which charge is transferred to the specific reaction of interest. Henceforth, it is assumed that η F = 1.
If the photo-system is driving multiple electrolysis systems that are electrically in series, then the solar to fuel efficiency is given by eq. (12) where i is the specific electrolyser and N elec is the total number of electrolysers.
In this work we will only consider electrolysers in series, hence for water splitting the solar to hydrogen efficiency can be written as:

Optical configuration
For a single photo-absorber, the light falling on the surface will be the complete input solar spectrum as defined in eq. (10). If there is more than one photo-absorber, consideration has to be made as to how they are optically configured, and Figure S2 shows the different configurations for up to 3 photo-absorbers. In order to model this, the configuration was represented as a directed acyclic graph made up of nodes and edges which connect nodes together. The first node is the input light spectrum (e.g AM1.5G) and each photo-absorber is placed at subsequent nodes and the edges represented the light path. The light before and after the photo-absorber at node i will be denoted by b − i and b + i respectively. The light at node 0 (input node) will therefore be b + 0 = b AM1.5G . For each connection between nodes, two quantities are specified. Firstly, the split fraction f i j from node i to j is defined in eq. (14) as the fraction of light reaching j versus the total amount of light leaving i. Here b i j is the photon flux from node i to node j.
Secondly, the area concentration between nodes is defined in eq. (15). This may differ from the 'light concentration' if the split fraction of that edge is less than unity.
In order to calculate the light at each node, the topological order of directed acyclic graph was computed as the graph is acyclic. Therefore before the light at a successor node was calculated, all the predecessor nodes had already been calculated. Light at a particular node i is given by the summation of the light from previous nodes, taking into account the transmission t i (E) of the photo-aborber at node i and the fraction of light from the previous node f i j .
The transmission function of the photo-absorber at i was assumed to be discontinuous corresponding to a sharp absorption edge.

5/9
In order to calculate the photocurrent of the photo-absorber corresponding to node i

Electrical configuration
From electrical connections between the components of the system (photo-absorbers and electrolysers), and the optical configuration of the photo-absorbers, we now can determine the currents and voltages within the circuit. Multiple electrolysers are considered in series and so the voltage required from the photo-absorber circuit, The circuit is represented as a directed cyclic graph where each branch is a component such as a PV or electrolyser. Unfortunately, as MATLAB does not permit multigraphs (multiple edges between nodes), we had to include a third type of branch, which was a simple electrical connection, in order to be able to represent all the permutations of design.
Due to the non-linear nature of the ideal diode equation used for the photo-absorbers, non-linear circuit analysis is required and the 'tableau matrix' methodology described in the book Linear and Nonlinear circuits by Chua, Desoer and Kuh 10 was used. Here, n is the number of nodes and b is the number of branches. φ i is the electric potential at node i. j j and v j are the branch currents and branch voltage respectively. The electrolyser connects the first (1) and last node (n) and is final branch (b).
1. Construct the directed graph of the circuit and choose the negative side of the electrolyser to be the datum node (φ = 0) 2. Formulate a reduced-incidence matrix, I I I, by removing the row associated with the datum node 3. Apply Kirchoff's current law to each node (except the datum node) (n − 1 equations)

Construct the branch equations (b equations)
• PV equations (b − 1 equations) where J PV is the current density of the photo-absorber calculated by eq. (8) and A PV is the illuminated area.
Non-linear root finding methods can be now used on the complete set of equations (KCL, KVL, branch eq.) in order to solve eq. (24).
As the numerical method used is an iterative one it requires an initial guess for j j j, v v v and φ φ φ . Through trial and error it was found that, whilst perhaps not optimal, the following initial guess appeared to be a relatively robust choice which found the solution in most cases. Furthermore, the fact that the initial guess is calculated rather than supplied by the user, serves to ensure the 'front end' graphical user interface is a simple as possible.
The initial guess for the voltage of the branches is given in eq. (25), and depends whether branch contains a photo-absorber, an electrolyser, or neither. The initial guess for the branch current j j j is given by eq. (26). Here, max( j mpp ) refers to the maximum current produced by any of the photo-absorbers operating at maximum power point.
The initial guess for the potential at each node, φ φ φ 0 , is calculated by finding the value of φ φ φ that minimises the residual when calculating Kirchhoff's Voltage Law (KVL) using a non-linear least-squares algorithm (see eq. (27) and eq. (21)).
Furthermore, in order to aid numerical solution, each quantity with a dimension was scaled before being passed to the numerical solver (fsolve). The scaling factors used were j * = max(J L ) and v * = e * = V electrolysis .