Possible influence of the Kuramoto length in a photo-catalytic water splitting reaction revealed by Poisson--Nernst--Planck equations involving ionization in a weak electrolyte
aa r X i v : . [ c ond - m a t . s o f t ] J a n Possible influence of the Kuramoto length in aphoto-catalytic water splitting reaction revealed byPoisson–Nernst–Planck equations involving ionization ina weak electrolyte
Yohichi Suzuki, Kazuhiko Seki ∗ Nanomaterials Research Institute(NMRI), National Institute of Advanced IndustrialScience and Technology (AIST)AIST Tsukuba Central 5, Higashi 1-1-1, Tsukuba, Ibaraki, Japan, 305-8565
Abstract
We studied ion concentration profiles and the charge density gradient causedby electrode reactions in weak electrolytes by using the Poisson–Nernst–Planck equations without assuming charge neutrality. In weak electrolytes,only a small fraction of molecules is ionized in bulk. Ion concentration profilesdepend on not only ion transport but also the ionization of molecules. Weconsidered the ionization of molecules and ion association in weak electrolytesand obtained analytical expressions for ion densities, electrostatic potentialprofiles, and ion currents. We found the case that the total ion densitygradient was given by the Kuramoto length which characterized the distanceover which an ion diffuses before association. The charge density gradientis characterized by the Debye length for 1:1 weak electrolytes. We discussthe role of these length scales for efficient water splitting reactions usingphoto-electrocatalytic electrodes.
Keywords:
Photo-catalytic water splitting reaction, Debye length,Kuramoto length, Poisson–Nernst–Planck equations ∗ Corresponding author
Email address: [email protected] (Kazuhiko Seki)
Preprint submitted to Elsevier August 2, 2018 . Introduction
The Poisson-Nernst-Planck (PNP) equations have been used to describea wide range of transport phenomena from electrons and holes in semicon-ductors to ions in electrolytes. [1–9] The conservation of charge for each ionspecies and electrostatic interactions among charge carriers are treated bythe PNP equations in a self-consistent manner. The PNP equations takeinto account the drift currents due to the electric fields generated by thedistribution of charge carriers. The effect is substantial e.g. when electrodereactions generate charge density gradients. The PNP equations can be ap-plied to obtain concentration profiles and an electro-static potential gener-ated by electrode reactions. Although the PNP equations are used for thestudy of coupled effects between electric fields and charge carrier transports,they are nonlinear and solved mainly by numerical methods. However, forcertain cases, approximate analytical solutions have also been obtained forstrong electrolytes. [2, 3, 8–13] Using the PNP equations, it has been shownthat the spatial dimensions of concentration gradient can be many ordersof magnitude larger than the characteristic length scale of charge densityprofiles given by the Debye length in strong electrolytes. [1–4, 13]In weak electrolytes, only a small fraction of molecules are ionized inbulk. Ion concentration profiles depend on ion transport and the ionizationof molecules. It should be noted that depleted ions in the vicinity of theelectrode can be replenished by the ionization of molecules and diffusionfrom the bulk in weak electrolytes. The current at the interface between theelectrolyte and electrode can be much larger than that in bulk due to themolecule ionization.The Nernst-Planck equation is the continuity equation representing theconservation of charge for each ion species, where dissociation and associationof ions are not considered. For weak electrolytes, the PNP equations havebeen extended to describe ionization and recombination of ions. [5, 14–25]In the extended PNP (e-PNP) equations, ion concentrations become non-conservative by local dissociation and association of ions in bulk phase.In this paper, we study charge transport induced by electrode reactionsin weak electrolytes using the e-PNP equations without assuming a prioricharge neutrality. Photo-electrochemical (PEC) conversion of water can bean example, but the fundamental results apply to other electrode reactions.We show that the gradient of total ion density (the sum of cation and anionconcentrations) is characterized by either the Kuramoto length [26–31] or the2
Anions
Nuetral Product
Electrode
Cations x Nuetral Product of cation reduction
Electrode
Cations x Nuetral Product of anion oxidation (a) (b)
Figure 1: Sketch of the model of ion discharge reactions at the electrode surface. Thedistance from the electrode is denoted by x . (a) Cations are reduced at the electrodesurface, and the neutral products are dissolved into the solution. Anions are reflected atthe electrode surface. (b) Cations are reduced, and anions are oxidized at the electrodesurface. The both neutral products are removed from the electrode surface. For photo-electrocatalytic (PEC) water splitting using a single photocatalyst, cations and the neutralproducts correspond to H O + and H , respectively; anions and the neutral products cor-respond to OH − and O , respectively. Debye length depending on the situation, while the gradient of charge density(the difference between cation and anion concentrations) is characterized bythe Debye length in binary monovalent weak electrolytes. The Kuramotolength characterizes the length scale of local density fluctuations around auniform concentration state. [26–31] In weak electrolytes, the Kuramotolength is given by the length scale of diffusive migration of ions within itslife-time; the life-time is determined by the association rate of ions. [26–31]Our results indicate that the ion density gradients can also be character-ized by the Kuramoto length when both cations and anions are dischargedat the electrode in binary monovalent weak electrolytes. For this case, iondensity drop caused by electrode reactions is recovered by ion density fluc-tuations localized within the Kuramoto length. We discuss the efficiency ofwater splitting reactions using photo-electrocatalytic electrodes in terms ofthe Kuramoto length and the Debye length. We also discuss the overpoten-tial related to a charge density gradient near the electrode and show that itcan be reduced when both cations and anions are discharged at the electrode.3 . Theory
As shown in Fig. 1, we considered the thermal motion of cation and anionunder Coulombic interactions. The one-dimensional (1D) x -coordinate wasintroduced by assuming a uniform distribution of ions in the plane parallelto the electrode surface. The origin of the x -coordinate was the electrodesurface, and the x -coordinate was normal to the electrode surface. The con-centration and diffusion constant of cations are denoted as n + and D + , re-spectively. The concentration and diffusion constant of anions are denoted as n − and D − , respectively, while the electric field is denoted as E . Each speciesmoves by electromigration and diffusion. For simplicity, we consider a binarymonovalent electrolyte (1:1 electrolyte). The concentration and the diffusionconstant of the undissociated neutral compound are denoted as m and D m ,respectively. At a sufficiently large distance L away from the electrode, theelectrolyte is neutral. The external electric field is not applied. We con-sider ion discharge reaction at the electrode followed by product formation.Products of electrode reactions are assumed to be removed. An exampleis hydrogen evolution reaction by photo-electrocatalytic (PEC) water split-ting. [32, 33] For PEC water splitting, cations, anions, and undissociatedmolecules correspond to H O + , OH − , and water molecules, respectively. Wemainly consider the case that one of the ion species is discharged at theelectrode and the other species are inert and reflected at the electrode asdescribed in Fig. 1 (a). In Sec. 6, we consider the special case that bothcation and anion species are discharged at the electrode as described in Fig.1 (b). The situation corresponds to overall water splitting into H and O using a single photocatalyst such as GaN-ZnO and ZnGeN -ZnO under lightillumination. [34, 35] For simplicity, we assume that reactive sites on theelectrode are structureless. In experiments, the electrode can be regarded ashomogeneous when the electrode is modified with molecular or metal complexcocatalyst to realize both oxidation and reduction reactions on the electrode.Some electrodes with nano-particulate photocatalysts could also be regardedas homogeneous.Ion transport was described by the Nernst–Planck equations and the elec-tric field satisfied the Poisson equation. By taking into account the associa-tion and dissociation reactions of ions with diffusion and ion migration underthe electric field E ( x ), the governing equations in the bulk phase are given4y [5, 14–21, 23–25] ∂∂t n + = ∂∂x D + (cid:18) ∂n + ∂x − eEn + k B T (cid:19) + k d m − k a n + n − , (1) ∂∂t n − = ∂∂x D − (cid:18) ∂n − ∂x + eEn − k B T (cid:19) + k d m − k a n + n − , (2) ∂∂t m = D m ∂ ∂x m − k d m + k a n + n − , (3) ∂∂x E ( x ) = 4 πeǫ ( n + − n − ) , (4)where k d and k a denote the dissociation rate constant and the association rateconstant, respectively. Equation (4) is the Poisson equation. − D + ( ∂n + ) / ( ∂x )and D + eEn + / ( k B T ) on the right-hand side of Eq. (1) represent the diffu-sive flux and the ion electro-migration by the electric field E , respectively. D + / ( k B T ) can be regarded as the mobility using the Einstein relation. Whenthe electrostatic potential is controllable by potentiostat operation, Eqs. (1)-(4) may be more conveniently expressed by using the electrostatic potentialrather than using the electric field. However, we study the ion currentsand the electrostatic potential induced by the ion discharge reactions at theelectrode without imposing the external electric field. In this case, the elec-trostatic potential is not a controllable parameter. To avoid a redundantintegration constant associated with the electrostatic potential we expressEqs. (1)-(4) using the electric field.In the following, we consider steady states. The left-hand sides of Eqs.(1) and (2) are zero in steady states. The electrostatic potential differencerelative to the electrostatic potential at L can be defined by V ( x ) = − Z xL dx E ( x ) , (5)and can be calculated by solving Eqs. (1)–(4) under proper boundary condi-tions. The total potential difference between the electrode surface and bulk(overpotential) is V (0). At a sufficiently large distance away from the electrode surface, the ionconcentrations are not affected by surface reactions. Correspondingly, we set5oundary conditions at x = L and assumed equal concentrations of cationsand anions given by n b . Boundary conditions at x = L are given by n + | x = L = n b , (6) n − | x = L = n b , (7) m | x = L = m b , (8)where k d m b = k a n b is satisfied in the bulk. We consider the situation thatthe external electric field is absent in the bulk E ( L ) = 0 . (9)In ideal situation of photo-electrochemical (PEC) water splitting, the exter-nal electric field is not applied. We also assumed a metal electrode, wherethe net image charge in the electrode is the same as the net charge in theelectrolyte system. Therefore, the electric fields at the boundaries enclosingboth the real charges in the electrolyte and the image charges are zero. Thissituation can also be viewed from another perspective. Because we considerneutral products from a neutral electrolyte by electrode reactions, the netcharge in the system including those induced in the electrode should be neu-tral. In this case, electric fields at some distance sufficiently far from theelectrodes should be zero.[6]We consider the case where cations are reduced at the interface betweenthe electrolyte and electrode and the anions are reflected as shown in Fig. 1(a). Boundary conditions at the cathode located at x = 0 can be expressedas n + (0) = n s , (10) D − (cid:18) ∂n − ∂x + eEn − k B T (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) x =0 = 0 , (11) ∂m∂x (cid:12)(cid:12)(cid:12)(cid:12) x =0 = 0 . (12)Equation (12) indicates that undissociated molecules are reflected by thecathode. The left-hand side of (11) represents the concentration current ofanions. Equation (11) indicates that anions are reflected by the cathode. InEq. (10), n s represents the net effect of surface reactions on n + (0). As aplausible boundary condition for irreversible surface reactions, we considered6hat n s = 0 on the right-hand side of Eq. (10). The boundary conditionis suitable to study the case when the overall reactions are limited by iontransport to the electrode surface. In general, we set n s being constant onthe right-hand side of Eq. (10). We can express concentrations of cations and anions in dimensionlessunits as N + = n + n b , N − = n − n b , N s = n s n b and N m = mm b . (13)For convenience, we introduce [12]¯ D = 2 D + D − D + + D − , ∆ = D + − D − D + + D − . (14)Using the above definitions, the original diffusion constants can be obtainedby D + = ¯ D/ (1 − ∆), and D − = ¯ D/ (1 + ∆). By using the Debye lengthgiven by [12, 13, 36–38] λ D = (cid:18) ǫk B T πe n b (cid:19) / , (15)we can express the characteristic electric field strength and diffusion time as[12] E = k B Teλ D = (cid:18) πn b k B Tǫ (cid:19) / , t = λ D ¯ D = ǫk B T πe n b ¯ D . (16)We also introduce dimensionless variables given by y = xλ D , τ = tt , ξ = EE . (17)Finally, we express the reaction rate constants in the dimensionless form,¯ k a = λ D k a n b ¯ D , ¯ k d = λ D k d ¯ D . (18)7ere, ¯ k a and ¯ k d represent the dimensionless ion association rate constantand ion dissociation rate constant, respectively. Equation (1)–(4) can berewritten using the dimensionless variables as ∂ N + ∂y − ∂∂y ( ξN + ) = (1 − ∆) (cid:20) ∂N + ∂τ − ¯ k a ( N m − N + N − ) (cid:21) , (19) ∂ N − ∂y + ∂∂y ( ξN − ) = (1 + ∆) (cid:20) ∂N − ∂τ − ¯ k a ( N m − N + N − ) (cid:21) , (20) ∂ N m ∂y = ¯ DD m (cid:20) ∂N − ∂τ + ¯ k d ( N m − N + N − ) (cid:21) , (21) ∂∂y ξ = 12 ( N + − N − ) . (22)Similarly, boundary conditions at y = L/λ D are given by N + | y = L/λ D = 1 , (23) N − | y = L/λ D = 1 , (24) N m | y = L/λ D = 1 , (25)and boundary conditions at y = 0 [Eqs. (10)-(12)] are given by N + (0) = N s , (26) (cid:18) ∂N − ∂y + ξN − (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) y =0 = 0 , (27) ∂N m ∂y (cid:12)(cid:12)(cid:12)(cid:12) y =0 = 0 . (28)Finally, the boundary condition for the electric field given by Eq. (9) can beexpressed as ξ | y = L/λ D = 0 . (29)The potential difference given by Eq. (5) can be expressed as V ( y ) = − E λ D Z yy = L/λ D dyξ ( y ) . (30)8elow we consider the case where the condition¯ k a n b m b ¯ DD m ≪ k d m b = k a n b , it can be shownthat the right-hand side of Eq. (21) can be ignored in steady states. Forthis case, the dimensionless concentration of undissociated molecules can beapproximated by N m = 1 from the solution of diffusion equation under theboundary conditions of Eqs. (25) and (28). This situation represents thatthe ratio of dissociated to undissociated neutral compound given by n b /m b is so small in weak electrolytes that the density of the neutral compoundcan be approximated as constant. Obviously, Eq. (31) is satisfied for water.For general weak electrolytes, additional condition for the concentration ofthe undissociated neutral compound may be needed. By using the degreeof ionization α , we have n b /m b = α/ (1 − α ), where the numerator and thedenominator indicate the increase of the amount of cation and the decrease inthe amount of undissociated molecules for the extent of ionization reaction,respectively. When we consider the dissociation of an initial amount n of theundissociated compound, the equilibrium dissociation constant is expressedas K a = n α / (1 − α ), where n is divided by the standard concentration 1mol/L. In weak electrolytes, n typically satisfies the relation n ≫ K a / α ≈ p K a /n . By substituting n b /m b ≈ α ≈ p K a /n , Eq. (31)can be rewritten as K a /n ≪ ( D m / ¯ D ) / ¯ k a . Using D m / ¯ D ∼
10 as a typicalvalues of the diffusion constants of Carboxylic acids [39] and 1 < ¯ k a <
10, wefind K a /n ≪ − ∼ − . For most weak electrolytes satisfying K a ≤ − ,the condition can be fulfilled under normal values of n . Below, we assumethat the condition is satisfied and put N m = 1 in Eqs. (19)-(20).
3. Analytical Results
When ion discharge reactions proceed at the electrode-electrolyte inter-face in a steady state, the time-derivatives in Eqs. (19)-(20) can be set zero.Here, we derive time independent approximate solutions of Eqs. (19)-(29) byassuming N m = 1 as explained in the previous section. As explained later,the steady–state analytical solution is obtained under a certain conditionrelevant to weak electrolytes. A separate consideration is needed for strongelectrolytes as shown in Appendix A.9ecause we are interested in the deviations from equilibrium concentra-tions induced by surface reactions, we introduce sum and change as [12] S = N + + N − , C = N + − N − , where S and C characterize the total massdensity and charge density, respectively. To investigate the deviations fromequilibrium states, we introduce δS = S − , δC = C. (32)In steady states, Eqs. (19), (20) and (22) can be rewritten using thesevariables as ∂ δS∂y − ∂∂y ( ξδC ) = 2¯ k a (cid:20) δS + 14 (cid:0) δS − δC (cid:1)(cid:21) , (33) ∂ δC∂y − ∂∂y [ ξ ( δS + 2)] = − k a (cid:20) δS + 14 (cid:0) δS − δC (cid:1)(cid:21) , (34) ∂∂y ξ = 12 δC, (35)where N m = 1 is set by assuming the condition Eq. (31) is satisfied. Ina weak electrolyte, the dissociation of a small amount of the solute can berepresented by the right-hand sides of Eqs. (33) and (34).The boundary conditions at y = L/λ D can be written as δS | y = L/λ D = 0 , (36) δC | y = L/λ D = 0 . (37)The boundary conditions at y = 0 can be expressed as δS + δC (cid:12)(cid:12)(cid:12)(cid:12) y =0 = N s − , (38) ∂∂y ( δS − δC ) + ξ (2 + δS − δC ) (cid:12)(cid:12)(cid:12)(cid:12) y =0 = 0 . (39)The boundary condition for the electric field is given by Eq. (29).By linearizing the above equations, we solved the following set of equa-tions, ∂ δS∂y = 2¯ k a δS, (40) ∂ δC∂y − ∂∂y ξ = − k a δS. (41)10y substituting Eq. (35) into Eq. (41), we obtain ∂ δC∂y − δC = − k a δS. (42)If the right-hand side is zero, Eq. (42) is equal to the equation derived byDebye and Falkenhagen in the steady state. [40] The right-hand side of Eq.(42) reflects ionization of molecules. When the ion concentrations are con-served, the Debye–Falkenhagen equation is obtained, and the length scale ofthe charge density variation is characterized by the Debye length. Equation(42) indicates that the charge density variation is still characterized by theDebye length even when the ion concentrations are non-conserved by ioniza-tion of molecules: the ionization effect is represented by the inhomogeneousterm in Eq. (42) when the linearization approximation is used. The validityof the linearization will be discussed later. Hereafter, we solve Eqs. (40),(42) and (35) using the boundary conditions given by Eqs. (36)-(39) and(29).The solution of Eq. (40) can be expressed as δS ( y ) = C s exp (cid:16) − p k a y (cid:17) (43)using the boundary condition given by Eq. (36) in the limit of L → ∞ ,where C s is an integration constant. Using the boundary condition given byEq. (37) in the limit of L → ∞ , we obtained from Eq. (42) δC ( y ) = − k a ∆ / (2¯ k a − C s exp (cid:16) − p k a y (cid:17) + C c exp ( − y ) , (44)where C c is an integration constant. By substituting them into Eqs. (38) and(39) and linearizing the equations, we obtain C s = 0 and C c = − − N s )up to the linear order in 1 − N s . By substituting δC and integrating Eq.(35), we found that ξ ( y ) = (1 − N s ) exp ( − y ) + C ξ , where C ξ is an integrationconstant. The boundary condition given by Eq. (29) can be expressed as ξ ( y ) → L → ∞ . Therefore, we find that C ξ = 0. Theseresults can be written as δC ( y ) = − − N s ) exp( − y ) , (45) ξ ( y ) = (1 − N s ) exp( − y ) . (46)11he electric field strength can be written as E ( y ) = k B Teλ D (1 − N s ) exp( − y ) . (47)From Eq. (30) and taking the limit of L → ∞ , the potential profile can beobtained as V ( y ) = k B Te (1 − N s ) exp ( − y ) . (48)Using N + = (2 + δC ) / N − = (2 − δC ) /
2, we obtained the ionic concen-tration profiles as N + ( y ) = 1 − (1 − N s ) exp ( − y ) , (49) N − ( y ) = 1 + (1 − N s ) exp ( − y ) . (50)Exponential screening of Eqs. (49)-(50) reminds us of the exponential de-crease of electric field formed near charged surfaces in strong electrolytesunder local equilibrium [13, 36–38, 41]; the solution for a certain bound-ary condition is known as the Debye–H¨uckel theory. However, the situa-tion considered here is different from that of the Debye–H¨uckel theory. TheDebye–H¨uckel theory is based on the Poisson–Boltzmann equation, whereion concentrations are assumed to obey local equilibrium. Here, we solvedthe e-PNP equations where the effect of ion currents is taken into accountand deviations from the local equilibrium distributions can be studied.The dimensionless concentration current density in Eq. (26) was obtainedfrom Eqs. (49) and (46) as (cid:18) ∂N + ∂y − ξN + (cid:19) = (1 − N s ) exp( − y ) . (51)The ion current density can be expressed as e D + λ D n b (cid:18) − n s n b (cid:19) exp( − x/λ D ) , (52)and its value at the electrode surface is given by e D + λ D n b (cid:18) − n s n b (cid:19) , (53)12sing a linearized approximation.The above results were derived by linearizing the Poisson–Nernst–Planckequations with the association and dissociation reactions. The linearizationis valid if | ∂ ( ξδC ) /∂y | on the left-hand side of Eq. (33) is smaller than2¯ k a | δS | on the right-hand side. The linearized solutions indicate δC∂ξ/∂y = ξ∂δC/∂y and | ∂ ( ξδC ) /∂y | = 2 δC∂ξ/∂y . The magnitude of charge density δC would decrease by increasing the distance from the electrode surface and δC∂ξ/∂y = δC / ≈ − N s ) near the electrode surface. If the electrodeis perfectly reactive, we have δC∂ξ/∂y ≈ | ∂ ( ξδC ) /∂y | = 2 δC∂ξ/∂y ≈ δS ≈ N − ≈ | δC | , we also obtain δS ≈ | ∂ ( ξδC ) /∂y | ≪ k a | δS | can be expressedas ¯ k a >
1. Otherwise, the coupling between ξ and δC in Eq. (33) cannot beignored. The linearization condition will be numerically studied in Sec. 5.
4. Interpretation of the linearization condition
As mentioned above, the linearization condition of the e-PNP equationsis given by ¯ k a >
1. The linearization condition can be interpreted in terms oftwo time scales using the definition ¯ k a = λ D k a n b / ¯ D . The value of 1 / ( k a n b )gives the time scale of recovering charge neutrality by bulk reactions. Thetime scale of diffusion over the Debye length is given by t = λ D / ¯ D ; theDebye length characterizes the spatial extent of the deviation from chargeneutrality. ¯ k a > k a > k a <
1, charge neutrality is mainly recovered by iondiffusion.As shown below, the linearization condition is satisfied when ion associ-ation is diffusion limited. The diffusion-limited ion association rate can berelevant for weak electrolytes, where the intrinsic association rate constantsare large, and the dissociation rate constants are small. For example, theassociation rate of water can be estimated as k a = 1 . × L/(mol s) us-ing the diffusion limited bulk rate constant 4 π ( D + + D − ) R eff together withthe effective reaction radius given by R eff = | r c | / [1 − exp( −| r c | /R )], [42–45] where | r c | = e / ( ǫk B T ) is the Onsager radius and the diffusion con-stants of H O + and OH − are 9 . × − m /s and 5 . × − m /s, respec-tively. [46] In the above estimation, we also used ǫ = 80 and R = 0 . k a value is close to the literature value givenby k a = 1 . × L/(mol s). [43, 44, 47] This indicates that the effect ofdiffusive interaction can be ignored owing to the low degree of ionization inweak electrolytes. [48] The similar estimation is also possible for other weakelectrolytes. [49] By noticing that the Debye length can be expressed usingthe Onsager length as λ D = (8 π | r c | n b ) − / , (54)¯ k a = λ D k a n b / ¯ D can be expressed as¯ k a = 14 ( D + + D − ) D + D − − exp( −| r c | /R ) . (55)Because the minimum value of ( D + + D − ) / ( D + D − ) is obtained when D + = D − , we find ¯ k a > − exp( −| r c | /R ) > ǫ and the reaction radius. For water, we find ¯ k a = 2 . k a > D + + D − ) D + D − − exp( −| r c | /R ) > . (57)using Eq. (55) and is always satisfied. The diffusion–limited association iscommonly adopted for weak electrolytes because the ion association rapidlyproceeds due to the Coulombic interaction without screening; [50, 51] thescreening by an ionic atmosphere only occurs in strong electrolytes.In calculating the diffusion–limited association rate, the concentrationof counterions at the distance orders of magnitude larger than the Onsagerradius is assumed to be homogeneous. [42, 43, 50] On the other hand, thedensity profiles obtained by solving the 1–dimensional e-PNP equations arecharacterized by the Debye length. If Onsager radius is orders of magnitudesmaller than the Debye length, the combined approach could be utilized. Theseparation of length scales is satisfied when the relative dielectric constant issufficiently high. For water, the Onsager length is 0 . µ m.14egarding the dissociation rate constant k d , we did not take into accountthe effect of electric fields on the ion dissociation rate. Strictly speaking, thedissociation rates of geminate ions have stronger electric field dependencecompared to the association rates in weak electrolytes. [50, 52–54] However,the field dependence of the dissociation rates is negligibly small in the absenceof external fields because the largest internal electrostatic potential differencein the spatial extent of the Debye length is on the same scale as thermalenergy at room temperature. The linearization condition will be numericallyjustified in the next section.
5. Numerical justification of the linearization condition
In this section, we confirm the linearization condition numerically. In asteady state, we numerically solved Eqs. (19)–(29) by the relaxation method.We present the numerical results obtained using the boundary conditionincluding the effect of the electrode.As a concrete example, we studied the photo-induced electrochemical re-duction of water at electrode surfaces yielding hydrogen. Water is a weakelectrolyte, and the dissociation of a small amount of water can be repre-sented by the association and dissociation reaction terms in Eqs. (33) and(34). The diffusion constants of H O + and OH − are 9 . × − m /s and5 . × − m /s, respectively. [46] By substituting D + = 9 . × − m /sand D − = 5 . × − m /s into Eq. (14), we obtained ¯ D = 6 . × − m /sand ∆ = 0 .
28. Using n b = 1 × − mol/L, the Debye length was obtainedas λ D = 9 . × − ∼ − m. The linearization condition of the e-PNPequations is given by ¯ k a >
1, where the rate of ion association given by k a n b is much higher than the inverse of the time required to diffuse over the De-bye length given by ¯ D/λ D . In this section, we study the case when ¯ k a = 2 . k a = 0 .
02; the linearized results shown in Sec. 3 could be applicableto the case of ¯ k a = 2 . k a = 0 .
02. For brevity, thecase of ¯ k a = 2 . k a = 0 .
02 isclassified as strong electrolytes. In weak electrolytes, the intrinsic associationrate constants are large, and the dissociation rate constants are small. Thelinearization condition is obtained when the association rate is diffusion lim-ited because the intrinsic association rate constant is large. The linearizationcondition is most likely satisfied for weak electrolytes. On the contrary, theassociation rate constants of strong electrolytes are small and can be ignored.15 -3 -2 E l e c t r o s t a t i c po t en t i a l ( V ) µ m) Weak electrolytesStrong electrolytes
Figure 2: (Color online) Electrostatic potential as a function of distance from the electrodewhen N s = 0 .
5. Thick lines indicate ¯ k a = 2 . k a = 0 .
02 (strong electrolytes). The solid lines represent the numerical solutions. (Red)dash-dotted lines represent the results of the Nernst equation given by Eq. (58). (Red)long dashed lines represent the results of Eq. (48). (Blue) short dashed lines represent theresults of the Henderson–Planck equation given by Eq. (62).
16n Fig. 2, we compare the results of the Nernst equation, the Henderson–Planck equation, and Eq. (48). The electrostatic potential profile of theNernst equation can be expressed as V N ( x ) = k B Te ln n b n + ( x ) , (58)which is equivalent to the local equilibrium assumption given by n + ( x ) n b = exp (cid:20) − eV N ( x ) k B T (cid:21) , (59)where ion mass flows are assumed to be absent. Intuitively, the Nernstequation could be thought of as applicable to weak electrolytes when the iondischarge reaction at the electrode is slow so that the resultant deviation fromcharge neutrality near the surface is recovered by bulk reactions rather thanion diffusion. Indeed, we can show that the electrostatic potential profilegiven by Eq. (48) is obtained by linearizing Eq. (58) by assuming [ n b − n + ( x )] /n b ≪ V N ( x ) = − k B Te ln (cid:18) n + ( x ) − n b n b (cid:19) (60) ≈ k B Te (cid:18) − n + ( x ) n b (cid:19) . (61)However, it should be reminded that Eq. (61) is obtained from the e-PNPequations when the linearization condition ¯ k a > n b − n + ( x )] /n b . In other words, if ¯ k a > n b − n + ( x )] /n b ∼ k a <
1. In the absence of associationand dissociation reactions, approximate analytical solutions of the Poisson–Nernst–Planck equations have been studied more than decades. As shown in17ppendix A and B, the potential difference known as the diffusion potentialis obtained, [3, 10–12] V HP ( x ) = − Z x =0 x = L dxE ( x ) = k B Te D + − D − D + + D − ln (cid:18) n b n + ( x ) (cid:19) . (62)Equation (62) is also called the Henderson–Planck equation. In strong elec-trolytes, the potential difference is reduced from the Nernst equation by thefactor given by ( D + − D − ) / ( D + + D − ). The Henderson–Planck equationrepresents the diffusion potential caused by different diffusion constants ofanions and cations. [3, 10–12]Coming back to Fig. 2, we notice that the result of the Nernst equationis closer to the numerical result than that of the Henderson–Planck equationwhen ¯ k a = 2 . k a >
1. The Nernst equation wasobtained by assuming local equilibrium, where mass currents were not takeninto account. The deviation found near the electrode could be caused bymass currents. When charge neutrality is broken near the electrode surface,both ionization and ion transport occur to recover charge neutrality. Becauseof ion transport, the local equilibrium assumption used to obtain the Nernstequation is violated, particularly near the electrode surface. The Henderson–Planck equation is derived for strong electrolytes, and its results cannot beapplied to weak electrolytes.When ¯ k a = 0 .
02, the linearization condition ¯ k a > k a = 0 .
02, the result of theHenderson–Planck equation is closer to the numerical result than that of theNernst equation. The Henderson–Planck equation is derived for strong elec-trolytes, where ionization is not taken into account. In contrast, the Nernstequation is derived by assuming local equilibrium, which can be violated bymass currents. We should note that breakdown of charge neutrality can berecovered only by mass currents in strong electrolytes.For ¯ k a = 0 .
2, the numerical results deviate from both the results of Eq.(48) and the Henderson-Planck equation. (Data not shown) The results areconsistent with the linearization condition given by ¯ k a > on c en t r a t i on Distance ( µ m) Strong electrolytesWeak electrolytes
Figure 3: (Color online) Concentrations of cations and anions as a function of distancefrom the electrode when N s = 0 .
5. Thick black lines indicate ¯ k a = 2 . k a = 0 .
02 (strong electrolytes). The solid lines representcation concentrations and dashed lines represent anion concentrations.
19n Fig. 3, we show concentrations of cations and anions as a function ofdistance from the electrode surface in strong and weak electrolytes. In weakelectrolytes, both cation and anion concentrations deviate from the bulkconcentration when the distance from the electrode is within a few Debyelengths. The result is consistent with the previous result obtained by as-suming weak external electric field. [25] Concentration profiles of cations arelong-ranged in strong electrolytes. Interestingly, the cation and anion con-centrations were different when the distance from the electrode was withina few Debye lengths in both strong and weak electrolytes. In other words,charge density caused by breakdown of charge neutrality was observed ap-proximately within the surface layer to a depth of a few Debye lengths. Theanion concentration profile is similar to that of cations outside this layer instrong electrolytes due to charge neutrality. The difference between the spa-tial range of charge density and that of mass density for strong electrolyteshas been pointed out. [4] Although small, the increase of anion concentra-tion near the electrode surface seen in Fig. 3 could be caused by ionizationof a few unionized molecules present in the strong electrolyte characterizedby ¯ k a = 0 .
02. We conclude that the linearization condition of the e-PNPequations is given by ¯ k a >
6. Both cations and anions are discharged at the electrode surface
As indicated in Eq. (43), the total charge density gradient could be scaledby the factor 1 / p k a from the Debye length when the total charge densitydrops by allowing anions to be discharged as well as cations at the electrodesurface as shown in Fig. 1 (b). In this section, we consider such situation.The boundary conditions at x = 0 ( y = 0) given by Eq. (11) changes into n − (0) = n s − . (63)As in Eq. (10) n s − represents the net effect of surface reactions on n − (0).The boundary condition at y = 0 obtained from Eq. (63) can be writtenas δS + δC + 2(1 − N s − ) | y =0 = 0 , (64)where N s − = n s − /n b as N s in Eq. (13). Similarly, we have δS + δC + 2(1 − N s ) | y =0 =0 from Eq. (38). By substituting Eqs. (43)-(44) into the above boundary20 Figure 4: Concentrations of cations and anions as a function of distance from the electrode.The solid lines indicate cations and the dashed lines indicate anions. The red lines representthe case when the perfectly reactive boundary condition ( k s → ∞ and n s = 0) is assumedfor cations with inert anions. The black lines represent the case when anions are alsoperfectly reactive at the electrode surface. The thick lines represent the exact numericalresults. The thin lines represent the approximate analytical results. The thin solid redline, dashed red line, solid black line and dashed black line is obtained from Eq. (49), Eq.(50), Eq. (65), and Eq. (66), respectively. C s = − N s + N s − and C c = N s − N s − − (2 − N s − N s − )2¯ k a ∆ / (2¯ k a − N s = N s − = 0, we obtain, N + = 1 − (cid:18) − k a ∆2¯ k a − (cid:19) exp (cid:16) − p k a y (cid:17) − k a ∆2¯ k a − − y ) , (65) N − = 1 − (cid:18) k a ∆2¯ k a − (cid:19) exp (cid:16) − p k a y (cid:17) + 2¯ k a ∆2¯ k a − − y ) , (66) ξ = − p k a ∆2¯ k a − (cid:16) − p k a y (cid:17) + 2¯ k a ∆2¯ k a − − y ) , (67) V = k B Te ∆2¯ k a − h k a exp ( − y ) − exp (cid:16) − p k a y (cid:17)i . (68)The factor p k a indicates the length scale given by the inverse of p k a n b / ¯ D .For diffusion limited ion association in weak electrolytes, ¯ k a > p ¯ D/ (2 k a n b ) might be smaller than the Debye lengthbecause ¯ k a > λ D > p ¯ D/ ( k a n b ). The concentration gra-dients are essentially characterized by the length scale given by p ¯ D/ (2 k a n b )rather than the Debye length when both cations and anions are discharged atthe electrode. Interestingly, the electrostatic potential difference is essentiallycharacterized by the Debye length because 2¯ k a > p ¯ D/ (2 k a n b ) characterizes theion concentration gradients. In an entirely different context but using thesimilar method, it was found that density profiles could be characterized bythe association reaction-diffusion length scale in certain multi-ionic solutions.[49] p ¯ D/ (2 k a n b ) represents the length scale of the diffusive migration withinthe life-time of ions which will disappear as a consequence of the ion associ-ation reaction. The length scale is called the Kuramoto length. [26–30] TheKuramoto length characterizes the spatial correlation of density fluctuationsaround a uniform concentration state. The density fluctuations beyond thelength scale given by the Kuramoto length are independent. Here, we showthat the ion density drop at the electrode surface can be localized withinthe Kuramoto length. The Kuramoto length can be expressed using the22iffusion-controlled rate shown in Sec. 4 as s D + D − [1 − exp( −| r c | /R )]4 π ( D + + D − ) | r c | n b . (69)In water, the Kuramoto length is given by λ D / ≈ . µ m.We compare the above results with the numerical solutions as in theprevious section. As an example, we consider photo-induced water splittingreactions by photocatalysts. We set ¯ k a = 2 .
7. Conclusion
The results of this paper are related to photo-electrochemical (PEC) wa-ter splitting. Although theoretically estimated PEC conversion efficiency23
Figure 5: (Color online) Electrostatic potential as a function of distance from the electrode.The thin (red) lines represent the case when the perfectly reactive boundary condition( n s = 0) is assumed for cations with inert anions. The thick (black) lines represent thecase when anions are also perfectly reactive at the electrode surface. The solid lines indicatethe exact numerical results. The dashed lines represent the approximate analytical results.The (red) dashed line represents V ( y ) = ( k B T /e )(1 − N + ) obtained from Eq. (48) usingthe numerical results for cation concentration. The (black) dashed line is obtained fromEq. (68). n b /λ D . By using D + ≈ − m /s, λ D ≈ − m, the ion current density can be expressed as eD + n b λ D × A/cm = 0 .
01 mA/cm , (70)when the boundary condition is given by n + (0) = 0 corresponding to the casewhere the back reaction can be ignored because of fast forward reactionsat the interface, where 96500 C/mol is approximated by 10 C/mol. Theabove result matches with that of linearization given by Eq. (53) and is30% smaller than the exact numerical result. We have also shown that theconcentration gradient is characterized by the Kuramoto length when bothcation reduction and anion oxidation reactions occur at the electrode for 1:1monovalent weak electrolytes. The Kuramoto length is about the half of theDebye length. As a result, the ion current can be twice larger than thatobtained when anions are not oxidized at the electrode surface. However, torealize 10% PEC conversion efficiency under 1 sun, a current density of 10 mAcm − is required. This result indicates that ion currents may limit the overallefficiency because of the small diffusion constant of H O + compared with thatof electrons. For photo-induced splitting of pure water, recent experimentalsetups have used internal convection caused by placing cathodic and anodicelectrodes in parallel configurations. [32, 33] In this configuration, an extraenergy source to generate convection was not needed. In principle, convectiveflows can be considered by introducing the hydrodynamic equation into theNernst–Planck equations. [63] An extension of the present work to includeconvective flows will be important to study the case when the current densityis high. 25he second limitation could be overpotential related to a charge densitygradient near the electrode. Our results indicate that the overpotential is atmost on the same scale as thermal energy at room temperature and can befurther reduced when both cation reduction and anion oxidation reactionsoccur at the electrode for 1:1 monovalent weak electrolytes. The overpotentialis caused by a charge density gradient whose spatial extent is characterizedby the Debye length. It should be noted that the charge density gradientcan be estimated from the pH gradient only for pure water. In general,charge density cannot be estimated from pH alone in the presence of otherion species.For simplicity, we considered the case that molecules could be ionizedto monovalent ions in the absence of buffer electrolytes composed of weakacid and/or weak base. We also did not consider supporting electrolytesto increase the conductivity of the bulk solution. For photo-induced water-splitting reactions, multi-ionic solutions are commonly used to control elec-trostatic potentials and engineer the pH of electrolytes. [20, 21, 64] There,large-scale pH gradients were observed. The observed large-scale pH gradi-ents could be related to the length-scale separation between charge densityprofiles and concentration profiles, which could be possible for multi-ionicsolutions. The effect of multiple ionic species on ion currents could be takeninto account using new numerical methods developed to calculate junctionpotentials,[11, 65] where the charge neutrality condition is unnecessary. Theanalytical method applied in this paper should be extended to multi-ionicsolutions in the future. Acknowledgement
This work was supported by the “Research Project for Future Develop-ment: Artificial Photosynthetic Chemical Process (ARPChem)” (METI/NEDO,Japan: 2012-2022).
Appendix A. Derivation of Eq. (62)
We performed a systematic expansion in 1 /τ when ¯ k a = 0 by applyingthe approach developed to calculate liquid junction potential to our problemof surface reactions. [12] When concentration change is restored only bydiffusion, a relevant variable is z = y/ √ τ . In terms of the new variable, Eq.2619)–(22) can be rewritten using δS = N + + N − − δC = N + − N − , and z by applying the chain rule for partials as1 τ ∂ δS∂z − ∂∂z (cid:18) ξ √ τ δC (cid:19) = ∂∂τ ( δS − ∆ δC ) − z τ ∂∂z ( δS − ∆ δC ) , (A.1)1 τ ∂ δC∂z − ∂∂z (cid:20) ξ √ τ ( δS + 2) (cid:21) = ∂∂τ ( δC − ∆ δS ) − z τ ∂∂z ( δC − ∆ δS ) , (A.2) ∂∂z ξ √ τ = 12 δC. (A.3)The diffusion may give rise to algebraic τ dependence, δS ( z, τ ) = δS ( z ) (0) + δS ( z ) (1) /τ + δS ( z ) (2) /τ + · · · , (A.4) δC ( z, τ ) = δC ( z ) (0) + δC ( z ) (1) /τ + δC ( z ) (2) /τ + · · · , (A.5) ξ ( z, τ ) √ τ = η (1) ( z ) τ + η (2) ( z ) τ + · · · . (A.6) η (0) ( z ) should be zero; otherwise, ξ ( z ) ∼ √ τ η (0) ( z ) diverges as τ → ∞ .Now, we study the lowest-order solution given by the order of 1 /τ . FromEq. (A.3), the charge neutrality condition is satisfied; i.e., δC ( z ) (0) = 0 . (A.7)By substituting Eqs. (A.4) and (A.5) into Eqs. (A.1)-(A.3), we obtained S (0) ( z ) and η (1) ( z ) of the order of 1 /τ as ∂ δS (0) ∂z = − z ∂∂z (cid:0) δS (0) (cid:1) , (A.8) − ∂∂z (cid:2) η (1) ( z )( δS (0) + 2) (cid:3) = ∆ z ∂∂z δS (0) . (A.9)Boundary conditions at z = L/ ( λ D √ τ ) are given by δS (0) (cid:12)(cid:12) z = L/ ( λ D √ τ ) = 0 , (A.10) δC (0) (cid:12)(cid:12) z = L/ ( λ D √ τ ) = 0 . (A.11)27he boundary condition at y = 0 is expressed as δS (0) (cid:12)(cid:12) z =0 ≈ N s − . (A.12)The potential difference can be obtained from V ( z ) = E λ D Z ∞ z dz η (1) ( z ) , (A.13)and the total potential difference between the electrode surface and bulk is V = V (0).First, we integrated Eq. (A.8) and obtained ∂∂z (cid:0) δS (0) (cid:1) = A exp (cid:18) − z (cid:19) , (A.14)where A is an integration constant. By integrating the above equation, wehave δS (0) ( z ) = A + A Z z dz exp (cid:18) − z (cid:19) , (A.15)where A is an integration constant.We note that R ∞ dz exp (cid:16) − z (cid:17) = √ π . Using the boundary conditionsgiven by Eqs. (A.10) and (A.12), we obtained A = 2 ( N s − , (A.16) A = 2 (1 − N s ) / √ π. (A.17)By substituting these equations into Eq. (A.15), it became δS (0) ( z ) = 2 (1 − N s ) [ − z/ , (A.18)where we used erf( z/
2) = (1 / √ π ) R z dz exp ( − z / z →∞ erf( z/
2) =1 and we have lim z →∞ δS (0) ( z ) = 0, as expected.Finally, using Eq. (A.14), we obtained from Eq. (A.9) ∂∂z (cid:2) η (1) ( z )( δS (0) ( z ) + 2) (cid:3) = 2∆ 1 − N s √ π ∂∂z exp (cid:18) − z (cid:19) . (A.19)28y further integration, we find η (1) ( z ) = ∆ − N s √ π exp (cid:16) − z (cid:17) + A N s + (1 − N s ) erf( z/ , (A.20)where we used lim z →∞ δS (0) = 0 and A is an integration constant.We consider the case that the electric field at the infinite distance fromthe electrode is zero. The L → ∞ limit of z is given by z → ∞ because wehave z = y/ √ τ , y = x/λ D , and x = L . Noting that the dimensionless fieldis denoted by η (1) , we obtained A = 0 from Eqs. (A.20) and (9).By substituting Eq. (A.20) into Eq. (A.13), the voltage difference wascalculated as V ( z ) = E λ D ∆ Z ∞ z dz − N s √ π exp (cid:16) − z (cid:17) N s + (1 − N s ) erf( z/
2) (A.21)= E λ D ∆ ln " N s + 1 − N s √ π Z z ′ dz exp (cid:18) − z (cid:19) z ′ = ∞ z ′ = z (A.22)= − E λ D ∆ ln h N s + (1 − N s ) erf (cid:16) z (cid:17)i . (A.23)The total voltage difference between the electrode surface and bulk was ob-tained from V (0) as V HP = − k B Te D + − D − D + + D − ln ( N s ) . (A.24)Although the concentration of cations should be lower than that of anionsnear the electrode because of the reduction of H O + at the electrode surface,the charge neutrality condition expressed by Eq. (A.7) was obtained as thelowest-order solution from the expansion of 1 /τ .Conventionally, Eq. (A.24) has been derived using the difference in chem-ical potentials at the boundaries of liquid junctions. [12] In the derivationin Appendix A, we used different boundary conditions, and the result indi-cates that the voltage change caused by electrolytic reactions at an electrodesurface can be approximated by the Henderson–Planck equation in strongelectrolytes. 29 ppendix B. The electrostatic potential under charge neutralitycondition Interestingly, the diffusion potential can be derived by assuming chargeneutrality regardless of presence of association and dissociation reactions.When the charges do not accumulate at the interface between the electrodeand electrolyte, the positive and negative currents should be equal. Theequality of the positive and negative currents in the steady state can beexpressed as D + (cid:20) ∂n + ∂x − (cid:18) eEn + k B T (cid:19)(cid:21) = D − (cid:20) ∂n − ∂x + (cid:18) eEn − k B T (cid:19)(cid:21) , (B.1)and the charge neutrality condition is given by n + = n − . By introducing thecharge neutrality condition, we obtain from Eq. (B.1) [12] E = k B Te D + − D − D + + D − n + ∂n + ∂x , (B.2)and Eq. (62) known as the Henderson–Planck equation for 1:1 monovalentelectrolytes even in the presence of association and dissociation reactions.The fact that Eq. (48) [or equivalently Eq. (61)] is not obtained by assuminglocal charge neutrality indicates the subtle aspects of charge neutrality tocalculate the electrostatic potential. ReferencesReferences [1] V. G. Levich, Physicochemical Hydrodynamics, 2nd Edition, Prentice-Hall, EnglewoodCliffs,NJ, 1962.[2] J. Newman, K. E. Thomas-Alyea, Electrochemical Systems, 3rd Edition,John Wiley & Sons, Inc., Honoken NJ, 2004.[3] A. J. Bard, L. R. Faulkner, Electrochemical Methods: Fundamentalsand Applications, John Wiley & Sons, Inc, 2001.[4] S. Maf´e, J. Pellicer, V. M. Aguilella, Ionic Transport and Space ChargeDensity in Electrolytic Solutions as Described by Nernst-Planck andPoisson Equations, J. Phys. Chem. 90 (22) (1986) 6045–6050.305] M. Bier, O. Palusinski, R. Mosher, D. Saville, Electrophoresis: Mathe-matical Modeling and Computer Simulation, Science 219 (4590) (1983)1281–1287.[6] R. P. Buck, Kinetics of Bulk and Interfacial Ionic Motion: MicroscopicBases and Limits for the Nernst–Planck Equation Applied to MembraneSystems, Journal of Membr. Sci. 17 (1) (1984) 1 – 62.[7] H. Cohen, J. W. Cooley, The Numerical Solution of the Time-DependentNernst-Planck Equations, Biophys. J 5 (2) (1965) 145–162.[8] I. Rubinstein, L. Shtilman, Voltage against Current Curves of CationExchange Membranes, J. Chem. Soc., Faraday Trans. 2 75 (1979) 231–246.[9] Y. Kim, W. S. Walker, D. F. Lawler, The Painlev´e Equationof the Second Kind for the Binary Ionic Transport in Diffu-sion Boundary Layers near Ion-Exchange Membranes at Overlim-iting Current, J. Electroanal. Chem. 639 (1) (2010) 59 – 66. doi:https://doi.org/10.1016/j.jelechem.2009.11.019 .[10] P. Henderson, Zur Thermodynamik der Fl¨ussigkeitsketten, Z Phys Chem59 (1907) 118–127.[11] E. J. F. Dickinson, J. G. Limon-Petersen, R. G. Compton, The Elec-troneutrality Approximation in Electrochemistry, J. Solid State Elec-trochem. 15 (7) (2011) 1335–1345.[12] J. L. Jackson, Charge Neutrality in Electrolytic Solutions and the LiquidJunction Potential, J. Phys. Chem. 78 (1974) 2060–2064.[13] M. Z. Bazant, K. Thornton, A. Ajdari, Diffuse-Charge Dynamics inElectrochemical Systems, Phys. Rev. E 70 (2004) 021506(1–24).[14] G. Grossman, Water Dissociation Effects in Ion Transport through Com-posite Membrane, J. Phys. Chem. 80 (14) (1976) 1616–1625.[15] V. V. Nikonenko, N. D. Pismenskaya, E. I. Belova, P. Sistat, P. Huguet,G. Pourcelly, C. Larchet, Intensive Current Transfer in Membrane Sys-tems: Modelling, Mechanisms and Application in Electrodialysis , Adv.Colloid and Interface Sci. 160 (1-2) (2010) 101–123.3116] L.-J. Cheng, H.-C. Chang, Microscale pH Regulation by Splitting Water,Biomicrofluidics 5 (4) (2011) 046502 (1–8).[17] V. I. Zabolotskii, V. V. Bugakov, M. V. Sharafan, R. K. Chermit, Trans-fer of Electrolyte ions and Water Dissociation in Anion-Exchange Mem-branes under Intense Current Conditions, Russ. J. Electrochem. 48 (6)(2012) 650–659.[18] Y. Kharkats, A. Sokirko, Theory of the Effect of Migration CurrentExaltation taking into Account Dissociation-Recombination Reactions,J. Electroanal. Chem. Interfacial Electrochem. 303 (1) (1991) 27 – 44.[19] C. P. Nielsen, H. Bruus, Transport-Limited Water Splitting at Ion-Selective Interfaces during Concentration Polarization, Phys. Rev. E 89(2014) 042405 (1–14).[20] J. Jin, K. Walczak, M. R. Singh, C. Karp, N. S. Lewis, C. Xiang, AnExperimental and Modeling/Simulation-Based Evaluation of the Effi-ciency and Operational Performance Characteristics of an Integrated,Membrane-Free, Neutral pH Solar-Driven Water-Splitting System, En-ergy Environ. Sci. 7 (2014) 3371–3380.[21] C. Xiang, A. Z. Weber, S. Ardo, A. Berger, Y. Chen, R. Coridan, K. T.Fountaine, S. Haussener, S. Hu, R. Liu, N. S. Lewis, M. A. Modestino,M. M. Shaner, M. R. Singh, J. C. Stevens, K. Sun, K. Walczak, Mod-eling, Simulation, and Implementation of Solar-Driven Water-SplittingDevices, Angew. Chem. Int. Ed. 55 (42) (2016) 12974–12988.[22] R. Femmer, A. Mani, M. Wessling, Ion Transport through Elec-trolyte/Polyelectrolyte Multi-layers, Sci. Rep. 5 (2015) 11583(1–12).[23] H. V. Parys, G. Telias, V. Nedashkivskyi, B. Mollay, I. Vandendael,S. V. Damme, J. Deconinck, A. Hubin, On the Modeling of Electro-chemical Systems with Simultaneous Gas Evolution. Case study: TheZinc Deposition Mechanism, Electrochim. Acta 55 (20) (2010) 5709 –5718. doi:https://doi.org/10.1016/j.electacta.2010.05.006 .[24] R. Kod´ym, V. F´ıla, D. ˇSnita, K. Bouzek, Poisson–Nernst–Planck Modelof Multiple Ion Transport across an Ion-Selective Membrane under Con-ditions close to Chlor-Alkali Electrolysis, J. Appl. Electrochem. 46 (6)(2016) 679–694. doi:10.1007/s10800-016-0945-1 .3225] K. J. Aoki, C. Li, T. Nishiumi, J. Chen, Electrolysis of Pure Water in aThin Layer Cell, J. Electroanal. Chem. 695 (Supplement C) (2013) 24 –29. doi:https://doi.org/10.1016/j.jelechem.2013.02.020 .[26] Y. Kuramoto, Fluctuations around Steady States in Chemi-cal Kinetics, Prog. Theor. Phys. 49 (5) (1973) 1782–1783. doi:10.1143/PTP.49.1782 .[27] Y. Kuramoto, Effects of Diffusion on the Fluctuations in OpenChemical Systems, Prog. Theor. Phys. 52 (2) (1974) 711–713. doi:10.1143/PTP.52.711 .[28] A. Nitzan, J. Ross, A Comment on Fluctuations around Nonequi-librium Steady States, J. Stat. Phys. 10 (5) (1974) 379–390. doi:10.1007/BF01008800 .[29] G. Nicolis, M. M. Mansour, Onset of Spatial Correlations in Nonequilib-rium Systems: A Master-Equation Description, Phys. Rev. A 29 (1984)2845–2853. doi:10.1103/PhysRevA.29.2845 .[30] K. Kitahara, K. Seki, S. Suzuki, On the Spatial Correlations inNonequilibrium Systems, J. Phys. Soc. Jpn. 59 (7) (1990) 2309–2311. doi:10.1143/JPSJ.59.2309 .[31] J. Wakou, K. Kitahara, M. Litniewski, J. Gorecki, Nonequilib-rium Spatial Correlations in Chemical Systems: Beyond Ornstein-Zernike, Int. J. Mod. Phys. C 13 (09) (2002) 1253–1262. doi:10.1142/S0129183102004108 .[32] Q. Wang, T. Hisatomi, Q. Jia, H. Tokudome, M. Zhong, C. Wang,Z. Pan, Z. Pan, T. Takata, M. Nakabayashi, N. Shibata, Y. Li, I. D.Sharp, A. Kudo, T. Yamada, K. Domen, Scalable Water Splitting onParticulate Photocatalyst Sheets with a Solar-to-Hydrogen Energy Con-version Efficiency Exceeding 1%, Nat. Mater. 15 (6) (2016) 611–615.[33] Q. Wang, T. Hisatomi, Y. Suzuki, Z. Pan, J. Seo, M. Katayama,T. Minegishi, H. Nishiyama, T. Takata, K. Seki, A. Kudo, T. Yamada,K. Domen, Particulate Photocatalyst Sheets Based on Carbon Con-ductor Layer for Efficient Z-Scheme Pure-Water Splitting at AmbientPressure, J. Am. Chem. Soc. 139 (4) (2017) 1675–1683.3334] Y. Lee, K. Teramura, M. Hara, K. Domen, Modification of(Zn1+xGe)(N2Ox) Solid Solution as a Visible Light Driven Photocata-lyst for Overall Water Splitting, Chem. Mater. 19 (8) (2007) 2120–2127. doi:10.1021/cm062980d .[35] K. Maeda, K. Teramura, D. Lu, T. Takata, N. Saito, Y. Inoue,K. Domen, Photocatalyst Releasing Hydrogen from Water, Nature440 (7082) (2006) 295. doi:http://dx.doi.org/10.1038/440295a .[36] P. W. Debye, E. H¨uckel, Zur Theorie der Elektrolyte. I. Gefrierpunkt-serniedrigung und Verwandte Erscheinungen, Phys. Z. 24 (1923) 185–206.[37] Y. Levin, Electrostatic Correlations: from Plasma to Biology, Rep. Prog.Phys. 65 (11) (2002) 1577–1632.[38] D. Andelman, Introduction to Electrostatics in Soft and Biological Mat-ter, in: W. Poon, D. Andelman (Eds.), Soft Condensed Matter Physicsin Molecular and Cell Biology, Taylor & Francis, New York, 2006, pp.97–122.[39] W. J. Albery, A. R. Greenwood, R. F. Kibble, Diffusion Coeffi-cients of Carboxylic Acids, Trans. Faraday Soc. 63 (1967) 360–368. doi:10.1039/TF9676300360 .[40] P. Debye, H. Falkenhagen, Dispersion der Leitf¨ahigkeit und der Dielektrizit¨atskonstante starker Elektrolyte. I., Phys. Z. 29 (1928) 121 – 132.[41] G. Barbero, D. Olivero, Ions and Nematic Surface Energy: Beyond theExponential Approximation for the Electric Field of Ionic Origin, Phys.Rev. E 65 (2002) 031701.[42] P. Debye, Reaction Rates in Ionic Solutions, J. Electrochem. Soc. 82 (1)(1942) 265–272. doi:10.1149/1.3071413 .[43] M. Eigen, L. de Maeyer, Self-Dissociation and Protonic Charge Trans-port in Water and Ice, Proc. R. Soc. (London) Ser. A 247 (1251) (1958)505–533. doi:10.1098/rspa.1958.0208 .[44] F. H. Stillinger, Proton Transfer Reactions and Kinetics in Water, in:H. EYRING, D. HENDERSON (Eds.), Theoretical Chemistry, Aca-demic Press, 1978, pp. 177 – 234.3445] W. C. Natzle, C. B. Moore, Recombination of Hydrogen Ion (H+) andHydroxide in Pure Liquid Water, J. Phys. Chem. 89 (12) (1985) 2605–2612. doi:10.1021/j100258a035 .[46] R. J. Silbey, R. A. Alberty, M. G. Bawendi, Physical Chemistry, 4thEdition, John Wiley & Sons, New Jersey, 2005.[47] D. A. McQuarrie, J. D. Simon, Physical Chemistry: a Molecular Ap-proach, University Science Books, California, 1997.[48] S. D. Traytak, M. Tachiya, Competition Effect in Diffusion-ControlledReactions between Ions, J. Phys.: Condens. Matter 19 (6) (2007) 065109.[49] A. T. Brown, W. C. K. Poon, C. Holm, J. de Graaf, Ionicscreening and dissociation are crucial for understanding chemicalself-propulsion in polar solvents, Soft Matter 13 (2017) 1200–1222. doi:10.1039/C6SM01867J .[50] L. Onsager, Deviations from Ohms Law in Weak Electrolytes, J. Chem.Phys. 2 (9) (1934) 599–615.[51] L. Onsager, Initial Recombination of Ions, Phys. Rev. 54 (1938) 554–557.[52] S. D. Traytak, The Rate Constant of Diffusion-Limited Reactions in aWeak External Electric Field, Chem. Phys. Lett. 181 (1991) 558–562. doi:10.1016/0009-2614(91)80313-M .[53] K. Isoda, N. Kouchi, Y. Hatano, M. Tachiya, The Effect of an ExternalElectric Field on Diffusion-Controlled Bulk ion Recombination, J. Chem.Phys. 100 (1994) 5874–5881. doi:10.1063/1.467099 .[54] M. Hilczer, M. Tachiya, Unified Theory of Geminate and Bulk Electron–Hole Recombination in Organic Solar Cells, J. Phys. Chem. C 114 (14)(2010) 6808–6813. doi:10.1021/jp912262h .[55] J. R. Bolton, S. J. Strickler, J. S. Connolly, Limiting and RealizableEfficiencies of Solar Photolysis of Water, Nature 316 (1985) 495 – 500.[56] M. C. Hanna, A. J. Nozik, Solar Conversion Efficiency of Photovoltaicand Photoelectrolysis Cells with Carrier Multiplication Absorbers, J.Appl. Phys. 100 (7) (2006) 74510(1–8).3557] M. F. Weber, M. J. Dignam, Efficiency of Splitting Water with Semicon-ducting Photoelectrodes, J. Electrochem. Soc. 131 (6) (1984) 1258–1265. doi:10.1149/1.2115797 .[58] M. Weber, M. Dignam, Splitting Water with SemiconductingPhotoelectrodes–Efficiency Considerations, Int. J. Hydrogen Energ.11 (4) (1986) 225 – 232.[59] L. C. Seitz, Z. Chen, A. J. Forman, B. A. Pinaud, J. D. Benck, T. F.Jaramillo, Modeling Practical Performance Limits of Photoelectrochem-ical Water Splitting Based on the Current State of Materials Research,ChemSusChem 7 (5) (2014) 1372–1385. doi:10.1002/cssc.201301030 .[60] R. E. Rocheleau, E. L. Miller, Photoelectrochemical Production of Hy-drogen: Engineering Loss Analysis, Int. J. Hydrogen Energ. 22 (8)(1997) 771 – 782.[61] C. C. L. McCrory, S. Jung, I. M. Ferrer, S. M. Chatman, J. C.Peters, T. F. Jaramillo, Benchmarking Hydrogen Evolving Reac-tion and Oxygen Evolving Reaction Electrocatalysts for Solar Wa-ter Splitting Devices, J. Am. Chem. Soc. 137 (13) (2015) 4347–4357. doi:10.1021/ja510442p .[62] C. R. Cox, J. Z. Lee, D. G. Nocera, T. Buonassisi, Ten-Percent Solar-to-Fuel Conversion with Nonprecious Materials, Proc. Natl. Acad. Sci.USA 111 (39) (2014) 14057–14061. doi:10.1073/pnas.1414290111doi:10.1073/pnas.1414290111