Band alignment at semiconductor-water interfaces using explicit and implicit descriptions for liquid water

In this work we study and contrast implicit solvation models against explicit atomistic, quantum mechanical models in the description of the band alignment of semiconductors in aqueous environment, using simulations based on density functional theory. We find consistent results for both methods for 9 different terminations across 6 different materials whenever the first solvation shell is treated explicitly, quantum mechanically. Interestingly this first layer of explicit water is more relevant when water is adsorbed but not dissociated, hinting at the importance of saturating the surface with quantum mechanical bonds. Furthermore, we provide absolute alignments by determining the position of the averaged electrostatic reference potential in the bulk region of explicit and implicit water with respect to vacuum. It is found that the absolute level alignments in explicit and implicit simulations agree within $\sim 0.1-0.2$ V if the implicit potential is assumed to lie 0.33 V below the vacuum reference level. By studying the interface between implicit and explicit water we are able to trace back the origin of this offset to the absence of a water surface dipole in the implicit model, as well as a small additional inherent polarization across the implicit-explicit interface.

In this work we study and contrast implicit solvation models against explicit atomistic, quantum mechanical models in the description of the band alignment of semiconductors in aqueous environment, using simulations based on density functional theory. We find consistent results for both methods for 9 different terminations across 6 different materials whenever the first solvation shell is treated explicitly, quantum mechanically. Interestingly this first layer of explicit water is more relevant when water is adsorbed but not dissociated, hinting at the importance of saturating the surface with quantum mechanical bonds. Furthermore, we provide absolute alignments by determining the position of the averaged electrostatic reference potential in the bulk region of explicit and implicit water with respect to vacuum. It is found that the absolute level alignments in explicit and implicit simulations agree within ∼ 0.1 − 0.2 V if the implicit potential is assumed to lie 0.33 V below the vacuum reference level. By studying the interface between implicit and explicit water we are able to trace back the origin of this offset to the absence of a water surface dipole in the implicit model, as well as a small additional inherent polarization across the implicit-explicit interface.

I. INTRODUCTION
Atomistic insight into the structure, composition and properties of solid-liquid interfaces is paramount for our understanding of the stability or performance of materials in many technological devices, such as chemical sensors, batteries or fuel cells. Very often, such detailed knowledge can only be obtained from atomistic simulations or from a combination of theory and experiments (e.g. interpretation of spectroscopic measurements). To date, the main challenges for such simulations are the enormous size of the compositional and configurational space, in particular for electrochemical solid-liquid interfaces, where the stability of a certain interface composition is not only influenced by the electrochemical potential of adsorbates but also by the energetics of the long-range space-charge layer, induced by the screening in the electrolyte over distances that are often out of reach for purely atomistic interface models due to the limits in cell size and number of atoms 1 .
Recently, several studies have shown that coupling density functional theory (DFT) to implicit solvation models can yield accurate results for the solvation energies [1][2][3][4][5][6][7] . In these cases, the explicit electrolyte solution (e.g. water with ions) is substituted by a meanfield description of the solvent, where the liquid is modelled via a polarizable continuum or via joint density-functional theory 2,3,5,6,8,9 . Computational costs are highly reduced not only because of smaller quantum mechanical systems but also because thermodynamic averaging, which is mainly necessary for the correct description of the liquid, is substituted by the response of the continuum model.
Reasonable agreement between explicit and implicit descriptions of water as well as with experiment also typically found for first principles simulations of metallic slabs with respect to interfacial structure, capacitance, potential of zero charge and interface energetics 1,3-5,10-16 . This points to rather unspecific interfacial water structures. Studies have shown that water ordering on metals depends mainly on the relative size of surfacewater and water-water interactions [17][18][19] , where significant ordering occurs only for very strong or very weak surface-water interactions, where the surface acts as a template or as a mere boundary supporting 2-dimensional intermolecular ordering via H-bonds.
For multicomponent semiconductor-water interfaces, on the other hand, rather strong surface-water interactions are more common 20-27 inducing very interface-specific properties of the solvation shell, where the applicability of implicit solvation is yet to be tested. As an example, we plot the time averaged density of O and H atoms as obtained from an ab-initio molecular dynamics (AIMD) trajectory of CdS(1010) in explicit water in Fig. 1.
The inhomogeneous accumulation of O and H (red and blue isosurfaces) at the interface are clear signatures of immobile, interfacial H 2 O molecules, whose properties are expected to be significantly different from those of bulk water molecules.
In this work we test protocols to simulate the band alignment of semiconductors accurately leveraging implicit solvation models, and validate them against explicit simulations for the following semiconductors in water: (rutile) r-TiO 2 (110) and CdS(1010) with molecularly adsorbed water, GaN(1010) with dissociatively adsorbed water and a-TiO 2 (101), GaAs (110) and GaP(110) with either of these two possibilities. All these terminations were previously found to be stable or meta-stable 28 . We will show that it is necessary to describe the first solvation shell explicitly in order to capture interfacial potential drops at semiconductorwater interfaces accurately within an implicit model, as similarly argued in Refs. 29 and 30.
Furthermore, we will provide an estimate for the position of the flat electrostatic potential in the bulk of the implicit model that is approximately 0.33 eV below the vacuum level (absolute potential scale) for the SCCS implicit solvation model 2 . This work provides best practices for simulations of semiconductor-water interfaces in implicit aqueous environment and also new insights in the effects and properties of interfacial water layers and to what extent they can be described by coarse grained continuum models, which we will tackle in the near future.

II. METHOD
The approaches adopted for the determination of band alignment at semiconductor-water interfaces are schematically illustrated in Fig. 2. For both interface models, explicit (e) and implicit (i), we simulate an explicit semiconductor slab with potential adsorbates (e.g. explicit water molecules) immersed in an explicit aqueous or an implicit environment. In order to determine the band alignment on an absolute scale (e.g. w.r.t. vacuum or SHE) three potential offsets need to be known: They are marked as e1-e3/i1-i3 in Fig. 2 and correspond to the position of the bands, e.g. the bottom of the conduction band c , w.r.t. the average bulk electrostatic potential V SC (e1/i1), the alignment of V SC w.r.t. the average bulk water potential V W (e2/i2) and the position of V W w.r.t. an absolute reference, most conveniently vacuum (e3/i3).
Whereas bulk simulations (e1/i1) are equivalent for explicit and implicit models, the alignment of V SC with respect to the water average V W (e2/i2) will be different. This is due to the fact that the main electrostatic contribution of the average potential in explicit water -the quadrupole contribution 31,32 -is missing in typical implicit models and due to different interface models used, e.g. with varying amount of explicitly treated interfacial water. In particular, it is yet unclear which and how many water molecules to treat explicitly, and whether similar results as in explicit simulations can be obtained. In order to answer both of these questions we analyse the properties of individual water molecules and interfacial water layers from all-explicit ab-initio molecular dynamics (AIMD) simulations for the 9 aforementioned water-semiconductor interfaces. Based on these results, we construct representative sampling sets of interface models with 0 to 3 explicit interfacial water layers, by removal of appropriate water molecules and using 100-150 random snapshot of the all-explicit AIMD trajectories per studied interface. Subsequently, the potential offsets are determined by SCF calculations based on density functional theory from the average electrostatic potential alignment V ex SC w.r.t. V ex W for the all-explicit snapshots and V im SC w.r.t. V im W for the interfaces with reduced number of molecules embedded in an implicit environment.
It will be shown that the potential offset V ex W −V im W is independent of the interface termination and amount of explicit water beyond the first solvation shell, supporting the suggestion that V im W can indeed be used as a universal reference level. For the final determination of the alignment with respect to an absolute level (vacuum) (e3/i3), V ex W is determined by simulation of the alignment of the potential of an explicit water slab in vacuum. For the implicit model we provide two pathways, namely by a fit to reproduce all-explicit calculations as well as a dedicated study of the potential alignment of explicit water slabs in implicit water, with both leading to consistent results.
It is important to note that all potential offsets (1-3) are not necessarily transferable between different DFT codes or pseudopotentials; however, the alignment of energy levels of the band edges w.r.t. vacuum remains unaltered, due to a cancellation in the differences (see also the comment below in section III).

III. COMPUTATIONAL SETUP
The fully equilibrated molecular dynamics trajectories used in this work were obtained in Ref. 28 and Ref. 32 by Car-Parrinello AIMD simulations in the NVT ensemble at 350 K with a time step of 0.5 fs and using a liquid water-adapted rVV10 functional 33 , that accounts for nonlocal van der Waals interactions 34,35 . Timeframes were carefully chosen for the 9 selected interface terminations such that no drifts of potential or change of interfacial composition are present in the trajectories. For each interfacial system a thermal distribution of structures was approximated by a random selection of 100-150 snapshots from a subset of structures corresponding to each 50th MD step (∆t = 25 fs). All reported results here are based on averages, obtained from SCF calculations of these snapshots and derived structures using the ENVIRON 36 module of Quantum ESPRESSO 37 and PBE 38 as exchange-correlation functional. As we are mainly interested here in demonstrating and testing the consistency of calculations for potential offsets at interfaces within implicit solvation models, we restrict the analysis to PBE, which underestimates band gaps but has been shown to reproduce potential offsets in good agreement with more advanced functionals and

IV. RESULTS
A.

B. Interfacial water properties and interface model construction
Surface-water interactions are generally very strong for semiconductors, such that theoretical understanding of hydration is typically separated into properties of primary and secondary hydration shells 22,24 , where the latter is dominated by generic water-water interactions more amenable to continuum models. In contrast, the primary hydration shell The distinct behaviour of molecules within the first solvation shell is evidenced in the systems studied here by the increased localization of water molecules in direct vicinity of surfaces, e.g. for CdS in Fig. 1. In order to obtain a more quantitative picture we performed a combined analysis of cumulative and molecule-specific descriptors, in particular the planar average of the water-related oxygen density (Fig. 3 a) and a molecular entropy descriptor S (Fig. 3 c), with: The probabilities p i describe the sampling of the degrees of freedom of the H 2 O molecule and are estimated from histograms with bins i in the 6-dimensional molecular configuration space as visualized in Fig. 3  with dΩ = sin θdΦdθ -leads to artificial anisotropies, which is why we decided instead to use as bins 320 equal, equilateral triangles of an accordingly decomposed spherical surface (see interfacial water molecules are not exchanging place with bulk-like water. Hence, molecules can be associated with individual water layers and used as probes for layer-specific properties. According to the molecular entropy-based analysis -as shown clearly for r-TiO 2 with molecularly adsorbed water in Fig. 3 c) -only the molecules of the first water layer -if associated with the first, pronounced oxygen density peak (Fig. 3 a) -show a distinctly different behavior than bulk-like water in the center of the water slab. H 2 O molecules in the additional layers, which become increasingly ambiguous due to reduced density oscillations, are found to be not significantly different from bulk ones (Fig. 3 a). Similar results are found for the other systems (see Figs. SF4 -SF8 in the SI), with the first density peak being significantly more pronounced for the unpassivated surfaces (i.e. the pristine surfaces with molecularly adsorbed water).
As a result, we chose to use the number of water molecules associated with the first density peak N l (for non-dissociated water) to define a water layer. In this setup, n water layers correspond to the collection of the n × N l water molecules closest to the surface. In most cases, but not always, this coincides with minima in the observed O-density and is associated with a "discontinuous jump" in the average distance to the surface between the n × N l th and the n × N l + 1st water molecule. Typical interfacial structures for r-TiO 2 as studied below are plotted in Fig. 4. We note in passing that -although well defined in the simulations here -a separation of the first water layer and the oxide surface can be ambiguous 25,26 . C.
Step e2/i2: Potential offsets at semiconductor/water interfaces As discussed in the methods part, the average electrostatic potential inside bulk explicit water is qualitatively and quantitatively different from an implicit electrostatic potential due to the presence/absence of electronic and nuclear charge density. The hypothesis that both models are equally suited to determine band alignment can be tested by studying the potential offset of bulk explicit water and implicit water. For each interfacial system studied,       It is worth noting that anatase TiO 2 with molecularly adsorbed water is a problematic case for the implicit procedure. The SCF does not converge for a single calculation using the purely implicit model with the ENVIRON code also after tweaking typical parameters as reported in Ref. 44 and in private communications with other researchers (see also Figs. SF16 and SF17 in the SI). In a few cases it converges for a single explicit water layer, based on which we report averages. We speculate that this is due to surface Ti ions, whose electronic structure is extremely sensitive to the environment, making self-consistency more challenging. The reported results need to be taken with care due to possible biases for the small converged subset of calculations. Whether such convergence problems might Table I. Position of the macroscopic average of the electrostatic potential V ex W in the center of explicit water slabs w.r.t. the vacuum or the implicit solvation potential (SCCS) (c.f. Fig. 7).
The system "real slab" refers to an AIMD trajectory a water slab in vacuum (Ref. 32). The other systems, denoted by "random set", refer to the random water slabs created according to the prescription in the text. "Polarization pot." denotes the potential drop due to the polarization charges of the implicit model as e.g. plotted in Fig. 8 a). universally hint towards the necessity of including explicit water beyond the first solvation shell is to be tested in the future.

D.
Step e3/i3: Relative and absolute alignment of explicit and implicit water It follows from the results in the previous section that V im W can be used in the same way as V ex W as a potential reference. In particular, if V ex W is known on an absolute scale V im W can be determined by a shift with ∆V W . The absolute position of V ex W is determined most conveniently from simulations of explicit water slabs in vacuum (trajectories already used in Ref. 32) and found to be -3.30 V here (Fig. 7 a and  The average ∆V W from the simulations of Fig. 6 a) with at least one explicit water layer is 2.97 ± 0.05 V; if the results with only one explicit layer are averaged (ignoring the results for a-TiO 2 (mol)) the result is 3.00 ± 0.10 V. At the same time the polarization potential contribution to the overall electrostatic potential offset varies much more strongly than ∆V W across all interface types and number of explicit water layers (see Fig. 6 a, b). This indicates The origin of the implicit polarization for random water interfaces seems related to the structure of a 'random' water surface, as model charge distributions inspired by the observed H and O densities (red solid and turquoise lines in Fig. 8 a) can explain the observed potential drop (Fig. 8 b). They derive from the more smeared out H distribution within the water molecule and the interface construction protocol that imposes a hard-wall for the O density, and implicit polarization charges that follow similar spatial features but with opposite charges. Fig. 8 a) also includes an analysis of force differences with respect to the periodic bulk calculations used for slab creation (red and green solid lines). The inclusion of mirror-symmetric results removes the dipolar potential drops across slabs in average, however, individually, most structures exhibit a macroscopic field across the slab inducing non-vanishing force differences |∆F | in the center of the water slab (note O-related forces are approximately twice as large as H-related ones in the center of the slab, in agreement with electrostatic considerations). The largest force differences, however, are observed for water molecules directly at the explicit/implicit interface. How far such differences might induce differences between molecular dynamics results of pure explicit and explicit/implicit models is unknown, and will be studied in the future but it goes beyond the scope of this work.

E. Semiconductor level alignment w.r.t. SHE
With the above relative alignments and the SHE potential set to the experimental value at -4.44 V w.r.t. vacuum, it is possible to compare conduction and valence band alignments for all interface models and all simulated materials on an absolute scale, yielding Fig. 9.
As mentioned before, apart from the case of a-TiO 2 (mol) where results might be unreliable, no significant increase in accuracy for the potential alignment is observed when including water layers beyond the first solvation shell. The differences to the explicit systems are of the order of 0.1-0.2 V which is comparable to the accuracy of the sampling statistics. On the other hand, simulations without inclusion of explicit water might very likely exhibit alignment errors of 1 V or more which should be kept in mind in future band alignment studies, as this has significant impact on whether a system can support e.g. water splitting reactions or not. Note that we count the explicitly treated dissociated water layer as water layer. We think the observed pattern reflects a general principle that interfacial electrostatic potential drops are captured correctly whenever the QM-implicit interface exhibits no unsaturated dangling bonds which can influence closeby solvent molecules in a very specific way, not captured in implicit models.

V. SUMMARY AND DISCUSSION
In this study we showed that the band alignment of semiconductor slabs in the SCCS implicit model reproduces explicit results up to an accuracy of ∼ 0.1 − 0.2 V, provided that the amount of explicitly simulated water molecules is chosen appropriately. In particular, all systems considered, agree very well with fully explicit simulations, when the chemically interacting interface molecules or their fragments are treated explicitly, which corresponds to the first dissociated or undissociated water layer above the unpassivated surface, in agreement with the findings in Refs. 29 and 30. This finding also agrees with the discussion in section C of the SI, where we observe that the unpassivated surface has interfacial water molecules with properties clearly different from bulk-like.
Furthermore, we studied the explicit-implicit water interface and showed that the electrostatic potential reference inside the implicit region does not correspond to the absolute potential reference (vacuum above the solution), but is shifted by approximately -0.33 V, which puts the SHE level at -4.11 eV w.r.t. the implicit level of the SCCS model. This offset is mainly related to the absence of an explicit water surface dipole contribution in the implicit model, and to a smaller extent to some generic polarization of the implicit model across the implicit/explicit interface. The relative insensitivity of the overall level alignment to the materials and amount of explicit water with at the same time strongly varying implicit-model-related polarization contributions ( Fig. 6 and Figs. SF9 -SF17 in the SI) suggests that the SCCS model is indeed able to mimic accurately all generic solvent-related screening properties at explicit/implicit water interfaces.
These results overall suggest that simulations of electrochemical interfaces can be performed with accuracies close to the all-explicit calculations at significantly reduced cost. We speculate that structural models with the QM-implicit boundary dominated by H, OH or H 2 O species, as obtained with included explicit water might also be favourable for consistent energetics across different systems due to chemical and structural similarity.
The final step to develop an accurate simulation protocol for solid-liquid interfaces with implicit models will be to assess and validate ab-initio molecular dynamics simulations for mixed explicit/implicit models, something that we will address in a follow up paper.