Quantum-continuum calculation of the surface states and electrical response of silicon in solution

TL;DR

This paper introduces a quantum-continuum computational method to analyze surface states and electrical responses of silicon electrodes in solution, offering microscopic insights into electrochemical interfaces relevant for practical applications.

Contribution

The study develops an embedded density-functional theory approach with SCCS to model semiconductor-electrolyte interfaces, integrating quantum surface descriptions with continuum band bending.

Findings

Reveals surface states influence silicon electrification at low voltages

Provides microscopic understanding of semiconductor-electrolyte interfaces

Demonstrates the method's applicability to electrochemical reactions

Abstract

A wide range of electrochemical reactions of practical importance occur at the interface between a semiconductor and an electrolyte. We present an embedded density-functional theory method using the recently released self-consistent continuum solvation (SCCS) approach to study these interfaces. In this model, a quantum description of the surface is incorporated into a continuum representation of the bending of the bands within the electrode. The model is applied to understand the electrical response of silicon electrodes in solution, providing microscopic insights into the low-voltage region, where surface states determine the electrification of the semiconductor electrode.

Tables

Table 1. Table 1: Fitted surface state properties for the 10 surface configurations tested. The percentage of the total potential drop that takes place within Region I, the bulk of the semiconductor, is represented by γ I subscript 𝛾 I \gamma_{\rm I} (%). The capacitance of the surface state, assuming a metal like distribution, is represented as 𝒞 II subscript 𝒞 II {\mathscr{C}}_{\rm II} .

	${\gamma }_I$ (%)	${𝒞}_{\mathrm{II}}$ ($\mu$F cm-2)
Si+O*	3.75	78.9
Si+O*+H*	5.39	54.1
Si+2O*	3.38	75.2
Si+2O*+H	5.92	38.9
Si+4O*+H*	11.7	25.3
Rutile (110)	0.03	1.28
Cristobalite (100)	1.69	786
Cristobalite (110)	8.23	1.89
Quartz (1000)	1.95	41.3
Quartz (11$\overline{2}$0)	4.71	30.1

Equations

$${\Phi }_{\mathrm{FB}}=({\phi }_0-{ϵ}_{\mathrm{F}})/e.$$

$$F=F_{\mathrm{I}}+F_{\mathrm{II}}+F_{\mathrm{III}}-\frac{1}{2}\int d\mathbf{r}{ϵ}_0ϵ(\mathbf{r}){\left|\nabla \phi (\mathbf{r})\right|}^2,$$

$$F_{\mathrm{I}}$$

$$s(f,f_{\mathrm{d}})=-k_{\mathrm{B}}\left[\left(f_{\mathrm{d}}-f\right)\ln \left(1-\frac{f}{f_{\mathrm{d}}}\right)+f\ln \left(\frac{f}{f_{\mathrm{d}}}\right)\right].$$

$$n_{\mathrm{d}}(\mathbf{r})=\frac{\mathcal{N}}{2}\left[\mathrm{erfc}\left(\frac{z-z_{\mathrm{I}}}{{\sigma }_{\mathrm{I}}}\right)+1\right]$$

$$p_{\mathrm{d}}(\mathbf{r})=\frac{\mathcal{P}}{2}\left[\mathrm{erfc}\left(\frac{z-z_{\mathrm{I}}}{{\sigma }_{\mathrm{I}}}\right)+1\right]$$

$$F_{\mathrm{III}}$$

$$\sigma (c_{\pm },c^{\circ })=-k_{\mathrm{B}}\left[f\ln \left(\frac{c_{\pm }}{c^{\circ }}\right)-c_{\pm }\right]$$

$$c^{\circ }(\mathbf{r})=\frac{\mathcal{C}}{2}\left[\mathrm{erfc}\left(\frac{z_{\mathrm{III}}-z}{{\sigma }_{\mathrm{III}}}\right)+1\right],$$

$$F_{\mathrm{II}}$$

$$\nabla \left({ϵ}_0ϵ(\mathbf{r})\nabla \phi (\mathbf{r})\right)$$

$$n(\mathbf{r})=n_{\mathrm{d}}(\mathbf{r}){\left[1+\exp \left(\frac{\phi (\mathbf{r})+{ϵ}_{\mathrm{C}}-{ϵ}_{\mathrm{F}}}{k_{\mathrm{B}}T}\right)\right]}^{-1}$$

$$p(\mathbf{r})=p_{\mathrm{d}}(\mathbf{r}){\left[1+\exp \left(\frac{{ϵ}_{\mathrm{F}}-{ϵ}_{\mathrm{V}}-\phi (\mathbf{r})}{k_{\mathrm{B}}T}\right)\right]}^{-1}$$

$$c_{\pm }(\mathbf{r})=c^{\circ }(\mathbf{r})\exp \left(\mp \frac{\phi (\mathbf{r})}{k_{\mathrm{B}}T}\right)$$

$$p(\mathbf{r})=p_{\mathrm{d}}(\mathbf{r})=0$$

$$n(\mathbf{r})=n_{\mathrm{d}}(\mathbf{r})\exp \left(\frac{{\phi }_0-\phi (\mathbf{r})}{k_{\mathrm{B}}T}\right),$$

$$n(\mathbf{r})=n_{\mathrm{d}}(\mathbf{r})=0$$

$$p(\mathbf{r})=p_{\mathrm{d}}(\mathbf{r})\exp \left(\frac{\phi (\mathbf{r})-{\phi }_0}{k_{\mathrm{B}}T}\right).$$

$$\frac{d^2\phi }{d{z}^2}=\frac{\mathcal{N}}{{ϵ}_0{ϵ}_{\mathrm{I}}}\left[1-\exp \left(\frac{{\phi }_0-\phi }{k_{\mathrm{B}}T}\right)\right]$$

$$\frac{d^2\phi }{d{z}^2}=\frac{\mathcal{P}}{{ϵ}_0{ϵ}_{\mathrm{I}}}\left[\exp \left(\frac{\phi -{\phi }_0}{k_{\mathrm{B}}T}\right)-1\right]$$

$${\left(\frac{d\phi }{dz}\right)}^2=\frac{2\mathcal{N}}{{ϵ}_0{ϵ}_{\mathrm{I}}}\left[\phi -{\phi }_0+k_{\mathrm{B}}T\left(e^{\frac{{\phi }_0-\phi }{k_{\mathrm{B}}T}}-1\right)\right]$$

$${\left(\frac{d\phi }{dz}\right)}^2=\frac{2\mathcal{P}}{{ϵ}_0{ϵ}_{\mathrm{I}}}\left[{\phi }_0-\phi +k_{\mathrm{B}}T\left(e^{\frac{\phi -{\phi }_0}{k_{\mathrm{B}}T}}-1\right)\right].$$

$$q_{\mathrm{I}}+q_{\mathrm{II}}+q_{\mathrm{III}}=0.$$

$${\epsilon }_{\mathrm{F},\mathrm{I}}={\bar{\phi }}_0-e{\Phi }_{\mathrm{FB}}$$

$${\bar{\phi }}_0=\bar{\phi }(z_{\mathrm{c}})-k_{\mathrm{B}}T-\frac{{ϵ}_0{ϵ}_{\mathrm{I}}}{2\mathcal{N}}{\left(\frac{d\bar{\phi }}{dz}(z_{\mathrm{c}})\right)}^2$$

$${\bar{\phi }}_0=\bar{\phi }(z_{\mathrm{c}})+k_{\mathrm{B}}T+\frac{{ϵ}_0{ϵ}_{\mathrm{I}}}{2\mathcal{P}}{\left(\frac{d\bar{\phi }}{dz}(z_{\mathrm{c}})\right)}^2$$

$${\epsilon }_{\mathrm{F},\mathrm{I}}={\epsilon }_{\mathrm{F},\mathrm{II}}.$$

Full text

I Introduction

Predicting the electrical response and stability of semiconductor–solution interfaces is of central relevance to a wide array of electrochemical and photoelectrochemical systems. These interfaces are involved in the photocatatytic splitting of water,Kudo and Miseki (2009) the photoreduction of carbon dioxide into hydrocarbons,Roy et al. (2010) the electrochemical etching of semiconductor surfaces,Tenne and Hodes (1983); Jin et al. (2004); Peng et al. (2006) the storage of energy at metal oxide electrodes,Xia et al. (2012) and the use of quantum dots as biological markers.Gill et al. (2008) The pivotal role of semiconductor electrodes at the frontier between solid state physics and electrochemistry provides a compelling motivation to study their behavior in solution.

Density-functional theory has been used to search for new photocatalysts,Bhatt and Lee (2015); Castelli et al. (2012, 2015) assess the alignment of the valence and conduction bands of semiconductor electrodes with the redox potentials of species in solution,Ping et al. (2015a, b); Cheng and Selloni (2010); Wu et al. (2011); Yang et al. (2012) determine reaction pathways for photoelectrochemical reactions,Setvin et al. (2016); Li and Selloni (2016); Cheng et al. (2016); Shen et al. (2010) and elucidate the dynamical interactions of the solvent molecules with the surface of the semiconductor.Oshikiri and Boero (2006); Rivest et al. (2014); Lee and Selloni (2014); Tritsaris et al. (2014); Kharche et al. (2014); Cheng and Sprik (2014) The same calculations can be applied to predict the electrical response of semiconductor–solution interfaces as long as they account for the long-range decay of the electrostatic potential within the semiconductor depletion region. However, electrostatic screening in doped semiconductors is much less effective than in metals,Rajeshwar (2007) causing the interfacial electric field to penetrate up to 10–103 nm into the electrode for typical dopant concentrations of $10^{16}$–$10^{18}$ cm*-3*. These length scales render the first-principles simulation of the interface computationally demanding.

Therefore, it is necessary to develop efficient models that will capture the essential features of a semiconductor–electrolyte interface at reduced computational cost. To this end, we exploit and further develop the self–consistent continuum solvation (SCCS) approach proposed by Andreussi et al.Andreussi et al. (2012) to simulate semiconductor electrodes under applied voltage in electrolytic media. While this method has been successful in modeling metal electrodes,Bonnet and Marzari (2013); Sementa et al. (2016); Montemore et al. (2016); Keilbart et al. (2017); Weitzner and Dabo (2017) no previous study has focused on describing band bending at semiconductor electrodes using a self–consistent continuum solvation approach.

To illustrate band bending, the electrostatic profile of a semiconductor–electrolyte interface is shown in Fig. 1. In order to reach equilibrium, the electrochemical potential of the semiconductor and the solution need to be equal. Since the excess charge is accommodated within the semiconductor by a low concentration of dopants, the electrostatic potential is seen to decay gradually across the extended depletion region (the Mott–Schottky layer). This ideal picture is made more complicated, however, by the presence of surface states, which result from the adsorption of ionic species. By trapping charge at the interface, these surface states lead to much stronger electrostatic screening within the semiconductor. It is the goal of this work to develop a quantum–continuum model for simulating the response of semiconductor–electrolyte interfaces under electrification, including the influence of the surface states. We apply the model to prototypical silicon electrodes in solution in an effort to elucidate the connection between their electrical response and surface structure at the molecular level.

II Method

II.1 Interface energy

The first step in constructing the model is to partition the system into the three regions shown in Fig. 2. Region I represents the bulk of the semiconductor, which will be modeled at the continuum level. Region II corresponds to the surface of the electrode; this region will be described quantum mechanically to represent the adsorbed species and resulting surface states. Region III denotes the electrolytic solution, which consists of a diffuse distribution of ions in a polarizable continuum. Also illustrated is the flatband potential $\Phi_{\rm FB}$, which corresponds to the difference between the asymptotic value $\varphi_{0}$ of the potential inside Region I and the Fermi energy of the neutral slab:

[TABLE]

Having defined the three regions, the free energy of the system is written as

[TABLE]

where $\epsilon$ denotes the space-dependent dielectric permittivity across the interface, and $F_{\rm I}$, $F_{\rm II}$, and $F_{\rm III}$ stand for the free energies of the continuum space charge, quantum surface slab, and continuum electrolyte, respectively.

This dielectric permittivity is defined using the self-consistent continuum solvation model (SCCS).Andreussi et al. (2012) In this model, a dielectric cavity is created at each lateral facet. On the semiconductor side, the local dielectric permittivity can be written as $\epsilon(\bm{r})=\exp[(\zeta(\bm{r})-\sin(2\pi\zeta(\bm{r}))/2\pi)\ln\epsilon_{\rm I}]$ where $\epsilon_{\rm I}$ is the dielectric constant of the bulk of the semiconductor and $\zeta(\bm{r})=(\ln\rho_{\rm max}-\ln\rho(\bm{r}))/(\ln\rho_{\rm max}-\ln\rho_{\rm min})$ is used as a smooth switching function, marking the transition between the quantum and continuum regions. Here, $\rho_{\rm min}$ and $\rho_{\rm max}$ serve as the density thresholds specifying the inner and outer isocontours of the dielectric cavity. We employ the same parametrization in Region III, replacing $\epsilon_{\rm I}$ with $\epsilon_{\rm III}$, the dielectric constant of the electrolyte. We specifically use $\rho_{\rm max}=$5\text{\times}{10}^{-3}$$ a.u. and $\rho_{\rm min}=10^{-4}$ a.u. for our calculations.

Focusing first on Region I, the contribution to the free energy can be expressed in terms of the local density of negative charge carriers $n$ and positive charge carriers $p$ as the sum of electrostatic and entropic terms:

[TABLE]

with

[TABLE]

In Eq. (3), $\epsilon_{\rm F}$ denotes the Fermi energy, $\epsilon_{\rm V}$ is the electronic energy at the top of the valence band, and $\epsilon_{\rm C}$ is the energy at the bottom of the conduction band. Furthermore, it is assumed that the donor and acceptor levels are shallow so that their energy levels sit at the band edges. Moreover, $s$ stands for the Fermi–Dirac entropy; it depends locally on the smooth switching functions

[TABLE]

with ${\mathscr{N}}$ and ${\mathscr{P}}$ being the concentrations of electron-donating and electron-accepting defects, and $z_{\rm I}$ and $\sigma_{\rm I}$ being the location and spatial extent of the transition between the semiconductor and surface.

Likewise, the free energy of the electrolyte can be expressed in terms of the concentrations $c_{+}$ and $c_{-}$ of the positive and negative ions as

[TABLE]

with

[TABLE]

In Eq. (5), we take a symmetric 1:1 ionic solution, in which the maximal ion concentration is defined as

[TABLE]

where $\mathscr{C}$ is the equilibrium ionic concentration inside the electrolyte, and $z_{\rm III}$ and $\sigma_{\rm III}$ are the location and spread of the frontier between the surface and electrolyte.

Finally, the free energy $F_{\rm II}$ is expressed as a Kohn–Sham functional of the density of the electrons $\rho_{-}$ and distribution of the atomic cores $\rho_{+}$:

[TABLE]

where $T_{\rm s}$ is the kinetic energy of the auxiliary system within the independent-electron mapping and $E_{\rm Hxc}$ is the sum of the Hartree and exchange-correlation energies. The electronic temperature and the entropy of the electronic smearing are denoted $\theta$ and ${\mathscr{S}}$, respectively.

II.2 Interface electrostatics

With the expression of the free energy in hand, the equilibrium charge density can be obtained by variations with respect to the occupations of the doping levels, ionic concentrations, and electrostatic potential, yielding the following electrostatic problem:

[TABLE]

where the source terms can be expressed as

[TABLE]

[TABLE]

[TABLE]

Here, it is understood that $\rho_{-}$ is obtained by solving the self-consistent Kohn-Sham equation for a given distribution $\rho_{+}$ of the atomic cores.

For a $n$-type semiconductor, using the Boltzman distribution, these equations become

[TABLE]

[TABLE]

where $\varphi_{0}$ stands for the asymptotic value of the potential in Region I. Conversely, for a $p$-type semiconductor, we can write

[TABLE]

[TABLE]

Furthermore, it is important to note that deep inside the semiconductor region, the electrostatic potential obeys the one-dimensional Poisson equations:

[TABLE]

[TABLE]

under conditions of $n$-type and $p$-type doping, respectively. In the long-range limit where $\varphi$ approaches $\varphi_{0}$, these equations imply that

[TABLE]

[TABLE]

As explained in Sec. II.3, these expressions are of central utility in describing the bending of the electronic bands and overcome the length scales that characterize electrostatic screening in the depletion region of the electrode.

II.3 Band bending

Equation (II.2) can be solved by implementing a fully self-consistent solution of the electrostatic problem. We plan to implement this method in the continuation of this study. For the moment, we use a simpler implementation to assess the model. The details of this approach are presented below.

To obtain the equilibrium charge-voltage distribution of the system, we start by specifying a total charge for the electrode, from which the potential of the system can be found. To this end, we place a plane of charge $q_{\rm I}$ in Region I, representing the defect charge in the bulk of the semiconductor, and another plane of countercharge $q_{\rm III}$ in Region III, representing the ionic charge of the electrolyte. Accordingly, an explicit charge $q_{\rm II}$ is added to the slab of Region II to fulfill charge neutrality:

[TABLE]

This planar setup provides an accurate approximation of the electrolytic side (Region III) as long as the ionic concentrations ($10^{20}$–$10^{21}$ cm*-3*) in the electrolyte are significantly larger than typical doping concentrations in the semiconductor ($10^{16}$–$10^{18}$ cm*-3*). This means that most of the potential drop takes place in the bulk of the semiconductor (Region I), making a plane of countercharge a reliable representation of the response of the electrolyte (Region III) (the Stern model).

Furthermore, on the semiconductor side, adding a plane of countercharge within Region I does not lead to any loss of generality in the solution of the problem within Region II and Region III. This can be seen by noting that once the Helmholtz plane of countercharges in Region III is set and the asymptotic boundary conditions of Poisson’s equation inside the electrolyte is fixed, the Fermi energy $\epsilon_{\rm F}$ and charge density $\rho_{-}$ at the surface are fully determined; they do not depend on the specific shape of the defect charge profile in virtue of Gauss’ law. In other words, the charge distribution and potential profile on the right hand side provide an accurate description of the interaction of the electrode with solution.

Although the electrostatic profile in Regions II and III is accurate, the potential in Region I is still not a reliable representation of the electrostatics of a semiconductor electrode where Region I consists of an infinite extension of Region II. To accurately describe the electrostatic potential in Region I, a cutoff plane is introduced within the slab as illustrated in Fig. 3, and the electrostatic potential to the left of the cutoff is set to follow the solution of the electrostatic equations of a continuum dielectric. The region to the right of the cutoff (but still within the explicit slab) will be taken as Region II, representing the surface states. The position of the cutoff plane is a user defined value which should correspond to around the inflection of the potential $\bar{\varphi}$. It should be far enough within the slab that a bulk like state has occurred, making a smooth transition from the surface states to the bulk of the semiconductor. From the value of the electrostatic potential and its derivative at the cutoff plane, the Fermi level of the bulk of the semiconductor in Region I can be easily determined from Eqs. (18) and (19):

[TABLE]

with

[TABLE]

[TABLE]

for $n$-type and $p$-type semiconductor, respectively. In these equations, the cutoff position $z_{\rm c}$ represents the location of the frontier between Region I and Region II in Fig. 3b and $\bar{\varphi}$ is the difference between the electrostatic potential of the charged slab and that of the neutral slab, corresponding to the dashed line in Fig. 3a. The bulk potential of the electrode can then be taken as $\varepsilon_{\rm F,I}$.

Finally, to find the equilibrium state of the charge electrode, we impose that the Fermi level of the bulk of the semiconductor must equal the Fermi level of the quantum slab:

[TABLE]

By satisfying this condition, the charge density of the electrode can be calculated as a function of voltage, and the surface state density can finally be obtained as the total charge on the right hand side of the frontier defined by $z_{\rm c}$ divided by the elementary charge. Different algorithms can be used to find the conditions of matching Fermi levels such as a dichotomy procedure. The procedure we use is defined in Sec. III.

This protocol enables us to determine how surface states and adsorption affect the potential profile and capacitance of the system. This approach is, however, limited in a few respects. First, the position of the frontier between Region I and Region II defined by $z_{\rm c}$ may affect the asymptotic value $\bar{\varphi}_{0}$ of the potential $\bar{\varphi}$ describing the overall trend of the potential $\varphi$ across the interface. This variation is, however, small and can be easily evaluated from $\Delta\bar{\varphi}_{0}=\frac{\Delta z_{\rm c}}{{\mathscr{L}}_{\rm I}}\frac{d\bar{\varphi}}{dz}(z_{\rm c})$ with ${\mathscr{L}}_{\rm I}$ being the electrostatic screening length of the semiconductor. Since ${\mathscr{L}}_{\rm I}$ is on the order of 10–103 nm, the sensitivity of $\bar{\varphi}_{0}$ to $z_{\rm c}$ is negligible under relevant doping conditions. Second, the range of charge that can be tested is dependent on the size of the slab. In fact, if the voltage drop between the two countercharge planes is larger than the band gap of the material, unphysical charge transfer by Zener tunneling will take place between the two sides of the slab. In these cases, smaller slabs serve to reduce the voltage drop for a constant Helmholtz charge density. It is important, however, to verify that the Fermi level converges with respect to the size of the slab used for the calculation. Third, the solution is essentially planar within Region I and III. This assumes that both the bulk of the semiconductor and the solution see no variation in the planar directions. For a few applications such as quantum dots, this assumption may need to be revisited. Nevertheless, this is a valid assumption within most applications of interest to first-principles surface electrochemistry.

III Computational details

Density-functional theory calculations are performed using the pw code of the Quantum-Espresso distribution.Giannozzi et al. (2009) As shown in Fig. 4, surface slabs of 1 $\times$ 1 Si(110) and rutile, cristobalite, and quartz ${\rm SiO}_{2}$ (110) and cristobalite and quartz ${\rm SiO}_{2}$ (100) are constructed with a slab width of 5 layers, which is sufficient to give converged Fermi levels within 0.01 eV. The slab is centered in the supercell with a vacuum height of 7 Å to ensure convergence of the atomic forces within a few meV/Å. We use ultrasoft pseudopotentials with the Perdew–Burke–Ernzerhof parameterization of exchange–correlation interactions.Perdew et al. (1996) The cutoffs of kinetic energy of the charge density and electrons are set at 50 Ry and 750 Ry, respectively. The Brillioun zone is sampled with a shifted 5 $\times$ 5 $\times$ 1 Monkhorst–Pack grid and 0.03 Ry of Marzari–Vanderbilt smearing.Marzari et al. (1997)

As explained above, the electrostatic response of the electrolyte interface is modeled using the environ module with the parameterizations developed for water.Andreussi et al. (2012) Several surface configurations of the Si and ${\rm SiO}_{2}$ surfaces are examined. In finding the equilibrium structure for these configurations, the three layers closest to Region I are frozen to create a bulk-like condition. The final relaxed positions are then used in the semiconductor–interface model discussed above. To find the equilibrium charge distribution between the surface states and the bulk of the semiconductor, several partitions of the charge $q_{\rm III}$ between $q_{\rm I}$ and $q_{\rm II}$ are tested for a fixed $q_{\rm III}$. In explicit terms, for each total electrode charge, 11 different partitions were considered: one where 0% of the electrode charge is in Region I and 100% in Region II, one where 10% of the electrode charge is Region I and 90% of the electrode charge is in Region II, and so on. Using a fixed cutoff position, $z_{\rm c}$, we find the charge distribution that minimizes the difference in Fermi level between Regions I and II from Eq. (23). A dopant concentration of $10^{18}$ cm*-3* was used along with a dielectric constant of silicon as $\epsilon_{\rm I}=11.7$ and a dielectric constant of water of $\epsilon_{\rm III}=78.3$. The results of these simulations are presented and discussed in Sec. IV.

IV Results and Discussion

Silicon electrodes represent an important test for the proposed method. Experiment shows that an amorphous layer of ${\rm SiO}_{2}$ forms on the surface of Si in contact with water.Runyan (2013); Gräf et al. (1989) However, the underlying driving force for this oxidation is not clearly understood and can vary with both pH and voltage.Leung et al. (2009); Li et al. (2004) In response, several different terminations of silicon were tested. The results of these calculations can be seen in Fig. 5.

When charge is placed on a solvated silicon slab with planes of countercharge for charge neutrality, all of the electronic charge accumulates at the semiconductor edge.Sup This contrasts with the expected electron distribution throughout a semiconductor electrode with surface states, implying that all the charge on a silicon semiconductor slab would accumulate deep inside the bulk of the semiconductor. As a result, there is no combination of charges $q_{\rm I}$ and $q_{\rm II}$ that equalizes the Fermi levels in Region I and Region II [Eq. (23)]. This shows that no equilibrium charge distribution exists between the surface of silicon and its bulk state, offering insight into the instability of pure silicon in water at low potentials.Nikolaychuk (2014) From this, we conclude that no surface states are likely to form on a pure silicon surface. Thus, some significant contribution from adsorption is expected for the electrification of surface states on a silicon electrode in water.

To gain insight into the role that surface states play on the charge-voltage curves semiconductor–solution interfaces, we use a simple model that consists of describing the charge–voltage behavior of a semiconductor electrode as an ideal Mott–Shottky semiconductor in series with a metal surface state. This gives an overall capacitance of the form $1/{\mathscr{C}}_{\rm I+II}(\Phi)=1/{\mathscr{C}}_{\rm I}(\Phi_{\rm B})+1/{\mathscr{C}}_{\rm II}$, where ${\mathscr{C}}_{\rm I}$ is the capacitance of Region I obtained from Eq. (18), ${\mathscr{C}}_{\rm II}$ is the capacitance of the interface region, and $\Phi_{\rm B}$ is the potential drop across Region I (the Schottky barrier). In this model, ${\mathscr{C}}_{\rm II}$ does not depend on the potential $\Phi$, whereas ${\mathscr{C}}_{\rm I}$ depends on it through $\Phi_{\rm B}$. Furthermore, we describe the relation between $\Phi_{\rm B}$ and $\Phi$ in the vicinity of the flatband potential, that is, for a small amount of charge accumulation at the electrode, as $\Phi_{\rm B}=\gamma_{\rm I}(\Phi-\Phi_{\rm FB})$, where $\gamma_{\rm I}$ is the fraction of the total potential drop that occurs within the semiconductor, limited to values between 0 and 1. The value of $\gamma_{\rm I}$ accounts for the contribution from surface states to the overall electrical response; the value of $\gamma_{\rm I}$ decreases with increasing charge buildup in the surface states. The parameters of the fitted curves, shown in Fig. 5, are reported in Table 1.

Simulations using different surface terminations of silicon show a strong buildup of charge upon increasing the potential, as can be seen in Fig. 5a. This reflects the fact that nearly all the charge is trapped in Region II of the material, as confirmed by the observation that the fraction $\gamma_{\rm I}$ of the potential drop taking place within the semiconductor is low for all of these structures. In general, adding hydroxyl groups to the surface causes the charge to be more distributed within the material than with just an oxygen added. This is supported by noting that the fraction of the potential drop within the semiconductor is higher for structures with hydrogen added than the corresponding structures without hydrogen. One likely explanation for this is that negatively charged adsorbates lead to a trapping of positive charge near the surface. Conversely, adding hydrogen to the surface reduces the electronegativity of the adsorbate and allows for more long-ranged charge distribution inside the semiconductor. This is further shown with the Si+4O*+H* structure, which has the most charge distribution within the electrode. Under typical conditions, however, the type of adsorbate at the silicon surface seems to have a moderate effect on the overall trend; in all these curves, a large initial buildup of charge due to surface states is followed by a much slower buildup dominated by the bulk of the semiconductor. All of the silicon adsorbate materials have a capacitance ${\mathscr{C}}_{\rm II}$ on the order of 10s of $\mu$F/cm2.

In contrast, simulations with ${\rm SiO}_{2}$ terminations show large differences in the resulting charge–voltage response. In general, these structures have a much more distributed charge profile. This is reflected in Fig. 5b by the lowered charge density in comparison to the silicon adsorbate structures. In particular, the rutile ${\rm SiO}_{2}$(110) structure and cristobalite ${\rm SiO}_{2}$(110) structure present a stark difference with the other shown structures. This can be attributed to surface states with a much lower capacitance. This leads to much shallower growth of the charge–voltage curve. It should be noted that the value of ${\mathscr{C}}_{\rm II}$ for Cristabolite (110) is unphysically high; this is not due to the quantum continuum model, but rather indicative of the limitations of a two parameter fitting function. Another important aspect to note is the change in electrode behavior is a function of the exposed surface. Changing from the (100) to (110) orientation for a cristobalite and quartz ${\rm SiO}_{2}$ structure leads to a lower accumulation of charge. This further underscores the importance of the exposed facet in determining the voltage-dependent charge distribution across the interface. In a traditional Mott–Schottky model, the specific surface termination would not change the charge–voltage response.

It should be noted that for surfaces that quickly grow positive with the application of a small amount of potential, it is expected that negatively charged species from the solution would adsorb at the surface. This would result in surface oxidation until a passivating oxide layer forms, providing insights into the experimentally observed formation of an oxide layer when silicon is in contact with water. Runyan (2013); Gräf et al. (1989) For a more complete comparison with experiment, it would be necessary to perform simulations on large-scale amorphous surface terminations under applied voltage and controlled pH. This will be the subject of a study in the continuation of this work.

The results presented here differ from the ideal Mott–Schottky picture by providing a detailed description of charge accumulation at low potential where the surface states dominate the electrochemical properties of the electrode.

V Conclusion

Semiconductor–electrolyte interfaces encompass numerous applications at the frontier of solid state physics and electrochemistry. We have presented a method to embed first-principles calculations of surface states between a Mott–Schottky description of band bending within the semiconductor and the Helmholtz representation of the surrounding electrolyte. We have applied the method on different surface terminations for silicon with a focus on Si and ${\rm SiO}_{2}$ structures. These simulations provide a comprehensive atom-level understanding of the experimentally observed electrification of silicon electrodes in water, suggesting the rapid accumulation of positive charge at the surface of solvated silicon electrodes leading to the formation of an oxide layer that shifts the potential of charge neutrality to more positive voltages and ultimately prevents further oxidation. This method is ideally positioned to examine the low-potential regime where surface state charges dominate the electrification of the electrode in a manner not captured by the Mott–Schottky theory alone. Future work will focus on the implementation and distribution of algorithms to determine the three-dimensional charge distribution between the bulk of the semiconductor, the surface states, and the electrolyte for predicting the structure and response of semiconductor–electrolyte interfaces under electrochemical conditions.

Acknowledgements.

The authors acknowledge primary support from the National Science Foundation under Grant DMR-1654625, and partial support from the 3M Graduate Fellowship and Penn State University Graduate Fellowship.