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

TL;DR

This study explores how the Kuramoto length influences ion concentration and charge density profiles in weak electrolytes during photo-catalytic water splitting, using Poisson--Nernst--Planck equations without assuming charge neutrality.

Contribution

The paper introduces analytical expressions linking ionization, ion association, and length scales like Kuramoto and Debye lengths in weak electrolytes.

Findings

Ion concentration profiles depend on ion transport and ionization.

Kuramoto length characterizes the ion density gradient.

Charge density gradient is described by the Debye length.

Abstract

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.

Equations

$$\frac{\partial }{\partial t}{n}_{+}$$

$$\frac{\partial }{\partial t}{n}_{-}$$

$$\frac{\partial }{\partial t}m$$

$$\frac{\partial }{\partial x}E(x)$$

$$V(x)=-{\int }_L^{x}d{x}_{1}E(x_1),$$

$${n_{+}|}_{x=L}=n_b,$$

$${n_{-}|}_{x=L}=n_b,$$

$${m|}_{x=L}=m_b,$$

$$E(L)=0.$$

$$n_{+}(0)$$

$${D_{-}\left(\frac{\partial n_{-}}{\partial x}+\frac{eE{n}_{-}}{k_{\mathrm{B}}T}\right)|}_{x=0}$$

$${\frac{\partial m}{\partial x}|}_{x=0}$$

$$N_{+}=\frac{n_{+}}{n_b},N_{-}=\frac{n_{-}}{n_b},N_s=\frac{n_s}{n_b}\text{mbox}and{N}_m=\frac{m}{m_b}.$$

$$\bar{D}$$

$${\lambda }_D={\left(\frac{ϵk_{\mathrm{B}}T}{8\pi e^2{n}_b}\right)}^{1/2},$$

$$E_0=\frac{k_{\mathrm{B}}T}{e{\lambda }_D}={\left(\frac{8\pi n_b{k}_{\mathrm{B}}T}{ϵ}\right)}^{1/2},t_0=\frac{{\lambda }_D^2}{\bar{D}}=\frac{ϵk_{\mathrm{B}}T}{8\pi e^2{n}_b\bar{D}}.$$

$$y=\frac{x}{{\lambda }_D},\tau =\frac{t}{t_0},\xi =\frac{E}{E_0}.$$

$${\bar{k}}_a=\frac{{\lambda }_D^2{k}_a{n}_b}{\bar{D}},{\bar{k}}_d=\frac{{\lambda }_D^2{k}_d}{\bar{D}}.$$

$$\frac{{\partial }^2{N}_{+}}{\partial y^2}-\frac{\partial }{\partial y}\left(\xi N_{+}\right)$$

$$\frac{{\partial }^2{N}_{-}}{\partial y^2}+\frac{\partial }{\partial y}\left(\xi N_{-}\right)$$

$$\frac{{\partial }^2{N}_m}{\partial y^2}$$

$$\frac{\partial }{\partial y}\xi$$

$${N_{+}|}_{y=L/{\lambda }_D}=1,$$

$${N_{-}|}_{y=L/{\lambda }_D}=1,$$

$${N_m|}_{y=L/{\lambda }_D}=1,$$

$$N_{+}(0)$$

$${\left(\frac{\partial N_{-}}{\partial y}+\xi N_{-}\right)|}_{y=0}$$

$${\frac{\partial N_m}{\partial y}|}_{y=0}$$

$${\xi |}_{y=L/{\lambda }_D}=0.$$

$$V(y)=-E_0{\lambda }_D{\int }_{y=L/{\lambda }_D}^{y}dy\xi (y).$$

$${\bar{k}}_a\frac{n_b}{m_b}\frac{\bar{D}}{D_m}\ll 1$$

$$\delta S=S-2,\delta C=C.$$

$$\frac{{\partial }^2\delta S}{\partial y^2}-\frac{\partial }{\partial y}\left(\xi \delta C\right)$$

$$\frac{{\partial }^2\delta C}{\partial y^2}-\frac{\partial }{\partial y}\left[\xi (\delta S+2)\right]$$

$$\frac{\partial }{\partial y}\xi$$

$${\delta S|}_{y=L/{\lambda }_D}=0,$$

$${\delta C|}_{y=L/{\lambda }_D}=0.$$

$${\frac{\delta S+\delta C}{2}|}_{y=0}$$

$${\frac{\partial }{\partial y}(\delta S-\delta C)+\xi (2+\delta S-\delta C)|}_{y=0}$$

Full text

1 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, 10, 3, 11, 12, 13, 8, 9] 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, 4, 3, 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.

The Nernst-Planck equation is the continuity equation representing the conservation of charge for each ion species, where dissociation and association of ions are not considered. For weak electrolytes, the PNP equations have been extended to describe ionization and recombination of ions. [5, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 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 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] or the 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.

2 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 H3O+, 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 H2 and O2 using a single photocatalyst such as GaN-ZnO and ZnGeN2-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]

[TABLE]

where $k_{d}$ and $k_{a}$ denote the dissociation rate constant and the association rate constant, respectively. Equation (4) is the Poisson equation. $-D_{+}(\partial n_{+})/(\partial x)$ and $D_{+}eEn_{+}/(k_{\rm 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_{\rm 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

[TABLE]

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)$.

2.1 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

[TABLE]

where $k_{d}m_{b}=k_{a}n_{b}^{2}$ is satisfied in the bulk. We consider the situation that the external electric field is absent in the bulk

[TABLE]

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

[TABLE]

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).

2.2 Dimensionless Equations

We can express concentrations of cations and anions in dimensionless units as

[TABLE]

For convenience, we introduce [12]

[TABLE]

Using the above definitions, the original diffusion constants can be obtained by $D_{+}=\bar{D}/\left(1-\Delta\right)$, and $D_{-}=\bar{D}/\left(1+\Delta\right)$. By using the Debye length given by [36, 12, 37, 13, 38]

[TABLE]

we can express the characteristic electric field strength and diffusion time as [12]

[TABLE]

We also introduce dimensionless variables given by

[TABLE]

Finally, we express the reaction rate constants in the dimensionless form,

[TABLE]

Here, $\bar{k}_{a}$ and $\bar{k}_{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

[TABLE]

Similarly, boundary conditions at $y=L/\lambda_{D}$ are given by

[TABLE]

and boundary conditions at $y=0$ [Eqs. (10)-(12)] are given by

[TABLE]

Finally, the boundary condition for the electric field given by Eq. (9) can be expressed as

[TABLE]

The potential difference given by Eq. (5) can be expressed as

[TABLE]

Below we consider the case where the condition

[TABLE]

is satisfied; under the condition and using $k_{d}m_{b}=k_{a}n_{b}^{2}$, 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 $\alpha$, we have $n_{b}/m_{b}=\alpha/(1-\alpha)$, 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}\alpha^{2}/(1-\alpha)$, where $n_{0}$ is divided by the standard concentration 1 mol/L. In weak electrolytes, $n_{0}$ typically satisfies the relation $n_{0}\gg K_{a}/4$ and we obtain $\alpha\approx\sqrt{K_{a}/n_{0}}$. By substituting $n_{b}/m_{b}\approx\alpha\approx\sqrt{K_{a}/n_{0}}$, Eq. (31) can be rewritten as $K_{a}/n_{0}\ll(D_{m}/\bar{D})^{2}/\bar{k}_{a}^{2}$. Using $D_{m}/\bar{D}\sim 10$ as a typical values of the diffusion constants of Carboxylic acids [39] and $1<\bar{k}_{a}<10$, we find $K_{a}/n_{0}\ll 10^{-2}\sim 10^{-4}$. For most weak electrolytes satisfying $K_{a}\leq 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).

3 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

[TABLE]

In steady states, Eqs. (19), (20) and (22) can be rewritten using these variables as

[TABLE]

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/\lambda_{D}$ can be written as

[TABLE]

The boundary conditions at $y=0$ can be expressed as

[TABLE]

The boundary condition for the electric field is given by Eq. (29).

By linearizing the above equations, we solved the following set of equations,

[TABLE]

By substituting Eq. (35) into Eq. (41), we obtain

[TABLE]

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

[TABLE]

using the boundary condition given by Eq. (36) in the limit of $L\rightarrow\infty$, where $C_{s}$ is an integration constant. Using the boundary condition given by Eq. (37) in the limit of $L\rightarrow\infty$, we obtained from Eq. (42)

[TABLE]

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 $\delta C$ and integrating Eq. (35), we found that $\xi(y)=(1-N_{s})\exp\left(-y\right)+C_{\xi}$, where $C_{\xi}$ is an integration constant. The boundary condition given by Eq. (29) can be expressed as $\xi(y)\rightarrow 0$ in the limit of $L\rightarrow\infty$. Therefore, we find that $C_{\xi}=0$. These results can be written as

[TABLE]

The electric field strength can be written as

[TABLE]

From Eq. (30) and taking the limit of $L\rightarrow\infty$, the potential profile can be obtained as

[TABLE]

Using $N_{+}=(2+\delta C)/2$ and $N_{-}=(2-\delta C)/2$, we obtained the ionic concentration profiles as

[TABLE]

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 [36, 37, 13, 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

[TABLE]

The ion current density can be expressed as

[TABLE]

and its value at the electrode surface is given by

[TABLE]

using a linearized approximation.

The above results were derived by linearizing the Poisson–Nernst–Planck equations with the association and dissociation reactions. The linearization is valid if $|\partial\left(\xi\delta C\right)/\partial y|$ on the left-hand side of Eq. (33) is smaller than $2\bar{k}_{a}|\delta S|$ on the right-hand side. The linearized solutions indicate $\delta C\partial\xi/\partial y=\xi\partial\delta C/\partial y$ and $|\partial\left(\xi\delta C\right)/\partial y|=2\delta C\partial\xi/\partial y$. The magnitude of charge density $\delta C$ would decrease by increasing the distance from the electrode surface and $\delta C\partial\xi/\partial y=\delta C^{2}/2\approx 2\left(1-N_{s}\right)^{2}$ near the electrode surface. If the electrode is perfectly reactive, we have $\delta C\partial\xi/\partial y\approx 2$ and $|\partial\left(\xi\delta C\right)/\partial y|=2\delta C\partial\xi/\partial y\approx 4$ near the electrode surface. By considering $\delta S\approx N_{-}\approx|\delta C|$, we also obtain $\delta S\approx 2$ near the electrode surface if the electrode is perfectly reactive. The condition for linearization given by $|\partial\left(\xi\delta C\right)/\partial y|\ll 2\bar{k}_{a}|\delta S|$ can be expressed as $\bar{k}_{a}>1$. Otherwise, the coupling between $\xi$ and $\delta C$ in Eq. (33) cannot be ignored. 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 equations is given by $\bar{k}_{a}>1$. The linearization condition can be interpreted in terms of two time scales using the definition $\bar{k}_{a}=\lambda_{D}^{2}k_{a}n_{b}/\bar{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}=\lambda_{D}^{2}/\bar{D}$; the Debye length characterizes the spatial extent of the deviation from charge neutrality. $\bar{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 when $\bar{k}_{a}>1$. In the opposite limit of $\bar{k}_{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\times 10^{11}\mbox{L/(mol s)}$ using the diffusion limited bulk rate constant $4\pi(D_{+}+D_{-})R_{\rm eff}$ together with the effective reaction radius given by $R_{\rm eff}=|r_{c}|/[1-\exp(-|r_{c}|/R)]$, [42, 43, 44, 45] where $|r_{c}|=e^{2}/(\epsilon k_{\rm B}T)$ is the Onsager radius and the diffusion constants of H3O+ and OH- are $9.4\times 10^{-9}$ m2/s and $5.3\times 10^{-9}$ m2/s, respectively. [46] In the above estimation, we also used $\epsilon=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\times 10^{11}\mbox{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

[TABLE]

$\bar{k}_{a}=\lambda_{D}^{2}k_{a}n_{b}/\bar{D}$ can be expressed as

[TABLE]

Because the minimum value of $(D_{+}+D_{-})^{2}/(D_{+}D_{-})$ is obtained when $D_{+}=D_{-}$, we find

[TABLE]

for given values of $\epsilon$ and the reaction radius. For water, we find $\bar{k}_{a}=2.0$ and the linearization condition $\bar{k}_{a}>1$ is satisfied. For diffusion-limited ion association, the linearization condition is rewritten as

[TABLE]

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$ $\mu$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.

5 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 H3O+ and OH- are $9.4\times 10^{-9}$ m2/s and $5.3\times 10^{-9}$ m2/s, respectively. [46] By substituting $D_{+}=9.4\times 10^{-9}$ m2/s and $D_{-}=5.3\times 10^{-9}$ m2/s into Eq. (14), we obtained $\bar{D}=6.78\times 10^{-9}$ m2/s and $\Delta=0.28$. Using $n_{b}=1\times 10^{-7}\mbox{mol/L}$, the Debye length was obtained as $\lambda_{D}=9.8\times 10^{-7}\sim 10^{-6}$ m. The linearization condition of the e-PNP equations is given by $\bar{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 Debye length given by $\bar{D}/\lambda_{D}^{2}$. In this section, we study the case when $\bar{k}_{a}=2.0$ and $\bar{k}_{a}=0.02$; the linearized results shown in Sec. 3 could be applicable to the case of $\bar{k}_{a}=2.0$ but not to the case of $\bar{k}_{a}=0.02$. For brevity, the case of $\bar{k}_{a}=2.0$ is classified as weak electrolytes and that of $\bar{k}_{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

[TABLE]

which is equivalent to the local equilibrium assumption given by

[TABLE]

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}\ll 1$ as follows:

[TABLE]

However, it should be reminded that Eq. (61) is obtained from the e-PNP equations when the linearization condition $\bar{k}_{a}>1$ is satisfied regardless of the condition on $[n_{b}-n_{+}(x)]/n_{b}$. In other words, if $\bar{k}_{a}>1$ is the linearization condition in the presence of diffusion, Eq. (61) should be satisfied even when $[n_{b}-n_{+}(x)]/n_{b}\sim 1$.

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, where $\bar{k}_{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, [10, 3, 11, 12]

[TABLE]

Equation (62) is also called the Henderson–Planck equation. In strong electrolytes, the potential difference is reduced from the Nernst equation by the factor given by $(D_{+}-D_{-})/(D_{+}+D_{-})$. The Henderson–Planck equation represents the diffusion potential caused by different diffusion constants of anions and cations. [10, 3, 11, 12]

Coming back to Fig. 2, we notice that the result of the Nernst equation is closer to the numerical result than that of the Henderson–Planck equation when $\bar{k}_{a}=2.0$ (weak electrolytes). However, the result of the Nernst equation deviates from the numerical result near the electrode. The potential profile obtained from Eq. (48) [or equivalently Eq. (61)] is close to the numerical result for all distances. Equation (48) was obtained by linearization of full transport equations and the Poisson equation. This result is valid under the linearization condition given by $\bar{k}_{a}>1$. The Nernst equation was obtained by assuming local equilibrium, where mass currents were not taken into account. The deviation found near the electrode could be caused by mass currents. When charge neutrality is broken near the electrode surface, both ionization and ion transport occur to recover charge neutrality. Because of ion transport, the local equilibrium assumption used to obtain the Nernst equation is violated, particularly near the electrode surface. The Henderson–Planck equation is derived for strong electrolytes, and its results cannot be applied to weak electrolytes.

When $\bar{k}_{a}=0.02$, the linearization condition $\bar{k}_{a}>1$ is not satisfied and the numerical results confirm that the results of Eq. (48) are not valid for this case. For strong electrolytes indicated by $\bar{k}_{a}=0.02$, the result of the Henderson–Planck equation is closer to the numerical result than that of the Nernst equation. The Henderson–Planck equation is derived for strong electrolytes, where ionization is not taken into account. In contrast, the Nernst equation is derived by assuming local equilibrium, which can be violated by mass currents. We should note that breakdown of charge neutrality can be recovered only by mass currents in strong electrolytes.

For $\bar{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 are consistent with the linearization condition given by $\bar{k}_{a}>1$.

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 by $\bar{k}_{a}=0.02$. We conclude that the linearization condition of the e-PNP equations is given by $\bar{k}_{a}>1$.

6 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/\sqrt{2\bar{k}_{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.

The boundary conditions at $x=0$ ($y=0$) given by Eq. (11) changes into

[TABLE]

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

[TABLE]

where $N_{s-}=n_{s-}/n_{b}$ as $N_{s}$ in Eq. (13). Similarly, we have $\left.\delta S+\delta C+2(1-N_{s})\right|_{y=0}=0$ from Eq. (38). By substituting Eqs. (43)-(44) into the above boundary conditions, we find $C_{s}=-2+N_{s}+N_{s-}$ and $C_{c}=N_{s}-N_{s-}-(2-N_{s}-N_{s-})2\bar{k}_{a}\Delta/(2\bar{k}_{a}-1)$. For the ideal case of $N_{s}=N_{s-}=0$, we obtain,

[TABLE]

The factor $\sqrt{2\bar{k}_{a}}$ indicates the length scale given by the inverse of $\sqrt{2k_{a}n_{b}/\bar{D}}$. For diffusion limited ion association in weak electrolytes, $\bar{k}_{a}>1$ holds and we can estimate that $\sqrt{\bar{D}/(2k_{a}n_{b})}$ might be smaller than the Debye length because $\bar{k}_{a}>1$ can be rewritten as $\lambda_{D}>\sqrt{\bar{D}/(k_{a}n_{b})}$. The concentration gradients are essentially characterized by the length scale given by $\sqrt{\bar{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 $2\bar{k}_{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 $\sqrt{\bar{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] $\sqrt{\bar{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

[TABLE]

In water, the Kuramoto length is given by $\lambda_{D}/2\approx 0.5$ $\mu$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 set $\bar{k}_{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.

7 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. [61, 55, 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}/\lambda_{D}$. By using $D_{+}\approx 10^{-8}$ m2/s, $\lambda_{D}\approx 10^{-6}$ m, the ion current density can be expressed as

[TABLE]

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 H3O+ 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 water-splitting reactions, multi-ionic solutions are commonly used to control electrostatic potentials and engineer the pH of electrolytes. [64, 21, 20] 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.

Acknowledgement

This work was supported by the “Research Project for Future Development: Artificial Photosynthetic Chemical Process (ARPChem)” (METI/NEDO, Japan: 2012-2022).

Appendix A Derivation of Eq. (62)

We performed a systematic expansion in $1/\tau$ when $\bar{k}_{a}=0$ by applying the approach developed to calculate liquid junction potential to our problem of surface reactions. [12] When concentration change is restored only by diffusion, a relevant variable is $z=y/\sqrt{\tau}$. In terms of the new variable, Eq. (19)–(22) can be rewritten using $\delta S=N_{+}+N_{-}-2$, $\delta C=N_{+}-N_{-}$, and $z$ by applying the chain rule for partials as

[TABLE]

The diffusion may give rise to algebraic $\tau$ dependence,

[TABLE]

$\eta^{(0)}(z)$ should be zero; otherwise, $\xi(z)\sim\sqrt{\tau}\eta^{(0)}(z)$ diverges as $\tau\rightarrow\infty$.

Now, we study the lowest-order solution given by the order of $1/\tau^{0}$. From Eq. (73), the charge neutrality condition is satisfied; i.e.,

[TABLE]

By substituting Eqs. (74) and (75) into Eqs. (71)-(73), we obtained $S^{(0)}(z)$ and $\eta^{(1)}(z)$ of the order of $1/\tau$ as

[TABLE]

Boundary conditions at $z=L/(\lambda_{D}\sqrt{\tau})$ are given by

[TABLE]

The boundary condition at $y=0$ is expressed as

[TABLE]

The potential difference can be obtained from

[TABLE]

and the total potential difference between the electrode surface and bulk is $V=V(0)$.

First, we integrated Eq. (78) and obtained

[TABLE]

where $A_{1}$ is an integration constant. By integrating the above equation, we have

[TABLE]

where $A_{0}$ is an integration constant.

We note that $\int_{0}^{\infty}dz_{1}\exp\left(-\frac{z_{1}^{2}}{4}\right)=\sqrt{\pi}$. Using the boundary conditions given by Eqs. (80) and (82), we obtained

[TABLE]

By substituting these equations into Eq. (85), it became

[TABLE]

where we used $\mbox{erf}(z/2)=\left(1/\sqrt{\pi}\right)\int_{0}^{z}dz_{1}\exp\left(-z_{1}^{2}/4\right)$. [66] Note that $\lim_{z\rightarrow\infty}\mbox{erf}(z/2)=1$ and we have $\lim_{z\rightarrow\infty}\delta S^{(0)}(z)=0$, as expected.

Finally, using Eq. (84), we obtained from Eq. (79)

[TABLE]

By further integration, we find

[TABLE]

where we used $\lim_{z\rightarrow\infty}\delta 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\rightarrow\infty$ limit of $z$ is given by $z\rightarrow\infty$ because we have $z=y/\sqrt{\tau}$, $y=x/\lambda_{D}$, and $x=L$. Noting that the dimensionless field is denoted by $\eta^{(1)}$, we obtained $A_{3}=0$ from Eqs. (90) and (9).

By substituting Eq. (90) into Eq. (83), the voltage difference was calculated as

[TABLE]

The total voltage difference between the electrode surface and bulk was obtained from $V(0)$ as

[TABLE]

Although the concentration of cations should be lower than that of anions near the electrode because of the reduction of H3O+ at the electrode surface, the charge neutrality condition expressed by Eq. (77) was obtained as the lowest-order solution from the expansion of $1/\tau$.

Conventionally, Eq. (94) 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

[TABLE]

and the charge neutrality condition is given by $n_{+}=n_{-}$. By introducing the charge neutrality condition, we obtain from Eq. (95) [12]

[TABLE]

and Eq. (62) known as the Henderson–Planck equation for 1:1 monovalent electrolytes even in the presence of association and dissociation reactions. The fact that Eq. (48) [or equivalently Eq. (61)] is not obtained by assuming local charge neutrality indicates the subtle aspects of charge neutrality to calculate the electrostatic potential.