# Terahertz Dynamics in the Glycerol-Water System

Johanna Kölbel  Department of Chemical Engineering, University of Cambridge,
Cambridge, CB3 0AS, UK  Walter Schirmacher  Department of Physics, Johannes
Gutenberg University, 55099 Mainz, Germany  Center for Life Nano Science
@Sapienza, Istituto Italiano di Tecnologia, 295 Viale Regina Elena, I-00161,
Roma, Italy  Evgenyi Shalaev  AbbVie, Irvine, CA, USA  J. Axel Zeitler
<EMAIL_ADDRESS>Department of Chemical Engineering, University of Cambridge,
Cambridge, CB3 0AS, UK

###### Abstract

The model glass-former glycerol and its aqueous mixtures were investigated
with terahertz-time domain spectroscopy (THz-TDS) in the frequency range of
$0.3\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}3.0\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$
at temperatures from $80\text{\,}\mathrm{K}305\text{\,}\mathrm{K}$. It was
shown that the infrared absorption coefficient measured with THz-TDS can be
theoretically related to the reduced Raman intensity
($\propto\alpha/\omega^{2}$) and the reduced density of states
($\propto\alpha/\omega^{3}$) and the agreement with experimental results
confirms this. The data were further used to investigate the behaviour of
model glasses in the harmonic (below the glass transition temperature
$T_{\text{g}}$), anharmonic (above $T_{\text{g}}$), and liquid regime. The
onset temperature of the molecular mobility as measured by the infrared active
dipoles, $T_{\text{g}}$, was found to correlate with the onset of anharmonic
effects, leading to an apparent shift of the boson peak and obscuring it at
elevated temperatures. The influence of clustered and unclustered water on the
dynamics, the boson peak, and the vibrational dynamics was also investigated.
A change in structural dynamics was observed at a water concentration of
approximately $5\text{\,}\,\mathrm{w}\mathrm{t}\,\%$, corresponding to a
transition from isolated water molecules distributed homogeneously throughout
the sample to the presence of small water clusters and an increased number of
water-water hydrogen bonds which lower the barriers on the potential energy
surface.

##  I Introduction

###  I.1 Dynamics of Glass-Formers

Upon cooling an organic molecular liquid to temperatures below its melting
point, it either forms a crystal or remains disordered, ultimately forming an
amorphous solid. Initially, at temperatures just below its freezing
temperature a liquid is referred to as supercooled. In such thermodynamically
metastable systems the Stokes-Einstein relation breaks down, and diffusion and
viscosity decouple [1]. The molecular arrangement remains similar to that of
liquids and long-range order is completely absent. Upon further cooling to
temperatures below the glass transition temperature ($T_{\text{g}}$), a sharp
change in the thermal expansion coefficient and specific heat is observed and
a glass is formed.

Motions that occur at temperatures at and above $T_{\text{g}}$ include
translational, cooperative molecular motions. At temperatures below
$T_{\text{g}}$ the molecular motions are restricted to local or so-called
secondary mobility of atoms rattling in the cages formed by their neighbours
[2]. Experimentally, the $T_{\text{g}}$ can be measured by a range of
techniques, such as differential scanning calorimetry (DSC), which detects the
change in heat capacity at the glass transition temperature.

Terahertz time-domain spectroscopy (THz-TDS) probes frequencies in the range
of $0.35\text{\,}\mathrm{THz}3\text{\,}\mathrm{THz}$ by coupling to dipole
moments at photon energies on the order of the hydrogen bond strength and
picosecond relaxation time. It hence provides a sensitive probe of mobility in
condensed phase organic molecular materials. Past studies using THz-TDS have
investigated the complex interplay between reorientation motion and vibration
dynamics in hydrogen-bonded liquids [3, 4, 5].

Due to the inherent disorder of glasses, no discrete phonon modes can be
sustained and the coupling of photons to the vibrational density of states
(VDOS) results in the relatively featureless so-called microscopical peak that
constitutes the main absorption mechanism in the terahertz range [6] with a
frequency-squared dependence as predicted by Debye theory [7]. With increasing
temperature, an overall increase in absorption is measured within the
accessible spectral bandwidth of THz-TDS instruments
($<$3\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$$). In addition to the
microscopical peak, THz-TDS also detects the excess density of states (DOS)
above the Debye level. This excess DOS is termed boson peak (BP).

While the exact origin of the BP is still not fully understood, it is
generally accepted that the BP is a harmonic phenomenon due to inherent
disorder, and should therefore be temperature independent [8, 9, 10].

Goldstein proposed that properties of glasses can be explained by the
characteristics of the potential energy surface (PES) [11]. The atomic
dynamics have been described as a consequence of the shape of the PES [12]. In
this context, at low temperatures, the system is trapped in a deep and steep
minimum which can be approximated as harmonic. Anharmonic effects play only a
small role at the low temperatures. The measurable increase in absorption
coefficient at single frequencies, however, is due to anharmonic effects [13].
At temperatures at which the apparent centre frequency of the BP is not
affected by anharmonicity, we will call the temperature range the harmonic
regime.

Once the temperature is increased sufficiently for the system to escape the
deep minimum, it has access to shallower regions of the PES that are
characterised by numerous local minima separated by smaller energy barriers.
The temperature at which the system contains sufficient thermal energy to
leave the deep minimum on the PES coincides with $T^{*}$, i.e. the onset of
local mobility, observed by THz-TDS as an increase in the rate of absorption
change with temperature. The beta-relaxation is the result of hopping
processes across low energy barriers. Anharmonic effects begin to dominate the
spectra above a certain temperature when non-harmonic parts of the PES minima
can be explored, resulting in an apparent change of the BP frequency.

A further temperature increase, ultimately resulting in $T_{\text{g}}$, opens
up more local minima on the PES while anharmonic effects continue to play an
important role. The larger the amplitude of the possible motions, the higher
the associated change in dipole moment and hence the measured absorption
coefficient. The shallower the minima, the more anharmonicity influences the
spectra.

Glycerol is a widely studied network material and model system for a glass
former [14] as it is non toxic and easily supercools and has been investigated
with experimental techniques like Raman scattering [15, 16], infrared
spectroscopy [17, 18], dielectric spectroscopy [19, 20], neutron and light
scattering [14, 21], DSC [22], terahertz spectroscopy [23], as well as
theoretical work with molecular dynamics simulations [17, 24]. Its
$T_{\text{g}}$ was found to lie between $185\text{\,}\mathrm{K}$ and
$194\text{\,}\mathrm{K}$ [22, 25, 26, 27]. In glycerol, the VDOS extends from
approximately $0.5\text{\,}\mathrm{THz}7\text{\,}\mathrm{THz}$ [25].

THz-TDS measurements of glycerol revealed that its $T_{\text{g,}\beta}$ is
more pronounced than the $T_{\text{g}}$[27, 28]. Glycerol is hence a promising
model system to investigate whether anharmonic effects can be detected at
temperatures even below $T_{\text{g}}$. Here we study the impact of
temperature changes on glycerol dynamics, especially anharmonic effects that
can obscure the BP.

The density of states of pure glycerol, including the elevated temperature
range, where the spectrum becomes temperature dependent, is well documented:
Chumakov  _et al._ [29] embedded a probe molecule in a glass matrix and
monitored its translational motions with nuclear inelastic scattering. Because
the probe followed the collective motions of the glass with a correlation
length larger than the probe size, the density of states of collective motions
of the glass matrix, i.e. the VDOS multiplied with the Debye-Waller factor,
was directly measured. It was reported that a significant part of the BP in
the model glass former glycerol is constituted of collective modes and that it
disappears close to $T_{\text{g}}$ due to increased sample mobility. They
further observed an exponential decrease in the reduced VDOS at energies above
the BP maximum [29].

Wuttke et al. found that the reduced density of states of glycerol is
temperature-independent at frequencies above the BP maximum [30]. At
frequencies below the maximum, a temperature-dependent change was observed
even at temperatures below $T_{\text{g}}$. These temperature-dependent changes
in the BP occur above the characteristic temperature $T^{*}$, which coincides
with the onset of localised motions and therefore increased anharmonicity. The
onset of localised motions at $T^{*}$ is typically too subtle to be detected
in DSC measurements [31]. Ruggiero et al. described $T^{*}$ as the
vitrification of the $\beta$-relaxation process as originally detailed by
Johari and Goldstein, also called $T_{\text{g,}\beta}$. [12]

###  I.2 Glycerol-Water Mixtures

Modern biopharmaceutical drugs are often developed into complex amorphous
formulations that are prepared by spray-drying or freeze-drying
(lyophilization). Given its resistance to crystallisation, glycerol is used
widely as a so-called cryoprotectant of cells and organs and also to stabilise
lyophilised proteins [32, 33]. Some pharmaceutical degradation mechanisms are
catalysed by water clusters and the influence of water content on the
stability of lyophilised products has been investigated previously [34]. In
macromolecules, water can both serve to stabilise the native structure while
simultaneously acting as a catalyst for destabilisation. Starciuc et al.
suggested that unclustered water molecules (that are commonly found in
products containing only a few wt. % water) are less catalytically active than
water clusters [35]. They hypothesised that proton transfer becomes possible
only once a water clustering threshold is exceeded, thereby supporting
pharmaceutical degradation reactions such as amide hydrolysis and deamidation.

The model system of glycerol-water mixtures has been studied widely [36, 17,
19, 37, 38, 39]: Towey et al. used a combination of neutron diffraction
experiments and computational modelling to investigate the structure and
hydrogen-bonding of pure glycerol [38], dilute aqueous glycerol solution [37],
and glycerol-water mixtures for glycerol mole fractions of
$x_{\text{g}}=0.05,0.10,0.25,0.50,0.80,\text{ and }1.00$ [39]. They found
bipercolating clusters of both water and glycerol in samples containing a mole
fraction between 0.75 and 0.5 water (i.e. between
$16.4\text{\,}\,\mathrm{w}\mathrm{t}\,\%37\text{\,}\,\mathrm{w}\mathrm{t}\,\%$
water). They further found that water maintained its full hydrogen-bonding
capacity independently of the glycerol concentration.

Murata and Tanaka found evidence for a water liquid-liquid transition (LLT)
without macroscopic phase separation in glycerol-water mixtures containing
between approximately
$45.5\text{\,}\,\mathrm{w}\mathrm{t}\,\%56\text{\,}\,\mathrm{w}\mathrm{t}\,\%$
water [36].

Recently, Starciuc et al. investigated glycerol-water mixtures containing
between
$0\text{\,}\,\mathrm{w}\mathrm{t}\,\%40\text{\,}\,\mathrm{w}\mathrm{t}\,\%$
water with low- and high-wavenumber Raman spectroscopy, and three water-
content regions were found [35]. This confirmed the existence of a threshold
for water clustering in glycerol. A transition from predominantly unclustered
water molecules to small water clusters was reported at $6.1\pm 0.7$ wt.%water
and a second threshold at $18.6\pm 4.4$ wt.%water, which was proposed to
correspond to the formation of large water clusters associated with the onset
of freezing.

And while glycerol and glycerol-water systems have been investigated
thoroughly with a range of methods, there has not yet been a comprehensive
study in the terahertz frequency range over a range of temperatures to
investigate the role of water, especially the formation of water clusters, in
glycerol-water mixtures with a water content below
$30\text{\,}\,\mathrm{w}\mathrm{t}\,\%\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}$,
where past studies have shown the existence of several water regions. Typical
specification limits for freeze-dried pharmaceuticals are below
$5\text{\,}\,\mathrm{w}\mathrm{t}\,\%$. By understanding the structural
dynamics of the glycerol-water system, such water content specifications in
the design of pharmaceutical products can be scientifically substantiated and
their design can be optimised.

The dynamics of the model system of glycerol-water mixtures are investigated
with THz-TDS by varying both the temperature and the water concentration. The
aim was to extract as much information as possible about both the glass and
the amorphous solid and to investigate the role of the BP and the shape of the
VDOS and their relationship with the glass transition temperatures.

In addition to the THz-TDS data at various (cryogenic) temperatures, high-
quality room temperature data of liquid glycerol-water mixtures with different
water concentrations were also acquired to further understand how the water
content influences the terahertz spectra of the liquids.

##  II Materials and Methods

###  II.1 Variable-Temperature THz-TDS Measurements

Glycerol was obtained from Sigma Aldrich (Poole, UK) and was mixed with
Milli-Q water (IQ 7000, Merck, Darmstadt, Germany, resistivity
$18.2\text{\,}\mathrm{M\SIUnitSymbolOhm}\text{\,}\mathrm{cm}$) in various
ratios
($0\text{\,}\,\mathrm{w}\mathrm{t}\,\%\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}30\text{\,}\,\mathrm{w}\mathrm{t}\,\%\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}$).
After mixing, the samples were sealed and left standing or briefly de-gassed
inside a desiccator attached to a vacuum pump. No experimental differences
between samples prepared by the two methods were observed. For subsequent
measurements, a drop of liquid sample was sealed between two z-cut quartz
windows (Crystran, UK) separated by a $100\text{\,}\mathrm{\SIUnitSymbolMicro
m}$ thick PTFE spacer in a liquid cell which was then fit into the cryostat
and placed under vacuum. The samples were analysed in transmission using a
commercial Terapulse 4000 spectrometer (TeraView Ltd., Cambridge, UK) with a
liquid nitrogen cryostat attached, as described by Sibik et al. [40]. Two
z-cut quartz windows (Crystran, UK) were used for reference measurements. The
spectral range the instrument was able to access was in the range of
$0.35\text{\,}\mathrm{THz}3\text{\,}\mathrm{THz}$, depending on sample
absorption.

Utilising the cryostat and liquid nitrogen cooling, a wide range of
temperatures were studied ($80\text{\,}\mathrm{K}305\text{\,}\mathrm{K}$). At
each temperature, a co-average of 1000 waveforms was acquired for both
reference and sample and used to calculate the optical constants.

To increase reproducibility and facilitate comparison between different water
concentrations, the measurements were controlled with a computer program that
kept the length of temperature steps and therefore the overall heating rate
constant at about $1\text{\,}\mathrm{K}\text{\,}{\mathrm{min}}^{-1}$.

###  II.2 THz-TDS Measurements at Room Temperature

Glycerol and its mixtures were prepared following the same steps as described
above. However, instead of sealing the sample into a liquid cell that fit into
the cryostat, a room temperature liquid cell was utilised which was easier to
assemble. This liquid cell also consisted of two z-cut quartz windows
separated by a $100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thick PTFE spacer.
This was subsequently placed into the measurement chamber which was purged
with dry nitrogen and samples were analysed in transmission using the
Terapulse 4000 spectrometer utilising the same settings as for measurements at
variable temperatures.

###  II.3 Relationship between reduced density of states measured with THz-
TDS, neutron, Raman, and nuclear inelastic scattering

In this section we develop theoretical ideas relating the infrared absorption
coefficient $\alpha(\omega)$ to the mechanical degrees of freedom. The
absorption coefficient is related to the dielectric loss function by
$\alpha(\omega)=\omega\epsilon^{\prime\prime}(\omega)/\left(n(\omega)c\right)$.
The loss function, in turn, can be related to the complex shear modulus
$G(\omega)=G^{\prime}(\omega)-iG^{\prime\prime}(\omega)$ by the formula of
Gemant [41, 42, 43], as utilised by Schirmacher, Ruocco and Mazzone [44].

$\epsilon^{\prime\prime}=\mathfrak{Im}\bigg{\\{}\frac{\textstyle 1}{\textstyle
1+\frac{VG}{kT}}\bigg{\\}}\propto\mathfrak{Im}\bigg{\\{}\frac{\textstyle
1}{\textstyle G}\bigg{\\}}$ (1)

$G(\omega)$ is the complex mechanical shear modulus, the frequency dependence
of which is mainly given by the sound attenuation of the shear waves
$\Gamma(\omega)$ [45]:

$G(\omega)=G^{\prime}(\omega)\left[1-i\Gamma(\omega)/\omega\right]$ (2)

Therefore $\alpha(\omega)$ is given by

$\alpha(\omega)\propto\omega\epsilon^{\prime\prime}(\omega)\propto\omega^{2}\frac{\Gamma(\omega)}{1+\Gamma(\omega)^{2}/\omega^{2}}$
(3)

This means that we can extract the transverse sound attenuation coefficient
from the infrared absorption, which is usually not accessible to experimental
investigation.

A very similar formula can be obtained for the Raman intensity
$I_{\text{AB}}(\omega)=\left[n(\omega)+1\right]\chi^{\prime\prime}_{\text{AB}}(\omega)$,
where A and B denote V or H polarisation and
$n(\omega)=(\exp{(\hbar\omega/(\text{k}T))}-1)^{-1}$. The Raman
susceptibilities can be written as [10]:

$\chi^{\prime\prime}_{\text{VV}}(\omega)=\frac{4}{3}\chi^{\prime\prime}_{\text{VH}}(\omega)+Af_{1}\chi^{\xi}_{\text{L}}(\omega)$
(4)

$\chi^{\prime\prime}_{\text{VH}}(\omega)=Af_{2}\frac{1}{30}\left(2\chi^{\xi}_{\text{L}}(\omega)+3\chi^{\xi}_{\text{T}}(\omega)\right)$
(5)

where

$\displaystyle\chi^{\xi}_{\text{L,T}}(\omega)$ $\displaystyle=$
$\displaystyle\mathfrak{Im}\bigg{\\{}\frac{3}{k^{3}_{\xi}}\int_{0}^{k_{\xi}}k^{2}dk\frac{\textstyle
k^{2}}{\textstyle-\omega^{2}+k^{2}c^{2}_{\text{L,T}}(\omega)}\bigg{\\}}$ (6)
$\displaystyle=$
$\displaystyle\mathfrak{Im}\bigg{\\{}\frac{1}{c^{2}_{\text{L,T}}(\omega)}+O(\omega^{2})\bigg{\\}}$

These equations have been derived under the assumption of fluctuating Pockels
constants [46, 47]. The upper cut-off $k_{\xi}$ is proportional to the inverse
correlation length $\xi$ of the Pockels constant fluctuations.
$c_{\text{L}}(\omega)$ and $c_{\text{T}}(\omega)$ are the complex longitudinal
and transverse sound velocities. $f_{1}$ and $f_{2}$ are the pre-factors of
the longitudinal and transverse Pockels correlation functions, respectively.
In terms of the longitudinal modulus $M(\omega)$ and the shear modulus
$G(\omega)$ the complex sound velocities are given by

$\rho c^{2}_{\text{L}}(\omega)=M(\omega)\quad\quad\rho
c^{2}_{\text{T}}(\omega)=G(\omega)$ (7)

where $\rho$ is the mass density.

Because the inverse of the transverse sound velocity is much larger than the
longitudinal one, the Raman scattering intensity will be dominated by the
transverse susceptibility. For low enough frequencies we essentially have:

$I_{\text{AB}}(\omega)\propto\left[n(\omega)+1\right]\mathfrak{Im}\bigg{\\{}\frac{\textstyle
1}{\textstyle G}\bigg{\\}}.$ (8)

So we expect the reduced Raman intensity to be proportional to the infrared
absorption coefficient divided by $\omega^{2}$:

$\tilde{I}(\omega)\equiv\frac{I(\omega)}{\omega(n(\omega)+1)}\propto\alpha(\omega)/\omega^{2}.$
(9)

The reduced Raman intensity is also commonly represented as: [47]

$\tilde{I}(\omega)\propto C(\omega)\frac{g(\omega)}{\omega^{2}}$ (10)

with the phenomenological frequency dependent coupling coefficient $C(\omega)$
and the reduced DOS $\frac{g(\omega)}{\omega^{2}}$. In most materials
$C(\omega)\propto\omega$ [48].

Thus, the reduced DOS $\frac{g(\omega)}{\omega^{2}}$ is related to the
infrared absorption coefficient by

$\frac{g(\omega)}{\omega^{2}}\propto\frac{\alpha(\omega)}{\omega^{3}}$ (11)

##  III Results and Discussion

###  III.1 Pure Glycerol Measurements and Analysis

####  III.1.1 The Boson Peak

The terahertz absorption spectra measured by THz-TDS of the investigated
glycerol-water samples are featureless. An example spectrum is shown in Figure
1a. Without further analysis it is not possible to distinguish between the
glass and the liquid states.

Figure 1: Comparison of terahertz to neutron, nuclear inelastic scattering,
and Raman data. (a) THz-TDS absorption spectrum of glycerol for temperature
steps of $10\text{\,}\mathrm{K}$ from
$80\text{\,}\mathrm{K}280\text{\,}\mathrm{K}$ and steps of
$5\text{\,}\mathrm{K}$ above. The sample is amorphous over the range of
temperatures measured and absorption increases with temperature at all
frequencies. (b) and (c) BP visualised of the same data by plotting
$\alpha/\nu^{2}$ (b) and $\alpha/\nu^{3}$ (c). A broad maximum is visible at
around $1.3\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$ at low temperatures that
disappears upon heating. (d) Comparison of the BP in the VDOS measured with
neutron scattering ([49], dots) and the BP measured with nuclear inelastic
scattering ([29], squares). The latter maximum is located at slightly lower
frequencies than in the neutron scattering data. (e) Comparison of the BP in
the VDOS measured with neutron scattering ([49], dots) and $\alpha/\nu^{3}$
measured with THz-TDS (lines). The maxima and shapes are very similar. (f)
Comparison of the $\alpha/\nu^{2}$ measured with THz-TDS (lines) and the VDOS
measured with Raman scattering ([50], dots). The maxima below
$260\text{\,}\mathrm{K}$ occur at similar frequencies and the shape is very
similar.

For comparison of the infrared data with neutron and Raman scattering data,
the absorption coefficient divided by frequency to the power of two and three
is plotted in Figures 1b and c, respectively. At low temperatures, the BP is
clearly visible at around $1\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$ and its
amplitude and centre frequency are mostly constant. With increasing
temperature, the maximum appears to shift to lower frequencies before features
of the BP can no longer be detected at even higher temperatures.

Chumakov et al. [29] used nuclear inelastic scattering to measure delocalised
collective motions in glycerol with a correlation length greater than
$20\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. In comparison to the measurements
by Wuttke et al. [14] the nuclear inelastic scattering data detects the BP at
slightly lower frequency, while exhibiting a similar shape overall (Figure
1d). In pure glycerol Wuttke et al. found a pronounced BP at
$170\text{\,}\mathrm{K}$ that disappeared at $210\text{\,}\mathrm{K}$. This
observation is confirmed by our THz-TDS measurements (Figure 1c). Even at very
low temperatures (below $120\text{\,}\mathrm{K}$), the BP, as measured with
THz-TDS, increases in intensity while the maximum frequency stays constant. In
the same temperature interval, the neutron data shows no temperature-dependent
change.

The BP was also visualised with low-frequency Raman scattering by Uchino and
Yoko [50]. The temperature range in their study was somewhat higher but the
peak was found situated at frequencies very close to the ones measured by THz-
TDS as shown in Figure 1f. At $170\text{\,}\mathrm{K}$, the intensity at
frequencies below the frequency of the maximum is very similar. It is notable
that while in the Raman data a distinct peak with a well-defined maximum was
observed at $260\text{\,}\mathrm{K}$, no distinct peak shape is obvious in the
terahertz data.

Light scattering data acquired by Cummins et al. [14] also show the BP located
at approximately $1.5\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$ at temperatures
above $173\text{\,}\mathrm{K}$. In their data, the BP does not disappear
completely, even at elevated temperatures in the liquid regime. The agreement
with terahertz data is good for low temperatures, as shown in Figure 2.

Figure 2: (a) $\chi^{\prime\prime}$ from Cummins’s light scattering data [14]
for various temperatures. (b) Comparison with terahertz data. Better agreement
with $\alpha/\nu^{2}$. The maximum in the light scattering data is located at
slightly higher frequencies than in the terahertz data.

While the BP itself is temperature-independent [8, 9, 10], additional
anharmonic effects that contribute to the absorption give the impression of an
apparent shift of the BP with temperature. By tracking the changes of the
apparent maximum we can observe when the anharmonic effects influence the
absolute maximum intensity. This process ultimately results in the intensity
of the BP being subsumed entirely by the anharmonic contributions to the
absorption intensity.

To track the apparent change in the peak frequency ($\nu_{\text{BP}}$) and its
intensity, the data are first smoothed with a moving average and the numerical
derivative is calculated. To determine whether the BP is still discernible at
elevated temperatures, the inflection point in the DOS is also calculated from
the derivative and is shown in Figure 3 for different temperatures.
$\nu_{\text{BP}}$ and the inflection point are plotted against temperature in
Figure 4. In pure glycerol, the decrease in the apparent centre frequency of
the BP is very sudden. For comparison, the more subtle change in some of the
glycerol/water mixtures is shown in Supplementary Figure S1.

Figure 3: The density of states and its derivative for different temperatures.
The inflection point is highlighted with a cross in the derivative. At
$295\text{\,}\mathrm{K}$, the BP is obscured by anharmonic effects and neither
an inflection point nor a maximum in the DOS can be found.  Figure 4:
$\nu_{\text{BP}}$ (circles) and inflection point (crosses) for glycerol for
different temperatures as determined from smoothed data. Above
$210\text{\,}\mathrm{K}$, no inflection points and maxima were found. The
frequency resolution (approximately
$30\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$) in the smoothed data limits the
accuracy when determining $T^{*}$ and $T_{\text{g}}$.

The behaviour of the BP is hence described by a set of three parameters: the
centre frequency at the lowest measured temperature (approximately
$0.9\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$), the temperature at which the
centre frequency is affected by anharmonic effects and appears to decrease
($T^{*}$, in glycerol approximately $170\text{\,}\mathrm{K}$), and the highest
temperature at which an inflection point is still present in the DOS ($T_{g}$,
in glycerol $210\text{\,}\mathrm{K}$, Figure 4).

####  III.1.2 The Absorption Coefficient at Cryogenic Temperatures

The temperature behaviour of glasses can be separated into 3 regimes: at low
temperatures in the harmonic regime, any mobility (e.g. dihedral angle
changes) is severely restricted and molecules vibrate at their equilibrium
position due to their thermal energy. An increase in thermal energy results in
an approximately linear increase in absorption that can be measured throughout
the accessible frequency range. This also affects the apparent intensity of
the (itself temperature-independent) BP. The onset of local mobility at
$T^{*}$ coincides with a transition region which is concluded by
$T_{\text{g}}$ and the anharmonic regime. With increasing temperature, the
VDOS is further enhanced until the liquid state (the third regime) is reached
where quasi-elastic scattering occurs and the BP is no longer defined.

The calorimetric glass transition temperature $T_{\text{g}}$ of glycerol is
hard to detect with THz-TDS, as has previously been reported by Sibik  _et
al._ [27] and Capaccioli  _et al._ [28]. Sibik et al. [27] found a mean
$T_{\text{g}}$ of $194\text{\,}\mathrm{K}$. Capaccioli  _et al._ [28] argued
that a secondary glass transition $T_{\text{g,}\beta}$ at
$161\text{\,}\mathrm{K}$ was consistent with the same $T_{\text{g}}$.

When the spectra can be described by a linear function, that function’s
gradient $a$ also yields information about the strength of anharmonic
processes.

Whether the absorption coefficient does indeed follow a linear frequency
dependence is first investigated by plotting the first derivative, as shown in
Figure 5a. This shows that the spectra can only be approximated with linear
functions below $0.5\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$ and above
$1.25\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$. In these regions, the average
gradients mirror the behaviour of the $\nu_{\text{BP}}$. Interestingly, lower
frequencies are stronger influenced by $T^{*}$ and higher frequencies by
$T_{g}$.

Figure 5: (a) First derivative of spectra. Above approximately
$1.25\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$ and below
$0.5\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$, the spectra can be described by
linear functions. The vertical line denotes $\nu_{\text{BP}}$. (b) Temperature
dependence of average gradient for frequencies between
$1.25\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}1.8\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$
(dots) and
$0.35\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}0.5\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$
(crosses), i.e. the anharmonicity parameter $a$. Shown by vertical lines are
also $T^{*}$ and $T_{\text{g}}$.

Below $T^{*}$, the gradients $a$ at low frequencies are almost constant,
implying that anharmonic effects do not dominate. Above $T^{*}$, however, $a$
increases noticeably. The gradient at higher frequencies is temperature-
dependent even below $T^{*}$. At $T_{\text{g}}$, a jump is observed, and $a$
decreases.

These observations and explanations are now used to examine the influence of
water content and temperature on the behaviour of glycerol-water mixtures.

###  III.2 The Influence of Water Concentration and Temperature on Glycerol-
Water Mixtures

When data at various concentrations are reported, one has to keep in mind that
Starciuc et al. showed that the transition from predominantly unclustered
water molecules to clusters occurs at around 6.1 $\pm$
$0.7\text{\,}\mathrm{w}\mathrm{t}\,\mathrm{\char 37\relax}\
\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}$. A second threshold was
found at 18.6 $\pm$ $4.4\text{\,}\mathrm{w}\mathrm{t}\,\mathrm{\char
37\relax}\ \mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}$ and linked to
larger, percolating water clusters [35]. At concentrations below the first
threshold, water molecules are predominantly unclustered, do not induce major
structural changes, and form mostly water-glycerol hydrogen bonds. The LLT
found by Murata and Tanaka occurs at higher water concentrations than
investigated here [36]. Any change in behaviour found at approximately 6 or
$19\text{\,}\mathrm{w}\mathrm{t}\,\mathrm{\char 37\relax}\
\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}$ can therefore be
attributed to the influence of water on the structural dynamics by increasing
the cluster size and connectivity, as well as a strengthening of the hydrogen
bonds between water molecules.

####  III.2.1 Harmonic Regime ($T\leq T_{g}$)

The BP is shown for all samples in Figure 6a at a temperature of
$90\text{\,}\mathrm{K}$. The dependence of its centre frequency and intensity
on water content are shown in Figure 6b. In most measurements, the BP is
centred at $0.95\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$. In measurements
where the BP is less pronounced, the centre frequency deviates from that
value. This seems to be more likely at water concentrations below
$5\text{\,}\,\mathrm{w}\mathrm{t}\,\%$.

Figure 6: (a) BP at $90\text{\,}\mathrm{K}$ for all samples. (b) Centre
frequency (black dots) and intensity (orange circles) at
$90\text{\,}\mathrm{K}$ for all measured concentrations.

In the harmonic regime, unclustered and clustered water molecules do not
influence the centre frequency of the BP. The Ioffe-Regel crossover is linked
to the centre frequency of the BP and denotes a crossover from wave-like to
random-matrix-like physics as the mean free path of transverse waves becomes
equal to their wavelength [29, 9, 10]. The results indicate that the Ioffe-
Regel crossover is not influenced by the presence of clustered or unclustered
water.

The onset of water clustering may decrease the intensity of the BP slightly,
indicating a higher level of disorder. While the centre frequency is largely
independent of water content, the peak intensity may exhibit a shallow minimum
around a concentration of about $5\text{\,}\,\mathrm{w}\mathrm{t}\,\%$ when
water molecules are homogeneously distributed throughout the sample.

In Figure 7a, an overview of the different transition temperatures is given
for all glycerol-water mixtures measured for transition temperatures based on
the different methods introduced when discussing pure glycerol. The
temperatures at which the BP appears to shift and dissolves are characterised
for all samples and the results are also shown in Figure 7 together with the
transition temperatures that are extracted from the anharmonicity parameter
$a$ in low and high frequency bands.

Figure 7: (a) Different transition temperatures characterising the onset of
mobility and anharmonicity and the temperature at which the BP dissolves. (b)
Mean $T^{*}$ for different water contents. (c) Rate of change of the
anharmonicity parameter with temperature in different temperature regimes.
Above $T_{\text{g}}$, a shallow minimum is observed around
$5\text{\,}\,\mathrm{w}\mathrm{t}\,\%$.

For clarity, the different $T^{*}$ have been combined into a mean $T^{*}$,
shown in Figure 7b, and it is apparent that it does not depend on water
content below approximately $5\text{\,}\,\mathrm{w}\mathrm{t}\,\%$ and
decreases at higher concentrations.

This means that unclustered water molecules do not change the temperature at
which anharmonic effects result in an apparent shift of the BP. The height of
the energy barrier separating the harmonic from the anharmonic regime is hence
not influenced by unclustered water molecules embedded in the glycerol matrix.

The shape of the deepest minimum of the PES is independent of water
concentration or the presence of water clusters, as otherwise the rate of
change of the anharmonicity parameter $a$ would be different, as shown in
Figure 7c. A change in $a$ with temperature corresponds to increased
anharmonic contributions.

At water concentrations above $5\text{\,}\,\mathrm{w}\mathrm{t}\,\%$,
parameters change their behaviour: the intensity of the BP increases slightly
with concentration, and anharmonic effects decrease the temperature at which
they influence the appearance of the BP.

A neutron scattering study by Towey et al. [38, 39, 37] has shown that at a
concentration of $5\text{\,}\,\mathrm{w}\mathrm{t}\,\%$, water monomers are
distributed homogenously throughout the material. As the water concentration
is increased, water clusters coexist with unclustered water molecules. Past
studies suggest that the onset of mobility in amorphous water may lie in the
range of $120\text{\,}\mathrm{K}150\text{\,}\mathrm{K}$[51] and the PES of
bulk-like water is hence expected to comprise shallower minima than that of
glycerol where the onset of mobility lies at higher temperatures. The THz-TDS
results indicate a shift in the structural dynamics once clusters form,
accompanied by a reduced height of the potential energy barrier separating the
harmonic from the anharmonic regime.

As the water concentration is increased up to
$16.4\text{\,}\,\mathrm{w}\mathrm{t}\,\%$ water, the percentage of water
molecules that are part of clusters increases to $80\text{\,}\mathrm{\char
37\relax}$ [38, 39, 37]. In our THz-TDS measurements we do not see a
difference between $T^{*}$ for mixtures containing
$15\text{\,}\,\mathrm{w}\mathrm{t}\,\%$ and
$20\text{\,}\,\mathrm{w}\mathrm{t}\,\%$. This indicates that the presence of
water-water hydrogen bonds lowers the barriers on the PES, but once clusters
dominate, the intensity of anharmonic effects does not depend on the number of
clusters, but on their mean size.

Water-water cooperative domains [52, 53] and percolating water clusters [35]
have been found at even higher water content. In our THz-TDS measurements, the
onset temperature of anharmonic effects is decreased noticeably in samples
with a water concentration of $30\text{\,}\,\mathrm{w}\mathrm{t}\,\%$. We can
hence infer that the mean size of clusters influences the PES and thereby
structural dynamics stronger than their number and that the height of the
lowest energy barrier is the limiting factor determining the onset temperature
of anharmonicity. Regions with lower energy barriers first exhibit an increase
in mobility which is measured by THz-TDS.

####  III.2.2 Anharmonic Regime ($T_{g}<T<T_{m}$)

Below $T^{*}$, the rate of absorption change does not strongly depend on water
content, as shown in Figure 7c. Above $T_{\text{g}}$, the rate of absorption
change starts to increase with water concentration above
$5\text{\,}\mathrm{w}\mathrm{t}\,\mathrm{\char 37\relax}\
\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}$ (see Figure 7c). There is
also evidence for a shallow minimum at around
$5\text{\,}\,\mathrm{w}\mathrm{t}\,\%$ which could be due to a very homogenous
sample including many unclustered water molecules.

Unclustered water reduces the rate of change of the anharmonicity parameter at
a water content below $30\text{\,}\,\mathrm{w}\mathrm{t}\,\%$ compared to pure
glycerol. This “plasticiser effect” could be because the unclustered water
molecules are surrounded by larger glycerol molecules, and any reorientation
of the glycerol molecules first requires the breaking of hydrogen bonds with
water. Once clusters are formed, however, the larger mobility of water
molecules in a bulk-like environment increases the rate again.

The higher the temperature and mobility, the more anharmonic effects obscure
the BP until it can no longer be resolved. $T_{\text{g}}$ does not depend on
water content below $10\text{\,}\,\mathrm{w}\mathrm{t}\,\%$, and decreases for
water concentrations up to $30\text{\,}\,\mathrm{w}\mathrm{t}\,\%$. This
indicates that the presence of a number of larger water clusters increases
anharmonic effects. Unclustered water and small isolated clusters have no
effect.

Raman scattering experiments have shown that an underlying BP continues to be
present above $T_{\text{g}}$ [50] which can no longer be resolved by THz-TDS
because it is obscured by anharmonic effects.

The calorimetric $T_{\text{g}}$ as measured with DSC shows a dependence on
water concentration even below $10\text{\,}\,\mathrm{w}\mathrm{t}\,\%$[35].
This can be attributed to the different criteria for determining the exact
value. While the calorimetric $T_{\text{g}}$ is related to a change in heat
flow, our definition based on anharmonic effects is therefore expected to
result in slightly different values.

At a temperature of approximately $30\text{\,}\mathrm{K}$ above $T_{\text{g}}$
lies the melting temperature of the crystalline system and by further heating
the supercooled liquid, it enters the liquid regime.

####  III.2.3 Liquid Regime ($T=$293\text{\,}\mathrm{K}$$)

A room temperature ($293\text{\,}\mathrm{K}$) set-up was utilised to measure
glycerol-water mixtures in the same concentration range
($0\text{\,}\mathrm{w}\mathrm{t}\,\mathrm{\char 37\relax}\
\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}30\text{\,}\mathrm{w}\mathrm{t}\,\mathrm{\char
37\relax}\ \mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}$) as discussed
above. With this setup it was easy to repeat measurements very rapidly and the
spectra shown in Figure S4a are the average of 18 measurements each. The
absolute error associated with the absorption coefficient measurement is hence
decreased and the data quality and spectral range are better than for
measurements performed at variable temperatures.

At room temperature, the mixtures are fully liquid and no BP is present. The
absorption coefficient at $1\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$ depends
linearly on water concentration $c$, as can be seen in Figure S4b:

$\alpha(c)=$1.47\text{\,}{\mathrm{cm}}^{-1}\text{\,}{\mathrm{\char
37\relax}}^{-1}$\cdot c[\,wt.\%]+$65.6\text{\,}{\mathrm{cm}}^{-1}$$ (12)

and the presence or absence of water clusters does not seem to influence this
relationship. To investigate how water clusters influence the THz-TDS spectra
of liquid glycerol-water mixtures, the frequency-dependent absorption is hence
evaluated. The high data quality allowed to extrapolate the absorption
coefficient measured in the range
$0.6\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}2.5\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$
according to a model developed by Schirmacher et al. [29]. However, the model
only makes predictions, and in a future experiment, the spectra could be
measured on a THz-TDS system with a higher spectral bandwidth.

The better data quality allows to fit a model to the data in the range
$0.6\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}2.5\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$
to extrapolate the spectra further and investigate the influence of
concentration:

$\alpha=A\cdot\nu^{2}\cdot\exp{(-\nu/(\nu_{\text{c}}/2))}+C$ (13)

Which predicts a broad peak centred at $\nu_{\text{c}}$. The model is fitted
to experimental data restricted to frequencies below
$2.5\text{\,}\mathrm{THz}$ because with increasing frequency the noise
increases [54]. Measured spectral data are reduced to the fit parameters and
the influence of water concentration is examined.

The parameters $\alpha_{\text{max}}$ (maximum absorption) and the full width
at half maximum (FWHM) describing the spectra are determined for different
concentrations as shown in Figure 8 with the concentration range below
$8.5\text{\,}\,\mathrm{w}\mathrm{t}\,\%$ water shown in detail in the insets.
As expected, an increase in water content leads to an increase in absorption,
as well as to a broadening of the peak. This also leads to the previously
discussed increase in the absorption at
$1\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$.

Figure 8: (a) Measured spectra (up to
$2.5\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$) and extrapolated spectra at
higher frequencies.
(b)-(d) Predicted centre frequency (a), maximum absorption (b), and FWHM of
the spectral shape for different water contents. Error bars are calculated
from 95% confidence interval of the fit. The insets show the concentration
range
$0\text{\,}\,\mathrm{w}\mathrm{t}\,\%8.5\text{\,}\,\mathrm{w}\mathrm{t}\,\%$
in more detail.

Unclustered water molecules are tightly integrated into the hydrogen-network,
with most glycerol hydrogen bonds formed between water and glycerol molecules
[35]. Above the glass transition temperature, these water molecules exhibit
different degrees of freedom which broadens the spectral feature at higher
frequencies. Increasing the amount of unclustered water leads to the most
noticeable change in those three parameters over the whole concentration range
studied.

At $6\text{\,}\,\mathrm{w}\mathrm{t}\,\%$, the rate of change in maximum
absorption with water content decreases slightly. Once water clusters form,
the number of water-glycerol hydrogen bonds decreases and the number of water-
water bonds increases. Once clusters are present, a higher water content
changes the dynamics less because the water would only increase the cluster
size and thereby disrupt the hydrogen-bonded network of glycerol molecules
less. Within the clusters, the water dynamics are markedly different from
unclustered water molecules that are hydrogen-bonded to glycerol.

While the second clustering threshold just below
$20\text{\,}\,\mathrm{w}\mathrm{t}\,\%$ was not clearly observed in the BP, it
influences the behaviour of the system at higher frequencies. All three
parameters describing the VDOS change less above a water content of
$20\text{\,}\,\mathrm{w}\mathrm{t}\,\%$. Most notably, the centre frequency
plateaus above $20\text{\,}\,\mathrm{w}\mathrm{t}\,\%$. At room temperature,
bulk water has a higher absorption at
$1\text{\,}\mathrm{T}\mathrm{H}\mathrm{z}$ than glycerol
($220\text{\,}{\mathrm{cm}}^{-1}$ compared to $66\text{\,}{\mathrm{cm}}^{-1}$)
and the larger the water clusters, the more they behave like bulk water. Once
percolating water clusters form, they dominate the absorption of liquid
glycerol-water mixtures at higher frequencies comprising the VDOS.

##  IV Conclusions

While the terahertz spectra of amorphous solids and liquids are featureless,
they provide considerable information about glass transition temperatures and
the related mobility changes. They also offer insight into the vibrational
density of states and the BP.

We show that the infrared absorption coefficient measured with THz-TDS can be
theoretically related to the reduced Raman intensity
($\propto\alpha/\omega^{2}$) and the reduced DOS ($\propto\alpha/\omega^{3}$)
and the agreement with experimental results confirms this.

The temperature behaviour is characterised by three regimes: below the glass
transition temperature, the system is trapped in a deep minimum of the
potential energy landscape and exhibits mostly harmonic vibrations. The onset
of local (and, at slightly higher temperatures), global mobility is
accompanied by anharmonic excitations, obscuring the BP and eventually leading
to its dissolution. Once the system crosses over into the liquid regime, the
absorption spectra are dominated by quasi-elastic scattering.

The implications of this model are examined in detail on the example of
glycerol and the insights are applied to study the influence of clustered and
unclustered water.

In the harmonic regime, the shape of the PES is unaltered in the presence of
clustered and unclustered water molecules. Unclustered water molecules are
embedded into the glycerol matrix and do not decrease the onset temperature of
anharmonicity. A change in structural dynamics is observed at a water
concentration of approximately $5\text{\,}\,\mathrm{w}\mathrm{t}\,\%$,
corresponding to a transition from isolated water molecules distributed
homogeneously throughout the sample to the presence of small water clusters
and an increased number of water-water hydrogen bonds which lower the barriers
on the PES.

Interestingly, the intensity of anharmonic effects does not depend on the
number of water clusters but on their mean size.

Data acquired at room temperature is extrapolated according to a model by
Schirmacher et al. which allows to investigate the spectral shape for
different concentrations for liquid glycerol-water mixtures and confirms a
change of dynamics once percolating water clusters formed between
$15\text{\,}\,\mathrm{w}\mathrm{t}\,\%20\text{\,}\,\mathrm{w}\mathrm{t}\,\%$.

This methodology has great potential to be applied to other systems to
investigate the change of mobility and dynamics.

###### Acknowledgements.

JK thanks the EPSRC Cambridge Centre for Doctoral Training in Sensor
Technologies and Applications (EP/L015889/1) and AstraZeneca for funding.
Disclosures: ES is an employee of AbbVie. He participated in conceptual design
of the study and in reviewing and approval of the publication.

## References

  * Corsaro  _et al._ [2019] C. Corsaro, E. Fazio, and D. Mallamace, The stokes-einstein relation in watermethanol solutions, The Journal of Chemical Physics 150, 234506 (2019).
  * Götze [1998] W. Götze, The essentials of the mode-coupling theory for glassy dynamics, Condensed Matter Physics (1998).
  * Yomogida  _et al._ [2010a] Y. Yomogida, Y. Sato, R. Nozaki, T. Mishina, and J. Nakahara, Comparative dielectric study of monohydric alcohols with terahertz time-domain spectroscopy, Journal of Molecular Structure 981, 173 (2010a).
  * Yomogida  _et al._ [2010b] Y. Yomogida, Y. Sato, R. Nozaki, T. Mishina, and J. Nakahara, Comparative study of boson peak in normal and secondary alcohols with terahertz time-domain spectroscopy, Physica B: Condensed Matter 405, 2208 (2010b).
  * Yomogida  _et al._ [2010c] Y. Yomogida, Y. Sato, K. Yamakawa, R. Nozaki, T. Mishina, and J. Nakahara, Comparative dielectric study of pentanol isomers with terahertz time-domain spectroscopy, Journal of Molecular Structure 970, 171 (2010c).
  * Taraskin  _et al._ [2006] S. Taraskin, S. Simdyankin, S. Elliott, J. Neilson, and T. Lo, Universal features of terahertz absorption in disordered materials, Physical Review Letters 97, 055504 (2006).
  * Sibik and Zeitler [2016] J. Sibik and J. A. Zeitler, Direct measurement of molecular mobility and crystallisation of amorphous pharmaceuticals using terahertz spectroscopy, Advanced Drug Delivery Reviews 100, 147 (2016).
  * Schirmacher  _et al._ [1998] W. Schirmacher, G. Diezemann, and C. Ganter, Harmonic vibrational excitations in disordered solids and the “boson peak”, Physical Review Letters 81, 136 (1998).
  * Marruzzo  _et al._ [2013a] A. Marruzzo, S. Koehler, A. Fratalocchi, G. Ruocco, and W. Schirmacher, Vibrational anomalies and marginal stability of glasses, The European Physical Journal Special Topics 216, 83 (2013a).
  * Schirmacher  _et al._ [2015] W. Schirmacher, T. Scopigno, and G. Ruocco, Theory of vibrational anomalies in glasses, Journal of Non-Crystalline Solids 407, 133 (2015).
  * Goldstein [1969] M. Goldstein, Viscous liquids and the glass transition: a potential energy barrier picture, The Journal of Chemical Physics 51, 3728 (1969).
  * Ruggiero  _et al._ [2017] M. T. Ruggiero, M. Krynski, E. O. Kissi, J. Sibik, D. Markl, N. Y. Tan, D. Arslanov, W. Van Der Zande, B. Redlich, T. M. Korter, _et al._ , The significance of the amorphous potential energy landscape for dictating glassy dynamics and driving solid-state crystallisation, Physical Chemistry Chemical Physics 19, 30039 (2017).
  * Markelz  _et al._ [2007] A. G. Markelz, J. R. Knab, J. Y. Chen, and Y. He, Protein dynamical transition in terahertz dielectric response, Chemical Physics Letters 442, 413 (2007).
  * Wuttke  _et al._ [1994] J. Wuttke, J. Hernandez, G. Li, G. Coddens, H. Cummins, F. Fujara, W. Petry, and H. Sillescu, Neutron and light scattering study of supercooled glycerol, Physical review letters 72, 3052 (1994).
  * Wang and Wright [1971a] C. Wang and R. Wright, Raman scattering studies of liquids and glasses. i. liquid and supercooled liquid glycerol, The Journal of Chemical Physics 55, 1617 (1971a).
  * Wang and Wright [1971b] C. Wang and R. Wright, Raman scattering studies in liquids and glasses. ii. liquid and supercooled liquid glycerol, The Journal of Chemical Physics 55, 3300 (1971b).
  * Dashnau  _et al._ [2006] J. L. Dashnau, N. V. Nucci, K. A. Sharp, and J. M. Vanderkooi, Hydrogen bonding and the cryoprotective properties of glycerol/water mixtures, The Journal of Physical Chemistry B 110, 13670 (2006).
  * Perova  _et al._ [1998] T. Perova, D. Christensen, U. Rasmussen, J. Vij, and O. Nielsen, Far-infrared spectra of highly viscous liquids: glycerol and triacetin (glycerol triacetate), Vibrational Spectroscopy 18, 149 (1998).
  * Murthy [2000] S. Murthy, Experimental study of the dynamics of water and the phase behavior of the supercooled aqueous solutions of propylene glycol, glycerol, poly (ethylene glycol) s, and poly (vinylpyrrolidone)., The Journal of Physical Chemistry B 104, 6955 (2000).
  * Huck  _et al._ [1988] J. Huck, G. Noyel, and L. Jorat, Dielectric properties of supercooled glycerol-water solutions, IEEE Transactions on Electrical Insulation 23, 627 (1988).
  * Dirama  _et al._ [2005] T. E. Dirama, G. A. Carri, and A. P. Sokolov, Role of hydrogen bonds in the fast dynamics of binary glasses of trehalose and glycerol: A molecular dynamics simulation study, The Journal of Chemical Physics 122, 114505 (2005).
  * Angell and Smith [1982] C. Angell and D. Smith, Test of the entropy basis of the vogel-tammann-fulcher equation. dielectric relaxation of polyalcohols near tg, The Journal of Physical Chemistry 86, 3845 (1982).
  * Sibik and Zeitler [2013] J. Sibik and J. A. Zeitler, Terahertz spectroscopy of hydrogen-bonded glass-forming liquids, in  _2013 38th International Conference on Infrared, Millimeter, and Terahertz Waves (IRMMW-THz)_ (IEEE, 2013) pp. 1–3.
  * Busselez  _et al._ [2009] R. Busselez, R. Lefort, M. Guendouz, B. Frick, O. Merdrignac-Conanec, and D. Morineau, Molecular dynamics of glycerol and glycerol-trehalose bioprotectant solutions nanoconfined in porous silicon, The Journal of Chemical Physics 130, 214502 (2009).
  * Schneider  _et al._ [1998] U. Schneider, P. Lunkenheimer, R. Brand, and A. Loidl, Dielectric and far-infrared spectroscopy of glycerol, Journal of Non-Crystalline Solids 235, 173 (1998).
  * Ryabov  _et al._ [2003] Y. E. Ryabov, Y. Hayashi, A. Gutina, and Y. Feldman, Features of supercooled glycerol dynamics, Physical Review B 67, 132202 (2003).
  * Sibik  _et al._ [2014] J. Sibik, S. R. Elliott, and J. A. Zeitler, Thermal decoupling of molecular-relaxation processes from the vibrational density of states at terahertz frequencies in supercooled hydrogen-bonded liquids, The Journal of Physical Chemistry Letters 5, 1968 (2014).
  * Capaccioli  _et al._ [2015] S. Capaccioli, K. Ngai, M. S. Thayyil, and D. Prevosto, Coupling of caged molecule dynamics to jg $\beta$-relaxation: I, The Journal of Physical Chemistry B 119, 8800 (2015).
  * Chumakov  _et al._ [2004] A. Chumakov, I. Sergueev, U. Van Buerck, W. Schirmacher, T. Asthalter, R. Rueffer, O. Leupold, and W. Petry, Collective nature of the boson peak and universal transboson dynamics of glasses, Physical Review Letters 92, 245508 (2004).
  * Wuttke  _et al._ [1995] J. Wuttke, W. Petry, G. Coddens, and F. Fujara, Fast dynamics of glass-forming glycerol, Physical Review E 52, 4026 (1995).
  * Sibik  _et al._ [2015] J. Sibik, K. Looebmann, T. Rades, and J. A. Zeitler, Predicting crystallization of amorphous drugs with terahertz spectroscopy, Molecular Pharmaceutics 12, 3062 (2015).
  * Sztein  _et al._ [2018] J. M. Sztein, T. Takeo, and N. Nakagata, History of cryobiology, with special emphasis in evolution of mouse sperm cryopreservation, Cryobiology 82, 57 (2018).
  * Cicerone  _et al._ [2015] M. T. Cicerone, M. J. Pikal, and K. K. Qian, Stabilization of proteins in solid form, Advanced Drug Delivery Reviews 93, 14 (2015).
  * Shalaev and Zografi [1996] E. Y. Shalaev and G. Zografi, How does residual water affect the solid-state degradation of drugs in the amorphous state?, Journal of Pharmaceutical Sciences 85, 1137 (1996).
  * Starciuc  _et al._ [2021] T. Starciuc, Y. Guinet, A. Hedoux, and E. Shalaev, Water content thresholds in glycerol/water system: Low-and high-wavenumber raman spectroscopy study, Journal of Molecular Liquids 321, 114678 (2021).
  * Murata and Tanaka [2012] K.-i. Murata and H. Tanaka, Liquid-liquid transition without macroscopic phase separation in a water-glycerol mixture, Nature Materials 11, 436 (2012).
  * Towey and Dougan [2012] J. Towey and L. Dougan, Structural examination of the impact of glycerol on water structure, The Journal of Physical Chemistry B 116, 1633 (2012).
  * Towey  _et al._ [2011] J. Towey, A. Soper, and L. Dougan, The structure of glycerol in the liquid state: a neutron diffraction study, Physical Chemistry Chemical Physics 13, 9397 (2011).
  * Towey  _et al._ [2012] J. Towey, A. Soper, and L. Dougan, Molecular insight into the hydrogen bonding and micro-segregation of a cryoprotectant molecule, The Journal of Physical Chemistry B 116, 13898 (2012).
  * Sibik  _et al._ [2013] J. Sibik, E. Y. Shalaev, and J. A. Zeitler, Glassy dynamics of sorbitol solutions at terahertz frequencies, Physical Chemistry Chemical Physics 15, 11931 (2013).
  * DiMarzio and Bishop [1974] E. A. DiMarzio and M. Bishop, Connection between the macroscopic electric and mechanical susceptibilities, The Journal of Chemical Physics 60, 3802 (1974).
  * Gemant [1935] A. Gemant, The conception of a complex viscosity and its application to dielectrics, Transactions of the Faraday Society 31, 1582 (1935).
  * Jakobsen  _et al._ [2005] B. Jakobsen, K. Niss, and N. B. Olsen, Dielectric and shear mechanical alpha and beta relaxations in seven glass-forming liquids, The Journal of Chemical Physics 123, 234511 (2005).
  * Schirmacher  _et al._ [2016] W. Schirmacher, G. Ruocco, and V. Mazzone, Theory of heterogeneous viscoelasticity, Philosophical Magazine 96, 620 (2016).
  * Marruzzo  _et al._ [2013b] A. Marruzzo, W. Schirmacher, A. Fratalocchi, and G. Ruocco, Heterogeneous shear elasticity of glasses: the origin of the boson peak, Scientific Reports 3, 1 (2013b).
  * Martin and Brenig [1974] A. Martin and W. Brenig, Model for brillouin scattering in amorphous solids, Physica Status Solidi (b) 64, 163 (1974).
  * Schmid and Schirmacher [2008] B. Schmid and W. Schirmacher, Raman scattering and the low-frequency vibrational spectrum of glasses, Physical Review Letters 100, 137402 (2008).
  * Surovtsev and Sokolov [2002] N. Surovtsev and A. Sokolov, Frequency behavior of raman coupling coefficient in glasses, Physical Review B 66, 054205 (2002).
  * Wuttke  _et al._ [1996] J. Wuttke, W. Petry, and S. Pouget, Structural relaxation in viscous glycerol: Coherent neutron scattering, The Journal of Chemical Physics 105, 5177 (1996).
  * Uchino and Yoko [1996] T. Uchino and T. Yoko, Low-frequency raman scattering and the fast relaxation process in glycerol, Science 273, 480 (1996).
  * Ewing  _et al._ [2015] S. Ewing, A. Hussain, G. Collins, C. Roberts, and E. Shalaev, Low-temperature mobility of water in sugar glasses: insights from thermally stimulated current study, in  _Water stress in biological, chemical, pharmaceutical and food systems_ (Springer, 2015) pp. 75–87.
  * Hayashi  _et al._ [2006] Y. Hayashi, A. Puzenko, and Y. Feldman, Slow and fast dynamics in glycerol–water mixtures, Journal of non-crystalline solids 352, 4696 (2006).
  * Sudo  _et al._ [2002] S. Sudo, M. Shimomura, N. Shinyashiki, and S. Yagihara, Broadband dielectric study of $\alpha$–$\beta$ separation for supercooled glycerol–water mixtures, Journal of non-crystalline solids 307, 356 (2002).
  * Jepsen and Fischer [2005] P. U. Jepsen and B. M. Fischer, Dynamic range in terahertz time-domain transmission and reflection spectroscopy, Optics Letters 30, 29 (2005).


# Inferring potential landscapes from noisy time series data

Pressé lab  Department of Physics, Arizona State University, USA  Bechhoefer
lab  Department of Physics, Simon Fraser University, CA

###### Abstract

While particle trajectories encode information on potentials dictating their
dynamics, the potentials themselves are often challenging to robustly
determine from such trajectories. Measurement errors corrupt the particle’s
position, sparse sampling of the potential landscape by the particle limits
data on in higher-energy regions such as barriers, and non-equilibrium
trajectories cannot be used in methods of obtaining potentials with an
equilibrium assumption. To overcome these difficulties, here we develop a
general framework to infer, simultaneously and self-consistently, a potential
of arbitrary shape along with the particle’s friction coefficient and
measurement noise properties that does not require the trajectories to be at
equilibrium. To accommodate any potential form, we invoke a Bayesian
nonparametric formulation. Within the Bayesian paradigm, we develop an
efficient Monte Carlo scheme to sample from the full posterior over all
candidate potentials. Our method is validated on experiments performed on
particles evolving within a feedback trap in both one and two dimensions.

††preprint: AIP/123-QED

##  I Introduction

Determining potentials governing particle dynamics is a fundamental question
across the biological science makarov2015shapes; wang2016nonlinear;
wang1997configurational; wang2003energy; chu2013quantifying, materials science
deringer2019machine; handle2018potential, and beyond perez2021forma;
dudko2006intrinsic; preisler2004modeling. The motivation for learning
potentials, often effective potentials or free energy landscapes, from a
trajectory is many-fold. For example, potentials provide reduced dimensional
descriptions of dynamics along a reaction coordinate wang2016nonlinear;
wang2003energy; chu2013quantifying and yield key estimates over thermodynamic
and kinetic quantities hanggi1990reaction; berezhkovskii2017mean;
bessarab2013potential.

… MENTION OPTICAL TWEEZERS AND FEEDBACK TRAP …

As inferring potentials is a general problem, many strategies have been
devised to deduce potentials from data reif2009fundamentals;
turkcan2012bayesian; perez2021forma; wang2019machine. These strategies,
however, are limited: some methods reif2009fundamentals; turkcan2012bayesian
necessarily rely on binned data, relating potential energies to Boltzmann
weights or average apparent force, thereby limiting the frequency of data in
each bin and requiring that equilibrium be reached before data acquisition.
Other methods perez2021forma assume harmonic analytic forms not otherwise
generalizable, while still others use neural networks wang2019machine to
deduce potentials for which the uncertainty derived from measurement error and
data sparsity is not easily propagated into uncertainty over the inferred
potential.

To address many of these challenges, in previous work bryan2020inferring, we
introduced a method starting from noiseless one-dimensional time series data
to infer effective potential landscapes without binning, nor assuming a
potential form, nor assuming equilibrium conditions, while admitting full
posterior inference (and thus error bars or, equivalently, credible intervals)
over any candidate potentials arising from sparse data.

However, our method was limited in three critical ways: 1) it was limited to
one dimension; 2) it ignored measurement error and implications for
uncertainty about the potential; 3) it could not be benchmarked on real data
since measurement error was ignored. In fact, no existing method can infer
potentials from trajectories corrupted by measurement error.

Here we introduce a method to learn potentials from noisy, multidimensional,
non-equilibrium data. By taking advantage of tools from Bayesian
nonparametrics, our method can learn the shapes of potentials from
trajectories while meeting all the following criteria simultaneously: 1) no
reliance on binning or pre-processing; 2) no assumed analytic potential form;
3) inferences drawn from a posteriors allowing for spatially nonuniform
uncertainties to be informed by local density of available data in specific
regions of the potential (e.g., fewer data points around barriers); 4)
treatment of multidimensional trajectories; 5) rigorous incorporation of
measurement noise through likelihoods; and 6) compatible with non-equilibrium
trajectories. No other existing method meets each of these 6 criteria
simultaneously.

… EXPLAIN USEFULNESS IN OPTICAL TRAPS …

##  II Methods

Here we first formulate the dynamical model (overdamped Langevin) and describe
the associated forward model for measurement noise.

Our goal is to use noisy positional measurements, $\bm{y}_{1:N}$, to learn the
potential at each point in space, $U(\cdot)$, the friction coefficient,
$\zeta$, the magnitude of the measurement noise, $\sigma^{2}$ (under a
Gaussian noise model), and the actual position at each time, $\bm{x}_{1:N}$.

We borrow functional notation $U(\cdot)$ to avoid confusion between the
function evaluated at a point $x$ ($U(x)$) versus the function as a whole over
the entire real axis.

Toward achieving our goal, after formulating the forward model, we develop a
strategy for constructing posterior probability distributions over all the
unknowns listed above. As our posterior does not admit an analytic form, we
devise an efficient Monte Carlo strategy to sample from this posterior.

###  II.1 Dynamics

We begin with an overdamped Langevin description of the particle dynamics
zwanzig2001nonequilibrium,

$\displaystyle\zeta\dot{\bm{x}}$ $\displaystyle=\bm{f}(\bm{x})+\bm{r}(t)$ (1a)
$\displaystyle\bm{f}(\bm{x})$ $\displaystyle=-\bm{\nabla}U(\bm{x}),$ (1b)

where $\bm{x}(t)$ is the possibly multidimensional position coordinate at time
$t$; $\dot{\bm{x}}(t)$ is the velocity; $\bm{f}(\bm{x})$ is the force at
position $\bm{x}(t)$; and $\zeta$ is the friction coefficient. The forces
acting on the particle include positional forces $\bm{f}(\bm{x})$ expressed as
the gradient of a conservative potential,
$\bm{f}(\bm{x})=-\bm{\nabla}U(\bm{x})$, and the stochastic (thermal) force,
$\bm{r}(t)$, defined as follows:

$\displaystyle\langle\bm{r}(t)\rangle$ $\displaystyle=\bm{0}$ (2a)
$\displaystyle\langle\bm{r}^{T}(t)\bm{r}(t^{\prime})\rangle$
$\displaystyle=2\zeta kT\delta(t-t^{\prime})$ (2b)

where $\langle\cdot\rangle$ denotes an ensemble average over realizations, $T$
is the temperature of the bath and $k$ is Boltzmann’s constant.

Under a forward Euler scheme leveque2007finite for Eq. 1a with time points
given by $t_{n}=n\Delta t$, each position, given its past realization, is
sampled from a normal distribution

$\displaystyle\bm{x}_{n+1}|\bm{x}_{n},\bm{f}(\cdot),\zeta\sim$
$\displaystyle\textbf{Normal}\left(\bm{x}_{n}+\frac{\Delta
t}{\zeta}\bm{f}(\bm{x}_{n}),~{}\frac{2\Delta tkT}{\zeta}\bm{I}\right).$ (3)

In words, “the position $\bm{x}_{n+1}$ given quantities
$\bm{x}_{n},\bm{f}(\cdot),$ and $\zeta$ is sampled from a Normal distribution
with mean $\bm{x}_{n}+\frac{\Delta t}{\zeta}\bm{f}(\bm{x}_{n})$ and variance
$\frac{2\Delta tkT}{\zeta}\bm{I}$.”

Next, we describe the measurement, $\bm{y}_{n}$, which is probabilistically
related to the position, $\bm{x}_{n}$. As is typical for experimental setups
(e.g. section III.1), we use a Gaussian noise model and write

$\displaystyle\bm{y}_{n}|\bm{x}_{n},\sigma^{2}\sim$
$\displaystyle\textbf{Normal}\left(\bm{x}_{n},~{}\sigma^{2}\bm{I}\right)$ (4)

where $\sigma^{2}$ is the measurement noise variance and the indexing notation
for the positions carries over for the measurements. In Eq. 4, the measurement
process is instantaneous, i.e. assumed to be faster than the dynamical time
scales.

Our choice of Gaussian measurement model here can be modified at minimal
computational cost (e.g., we can select more complex noise models
hirsch2013stochastic if warranted by the data).

###  II.2 Probabilities

From Eqs. 3 and 4, we wish to obtain the probability distribution over the
potential, $U(\cdot)$, the friction coefficient, $\zeta$, and the positions at
each time, $\bm{x}_{1:N}$, given the data, i.e., the sequence of noisy
measurements, $\bm{y}_{1:N}$.

In a Bayesian formulation, the probability of the model parameters is given by
the product of the likelihood
($\mathcal{P}\left(\bm{y}_{1:N}|U(\cdot),\zeta,\bm{x}_{1:N},\sigma^{2}\right)$)
and $\mathcal{P}\left(x_{1:N}|U(\cdot),\zeta\right)$) and the priors
($\mathcal{P}\left(U(\cdot)\right)$,
$\mathcal{P}\left(x_{1}|U(\cdot),\zeta\right)$,
$\mathcal{P}\left(\zeta\right)$, and $\mathcal{P}\left(\sigma^{2}\right)$)

$\displaystyle\mathcal{P}\left(U(\cdot),\zeta,\bm{x}_{1:N},\sigma^{2}|\bm{y}_{1:N}\right)$
$\displaystyle\propto\mathcal{P}\left(\bm{y}_{1:N}|U(\cdot),\zeta,\bm{x}_{1:N},\sigma^{2}\right)\mathcal{P}\left(U(\cdot),\zeta,\bm{x}_{1:N},\sigma^{2}\right).$
(5)

To simplify our presentation, we have excluded the normalization,
$\left[\mathcal{P}\left(\bm{y}_{1:N}\right)\right]^{-1}$.

Using Eq. 4, we can write the measurement portion of the likelihood as a
product over independent measurements,

$\displaystyle\mathcal{P}\left(\bm{y}_{1:N}|U(\cdot),\zeta,\bm{x}_{1:N},\sigma^{2}\right)$
$\displaystyle=\prod_{n=1}^{N}\textbf{Normal}\left(\bm{y}_{n};~{}\bm{x}_{n},~{}\sigma^{2}\bm{I}\right).$
(6)

We decompose the remaining probabilities as follows:

$\displaystyle\mathcal{P}\left(U(\cdot),\zeta,\bm{x}_{1:N},\sigma^{2}\right)=$
$\displaystyle\mathcal{P}\left(\bm{x}_{1:N}|U(\cdot),\zeta\right)\mathcal{P}\left(U(\cdot)\right)\mathcal{P}\left(\zeta\right)\mathcal{P}\left(\sigma^{2}\right).$
(7)

The first term on the right-hand side of Eq. 7 can be further decomposed since
the trajectory, $\bm{x}_{1:N}$, results from the Langevin dynamics, Eq 1a,

$\displaystyle\mathcal{P}\left(\bm{x}_{1:N}|U(\cdot),\zeta\right)=$
$\displaystyle\mathcal{P}\left(\bm{x}_{1}|U(\cdot),\zeta\right)\prod_{n=2}^{N}\textbf{Normal}\left(\bm{x}_{n};~{}\bm{x}_{n-1}+\frac{\Delta
t}{\zeta}\bm{f}(\bm{x}_{n-1}),\frac{2\Delta tkT}{\zeta}\bm{I}\right).$ (8)

We are free to choose the remaining priors,
$\mathcal{P}\left(U(\cdot)\right)$,
$\mathcal{P}\left(x_{1}|U(\cdot),\zeta\right)$,
$\mathcal{P}\left(\zeta\right)$, and $\mathcal{P}\left(\sigma^{2}\right)$.

In particular, important considerations dictate the form of the prior on the
potential. First, the potential have any shape (so it must be nonparametric);
however, the potential at points that are close to each other should be
correlated. A Gaussian-process prior williams2006gaussian allows us to sample
continuous curves with covariance (correlation of points) given by a
prespecified kernel function. However, naive implementation of Gaussian
process priors are computationally expensive, with time and memory
requirements scaling like the number of data points cubed bryan2020inferring;
williams2006gaussian.

To allow for efficient computations when the Gaussian process is invoked for
larger number of data points, we use a structured-kernel-interpolation
Gaussian process (SKI-GP) wilson2015kernel; wilson2015thoughts;
titsias2009variational; gal2014variational prior for the potential,
$U(\cdot)$. The SKI-GP prior is a hierarchical structure where the potential
at all points is interpolated according to $M$ chosen inducing points at fixed
locations, $\bm{x}^{u}_{1:M}$, and where the values of the potential at the
inducing points are drawn from a Gaussian process (or a multivariate
Gaussian). That is,

$\displaystyle U(\bm{x})=$
$\displaystyle\bm{K^{\dagger}}\bm{K}^{-1}\bm{u}_{1:M}$ (9a)
$\displaystyle\bm{u}_{1:M}\sim$
$\displaystyle\textbf{Normal}\left(\bm{0},~{}\bm{K}\right)$ (9b)
$\displaystyle K^{\dagger}_{m}=$ $\displaystyle K(\bm{x},\bm{x}^{u}_{m})$ (9c)
$\displaystyle K_{ij}=$ $\displaystyle K(\bm{x}^{u}_{i},\bm{x}^{u}_{j})$ (9d)

where $\bm{K^{\dagger}}$ is the covariance between inducing points and the
point at which the potential is being evaluated (the $m$th element of
$\bm{K}^{\dagger}$ is $K^{\dagger}_{m}$), $\bm{K}$ is the covariance between
inducing points, and $K(\cdot,\cdot)$ is a kernel function (which is often
chosen to be quadratic exponential in form williams2006gaussian).

Notice that under this model, we shift the focus from learning $U(\cdot)$ to
learning $\bm{u}_{1:M}$. With $\bm{u}_{1:M}$, we can recover $U(\cdot)$ and
$\bm{f}(\cdot)$ using a modified kernel matrix (See SI section V.1).

For $\mathcal{P}\left(x_{1}|U(\cdot),\zeta\right)$, we are free to choose any
prior and select

$\displaystyle\mathcal{P}\left(\bm{x}_{1}|U(\cdot),\zeta\right)=\textbf{Normal}\left(\bm{x}_{1};0,\Theta^{2}\right).$
(10)

For the friction, $\zeta$, we select a gamma distribution for which support on
the positive real axis is assured. That is,

$\displaystyle\zeta\sim$
$\displaystyle\textbf{Gamma}\left(\alpha_{\zeta},\beta_{\zeta}\right).$ (11)

Lastly, for the noise variance, $\sigma^{2}$, we place a inverse gamma prior

$\displaystyle\sigma^{2}\sim$
$\displaystyle\textbf{InvGamma}\left(\alpha_{\sigma^{2}},\beta_{\sigma^{2}}\right).$
(12)

The choice of inverse gamma prior is conjugate to the likelihood (Eq. 6)
meaning that the posterior constructed by the prior multiplied by the
likelihood can be directly sampled from bishop2006pattern. The variables,
$\Theta$, $\alpha_{\zeta}$, $\beta_{\zeta}$, $\alpha_{\sigma^{2}}$, and
$\beta_{\sigma^{2}}$ are hyperparameters that we are free to choose, and whose
impact on the ultimate shape of the posterior reduces as more data are
collected gelman2013bayesian; sivia2006data).

###  II.3 Inference

Using Eq. 5, the goal is to learn the potential, $U(\cdot)$, the friction
coefficient, $\zeta$, the variance of the measurement noise, $\sigma^{2}$, and
the position at each time level, $\bm{x}_{1:N}$.

As our posterior does not have an analytical expression, we devise a Gibbs
sampling scheme. Within this scheme, we start with an initial set of values
for the parameters and then iteratively sample each variable holding all
others fixed geman1984stochastic. We summarize the algorithm below:

  * •

Step 1: Choose initial
$\zeta^{(0)},\sigma^{(0)},\bm{x}_{1:N}^{(0)},U(\cdot)^{(0)}$.

  * •

Step 2: For all iterations, $i$:

    * –

Sample a new friction,
$\zeta^{(i)}\sim\mathcal{P}\left(\zeta^{(i)}|\sigma^{(i-1)},\bm{x}_{1:N}^{(i-1)},U(\cdot)^{(i-1)}\right)$.

    * –

Sample a new noise,
$\sigma^{(i)}\sim\mathcal{P}\left(\sigma^{(i)}|\zeta^{(i)},\bm{x}_{1:N}^{(i-1)},U(\cdot)^{(i-1)}\right)$.

    * –

Sample a new trajectory,
$\bm{x}_{1:N}^{(i)}\sim\mathcal{P}\left(\bm{x}_{1:N}^{(i)}|\zeta^{(i)},\sigma^{(i)},U(\cdot)^{(i-1)}\right)$.

    * –

Sample a new potential,
$U(\cdot)^{(i)}\sim\mathcal{P}\left(U(\cdot)^{(i)}|\zeta^{(i)},\sigma^{(i)},\bm{x}_{1:N}^{(i)}\right)$.

See SI section V.2 for how we derive and sample from the required conditional
probabilities.

##  III Results

Here we benchmark our method on experimental data. We start with data acquired
from a simple harmonic potential and show that we can accurately learn the
shape of the potential along with the friction coefficient, $\zeta$, and
measurement noise, $\sigma^{2}$. We show that increasing the width of the
measurement noise increases the uncertainty of all learned quantities, but
does not prohibit us from learning the potential accurately especially in
regions where it is well sampled. We then introduce a particle into a more
complicated double-well potential and show that our method can still learn the
shape of this potential. Next, we look at two-dimensional time series data and
demonstrate that we can learn multidimensional potentials. A robustness test
over parameters of interest using simulated data can be found in the SI.

Ground truth, with which we can validate our method, is known as the particle
can be evolved in a feedback trap with a known, sculpted, potential and
calibrated noise as described in the next section.

… NOTE THAT THE FIGURE TITLES ARE NAMED AFTER THE DATA SET AND NOT THE
MAGNITUDE OF THE MEASUREMENT NOISE. WE WILL NEED TO CHANGE THIS LATER …

###  III.1 Data acquisition

We collected data using a feedback trap

… EXPLAIN DATA ACQUISITION …

… EXPLAIN ESTIMATION OF GROUND TRUTH TRAJECTORY …

… CAN WE JUSTIFY GAUSSIAN NOISE MODEL? …

###  III.2 Single well

Figure 1: Demonstration on data from a harmonic potential. Here we analyze 5
data sets with increasing measurement noise. Each column shows results for a
different data set. For each data set we plot the learned potential in the top
row along with the ground truth and the results of the Boltzmann method. We
plot the learned trajectory with uncertainty against the ground truth in the
bottom row. For clarity we zoom into a region of the trajectory (200ms to
201ms).

We first tested our method for data acquired from a particle evolving in an
approximately single harmonic well potential. Results are shown in Fig. 1.
Each column shows the learned potential (top row) and learned trajectory
(bottom row) for each data set analyzed. We provide uncertainties and ground
truth estimates for both the potential and trajectory. Additionally, for sake
of comparison, we also show the potential estimated using the Boltzmann method
bryan2020inferring; reif2009fundamentals (see SI section V.3). We highlight
that the Boltzmann method does not provide trajectory estimates. By contrast,
our method infers those positions obscured by noise.

Figure 1 shows that the ground truth potential and trajectory fall within our
error bars (credible interval) for all data sets analyzed. On the other hand,
the Boltzmann method (see details in SI section V.3) does well in the case of
minimal noise (Fig. 1 left column), but fails as the noise introduced grows
(Fig. 1 right column). This is expected as the Boltzmann method does not
differentiate between measurement noise and inherent Brownian motion of the
trajectory, thus falling victim to overestimating the potential in those
regions around the trajectory and underestimating the potential along the
trajectory because when the particle is measured in the wrong bin this both
places a measurement in the wrong bin and takes away a measurement from the
true bin.

We note that the estimate of the potential falls off at the spatial edges of
the region sampled by the particle. This is because at the edges with high
potential, rarely visited by the particle, we do not have sufficient
information through the likelihood to inform those regions. As such the
potential inferred reverts back to the prior (set at 0 as described in SI
section II.2).

###  III.3 Double well

Figure 2: Demonstration on data from a double well potential. Here we analyze
5 data sets with increasing measurement noise. Each column shows results for a
different data set. For each data set, we plot the learned potential in the
top row alongside the ground truth and results of the Boltzmann method. We
plot the learned trajectory against the ground truth in the bottom row. For
clarity we zoom into a region of the trajectory (200ms to 201ms).

Next, we analyzed data from a particle in a double well potential. Results are
shown in Fig. 2. Again, the ground truth potential and trajectory fall withing
the estimated range even when the measurement noise is so large that the
particle is often seen in the wrong well (Fig. 2 top right panel). We note
that both our method and the Boltzmann method slightly overestimate the
potential of the left well at the lowest noise level. This is because, due to
finiteness of data, the trajectory spent a larger than expected amount of time
in the right well leaving the left well slightly undersampled (with both
methods concluding that it must be higher in energy than it is).

###  III.4 2D single well

Figure 3: Demonstration on data from a 2D harmonic potential. Here we analyze
4 data sets with increasing measurement noise. Each column shows the learned
potential results along with the ground truth potential for a different data
set. In the top we show the learned potential and ground truth plotted in 2D.
On the bottom, we show a 1D slice taken through the middle of the potential.

Lastly, we analyzed data from a particle in a 2D harmonic potential. Results
are shown in Fig. 3. For clarity we do not show uncertainties, trajectories,
nor Boltzmann method estimates for the 2D plot, but we do show them for a 1D
slice of the potential. Despite the added complexity of learning the potential
in 2D at once, our estimates qualitatively match with the ground truth even at
high measurement noise.

###  III.5 Non-equilibrium trajectories

Figure 4: Demonstration on data from non-equilibrium trajectories. Here we
analyze 4 data sets with increasing barrier heights taken from non-equilibrium
trajectories. For each data set we analyze multiple trajectories in which we
start the particle at the top of the well and watch it ”roll off” for a short
duration before starting a new trajectory. We are able to faithfully
reconstruct the shape of each potential despite the fact that the Boltzmann
method completely fails. For comparison, the learned and ground truth energy
landscapes of each panel were shifted so that the energy is 0 at the origin.

One advantage to our method is that it does not rely on an equilibrium
assumption. As such, we can analyze trajectories starting from non-equilibrium
conditions. To demonstrate this we created datasets in which we start the
particle at the top of the potential well and watch it ”roll off” to either
side. The trajectories are short (1$ms$) so that the particle does not reach
equilibrium during the time trace. We record 50 trajectories for each barrier
height. By including the likelihoods from multiple trajectories into our
likelihood we can gain information on either side of the well and recreate the
full potential.

Figure 4 shows the results. Note that in all four experiments, our method
learns that there is a barrier located at the origin, whereas the Boltzmann
histogram method predicts a well at the origin. This is because, as all
trajectories start at the origin, the histogram method will impart a low
energy at the origin. On the other hand, our method does not rely on
equilibrium assumption and learns the potential based off of the size of the
thermal kicks, thus we do not fall victim to the limitation imposed on the
Boltzmann method.

For the low barrier (figure 4 left panel), our method accurately predicts the
shape of the right barrier, but underestimates the depth of the left well
because the particle rolled to the left side few times and stochastic events
in the few trajectories had a large influence on the total inference. On the
high potential barrier (figure 4 right panel) we underestimate the depth of
the well because most of the trajectories ended before they could explore the
bottom of the well and thus we do not have much information about the far side
of the well. For the middle two experiments of figure 4, the ground truth is
within our uncertainty estimate.

##  IV Discussion

Inferring potential landscapes is a key step toward providing a reduced
dimensional description of complex systems wang2019machine;
chmiela2017machine; espanol2011obtaining; izvekov2005multiscale;
manzhos2015neural. Here, we go beyond existing methods by providing a means to
obtain potentials from multidimensional time series data corrupted by
measurement noise. Our method is capable of learning potentials along with the
positions (treated as latent variables hidden by measurement noise), friction
coefficient, and measurement noise variance simultaneously (assuming a
Gaussian noise model for simplicity alone).

We showed that we could accurately learn the potential landscape even if the
noise was of the same order of magnitude as the size of the confining
potential well and if the trajectories started away from equilibrium. We
benchmarked our method on data from a feedback trap. Additional synthetic data
sets is provided in the SI section V.4.

Our method provides an added advantage over other methods in treating
measurement noise without pre-processing (e.g., binning), without assuming an
analytical potential form, without requiring equilibrium conditions, and that
allowing for rigorous error propagation through a likelihood.

In particular, our method is advantageous over methods perez2021forma assuming
analytic potential forms and methods based on neural nets manzhos2015neural
that cannot currently propagate known error. Importantly, unlike the Boltzmann
method, our method does not invoke any equilibrium assumption.

… DISCUSSION RELATED TO OPTICAL TRAP? …

Our method, as it stands, uses a Gaussian noise model. We could, in principle,
use different measurement likelihoods by modifying Eq. 6. Additionally, our
model could be generalized to likelihoods beyond noisy measurements. In fact,
moving forward it would be interesting to analyze data from, for example,
smFRET (single molecule fluorescence resonance energy transfer) experiments,
where the position of the particle at each point in time is indirectly related
to the smFRET efficiency over time. This application would be relevant to
inferring protein conformational dynamics or binding kinetics
schuler2008protein; chung2018protein; sturzenegger2018transition.

… APPLICATION TO FORCE SPECTROSCOPY? …

##  V Supplemental information

###  V.1 Modified kernel matrix

Here we show how to calculate the force at an arbitrary test location using
the structured kernel interpolation Gaussian process. We first derive the
force for one-dimensional potentials exactly and then briefly explain how to
generalize to higher dimensions. In order to calculate the force given
$\bm{u}_{1:M}$ we write,

$\displaystyle f(x)=$ $\displaystyle-\frac{d}{dx}U(x)$ (13a) $\displaystyle=$
$\displaystyle-\frac{d}{dx}\bm{K^{\dagger}}\bm{K}^{-1}\bm{u}_{1:M}$ (13b)
$\displaystyle=$ $\displaystyle\bm{K^{*}}\bm{K}^{-1}\bm{u}_{1:M}$ (13c)
$\displaystyle\bm{K^{*}}=$ $\displaystyle-\frac{d}{dx}\bm{K^{\dagger}}$ (13d)
$\displaystyle K^{*}_{ij}=$ $\displaystyle-\frac{d}{dx}K(x,x^{u}_{m})$ (13e)
$\displaystyle=$
$\displaystyle\frac{h^{2}}{\ell^{2}}\left(x-x^{u}_{m}\right)\exp\left(-\frac{1}{2}\left(\frac{x-x^{u}_{m}}{\ell}\right)^{2}\right)$
(13f)

where $\bm{K^{*}}$ is now the covariance between the potential evaluated at
the inducing points and the force evaluated at a test location, and we have
assumed the kernel follows the familiar squared exponential form
[williams2006gaussian].

To generalize to higher dimensions, we will need to calculate the force in
each dimension separately. For example, to find the force in the $k$ direction
at a test point, $\bm{x}$,

$\displaystyle f^{k}(\bm{x})=$ $\displaystyle-\frac{d}{dx^{k}}U(\bm{x})$ (14)
$\displaystyle=$
$\displaystyle-\frac{d}{dx^{k}}\bm{K^{\dagger}}\bm{K}^{-1}\bm{u}_{1:M}$ (15)
$\displaystyle=$ $\displaystyle\bm{K^{*k}}\bm{K}^{-1}\bm{u}_{1:M}$ (16)
$\displaystyle K^{*k}_{ij}=$
$\displaystyle\frac{h^{2}}{\ell^{2}}\left(x^{k}-x^{u,k}_{m}\right)\exp\left(-\frac{1}{2}\left|\frac{\bm{x}-\bm{x}^{u}_{m}}{\ell}\right|^{2}\right).$
(17)

###  V.2 Conditional probabilities for Gibbs sampling

Here, we show the conditional probabilities used for our Gibbs sampling
algorithm.

####  V.2.1 Positions

The distribution for $\bm{x}_{1:N}$ is simple. If we sample each $x_{n}$ one
at a time,

$\displaystyle\mathcal{P}\left(x_{1}|\bm{u}_{1:M},\bm{x}_{2:N}\right)\propto$
$\displaystyle\textbf{Normal}\left(x_{1};~{}0,~{}\frac{2\tau
kT}{\zeta}\right)$ (18)
$\displaystyle\times\textbf{Normal}\left(x_{2};~{}x_{1}+\frac{\tau}{\zeta}f(x_{1}),~{}\frac{2\tau
kT}{\zeta}\right)$
$\displaystyle\times\textbf{Normal}\left(y_{1};~{}x_{1},~{}\sigma^{2}\right)$
$\displaystyle\mathcal{P}\left(x_{n}|\bm{u}_{1:M},\bm{x}_{1:n-1}\right)\propto$
$\displaystyle\textbf{Normal}\left(x_{n};~{}x_{n-1}+\frac{\tau}{\zeta}f(x_{n-1}),~{}\frac{2\tau
kT}{\zeta}\right)$ (19)
$\displaystyle\times\textbf{Normal}\left(x_{n+1};~{}x_{n}+\frac{\tau}{\zeta}f(x_{n}),~{}\frac{2\tau
kT}{\zeta}\right)$
$\displaystyle\times\textbf{Normal}\left(y_{n};~{}x_{n},~{}\sigma^{2}\right)$
$\displaystyle\mathcal{P}\left(x_{N}|\bm{u}_{1:M},\bm{x}_{1:N-1}\right)\propto$
$\displaystyle\textbf{Normal}\left(x_{N};~{}x_{N-1}+\frac{\tau}{\zeta}f(x_{N-1}),~{}\frac{2\tau
kT}{\zeta}\right)$ (20)
$\displaystyle\times\textbf{Normal}\left(y_{N};~{}x_{N},~{}\sigma^{2}\right),$

however, although these equations look Gaussian, they are not (except the last
one). We have to remember that $f_{n}=f(x_{n})$ is a function of $x_{n}$. Thus
we can only sample position using a Metropolis-Hastings algorithm
[bishop2006pattern].

####  V.2.2 Potential

Following the logic outlined in our previous manuscript [bryan2020inferring],
we can learn the potential. For now let us focus on one-dimensional data. Let
us fully write the prior on positions, $\bm{x}_{2:N}|U(\cdot),x_{1},\zeta$ as

$\displaystyle\mathcal{P}\left(\bm{x}_{2:N}|f(\cdot),x_{1}\right)=$
$\displaystyle\textbf{Normal}\left(\tau\bm{v}_{1:N-1};~{}\frac{\tau}{\zeta}\bm{f}_{1:N-1},~{}\frac{2\tau
kT}{\zeta}\bm{I}\right)$ (21) $\displaystyle\propto$
$\displaystyle\textbf{Normal}\left(\bm{f}_{1:N-1};~{}\zeta\bm{v}_{1:N-1},~{}\frac{2\zeta
kT}{\tau}\bm{I}\right)$ (22)

where $f_{n}=f(x_{n})$ and $v_{n}=(x_{n+1}-x_{n})/\tau$. Substituting in Eq.
9a, we get,

$\displaystyle=$
$\displaystyle\textbf{Normal}\left(\bm{K^{*}}\bm{K}^{-1}\bm{u}_{1:M};~{}\zeta\bm{v}_{1:N-1},~{}\frac{2\zeta
kT}{\tau}\bm{I}\right)$ (23) $\displaystyle\propto$
$\displaystyle\exp\left(-\frac{1}{2}\left(\zeta\bm{v}_{1:N-1}-\bm{K^{*}}\bm{K}^{-1}\bm{u}_{1:M}\right)^{T}\frac{\tau}{2\zeta
kT}\bm{I}\left(\zeta\bm{v}_{1:N-1}-\bm{K^{*}}\bm{K}^{-1}\bm{u}_{1:M}\right)\right)$
(24) $\displaystyle\propto$
$\displaystyle\exp\left(-\frac{1}{2}\frac{\tau}{2\zeta
kT}\bm{u}_{1:M}^{T}\bm{K}^{-1}\bm{K^{*}}{}^{T}\bm{K^{*}}\bm{K}^{-1}\bm{u}_{1:M}+\frac{1}{2}\frac{2\tau}{2kT}\bm{v}_{1:N-1}^{T}\bm{K^{*}}\bm{K}^{-1}\bm{u}_{1:M}\right)$
(25) $\displaystyle\propto$
$\displaystyle\exp\left(-\frac{1}{2}\left(\bm{u}_{1:M}-\bm{M}^{-1}\bm{b}\right)^{T}\bm{M}\left(\bm{u}_{1:M}-\bm{M}^{-1}\bm{b}\right)\right)$
(26) $\displaystyle\propto$
$\displaystyle\textbf{Normal}\left(\bm{u}_{1:M};~{}\bm{M}^{-1}\bm{b},~{}\bm{M}^{-1}\right)$
(27)

where we have used matrix completing the square with

$\displaystyle\bm{b}=$
$\displaystyle\frac{\tau}{2kT}\bm{K}^{-1}\bm{K^{*}}{}^{T}\bm{v}_{1:N-1}$ (28)
$\displaystyle\bm{M}=$ $\displaystyle\frac{\tau}{2\zeta
kT}\bm{K}^{-1}\bm{K^{*}}{}^{T}\bm{K^{*}}\bm{K}^{-1}.$ (29)

Combining with the prior, we get

$\displaystyle\mathcal{P}\left(u_{1:M}|\bm{x}_{1:N}\right)\propto$
$\displaystyle\textbf{Normal}\left(\bm{u}_{1:M};~{}\bm{M}^{-1}\bm{b},~{}\bm{M}^{-1}\right)\textbf{Normal}\left(\bm{u}_{1:M};~{}\bm{0},~{}\bm{K}\right)$
(30) $\displaystyle\propto$
$\displaystyle\textbf{Normal}\left(\bm{u}_{1:M};~{}\left(\bm{K}^{-1}+\bm{M}\right)^{-1}\left(\bm{K}^{-1}\bm{0}+\bm{M}\bm{M}^{-1}\bm{b}\right),~{}\bm{K}^{-1}+\bm{M}\right)$
(31) $\displaystyle=$
$\displaystyle\textbf{Normal}\left(\bm{u}_{1:M};~{}\tilde{\bm{\mu}},~{}\tilde{\bm{K}}\right)$
(32) $\displaystyle\tilde{\bm{\mu}}=$
$\displaystyle\frac{\tau}{2kT}\left(\bm{K}^{-1}+\frac{\tau}{2\zeta
kT}\bm{K}^{-1}\bm{K^{*}}{}^{T}\bm{K^{*}}\bm{K}^{-1}\right)^{-1}\bm{K}^{-1}\bm{K^{*}}{}^{T}\bm{v}_{1:N-1}$
(33) $\displaystyle\tilde{\bm{K}}=$
$\displaystyle\left(\bm{K}^{-1}+\frac{\tau}{2\zeta
kT}\bm{K}^{-1}\bm{K^{*}}{}^{T}\bm{K^{*}}\bm{K}^{-1}\right)^{-1}$ (34)

For multidimensional trajectories, we separate the likelihood into the forces
experienced along each trajectory,

$\displaystyle\mathcal{P}\left(\bm{x}_{2:N}|f(\cdot),x_{1}\right)\propto$
$\displaystyle\prod_{d}^{D}\textbf{Normal}\left(\bm{f}^{d}_{1:N-1};~{}\zeta\bm{v}^{d}_{1:N-1},~{}\frac{2\zeta
kT}{\tau}\bm{I}\right),$ (35)

where $d$ is the dimension index and $D$ is the number of dimensions. A clever
choice of indexing allows us to sample the potential directly even with
multidimensional trajectories. It consists of “flattening out” our
multidimensional force, velocity and kernel arrays into one dimensional
vectors, and “flattening out” the multidimensional kernel matrices into a
matrix.

The clever choice is as follows: Let the first $N$ indices of all variables be
reserved for values concerning the $N$ measurements in the first dimension,
then the $N+1$ through $2N$ indices be reserved for values concerning the $N$
measurements in the second dimension, and continue this pattern until all
$D\times N$ measurements have been accounted for. For example, define $\bm{f}$
to be a vector such that $f_{1}$ is the force at the first time level along
the first dimension, $f_{2}$ is the force at the first time level along the
second dimension, …, $f_{D+1}$ is the force at the first time level along the
second dimension, $f_{D+2}$ is the force at the second time level along the
second dimension, and so forth. By utilizing this pattern for $\bm{f}$,
$\bm{v}$, and $\bm{K}^{*}$ we can sample the potential directly using Eq. 32.

####  V.2.3 Friction coefficient

The conditional probability for the friction coefficient is the product of
Eq.s 3, 10 and 11

$\displaystyle\mathcal{P}\left(\zeta|\bm{x}_{1:N},U(\cdot)\right)\propto$
$\displaystyle\textbf{Gamma}\left(\zeta;\alpha_{\zeta},\beta_{\zeta}\right)\textbf{Normal}\left(x_{1};~{}0,~{}\frac{2\tau
kT}{\zeta}\right)$
$\displaystyle\times\prod_{n=2}^{N}\textbf{Normal}\left(x_{n};~{}x_{n-1}+\frac{\tau}{\zeta}f(x_{n-1}),~{}\frac{2\tau
kT}{\zeta}\right)$ (36)

which cannot be simplified into any known elementary distribution. We
therefore sample, $\zeta$ using a Metopolis Hastings algorithm
[bishop2006pattern].

####  V.2.4 Measurement noise

As our inverse gamma prior for $\sigma^{2}$, Eq. 12 is conjugate to our
likelihood, Eq. 6, our conditional probability for $\sigma^{2}$ can be
simplified to an inverse gamma distribution [gelman2013bayesian,
sivia2006data]

$\displaystyle\mathcal{P}\left(\sigma^{2}|\bm{y}_{1:N},\bm{x}_{1:N}\right)=$
$\displaystyle\textbf{InvGamma}\left(\alpha_{\sigma^{2}}+\frac{N}{2},~{}\beta_{\sigma^{2}}+\frac{1}{2}\sum|\bm{x}_{n}-\bm{y}_{n}|^{2}\right).$
(37)

###  V.3 Boltzmann method

In the Results section, we will compare to the Boltzmann method. The Boltzmann
method, as opposed to the other existing methods, has the advantage that it is
physically intuitive (it is derived from thermodynamics), it does not assume a
potential shape a priori and therefore can be used for non-harmonic
potentials. It’s main limitation, that it does not treat measurement noise, is
a limitation that all other competing methods posses.

The Boltzmann method uses the Boltzmann distribution from thermodynamics,
which relates the potential in a region of space to the fraction of time that
the particle will be seen in that region,

$\displaystyle U_{i}$ $\displaystyle=-kT\log(p_{i})$ (38)

where $U_{i}$ is the potential of region $i$ and $p_{i}$ is the fraction of
the time that the particle spent in region $i$.

The limitation of the Boltzmann method is that when measurement noise is
introduced then the fraction of time that a particle was seen in a region,
$p_{i}$, is not equivalent to the fraction of time that the particle was in
the region. In particular, measurement noise will cause the learned potential
to smear out over the range of the measurement noise (see Fig. 3).

###  V.4 Robustness test using simulated data

In this section we benchmark our method using simulated data and demonstrate
its robustness on a range of model parameters. We start by showing results for
a base set of parameters, then one by one change to parameters to show how the
results vary.

####  V.4.1 Demonstration

…

####  V.4.2 Varying the number of data points

…

####  V.4.3 Varying the noise

…