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.

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

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 (δχsm) = (Δs)(Δm), where (Δs) and (Δm) are the differences between the Fermi levels and average electrostatic potentials for the semiconductor and the surrounding medium,18 respectively.

FIG. 1.

(a) A pristine semiconductor in contact with a chemically inert medium (here, a solvent of wide redox stability window) has all of its potential drop located at the interface, leading to the formation of the surface dipole (δχsm). (b) When the Fermi level of an extrinsic semiconductor and that of a reactive medium (here, the same solvent with electron-accepting ionic species) are initially shifted from one another, the potential drop to realign the energy levels is distributed between the surface dipole δχsm, the Schottky barrier φs within the semiconductor, and the potential shift φm within the medium.

FIG. 1.

(a) A pristine semiconductor in contact with a chemically inert medium (here, a solvent of wide redox stability window) has all of its potential drop located at the interface, leading to the formation of the surface dipole (δχsm). (b) When the Fermi level of an extrinsic semiconductor and that of a reactive medium (here, the same solvent with electron-accepting ionic species) are initially shifted from one another, the potential drop to realign the energy levels is distributed between the surface dipole δχsm, the Schottky barrier φs within the semiconductor, and the potential shift φm within the medium.

Close modal

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 δχ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

φs=ΔsΔmδχsmφm,
(1)

where Δs and Δm are the differences between the Fermi level (divided by the electron charge e0) 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 1016 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 (δχsm) 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 barrier15,21–30

φsΔsΔm(δχsm).
(2)

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

φs=[ΔsΔm(δχsm)][δχsm(δχsm)].
(3)

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 δχsm(δχsm), 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.

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

ϵr=expζrsin2πζr/2πlnϵm,
(4)

where ϵm is the dielectric constant of the electrolytic medium and ζ(r) = (ln ρmax − ln ρ(r))/(ln ρ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 Gcav = γS and the dispersion and repulsion effects by Gdis+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 = ∫dr(dΘ/)|∇ρ| and V = ∫drΘ(ρ), 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 species34 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 of the neutral surface through

φFB=EFe0,
(5)

where e0 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 qsurf and the bulk semiconductor qbulk such that q = qbulk + qsurf. (Here, the charge qbulk 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 = zc 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),

n(r)=Nexp(φ0φ(r))/(kBT),
(6)

where N is the concentration of electron-donating defects, kB 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 

φz2=2Nϵ0φφ0+kBTexpφ0φkBT1.
(7)

Then, noting that the exponential term is typically small (i.e., φφ0kBT), the asymptotic value φ¯0 of the difference between the electrostatic potential of the charged and neutral slabs φ¯ can be obtained from

φ¯0=φ¯(zc)kBTϵ02Ndφ¯dz(zc)2.
(8)

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 φ¯0 [Eq. (8)] to the Fermi level of the neutral slab,

EFbulk=e0φ¯0+EF.
(9)

Once the equilibrium Fermi level is known, the voltage φ = −EF/e0 and charge q = qbulk + qsurf 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.

FIG. 2.

The equilibrium state of the silicon-water interface is found by computing the Fermi level of the bulk semiconductor EFbulk and that of the surface region EFsurf as a function of the surface charge. At the intersections of these two electronic levels, the chemical potential of the electrons is constant across the entire system, as shown by the horizontal dashed line on the electrostatic profile of the inset. Repeating this procedure for multiple net charges yields the charge-voltage response of the silicon electrode.

FIG. 2.

The equilibrium state of the silicon-water interface is found by computing the Fermi level of the bulk semiconductor EFbulk and that of the surface region EFsurf as a function of the surface charge. At the intersections of these two electronic levels, the chemical potential of the electrons is constant across the entire system, as shown by the horizontal dashed line on the electrostatic profile of the inset. Repeating this procedure for multiple net charges yields the charge-voltage response of the silicon electrode.

Close modal

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 correlation39 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 

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

FIG. 3.

(a) Surface terminations of the silicon (110) electrodes, showing the configuration of the hydrogen and oxygen adsorbates. (b) Charge-voltage response of Si (110) structures with these adsorbates. Data points closely follow the proposed analytic model. (c) Surface free energy γ of each silicon termination at pH = 7. The termination with the lowest free energy is the most stable at that potential. (d) Space-charge fraction θ of each electrode, the fraction of applied bias that develop in the depletion region.

FIG. 3.

(a) Surface terminations of the silicon (110) electrodes, showing the configuration of the hydrogen and oxygen adsorbates. (b) Charge-voltage response of Si (110) structures with these adsorbates. Data points closely follow the proposed analytic model. (c) Surface free energy γ of each silicon termination at pH = 7. The termination with the lowest free energy is the most stable at that potential. (d) Space-charge fraction θ of each electrode, the fraction of applied bias that develop in the depletion region.

Close modal

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.

FIG. 4.

Accumulation of electronic charge in the explicit quantum mechanical portion of the electrode under bias. The isocontours indicate the positions of 90% of the charge. When the majority of charge is at the surface, a significant surface dipole is formed, corresponding to a lower Schottky barrier. When charge is distributed across several layers (e.g., 2H), a smaller surface dipole is formed and a larger potential drop takes place across the space-charge region, inducing a higher Schottky barrier.

FIG. 4.

Accumulation of electronic charge in the explicit quantum mechanical portion of the electrode under bias. The isocontours indicate the positions of 90% of the charge. When the majority of charge is at the surface, a significant surface dipole is formed, corresponding to a lower Schottky barrier. When charge is distributed across several layers (e.g., 2H), a smaller surface dipole is formed and a larger potential drop takes place across the space-charge region, inducing a higher Schottky barrier.

Close modal

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

Q=2ϵ0ϵse0Nθφ+C(1θ)φ,
(10)

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 qbulk is typically a constant fraction of the surface charge qsurf for all the terminations that we have studied (see the supplementary material). We thus introduce the coefficient η = qbulk/qsurf that describes the fraction of the charge that resides in the semiconductor section, allowing us to derive the following expression for θ:

θ=1+ϕ/φϕ2/φ2+2ϕ/φ,
(11)

where ϕ=(ϵsϵ0e0N)/(η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.

TABLE I.

Fitted surface state properties for the seven surface configurations tested. The fraction of the charge on the surface states that is on the bulk of the semiconductor is represented by η = qbulk/qsurf. The capacitance of the surface state, assuming a metal like distribution, is denoted as C.

Charge ratio ηCapacitance C (μF/cm2)
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 
0.087 14.8 
2H 0.173 10.5 
Charge ratio ηCapacitance C (μF/cm2)
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 
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 γ=γ0φFBφσ(Φ)dΦ. 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) method42–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 φs=φ¯0φ¯(zc), where φ¯0 is computed from Eq. (8) and φ¯(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 δχ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

φs=S(χmχs),
(12)

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

FIG. 5.

The equilibrium barrier height as a function of the difference between the hydrogen evolution potential and the flatband potential. An ideal semiconductor junction would show a unit slope, S=1 (dashed line). Due to charge trapping by surface states, the Schottky barrier heights are significantly reduced, corresponding to a charge-pinning fraction equal to S=0.7±0.005 for the silicon-water system.

FIG. 5.

The equilibrium barrier height as a function of the difference between the hydrogen evolution potential and the flatband potential. An ideal semiconductor junction would show a unit slope, S=1 (dashed line). Due to charge trapping by surface states, the Schottky barrier heights are significantly reduced, corresponding to a charge-pinning fraction equal to S=0.7±0.005 for the silicon-water system.

Close modal

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.

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.

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

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.

1.
Report of the Hydrogen Production Expert Panel: A Subcommittee of the Hydrogen & Fuel Cell Technical Advisory Committee, Technical Report, United States Department of Energy, 2013.
2.
N. S.
Lewis
, “
Powering the planet
,”
MRS Bull.
32
,
808
820
(
2007
).
3.
A.
Fujishima
and
K.
Honda
, “
Electrochemical photolysis of water at a semiconductor electrode
,”
Nature
238
,
37
38
(
1972
); e-print arXiv:1011.1669v3.
4.
T.
Jafari
,
E.
Moharreri
,
A.
Amin
,
R.
Miao
,
W.
Song
, and
S.
Suib
, “
Photocatalytic water splitting—The untamed dream: A review of recent advances
,”
Molecules
21
,
900
(
2016
); e-print arXiv:0307014 [cond-mat].
5.
P.
Leempoel
,
M.
Castro-Acuna
,
F.-R. F.
Fan
, and
A. J.
Bard
, “
The effect of light intensity and iodine doping on the stabilization of n-silicon by phthalocyanine films
,”
J. Phys. Chem.
86
,
1396
1400
(
1982
).
6.
Y. W.
Chen
,
J. D.
Prange
,
S.
Duhnen
,
Y.
Park
,
M.
Gunji
,
C. E. D.
Chidsey
, and
P. M.
McIntyre
, “
Atomic layer-deposited tunnel oxide stabilizes silicon photoanodes for water oxidation
,”
Nat. Mater.
10
,
539
544
(
2011
).
7.
J. M.
Bolts
,
A. B.
Bocarsly
,
M. C.
Palazzotto
,
E. G.
Walton
,
N. S.
Lewis
, and
M. S.
Wrighton
, “
Chemically derivatized n-type silicon photoelectrodes. Stabilization to surface corrosion in aqueous electrolyte solutions and mediation of oxidation reactions by surface-attached electroactive ferrocene reagents
,”
J. Am. Chem. Soc.
101
,
1378
1385
(
1979
).
8.
W. R.
Runyan
, “
Silicon
,” in
Kirk-Othmer Encyclopedia of Chemical Technology
(
John Wiley & Sons, Inc.
,
2013
), pp.
1
22
.
9.
D.
Graf
,
M.
Grundner
, and
R.
Schulz
, “
Reaction of water with hydrofluoric acid treated silicon(111) and (100) surfaces
,”
J. Vac. Sci. Technol., A
7
,
808
(
1989
).
10.
J.
Basquert
, “
Planar and nanostructured semiconductor junctions
,” in
Nanostructured Energy Devices: Equilibrium Concepts and Kinetics
(
CRC Press
,
2014
), Chap. 9, pp.
275
320
.
11.
K.
Rajeshwar
, “
Fundamentals of semiconductors electrochemistry and photoelectrochemistry
,” in
Encyclopedia of Electrochemistry
(
Wiley
,
2007
), Chap. 1, pp.
1
51
.
12.
L. J.
Brillson
,
An Essential Guide to Electronic Material
(
John Wiley & Sons, Inc.
,
Hoboken
,
2016
).
13.
N.
Kharche
,
J. T.
Muckerman
, and
M. S.
Hybertsen
, “
First-principles approach to calculating energy level alignment at aqueous semiconductor interfaces
,”
Phys. Rev. Lett.
113
,
176802
(
2014
).
14.
J.
Cheng
and
M.
Sprik
, “
Alignment of electronic energy levels at electrochemical interfaces
,”
Phys. Chem. Chem. Phys.
14
,
11245
11267
(
2012
).
15.
Y.
Ping
,
W. A.
Goddard
, and
G. A.
Galli
, “
Energetics and solvation effects at the photoanode/catalyst interface: Ohmic contact versus Schottky barrier
,”
J. Am. Chem. Soc.
137
,
5264
5267
(
2015
).
16.
Y.
Wu
,
M. K. Y.
Chan
, and
G.
Ceder
, “
Prediction of semiconductor band edge positions in aqueous environments from first principles
,”
Phys. Rev. B
83
,
235301
(
2011
).
17.
Q.
Campbell
and
I.
Dabo
, “
Quantum-continuum calculation of the surface states and electrical response of silicon in solution
,”
Phys. Rev. B
95
,
205308
(
2017
).
18.

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.

19.
J. O.
Bockris
and
A.
Reddy
,
Modern Electrochemistry 2B
, 2nd ed. (
Kluwer Academic Publishers
,
2000
).
20.
W.
Schmickler
and
E.
Santos
,
Interfacial Electrochemistry
, 2nd ed. (
Springer
,
2010
).
21.
Y.
Jiao
,
A.
Hellman
,
Y.
Fang
,
S.
Gao
, and
M.
Käll
, “
Schottky barrier formation and band bending revealed by first- principles calculations
,”
Sci. Rep.
5
,
11374
(
2015
).
22.
M.
Farmanbar
and
G.
Brocks
, “
Controlling the Schottky barrier at MoS2/metal contacts by inserting a BN monolayer
,”
Phys. Rev. B
91
,
161304(R)
(
2015
).
23.
J.
Goniakowski
and
C.
Noguera
, “
Electronic states and Schottky barrier height at metal/MgO(100) interfaces
,”
Interface Sci.
12
,
93
103
(
2004
).
24.
M.
Stengel
,
P.
Aguado-Puente
,
N. A.
Spaldin
, and
J.
Junquera
, “
Band alignment at metal/ferroelectric interfaces: Insights and artifacts from first principles
,”
Phys. Rev. B
83
,
235112
(
2011
).
25.
P. A.
Khomyakov
,
G.
Giovannetti
,
P. C.
Rusu
,
G.
Brocks
,
J.
Van Den Brink
, and
P. J.
Kelly
, “
First-principles study of the interaction and charge transfer between graphene and metals
,”
Phys. Rev. B
79
,
195425
(
2009
).
26.
Y.
Dong
and
L. J.
Brillson
, “
First-principles studies of metal (111)/ZnO{0001} interfaces
,”
J. Electron. Mater.
37
,
743
748
(
2008
).
27.
R. A.
McKee
,
F. J.
Walker
,
M. B.
Nardelli
,
W. A.
Shelton
, and
G. M.
Stocks
, “
The interface phase and the Schottky barrier for a crystalline dielectric on silicon
,”
Science
300
,
1726
1730
(
2003
).
28.
J. E.
Padilha
,
A.
Fazzio
, and
A. J.
Da Silva
, “
van der Waals heterostructure of phosphorene and graphene: Tuning the Schottky barrier and doping by electrostatic gating
,”
Phys. Rev. Lett.
114
,
066803
(
2015
).
29.
S.
Picozzi
,
A.
Continenza
, and
S.
Massidda
, “
Structural and electronic properties of ideal nitride/Al interfaces
,”
Phys. Rev. B
57
,
4849
4856
(
1998
).
30.
S. H.
Jeon
,
B. H.
Park
,
J.
Lee
,
B.
Lee
, and
S.
Han
, “
First-principles modeling of resistance switching in perovskite oxide material
,”
Appl. Phys. Lett.
89
,
042904
(
2006
).
31.
O.
Andreussi
,
I.
Dabo
, and
N.
Marzari
, “
Revised self-consistent continuum solvation in electronic-structure calculations
,”
J. Chem. Phys.
136
,
064102
(
2012
).
32.
G.
Fisicaro
,
L.
Genovese
,
O.
Andreussi
,
S.
Mandal
,
N. N.
Nair
,
N.
Marzari
, and
S.
Goedecker
, “
Soft-sphere continuum solvation in electronic-structure calculations
,”
J. Chem. Theory Comput.
13
,
3829
3845
(
2017
).
33.
O.
Andreussi
and
G.
Fisicaro
, “
Continuum embeddings in condensed-matter simulations
,”
Int. J. Quantum Chem.
119
,
e25725
(
2018
).
34.
C.
Dupont
,
O.
Andreussi
, and
N.
Marzari
, “
Self-consistent continuum solvation (SCCS): The case of charged systems
,”
J. Chem. Phys.
139
,
214110
(
2013
).
35.
N.
Hörmann
,
O.
Andreussi
, and
N.
Marzari
, “
Grand canonical simulations of electrochemical interfaces in implicit solvation models
,”
J. Chem. Phys.
150
,
041730
(
2019
).
36.
F.
Nattino
,
M.
Truscott
,
N.
Marzari
, and
O.
Andreussi
, “
Continuum models of the electrochemical diffuse layer in electronic-structure calculations
,”
J. Chem. Phys.
150
,
041722
(
2019
); e-print arXiv:1810.09797.
37.
K.
Gelderman
,
L.
Lee
, and
S. W.
Donne
, “
Flat-Band potential of a semiconductor: Using the Mott–Schottky equation
,”
J. Chem. Educ.
84
,
685
(
2007
).
38.
P.
Giannozzi
,
S.
Baroni
,
N.
Bonini
,
M.
Calandra
,
R.
Car
,
C.
Cavazzoni
,
D.
Ceresoli
,
G. L.
Chiarotti
,
M.
Cococcioni
,
I.
Dabo
,
A.
Dal Corso
,
S.
de Gironcoli
,
S.
Fabris
,
G.
Fratesi
,
R.
Gebauer
,
U.
Gerstmann
,
C.
Gougoussis
,
A.
Kokalj
,
M.
Lazzeri
,
L.
Martin-Samos
,
N.
Marzari
,
F.
Mauri
,
R.
Mazzarello
,
S.
Paolini
,
A.
Pasquarello
,
L.
Paulatto
,
C.
Sbraccia
,
S.
Scandolo
,
G.
Sclauzero
,
A. P.
Seitsonen
,
A.
Smogunov
,
P.
Umari
, and
R. M.
Wentzcovitch
, “
QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials
,”
J. Phys.: Condens. Matter
21
,
395502
(
2009
).
39.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
,
3865
3868
(
1996
).
40.
K.
Lejaeghere
,
G.
Bihlmayer
,
T.
Björkman
,
P.
Blaha
,
S.
Blügel
,
V.
Blum
,
D.
Caliste
,
I. E.
Castelli
,
S. J.
Clark
,
A.
Dal Corso
,
S.
De Gironcoli
,
T.
Deutsch
,
J. K.
Dewhurst
,
I.
Di Marco
,
C.
Draxl
,
M.
Dułak
,
O.
Eriksson
,
J. A.
Flores-Livas
,
K. F.
Garrity
,
L.
Genovese
,
P.
Giannozzi
,
M.
Giantomassi
,
S.
Goedecker
,
X.
Gonze
,
O.
Grånäs
,
E. K.
Gross
,
A.
Gulans
,
F.
Gygi
,
D. R.
Hamann
,
P. J.
Hasnip
,
N. A.
Holzwarth
,
D.
Iuşan
,
D. B.
Jochym
,
F.
Jollet
,
D.
Jones
,
G.
Kresse
,
K.
Koepernik
,
E.
Küçükbenli
,
Y. O.
Kvashnin
,
I. L.
Locht
,
S.
Lubeck
,
M.
Marsman
,
N.
Marzari
,
U.
Nitzsche
,
L.
Nordström
,
T.
Ozaki
,
L.
Paulatto
,
C. J.
Pickard
,
W.
Poelmans
,
M. I.
Probert
,
K.
Refson
,
M.
Richter
,
G. M.
Rignanese
,
S.
Saha
,
M.
Scheffler
,
M.
Schlipf
,
K.
Schwarz
,
S.
Sharma
,
F.
Tavazza
,
P.
Thunström
,
A.
Tkatchenko
,
M.
Torrent
,
D.
Vanderbilt
,
M. J.
Van Setten
,
V.
Van Speybroeck
,
J. M.
Wills
,
J. R.
Yates
,
G. X.
Zhang
, and
S.
Cottenier
, “
Reproducibility in density functional theory calculations of solids
,”
Science
351
,
aad3000
(
2016
).
41.
N.
Marzari
,
D.
Vanderbilt
, and
M. C.
Payne
, “
Ensemble density-functional theory for ab-initio molecular dynamics of metals and finite-temperature insulators
,”
Phys. Rev. Lett.
79
,
1337
1340
(
1997
).
42.
I. C.
Man
,
H.-Y.
Su
,
F.
Calle-Vallejo
,
H. A.
Hansen
,
J. I.
Martínez
,
N. G.
Inoglu
,
J.
Kitchin
,
T. F.
Jaramillo
,
J. K.
Nørskov
, and
J.
Rossmeisl
, “
Universality in oxygen evolution electrocatalysis on oxide surfaces
,”
ChemCatChem
3
,
1085
(
2011
).
43.
J. K.
Norskov
,
J.
Rossmeisl
,
A.
Logadottir
,
L.
Lindqvist
,
J. R.
Kitchin
,
T.
Bligaard
, and
H.
Jonsson
, “
Origin of the overpotential for oxygen reduction at a fuel-cell cathode
,”
J. Phys. Chem. B
108
,
17886
17892
(
2004
).
44.
J.
Rossmeisl
,
Z. W.
Qu
,
H.
Zhu
,
G. J.
Kroes
, and
J. K.
Nørskov
, “
Electrolysis of water on oxide surfaces
,”
J. Electroanal. Chem.
607
,
83
89
(
2007
).
45.
S.
Kurtin
,
T. C.
McGill
, and
C. A.
Mead
, “
Fundamental transition in the electronic nature of solids
,”
Phys. Rev. Lett.
22
,
1433
1436
(
1969
).
46.
L.
Brillson
, “
Transition in Schottky barrier formation with chemical reactivity
,”
Phys. Rev. Lett.
40
,
260
263
(
1978
).
47.
R. T.
Tung
, “
The physics and chemistry of the Schottky barrier height
,”
Appl. Phys. Rev.
1
,
11304
(
2014
).
48.
A.
Natarajan
,
G.
Oskam
, and
P. C.
Searson
, “
The potential distribution at the semiconductor/solution interface
,”
J. Phys. Chem. B
102
,
7793
7799
(
1998
).
49.
A. J.
Bard
,
A. B.
Bocarsly
,
F. R. F.
Fan
,
E. G.
Walton
, and
M. S.
Wrighton
, “
The concept of Fermi level pinning at semiconductor/liquid junctions. Consequences for energy conversion efficiency and selection of useful solution redox couples in solar devices
,”
J. Am. Chem. Soc.
102
,
3671
3677
(
1980
).
50.
A. M.
Cowly
and
S. M.
Sze
, “
Surface states and barrier height of metal-semiconductor systems
,”
J. Appl. Phys.
36
,
3212
(
1965
).
51.
S. E.
Weitzner
and
I.
Dabo
, “
Quantum-continuum simulation of underpotential deposition at electrified metal-solution interfaces
,”
npj Comput. Mater.
3
,
1
(
2017
).
52.
N.
Keilbart
,
Y.
Okada
,
A.
Feehan
,
S.
Higai
, and
I.
Dabo
, “
Quantum-continuum simulation of the electrochemical response of pseudocapacitor electrodes under realistic conditions
,”
Phys. Rev. B
95
,
115423
(
2017
).
53.
S. E.
Weitzner
and
I.
Dabo
, “
Voltage effects on the stability of Pd ensembles in PdAu/Au(111) surface alloys
,”
J. Chem. Phys.
150
,
041715
(
2019
).
54.
N. L.
Nguyen
,
N.
Colonna
,
A.
Ferretti
, and
N.
Marzari
, “
Koopmans-compliant spectral functionals for extended systems
,”
Phys. Rev. X
8
,
021051
(
2018
).
55.
I.
Timrov
,
N.
Marzari
, and
M.
Cococcioni
, “
Hubbard parameters from density-functional perturbation theory
,”
Phys. Rev. B
98
,
085127
(
2018
).

Supplementary Material