The relevance of structural variability in the time-domain for computational reflection anisotropy spectroscopy at solid–liquid interfaces

Solid–liquid interfaces are crucial for many technologically relevant processes, ranging from energy storage in batteries to hydrogen production via water splitting [1]. Charges need to be transferred over this interface to drive an electrochemical reaction, yet the structure and the potential distribution can vary significantly both in the time and space domain, but also as a function of an applied external potential [2]. Hence, a detailed understanding of these interfaces is of utmost importance for systematic improvements of electrochemical devices. Significant experimental and computational efforts are therefore dedicated to unravelling the structure of electrochemical interfaces [2–4]. However, most experimental techniques applicable to solid–liquid interfaces under operando conditions are either restricted in structural or temporal resolution, which means that the availability of information on the atomistic scale under realistic electrochemical conditions is limited. Here, electrochemical reflection anisotropy spectroscopy (RAS) is an emerging optical method in the field, allowing for a non–destructive investigation of crystalline surfaces and interfaces providing insights in their atomistic structure [5].

RAS measures the difference in reflectivity for two orthogonal polarisation directions, thus yielding information on the anisotropy of the investigated surface. This anisotropy can originate (almost) exclusively from the investigated crystal surface, as is for instance the case for the (110) surface of fcc metals or the (100) surface of III–V semiconductors in the zincblende structure. The combination of structural and temporal resolution of RAS has made it a popular tool to follow and control epitaxial growth processes [6]. Yet due to its near-normal incidence reflection geometry, RAS can also be straightforwardly applied in electrochemical environments to probe the anisotropy and hence the atomistic structure at the electrochemical solid–liquid interface under applied potentials. This means that electrodes with anisotropic crystal surfaces can be investigated under operating conditions, allowing to relate particular features in a cyclic voltammogram to an in operando measured spectroscopic response. The observed spectra can then ideally be connected to particular structural modifications such as reconstructions or the adsorption of certain atoms and molecules on the surface [5, 7]. Early electrochemical RAS studies used gold as electrode in aqueous electrolytes [8–10]. As Au is not prone to corrosion in mild electrolyte concentrations and moderate potential windows, structural changes can be assumed to be driven by the applied potential in a reversible manner [11, 12].

While RA spectroscopy allows to monitor changes on surfaces and interfaces, the interpretation of such changes in terms of the underlying atomic structure is non–trivial. In fact, such an interpretation can only be obtained by comparison between experimental and computational spectra or the tight correlation with complementary experimental methods. While the latter is often challenging or not possible at all for electrochemical systems, the former requires a rather evolved computational effort, depending on the level of theory. Before density functional theory (DFT) calculations of large supercells became affordable, early studies used perturbative approaches [13]. The currently most well-established approach for such a calculation is an excited state calculation on a random phase approximation (RPA) or even Bethe–Salpeter equation level on top of a DFT ground state [14]. The existing computational studies on RAS so far focused on fixed structures after geometry optimisation [7, 15, 16], thus lacking any information on the structural interface dynamics inherent at an electrochemical interface. These dynamics in RAS are expected to qualitatively differ from infrared spectroscopy probing molecular vibrations, where methods are already established to derive experimental signatures from molecular dynamics (MD) simulations [17, 18]: For most systems it is safe to assume that the bulk electrolyte is isotropic and hence the anisotropy confined to the interface and furthermore, experimental realisations of RAS are typically operating at different energies, in the visible regime ($\gt\!\!1.5$ eV).

In this paper, we therefore address the question on how structural dynamics at the solid–liquid interface impact RAS. As established model system, we investigate the missing row reconstruction of the Au(110) surface in contact with water. We find significant fluctuations of spectral features over the time of the MD trajectory which should be taken into account to realistically model electrochemical interfaces for comparison with experiment.

2.1. Reflection anisotropy spectroscopy

In experimental RA spectroscopy, linearly polarised light at near-normal incidence is used to investigate the surface/interface of interest. The difference in reflectivity along two orthogonal crystal directions in the surface plane is measured and analyzed as a function of the energy of the incoming photons. The observed reflectivity difference is normalised by the overall reflectivity, thus resulting in the following expression:

$$\begin{align} \frac{\Delta r}{r} = 2\frac{r_x - r_y}{r_x + r_y}; r \in \mathbb{C}. \end{align} \tag{ 1 }$$

Here, it has to be noted that the typically determined quantity in experiment, $\Delta r/r$, refers to the reflectivity, which can be related to the reflectance via $\mbox{Re}(\Delta r/r)\approx 1/2(\Delta R/R)$ [19]. In the case of optically isotropic bulk materials, the signal originates exclusively from the (near) surface region. Especially then, RAS is a highly surface-/interface-sensitive probe that allows to identify subtle structural changes, for instance during an electrochemical process.

From a theoretical perspective, the optical properties of a given system are determined by its dielectric function, such that accessing RA spectra by computational approaches necessitates the determination of the surface dielectric function, as described in the following equation [20, 21]:

$$\begin{align} {{\frac{\Delta R}{R}}} = \frac{4\pi d}{\lambda}\left[A\Delta\varepsilon_{s}^{^{\prime\prime}}-B\Delta\varepsilon_{s}^{^{\prime}}\right], \end{align} \tag{ 2 }$$

with $\Delta\varepsilon_{s}^{^{\prime\prime}}$ and $\Delta\varepsilon_{s}^{^{\prime}}$ corresponding to real and imaginary part of the surface dielectric anisotropy, respectively, d to the layer thickness, and λ to the photon wavelength.

Finally, A and B are related to real and imaginary part of the complex bulk dielectric function, $\varepsilon_{b} = \varepsilon_{b}^{^{\prime}}+ i\varepsilon_{b}^{^{\prime\prime}}$:

$$\begin{align} A = \Re\Bigl\lbrack\frac{1}{\varepsilon_{b} - 1}\Bigr\rbrack = \frac{\varepsilon_{b}^{^{\prime}}-1}{\left(1-\varepsilon_{b}^{^{\prime}}\right)^{2}+\left(\varepsilon_{b}^{^{\prime\prime} 2}\right)}, \end{align} \tag{ 3 }$$

and

$$\begin{align} B = -\Im\Bigl\lbrack\frac{1}{\varepsilon_{b} - 1}\Bigr\rbrack = \frac{\varepsilon_{b}^{^{\prime\prime}}}{\left(1-\varepsilon_{b}^{^{\prime}}\right)^{2}+\left(\varepsilon_{b}^{^{\prime\prime} 2}\right)}.\end{align} \tag{ 4 }$$

These quantities need to be computed starting from a ground-state electronic structure which is typically provided on a DFT-level of theory.

2.2. Computational details

The Au(110) surface and its interaction with water was modelled by periodic DFT using the CP2K code package [22]. For all calculations, the Gaussian-and-plane-waves approach was employed, making use of the Goedecker–Teter–Hutter pseudopotentials for describing the electron core interaction [23]. A cutoff of 400 Ry and a relative cutoff of 60 Ry were selected. Exchange and correlation were described by the generalised gradient approximation in the formulation of Perdew, Burke and Ernzerhof (PBE) [24]. Additionally, van–der–Waals interactions were accounted for by the Grimme D3 correction [25].

Before constructing the Au surface, the unit cell of bulk Au was optimised to determine the lattice parameter. Then, a $2 \times4$ Au(110) surface with 6 layers was constructed and a row of Au atoms in the top layer was removed to create the missing row-reconstruction (see figure 1) [13]. After adding 15 Å of vacuum on top of the surface, the slab geometry was optimised with the bottom 2 layers fixed. The calculations were performed under periodic boundary conditions, using a $5 \times 5 \times 1$ k-point mesh constructed by the Monkhorst–Pack scheme. To allow for the investigation of the solid–liquid interface, 5 water layers—corresponding to 30 water molecules—were added on top of the optimised surface to fill the vacuum, resulting in a density close to 1 g cm⁻³.

Zoom In Zoom Out Reset image size

To obtain the corresponding minimum energy configuration, the geometry of the Au-water system was optimised as well. This configuration was then used as starting structure for subsequent ab initio MD (AIMD) simulations. For this purpose a canonical ensemble with a temperature of 300 K, controlled by a Nose–Hoover thermostat, was set up. The AIMD simulations were run for 5000 steps with a timestep of 0.5 fs to equilibrate the system. From the following 5000 steps, 20 configurations (i.e. every 250 steps) were extracted and selected for excited state calculation. To determine the corresponding RA spectra, the RAS module as implemented in the Yambo code [26, 27] was applied.

Yambo requires a highly accurate calculation of the electronic ground state of the system and is directly interfaced with Quantum Espresso (QE) [28]. Therefore, single-point calculations for the selected AIMD snapshots were performed using the plane-wave DFT code QE. As for the previous DFT calculations, the PBE exchange-correlation functional was employed, whereas Optimised Norm-Conserving Vanderbilt pseudopotentials were used. A plane-wave cutoff of 60 Ry with an increased $7 \times 7 \times 1$ k-point mesh using the Monkhorst–Pack scheme was applied. The converged single-point calculations and the corresponding single-particle wave functions were then used as starting point to perform a non-self-consistent field calculation, including an increased number of unoccupied bands. For the water-covered Au(110) surface, which contained altogether 1076 electrons, a total of 750 bands was considered. This number of bands ensures that enough states above the Fermi level are included, which is necessary for the subsequent RAS calculation via the Yambo code.

Finally, Yambo calculations were performed to obtain the RA spectrum in an energy range of 0–5 eV. A real-space cutoff was selected such that only the topmost four layers of the Au slab (including the reconstructed top-layer) contributed to the RA spectrum. As evident from equations (3) and (4), the computational RA spectra have to be normalised by the dielectric function of the bulk material. For this purpose, an experimental bulk dielectric function was applied [29].

It has to be pointed out that RAS calculation for the here presented system sizes became only possible by using a current Yambo code version with enhanced parallelisation. We therefore ported an optimised RAS module to the public development branch of the latest Yambo version (Yambo 5.1) [30]. This module currently features the RPA without local field effects (IP-RPA) to calculate dielectric functions. In the IP-RPA method, the imaginary part of the dielectric function of a slab can be computed as

$$\begin{align} \begin{aligned} \mathrm{Im}\left[\varepsilon_{s}^{xx}\right]& = \frac{8\pi^2 e^2}{m^2\omega^2AL}\sum\limits_{{\mathbf{k}}}\sum\limits_{v,c}|P_{v{\mathbf{k}},c\mathbf{k}}^{x}|^{2}\\ & \quad \times \delta\left(E_{c{\mathbf{k}}}-E_{v{\mathbf{k}}}-\hslash \omega\right), \end{aligned} \end{align} \tag{ 5 }$$

where $P_{v{\mathbf{k}},c{\mathbf{k}}}^{x}$ are transition matrix elements, ω the angular frequency of the light, while A and L correspond to area and height of the cell. These matrix elements consist of the momentum operator $\textbf{{p}}$ and the commutator containing the non-local elements of the pseudo-potentials V^nl , i.e. $\textbf{{p}}+i[V^{nl},\,\textbf{{r}}]$. To calculate the RA spectrum stemming from a specific surface structure on a slab, rather than including the whole slab, a so-called real-space cutoff can be introduced [31]. By employing such a cutoff θ, the transition matrix elements transform into

$$\begin{align} \tilde{P}_{v{\mathbf{k}},c{\mathbf{k}}}^{x} = \langle v{\mathbf{k}} |\theta\left(z\right)\left[{\textbf{p}}^{x}+i\left[V^{nl},{\textbf{r}}^{x}\right]\right] |c{\mathbf{k}} \rangle, \end{align} \tag{ 6 }$$

and we can rewrite equation (5) as

$$\begin{align} \begin{aligned} \mathrm{Im}\left[\varepsilon_{s}^{xx}\right]& = \frac{8\pi^2 e^2}{m^2\omega^2AL}\sum\limits_{{\mathbf{k}}}\sum\limits_{v,c}\left[P_{v{\mathbf{k}},c{\mathbf{k}}}^{x}\right]^{\ast}\tilde{P}_{v{\mathbf{k}},c{\mathbf{k}}}^{x} \\ & \quad \times \delta\left(E_{c{\mathbf{k}}}-E_{v{\mathbf{k}}}- \hslash \omega\right). \end{aligned} \end{align} \tag{ 7 }$$

This updated implementation, allowing for a significant speedup and larger system sizes, will become part of one of the next Yambo releases.

For a typical RAS experiment with a commercial spectrometer, the measured spectrum corresponds to the anisotropy signal averaged over a macroscopic surface area (several mm squared) and integrated over a time of at least 10 ms. To investigate how the dynamics of the solid–liquid interface of a nm-sized supercell impact the spectroscopic response, we first computed RA spectra along a 2.5 ps-long part of an MD trajectory of the Au surface in contact with water. The spectral response corresponding to the different MD snapshots and the resulting average are depicted in figure 2. All spectra show elements of the experimentally observed main features of the reconstructed Au(110) surface, in particular the pronounced negative intensity at 2.7 eV [5, 8, 9]. Interestingly, the single spectra show certain differences to the time-average, especially in the high-energy part beyond 2.5 eV (see figure 2). It can also be observed that larger deviations are conserved over a longer time period, here for frames 9500–10 000, i.e. over a time-period of at least 250 fs. To correlate changes in surface structure and RA spectra, the side view of the corresponding MD configuration and the resulting electrostatic potential are shown below each RA spectrum (see figure 2). From this comparison, it becomes evident that the frames that show an increased deviation from the average RA signal (especially frames 9500–10 000), correspond to those configurations with larger deviations from the average potential. This makes it clear that the structural changes that are at the origin of the variation in electrostatic potential affect the resulting RA spectra. However, it also has to be pointed out that it is not possible to directly derive the evolution of the RA spectrum from observed changes in the surface structure (see also the discussion below for the water-free case).

Zoom In Zoom Out Reset image size

To further elaborate on the fluctuations induced by the water layer, a comparative MD simulation without water was performed and the resulting RA spectra were investigated. Here, the MD trajectory [32] of the plain Au(110) missing-row reconstruction was also analyzed with respect to the time evolution of the RA spectra (see figure 3). As in the case of the water containing surface, the electrostatic potential for the respective snapshots was determined (see figure S2 in the SI). Interestingly, for the plain surface the differences of the electrostatic potential from the average are generally more pronounced than in the presence of water. However, for the corresponding RA spectra this does not necessarily translate in strong variations. Also here, a persistence of larger deviations from the average RA signal can be observed, in this case from frames 8250–9000, which translates to a time-span of 375 fs. The reason for these deviations is most likely the characteristic frequencies of the Au surface represented in the MD trajectory. This emphasises that computational RA spectra from single arbitrary snapshots of an MD trajectory will typically not capture effects due to surface dynamics. A comparison of MD-averaged spectra with and without water indicates rather small differences, which are most pronounced in the intermediate energy range from about 2.5–3.5 eV (see figure 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A comparison to the RA spectrum of the relaxed, water-free equilibrium structure of the Au(110) missing row reconstruction in vacuum with both average spectra shows significant differences. Overall, the intensity of the averaged RA spectra amounts to only about $1/4$ of the intensity obtained for the geometry-optimised, minimum energy surface displayed in figure 4. Interestingly, the low-energy, positive anisotropy feature between 1.0 and 2.2 eV is much more pronounced in the case of the average spectra. A comparison of the averaged spectra with and without water shows that the water-decorated surface leads to a spectrum with an intensity further reduced by about 20%, indicating a slightly reduced ordering of the Au surface. While most spectral features are almost conserved, the spectrum from the slab with water shows more structure in the energy range around 4 eV.

The observed intensity decrease for the MD-averaged spectra can be understood as a consequence of the surface dynamics, which makes the originally present anisotropy less pronounced, thus resulting in a reduced RA signal. This is an important finding, as differences between experimental and calculated spectra are often interpreted as stemming from not fully-covered areas on the surface. Our findings, however, strongly indicate that such differences may also be a consequence of neglecting the surface dynamics: While the geometry-optimised surface structure at 0 K represents ideally the absolute minimum energy configuration of the energy hypersurface with respect to the atomic coördinates, the kinetic energy in the MD simulation will drive the system out of this minimum and can therefore reduce anisotropies of the structure.

A remaining question is now the origin of the differences between the averaged RA spectra of the Au surface with and without water. These may either be due to a direct contribution of the water molecules to the RA spectra or due to an indirect impact via a water-induced alteration of the Au surface, i.e. via changes in structure and dynamics of the interface. To resolve this issue, the RA spectra of the water containing structure were recalculated by applying the real space cutoff such that potential contributions of water molecules were excluded.

As can be inferred from figure 5, the average RA spectra obtained for the MD simulation including water is almost independent of the direct water contribution, similar to what was found for a static calculation of a semiconducting system in contact with water [16]. This does, on the other hand, mean that the differences between the average spectra obtained for the MD simulation with and without water are an evidence of the impact of the water molecules on real-space and electronic structure of the gold surface. The latter one can also originate from the electric field inducing a linear electro-optic effect [16]. In other words, the difference between the average spectra of the two MD simulations appear to be a signature of the solid–liquid interface.

Zoom In Zoom Out Reset image size

Finally, it has to be pointed out that the here presented model system does not perfectly agree with experimentally observed data of the full electrochemical system [5, 8–10]. This, however, is likely to be a consequence of the limited system size as well as the absence of electrolyte ions. For a full comparison to experiment, an increased number of Au layers might have to be investigated (see figures S1 in the SI), while a better statistical averaging (more frames/longer simulation times) could also further improve the comparison. Yet, the focus of the present work lies in emphasising the impact of the surface dynamics on the RA signal, which is already evident for the here presented system size.

In this paper, we investigated the impact of structural variability in the time-domain on RAS by computing RA spectra along an MD trajectory of the Au(110) surface, featuring the frequently observed missing-row reconstruction with and without considering the presence of water. In both cases the dynamically changing structure of the Au surface is found to result in a significant reduction of the observed optical anisotropy, while at the same time the main RAS peak becomes more pronounced. This finding indicates that surface and interface dynamics are important when RA spectra are analysed with respect to computational data. Indeed, differences in overall signal intensity between experimental and computed RA spectra are frequently attributed to macroscopic properties of the surface, whereas a microscopic origin in terms of surface/interface variability was so far not considered. Furthermore, the observed differences in the RA spectra between the plain surface and the water-containing setup show that the presence of water molecules changes the surface/interface structure in a well-ordered manner. This points to the high sensitivity of RAS with respect to the solid–liquid or solid-electrolyte interface (SEI), while at the same time it becomes evident that for the interpretation of such interfaces, static calculations that only consider the geometry-optimised, minimum energy structure will not necessarily give the full picture. Finally, we want to emphasise the potential of RAS as a method to study long-standing problems in applied electrochemistry, such as the SEI growth in batteries. When single-crystalline model systems can be prepared, RAS is capable of deciphering the processes that are determininig the initial state of the SEI growth, which are still not well understood in many battery systems. With the more complex electrolytes including large molecules, however, this will be a challenging undertaking. Yet, also for aqueous electrochemistry, a closer approximation of real electrochemical systems will require the analysis of MD trajectories that include both ions and applied potentials [4].

This work was funded by the German Research Foundation (DFG) under Project No 434023472 and the Cluster of Excellence EXC 2154/'Post-Li Storage', Project No. 390874152. The German Bundesministerium für Bildung and Forschung (BMBF), supported this work with the Project 'NETPEC' (No. 01LS2103A). The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through Grant No. INST 40/575-1 FUGG (JUSTUS 2 cluster). Part of the simulations were performed on the national supercomputer Hawk at the High Performance Computing Center Stuttgart (HLRS) under the Grant No. SPECSY/44227. The authors thank Conor Hogan and Davide Sangalli for assistance with re-implementing the RAS module in Yambo.

The AIMD trajectories that support this study are openly available at the following URL/DOI: https://doi.org/10.5281/zenodo.10455322.