We study the factors that affect the photoactivity of silicon electrodes for the water-splitting reaction using a self-consistent continuum solvation model of the solid-liquid interface. This model allows us to calculate the charge-voltage response, Schottky barriers, and surface stability of different terminations while accounting for the interactions between the charge-pinning centers at the surface and the depletion region of the semiconductor. We predict that the most stable oxidized surface does not have a favorable Schottky barrier, which further explains the low solar-to-hydrogen performance of passivated silicon electrodes.

## I. INTRODUCTION

Hydrogen is a sustainable energy carrier that can be produced by splitting water at the surface of a photocatalytic semiconductor.^{1–4} Yet, most of the semiconductors that are presently well developed for photovoltaics transfer poorly to photocatalysis.^{5–7} In the case of silicon, for instance, multiple factors contribute to the low solar-to-hydrogen efficiency. Notably, the electrochemical corrosion and surface restructuring of silicon electrodes are primary limitations to their use as photocatalysts. In fact, silicon is prone to oxidize in water, forming a passivating silica deposit at the solid-liquid interface.^{8,9}

While this process is known, the mechanisms that limit the photocatalytic efficiency have not been fully explored at the molecular level. The oxidation alters several properties of the interface, including the band gap, the surface states, and the Schottky barrier height. The latter is the potential difference between the bulk of the semiconductor and the surface; it provides the motive force for separating the photogenerated electrons and holes and transporting the excited charges from the electrode to the electrolyte.^{10–12} Therefore, the Schottky barrier is a central indicator of the light-harvesting ability of semiconductor electrodes.

First-principles methods, with their ability to probe atomic length scales, provide an ideal approach for examining the mechanisms that underlie the activity of photoelectrodes as a function of the surface termination. While several electronic-structure methods have been developed to understand neutral photocatalytic interfaces,^{13–16} until recently, no work has included the effects of applying a potential and controlling the hydrogen activity. Here, we exploit our newly developed semiconductor-continuum methodology^{17} to understand the charge-voltage response of different structures of the silicon-water interface. We extend this methodology to predict the stability and Schottky barriers of the various terminations, showing that while the most oxidized surface is typically the most stable, it also exhibits a low Schottky barrier, which helps understand the limited efficiency of silicon for photocatalytic hydrogen generation. Our work illustrates the capabilities of first-principles methods to identify the molecular factors controlling the photoactivity of semiconductor electrodes.

## II. BACKGROUND

To predict the influence of surface termination on charge separation, we first examine the microscopic mechanisms that lead to the formation of the Schottky barrier. We start by focusing on the interface between an intrinsic semiconductor and a chemically inert medium; in specific terms, we consider a semiconducting electrode in contact with an ideal electrolyte, i.e., one that does not interact chemically with the electrode and is stable over a wide range of applied voltages. Since the electrolyte is, in this ideal case, insensitive to the applied voltage, the equilibration of the system will take place without constraint on the Fermi energy. Consequently, no electronic charge will be injected or withdrawn from the semiconductor. Figure 1(a) depicts the resulting equilibrium state, where a surface dipole forms due to the reorientation of the solvent molecules. This surface dipole can be expressed as $(\delta \chi sm)\u25e6$ = $(\Delta s)\u25e6$ − $(\Delta m)\u25e6$, where $(\Delta s)\u25e6$ and $(\Delta m)\u25e6$ are the differences between the Fermi levels and average electrostatic potentials for the semiconductor and the surrounding medium,^{18} respectively.

Then, if defects are introduced in the semiconductor and chemically active ionic species are added to the embedding electrolyte, the chemical window of the solution will be reduced, causing the Fermi level of the electrode to be pinned by the chemical potential of the reaction that limits the stability of the reactive medium. This constraint leads to a different equilibrium state [Fig. 1(b)], where defect charge builds up within the depletion layer of the electrode and compensating ions accumulate within the double layer of the electrolyte. As a result, a large drop in the electrostatic potential is observed, corresponding to an increased surface dipole $\delta \chi sm$, a potential drop within the depletion layer of the electrode, namely, the Schottky barrier *φ*_{s}, and another electrostatic shift in the electrical double layer *φ*_{m}.^{19,20} In this equilibrated state, the Schottky barrier *φ*_{s} is simply related to the surface dipole through

where Δ_{s} and Δ_{m} are the differences between the Fermi level (divided by the electron charge *e*_{0}) and the bulk potential of the doped semiconductor and the chemically active medium, respectively. An electrostatic shift *φ*_{m} will also take place in the solvating medium. Since the ionic concentration of this medium is typically orders of magnitude greater than the carrier concentration of the semiconductor, the potential drop in the medium will be negligible: *φ*_{m} ≪ *φ*_{s}. (In practical experiments, *φ*_{m} can even be eliminated by varying the potential and pH of the solution to reach the “zeta potential” at which there is no interfacial charge and no accumulation of H^{+} and OH^{−} species.^{13,14})

Predicting the Schottky barrier of a semiconductor interface from first principles is challenging; since a barrier height of 0.5 V across a semiconductor with a dielectric constant of ∼10 and a carrier density of 10^{16} cm^{−3} would extend as far as ∼250 nm into the semiconductor, the direct quantum-mechanical modeling of Schottky barriers is prohibitively demanding. Hence, a common practice in the first-principles density-functional theory literature has been to compute Schottky barriers at electrically neutral surfaces with a simulation range on the order of ∼10 nm. In this approach, the neutral-interface dipole $(\delta \chi sm)\u25e6$ between an intrinsic semiconducting surface and an electrically neutral embedding medium is found and any remaining difference between the Fermi levels is equated to the Schottky barrier^{15,21–30}

This commonly used approach overestimates the barrier height, as can be seen by rearranging Eq. (1) into

In many cases, calculating only the first half of the right-hand side of Eq. (3), as done in this conventional approach, is a reasonable approximation; the overestimation of the Schottky barrier height may be compensated to some degree by the underestimated band gap within local and semilocal density-functional theory. Yet, this approach is applicable only when $\delta \chi sm\u2248(\delta \chi sm)\u25e6$, i.e., when the surface dipole of the charged interface under bias is close to that of the neutral interface, which is only true for interfaces that remain moderately polarized under bias. Therefore, there is a need for a computational model that would incorporate the effects of surface termination and external bias and would predict Schottky barrier heights accurately.

To fill this gap, we apply and further develop a first-principles approach that integrates an implicit semiclassical description of the bulk semiconductor with an explicit electronic-structure treatment of the semiconductor surface. This method, detailed in Sec. III, allows us to determine the potential-dependent interfacial dipole and the resulting equilibrium Schottky barriers.

## III. MODEL

### A. Quantum-continuum embedding

The proposed computational method consists of embedding an explicit quantum-mechanical description of the semiconductor surface layers into an implicit continuum model of the semiconducting and electrolytic environments. As described in Ref. 17, this procedure exploits the self-consistent continuum solvation (SCCS) approach,^{31} which introduces dielectric cavities around each facet of the slab. The dielectric permittivity is expressed as

where *ϵ*_{m} is the dielectric constant of the electrolytic medium and *ζ*(** r**) = (ln

*ρ*

_{max}− ln

*ρ*(

**))/(ln**

*r**ρ*

_{max}− ln

*ρ*

_{min}) is used as a smooth switching function, marking the transition between the quantum-mechanical and semiclassical continuum regions. Here,

*ρ*

_{min}and

*ρ*

_{max}serve as the density thresholds specifying the inner and outer isocontours of the dielectric cavity. The SCCS model also includes contributions from the external pressure, solvent surface tension, and solvent dispersion and repulsion. The surface tension is described by

*G*

_{cav}=

*γS*and the dispersion and repulsion effects by

*G*

_{dis+rep}=

*αS*+

*βV*. Here,

*γ*is the solvent surface tension, taken from experiment,

*α*and

*β*are the fitted parameters, and

*S*and

*V*are the quantum surface and volume of the solute, defined as

*S*=

*∫d*

**r**(

*d*Θ/

*dρ*)|∇

*ρ*| and

*V*=

*∫d*

**r**Θ(

*ρ*), where Θ is another smooth switching function, defined by Θ(

*ρ*) = (

*ϵ*

_{s}−

*ϵ*(

*ρ*))/(

*ϵ*

_{s}− 1). We utilize the parameterization of Andreussi

*et al.*, where

*ρ*

_{max}= 5 × 10

^{−3}a.u.,

*ρ*

_{min}= 1 × 10

^{−4}a.u.,

*γ*= 72.0 dyn/cm,

*α*= −22 dyn/cm, and

*β*= −0.35 GPa.

^{31}

A recent study has shown that the volume contribution *βV* to the total energy of the interface is physically inappropriate for a slab model, as it leads to a spurious dependence of this energy on the thickness size of the slab.^{32,33} However, because the simulated slabs are of comparable volume, the spurious term cancels out and does not affect the relative energies.

We also note that the parameterization of *ρ*_{max} and *ρ*_{min} was initially developed for neutral molecular systems in water. Adopting the alternative parameterization for cation species^{34} did not lead to significant changes in the final energy and charge-voltage curves as shown in the supplementary material. This is due to the fact that at a semiconductor surface, small changes in the surface charge lead to large potential shifts due to poor electrostatic screening, which implies that the voltage-dependent charge per atom remains small (on the order of a few tenths of the elementary charge), suggesting that the neutral-atom parameterization is more appropriate. The anion parameterization of Dupont *et al.*^{34} did lead to larger changes, but the qualitative trend remained the same. In the same vein, another parameterization has been developed by Hörmann *et al.*^{35} and tested by Nattino *et al.*, demonstrating improved agreement between the calculated differential capacitance and experimental data.^{36} As shown in the supplementary material, using the SCCS parameterization of Hörmann *et al.* does change the absolute voltages, but does not alter the voltage differences and thermodynamics trends significantly.

Using the selected parameterization, we first simulate an electrically neutral system with different dielectric responses on each side of the slab. These dielectric responses are described by the experimental permittivities of silicon (*ϵ*_{s} = 11.7) and water (*ϵ*_{m} = 78.3). By aligning the potential inside the solution to the electrostatic reference of the continuum solvent, the flatband potential *φ*_{FB} can be directly related to the Fermi level $EF\u25e6$ of the neutral surface through

where *e*_{0} denotes the elementary charge.

The charged electrode is then simulated by placing planar countercharges on each side of the slab and by assigning to the electrode a total charge *q*, split between the quantum-mechanical region *q*^{surf} and the bulk semiconductor *q*^{bulk} such that *q* = *q*^{bulk} + *q*^{surf}. (Here, the charge *q*^{bulk} is placed on the plane of the semiconductor side, while the charge −*q* is assigned to the Helmholtz plane in the solution.)

Finally, the long-range polarization of the depletion region is included by setting a frontier between the explicit and implicit semiconductor regions at a cutoff coordinate *z* = *z*_{c} along the transverse *z*-axis, typically two layers within the semiconductor electrode, corresponding to the inflection point of the potential profile where the net charge density vanishes (cf. Poisson’s equation). On the side of the cutoff bordering the solution, the electrode surface is described quantum mechanically, whereas on the other side, the Mott-Schottky model is used to describe the screening of the electrostatic field that arises from the dopant distribution (here, chosen to be *n*-type),

where *N* is the concentration of electron-donating defects, *k*_{B} is the Boltzmann constant, *T* is the ambient temperature, and *φ*_{0} is the asymptotic value of the potential in the bulk semiconductor. (A similar equation with opposite charge applies for *p*-doped semiconductors.) It is important to note that this Mott-Schottky model circumvents the need to reparameterize the SCCS model on the semiconductor side, since the embedding contribution from the continuum semiconductor only serves as an intermediate step that is ultimately replaced by the Mott-Schottky medium.

Using Poisson’s equation with a Boltzmann distribution of electronic charges in the depletion layer, we can then derive the Mott-Schottky asymptotic trend of the electrostatic potential in the continuum region,^{37}

Then, noting that the exponential term is typically small (i.e., *φ* − *φ*_{0} ≫ *k*_{B}*T*), the asymptotic value $\phi \xaf0$ of the difference between the electrostatic potential of the charged and neutral slabs $\phi \xaf$ can be obtained from

As shown in Fig. 2, the last step of the calculation consists of finding the charge distribution that self-consistently aligns the Fermi level of the explicit electrode surface ($EFsurf$) to that of the bulk material ($EFbulk$). The latter is calculated by adding the asymptotic potential difference $\phi \xaf0$ [Eq. (8)] to the Fermi level of the neutral slab,

Once the equilibrium Fermi level is known, the voltage *φ* = −*E*_{F}/*e*_{0} and charge *q* = *q*^{bulk} + *q*^{surf} can be calculated, allowing us to derive the equilibrium Schottky barrier height as a function of the applied voltage, as discussed in Sec. IV. By dividing the total charge of the electrode *q* by the surface area of the slab, we determine the charge density of each electrode surface.

### B. Computational details

Electronic-structure calculations are performed using the QUANTUM-ESPRESSO software.^{38} A slab of five layers, with three of these layers geometrically constrained on the semiconductor side, is found to be sufficient to converge the Fermi level of the semiconductor-solution system within 50 meV. We center the slab in the supercell with a separation of 14 Å between periodic slabs and optimize the geometry of the structure until interatomic forces are lower than 50 meV/Å. We use pseudopotentials generated with Perdew-Burke-Ernzerhof exchange correlation^{39} with the projector augmented wave (PAW) method from the SSSP library.^{40} The kinetic and charge density cutoffs are of 50 Ry and 750 Ry, respectively. The Brillouin zone is sampled with a shifted 5 × 5 × 1 Monkhorst-Pack grid and 0.03 Ry of Marzari-Vanderbilt smearing.^{41} We exploit the ENVIRON module for the continuum solvent.^{31}

## IV. RESULTS

To understand the influence of the applied potential on the silicon-water interface, we simulate seven representative surfaces terminated by different combinations of oxygen and hydrogen species, namely, O, (O, H), 2O, (2O, H), (4O, H), H, and 2H, as shown in Fig. 3(a), with a semiconductor carrier concentration of 10^{18} cm^{−3} for the semiclassical Mott-Schottky model. The calculated charge-voltage responses are shown in Fig. 3(b). Here, the potentials are measured with respect to the flatband potential *φ*_{FB}. These charge-voltage characteristics reveal that substantial differences arise from changing the surface termination due to charge pinning on the surface states associated with the dangling bonds that are present at the interface.

This observation is consistent with the fact that the charge of the depletion region causes a potential drop orders of magnitude larger than the surface dipole, due to poor electrostatic screening in the doped semiconductor. In contrast, surface terminations with significant charge accumulation (e.g., for the O and 2O covered surfaces) have most of their charge pinned at the interface. Figure 4 confirms this interpretation by showing the redistribution of the surface charge under applied voltage; for most of the surfaces, charge accumulates almost entirely on the adsorbate. For the 2H-terminated surface, however, the charge is distributed through all five layers of the simulated semiconductor, indicating that the charge extends deep into its depletion region.

We then interpret the electrical response of the semiconductor-solution interface analytically by developing a model for predicting the charge-voltage relationship. This analytic model is premised on the fact that surface states contribute a metallic capacitance, as evidenced by examining the interfacial density of states. This surface-state capacitive contribution is then connected to that of the ideal Mott-Schottky depletion layer, yielding

where the first term is the contribution from the bulk semiconductor to the total charge and the second term is that of the surface states. Here, $C$ is the capacitance of the surface states (which will be obtained by fitting the calculated charge-voltage response), *ϵ*_{s} is the dielectric constant of the semiconductor, and *φ* is the total potential drop across the electrode. We further introduce the space charge fraction, *θ* = *φ*_{s}/*φ*, to quantify the extent to which surface states dominate the electrical response of the photoelectrode. This ratio represents the fraction of potential on the electrode that falls within the bulk of the semiconductor; *θ* is not a constant and varies as a function of the total potential across the electrode.

An analytic expression for *θ* can be derived by noting that the amount of charge in the bulk semiconductor *q*^{bulk} is typically a constant fraction of the surface charge *q*^{surf} for all the terminations that we have studied (see the supplementary material). We thus introduce the coefficient *η* = *q*^{bulk}/*q*^{surf} that describes the fraction of the charge that resides in the semiconductor section, allowing us to derive the following expression for *θ*:

where $\varphi =(\u03f5s\u03f50e0N)/(\eta 2C2)$ is in units of volts and can be thought of as a switching potential, representing the point at which ∼25% of the total potential drop takes place across the bulk semiconductor. At potentials above the switching potential, the fraction shifts rapidly such that the potential drop in the depletion layer makes up the majority of the total electrostatic drop. The charge-voltage model for the semiconductor-solution interface can now be fully described analytically by fitting the two constants *η* and $C$. Using this fitting procedure, Eq. (10) can serve as an accurate analytic model of the charge-voltage response.

It should be noted that this model is only valid at potentials less than the band gap of the material; if the potential exceeds this band gap, Zener tunneling will cause the semiconductor to act as a metal. The fitted values of *η* and $C$ for each adsorbate are reported in Table I, showing that the surface charge is more delocalized into the depletion region for the 2H termination (*η* = 17.3%) than it is for the O termination (*η* = 4.4%), causing the capacitance $C$ to be higher in the former case.

. | Charge ratio η
. | Capacitance $C$ (μF/cm^{2})
. |
---|---|---|

O | 0.044 | 17.3 |

(O, H) | 0.060 | 14.0 |

2O | 0.040 | 17.6 |

(2O, H) | 0.063 | 13.1 |

(4O, H) | 0.067 | 8.1 |

H | 0.087 | 14.8 |

2H | 0.173 | 10.5 |

. | Charge ratio η
. | Capacitance $C$ (μF/cm^{2})
. |
---|---|---|

O | 0.044 | 17.3 |

(O, H) | 0.060 | 14.0 |

2O | 0.040 | 17.6 |

(2O, H) | 0.063 | 13.1 |

(4O, H) | 0.067 | 8.1 |

H | 0.087 | 14.8 |

2H | 0.173 | 10.5 |

With the charge-voltage relation in hand, we can calculate the surface free energy of each termination as a function of the applied potential using the Lippmann electrocapillary equation $\gamma =\gamma 0\u2212\u2009\u222b\phi FB\phi \sigma (\Phi )d\Phi $. Here, *σ* is the charge per surface area of the electrode, *γ* is the surface free energy of the charged slab at a certain potential *φ*, and *γ*_{0} is the surface free energy of the slab under neutral charge conditions at the flatband potential. We determine the surface free energy of a neutral surface following the computational standard hydrogen electrode (SHE) method^{42–44} by subtracting the energy of a surface terminated with adsorbates from the energy of the same surface without adsorbates and further subtracting the energy required to pull out a given adsorbate from the surrounding solution (see Sec. S2 of the supplementary material). We then calculate the free energy curves for each surface, as shown in Fig. 3(c). The structure with the lowest free energy at a given electrode potential is the thermodynamically stable configuration. Under most potential and pH conditions within the stability window of water, the (4O, H) configuration (which is the most oxidized termination tested) is the most stable, in agreement with the known tendency of silicon to oxidize in contact with water.^{8,9}

Next, we calculate the Schottky barrier for each termination by determining the potential drop across the bulk semiconductor within our simulation, determined as $\phi s=\phi \xaf0\u2212\phi \xaf(zc)$, where $\phi \xaf0$ is computed from Eq. (8) and $\phi \xaf(zc)$ is obtained from our calculations. As shown in Fig. 3(d), at low bias, the majority of the potential drop across the electrode is associated with the surface dipole $\delta \chi sm$ instead of the Schottky barrier, leading to a low space-charge fraction *θ*. The fraction of the potential drop across the semiconductor increases rapidly with applied bias and then shows a much slower increase after ∼0.3 V. The surface termination with the least surface charge pinning (2H) unsurprisingly exhibits the largest fraction of the potential drop inside the semiconductor in the calculated voltage range.

Having computed the Schottky barrier as a function of external bias, we can now calculate an important parameter, the charge-pinning fraction $S$, which describes the reduction of the Schottky barrier from its theoretical maximum due to charge pinning.^{45–50} Explicitly, the charge-pinning fraction is defined (for *n*-type semiconductors) as

where $S=1$ corresponds to an ideal semiconductor junction with no charge trapping, while a value of $S=0$ indicates that all the charge accumulates in surface states, inducing complete Fermi-level pinning. Here, *χ*_{m} and *χ*_{s} are the vacuum referenced electronegativity of the embedding medium and the semiconductor, respectively, which can be calculated as the opposite of the Fermi level when the vacuum section of the slab is aligned to zero.^{51–53} The value of $S$ is typically obtained experimentally by measuring the Schottky barrier of a semiconductor against different chemical environments and calculating the slope of the resulting line of best fit.

To calculate the equilibrium Schottky barrier *φ*_{s} and charge-pinning factor $S$, we first find the applied bias that sets the voltage of the electrodes to the hydrogen evolution potential. We then directly extract the Schottky barrier height at this applied bias by calculating the potential drop within the bulk semiconductor. By plotting the Schottky barrier height as a function of the applied bias to bring the interface into equilibrium with the hydrogen evolution redox level, we can find the charge-pinning fraction from the slope, as shown in Fig. 5. By linear regression, we measure a slope of $S\u22480.7$, which corresponds to a significant deviation from the ideal trend of $S=1$, reflecting the contribution from the surface dipole to the renormalization of the Schottky barrier height.

Our work highlights some of the contradictory requirements that limit the photocatalytic activity of silicon photoelectrodes. For the hydrogen evolution reaction, the ideal Schottky barrier would be positive—driving electrons from the bulk of the electrode to the surface. However, the two terminations with the most positive Schottky barriers (one oxygen and two oxygens adsorbed onto the surface) are both unstable in the redox window of water as shown in Fig. 3(c). The difficulty of simultaneously achieving surface stability and effective charge transfer across the silicon-water interface provides quantitative insights into the limited photocatalytic activity of silicon photoelectrodes in an aqueous environment. It is interesting to note that this interpretation requires only a single layer of oxide to form on the silicon surface and is not predicated on the formation of a thick deposit of silica. Instead, the low photocatalytic performance of silicon is here possibly explained by the surface states induced by an atomically thin oxide layer, suggesting a much more drastic influence on the solar-to-hydrogen conversion efficiency.

## V. CONCLUSION

We have examined the performance of silicon for water splitting as a function of the exposed termination of the semiconductor electrode from first principles. We have developed a methodology for predicting Schottky barriers accurately, taking into account the interactions between the charge pinned in the interfacial region and the charge accumulated within the depletion layer of the semiconductor. We have applied this methodology to predict the stability, Schottky barrier, and charge-pinning fraction of different surface terminations for silicon in water. Our study shows that the structures with the most favorable Schottky barriers for water splitting are electrochemically unstable, shedding light on the physical origins of their low solar-to-hydrogen conversion efficiency. Our work demonstrates the broad capabilities of the SCCS model, which can be used in conjunction with recent developments in predicting accurate electronic structures,^{54,55} to predict and understand the photocatalytic activity of passivated semiconductor electrodes from first principles.

## SUPPLEMENTARY MATERIAL

See supplementary material for a discussion of electrostatic potential references, further details on the model, and further explanations on the calculation of surface stability.

## ACKNOWLEDGMENTS

The authors acknowledge primary support from the National Science Foundation under Grant No. DMR-1654625 and partial support from the 3M Graduate Fellowship and the Penn State University Graduate Fellowship. The authors gratefully acknowledge Giulia Galli, Héctor D. Abruña, Roman Engel-Herbert, and Suzanne Mohney for fruitful discussions.

## REFERENCES

In this electronic-structure study, we describe the electronic properties of the junction in terms of differences between the Fermi level and the average electrostatic potential of the two terminals, instead of electronegativity differences as is common in the experimental literature. These conventions are equivalent under an appropriate shift of the relative energy scales, as shown in Sec. S1 of the supplementary material.

_{2}/metal contacts by inserting a BN monolayer