The origin of the Debye relaxation in liquid water and fitting the high frequency excess response

TL;DR

This study investigates the microscopic origin of Debye relaxation in liquid water through simulations and spectral fitting, challenging local process models and supporting a defect migration theory.

Contribution

The paper introduces a new spectral fitting package, analyzes the high frequency excess response, and tests the Debye relaxation mechanism against experimental data and theoretical models.

Findings

Debye relaxation involves processes at 1.5-2.0 nm scales.

Debye relaxation is highly dispersive with respect to wavevector.

Traditional multi-Debye models are inconsistent with the gLST constraint.

Abstract

We critically review the literature on the Debye absorption peak of liquid water and the excess response found on the high frequency side of the Debye peak. We find a lack of agreement on the microscopic phenomena underlying both of these features. To better understand the molecular origin of Debye peak we ran large scale molecular dynamics simulations and performed several different distance-dependent decompositions of the low frequency dielectric spectra, finding that it involves processes that take place on scales of 1.5-2.0 nm. We also calculated the $k$-dependence of the Debye relaxation, finding it to be highly dispersive. These findings are inconsistent with models that relate Debye relaxation to local processes such as the rotation/translation of molecules after H-bond breaking. We introduce the spectrumfitter Python package for fitting dielectric spectra and analyze different ways of fitting the high frequency excess, such as including one or two additional Debye peaks. We propose using the generalized Lydanne-Sachs-Teller (gLST) equation as a way of testing the physicality of model dielectric functions. Our attempts at fitting the experimental spectrum using the gLST relation as a constraint indicate that the traditional way of fitting the excess response with secondary and tertiary Debye relaxations is problematic. All of our work is consistent with the recent theory of Popov et al. (2016) that Debye relaxation is due to the migration of Bjerrum-like defects in the hydrogen bond network. Under this theory, the mechanism of Debye relaxation in liquid water is similar to the mechanism in ice, but the heterogeneity and power-law dynamics of the H-bond network in water results in excess response on the high frequency side of the peak.

Tables

Table 1. Table 1: Reported two-Debye and three-Debye fits for experimental data taken at 298 K (25 C) over the last 30 years. DRS = microwave dielectric relaxation spectroscopy, ATR = THz attenuated total reflectance spectroscopy, TDS = THz time domain reflection spectroscopy, fLS - femotosecond laser spectroscopy, dFTS = dispersive Fourier Transform Spectroscopy ∗ HN model for τ 2 subscript 𝜏 2 \tau_{2} , α = 1 , β = .77 formulae-sequence 𝛼 1 𝛽 .77 \alpha=1,\beta=.77   † † {}^{\ \dagger} HN model for τ 3 subscript 𝜏 3 \tau_{3} , α = .9 , β = .8 formulae-sequence 𝛼 .9 𝛽 .8 \alpha=.9,\beta=.8

${\tau }_D$(ps)	${\tau }_2$(ps)	${\tau }_3\left(ps\right)$	$f_1$	$f_2$	$f_3$	range (cm-1)	method	ref
8.3	1.0		72	1.69(3)		.03 - 3	DRS	Barthel, 1990Barthel et al. (1990)
8.4	1.1		72	1.75		.006 - 14	DRS	Buchner, 1997Buchner et al. (1999)
8.4	0.91		72(1)	1.77(6)		.075 -10	DRS	Peacock, 2009Peacock (2009)
8.3	0.39		75	1.67(3)		.2 - 4	DRS	Sato, 2008Sato and Buchner (2004)
8.3	0.36		72	2.12		6 - 83	ATR	Möller,2009 Moller et al. (2009)
	0.248(8)		75(1)	1.67(3)			ATR	Yada, 2008Yadaa et al. (2008)
7.0(3)	0.92(6)		70(1)	2.0(3)		2 - 66	TDS	Ronne, 1997Ronne et al. (1997)
8.3	0.42		73	2		.001 - 3	TDS	Fukasawa, 2005Fukasawa et al. (2005)
8.24(4)	0.18(14)		73	1.9(5)		2-50	fLS	Kindt, 1996Kindt and Schmuttenmaer (1996)
8.8(6)	0.21(6)		73	1.5(8)		3-55	fLS	Venables, 1998Venables and Schmuttenmaer (1998)
7.8	0.2		73	1.6		.16-33	var	Liebe, 1991 Liebe et al. (1991)
8.21	0.39(5)		73	2.5(2)		.1 - 33	var	Benduci, 2007Beneduci (2008)
8.31	1.0∗	0.10†	71.5	2.8	1.6	50-220	dFTS	Vij, et al. 2004Vij (2004)
8.26(3)	1.1(5)	0.14(4)	73	2.2(2)	1.3(3)	.1 - 33	var	Benduci, 2007Beneduci (2008)
8.5	0.93	0.08				.03-800	var	Ellison, 2007Ellison (2007)
8.4(3)	1.05(15)	0.18(5)				0.02-37	var	Vinh, et al., 2015Vinh et al. (2015)

Table 2. Table 2: Some of the H-bond network modes in liquid water. All of these modes are IR active and should appear in the dielectric response.

approx freq. (cm-1)	description	ref
50-65	H-bond bending (in plane)	Agmon (1996)
70	H-bond torsion	Agmon (1996)
150	H-bond sym. stretch (“breathing”)	Agmon (1996)
180	H-bond asym. stretch	Agmon (1996)
80-150	assymetric umbrella mode	Heyden et al. (2010)

Equations

$$\frac{\epsilon (\omega )-{\epsilon }_{\infty }}{\epsilon (0)-{\epsilon }_{\infty }}=\frac{A_{\text{mbox}HN}}{[1+(i\omega {\tau }_{\text{mbox}HN}{)}^{\alpha }{]}^{\beta }}$$

$$\epsilon (\omega )=\frac{\epsilon (0)-{\epsilon }_{\infty }}{1+i\omega {\tau }_D}+{\epsilon }_{\infty }$$

$${\tau }_{\text{mbox}DS}(T)=\frac{4\pi \eta (T)R^3}{k_{B}T}$$

$$\frac{1}{D(T)}=\frac{6\pi \eta (T)R}{k_{B}T}\propto {\tau }_{\text{mbox}DS}(T)$$

$$\frac{{\tau }_D}{{\tau }_s}=\frac{\epsilon (0)+2}{{\epsilon }_{\infty }+2}\approx 11$$

$$\frac{{\tau }_D}{{\tau }_s}=\frac{3\epsilon (0)}{2\epsilon (0)+{\epsilon }_{\infty }}=1.46$$

$$\frac{{\tau }_D}{{\tau }_s}=\frac{3\epsilon (0)G_K}{2\epsilon (0)+{\epsilon }_{\infty }}$$

$$\frac{{\tau }_D}{{\tau }_s}=\frac{2\epsilon (0)+{\epsilon }_{\infty }}{\epsilon (0)+2{\epsilon }_{\infty }}=1.8$$

$$\tau (T)$$

$$\varphi (t)=⟨\sum _i{\mathbf{\mu }}_i(0)\cdot \sum _j{\mathbf{\mu }}_j(t)⟩$$

$${\varphi }^{\text{mbox}ds}(t,R)=⟨\sum _i{\mathbf{\mu }}_i(0)\cdot \sum _{j\in R_i}{\mathbf{\mu }}_j(t)⟩$$

$${\varphi }^{\text{mbox}ss}(t,R)=\sum _i⟨{\mathbf{\mu }}_i^s(0)\cdot {\mathbf{\mu }}_i^s(t)⟩$$

$${\mathbf{\mu }}_s(t)={\mathcal{N}}_i(t)\sum _{j\in R_i}{\mathbf{\mu }}_j(t)$$

$$\epsilon (\omega )$$

$$f(\omega )$$

$$\epsilon (\omega )=\frac{f}{1+{\left[(i\omega {\tau }_{\text{mbox}defect}{)}^{-1}+(i\omega {\tau }_{\text{mbox}osc}{)}^{-\delta }\right]}^{-1}}+{\epsilon }_{\infty }$$

$$\epsilon (\omega )=\sum _i^N\frac{f_i}{1+i\omega {\tau }_{Di}}+\sum _j^M\frac{f_j{\omega }_{Tj}^2}{{\omega }_{Tj}^2-{\omega }^2-i\omega {\gamma }_j}+{\epsilon }_{\infty }$$

$$\epsilon (0)-{\epsilon }_{\infty }=\sum _i{f}_i$$

$${\epsilon }^{\prime }(\omega )$$

$${\epsilon }^{\prime \prime }(\omega )$$

$$\text{mbox}Cost=\sum _i{\left(\frac{\text{mbox}{fit}_i-\text{mbox}{data}_i}{\text{mbox}{data}_i}\right)}^2$$

$$\frac{({\omega }_L^2+{\gamma }_L^2)}{{\omega }_T^2}=\frac{\epsilon (0)}{{\epsilon }_{\infty }}$$

$$\sum _i^N\frac{{\tau }_{Ti}}{{\tau }_{Li}}\sum _j^M\frac{({\omega }_{Lj}^2+{\gamma }_{Lj}^2)}{{\omega }_{Tj}^2}=\frac{\epsilon (0)}{{\epsilon }_{\infty }}$$

$$\epsilon (\omega )={\epsilon }_{\infty }\prod _i^N\frac{(\omega -i{\omega }_{Li})}{(\omega -i{\omega }_{Ti})}\prod _j^M\frac{{\omega }_{Lj}^2-{\omega }^2-i\omega {\gamma }_j}{{\omega }_{Tj}^2-{\omega }^2-i\omega {\gamma }_j}$$

$$\frac{1}{\epsilon (\omega )}=\frac{1}{{\epsilon }_{\infty }}\prod _i^N\frac{(\omega -i{\omega }_{Ti})}{(\omega -i{\omega }_{Li})}\prod _j^M\frac{{\omega }_{Tj}^2-{\omega }^2-i\omega {\gamma }_j}{{\omega }_{Lj}^2-{\omega }^2-i\omega {\gamma }_j}$$

$$\epsilon (\omega )=\frac{1}{\sqrt{2\pi }\sigma }{\int }_{-\infty }^{\infty }{e}^{-\frac{(x-{\omega }_0{)}^2}{2{\sigma }^2}}\frac{{\omega }_p^2}{x^2-{\omega }^2+i\omega \gamma }+{\epsilon }_{\infty }$$

Full text

I Introduction

The Debye relaxation peak dominates the dielectric absorption spectrum of water. The peak is centered at $\approx$ 20 GHz (0.66 cm*-1*) and spans six decades of frequency. The large oscillator strength of the Debye peak (73 at 25 C*∘*) can be viewed as the main contributor to water’s anomalously high static dielectric constant through the $f$-sum rule. The Debye absorption peak is of immense practical importance as it is used in microwave ovens and in satellite-based microwave radar sensing of water and ice.Rosenkranz (2015); Wang and Schmugge (1980); Kneifel et al. (2014)

Water is an anomalous, complex liquid that serves as the arena for all life on planet Earth. As such, there is a long legacy of research on water which continues to be built upon today. Current topics of interest among researchers include the possible liquid-liquid phase transition in water,Nilsson and Pettersson (2016); Limmer and Chandler (2011) understanding water’s behaviour under extreme thermodynamic conditions and confinement, and understanding the important role of nuclear quantum effects in water.Ceriotti et al. (2016); PhysRevB.92.134105 In this paper we focus on the dielectric relaxation of water, about which an enormous literature already exists, including several recent detailed theoretical studies.Popov et al. (2016); Ben Ishai et al. (2015); Hansen et al. (2016); Arbe et al. (2016); Lunkenheimer et al. (2016) In our review of this literature, we found a lack of agreement on the molecular origins of the Debye relaxation and further disagreement on how to fit excess response on the high frequency side of the Debye relaxation, which historically had been fit with a secondary Debye relaxation. Some authors relate Debye relaxation to particular translational and/or rotational motions of a single molecules after breaking one or more hydrogen bonds.Ronne et al. (1997); Vij (2004); Agmon (1996) Other authors propose that Debye relaxation is due to the movement of “free” molecules which only have one or two hydrogen bonds.Buchner et al. (1999); Nabokov and Lubimov (1988); Stanley and Teixeira (1980); J.B. Hasted (1985); Zasetsky and Buchner (2011); von Hippel (1988) In contrast to these theories, other authors describe Debye relaxation as a collective relaxation of a cluster or large collection of molecules,Arkhipov (2002); Arkhipov and Agmon (2003) which is more in line with dielectric theories that establish how the collective relaxation time increases due to dipole-dipole correlation. Finally, recently a few authors have also proposed the hopping of defects allow the rearrangement of the hydrogen bond network and Debye relaxation,Popov et al. (2016); Ben Ishai et al. (2015); Artemov and Volkov (2014) in analogy to the process underlying Debye relaxation in ice.Hobbs (2010) As a macroscopic measurement, dielectric relaxation spectroscopy alone does not provide enough information to distinguish between the various mechanisms that have been proposed. In this work, we analyze large scale molecular dynamics simulations to gain insight into the collective nature of Debye relaxation.

It is worth dwelling on the fact that the Debye peak in water is nearly perfectly Debye. Most other dipolar liquids and virtually all polymeric liquids exhibit more complex relaxation.Kindt and Schmuttenmaer (1996) In such liquids, a phenomenological equation called the Havriliak-Negami (HN) formula is often used:

[TABLE]

The case $\alpha\neq 1,\beta=1$ is known as Cole-Cole relaxation, and corresponds to a symmetric distribution of exponential relaxations. The applicability of the HN equation to water has been tested several times before. For instance Kaatze (1993) found $\alpha=.989(2)$ and $\beta=.959(4)$.Kaatze (1993) Vij et al. (2004) report $\alpha=1$, $\beta=1$.Vij (2004) Mason’s careful fit of the Cole-Cole equation yields $\alpha=0.988\pm.008$.Mason et al. (1974)

The fact that the Debye relaxation peak is so well described by a single relaxation process clashes with our understanding of the structure and dynamics of the hydrogen bond network. Saito and Ohme Saito and Ohmine (1994) used molecular dynamics simulation to find that the relaxation of the polarization vector of a large water cluster exhibits a long-time tail, which can be fit with either a stretched exponential (also called a Kohlraush function, $P(t)\propto\exp[-(t/\tau)^{\beta}]$) or a $1/f^{\alpha}$ power law, both of which correspond to a broad distribution of relaxation times.Baba et al. (1997) Similarly, molecular dynamics studies show that both the relaxation of the polarization vector of single molecules and hydrogen bond auto-correlation functions are well described by a stretched exponential.Luzar and Chandler (1996); Kumar et al. (2006) Ohmine and Tanaka, in a detailed review, present evidence from molecular dynamics simulations that hydrogen-bond network rearrangement dynamics are complex and highly collective in nature.Ohmine and Tanaka (1993) They find evidence that the hydrogen bond network “contains many relaxation processes, with many time scales”. This is reflected in the fact that the low frequency Raman spectrum of water between $1-20$ cm*-1* is very diffuse and can be fit with a power law.Ohmine and Tanaka (1993) Reconciling the complex heterogeneous dynamics of the H-bond network with the pure exponential character of the Deybe relaxation has been recognized as an important unresolved issue.Shiratani and Sasai (1996); Yamaguchi et al. (2003); Ohmine and Tanaka (1993)A possible clue to solving this issue comes from the fact that if Coulomb interactions are smoothly truncated at 9 Å, dielectric relaxation decreases from $\approx 9$ ps to only 1 ps and assumes a $1/f^{\alpha}$ character.Ohmine and Tanaka (1993) The structure of water in such situations may be highly non-physical as well.Yonetani (2005); van der Spoel and van Maaren (2006) These results suggest that long range dipole-dipole interactions and/or the long range structure of the H-bond network are necessary to recover the exponential character of Debye relaxation. Further evidence comes from studies with salt, which show that Debye relaxation does not change very much with increasing salt concentration.Vinh et al. (2015); Turton et al. (2008)

II Critical analysis of previous ideas about Debye relaxation

II.1 Models based off Debye-Stokes theory

We first consider Debye’s original theory from 1929.Debye (1929) Debye considers a thermal ensemble of non-interacting molecules in an applied electric field, and considers what happens when the field is turned off. His starting assumption, which is now known to be incorrect, is that each individual water molecule undergoes Brownian rotational motion. Molecular dynamics simulations show that molecular relaxation in water actually occurs in highly discontinuous “jumps” due to the breaking of H-bonds, rather than small angle Brownian diffusion.Luzar and Chandler (1996); Laage and Hynes (2008, 2006); Ludwig (2007); Geiger (2003) Debye then takes the diffusion equation in spherical coordinates, linearizes it, and finds that the average moment decays exponentially with time, leading directly to the Debye equation for dielectric relaxation:

[TABLE]

Under this model, the Debye relaxation time $\tau_{D}$ is related to the rotational friction constant $\zeta$ via $\tau_{D}=\zeta/2k_{B}T$. Stokes showed that for a sphere of radius $a$ rotating in a medium with shear viscosity $\eta$, the rotational friction is given by $\zeta=8\pi\eta R^{3}$. This leads to the Debye-Stokes model for $\tau_{D}(T)$:

[TABLE]

This equation fits the experimental data for $\eta(T)$ and $\tau_{D}(T)$ remarkably well, with a value of $R=1.44\mbox{\AA}$ at 0 C∘ which is about the right radius for a single water molecule.Ronne et al. (1997); Agmon (1996) Many authors have noted this agreement and concluded that Debye’s model is essentially correct, and that Debye relaxation is due to the rotational relaxation of single molecules. Furthermore, Stoke’s model leads to the Stokes-Einstein relation, which says that the translational diffusion constant is given by:

[TABLE]

The Stokes-Einstein relation is borne out experimentally, as is the proportionality between $\tau_{D}(T)$ and $1/D(T)$ – both follow Arrhenius-like temperature dependencies with very similar rate coefficients.Nabokov and Lubimov (1988) Equation 4 only breaks down when water comes supercooled.Bertolini et al. (1982); Agmon (1996); Arkhipov (2002) Bertolini argues that the Arrhenius temperature dependence must be due to a barrier hoping process.Bertolini et al. (1982) Agmon builds his theory for Debye relaxation in water on equations 3 and 4, suggesting that it is due to translational hopping. In particular, he proposes that Debye relaxation is due to a hopping process called “tetrahedral displacement”.Agmon (1996) Tetrahedral displacement has a hopping distance of 3.3 Å, which is “the separation between an occupied and unoccupied corners of a cube binding the pentawater tetrahedron”.Agmon (1996) We will argue later that Agmon’s model is either incorrect or incomplete, because it does not explain the collective nature Debye relaxation that we find.

Hansen et al. note that both the Debye-Stokes and Stokes-Einstein equations fail for ordinary liquids, yielding a value $R_{H}$ which is much smaller than the molecular $R$.Hansen et al. (2016) They note that for most molecular liquids, empirically it is found that $R_{H}\approx R/2$.Hansen et al. (2016) If one uses this empirical relation, one finds $\tau_{\mbox{\footnotesize{DS}}}\approx\tau_{D}/8\approx\tau_{2}$. In other words, the agreement of $R_{H}\approx R$ when equation 3 is applied to water’s Debye relaxation should be taken as coincidental. This makes sense if one understands $\tau_{D}$ as a collective phenomena, and $\tau_{2}$ as being related to single molecule relaxation.

II.2 Mean-field theories

The Debye-Stokes model is based on the flawed assumption that interactions between molecules can be ignored - notably dipole-dipole forces and H-bonding interactions. The effects of dipole-dipole interaction can be approximately accounted for by using a mean-field theory. The simplest mean-field theory is obtained by taking Debye’s molecular dipole undergoing Brownian diffusion and moving it into a homogeneous medium described by a frequency dependent dielectric function $\varepsilon(\omega)$. The dipole then feels an additional “internal field” due to the polarization response of the medium. After solving the system self-consistently, one again obtains the Debye equations for $\varepsilon(\omega)$, but now:Arkhipov (2002)

[TABLE]

In other words, the relaxation time for the polarization of the entire system, as measured through $\varepsilon(\omega)$, is greater than the relaxation time for the single dipole $\tau_{S}$. Here $\varepsilon_{\infty}$ can be understood as the excess oscillator strength not described by Debye relaxation. We use the experimental values $\varepsilon(0)=78.6$ and $\varepsilon_{\infty}=5.4$ to obtain a ratio of $10.8$. The ratio experimentally varies considerably with experiment (table I) between 8.3 - 34, with an average of $\approx 13$. We know that the Debye model is wrong, though, because the same mean field theory framework gives a completely wrong estimate of $\varepsilon(0)$, and predicts that water exists in a ferroelectric phase at room temperature.

The Onsager mean-field model for $\varepsilon(0)$,Onsager (1936) consisting of a dipole in a cavity, is considered a significant improvement over Debye’s model. When the dipole moment is in a cavity, a “reaction field” field appears in addition to the internal field. Cole extended the Onsager mean-field model to the time dependent case, yielding a complicated dielectric function which corresponds to a non-exponential relaxation. Glarum extended Onsager’s model to the time-dependent case,Glarum (1960) recovering the Debye equations with:

[TABLE]

A similar theory by Powles yields:Powles (1953)

[TABLE]

Here the Kirkwood factor $G_{K}$ includes the effects of local dipole-dipole correlation.

[TABLE]

A variety of more sophisticated mean field theories have been developed, the details of which we will not recount here.Bagchi (2012); Fatuzzo and Mason (1967); Nee and Zwanzig (1970); Arkhipov (2002); Bagchi and Chandra (1990); Arkhipov (2002) A common outcome of these models is that dipole-dipole correlations lead to the macroscopic dipole relaxation time of being longer than the molecular one. The same principle holds for clusters relaxing in a dielectric environment.

II.3 Propagation of defects

Finally we come to the recent idea that Debye relaxation is due to the propagation of defects, as is the case in ice. There are four main defects in ice – two charged ($H_{3}O^{+}$ & $OH^{-}$) and two uncharged (Bjerrum L & D defects). The propagation of these defects are responsible for Debye relaxation in ice, with the defect with the lowest energy barrier determining the timescale of the relaxation.Petrenko and Whitworth (1999)

Popov et al. show that in liquid water the activation energies of charged and uncharged defects should be roughly the same, so distinguishing which may underlie Debye relaxation requires some additional analysis. Artemov & Volkov propose that Debye relaxation is entirely due to the conduction of the charged defects.Artemov and Volkov (2014) Their model makes the claim that nearly 10% of water molecules are ionized, which differs by six orders of magnitude from the accepted value of $K_{w}=10^{-7}$.Artemov and Volkov (2014); Volkov et al. (2014) This idea is inconsistent with classical molecular dynamics simulations, which can satisfactorily reproduce the dielectric response of water, even though they do not contain charged defects.Elton and Fernández-Serra (2014); Sega and Schröder (2015) Furthermore, Popov et al. note that the dielectric relaxation does not depend on pH, as it would in Artemov & Volkov’s model.Popov et al. (2016) Popov et al. propose Debye relaxation is due entirely to Bjerrum-like defects, which carry an effective charge. The bifuricated hydrogen defect (Bjerrium D-like) results in excess positive charge locally and the bifuricated Oxygen defect (Bjerrium L-like) results in excess negative charge locally. If defects follow ordinary diffusive behaviour ($\langle r^{2}(t)\rangle=6D_{\mbox{\footnotesize{defect}}}t$), then it is easy to show that the Debye equation for the dielectric response results.Popov et al. (2016) The existence of defects Bjerrium L & D defects (biffuricated bonds) in liquid water is supported the analysis of x-ray scattering data by Sciortino et al.Sciortino et al. (1990a) It is also supported by the molecular dynamics simulations of Laage & Hynes which found bifuricated bonds lead to jump relaxation.Laage and Hynes (2006) We find the theory of Popov et al. attractive as previously we showed that the librational and OH-stretching dynamics of water are very similar to that of ice and that both originate from propagating phonon-like modes which travel through the hydrogen bond network.Elton and Fernández-Serra (2016); Elton (2016)

III The secondary Debye process and high frequency excess

We also seek to address an ongoing controversy on how to fit the high frequency side of the Debye relaxation, between 1 - 100 cm*-1*, where there is unaccounted for excess response. Ishai et al. and Popov et al. propose that excess is also due to the dynamics of defects propagating through the hydrogen bond network which was discussed in the previous section.Ben Ishai et al. (2015); Popov et al. (2016) Traditionally however, this excess response has been fit by introducing a second Debye mode, characterized by a time constant $\tau_{2}$. However, while the value of $\tau_{D}$ is very consistent among experiments, values for $\tau_{2}$ vary considerably in the experimental literature, as shown in table I. An additional problem, noted by Beneduci, is that the data on the temperature dependence of $\tau_{2}$ is contradictory - Barthel et al. find it increasing with temperature while Ronne, et al. find it decreasing.Beneduci (2008) To fix these issues, recently it was proposed that a 3rd Debye relaxation ($\tau_{3}$) is also required to properly fit the high frequency excess.Beneduci (2008); Vinh et al. (2015)

As with the primary Debye mode, different authors have different hypotheses about the microscopic mechanisms that underlie the secondary and putative tertiary Debye modes. $\tau_{2}$ is usually associated with the rotational relaxation of single molecules and hydrogen bond breaking. Some attribute it to the relaxation of weakly bound molecules.Buchner et al. (1999) Molecular dynamics show the average hydrogen bond lifetime to be around 0.5 - 1.0 ps at 300 K,Ohmine (1995) with a very broad distribution.Sciortino et al. (1990b) Others propose that the excess response is either the $\alpha$ or (more oftenly) $\beta$ relaxation found in supercooled liquids.Lunkenheimer et al. (2016)

In light of the heterogeneous and non-exponential dynamics of the hydrogen bond network, the use of one or two additional exponentials to fit the excess response seems ad-hoc. The ad-hoc nature of this fit can also be seen by considering some of the infrared active H-bond vibrational modes that have been shown to exist in the region of $10-200$ cm*-1*, as shown in table II. Heyden et al. have shown that all of these modes span a broad frequency range between $20-200$ cm*-1* due to the inhomogeneous nature of the H-bond network.Heyden et al. (2010) Instantaneous normal mode analysis of water shows a broad spectrum of translational modes, extending from very 1-400 cm*-1* and peaked around 100 cm*-1*.Cho et al. (1994); Ohmine (1995); Kindt and Schmuttenmaer (1996)

Thus it appears that a large number of Debye and resonance processes contribute to the excess response. This can be modeled by fitting a distribution of modes. Unfortunately, fitting a distribution of modes is a mathematically ill-posed problem - many distributions and combinations of modes may be consistent to the experimental data within the uncertainty of the data. The wide variation in experimental fits is at least partially explained by the fact that the approximation of one or two additional Debye relaxations to fit the excess is crude and that each experiment only measures a certain window of frequencies. Since the fit is approximate, it varies depending on the particular range of the experiment. Furthermore, some authors do not unbias their fitting, so the fit can be biased toward the side of the spectrum that is higher in magnitude. Dielectric relaxation spectroscopy (DRS) experiments can only probe the low frequency part (.0001 - 2 cm*-1*), while THz or microwave time-domain reflection spectroscopy (TDS) probe the“middle” frequencies (1 - 10 cm*-1*), as does a variation known as attenuated total reflectance spectroscopy (ATR).Moller et al. (2009) Finally, Fourier transform infrared spectroscopy (FTIR) covers the region above 10 cm*-1*.

A side point is that an excess response on the high frequency side of Debye relaxation is a general feature found in many dipolar liquids, including non H-bonding liquids. It seems largely forgotten that this was first pointed out in 1955, when Poley noticed that an excess around 10 cm*-1* (0.3 THz) exists in many dipolar liquids.Poley (1955) To explain this phenomena, which was called “Poley absorption” at the time, Hill and others proposed that it was due to inertial motion.Hill (1971); Nora E. Hill (1963) Physically, inertial motion can be pictured as either fast “rattling” of molecules within sharply defined potential energy wells or as nearly-free rotations over small angles.Bagchi and Chandra (1993); Ferraro and Basile (1978) In water, inertial absorption (if relevant) would overlap with hydrogen bond network vibrations and modes. Inertial relaxation results in approximately a Gaussian form for $\phi(t)$ near $t=0$.Hill (1971) In our previous workElton and Fernández-Serra (2016) we did not find any evidence of such Gaussian relaxation except in $\phi(k,t)$ for $k>3\mbox{\AA}^{-1}$.

IV Fitting the temperature dependence

The temperature dependence of Debye relaxation is usually cited as being Arrhenius.Agmon (1996) According to transition state theory, Arrhenius temperature dependence implies an free energy barrier $\Delta H$, and in water one finds $\Delta H$ is roughly the hydrogen bond energy. On the basis of this, Buchner proposes that the Arrhenius temperature dependence of $\tau_{D}(T)$ is due to the production rate of “free molecules”, which he defines as molecules having 1 or 0 H-bonds.Buchner et al. (1999) The percolation model of Stanley & Tiexera indicates that $\tau_{D}(T)$ is due to reorientation of molecules having only one hydrogen bond.Stanley and Teixeira (1980) Nabokov and Lubinov adapt a similar type of model, but argue that molecules with both one and two hydrogen bonds are mobile and thus contribute.Nabokov and Lubimov (1988)

Nabokov et al. note that $\tau_{D}(T)$ deviates from Arrhenius behaviour at low temperature.Nabokov and Lubimov (1988) To our knowledge a test of different fit functions over the entire experimentally accessible range has not yet been published. We combined data from Ellison (2007)Ellison (2007) from between 273 - 373 K and the data collected by Nabokov (1988)Nabokov and Lubimov (1988) from between 255 - 533 K. The data taken above 373 K was measured along the gas-liquid co-existence curve. We fit the log normalized data with the following fit functions:

[TABLE]

The results are shown in an Arrhenius plot in figure 1. Clear deviations from Arrhenius behaviour (which would be a straight line in this plot) are observed. The behaviour of $\tau_{D}(T)$ is very nearly Arrhenius between 273-313 K, but the only model which fits the data through the entire temperature range is the Vogel-Fulcher-Tammann relation, which is used to fit the “$\alpha$-relaxation” in glasses. We found $\tau_{\infty}=0.14$ ps, $D=4.52$ and $T_{\mbox{\footnotesize{VFT}}}=141$ K. For the Arrhenius model we found $\Delta H=0.16$ eV (3.7 kcal/mol), which is approximately the hydrogen bond energy. The mode-coupling theory power law was previously shown to fit $\tau_{D}(T)$Kumar et al. (2006) in simulations of SPC/E water.Starr et al. (1999) The Speedy equation was included because it is used to fit other response functions for water, especially in the supercooled region.Speedy and Angell (1976)

VFT temperature dependence is a universal feature of both relaxor ferroelectrics and dipolar glasses,Pirc and Blinc (2007); Bokov et al. (1999) and indicates the presence of spatial heterogeneity.Vilgis (1990); Giovambattista et al. (2003) A phenomenological theory by Tagantsev shows that near $T_{\mbox{\footnotesize{VFT}}}$ VFT temperature dependence is a consequence of a very wide distribution of relaxation times in the system.Tagantsev (1994) A very general theory for the VFT equation is the Adam-Gibbs model, which assumes the existence of cooperatively rearranging regions (CRRs), which are clusters that relax independently from each other, similar to the polar nanoregions concept we introduced in a previous work.Giovambattista et al. (2003)

V Molecular dynamics simulations

We performed a number of molecular dynamics simulations using the TIP4P/2005 water modelAbascal and Vega (2005) and the GROMACS 4.5.5 molecular dynamics package.Abascal and Vega (2005) Our simulations used a Nosé-Hoover thermostat with $\tau=0.5$ ps and a timestep of 2 fs. Long range Coulomb interactions were handled with the particle mesh Ewald method. All simulations were equilibrated for at least 100 ps before trajectory output. Our 10,000 molecule (box size $L=6.68$ nm) simulation was 4 ns long and used a Coulomb cutoff of $r=3.34$ nm. We analysed the simulation trajectories using our epskw Fortran code for molecular liquids. The code is open source and available online at www.github.com/delton137/epskw.

VI $k$-dependence

Figure 2 shows the $k$ dependence of the Debye relaxation, calculated using the method we describe in our previous work.Elton (2016); Elton and Fernández-Serra (2016) By fitting the underlying $k$ dependent correlation functions we produced a plot of $\tau_{D}(k)$ which is shown in figure 3. At larger spatial scales (smaller $k=2\pi/\lambda$) the relaxation is slower than at shorter scales (larger $k$).

Recently, Arbe et al. have measured the dynamic structure factor of water, $S(Q,\omega)$ using incoherent and coherent neutron scattering in the GHz-THz frequency range and the intermediate $Q$ range (0.3 - 2.0 $\mbox{\AA}^{-1}$).Arbe et al. (2016) They observe a structural relaxation process at same frequency as the Debye peak in the dielectric spectra. Assuming the two processes correspond, our result (fig. 3) is qualitatively consistent with their measurement of $\tau_{D}(k)$ they derived from fitting a double exponential process to $F(Q,t)$, the intermediate incoherent scattering function for $H$ nuclei. Arbe et al. conclude that the $k$ dependence indicates the presence of translational diffusion, but is not consistent with an interpretation in terms of a single diffusive motion, as in Agmon’s model. They label the Debye mode as the “diffusive” model, and the region of excess response is fit with a “local” Debye relaxation. Arbe et al. also found that the “local” Debye relaxation does not exhibit dispersion with $k$, (hence the name “local”).Arbe et al. (2016)

VI.1 Distance decomposition of the Debye relaxation

There are several ways to perform distance-dependent decomposition of the Debye relaxation. The most intuitive way is to break the simulation cell into sub-boxes of different sizes and compute the total dipole moment for each sub-box. Dipole time correlation functions are then computed separately for each sub-box and averaged. The process is repeated for sub-boxes of different sizes (fig. 4). We find that the dipole relaxation time does not converge to the bulk value until the box size is increased to $\approx 2.0$ nm.

We investigated two other methods of distance-decomposition. In the first method, which we call the “dip-sphere” method, one starts with the dipole-dipole time-correlation function:

[TABLE]

Then, one limits the molecules around each molecule $i$ to those in a sphere of radius $R$:

[TABLE]

The other method, which we call the “sphere-sphere” method, was introduced by Heyden, et al.Heyden et al. (2010) The sphere-sphere method is so-called because one calculates the autocorrelation of the total dipole moment of a sphere of radius $R$ centered around a reference molecule, and then averages this over each molecule in the system:

[TABLE]

where:

[TABLE]

Heyden et al. recommend the normalization factor $\mathcal{N}_{i}(t)=(1+N_{ij}^{2})^{-1/2}$ to normalize for number of molecules in each sphere. This normalization factor is chosen so that in the bulk limit ($R\rightarrow\infty$) the original full response function is obtained. Heyden et al. also propose introducing a smoothing function to weight molecules within the sphere, which we neglect here for simplicity.

The results of these two types of decomposition are shown in figure 5. The dielectric relaxation time increases up to about 1.5 nm. Interestingly, this is the maximum length scale we found feasible for clusters in liquid water – beyond 1.5 nm water-water interactions behave as in a dielectric continuum.Elton and Fernández-Serra (2014) A sphere of 1.5 nm contains approximately 470 water molecules. Our results indicate that complexes of 100s of molecules participate in Debye relaxation. The two models that are most consistent with this finding are the propagation of defects model and the model of relaxing clusters.

VII Fitting the high frequency excess

As mentioned, there is uncertainty and confusion on how to fit and interpret the excess response on the high frequency side of the Debye relaxation. There are many pitfalls when fitting a dielectric spectra, especially when there there are many overlapping modes. When the underlying relaxation is described by several overlapping modes or a distribution of relaxations and/or damped harmonic oscillators, the number of parameters becomes cumbersome or impossible for many fitting algorithms. More importantly, the fitting becomes fundamentally under-defined – many fits become consistent with the data, within the noise of the data. Trying to fit a distribution of relaxation times $G(\tau)$ to a spectrum is an ill-defined problem, and attempts to do so can yield peculiar results.Zasetsky and Buchner (2011) Still, this type of fitting can be done using Tichonov regularization, which enforces the distribution to be smooth and prevents overfitting.Schäfer et al. (1996); Dias (1996); Weese (1992); Tuncer and Macdonald (2006); Ktitorov (2003)

Distinguishing the validity of various proposed dielectric functions is difficult. Sheppard and Grant found that “.. data represented by a small departure from single relaxation time behaviour ($0.9<\alpha<1$ in eqn. 1) could be equally well interpreted as being due to two distinct kinds of relaxation process with relaxation times separated by a factor as high as three.”Sheppard and Grant (1974) A similar point is made by Barker, who notes that a single Debye mode is nearly indistinguishable from two closely overlapping modes.Barker (1975)

Puzenko, et al. introduced the following dielectric function,Puzenko et al. (2005)

[TABLE]

This function was recently used by Ishasi et al. to fit the high frequency excess in water.Ben Ishai et al. (2015) A cutoff frequency $\omega_{c}$ is required, since otherwise this expression violates the Kramers-Kronig relations. In our experiments with this fit function, we found this cutoff could be ignored. Putting the cutoff on at the center of the Debye peak results in an unphysical discontinuity in the fit function. As a power law, this fit corresponds to a flat distribution of relaxation times in the hydrogen bond network. The interpretation of $A$ and $q$ in this fit function are not clear.

Under the defect-migration theory for Debye relaxation introduced by Popov, et al., the Debye peak assumes the following form:Popov et al. (2016)

[TABLE]

This fit takes into account the high frequency wing, without the need for the addition of one or two extra Debye relaxation processes. Hidden in this dielectric function is a power law of the form $\omega^{-\delta}$, where $0<\delta<1$.

A custom Python package called spectrumfitter was developed, which is available at github.com/delton137/spectrumfitter or from the Python Package Index. We started with the assumption that dielectric function has $N$ Debye relaxation processes and $M$ damped harmonic oscillator processes, and thus has the following form:

[TABLE]

Fitting is performed with the “$f$-sum rule” enforced as a constraint:

[TABLE]

where $f_{i}$ is the oscillator strength of the $i$th mode. We fit the experimental refractive index data compiled by Segelstein (1981), which comes from a compilation of all experimental sources of dielectric function and index data that were available at the time.Segelstein (1981) While dated, the Segelstein dataset has the benefit that it covers the entire frequency range from $0.001-200,000$ cm*-1* The dielectric permittivity is obtained from the complex index of refraction ($n,k$) data using:

[TABLE]

We first did an interpolation of the Segelstein data on a logarithmic grid up to 18 cm*-1* and a linear grid from 18-4000 cm*-1*, to prevent biasing the fitting towards the Debye mode. We performed the fitting in an unbiased manor by minimizing the sum of the relative errors squared:

[TABLE]

VII.1 Using the gLST relation as a novel constraint

Any dielectric function should obey the Kramers-Kronig relations, which is derived from the basic principle of causality (causes must proceed effects). The $f$-sum rule (eqn. 17) comes from taking $\omega=0$ in the Kramers-Kronig relations. Barker shows how a generalized Lydanne-Sachs-Teller (LST) relation can also be derived from the Kramers-Kronig relations.Barker (1975) For a single damped harmonic oscillator mode, the gLST relation is:

[TABLE]

The dampening factor appears in the numerator since the longitudinal frequency is complex $\bar{\omega}_{L}=\omega_{L}+i\gamma_{L}$. The generalized LST relation for $N$ Debye modes and $M$ damped harmonic oscillator modes reads:Barker (1975)

[TABLE]

Barker’s exposition on the gLST equation suggests it can help distinguish the validity of various fit functions, for instance, whether one should fit with a single Debye relaxation or two overlapping Debye relaxations.Barker (1975) However, to do this, one needs to work with the entire spectrum, so that all the modes that contribute to $\varepsilon(0)$ can be accounted for.

We attempt to do this by first rearranging the dielectric function given in eqn. 16 into the following form:

[TABLE]

The longitudinal frequencies correspond to points where $\varepsilon(\omega)=0$. Solving for the longitudinal frequencies requires factoring a messy $(N+2M)$th degree polynomial, which can be done numerically. An easier alternative method is to simply invert eqn. 22

[TABLE]

By symmetry it should be easy to see that the dielectric function for $\frac{1}{\varepsilon(\omega)}$ has the same form as for $\varepsilon(\omega)$, but with the transverse frequencies changed to the longitudinal frequencies. Thus, we can do a separate fit of $\frac{1}{\varepsilon(\omega)}$ using the same dielectric function given in eqn. 16 (with the change of $\varepsilon_{\infty}\rightarrow\frac{1}{\varepsilon_{\infty}}$) to obtain the longitudinal frequencies.

This method of fitting allows us to try fitting with the gLST equation as a novel constraint. Alternatively, one could use the Kramers-Kronig relations directly as a constraint, but this requires one has the entire spectrum available. Sometimes the optimization process is not able to satisfy the gLST constraint. Since the spectrumfitter code shows the contribution to the left hand side of gLST equation (eqn. 21) from each lineshape, one can often pinpoint features in the spectrum are problematic to meeting the constraint.

VII.2 Results of fitting $f$-sum and gLST constraints

We performed fitting with the $f$-sum rule and gLST equation as constraints. The inclusion of the gLST equation as a constraint did not significantly improve the quality of fit, but it did yield some insights into the physicality of the fit functions being used. Figure 7 shows a fit with 3 Debye relaxations, 1 DHO peak for H-bond stretching and 3 DHO peaks for the librational region. Under this model, the LHS of the gLST relation here is 284 while the RHS is 48.15. The 2nd Debye relaxation contributes a factor of 7.11 and the 3rd Debye relaxation contributes a factor of 4.2 to the gLST relation. Next we gave the fitting fitting procedure the option of suppressing the 3rd Debye relaxation. It was completely suppressed in the transverse case ($f\rightarrow 0$), suggesting that it is not physical. We found that if we tried to include an additional DHO lineshape between 10 and 100 cm*-1* for the H-bond bending peak it was also completely suppressed by the fitting procedure.

We found a better fit can be obtained using the “Brendel” lineshape instead of a damped harmonic oscillator for the librational and stretching peaks. The Brendel lineshape is a convolution of the DHO lineshape with Gaussian:Brendel and Bormann (1992)

[TABLE]

This lineshape can be calculated analytically in terms of complex error functions, making it easy to evaluate numerically.Brendel and Bormann (1992) We found the Brendel lineshape especially helpful for fitting the H-bond stretching peak at $\approx 150$ cm*-1*. Figure 6 shows a fit with 3 Debye relaxations, 2 Brendel peaks for H-bond stretching and H-bond bending, and 3 Brendel peaks for the librational region. To perform the gLST analysis with the Brendel peaks, we consider the Brendel lineshape can be approximately factored into the gLST equation as a single DHO. Then, the LHS of the gLST relation is 456 while the RHS is 48.15. The 1st Debye relaxation contributes a factor of 1.55 to the gLST relation, while the 2nd Debye relaxation contributes a factor of 5.49 and the 3rd Debye relaxation contributes a factor of 6.48.

In every type of fit we tried, inclusion of the 3rd Debye relaxation resulted in a significantly larger departure from the the gLST equation. If only 2 Debye relaxations are used in this model, the gLST LHS becomes 21.921, while the RMS error increases slightly from 0.120 to 0.125.

VIII Conclusion

We critically reviewed the literature on the Debye relaxation and the high frequency excess response and found widespread disagreement on the molecular mechanisms underlying these processes. We studied the distance decomposition of the Debye relaxation, finding that it involves significant correlations between molecules on scales of 1.5 - 2 nm. The $k$-dependence of Debye relaxation suggests that Debye relaxation is a collective process, and the dispersion relation of the Debye peak suggests a propagating mode. Our findings call into question models of Debye relaxation such as Agmon’s model,Agmon (1996) which related Debye relaxation to rotation/translation after a hydrogen bond breaks, and Buchner’s model, which related it to the rotation of nearly free molecules.Buchner et al. (1999)

Using our spectrumfitter Python package we attempted fitting experimental spectra with secondary and tertiary Debye processes to model the high frequency excess. When we tried to apply the gLST relation to our fitting we found that both the fits with one secondary and an additional tertiary Debye process are problematic. When the gLST relation is used as a constraint with such fits, large violations must be incurred, especially when the tertiary Debye process is included. Given the understanding of H-bond network dynamics that has been elucidated by the work of many authors during the last few decadesSciortino et al. (1990a); Ohmine (1995); Heyden et al. (2010); Kindt and Schmuttenmaer (1996) as well as in our recent work,Elton and Fernández-Serra (2016) we believe there is a wide distribution of H-bond network modes which contribute in this region. In light of this, the lack of agreement of values for $\tau_{2}$ and $\tau_{3}$ reported in the previous literature (table I) is not surprising. The narrow frequency windows of most experiments and failure to unbias when fitting further contribute to the dispersion in values.

The totality of our findings are consistent with the recent model of Debye relaxation by Popov et al.,Popov et al. (2016) which posits that it is caused by the movement of defects in the hydrogen bond network. These defects are similar to the L and D Bjerrum defects found in ice, but likely more smeared out.Popov et al. (2016) The existence of such defects in liquid water, under the name of “bifuricated” H-bonds has already been established from x-ray scatteringSciortino et al. (1990a) and simulation.Laage and Hynes (2006) The migration of these defects through the H-bond network explains the semi-long range dependence we found ($\approx$ 15 Å) as well as the dispersion of the Debye peak we found from simulation and that was found in X-ray scattering experiments.Arbe et al. (2016) The migration of these defects is not purely diffusive, but occurs due to hopping motions. The idea is attractive in light of our recent work showing ice-like phonons that propagate through the H-bond network.Elton and Fernández-Serra (2016) If the dynamics of water are ice-like at the librational frequencies $400-800$ cm*-1*, then it is not surprising they are also ice-like at lower frequencies. Dielectric relaxation in ice is almost perfectly described by a single Debye relaxation, at least at high temperatures.Hobbs (2010)

According to the theory of Popov et al., the hopping motion of a defect results in both translational and oscillatory motions of water molecules, which cause vibrations in the H-bond network. Under this assumption, and taking into account the power-law vibrational dynamics of the H-bond network, Popov et al. are able to derive the excess wing of the Debye relaxation.Popov et al. (2016) Thus we also see that the theory of Popov et al. helps resolve the mystery of how a single Debye peak arises from the heterogeneous (and fractal-like) H-bond network. The heterogeneous and fractal-like nature of the H-bond network does in fact appear in the spectrum, giving rise to the excess response. In other words, the Debye peak in liquid water is not a standard pure Debye peak, but is a Debye peak whose high frequency side is modified by the presence of the complex H-bond network.

IX Acknowledgements

This work was supported by DOE Award No. DE-FG02-09ER16052. The author would like to thank his Ph.D. adviser, Prof. Marivi Fernández-Serra for allowing him to pursue this research and providing guidance and critical proofreading of the manuscript. The author also acknowledges Dr. Kevin Hauser, Dr. Ivan Popov, Prof. Marcello Sega, and Prof. John Watterson for providing useful feedback.