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

We studied ion concentration profiles and the charge density gradient caused by 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 profiles depend on not only ion transport but also the ionization of molecules. We considered the ionization of molecules and ion association in weak electrolytes and obtained analytical expressions for ion densities, electrostatic potential profiles, and ion currents. We found the case that the total ion density gradient was given by the Kuramoto length which characterized the distance over which an ion diffuses before association. The charge density gradient is characterized by the Debye length for 1:1 weak electrolytes. We discuss the role of these length scales for efficient water splitting reactions using photo-electrocatalytic electrodes.


Introduction
The Poisson-Nernst-Planck (PNP) equations have been used to describe a wide range of transport phenomena from electrons and holes in semiconductors to ions in electrolytes. [1][2][3][4][5][6][7][8][9] The conservation of charge for each ion species and electrostatic interactions among charge carriers are treated by the PNP equations in a self-consistent manner. The PNP equations take into account the drift currents due to the electric fields generated by the distribution of charge carriers. The effect is substantial e.g. when electrode reactions generate charge density gradients. The PNP equations can be applied to obtain concentration profiles and an electro-static potential generated by electrode reactions. Although the PNP equations are used for the study of coupled effects between electric fields and charge carrier transports, they are nonlinear and solved mainly by numerical methods. However, for certain cases, approximate analytical solutions have also been obtained for strong electrolytes. [2,3,[8][9][10][11][12][13] Using the PNP equations, it has been shown that the spatial dimensions of concentration gradient can be many orders of magnitude larger than the characteristic length scale of charge density profiles given by the Debye length in strong electrolytes. [1][2][3][4]13] In weak electrolytes, only a small fraction of molecules are ionized in bulk. Ion concentration profiles depend on ion transport and the ionization of molecules. It should be noted that depleted ions in the vicinity of the electrode can be replenished by the ionization of molecules and diffusion from the bulk in weak electrolytes. The current at the interface between the electrolyte and electrode can be much larger than that in bulk due to the molecule ionization.
In this paper, we study charge transport induced by electrode reactions in weak electrolytes using the e-PNP equations without assuming a priori charge neutrality. Photo-electrochemical (PEC) conversion of water can be an example, but the fundamental results apply to other electrode reactions. We show that the gradient of total ion density (the sum of cation and anion concentrations) is characterized by either the Kuramoto length [26][27][28][29][30][31]  Debye length depending on the situation, while the gradient of charge density (the difference between cation and anion concentrations) is characterized by the Debye length in binary monovalent weak electrolytes. The Kuramoto length characterizes the length scale of local density fluctuations around a uniform concentration state. [26][27][28][29][30][31] In weak electrolytes, the Kuramoto length is given by the length scale of diffusive migration of ions within its life-time; the life-time is determined by the association rate of ions. [26][27][28][29][30][31] Our results indicate that the ion density gradients can also be characterized by the Kuramoto length when both cations and anions are discharged at the electrode in binary monovalent weak electrolytes. For this case, ion density drop caused by electrode reactions is recovered by ion density fluctuations localized within the Kuramoto length. We discuss the efficiency of water splitting reactions using photo-electrocatalytic electrodes in terms of the Kuramoto length and the Debye length. We also discuss the overpotential related to a charge density gradient near the electrode and show that it can be reduced when both cations and anions are discharged at the electrode.

Theory
As shown in Fig. 1, we considered the thermal motion of cation and anion under Coulombic interactions. The one-dimensional (1D) x-coordinate was introduced by assuming a uniform distribution of ions in the plane parallel to the electrode surface. The origin of the x-coordinate was the electrode surface, and the x-coordinate was normal to the electrode surface. The concentration and diffusion constant of cations are denoted as n + and D + , respectively. The concentration and diffusion constant of anions are denoted as n − and D − , respectively, while the electric field is denoted as E. Each species moves by electromigration and diffusion. For simplicity, we consider a binary monovalent electrolyte (1:1 electrolyte). The concentration and the diffusion constant of the undissociated neutral compound are denoted as m and D m , respectively. At a sufficiently large distance L away from the electrode, the electrolyte is neutral. The external electric field is not applied. We consider ion discharge reaction at the electrode followed by product formation. Products of electrode reactions are assumed to be removed. An example is hydrogen evolution reaction by photo-electrocatalytic (PEC) water splitting. [32,33] For PEC water splitting, cations, anions, and undissociated molecules correspond to H 3 O + , OH − , and water molecules, respectively. We mainly consider the case that one of the ion species is discharged at the electrode and the other species are inert and reflected at the electrode as described in Fig. 1 (a). In Sec. 6, we consider the special case that both cation and anion species are discharged at the electrode as described in Fig.  1 (b). The situation corresponds to overall water splitting into H 2 and O 2 using a single photocatalyst such as GaN-ZnO and ZnGeN 2 -ZnO under light illumination. [34,35] For simplicity, we assume that reactive sites on the electrode are structureless. In experiments, the electrode can be regarded as homogeneous when the electrode is modified with molecular or metal complex cocatalyst to realize both oxidation and reduction reactions on the electrode. Some electrodes with nano-particulate photocatalysts could also be regarded as homogeneous.
Ion transport was described by the Nernst-Planck equations and the electric field satisfied the Poisson equation. By taking into account the association and dissociation reactions of ions with diffusion and ion migration under the electric field E(x), the governing equations in the bulk phase are given by [5,[14][15][16][17][18][19][20][21][23][24][25] where k d and k a denote the dissociation rate constant and the association rate constant, 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 diffusive 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. When the electrostatic potential is controllable by potentiostat operation, Eqs. (1)-(4) may be more conveniently expressed by using the electrostatic potential rather than using the electric field. However, we study the ion currents and the electrostatic potential induced by the ion discharge reactions at the electrode without imposing the external electric field. In this case, the electrostatic potential is not a controllable parameter. To avoid a redundant integration constant associated with the electrostatic potential we express Eqs. (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 difference relative to the electrostatic potential at L can be defined by and can be calculated by solving Eqs. (1)-(4) under proper boundary conditions. The total potential difference between the electrode surface and bulk (overpotential) is V (0).

Boundary Conditions
At a sufficiently large distance away from the electrode surface, the ion concentrations are not affected by surface reactions. Correspondingly, we set boundary conditions at x = L and assumed equal concentrations of cations and anions given by n b . Boundary conditions at x = L are given by where k d m b = k a n 2 b is satisfied in the bulk. We consider the situation that the external electric field is absent in the bulk In ideal situation of photo-electrochemical (PEC) water splitting, the external electric field is not applied. We also assumed a metal electrode, where the net image charge in the electrode is the same as the net charge in the electrolyte system. Therefore, the electric fields at the boundaries enclosing both the real charges in the electrolyte and the image charges are zero. This situation can also be viewed from another perspective. Because we consider neutral products from a neutral electrolyte by electrode reactions, the net charge in the system including those induced in the electrode should be neutral. In this case, electric fields at some distance sufficiently far from the electrodes should be zero. [6] We consider the case where cations are reduced at the interface between the 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 expressed as Equation (12) indicates that undissociated molecules are reflected by the cathode. The left-hand side of (11) represents the concentration current of anions. Equation (11) indicates that anions are reflected by the cathode. In Eq. (10), n s represents the net effect of surface reactions on n + (0). As a plausible boundary condition for irreversible surface reactions, we considered that n s = 0 on the right-hand side of Eq. (10). The boundary condition is suitable to study the case when the overall reactions are limited by ion transport to the electrode surface. In general, we set n s being constant on the right-hand side of Eq. (10).

Dimensionless Equations
We can express concentrations of cations and anions in dimensionless units as For convenience, we introduce [12] Using the above definitions, the original diffusion constants can be obtained by D + =D/ (1 − ∆), and D − =D/ (1 + ∆). By using the Debye length given by [12,13,[36][37][38] we can express the characteristic electric field strength and diffusion time as [12] We also introduce dimensionless variables given by Finally, we express the reaction rate constants in the dimensionless form, Here,k a andk d represent the dimensionless ion association rate constant and ion dissociation rate constant, respectively. Equation (1)-(4) can be rewritten using the dimensionless variables as Similarly, boundary conditions at y = L/λ D are given by and boundary conditions at y = 0 [Eqs. (10)- (12)] are given by Finally, the boundary condition for the electric field given by Eq. (9) can be expressed as The potential difference given by Eq. (5) can be expressed as Below we consider the case where the condition is satisfied; under the condition and using k d m b = k a n 2 b , it can be shown that the right-hand side of Eq. (21) can be ignored in steady states. For this case, the dimensionless concentration of undissociated molecules can be approximated by N m = 1 from the solution of diffusion equation under the boundary conditions of Eqs. (25) and (28). This situation represents that the 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 compound can be approximated as constant. Obviously, Eq. (31) is satisfied for water. For general weak electrolytes, additional condition for the concentration of the undissociated neutral compound may be needed. By using the degree of ionization α, we have n b /m b = α/(1 − α), where the numerator and the denominator indicate the increase of the amount of cation and the decrease in the amount of undissociated molecules for the extent of ionization reaction, respectively. When we consider the dissociation of an initial amount n 0 of the undissociated compound, the equilibrium dissociation constant is expressed as K a = n 0 α 2 /(1 − α), where n 0 is divided by the standard concentration 1 mol/L. In weak electrolytes, n 0 typically satisfies the relation n 0 ≫ K a /4 and we obtain α ≈ K a /n 0 . By substituting n b /m b ≈ α ≈ K a /n 0 , Eq. (31) can be rewritten as K a /n 0 ≪ (D m /D) 2 /k 2 a . Using D m /D ∼ 10 as a typical values of the diffusion constants of Carboxylic acids [39] and 1 <k a < 10, we find K a /n 0 ≪ 10 −2 ∼ 10 −4 . For most weak electrolytes satisfying K a ≤ 10 −4 , the condition can be fulfilled under normal values of n 0 . Below, we assume that the condition is satisfied and put N m = 1 in Eqs. (19)- (20).

Analytical Results
When ion discharge reactions proceed at the electrode-electrolyte interface 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) by assuming N m = 1 as explained in the previous section. As explained later, the steady-state analytical solution is obtained under a certain condition relevant to weak electrolytes. A separate consideration is needed for strong electrolytes as shown in Appendix A.
Because we are interested in the deviations from equilibrium concentrations induced by surface reactions, we introduce sum and change as [12] S = N + + N − , C = N + − N − , where S and C characterize the total mass density and charge density, respectively. To investigate the deviations from equilibrium states, we introduce In steady states, Eqs. (19), (20) and (22) can be rewritten using these variables as where N m = 1 is set by assuming the condition Eq. (31) is satisfied. In a weak electrolyte, the dissociation of a small amount of the solute can be represented by the right-hand sides of Eqs. (33) and (34).
The boundary conditions at y = L/λ D can be written as The boundary conditions at y = 0 can be expressed as The boundary condition for the electric field is given by Eq. (29). By linearizing the above equations, we solved the following set of equations, By substituting Eq. (35) into Eq. (41), we obtain If the right-hand side is zero, Eq. (42) is equal to the equation derived by Debye and Falkenhagen in the steady state. [40] The right-hand side of Eq. (42) reflects ionization of molecules. When the ion concentrations are conserved, the Debye-Falkenhagen equation is obtained, and the length scale of the charge density variation is characterized by the Debye length. Equation (42) indicates that the charge density variation is still characterized by the Debye length even when the ion concentrations are non-conserved by ionization of molecules: the ionization effect is represented by the inhomogeneous term in Eq. (42) when the linearization approximation is used. The validity of 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 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 by Eq. (37) in the limit of L → ∞, we obtained from Eq. (42) 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 = −2(1 − N s ) up to the linear order in 1 − N s . By substituting δC and integrating Eq.
The boundary condition given by Eq. (29) can be expressed as ξ(y) → 0 in the limit of L → ∞. Therefore, we find that C ξ = 0. These results can be written as The electric field strength can be written as From Eq. (30) and taking the limit of L → ∞, the potential profile can be obtained as Using N + = (2 + δC)/2 and N − = (2 − δC)/2, we obtained the ionic concentration profiles as Exponential screening of Eqs. (49)-(50) reminds us of the exponential decrease of electric field formed near charged surfaces in strong electrolytes under local equilibrium [13,[36][37][38]41]; the solution for a certain boundary condition is known as the Debye-Hückel theory. However, the situation considered here is different from that of the Debye-Hückel theory. The Debye-Hückel theory is based on the Poisson-Boltzmann equation, where ion concentrations are assumed to obey local equilibrium. Here, we solved the e-PNP equations where the effect of ion currents is taken into account and deviations from the local equilibrium distributions can be studied. The dimensionless concentration current density in Eq. (26) was obtained from Eqs. (49) and (46) as The ion current density can be expressed as and its value at the electrode surface is given by using a linearized approximation.

Interpretation of the linearization condition
As mentioned above, the linearization condition of the e-PNP equations is given byk a > 1. The linearization condition can be interpreted in terms of two time scales using the definitionk a = λ 2 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. The time scale of diffusion over the Debye length is given by t 0 = λ 2 D /D; the Debye length characterizes the spatial extent of the deviation from charge neutrality.k a > 1 indicates that the time scale of bulk reactions is shorter than that of diffusion over the Debye length. Therefore, charge neutrality is mainly recovered by ion association and dissociation reactions whenk a > 1. In the opposite limit ofk a < 1, charge neutrality is mainly recovered by ion diffusion.
As shown below, the linearization condition is satisfied when ion association is diffusion limited. The diffusion-limited ion association rate can be relevant for weak electrolytes, where the intrinsic association rate constants are large, and the dissociation rate constants are small. For example, the association rate of water can be estimated as k a = 1.3 × 10 11 L/(mol s) using the diffusion limited bulk rate constant 4π(D + + D − )R eff together with the effective reaction radius given by R eff = |r c |/[1 − exp(−|r c |/R)], [42][43][44][45] where |r c | = e 2 /(ǫk B T ) is the Onsager radius and the diffusion constants of H 3 O + and OH − are 9.4 × 10 −9 m 2 /s and 5.3 × 10 −9 m 2 /s, respectively. [46] In the above estimation, we also used ǫ = 80 and R = 0.75 nm. [43,44] The estimated k a value is close to the literature value given by k a = 1.4 × 10 11 L/(mol s). [43,44,47] This indicates that the effect of diffusive interaction can be ignored owing to the low degree of ionization in weak electrolytes. [48] The similar estimation is also possible for other weak electrolytes. [49] By noticing that the Debye length can be expressed using the Onsager length as k a = λ 2 D k a n b /D can be expressed as Because the minimum value of ( for given values of ǫ and the reaction radius. For water, we findk a = 2.0 and the linearization conditionk a > 1 is satisfied. For diffusion-limited ion association, the linearization condition is rewritten as using Eq. (55) and is always satisfied. The diffusion-limited association is commonly adopted for weak electrolytes because the ion association rapidly proceeds due to the Coulombic interaction without screening; [50,51] the screening by an ionic atmosphere only occurs in strong electrolytes.
In calculating the diffusion-limited association rate, the concentration of counterions at the distance orders of magnitude larger than the Onsager radius is assumed to be homogeneous. [42,43,50] On the other hand, the density profiles obtained by solving the 1-dimensional e-PNP equations are characterized by the Debye length. If Onsager radius is orders of magnitude smaller than the Debye length, the combined approach could be utilized. The separation of length scales is satisfied when the relative dielectric constant is sufficiently high. For water, the Onsager length is 0.7 nm which is orders of magnitude smaller than the Debye length given by 1 µm.
Regarding the dissociation rate constant k d , we did not take into account the effect of electric fields on the ion dissociation rate. Strictly speaking, the dissociation rates of geminate ions have stronger electric field dependence compared to the association rates in weak electrolytes. [50,[52][53][54] However, the field dependence of the dissociation rates is negligibly small in the absence of external fields because the largest internal electrostatic potential difference in the spatial extent of the Debye length is on the same scale as thermal energy at room temperature. The linearization condition will be numerically justified in the next section.

Numerical justification of the linearization condition
In this section, we confirm the linearization condition numerically. In a steady state, we numerically solved Eqs. (19)-(29) by the relaxation method. We present the numerical results obtained using the boundary condition including the effect of the electrode.
As a concrete example, we studied the photo-induced electrochemical reduction of water at electrode surfaces yielding hydrogen. Water is a weak electrolyte, and the dissociation of a small amount of water can be represented by the association and dissociation reaction terms in Eqs. (33) and (34). The diffusion constants of H 3 O + and OH − are 9.4 × 10 −9 m 2 /s and 5.3 × 10 −9 m 2 /s, respectively. [46] By substituting D + = 9.4 × 10 −9 m 2 /s and D − = 5.3 × 10 −9 m 2 /s into Eq. (14), we obtainedD = 6.78 × 10 −9 m 2 /s and ∆ = 0.28. Using n b = 1 × 10 −7 mol/L, the Debye length was obtained as λ D = 9.8 × 10 −7 ∼ 10 −6 m. The linearization condition of the e-PNP equations is given byk 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 Debye length given byD/λ 2 D . In this section, we study the case whenk a = 2.0 andk a = 0.02; the linearized results shown in Sec. 3 could be applicable to the case ofk a = 2.0 but not to the case ofk a = 0.02. For brevity, the case ofk a = 2.0 is classified as weak electrolytes and that ofk a = 0.02 is classified as strong electrolytes. In weak electrolytes, the intrinsic association rate constants are large, and the dissociation rate constants are small. The linearization condition is obtained when the association rate is diffusion limited because the intrinsic association rate constant is large. The linearization condition is most likely satisfied for weak electrolytes. On the contrary, the association rate constants of strong electrolytes are small and can be ignored.  In Fig. 2, we compare the results of the Nernst equation, the Henderson-Planck equation, and Eq. (48). The electrostatic potential profile of the Nernst equation can be expressed as which is equivalent to the local equilibrium assumption given by where ion mass flows are assumed to be absent. Intuitively, the Nernst equation could be thought of as applicable to weak electrolytes when the ion discharge reaction at the electrode is slow so that the resultant deviation from charge neutrality near the surface is recovered by bulk reactions rather than ion diffusion. Indeed, we can show that the electrostatic potential profile given by Eq. (48) is obtained by linearizing Eq. (58) by assuming [n b − n + (x)]/n b ≪ 1 as follows: However, it should be reminded that Eq. (61) is obtained from the e-PNP equations when the linearization conditionk a > 1 is satisfied regardless of the condition on [n b −n + (x)]/n b . In other words, ifk a > 1 is the linearization condition in the presence of diffusion, Eq. (61) should be satisfied even when The Nernst equation is obtained by assuming that ion mass flows are absent. The situation is irrelevant for strong electrolytes because electrode reactions induce ion mass flows in particular in the absence of ionization and association of ions. In strong electrolytes, ions are always dissociated and the terms that represent association and dissociation reactions on the right-hand sides of Eqs. (1) and (2) should be absent. As mentioned above, linearization fails for strong electrolytes, wherek a < 1. In the absence of association and dissociation reactions, approximate analytical solutions of the Poisson-Nernst-Planck equations have been studied more than decades. As shown in Appendix A and B, the potential difference known as the diffusion potential is obtained, [3,[10][11][12] Equation (62)   In Fig. 3, we show concentrations of cations and anions as a function of distance from the electrode surface in strong and weak electrolytes. In weak electrolytes, both cation and anion concentrations deviate from the bulk concentration when the distance from the electrode is within a few Debye lengths. The result is consistent with the previous result obtained by assuming weak external electric field. [25] Concentration profiles of cations are long-ranged in strong electrolytes. Interestingly, the cation and anion concentrations were different when the distance from the electrode was within a few Debye lengths in both strong and weak electrolytes. In other words, charge density caused by breakdown of charge neutrality was observed approximately within the surface layer to a depth of a few Debye lengths. The anion concentration profile is similar to that of cations outside this layer in strong electrolytes due to charge neutrality. The difference between the spatial range of charge density and that of mass density for strong electrolytes has been pointed out. [4] Although small, the increase of anion concentration near the electrode surface seen in Fig. 3 could be caused by ionization of a few unionized molecules present in the strong electrolyte characterized byk a = 0.02. We conclude that the linearization condition of the e-PNP equations is given byk a > 1.

Both cations and anions are discharged at the electrode surface
As indicated in Eq. (43), the total charge density gradient could be scaled by the factor 1/ 2k a from the Debye length when the total charge density drops by allowing anions to be discharged as well as cations at the electrode surface as shown in Fig. 1 (b). In this section, we consider such situation.
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 written as 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. conditions, we find 1). For the ideal case of N s = N s− = 0, we obtain, The factor 2k a indicates the length scale given by the inverse of 2k a n b /D. For diffusion limited ion association in weak electrolytes,k a > 1 holds and we can estimate that D /(2k a n b ) might be smaller than the Debye length becausek a > 1 can be rewritten as λ D > D /(k a n b ). The concentration gradients are essentially characterized by the length scale given by D /(2k a n b ) rather than the Debye length when both cations and anions are discharged at the electrode. Interestingly, the electrostatic potential difference is essentially characterized by the Debye length because 2k a > 1 holds in Eq. (68). The result is reasonable because the electrostatic potential difference originates from the capacitance effect caused by the charge density near the electrode; the charge density gradients are characterized by the Debye length.
As shown above, the length scale given by D /(2k a n b ) characterizes the ion concentration gradients. In an entirely different context but using the similar method, it was found that density profiles could be characterized by the association reaction-diffusion length scale in certain multi-ionic solutions. [49] D /(2k a n b ) represents the length scale of the diffusive migration within the life-time of ions which will disappear as a consequence of the ion association reaction. The length scale is called the Kuramoto length. [26][27][28][29][30] The Kuramoto length characterizes the spatial correlation of density fluctuations around a uniform concentration state. The density fluctuations beyond the length scale given by the Kuramoto length are independent. Here, we show that the ion density drop at the electrode surface can be localized within the Kuramoto length. The Kuramoto length can be expressed using the diffusion-controlled rate shown in Sec. 4 as In water, the Kuramoto length is given by λ D /2 ≈ 0.5 µm. We compare the above results with the numerical solutions as in the previous section. As an example, we consider photo-induced water splitting reactions by photocatalysts. We setk a = 2.0 and the other parameters are the same as those in the previous section. In Fig. 5, we show the cation and anion concentration profiles as a function of distance from the electrode. We can see that the above analytical solutions satisfactory reproduce the numerical results. Note that the concentration gradients are short ranged when both cations and anions are reactive at the electrode surface compared to the case when only cations are reactive. The difference reflects the fact that the concentration gradients are characterized by the Kuramoto length for the former and the Debye length for the latter. The steeper concentration gradients generated the larger ion currents by diffusion. In this respect, as long as side effects such as back reactions can be ignored, the photo-catalysts that are reactive to both cations and anions are desirable.
In Fig. 5, we show electrostatic potential as a function of distance from the electrode. Again, the analytical expressions reproduce the numerical results. When only cations are reduced, and anions are not reactive at the electrode, the electrode potential is higher than the electrostatic potential in the bulk water, and the result can be interpreted as overpotential. When both cations and anions are reactive at the electrode, the difference in the electrode potential is reduced and is proportional to the difference between the cation diffusion constant and the anion diffusion constant as shown in Eq. (68). Depending on the difference between the diffusion constants, the electrostatic potential difference can be negative. For any case, the electrostatic potential gradient is characterized by the same length scale given by the Debye length. The situation is very different from the ion concentration profiles, where the length scale is characterized by the Kuramoto length when both cations and anions are discharged.

Conclusion
The results of this paper are related to photo-electrochemical (PEC) water splitting. Although theoretically estimated PEC conversion efficiency exceeds 10%, [55][56][57][58][59][60] such high efficiency has not yet been achieved without an external energy supply. [55,61,62] Unlike photovoltaics, the PEC conversion efficiency can be limited by electrochemical processes occurring at the interface between an electrode and electrolyte even if a semiconductor with a suitable band gap is used. There could be two main limiting factors for PEC conversion efficiency related to the electrode-electrolyte interface.
The first limitation of PEC conversion efficiency could be that ion mobilities can be lower than electron mobilities in semiconductors. Steep ion concentration gradients would be preferable to induce ion currents to resolve the inhomogeneous ion distribution. We have shown that the concentration gradients are characterized by the Debye length when cations are reduced at the electrode, and the anions are inert for 1:1 monovalent weak electrolytes. The concentration gradient can be approximately given by n b /λ D . By using D + ≈ 10 −8 m 2 /s, λ D ≈ 10 −6 m, the ion current density can be expressed as when the boundary condition is given by n + (0) = 0 corresponding to the case where the back reaction can be ignored because of fast forward reactions at the interface, where 96500 C/mol is approximated by 10 5 C/mol. The above result matches with that of linearization given by Eq. (53) and is 30% smaller than the exact numerical result. We have also shown that the concentration gradient is characterized by the Kuramoto length when both cation reduction and anion oxidation reactions occur at the electrode for 1:1 monovalent weak electrolytes. The Kuramoto length is about the half of the Debye length. As a result, the ion current can be twice larger than that obtained when anions are not oxidized at the electrode surface. However, to realize 10% PEC conversion efficiency under 1 sun, a current density of 10 mA cm −2 is required. This result indicates that ion currents may limit the overall efficiency because of the small diffusion constant of H 3 O + compared with that of electrons. For photo-induced splitting of pure water, recent experimental setups have used internal convection caused by placing cathodic and anodic electrodes in parallel configurations. [32,33] In this configuration, an extra energy source to generate convection was not needed. In principle, convective flows can be considered by introducing the hydrodynamic equation into the Nernst-Planck equations. [63] An extension of the present work to include convective flows will be important to study the case when the current density is high.
The second limitation could be overpotential related to a charge density gradient near the electrode. Our results indicate that the overpotential is at most on the same scale as thermal energy at room temperature and can be further reduced when both cation reduction and anion oxidation reactions occur at the electrode for 1:1 monovalent weak electrolytes. The overpotential is caused by a charge density gradient whose spatial extent is characterized by the Debye length. It should be noted that the charge density gradient can be estimated from the pH gradient only for pure water. In general, charge density cannot be estimated from pH alone in the presence of other ion species.
For simplicity, we considered the case that molecules could be ionized to monovalent ions in the absence of buffer electrolytes composed of weak acid and/or weak base. We also did not consider supporting electrolytes to increase the conductivity of the bulk solution. For photo-induced watersplitting reactions, multi-ionic solutions are commonly used to control electrostatic potentials and engineer the pH of electrolytes. [20,21,64] There, large-scale pH gradients were observed. The observed large-scale pH gradients could be related to the length-scale separation between charge density profiles and concentration profiles, which could be possible for multi-ionic solutions. The effect of multiple ionic species on ion currents could be taken into account using new numerical methods developed to calculate junction potentials, [11,65] where the charge neutrality condition is unnecessary. The analytical method applied in this paper should be extended to multi-ionic solutions in the future.
(19)- (22) can be rewritten using δS = N + + N − − 2, δC = N + − N − , and z by applying the chain rule for partials as The diffusion may give rise to algebraic τ dependence, The boundary condition at y = 0 is expressed as The potential difference can be obtained from where A 1 is an integration constant. By integrating the above equation, we have where A 0 is an integration constant. We note that where we used erf(z/2) = (1/ √ π) z 0 dz 1 exp (−z 2 1 /4). [66] Note that lim 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) By further integration, we find where we used lim z→∞ δS (0) = 0 and A 3 is an integration constant. We consider the case that the electric field at the infinite distance from the electrode is zero. The L → ∞ limit of z is given by z → ∞ because we have z = y/ √ τ , y = x/λ D , and x = L. Noting that the dimensionless field is denoted by η (1) , we obtained Although the concentration of cations should be lower than that of anions near the electrode because of the reduction of H 3 O + at the electrode surface, the charge neutrality condition expressed by Eq. (A.7) was obtained as the lowest-order solution from the expansion of 1/τ . Conventionally, Eq. (A.24) has been derived using the difference in chemical potentials at the boundaries of liquid junctions. [12] In the derivation in Appendix A, we used different boundary conditions, and the result indicates that the voltage change caused by electrolytic reactions at an electrode surface can be approximated by the Henderson-Planck equation in strong electrolytes.

Appendix B. The electrostatic potential under charge neutrality condition
Interestingly, the diffusion potential can be derived by assuming charge neutrality regardless of presence of association and dissociation reactions. When the charges do not accumulate at the interface between the electrode and electrolyte, the positive and negative currents should be equal. The equality of the positive and negative currents in the steady state can be expressed as