††footnotetext: *Equal contributions.

# Excess wings and asymmetric relaxation spectra in a facilitated trap model

Camille Scalliet  Department of Applied Mathematics and Theoretical Physics,
University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom
<EMAIL_ADDRESS>Benjamin Guiselin  Laboratoire Charles Coulomb (L2C),
Université de Montpellier, CNRS, 34095 Montpellier, France  Ludovic Berthier
Laboratoire Charles Coulomb (L2C), Université de Montpellier, CNRS, 34095
Montpellier, France  Yusuf Hamied Department of Chemistry, University of
Cambridge, Lensfield Road, Cambridge CB2 1EW, United Kingdom

###### Abstract

In a recent computer study, we have shown that the combination of spatially
heterogeneous dynamics and kinetic facilitation provides a microscopic
explanation for the emergence of excess wings in deeply supercooled liquids.
Motivated by these findings, we construct a minimal empirical model to
describe this physics and introduce dynamic facilitation in the trap model,
which was initially developed to capture the thermally-activated dynamics of
glassy systems. We fully characterise the relaxation dynamics of this
facilitated trap model varying the functional form of energy distributions and
the strength of dynamic facilitation, combining numerical results and analytic
arguments. Dynamic facilitation generically accelerates the relaxation of the
deepest traps, thus making relaxation spectra strongly asymmetric, with an
apparent “excess” signal at high frequencies. For well-chosen values of the
parameters, the obtained spectra mimic experimental results for organic
liquids displaying an excess wing. Overall, our results identify the minimal
physical ingredients needed to describe excess processes in relaxation spectra
of supercooled liquids.

##  I Introduction

The large body of experimental work Ediger  _et al._ (1996); Angell (1995)
characterising structural relaxation in supercooled liquids approaching the
glass transition reveals that the slow molecular motion is broadly distributed
in time Böhmer  _et al._ (1993). Correlation functions (in the time domain) or
relaxation spectra (in the frequency domain) demonstrate the existence of
molecular dynamics occurring over all timescales separating microscopic motion
(in the picosecond range) and the average structural relaxation time
$\tau_{\alpha}$ (about hundreds of seconds at the experimental glass
transition temperature $T_{g}$) Böhmer  _et al._ (1993); Lunkenheimer  _et
al._ (2002); Roland  _et al._ (2005). A central theme in glass transition
studies concerns the physical origin of the rapid growth of $\tau_{\alpha}(T)$
near $T_{g}$. However, any physical explanation should also be able to
rationalise the main features of the broad distribution of relaxation times
characterising the dynamics Ediger (2000); Tarjus (2011); Berthier and Biroli
(2011).

Experimentally, various techniques are able to probe the orientational
dynamics of molecules over a broad frequency range, such as dielectric
spectroscopy Lunkenheimer  _et al._ (2000), magnetic resonance Flämig  _et
al._ (2020), and dynamic light scattering Schmidtke  _et al._ (2013). In most
molecular liquids, the structural relaxation is described in the time domain
by a stretched (rather than simple) exponential form. In the Fourier domain,
this becomes a peak that is broader than a simple Lorentzian, and various
empirical functions are used Böhmer  _et al._ (1993). In a large number of
liquids, an “excess” relaxation signal appears in addition to the main
$\alpha$-peak on its high-frequency flank Johari and Goldstein (1970); Dixon
_et al._ (1990); Menon and Nagel (1995); Menon  _et al._ (1992); Leheny and
Nagel (1997, 1998); Adichtchev  _et al._ (2003); Blochowicz  _et al._ (2003);
Gainaru  _et al._ (2009); Caporaletti  _et al._ (2021); Blochowicz  _et al._
(2006); Gainaru (2019). It is then customary to fit this signal as an
additional contribution, or “process”, to the dynamics. Such secondary
relaxation may also take different functional forms depending on the studied
liquid Ngai and Paluch (2004); Wu (1991), but quite often it appears as an
extremely broad peak (in the frequency domain), or even a pure power law
covering many decades Menon and Nagel (1995). This excess signal then
resembles a “wing” in a log-log representation of the spectra, which are then
very asymmetric. The presence of such signal is observed across a wide variety
of materials comprising metallic, organic, ionic, polymeric and atomic
glasses. The emergence of excess wings is therefore a reasonably universal
feature of deeply supercooled liquids.

There are several important questions related to secondary processes in
supercooled liquids. Given that the strength of the signal is much smaller
than the main $\alpha$-peak (typically 1 %), two generic explanations may be
considered. Either all molecules contribute by performing motions that are
about 100 times smaller than during structural relaxation itself Vogel  _et
al._ (2002), or instead a small fraction of about 1 % of the molecules undergo
structural relaxation much before the bulk. A third alternative is that
secondary relaxation is caused by a slow degree of freedom (say, some intra-
molecular motion) that is distinct from the ones contributing to the
$\alpha$-peak, which then raises further questions about the inter-relation
between these processes. Previous work also tried to identify how universal
the shape of the excess signal can be depending on experimental probes Petzold
_et al._ (2013); Körber  _et al._ (2020), the type of molecules Ngai and
Paluch (2004); Buchenau  _et al._ (2020), and thermodynamic conditions (by,
say, varying pressure Roland  _et al._ (2005)).

A variety of physical interpretations and empirical models have been proposed
to explain the asymmetric shape of relaxation spectra. In the simplest
empirical descriptions, the structural and secondary relaxations are described
by the linear superposition of two independent contributions Kudlik  _et al._
(1998); Gainaru  _et al._ (2010). This is of course a valid mathematical
option, but this does not address the nature of the two relaxation processes,
and even implicitly suggests that they can be considered as stemming from
independent and distinct molecular motion, which is a questionable hypothesis.

In most models proposed to describe the glass transition of supercooled
liquids, secondary processes do not necessarily appear at the level of the
basic formulation Berthier and Biroli (2011). Still, there now exists a good
number of proposals regarding the microscopic origin and physical
interpretation of excess wings. In a broad, first category of models,
asymmetric relaxation spectra emerge when the “main” glassy degree of freedom
responsible for the $\alpha$-relaxation is coupled to an “additional” degree
of freedom. In Refs. Diezemann  _et al._ (1999); Mohanty  _et al._ (2000), for
instance, a trap model describing the structural relaxation is passively
coupled to an independent relaxation process, possibly representing a local
molecular process. In mode-coupling approaches, the translational degrees of
freedom responsible for dynamic arrest are also passively coupled to a second
glassy process (for instance mimicking the rotation of molecules) to produce
more complex spectra Götze and Sjögren (1989); Götze and Sperl (2002);
Domschke  _et al._ (2011); Chong and Götze (2002); Cummins (2005); Sjögren
(1986). In a very similar spirit, the coupling of two schematic glass models
was shown to produce spectra with several slow processes Crisanti  _et al._
(2011); Crisanti and Leuzzi (2015), a notable difference with mode-coupling
approaches being the more reciprocal interaction between the two glassy
processes.

In a second broad family of approaches, asymmetric spectra are directly
obtained from an underlying broad distribution of relaxation times with an
asymmetric shape. This broad distribution is meant to describe the
heterogeneous nature of glassy relaxation, which is ascribed to static Sethna
_et al._ (1991); Stevenson and Wolynes (2010); Viot  _et al._ (2000);
Chamberlin (1999); Dyre and Schrøder (2000) or kinetic Berthier and Garrahan
(2005) underlying mechanisms. To account for the winged asymmetric shapes of
spectra in such approaches, one needs to invoke at least one physical
ingredient. In the approach based on geometric frustration Tarjus  _et al._
(2005), large correlated domains are more energetically suppressed than small
ones due to frustration while they all relax via thermal activation, which
accounts for excess wings Viot  _et al._ (2000). In the model developed by
Chamberlin Chamberlin (1993, 1998, 1999) the distribution of domain sizes is
symmetric, but large and small domains relax differently Chamberlin (1999),
thus resulting in asymmetric time distributions. In some dynamically
facilitated models Garrahan  _et al._ (2011); Chandler and Garrahan (2010),
the structural relaxation emerges as a hierarchical process involving a
distribution of timescales which is skewed towards small times, giving rise to
asymmetric spectra in the Fourier domain Berthier and Garrahan (2005). In the
mosaic (droplet) picture of the random first order transition theory
Kirkpatrick  _et al._ (1989); Lubchenko and Wolynes (2007), secondary
processes result from fluctuations of the droplet sizes and shapes Stevenson
and Wolynes (2010). They thus induce an asymmetric distribution of free energy
barriers overcome by thermal activation. We note that the impact of dynamic
facilitation on the dynamics of the largest droplets was also invoked in this
context Xia and Wolynes (2001); Bhattacharyya  _et al._ (2008); Biroli and
Bouchaud (2012); Berthier  _et al._ (2019a), but this was not related to the
existence of excess wings.

In a recent computational study Guiselin  _et al._ (2021), we investigated the
equilibrium relaxation dynamics of a simple supercooled liquid thermalised
down to the experimental glass transition temperature $T_{g}$ using the swap
Monte Carlo algorithm Berthier  _et al._ (2016); Ninarello  _et al._ (2017);
Berthier  _et al._ (2019b), and we were able to follow the relaxation dynamics
of the particles over a time window of about 10 orders of magnitude. In
particular, we could access the temperature regime where excess wings are
observed experimentally and record the dynamics over the relevant time window
as well. We demonstrated the emergence of an excess wing in this deeply
supercooled regime, and observed the microscopic particle motion responsible
for the asymmetric shape of the relaxation spectra. Our main conclusion was
that relaxation close to $T_{g}$ is the result of two central physical
ingredients, dynamic heterogeneity and dynamic facilitation. First, relaxation
is initiated over timescales much shorter than $\tau_{\alpha}$ at a sparse
population of localised regions, these relaxation events being extremely
broadly distributed. The dynamics is therefore spatially and temporally
heterogeneous in this regime. These sparse relaxed regions were found to
trigger the relaxation of neighbouring regions by kinetic facilitation. We
argued that since dynamic facilitation accelerates the relaxation of the
slowest regions in the liquid, the broad relaxation spectrum observed at high
frequencies becomes compressed at low frequencies, thus explaining the
asymmetric shape of the observed spectra.

In our previous work Guiselin  _et al._ (2021), we suggested that a simple
facilitated trap model combining these two ingredients (heterogeneity and
facilitation) should generically lead to asymmetric relaxation spectra. Here,
we provide a complete study of the relaxation dynamics of the proposed model.
We describe how this facilitated trap model is actually constructed, by making
the simplest possible assumptions. We then extensively explore the parameter
space. We show that, despite its extreme simplicity, the model naturally and
generically gives rise to asymmetric, winged relaxation spectra in the
presence of dynamic facilitation. We also propose a microscopic analysis of
the model itself, combining detailed numerical simulations to an approximate
analytic solution. Combined to our numerical simulations Guiselin  _et al._
(2021), these results provide a simple, physically-motivated microscopic
picture for the shape of relaxation spectra in deeply supercooled liquids.

Our article is organised as follows. In Sec. II, we introduce the facilitated
trap model and dynamic observables. In Sec. III, we obtain the relaxation
spectra. We present in Sec. IV our analysis of the individual trap dynamics to
rationalise the spectral shapes. In Sec. V we summarise our results and
discuss the microscopic picture that explains the asymmetric, winged
relaxation spectra of deeply supercooled liquids which emerges from our study,
in relation to both experiments and previous models.

##  II Facilitated trap model

In this section we define the facilitated trap model, along with observables
employed to investigate its dynamics. We also describe our numerical
simulations of the model.

###  II.1 Original trap model

The dynamics of deeply supercooled liquids exhibits very strong dynamic
heterogeneities Ediger (2000); Berthier  _et al._ (2011); Berthier (2011).
This suggests that the supercooled liquid can be coarse-grained into
independently relaxing domains characterised by a local relaxation time that
is broadly distributed, as illustrated in Fig. 1(a). This spatially
heterogeneous viewpoint is mathematically captured by the original trap model
Brawer (1984); Dyre (1987); Bässler (1987); Bouchaud (1992); Monthus and
Bouchaud (1996); Jean-Philippe Bouchaud  _et al._ (1995); Denny  _et al._
(2003); Heuer  _et al._ (2005); Diezemann and Heuer (2011). Our computer
simulations of a simple atomistic glass-forming model has revealed the
heterogeneous nature of the particle motion that gives rise to the high-
frequency part of the relaxation spectra Guiselin  _et al._ (2021). The trap
model is therefore a natural starting point.

Figure 1: Facilitated trap model. (a) The liquid is coarse-grained as a
collection of traps of depth $E$ distributed according to $\rho(E)$. Hopping
of a mesoscopic domain is a thermally-activated process. (b) Whenever a trap
relaxes, it perturbs the energy of all other traps by an amount proportional
to $\Delta$.

The original trap model describes the thermally-activated dynamics of $N$
mesoscopic domains hopping between energy traps of depth $E\geq 0$, with or
without any spatial structure. This description makes it clear that a direct
connection to detailed molecular motion in the original system is difficult.
At a temperature $T$, the domains escape the trap they occupy after a Poisson-
distributed relaxation time $\tau$ of mean $\tau_{0}e^{\beta E}$, with
$\beta=1/T$ the inverse temperature (the Boltzmann constant is set to unity)
and $\tau_{0}$ a microscopic time constant, as illustrated in Fig. 1(b).
Dynamic heterogeneity is introduced through the distribution $\rho(E)$ of trap
depths. We consider a general functional form for this distribution

$\rho(E)=\frac{\alpha}{\Gamma(1/\alpha)E_{0}}e^{-(E/E_{0})^{\alpha}},$ (1)

with $E_{0}$ a constant energy scale, and an exponent $1<\alpha\leq 2$ to
smoothly interpolate between the exponential Bouchaud (1992); Monthus and
Bouchaud (1996); Jean-Philippe Bouchaud  _et al._ (1995) and Gaussian Brawer
(1984); Dyre (1987); Bässler (1987); Diezemann (1997); Rehwald  _et al._
(2010) distributions. Both these models have been widely studied in various
contexts. Bässler and coworkers Arkhipov and Baessler (1994); Arkhipov  _et
al._ (1996) have previously introduced the idea of an exponent $\alpha$ which
varies continuously, in order to perform quantitative comparisons to some
experimental results.

When a given domain relaxes and jumps out of its trap, it chooses a new trap
depth $E^{\prime}$ with probability $\rho(E^{\prime})$ without memory of the
previous energy $E$. The resulting dynamics of the model thus depends on two
parameters: $\alpha$, which specifies the distribution $\rho(E)$, and the
temperature $T$. We express energies and times in units of $E_{0}$ and
$\tau_{0}$, respectively, and consider a mean-field version with no spatial
structure or interaction between traps.

In thermal equilibrium, the probability distribution $P_{\text{eq}}(E)$ for a
sub-system to be in a trap of depth $E$ at the temperature $T$ is given by
Monthus and Bouchaud (1996)

$P_{\text{eq}}(E)=\frac{\rho(E)e^{\beta E}}{Z(\beta)},$ (2)

with $Z(\beta)=\int_{0}^{\infty}\rho(E)e^{\beta E}\mathrm{d}E$ a normalisation
factor which is finite and strictly positive at all $T>0$ as long as
$\alpha>1$ (the exponential trap model with $\alpha=1$ has instead a finite
temperature phase transition to a non-ergodic low-temperature phase Bouchaud
(1992)). We define the average trap depth $\overline{E}(T)$ via
$\overline{E}=\int_{0}^{\infty}P_{\text{eq}}(E)E\mathrm{d}E$.

By considering the thermally-activated relaxation over distributed energy
barriers, the trap model naturally gives rise to dynamic heterogeneities. Sub-
systems in shallow traps $E\ll\overline{E}$ relax much faster than the ones
stuck in deep traps with $E\gg\overline{E}$. The parameter $\alpha$ directly
controls the breadth of dynamic heterogeneities: smaller $\alpha$ corresponds
to broader distribution $\rho(E)$, and thus to much wider distribution of
local relaxation times.

###  II.2 Facilitated trap model

The above model considers that sub-systems are dynamically independent, and
that the motion out of a given trap does not affect the state of the other
domains. This is not true in realistic glass-formers, in which a relaxation
event corresponds to a local rearrangement of particles, which inevitably
modifies the environment of the neighbouring particles Chandler and Garrahan
(2010). Computer simulations reveal that this effect facilitates the gradual
spreading of relaxation events through the liquid Bergroth  _et al._ (2005);
Vogel and Glotzer (2004); Candelier  _et al._ (2009, 2010); Keys  _et al._
(2011); Schoenholz  _et al._ (2016), and plays an increasingly important role
at low temperatures Chacko  _et al._ (2021) (an opposite conclusion was
reached in Ref. Candelier  _et al._ (2009)). In our numerical investigation of
excess wings, we also found that the initial relaxation in isolated locations
facilitates the relaxation of nearby regions Guiselin  _et al._ (2021).

While most studies of the trap model ignore the effect of dynamic
facilitation, we note that the initial studies of the exponential trap model
in Refs. Jean-Philippe Bouchaud  _et al._ (1995); Monthus and Bouchaud (1996)
had actually introduced such dynamic correlation between traps. A Gaussian
trap model with dynamic exchanges between fast and slow regions was studied in
Ref. Diezemann (2005) to analyse the consequences of spatially heterogeneous
dynamics in supercooled liquids. More recently, Heuer and coworkers have also
included a similar effect in a lattice version of the model, in order to
describe the finite-size dependence of the dynamics of supercooled liquids
Rehwald  _et al._ (2010); Rehwald and Heuer (2012). Here we follow a similar
philosophy and also assume that every hopping event induces a small shift in
the depth of the neighbouring traps, as illustrated in Fig. 1(b).

To simplify the analytic description of this facilitation effect we assume,
again in a mean-field spirit, that all traps are equally affected by hopping
processes. In practice, whenever a trap relaxes, we attempt to shift the
energy of all other traps by a random amount $\delta E$, different for each
trap, uniformly distributed in the interval
$\left[-\frac{\Delta}{2\sqrt{N}},\frac{\Delta}{2\sqrt{N}}\right]$:
$E\rightarrow E^{\prime}=E+\delta E$. The energy scale $\Delta$ controls the
strength of dynamic facilitation ($\Delta=0$ gives back the original trap
model) and the scaling with $N$ ensures that the resulting dynamics scales
correctly with the system size. The change in energy is then accepted or
rejected in order to leave the equilibrium probability distribution
$P_{\text{eq}}(E)$ unchanged. This facilitation effect thus corresponds to a
drift-diffusion process of the trap depth in energy space in an effective
confining potential (in appropriate units) $V_{\rm eff}=-\log P_{\mathrm{eq}}$
(see Sec. IV.2 for an analytic description). If the new trap energy
$E^{\prime}$ is accepted, we assign a new Poisson-distributed relaxation time
of mean $e^{\beta E^{\prime}}$ to that site. At fixed $\alpha$ and $T$,
increasing $\Delta$ corresponds to increasing the strength of dynamic
facilitation. Notice that due to dynamic facilitation, the energy of a given
trap can both increase or decrease in a stochastic manner. However, since the
energy diffusion occurs in a potential $V_{\rm eff}$, the energy of traps with
a large initial depth will slowly return back to the average energy value. We
shall show below that such dynamic facilitation leads to an average speedup of
the global dynamics, which is not introduced by hand in the formulation of the
model.

###  II.3 Dynamic observables

By construction, the equilibrium distribution of energies $P_{\text{eq}}(E)$
does not depend on the presence of dynamic facilitation and is independent of
$\Delta$. Although thermodynamic properties do not depend on $\Delta$ we
expect the dynamic relaxation to be strongly affected by the presence of
dynamic facilitation.

To investigate this effect, we define a time correlation function to describe
the dynamics of the traps. To this end, we introduce the mean persistence
$p(t)$, which quantifies the fraction of sub-systems which, at time $t$, have
not escaped the trap they occupied at time 0. The persistence function
monotonically decreases from one at $t=0$ to zero at long times, when all sub-
systems have relaxed. It is physically equivalent to the self-intermediate
scattering function at large wavevectors, which is more typically used in off-
lattice glass models. The mean persistence $p(t)$ can be expressed as an
average of the persistence $p(t;E_{i})$ of a single sub-system starting from
an initial trap depth $E_{i}$ at $t=0$ over the equilibrium distribution
$P_{\mathrm{eq}}(E_{i})$,

$p(t)=\int_{0}^{\infty}\mathrm{P}_{\text{eq}}(E_{i})p(t;E_{i})\mathrm{d}E_{i}.$
(3)

In the absence of dynamic facilitation, the persistence of a single sub-system
can be derived exactly from the exponential distribution of relaxation times
and reads $p(t;E_{i})=e^{-te^{-\beta E_{i}}}$. The global persistence can thus
be computed analytically from Eq. (3). In the presence of dynamic facilitation
$\Delta>0$, the persistence of a single sub-system $p(t;E_{i})$ stems from a
non-trivial interplay between the activated jump dynamics and the facilitated,
diffusive, dynamics of the trap depth in energy space, and cannot be computed
analytically exactly. In Sec. IV.2 we provide an approximate analytic solution
for $p(t)$ in the presence of facilitation. The exact persistence of the
facilitated trap model is thus measured by means of numerical simulations
using

$p(t)=\left\langle\frac{1}{N}\sum_{i=1}^{N}p_{i}(t)\right\rangle,$ (4)

where $p_{i}(t)=1$ if the sub-system $i$ has not escaped at time $t$ from the
trap it occupied at time 0, and $p_{i}(t)=0$ otherwise. The brackets indicate
an average over independent initial conditions (trap energies) drawn from the
equilibrium distribution $P_{\rm eq}(E)$.

In experiments, glassy features are often investigated thanks to spectroscopic
techniques. In order to allow for a direct comparison, we introduce a dynamic
susceptibility $\chi^{\prime\prime}(\omega)$ in the frequency domain.
Following the experimental literature, we assume that the dynamics can be
described by a distribution of relaxation times $G(\log\tau)$ such that Böhmer
_et al._ (1998); Blochowicz  _et al._ (2003); Berthier and Garrahan (2005)

$\chi^{\prime\prime}(\omega)=\int_{-\infty}^{\infty}G(\log\tau)\frac{\omega\tau}{1+(\omega\tau)^{2}}\mathrm{d}\log\tau,$
(5)

with $\omega$ the frequency (up to a factor $2\pi$), and where the double
prime indicates that we consider the imaginary part of the total
susceptibility. The global persistence corresponds to the superposition of
step functions decaying from 1 to 0 when a sub-system jumps. As a consequence,
the relation $G=-\mathrm{d}p/\mathrm{d}\log t$ holds Berthier and Garrahan
(2005) and the susceptibility spectrum $\chi^{\prime\prime}(\omega)$ can be
directly computed from the persistence data:

$\chi^{\prime\prime}(\omega)=-\int_{-\infty}^{+\infty}\frac{\mathrm{d}p}{\mathrm{d}\log
t}(\log\tau)\frac{\omega\tau}{1+(\omega\tau)^{2}}\mathrm{d}\log\tau,$ (6)

for any value $\Delta\geq 0$. In the case $\Delta=0$, we directly evaluate the
persistence given by Eq. (3). For $\Delta>0$, we instead use the persistence
measured in the simulations via Eq. (4). Both are then used to compute
$\chi^{\prime\prime}(\omega)$ via Eq. (6).

###  II.4 Simulations of the facilitated trap model

When $\Delta>0$, the persistence cannot be computed analytically exactly from
Eq. (3) and we instead use simulations and Eq. (4) to measure the time
correlation function. We provide the details of these simulations.

First we initialize the $N$ trap energies by directly sampling the equilibrium
distribution $P_{\text{eq}}(E)$. To do so, we numerically evaluate the
cumulative probability distribution of energies
$\mathcal{C}_{\text{eq}}(E)=\int_{0}^{E}P_{\mathrm{eq}}(E^{\prime})\mathrm{d}E^{\prime}$
on a grid of a thousand points from $E_{\mathrm{min}}$ where
$\mathcal{C}_{\mathrm{eq}}(E_{\mathrm{min}})=10^{-10}$, to $E_{\mathrm{max}}$
where $\mathcal{C}_{\mathrm{eq}}(E_{\mathrm{max}})=1-10^{-10}$, as ${\cal
C}_{\rm eq}$ cannot be computed analytically for arbitrary values of $\alpha$.
We use a cubic spline interpolation to construct numerically the reciprocal
function $\mathcal{E}={{\mathcal{C}}_{\rm eq}}^{-1}$. For each trap, we
generate a random variable $X$ from a uniform distribution in the range
$\left[0,1\right]$, and we assign it an initial energy $E_{i}=\mathcal{E}(X)$.

The relaxation time of a trap of depth $E$ is drawn from an exponential
distribution of mean $e^{\beta E}$. To accomplish this, we generate a random
number $X$ from a flat distribution in the range $\left[0,1\right]$ and
construct the relaxation time $\tau=-e^{\beta E}\log(X)$.

In the course of the simulation, when a trap relaxes, we sample the
probability distribution $\rho(E)$ defined in Eq. (1) to select its new
energy. For this purpose, we generate two random variables: $X$ from a flat
distribution in the range $\left[0,1\right]$ and $Y$ from the Gamma
distribution $\Gamma(1+1/\alpha,1)$ with shape parameter 1 + $1/\alpha$ and
rate parameter 1. The resulting energy $E=XY^{1/\alpha}$ is then distributed
according to $\rho(E)$ Marsaglia and Tsang (2000).

Facilitation in the model corresponds to a diffusion process in energy space
that preserves the equilibrium distribution of trap energies, hence which is
confined to an effective potential $V_{\rm eff}=-\log P_{\text{eq}}$ (in
appropriate units). To this end, we compute the difference in effective
potential that a proposed shift of energy would imply, $\delta V_{\rm
eff}=\log(P_{\rm eq}(E)/P_{\rm eq}(E^{\prime}))$ and use a Metropolis filter.
If $\delta V_{\rm eff}<0$, we accept the new trap energy $E^{\prime}$. If
$\delta V_{\rm eff}>0$, the change is accepted with probability $e^{-\delta
V_{\rm eff}}$.

To measure the global persistence given by Eq. (4), we perform numerical
simulations for $N=100$ until full decorrelation is reached, using $\Delta$
values in the range $\left[10^{-4},1\right]$. Since the persistence decreases
by steps of size $1/N$, large systems of size $N=10^{4}$ are used to resolve
the persistence function at very short times, when needed. We then combine the
persistence data measured at short times in very large systems with the data
measured in smaller systems at longer times. We typically average the
persistence over 100 independent simulations for $N=100$, and 10 runs for
$N=10^{4}$ in order to obtain sufficient statistics.

##  III Relaxation spectra

In this section, we first report the relaxation spectra of the trap model in
the absence of dynamic facilitation for $\Delta=0$. We present results for
various underlying distributions of energy barriers $\rho(E)$, parametrized by
$\alpha$, and for different temperatures $T$. We then present the effect of
dynamic facilitation, characterised by $\Delta$, on the relaxation spectra.

###  III.1 Original trap model without facilitation

Figure 2: Relaxation spectra of the original trap model in the absence of
facilitation for different distributions of trap energies $\rho(E)\sim
e^{-E^{\alpha}}$: (a) $\alpha=2$, (b) $\alpha=1.5$ and (c) $\alpha=1.1$, and
various temperatures $T=1/\beta$. For a given value of $\alpha$, the spectra
broaden and shift to lower frequencies with decreasing temperature. In
addition, the spectra also broaden moving from a Gaussian (a) to an almost-
exponential (c) distribution $\rho(E)$ of trap energies.

We consider the relaxation spectra of the trap model in the absence of dynamic
facilitation ($\Delta=0$) to illustrate the variety of spectral shapes
exhibited by the model.

We present in Fig. 2 the relaxation spectra obtained for $\alpha=2$, 1.5, and
1.1 at various temperatures. The former corresponds to the Gaussian model,
which has been studied in the past, while the latter is very close to an
exponential distribution for $\rho(E)$. In order to highlight quantitative
differences, we show the data over the same frequency and susceptibility
ranges in all cases. We also restrict the frequencies to an experimentally-
relevant regime, extending over about 14 decades. We observe that all spectra,
at various $\alpha$ and temperatures $T$, share a qualitatively similar shape:
they are broad, spanning many orders of magnitude in the frequency domain, and
are relatively symmetric around their maximum value. The broadness of the
spectra is a direct signature of dynamic heterogeneity, and reflects the fact
that the sub-systems relax over a broad distribution of timescales.

Yet, we notice from Fig. 2 a clear quantitative evolution of the spectra as a
function of $\alpha$ and $\beta=1/T$. At fixed $\alpha$, the spectra shift to
lower frequencies and broaden as the temperature decreases. This is due to
thermal activation, which gives rise to longer relaxation times and
simultaneously enhances dynamic heterogeneity at low temperatures. The
exponent $\alpha$ directly controls how broad the underlying distribution of
trap energies $\rho(E)$ is, and therefore how strong dynamic heterogeneities
are, irrespective of the temperature. While a Gaussian distribution $\rho(E)$
for $\alpha=2$ yields relatively narrow spectra, see Fig. 2(a), the almost-
exponential distribution $\rho(E)$ for $\alpha=1.1$ generates the extremely
broad spectra shown in Fig. 2(c), which extend to physically unreachable low
frequencies at low temperatures, for example for $\beta=1.59$, where the
entire spectrum does not fit the experimental frequency window.

These results confirm the salient and well-known features of the trap model
that generically leads to slow dynamics and broad relaxation spectra, which
appear, however relatively featureless with a main relaxation peak in the
Fourier domain. Notice that the spectra found within the trap model are very
broad (the behaviour at low frequencies is typically slower than a power law
and depends on both $\alpha$ and $T$) and temperature dependent (so that time
temperature superposition is not obeyed).

###  III.2 Effect of dynamic facilitation

We now investigate how the introduction of facilitation, parametrized by
$\Delta>0$, affects the spectral shapes. We have systematically investigated
the effect of $\Delta$ on the dynamics for $\alpha\in]1,2]$ and various
temperatures. We present results for well-chosen values of $\alpha$, $\beta$,
and varying $\Delta$, to illustrate the variety of observed behaviours.

We present in Fig. 3 the relaxation spectra of the facilitated trap model for
$\alpha=2$, $1.5$, and $1.1$. For each $\alpha$, we choose a relatively low
temperature, for which the spectrum without facilitation is very broad and
thus for which facilitation is expected to have a larger impact. In Fig. 3(a),
we show the spectra for the Gaussian model ($\alpha=2$) at an inverse
temperature $\beta=8$. The spectrum without facilitation $\Delta=0$ (solid
line) is broad, relatively symmetric, and reaches a maximum at
$\log_{10}\omega_{\alpha}\sim-14$. We show the data for $\Delta$ increasing
logarithmically from $0.001$ to $1$ (left to right). At the lowest
$\Delta=0.001$, the data coincides with $\Delta=0$ at high frequencies, but
deviates from it below $\log_{10}\omega\sim-12$. We see that low frequencies
are suppressed by $\Delta>0$, causing a compression of the spectrum into a
peak around $\log_{10}\omega_{p}\sim-14$. The resulting peak is much sharper
than the broad hump of the underlying $\Delta=0$ curve. As $\Delta$ increases,
the low-frequencies are increasingly suppressed, and the resulting peak
increasingly sharper, shifting to a higher $\omega_{p}$. The frequency range
over which the data coincide with $\Delta=0$ decreases with increasing
$\Delta$.

Figure 3: Effect of dynamic facilitation, parametrized by $\Delta$, on the
relaxation spectra with parameters: (a) $\alpha=2$, $\beta=8$, (b)
$\alpha=1.5$, $\beta=4$, (c) $\alpha=1.1$, $\beta=1.55$, (d) $\alpha=1.5$,
$\beta=3.5$.

The suppression of low frequencies and the emergence of a sharper peak is a
generic effect of imposing $\Delta>0$ in the facilitated trap model. This is
confirmed by the results in Figs. 3(b,c) obtained with $\alpha=1.5$, $\beta=4$
and $\alpha=1.1$, $\beta=1.55$, respectively. They show a similar trend to the
Gaussian case in Fig. 3(a). In the three first panels, facilitation gives rise
to a sharp peak at $\omega_{p}>\omega_{\alpha}$ (the latter being the location
of the maximum in the $\Delta=0$ spectrum). This peak shifts to higher
frequencies as the strength of facilitation is raised. On the other hand,
spectra follow the underlying $\Delta=0$ signal at high frequencies, resulting
in a highly asymmetric spectrum, in great contrast with the symmetric spectra
without facilitation shown in Fig. 2. Finally, we note that the effect of
facilitation is more spectacular at small $\alpha$, where the $\Delta=0$
spectrum is extremely broad, extending to unphysically small frequencies. The
relative amplitude of the peak is almost two orders of magnitude larger than
the underlying signal at $\Delta=0$ in this case, while it represents about
one order of magnitude for $\alpha=1.5$, and even less for $\alpha=2$.

Finally, we show in Fig. 3(d) the spectra obtained for $\alpha=1.5$,
$\beta=3.5$. The exponent $\alpha$ is the same as in panel (b), but the
temperature is higher, explaining why the $\Delta=0$ spectrum is centered
around a larger frequency $\log_{10}\omega_{\alpha}\sim-8$. Here, a small
$\Delta=0.001$ also compresses the low-frequency part of the spectrum, giving
rise to a secondary peak but the effect is however extremely weak. As $\Delta$
increases, the peak clearly emerges, and shifts to higher frequencies as for
the lower temperature. The comparison between panels (b) and (d) reveals that
for a given value of $\Delta$, the effect of facilitation is more pronounced
at lower temperature and depends on both $\Delta$ and $T$. This point is
further addressed in Sec. IV.2.

Overall, this analysis demonstrates that dynamic facilitation compresses the
low-frequency part of the relaxation spectra, but leaves the high-frequency
regime unaffected. This generically gives rise to highly asymmetric shapes
composed of a relatively sharp peak at a low frequency $\omega_{p}$ mainly
controlled by $\Delta$, and a broader shape on the high-frequency flank,
controlled by the underlying distribution $P_{\rm eq}(E)$, and thus by the
exponent $\alpha$. Whereas the high-frequency regime of the spectra shown in
Fig. 3 depends on the control parameters, the low-frequency regime is
identical with a linear frequency dependence, expected when the underlying
distribution of relaxation times is bounded at large times. The linear
frequency dependence is not obeyed in the original trap model for $\Delta=0$.
We have not systematically investigated whether time temperature superposition
holds for the facilitated trap model, but we recall that some models with
kinetic facilitation do obey superposition while some others do not, the
difference lying in the details of the kinetic facilitation rules Garrahan
_et al._ (2011).

##  IV Microscopic analysis of the facilitated trap model

In this section we investigate the dynamics of individual traps in the
presence of dynamic facilitation to rationalise the evolution of the
relaxation spectra characterised by an emerging peak frequency $\omega_{p}$
and an asymmetric winged shape.

###  IV.1 Origin of asymmetric spectra

Figure 4: Representative time series of a given trap depth $E(t)$ in the
absence (a) and presence (b-d) of dynamic facilitation. The vertical axis is
reversed so that large $E$ values appear as “deep” traps. While the time
series (a, b) look similar, the inset in (b) reveals the slow diffusion of the
energy in deep traps. (c,d) Time evolution of a trap starting with an energy
$E_{i}$ at $t=0$ in the presence of facilitation. The effect of $\Delta>0$ for
$E_{i}=6$ is little (c), but accelerates greatly the dynamics of a deep trap
with $E_{i}=12$ (d), which undergoes its next jump when its energy is $E_{\rm
jump}$ at a time $\tau_{\rm jump}\ll e^{\beta E_{i}}$. Data are obtained for
$\alpha=1.5$, $\beta=4$, $\Delta=0.01$.

We illustrate the dynamics of a single trap with and without dynamic
facilitation in Fig. 4, for $\alpha=1.5$, $\beta=4$. In the original trap
model, see Fig. 4(a), the time series is a succession of plateaus
corresponding to the residence times at each given depth. As expected, shallow
traps are short-lived, while deep ones are much longer-lived.

When $\Delta>0$, the trap energy constantly receives energy kicks when other
traps hop. The time series is therefore composed of sharp jumps when the
energy of the probed trap is renewed, but its energy now slowly evolves in
between the jumps. The inset in Fig. 4(b) magnifies these small fluctuations.
The effect of a small $\Delta>0$ on the dynamics of shallow traps is
negligible.

We display in Fig. 4(c) the evolution of the energy for a trap starting from a
large initial depth $E_{i}$ for the same parameters. In particular the value
$E_{i}\sim 12$ is larger than the mean energy value $\overline{E}\sim 7.2$.
For $\Delta=0$, the energy would remain constant over a typical time $e^{\beta
E_{i}}\sim 10^{20}$. When $\Delta=0.01$, the energy slowly drifts towards
lower values, closer to the mean. From this time evolution, we define the
energy $E_{\text{jump}}$ of the trap when it jumps, along with the time
$\tau_{\text{jump}}$ before it jumps, see Fig. 4(c). Remarkably, we have that
$\tau_{\text{jump}}\ll e^{\beta E_{i}}$. Therefore, adding facilitation
dramatically accelerates the relaxation of deep traps. In this specific
example, we also note that $\tau_{\text{jump}}\neq e^{\beta E_{\text{jump}}}$.
In general we find that for large $E_{i}$, we have $\tau_{\text{jump}}\gg
e^{\beta E_{\text{jump}}}(\sim 10^{8})$. This suggests that the relaxation
mechanism of the deepest traps is ruled by the diffusion of the trap depth in
energy space: the trap energy slowly diffuses until it explores low values
which allow for its relaxation.

We present in Fig. 5 scatter plots of the energy $E_{\text{jump}}$ and the
time $\tau_{\text{jump}}$ at which a trap of initial energy $E_{i}$ jumps. The
data shown in Fig. 5 are obtained for $\alpha=1.5$ and $\beta=3.5$ and various
$\Delta$ values. The initial energies $E_{i}$ were drawn from the equilibrium
distribution $P_{\mathrm{eq}}(E)$ which is centered, for these parameters,
around $\overline{E}\sim 5.5$, and becomes smaller than $10^{-2}$ for $E_{i}$
values above 10. Starting with energies, we observe that for $\Delta=0$, the
data trivially collapse on the $E_{\text{jump}}=E_{i}$ line. For very small
$\Delta=10^{-4}$, the data points follow the $E_{\text{jump}}=E_{i}$ line at
low $E_{i}<6$, and depart from it above. For the majority of large $E_{i}$, we
see that $E_{\text{jump}}<E_{i}$, confirming the trend of Fig. 4. With
increasing $\Delta$, the data depart from the $E_{\text{jump}}=E_{i}$ line at
smaller energies, and saturate to smaller $E_{\text{jump}}$ values at large
$E_{i}$.

We present in Fig. 5(b) the scatter plot of the escape time
$\tau_{\text{jump}}$ from a trap of initial energy $E_{i}$. For $\Delta=0$,
the relaxation time $\tau_{\text{jump}}$ for a given $E_{i}$ is Poisson-
distributed with mean $e^{\beta E_{i}}$. As $\Delta$ increases, we see that
large $\tau_{\text{jump}}$ are suppressed, with a decreasing cutoff value,
mimicking the behaviour of $E_{\rm jump}$. The suppression of large escape
times is better appreciated in the distributions shown in the inset in Fig.
5(b).

Figure 5: Scatter plot of (a) the energy $E_{\text{jump}}$ and (b) time
$\tau_{\text{jump}}$ at which traps of initial energy $E_{i}$ jump. Data are
for $\alpha=1.5$, $\beta=3.5$, and from top to bottom $\Delta=0,\ 10^{-4},\
10^{-3},\ 10^{-2},\ 10^{-1}$. The traps relax at smaller energies and
timescales with increasing $\Delta$, on average. The inset in (b) collects the
data for $\tau_{\rm jump}$ in the distribution $P(\log_{10}(\tau_{\rm
jump}))$, which confirms that large relaxation times are suppressed by dynamic
facilitation.

The observations in Figs. 4 and 5 qualitatively explain the shape of the
relaxation spectra presented in Fig. 3. The dynamics of shallow traps is
essentially unaffected by facilitation. On the contrary, the relaxation of
deep traps is greatly accelerated with $\Delta>0$, as illustrated in Fig.
4(c), which effectively suppresses the long relaxation times. Facilitation
thus leaves unchanged the distribution of relaxation times $G(\log\tau)$ at
small $\tau$ but results in a compression and a rather sharp cutoff at large
$\tau$ [see the inset in Fig. 5(b)]. At fixed frequency $\omega$, the integral
in Eq. (5) providing $\chi^{\prime\prime}(\omega)$ is dominated by the
behaviour of the distribution $G(\log\tau)$ for $\tau\sim 1/\omega$. This
explains why the spectra with and without facilitation coincide at high
frequencies. Instead, at low frequencies, the compression of $G(\log\tau)$ at
large times gives rise to a similar compression at low $\omega$ in the
frequency domain. For example, we have $\log_{10}\tau_{\text{jump}}\sim 7$ for
the largest $\Delta=0.1$. This translates into a sharp peak around
$\log_{10}\omega_{p}\sim-7$ in the corresponding spectrum, see Fig. 3(d).

###  IV.2 Approximate analytic description of the spectra

Next, we want to capture more quantitatively the effect of $\Delta$ on the
distribution of timescales which is ruled by a competition between diffusion
of the trap depth in energy space and thermally-activated jumps. We propose an
approximate analytic treatment of the facilitated trap model, as the complete
dynamics cannot be solved exactly.

We seek an approximate description of the single-trap dynamics before it is
escaped and study the diffusion of one trap in energy space before its escape.
We define the probability distribution $P(E,t;E_{i})$ for the energy $E$ of a
single trap at time $t$, starting from the energy $E_{i}$ at $t=0$. This
distribution follows a Fokker-Planck equation which can be derived from the
master equation with the appropriate Metropolis rule:

$\displaystyle\frac{\partial P}{\partial t}(E,t;E_{i})=$ (7) $\displaystyle
D_{E}\frac{\partial}{\partial E}\left[-\frac{\mathrm{d}\log
P_{\mathrm{eq}}}{\mathrm{d}E}(E)P(E,t;E_{i})+\frac{\partial P}{\partial
E}(E,t;E_{i})\right].$

This describes the drift-diffusion of a single trap in the confining potential
$V_{\rm eff}=-\log P_{\mathrm{eq}}$ (in units of $D_{E}$) and involves a
diffusion constant $D_{E}\sim\Delta^{2}\Gamma_{\mathrm{eq}}$ (up to a
numerical constant), where $\Gamma_{\mathrm{eq}}(T)$ represents the average
jump rate at temperature $T$. The latter only depends on the equilibrium
distribution and reads
$\Gamma_{\mathrm{eq}}=\int_{0}^{+\infty}\mathrm{d}EP_{\mathrm{eq}}(E)e^{-\beta
E}$. To ensure the probability conservation, we impose a hard-wall boundary
condition

$\left.\left[-\frac{\mathrm{d}\log
P_{\mathrm{eq}}}{\mathrm{d}E}(E)P(E,t;E_{i})+\frac{\partial P}{\partial
E}(E,t;E_{i})\right]\right|_{E=0}=0$ (8)

at the origin, while the initial condition is given by
$P(E,t;E_{i})|_{t=0}=\delta(E-E_{i})$.

We then come to the main approximation, which amounts to expressing the
persistence function for this process as

$p(t;E_{i})=e^{-\int_{0}^{t}\mathrm{d}t^{\prime}\Gamma(t^{\prime};E_{i})},$
(9)

with $\Gamma(t;E_{i})$ the average jump rate at time $t$, given by

$\Gamma(t;E_{i})=\int_{0}^{+\infty}P(E,t;E_{i})e^{-\beta E}\mathrm{d}E.$ (10)

We note that the approximation in Eq. (9) yields the exact solution when
$\Delta=0$, as one recovers that $p(t;E_{i})=e^{-te^{-\beta E_{i}}}$ in this
limit. For $\Delta>0$ it is only approximate as we have not described the full
distribution of jumping times of a single trap by combining energy diffusion
with Poisson-distributed relaxation times at fixed energy. Instead we have
only characterised its dynamics by the average jump rate $\Gamma(t;E_{i})$.

Still, our approximation captures the broad features observed in the
simulations of Sec. III, namely, the fact that the dynamics of shallow traps
is unaffected while the one of deep traps is strongly accelerated because they
diffuse towards smaller energies. Indeed, Eq. (9) yields $p(t;E_{i})\sim
e^{-te^{-\beta E_{i}}}$ in the short-time limit, and $p(t;E_{i})\sim
e^{-\Gamma_{\mathrm{eq}}t}$ in the long-time limit, because the distribution
of single-trap energies tends to the equilibrium distribution
$P_{\mathrm{eq}}$ when $t\to+\infty$, by virtue of Eq. (7). The characteristic
timescale for this crossover is given by the typical diffusion timescale over
an energy range of order 1, namely,

$\tau_{\Delta}\sim 1/D_{\mathrm{E}}\sim\Delta^{-2}\Gamma_{\rm eq}^{-1}(T),$
(11)

which depends both on $\Delta$ and $T$. This means that traps with $e^{\beta
E_{i}}\ll\tau_{\Delta}$ are unaffected by the diffusion process and relax as
when $\Delta=0$. For these traps, the typical diffusion timescale is too large
and they relax by thermal activation before their energy can significantly
diffuse. Instead, deep traps with $e^{\beta E_{i}}\gg\tau_{\Delta}$ have ample
time to explore lower energies and thus jump faster, with a typical timescale
controlled by $\tau_{\Delta}$. As $\tau_{\Delta}\sim\Delta^{-2}$, the
crossover energy between these two limits becomes smaller for larger $\Delta$,
in agreement with the scatter plots in Fig. 5.

Figure 6: Relaxation spectra for the facilitated Gaussian trap model
($\alpha=2$ and $\beta=8$) obtained from our approximate solution of the
dynamics, solving Eqs. (3, 6, 9). The broad feature of the direct simulations
in Fig. 3(a) are correctly reproduced (see the direct comparison for
$\Delta=0.01$).

Moving to the average persistence in Eq. (3), it follows that its time decay
is unaffected by facilitation when $t\ll\tau_{\Delta}$ while it decays very
sharply towards $0$ for $t\gg\tau_{\Delta}$ for $\Delta<1$. Physically, the
ratio $\tau_{\mathrm{eq}}/\tau_{\Delta}$ quantifies the effective strength of
the dynamic facilitation, with
$\tau_{\mathrm{eq}}=\int_{0}^{+\infty}P_{\mathrm{eq}}(E)e^{\beta
E}\mathrm{d}E$ the auto-correlation time of the $\Delta=0$ trap model. This
ratio trivially increases with $\Delta$ at constant temperature. More
interestingly, for a fixed value of $\Delta$, this ratio also increases with
decreasing temperature. This shows that $\Delta$ itself does not uniquely
characterise the strength of facilitation in the system.

To assess the quality of our approximate description of the dynamics, we solve
Eq. (7) analytically for the Gaussian case $\alpha=2$ for which a closed
formula for the jump rate $\Gamma(t;E_{i})$ in Eq. (10) can be derived, as
this represents a solvable Ornstein-Uhlenbeck process in energy space Gardiner
_et al._ (1985). Fig. 6 shows the resulting spectra for different values of
$\Delta$ at inverse temperature $\beta=8$. As $\Delta$ increases, we observe
the same trend as the simulation results presented in Fig. 3(a), even though
the quantitative agreement is not exact, due to the approximation involved in
Eq. (9). This suggests that our approximation correctly captures the interplay
between dynamic heterogeneity and dynamic facilitation which accounts for the
shape of the relaxation spectra in the facilitated trap model, and their
evolution with the control parameters.

Figure 7: Evolution with $\Delta$ of the relaxation time
$\tau_{p}=1/\omega_{p}$, with $\omega_{p}$ the location of the peak in the
relaxation spectra $\chi^{\prime\prime}(\omega)$ for various sets of
parameters: ($\alpha,\beta$)=(2,8), (1.5,3.5), (1.1,1.55). Closed symbols are
from direct simulations, open symbols are from the analytic approximation. The
data are correctly described by an effective power law dependence:
$\tau_{p}\sim\Delta^{-\gamma}$, with $\gamma\in\left[1.45,1.85\right]$.

Our approximate solution reveals the existence of a diffusive timescale
$\tau_{\Delta}$ in Eq. (11) which controls the relaxation of the deepest traps
and provides a terminal cutoff to the distribution of relaxation times, which
controls the linear behaviour of $\chi^{\prime\prime}(\omega)$ at low
frequency. We thus analyse in Fig. 7 the dependence of the relaxation time
defined as $\tau_{p}=\omega_{p}^{-1}$, as a function of $\Delta$ in the
simulations of the trap model. The open symbols correspond to the approximate
analytic solution. In the range of $\Delta$ investigated, we find that
$\tau_{p}\sim\Delta^{-\gamma}$ with a fitted exponent
$\gamma\in\left[1.45,1.85\right]$, depending on the parameters ($\alpha$,
$\beta$) of the model. The measured exponent $\gamma$ increases slightly when
moving from the Gaussian model $\alpha=2$ to $\alpha=1.1$. The exponent
remains slightly smaller than prediction $\gamma=2$ emerging from the
diffusive timescale $\tau_{\Delta}$, see Eq. (11). The difference probably
stems from the fact that the spectrum is an averaged quantity receiving
contributions from all traps, including the ones that are less affected by
facilitation. This interpretation is confirmed by our approximate analytic
solution where the clear presence of the diffusive cutoff nevertheless results
in an effective exponent $\gamma\sim 1.7$, different from 2.

##  V Discussion

We constructed and analysed a simple version of a facilitated trap model,
which combines dynamic heterogeneity and kinetic facilitation, recently
identified as the two key ingredients explaining the asymmetric relaxation
spectra in deeply supercooled liquids. Guiselin  _et al._ (2021). We showed
that the model generates asymmetric relaxation spectra, composed of a sharp
peak at low frequencies, and a much flatter and broader high-frequency signal.
This shape qualitatively resembles the dielectric spectra of deeply
supercooled liquids. Because there is no internal microscopic degrees of
freedom inside traps, the obtained spectra lack the microscopic peak exhibited
by experimental data in the THz region.

Figure 8: Two interpretations of the same relaxation spectrum. (a) In the
common picture, an empirical fit to the $\alpha$-peak reveals an “excess wing”
in the high-frequency range (shaded), highlighted by the dashed line. (b) In
our picture, a broad underlying distribution of timescales (blue line) becomes
compressed at low frequencies by dynamic facilitation (arrow). Some of the
$\alpha$-peak then appears in “excess” (shaded) of the intrinsic spectrum.
Symbols correspond to $\alpha=1.1$, $\beta=1.59$ and $\Delta=0.1$.

In Fig. 8(a) we reproduce the common interpretation of relaxation data
generated from the present facilitated trap, where a fit to the main
$\alpha$-peak reveals the existence of an additional “excess” signal at high
frequencies – the excess wing. Based on our results, we propose a different
interpretation of the same spectrum, shown in Fig. 8(b). The self-induced
heterogeneity Berthier (2020) of the supercooled liquid produces a broad
distribution of activation timescales, or, equivalently, a broad spectrum
(blue line). The measured spectrum does not correspond to this underlying
distribution because the relaxation of the fastest regions facilitates the
relaxation of the slow regions. As a result, the low-frequency part of the
underlying spectrum becomes compressed in the presence of facilitation, thus
producing an asymmetric shape, see Fig. 8(b). Somewhat provocatively, we might
say that in this picture, the $\alpha$-peak appears in “excess” of a much
broader underlying time distribution, while the excess wing in fact reveals
the intrinsic shape of the relaxation spectrum.

Figure 9: (a) Frequency dependence of the dielectric loss (relaxation
spectrum) in glycerol at various temperatures, data from Ref. Schneider  _et
al._ (2000). (b) Relaxation spectra of the facilitated trap model for
$\rho(E)\sim e^{-E^{\alpha}}$, $\alpha=1.1$, various temperatures $1/\beta$
and values of $\Delta$. By combining dynamic heterogeneities and dynamic
facilitation, the model qualitatively reproduces the experimentally-observed
behaviour.

The facilitated trap model is an empirical model based on simple, physically-
motivated, ingredients. This is not a microscopic model, it even lacks any
spatial structure. As such, the scope of the model is to describe the essence
of the asymmetric spectra, rather than to provide a quantitative description
of any kind of data.

It is however tempting to push the model to its limits to generate relaxation
spectra which resemble experimental data as closely as possible. We report in
Fig. 9(a) the dielectric loss of supercooled glycerol at various temperatures
down to the glass transition temperature $T_{g}$, reproduced from Ref.
Schneider  _et al._ (2000). We clearly identify a main $\alpha$-peak, which
shifts to lower frequencies with decreasing temperature, as well as an “excess
wing”, which seems to flatten as the temperature decreases. We report in Fig.
9(b) the relaxation spectra of the facilitated trap model for a fixed
$\alpha=1.1$, obtained for various values of $\beta$ and $\Delta$. The model
yields spectra which are in semi-quantitative agreement with experimental
measurements. The temperature variation of the excess wing is in particular
simply reproduced by changing $T$ at constant $\alpha$. The value of $\Delta$,
on the other hand, was adjusted to obtain the correct location of the main
$\alpha$-peak. Although the obtained $\Delta$ values decrease as $T$
decreases, the effect of facilitation, as quantified by the ratio
$\tau_{\mathrm{eq}}/\tau_{\Delta}$ (see Sec. IV.2) steadily increases towards
low temperatures by several orders of magnitude. Therefore, by constraining
our model to yield a semi-quantitative agreement with experimental data, we
conclude that the effect of dynamic facilitation becomes stronger at lower
temperatures, which is fully consistent with our recent atomistic simulations
Guiselin  _et al._ (2021).

There are some differences between the experimental data and results from the
facilitated trap model in Fig. 9. First, the excess wing obtained with the
trap model is not a pure power law, as reported in experiments Menon and Nagel
(1995) (but this is also a debated issue Schneider  _et al._ (2000);
Lunkenheimer  _et al._ (2002)). Strictly speaking, it corresponds to the high-
frequency part of a very broad and very flat relaxation spectrum. However,
since the observed wing extends over at most a few decades, it may be well-
fitted by a power law. Another difference lies in the crossover between the
relaxation peak and the high-frequency spectrum, which is quite marked in the
model, while it appears much smoother in experiments. This sharp crossover may
be partially explained by the lack of spatial resolution of the model, which
may miss some of the heterogeneity present in real liquids.

Our work demonstrates that a simple model combining dynamic heterogeneity and
kinetic facilitation generically produces relaxation spectra that compare well
with experimental and numerical findings. Since the ingredients of the model
are directly motivated by a microscopic analysis of simulation results
Guiselin  _et al._ (2021), we can critically revisit the assumptions made in
alternative models. First, the linear superposition of two processes provides
the misleading picture of distinct and independent molecular motions in the
$\alpha$-peak and in the excess wing. In our model, there are not two distinct
types of dynamics as the asymmetric spectrum emerges from a single type of
relaxation events. Second, our approach suggests that a coupling between
distinct glassy degrees of freedom is not needed to produce excess wings and
complex spectra. Excess wings instead emerge in the facilitated trap model and
in the atomistic simulations from a single type of glassy degree of freedom.
As a consequence, excess wings appear as an intrinsic feature of slow dynamics
near $T_{g}$. Third, our results clearly favour the interpretation of
asymmetric spectra in terms of a single, asymmetric distribution of
timescales. The key role played by kinetic facilitation in our argument is
difficult to reconcile with the thermodynamic nature of the geometric
frustration argument proposed in Ref. Viot  _et al._ (2000). Our model appears
closer in spirit to the one proposed by Chamberlin Chamberlin (1993, 1998,
1999) in which a symmetric underlying distribution (of domain sizes Chamberlin
(1999)) relaxes in a peculiar manner whereby the relaxation time of the large
domains saturates to a finite limit. This is mathematically consistent with
the saturation effect reported in the scatter plots in Fig. 5 for the
facilitated trap model, even though the kinetic argument used by Chamberlin
appears physically unrelated from dynamic facilitation. Finally, our
conclusions regarding the central role played by dynamic facilitation to
account for excess wings are obviously consistent with results obtained with
kinetically constrained models Berthier and Garrahan (2005), and the
qualitative argument proposed in Ref. Xia and Wolynes (2001) to account for
the relative narrowness of the $\alpha$-peak. In all these pictures, the
excess wing at high-frequency appears to result from the full decorrelation of
a small fraction of the system, rather than the small motion of the entire
system advocated in some NMR studies Vogel  _et al._ (2002).

In future work, it would be interesting to develop more realistic coarse-
grained models of structural relaxation in supercooled liquids. An obvious
step would be develop a version of the facilitated trap model with a spatial
structure, to understand better the nature of the spatially heterogeneous
dynamics in the presence of dynamic facilitation, in the spirit of Refs.
Rehwald  _et al._ (2010); Rehwald and Heuer (2012). It would also be
interesting to characterise better the specific nature and geometry of dynamic
facilitation in computer simulations Guiselin  _et al._ (2021); Chacko  _et
al._ (2021), especially at very low temperatures that can now be simulated
more easily thanks to the swap Monte Carlo algorithm.

###### Acknowledgements.

We thank E. Bertin, G. Biroli, U. Buchenau, M. D. Ediger, J. Kurchan, and E.
Rössler for useful discussions and correspondance, and S. R. Nagel for
providing us wth detailed explanations about the discovery of excess wings in
experiments. This work was supported by a grant from the Simons Foundation
(#454933, LB), the European Research Council under the EU’s Horizon 2020
Program, Grant No. 740269, a Herchel Smith Postdoctoral Research Fellowship
(CS), Sidney Sussex College, Cambridge (Ramon Jenkins Research Fellowship to
CS) and Capital Fund Management – Fondation pour la Recherche (BG).

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Ediger  _et al._ (1996) M. D. Ediger, C. A. Angell, and S. R. Nagel, The Journal of Physical Chemistry 100, 13200 (1996).
  * Angell (1995) C. A. Angell, Science 267, 1924 (1995).
  * Böhmer  _et al._ (1993) R. Böhmer, K. L. Ngai, C. A. Angell, and D. J. Plazek, The Journal of Chemical Physics 99, 4201 (1993).
  * Lunkenheimer  _et al._ (2002) P. Lunkenheimer, R. Wehn, T. Riegger, and A. Loidl, Journal of Non-Crystalline Solids 307-310, 336 (2002).
  * Roland  _et al._ (2005) C. M. Roland, S. Hensel-Bielowka, M. Paluch, and R. Casalini, Reports on Progress in Physics 68, 1405 (2005).
  * Ediger (2000) M. D. Ediger, Annual Review of Physical Chemistry 51, 99 (2000).
  * Tarjus (2011) G. Tarjus, Dynamical Heterogeneities in Glasses, Colloids, and Granular Media 150, 39 (2011).
  * Berthier and Biroli (2011) L. Berthier and G. Biroli, Reviews of modern physics 83, 587 (2011).
  * Lunkenheimer  _et al._ (2000) P. Lunkenheimer, U. Schneider, R. Brand, and A. Loid, Contemporary Physics 41, 15 (2000).
  * Flämig  _et al._ (2020) M. Flämig, M. Hofmann, N. Fatkullin, and E. Rössler, The Journal of Physical Chemistry B 124, 1557 (2020).
  * Schmidtke  _et al._ (2013) B. Schmidtke, N. Petzold, R. Kahlau, and E. Rössler, The Journal of chemical physics 139, 084504 (2013).
  * Johari and Goldstein (1970) G. P. Johari and M. Goldstein, The Journal of chemical physics 53, 2372 (1970).
  * Dixon  _et al._ (1990) P. K. Dixon, L. Wu, S. R. Nagel, B. D. Williams, and J. P. Carini, Phys. Rev. Lett. 65, 1108 (1990).
  * Menon and Nagel (1995) N. Menon and S. R. Nagel, Phys. Rev. Lett. 74, 1230 (1995).
  * Menon  _et al._ (1992) N. Menon, K. P. O’Brien, P. K. Dixon, L. Wu, S. R. Nagel, B. D. Williams, and J. P. Carini, Journal of Non-Crystalline Solids 141, 61 (1992).
  * Leheny and Nagel (1997) R. Leheny and S. Nagel, EPL (Europhysics Letters) 39, 447 (1997).
  * Leheny and Nagel (1998) R. L. Leheny and S. R. Nagel, Journal of Non-Crystalline Solids 235-237, 278 (1998).
  * Adichtchev  _et al._ (2003) S. Adichtchev, T. Blochowicz, C. Tschirwitz, V. N. Novikov, and E. A. Rössler, Phys. Rev. E 68, 011504 (2003).
  * Blochowicz  _et al._ (2003) T. Blochowicz, C. Tschirwitz, S. Benkhof, and E. Rössler, The Journal of chemical physics 118, 7544 (2003).
  * Gainaru  _et al._ (2009) C. Gainaru, R. Kahlau, E. A. Rössler, and R. Böhmer, The Journal of chemical physics 131, 184510 (2009).
  * Caporaletti  _et al._ (2021) F. Caporaletti, S. Capaccioli, S. Valenti, M. Mikolasek, A. Chumakov, and G. Monaco, Nature communications 12, 1 (2021).
  * Blochowicz  _et al._ (2006) T. Blochowicz, A. Brodin, and E. A. Rössler, Fractals, Diffusion, and Relaxation in Disordered Complex Systems: Advances in Chemical Physics, Part A 133, 127 (2006).
  * Gainaru (2019) C. Gainaru, Physical Review E 100, 020601 (2019).
  * Ngai and Paluch (2004) K. L. Ngai and M. Paluch, The Journal of Chemical Physics 120, 857 (2004).
  * Wu (1991) L. Wu, Phys. Rev. B 43, 9906 (1991).
  * Vogel  _et al._ (2002) M. Vogel, C. Tschirwitz, G. Schneider, C. Koplin, P. Medick, and E. Rössler, Journal of non-crystalline solids 307, 326 (2002).
  * Petzold  _et al._ (2013) N. Petzold, B. Schmidtke, R. Kahlau, D. Bock, R. Meier, B. Micko, D. Kruk, and E. Rössler, The Journal of chemical physics 138, 12A510 (2013).
  * Körber  _et al._ (2020) T. Körber, R. Stäglich, C. Gainaru, R. Böhmer, and E. A. Rössler, The Journal of Chemical Physics 153, 124510 (2020).
  * Buchenau  _et al._ (2020) U. Buchenau, G. D’Angelo, G. Carini, X. Liu, and M. Ramos, arXiv preprint arXiv:2012.10139 (2020).
  * Kudlik  _et al._ (1998) A. Kudlik, C. Tschirwitz, T. Blochowicz, S. Benkhof, and E. Rössler, Journal of non-crystalline solids 235, 406 (1998).
  * Gainaru  _et al._ (2010) C. Gainaru, R. Böhmer, R. Kahlau, and E. Rössler, Physical Review B 82, 104205 (2010).
  * Diezemann  _et al._ (1999) G. Diezemann, U. Mohanty, and I. Oppenheim, Phys. Rev. E 59, 2067 (1999).
  * Mohanty  _et al._ (2000) U. Mohanty, G. Diezemann, and I. Oppenheim, Journal of Physics: Condensed Matter 12, 6431 (2000).
  * Götze and Sjögren (1989) W. Götze and L. Sjögren, Journal of Physics: Condensed Matter 1, 4183 (1989).
  * Götze and Sperl (2002) W. Götze and M. Sperl, Phys. Rev. E 66, 011405 (2002).
  * Domschke  _et al._ (2011) M. Domschke, M. Marsilius, T. Blochowicz, and T. Voigtmann, Phys. Rev. E 84, 031506 (2011).
  * Chong and Götze (2002) S.-H. Chong and W. Götze, Phys. Rev. E 65, 051201 (2002).
  * Cummins (2005) H. Cummins, Journal of Physics: Condensed Matter 17, 1457 (2005).
  * Sjögren (1986) L. Sjögren, Physical Review A 33, 1254 (1986).
  * Crisanti  _et al._ (2011) A. Crisanti, L. Leuzzi, and M. Paoluzzi, The European Physical Journal E 34, 1 (2011).
  * Crisanti and Leuzzi (2015) A. Crisanti and L. Leuzzi, Journal of Non-Crystalline Solids 407, 110 (2015).
  * Sethna  _et al._ (1991) J. P. Sethna, J. D. Shore, and M. Huang, Phys. Rev. B 44, 4943 (1991).
  * Stevenson and Wolynes (2010) J. D. Stevenson and P. G. Wolynes, Nature physics 6, 62 (2010).
  * Viot  _et al._ (2000) P. Viot, G. Tarjus, and D. Kivelson, The Journal of Chemical Physics 112, 10368 (2000).
  * Chamberlin (1999) R. V. Chamberlin, Phys. Rev. Lett. 82, 2520 (1999).
  * Dyre and Schrøder (2000) J. C. Dyre and T. B. Schrøder, Rev. Mod. Phys. 72, 873 (2000).
  * Berthier and Garrahan (2005) L. Berthier and J. P. Garrahan, The Journal of Physical Chemistry B 109, 3578 (2005).
  * Tarjus  _et al._ (2005) G. Tarjus, S. A. Kivelson, Z. Nussinov, and P. Viot, Journal of Physics: Condensed Matter 17, R1143 (2005).
  * Chamberlin (1993) R. V. Chamberlin, Phys. Rev. B 48, 15638 (1993).
  * Chamberlin (1998) R. V. Chamberlin, Phase Transitions: A Multinational Journal 65, 169 (1998).
  * Garrahan  _et al._ (2011) J. P. Garrahan, P. Sollich, and C. Toninelli, Dynamical heterogeneities in glasses, colloids, and granular media 150, 111 (2011).
  * Chandler and Garrahan (2010) D. Chandler and J. P. Garrahan, Annual Review of Physical Chemistry 61, 191 (2010).
  * Kirkpatrick  _et al._ (1989) T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, Phys. Rev. A 40, 1045 (1989).
  * Lubchenko and Wolynes (2007) V. Lubchenko and P. G. Wolynes, Annual Review of Physical Chemistry 58, 235 (2007).
  * Xia and Wolynes (2001) X. Xia and P. G. Wolynes, Phys. Rev. Lett. 86, 5526 (2001).
  * Bhattacharyya  _et al._ (2008) S. M. Bhattacharyya, B. Bagchi, and P. G. Wolynes, Proceedings of the National Academy of Sciences 105, 16077 (2008).
  * Biroli and Bouchaud (2012) G. Biroli and J.-P. Bouchaud, Structural Glasses and Supercooled Liquids: Theory, Experiment, and Applications , 31 (2012).
  * Berthier  _et al._ (2019a) L. Berthier, G. Biroli, J.-P. Bouchaud, and G. Tarjus, The Journal of Chemical Physics 150, 094501 (2019a).
  * Guiselin  _et al._ (2021) B. Guiselin, C. Scalliet, and L. Berthier, arXiv preprint arXiv:2103.01569 (2021).
  * Berthier  _et al._ (2016) L. Berthier, D. Coslovich, A. Ninarello, and M. Ozawa, Phys. Rev. Lett. 116, 238002 (2016).
  * Ninarello  _et al._ (2017) A. Ninarello, L. Berthier, and D. Coslovich, Phys. Rev. X 7, 021039 (2017).
  * Berthier  _et al._ (2019b) L. Berthier, E. Flenner, C. J. Fullerton, C. Scalliet, and M. Singh, Journal of Statistical Mechanics: Theory and Experiment 2019, 064004 (2019b).
  * Berthier  _et al._ (2011) L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, and W. van Saarloos,  _Dynamical heterogeneities in glasses, colloids, and granular media_ , Vol. 150 (OUP Oxford, 2011).
  * Berthier (2011) L. Berthier, Physics 4, 42 (2011).
  * Brawer (1984) S. Brawer, The Journal of chemical physics 81, 954 (1984).
  * Dyre (1987) J. C. Dyre, Phys. Rev. Lett. 58, 792 (1987).
  * Bässler (1987) H. Bässler, Phys. Rev. Lett. 58, 767 (1987).
  * Bouchaud (1992) J.-P. Bouchaud, Journal de Physique I 2, 1705 (1992).
  * Monthus and Bouchaud (1996) C. Monthus and J.-P. Bouchaud, Journal of Physics A: Mathematical and General 29, 3847 (1996).
  * Jean-Philippe Bouchaud  _et al._ (1995) Jean-Philippe Bouchaud, Alain Comtet, and Cécile Monthus, J. Phys. I France 5, 1521 (1995).
  * Denny  _et al._ (2003) R. A. Denny, D. R. Reichman, and J.-P. Bouchaud, Phys. Rev. Lett. 90, 025503 (2003).
  * Heuer  _et al._ (2005) A. Heuer, B. Doliwa, and A. Saksaengwijit, Phys. Rev. E 72, 021503 (2005).
  * Diezemann and Heuer (2011) G. Diezemann and A. Heuer, Phys. Rev. E 83, 031505 (2011).
  * Diezemann (1997) G. Diezemann, The Journal of Chemical Physics 107, 10112 (1997).
  * Rehwald  _et al._ (2010) C. Rehwald, O. Rubner, and A. Heuer, Phys. Rev. Lett. 105, 117801 (2010).
  * Arkhipov and Baessler (1994) V. I. Arkhipov and H. Baessler, The Journal of Physical Chemistry 98, 662 (1994).
  * Arkhipov  _et al._ (1996) V. Arkhipov, H. Bässler, and D. Khramtchenkov, The Journal of Physical Chemistry 100, 5118 (1996).
  * Bergroth  _et al._ (2005) M. N. Bergroth, M. Vogel, and S. C. Glotzer, The Journal of Physical Chemistry B 109, 6748 (2005).
  * Vogel and Glotzer (2004) M. Vogel and S. C. Glotzer, Phys. Rev. Lett. 92, 255901 (2004).
  * Candelier  _et al._ (2009) R. Candelier, O. Dauchot, and G. Biroli, Phys. Rev. Lett. 102, 088001 (2009).
  * Candelier  _et al._ (2010) R. Candelier, A. Widmer-Cooper, J. K. Kummerfeld, O. Dauchot, G. Biroli, P. Harrowell, and D. R. Reichman, Phys. Rev. Lett. 105, 135702 (2010).
  * Keys  _et al._ (2011) A. S. Keys, L. O. Hedges, J. P. Garrahan, S. C. Glotzer, and D. Chandler, Phys. Rev. X 1, 021013 (2011).
  * Schoenholz  _et al._ (2016) S. S. Schoenholz, E. D. Cubuk, D. M. Sussman, E. Kaxiras, and A. J. Liu, Nature Physics 12, 469 (2016).
  * Chacko  _et al._ (2021) R. N. Chacko, F. P. Landes, G. Biroli, O. Dauchot, A. J. Liu, and D. R. Reichman, arXiv preprint arXiv:2103.01852 (2021).
  * Diezemann (2005) G. Diezemann, The Journal of chemical physics 123, 204510 (2005).
  * Rehwald and Heuer (2012) C. Rehwald and A. Heuer, Phys. Rev. E 86, 051504 (2012).
  * Böhmer  _et al._ (1998) R. Böhmer, R. Chamberlin, G. Diezemann, B. Geil, A. Heuer, G. Hinze, S. Kuebler, R. Richert, B. Schiener, H. Sillescu, _et al._ , Journal of non-crystalline solids 235, 1 (1998).
  * Marsaglia and Tsang (2000) G. Marsaglia and W. W. Tsang, ACM Trans. Math. Softw. 26, 363 (2000).
  * Gardiner  _et al._ (1985) C. W. Gardiner _et al._ ,  _Handbook of stochastic methods_ , Vol. 3 (springer Berlin, 1985).
  * Berthier (2020) L. Berthier, arXiv preprint arXiv:2010.12244 (2020).
  * Schneider  _et al._ (2000) U. Schneider, R. Brand, P. Lunkenheimer, and A. Loidl, Phys. Rev. Lett. 84, 5560 (2000).