Cooperation and Environment Characterize the Low-Lying Optical Spectrum of Liquid Water

TL;DR

This study uses subsystem time-dependent density functional theory to analyze liquid water's optical spectrum, explaining features like broadening, shifts, and excitonic effects due to environmental interactions.

Contribution

It provides a detailed theoretical explanation for the optical spectral features of liquid water, including environmental and many-body effects, which were previously elusive.

Findings

Broader joint density of states causes red-shifted Urbach tail.

First solvation shell induces blue shift of the absorption peak.

Many-body excitonic effects influence low-frequency spectral weights.

Abstract

The optical spectrum of liquid water is analyzed by subsystem time-dependent density functional theory. We provide simple explanations for several important (and so far elusive) features. Due to the disordered environment surrounding each water molecule, the joint density of states of the liquid is much broader than that of the vapor. This results in a red shifted Urbach tail. Confinement effects provided by the first solvation shell are responsible for the blue shift of the first absorption peak compared to the vapor. In addition, we also characterize many-body excitonic effects. These dramatically affect the spectral weights at low frequencies, contributing to the refractive index by a small but significant amount.

Equations

$$\chi =\sum _I^{N_S}{\chi }_I^c,$$

$${\chi }_I^c=\underset{\text{one body}}{\underbrace{{\chi }_I^u}}+\underset{\text{many body}}{\underbrace{\sum _{J\ne I}^{N_S}{\chi }_I^u{K}_{IJ}{\chi }_J^c,}}$$

Full text

DFT, embedding, liquid water, optical spectrum

Water is the most important liquid on Earth. Thus, understanding its optical spectrum is of pivotal importance. Although there has been tremendous progress in this area, difficulties arise in the ab-initio models of the liquid because large simulation cells need to be employed due to its disordered nature. In turn, this either forces the use of approximate methods on large cells, or accurate methods on cells that are too small. As a result of this limitation, there are several open questions and interesting features of the optical absorption spectrum of water that are yet to be fully explained. In this work, we explore two themes: (1) Many-body excitonic interactions between the water molecules, and (2) coupling of the first absorption band to the nuclear degrees of freedom describing the liquid structure.

By many-body effects, we mean the effects that arise when single water molecules (single bodies) in the liquid interact with each other (other bodies) both in the ground state as well as in their excited states. Although related, this definition is different in spirit from the kinds of many-body effects that a Bethe-Salpeter Equation (BSE) treatment would recover. In the latter, many-body refers to electron–hole interactions. We set out to investigate how many-body interactions affect the optical spectrum and other related quantities (such as the refractive index). These can be cooperative or anticooperative in nature.

Many-body effects have been discussed before in terms of the screening properties of the bulk in the computation of self-energies for GW calculations Garbuio et al. (2006); Ziaei and Bredow (2016). It was found that screening is independent of the particular configuration of water considered. Thus, it can be inferred that it is not affected by the structure of the environment surrounding the water molecules. It has also been shown that for ice, cooperative many-body effects in the form of excitonic couplings increase the oscillator strength of low-lying excitations and are responsible for an increase of the index of refraction with increasing pressure Pan et al. (2014).

In addition, we also set out to investigate coupling between the first absorption band and the structure of the liquid. This helps us understand what influences the peak position and shape (broadening). The underlying reasons for a blue shift of the first absorption band of water from 7.4 eV in gas phase Chan et al. (1993) to 8.3–8.6 eV in the liquid phase Elles et al. (2009); Heller (1974); Hayashi and Hiraoka (2015) are still matter of debate. G0W0 calculations based on PBE orbitals predict a shift of about 0.7 eV and attribute it to excitonic effects Ziaei and Bredow (2016); Hahn et al. (2005) (i.e., the isolated water molecule features a more localized exciton than that of the liquid phase). Calculations based on clusters, instead, define the blue shift as a result of electrostatic embedding Hermann et al. (2008) (a sort of confinement effect). This view is corroborated by excitonically coupled coupled cluster (CCSD) and semilocal time-dependent Density Functional Theory (TDDFT) calculations showing that the first absorption band is composed of localized states Mata et al. (2009); Lu et al. (2008). A recent analysis of the non self–consistency of the previously mentioned G0W0 results uncovered some possible biases that are carried over from using DFT orbitals in the GW+BSE part of the calculation Blase et al. (2016) putting into question the previously predicted exciton sizes of liquid and vapor.

We set out to attack these open questions with simulations based on subsystem DFT Jacob and Neugebauer (2014); Wesolowski et al. (2015); Krishtal et al. (2015a). The electron density is decomposed into subsystem contributions, $\rho({\mathbf{r}})={\sum_{I}^{N_{S}}}{\rho_{I}({\bf r})}$, where $N_{S}$ is the number of subsystems. The subsystem densities are recovered variationally solving subsystem-coupled Kohn–Sham equations, $\left(-\frac{1}{2}\nabla^{2}+v_{s}^{I}({\mathbf{r}})+v_{\rm emb}^{I}({\mathbf{r}})\right)\phi_{i}^{I}({\mathbf{r}})=\varepsilon_{i}^{I}\phi_{i}^{I}({\mathbf{r}})$, where $v_{s}^{I}$ is the Kohn–Sham potential of the isolated subsystem evaluated at ${\rho_{I}({\bf r})}$. The embedding potential, $v^{I}_{\rm emb}$, contains exact Coulomb interactions with surrounding subsystems, as well as functional derivatives of nonadditive exchange–correlation (xc) and nonadditive noninteracting kinetic energy functional ($T_{s}$) Jacob and Neugebauer (2014); Wesolowski et al. (2015); Krishtal et al. (2015a). To access information about excited states, subsystem TDDFT can be formulated either in frequency Casida and Wesolowski (2004) or time domain Krishtal et al. (2015a).

Subsystem TDDFT allows us to approach larger supercells than before and naturally gives us access to many-body effects without compromising the accuracy of the model. Many-body effects arise as inter-subsystem static and time-dependent interactions Krishtal et al. (2015a); Neugebauer (2010); Krishtal et al. (2015b). These perks are further complemented by the ability of subsystem DFT to naturally localize the electronic structure of the subsystems while still allowing electron density overlap. We have showed, for example, that it can quantitatively reproduce the dynamics and structure of liquid water because the too strong hydrogen bond resulting from self-interaction in semilocal xc functionals, is offset by the subsystem-local electronic structure Genova et al. (2016). While we do not expect that subsystem electron density localization be a generally beneficial feature of the model, it is beneficial in the limit of simulating systems composed of noncovalently bound subunits. Such as molecular liquids, crystals, and layered periodic systems Genova et al. (2017).

Simulations were carried out with the package embedded Quantum-ESPRESSO (eQE) Genova et al. (2017), employing the PBE xc functional and ultrasoft pseudopotentials (O.pbe-rrkjus.UPF/H.pbe-rrkjus.UPF from the Quantum-ESPRESSO PP Library Rappe et al. (1990)). We employ plane wave kinetic energy cutoffs of 50 and 500 Ry for the waves and charge densities, respectively. Real-time subsystem TDDFT was implemented in the linear response regime, applying a $\delta$ electric field perturbation Yabana and Bertsch (1996) of 0.02 Ry/Å, a time step of 1 as, and assuming the adiabatic approximation for all the density-dependent potentials (xc as well as nonadditive $T_{s}$). The time-dependent KS orbitals were propagated for 50,000 steps with a Crank-Nicholson propagator, totalling 50 fs of electron dynamics carried out for each of the three Cartesian directions. We considered a total of 8 snapshots of a subsystem DFT based Born-Oppenheimer dynamics of bulk liquid water modelled by 64 independent water molecules in a cubic box ($a=12.4278$Å, presented in a previous work Genova et al. (2016)). We employ the LC94 Lembarki and Chermette (1994) nonadditive $T_{s}$ functional, which was shown to satisfactorily reproduce energy surfaces of CCSD(T) benchmarks for water dimer as well as the structural parameters of the liquid (e.g., radial and angular distribution functions) Genova et al. (2016).

The real-time subsystem dipole change, $\delta\mu_{I}(t)=\int{\mathbf{r}}\left(\rho({\mathbf{r}},t)-\rho({\mathbf{r}},t=0)\right){\mathrm{d}}{\mathbf{r}}$, is Fourier transformed to frequency space to yield the isotropic dipole polarizability, $\alpha_{I}(\omega)$ which is related to the frequency-dependent dielectric constant Hayashi and Hiraoka (2015); Bertsch et al. (2000), $\epsilon(\omega)=\left[1+4\pi\langle\alpha(\omega)\rangle\right]^{-1}$, where the average is carried out over all water molecules for all snapshots, $\langle\alpha(\omega)\rangle=\langle\frac{1}{N_{S}}{\sum_{I}^{N_{S}}}\alpha_{I}(\omega)\rangle_{\text{Snapshots}}$. The Fourier transformation was carried out including an artificial peak broadening (see supplementary information), and the sum rule for the polarizability was normalized to the experimental value in the range 0-25 eV.

In Figure 1, we show the comparison between our calculation and the most recent experiment for the real ($\epsilon_{1}$) and imaginary ($\epsilon_{2}$) parts of the frequency-dependent dielectric constant. Although differences are evident, our subsystem TDDFT calculation recovers the overall trend and also reproduces some interesting features. The computed $\epsilon_{2}$ is consistently red shifted, however the overall multimodal shape is reproduced. For $\epsilon_{1}$, we satisfactorily reproduce the $\omega=0$ limit, as well as the peak maximum.

Subsystem TDDFT allows us to inspect the dynamical interactions between subsystems, and allows us to glimpse into the cooperative interactions that arise when the liquid responds to an external perturbation. A formally exact decomposition of the interacting electronic response function of a system into subsystem contribution Casida and Wesolowski (2004); Pavanello (2013); Neugebauer (2007) can be implemented. Namely,

[TABLE]

with each interacting subsystem response function given by a local interacting one-body term ($\chi_{I}^{u}$) and a nonlocal many-body term,

[TABLE]

where hereafter superscripts $c$ and $u$ stand for “coupled” (i.e., employing the full many-body response) and “uncoupled” (i.e., employing only the one-body response).

The subsystem TDDFT kernel for $I\neq J$ is given by $K_{IJ}({\mathbf{r}}_{1},{\mathbf{r}}_{2},t-t^{\prime})=\frac{\delta(t-t^{\prime})}{|{\mathbf{r}}_{1}-{\mathbf{r}}_{2}|}+\frac{\delta^{2}E_{\rm xc}}{\delta\rho({\mathbf{r}}_{1},t)\delta\rho({\mathbf{r}}_{2},t^{\prime})}+\frac{\delta^{2}T_{\rm s}}{\delta\rho({\mathbf{r}}_{1},t)\delta\rho({\mathbf{r}}_{2},t^{\prime})}$ Casida and Wesolowski (2004); Neugebauer (2010). The dipole polarizability is related to the response function by $\alpha_{xy}(\omega)=-\int x^{\prime}\chi({\mathbf{r}^{\prime}},{\mathbf{r}},\omega)y{\mathrm{d}}{\mathbf{r}}{\mathrm{d}}{\mathbf{r}^{\prime}}$. The uncoupled subsystem-local response function is computed with response equations that include occupied-virtual orbital excitations of only one subsystem. Or alternatively, in the real-time approach the embedding potential is computed at each time step only updating the subsystem time-dependent density and leaving the density of the surrounding subsystems frozen at $t=0$. The many-body terms include couplings between subsystem excitations at the level of two and higher bodies Krishtal and Pavanello (2016), $\chi_{I}^{u}K_{IJ}\chi_{J}^{u}+\chi_{I}^{u}K_{IJ}\chi_{J}^{u}K_{JK}\chi_{K}^{u}\ldots$

It is evident from Figure 2 that the many-body contributions are cooperative at low frequencies, and anticooperative at high frequencies. Although it is difficult to pinpoint the specific reasons for this behavior, we note that increased spectral weights at low frequencies are typical of excitonically coupled systems Onida et al. (2002). We should remark that the many-body term in Eq. (2) must be associated with an overall zero sum rule Krishtal and Pavanello (2016). Thus, if there is a many-body enhancement in one region of the spectrum there must be a many-body depletion in another region of the spectrum. In addition, the predicted many-body enhancement at low frequencies is consistent with the finding that the refractive index in water increases with increasing pressure, ascribed to an increase of the oscillator strengths with pressure Pan et al. (2014). From our results, we estimate the refractive index as $n=\sqrt{\epsilon_{1}(\omega=0)}$, and obtain $n=1.25$ and $1.30$ for the one-body and the full many-body results. Thus, many-body effects increase $n$ by about 4%.

After having characterized the role of many-body dynamical interactions between water molecules in the liquid, we now analyze the first absorption band. Subsystem TDDFT places the average peak position of the first absorption band at $\omega_{1}^{u}=6.4$ eV. Compared to the experimental 8.3–8.6 eV, we underestimate it by $\sim$2 eV. This is expected Garbuio et al. (2006), as we employ the PBE semilocal xc functional and the adiabatic approximation. Due to the more localized electronic structure compared to semilocal TDDFT of the supersystem, it is also expected that subsystem TDDFT yields a larger gap (semilocal TDDFT of the supersystem finds a gap of $\sim 5$ eV Tavernelli (2006)). For the same reasons, the subsystem TDDFT value for $\epsilon_{1}(\omega=0)=1.68$, slightly underestimates the DFT-PBE value of 1.72 Lu et al. (2008).

We find that the many-body terms have a negligible effect on the peak position of the first absorption band, which only shifts by about 0.02 eV when they are included. In other words, $\omega_{1}^{c}-\omega_{1}^{u}=-0.02$ eV. This points to a unique feature of liquid water, that is the many-body terms affect the magnitude of the dielectric constant but do not change the position of the first peak in the optical spectrum.

In Figure 3, we summarize our results pertaining to the first absorption band. Inset (a) of the figure displays histograms of the first excitation energy on the $x$-axis and the subsystem count on the $y$-axis, for the liquid (in blue) and for the vapor (dotted line). These histograms can be interpreted as vibrationally modulated joint density of states (JDOS). The gas-phase data is generated with water molecules having the same geometry as the liquid, but their electronic structure is computed in the absence of environment with the ADF te Velde et al. (2001) software employing the PBE xc functional and a quadruple-$\zeta$ Slater-Type Orbital basis set. From the JDOS, we evince that the gas-phase first excitation energy is sharply centered around the average (there are 512 data points in the histogram) at $\omega^{iso}_{1}=6.0$ eV ($\sigma^{iso}_{1}=0.07$ eV, onset at 5.8 eV), while the subsystem TDDFT places the band maximum at $\omega^{u}_{1}=6.4$ eV featuring a very broad range of excitation energies ($\sigma^{u}_{1}=0.4$ eV, onset at 5.4 eV). Such a broadening is entirely due to interactions of the water molecules with their environment. Because the environment is heterogenous and dynamic in the liquid, a broad distribution of excitation energies emerges. Additionally, Figure 3a shows that the onset of the absorption band (Urbach tail) is red shifted compared to the isolated water molecule case by 0.4 eV, again due to vibronic coupling to environmental degrees of freedom.

To justify the above claim, we provide in Figure 3b and 3c, the scatter plot of the computed excitation energy for each of the isolated and embedded (liquid phase) water molecules against the symmetric stretching degree of freedom. We find that the excitation energy for the isolated water molecules almost perfectly anticorrelates to the stretching mode. Conversely, the first excitation of the embedded water molecules does not correlate at all with this internal degree of freedom. Inspection of Figure 3b and 3c unequivocally determines that the symmetric stretching, which was determinant for the spectrum of gas phase water, no longer plays a role in the spectrum of the liquid. The other local degrees of freedom (asymmetric stretch and the bond angle) do not correlate to the absorption peak of either gas or liquid phases.

Thus, we set out to find the degree of freedom describing the environment surrounding a water molecule in the liquid that correlates the most to the first absorption peak of each subsystem. Among the ubiquitous order parameters (such as number of donated/accepted hydrogen bonds, coordination number, and tetrahedral order parameter DiStasio et al. (2014)), we found that only the number of accepted hydrogen bonds features a nonzero correlation to the first peak position. However, the computed correlation value ($\text{correlation}=\frac{\sigma(X,Y)}{\sqrt{\sigma(X,X)\sigma(Y,Y)}}$, with $\sigma$ being the variance and $X=\text{excitation energy}$, $Y=\text{degree of freedom}$) is only 0.3, and thus it is not large enough to confidently attribute the presence of a correlation. In Figure 3d we show the scatter plot of the first excitation energy of each subsystem against a composite order parameter which we call “Environment Order Parameter” (EOP) that includes degrees of freedom of molecules in the environment. EOP was constructed from 33 independent descriptors (O–O distances with the closest 6 waters and 12 respective angles with the hydrogen atoms, the O1–O2–H2 and O2–O1–H1 angles, $\theta$, the tetrahedral order parameter DiStasio et al. (2014), accepted and donated number of hydrogen bonds). These initial descriptors were combined in nonlinear ways (for example, slicing the O–O distance in this way: $\text{O--O distance}\times(\theta-\theta_{0})$, with $\theta_{0}$ being a constant angle value) to generate 268 linearly dependent environmental degrees of freedom. These were reduced by singular value decomposition of their covariance matrix (i.e., a principal component analysis) to only 29. This is a massive reduction in dimensionality, especially in view of the fact that the total number of possible binary nonlinear combinations of order parameters with 3 slices per set (i.e., 3 values of $\theta_{0}$) with 63 molecules in the environment equals to $63\times 3\times 62\times 3=35154$ parameters. The 29 linearly independent order parameters resulting from the principal component analysis where linearly combined to yield EOP in such a way to maximize the correlation to the first excitation energy of each subsystem. The maximization was carried out with a multivariable optimization based on the BFGS algorithm. Correlation of the first absorption band peak position to EOP is 0.6. This is double the original best value from simple order parameters.

Analysis of the major components of EOP highlights four order parameters. In order of contribution: the distances to the 4-th (negative correlation), 6-th (negative correlation), 1-st (positive correlation) and 2-nd (negative correlation) oxygen atoms in the environment, all multiplied by $\theta-60^{o}$ (i.e., the O1–O2–H2 angle). This explains the not large enough correlation to the number of hydrogen bonds, because the hydrogen bond definition is too simplistic and only combines a distance cut-off criterion in conjunction with an angle cut-off therefore missing the complexity of the interaction. In the supplementary information, we provide further details about the procedure to obtain EOP.

In conclusion, we carried out real-time subsystem TDDFT simulations of liquid water and compared the results to the gas-phase water. We show that, although many-body excitonic effects have no effect on the position of the first absorption band peak, they dramatically enhance the spectral intensities in the low frequency range, and deplete the intensities at high frequencies. The simulations also predict the liquid’s Urbach tail to be 0.4 eV red shifted compared to the gas phase, due to the strong coupling of the first absorption band with the environment within the first solvation shell. Interaction with the environment thus dramatically broadens the joint density of states and blue shifts the band peak by 0.4 eV. The results reproduce semiquantitatively the experimental red shift of the Urbach tail and the blue shift of the band peak.

We thank Robert A. DiStasio, Jr. for helpful discussions, particularly for pointing out the need to slice the O–O distances to enhance correlations. We also thank Alessandro Genova for helping us with the dipole Fourier transform script. This material is based upon work supported by the U.S. Department of Energy, Office of Basic Energy Sciences, under Award Number [award number].