††thanks: These two authors contributed equally††thanks: These two authors
contributed equally

# Local barriers control relaxation near the glass transition up to
millisecond time scales

Massimo Pica Ciamarra<EMAIL_ADDRESS>Division of Physics and Applied
Physics, School of Physical and Mathematical Sciences, Nanyang Technological
University, 21 Nanyang Link, 637371, Singapore  Consiglio Nazionale delle
Ricerce, CNR-SPIN, Napoli, I-80126, Italy  Wencheng Ji
<EMAIL_ADDRESS>Department of Physics of Complex Systems, Weizmann
Institute of Science, Rehovot, 234 Hertzl St., Israel  John A. Paulson School
of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138,
USA  Matthieu Wyart<EMAIL_ADDRESS>Institute of Physics, École
Polytechnique Fédérale de Lausanne, Lausanne, CH-1015, Switzerland

###### Abstract

Which phenomenon slows down the dynamics in super-cooled liquids and turns
them into glasses is a long-standing question of condensed-matter. Most
popular theories posit that the dynamics becomes cooperative as the
temperature decreases: many events must occur in a coordinated fashion on a
growing length scale for relaxation to occur. Instead, other approaches
consider that local barriers associated with the elementary rearrangement of a
few particles or ‘excitations’ govern the dynamics. Here, we introduce a new
algorithm, SEER, to resolve this conundrum. SEER can systematically extract
hundreds of excitations and their energy from any given configuration.
Remarkably, we find that the excitation density of states is essentially
shifted to higher energy under cooling. This observation directly predicts how
local barriers slow down the dynamics. We compare this prediction with the
observed dynamics in liquids that can be equilibrated by swap algorithms up to
millisecond time scales. The agreement is quantitative, revealing that
cooperative effects are not controlling the fragility of the liquid, and
suggesting new perspectives on the glass transition.

Glass, Excitation, Landscape

##  I Introduction

Although glass making appeared four thousand years ago [1], and was already
used to make windows in Roman times [2], the physical process preventing
glasses from flowing remains an enduring mystery [3, 4]. As a liquid is
cooled, the relaxation time $\tau$ below which it acts as a solid continuously
grows from picoseconds at high temperatures up to minutes at the glass
transition temperature $T_{g}$. In liquids called strong such as silica [5],
$\tau$ follows an Arrhenius law $\tau=t_{0}\exp(E_{a}/T)$, where $t_{0}$ is
some microscopic time scale, the activation energy $E_{a}$ is constant and $T$
is the temperature, in units where the Boltzmann constant is unity. By
contrast, in liquids called fragile, $E_{a}(T)$ can increase by up to several
folds as $T$ decreases. In both cases, as the dynamics slows down it also
becomes correlated on a growing length scale $\xi(T)$ [6, 7]. These phenomena,
challenging to explain, occur while the static structure of the material
displays limited changes under cooling.

This work focuses on the mechanism inducing the growth of the activation
energy in liquids, which is the central problem of the glass transition. For
our purpose, it is useful to classify theories into two groups. In most
popular approaches, activation is cooperative: the activation energy $E_{a}$
increase stems from the growth of a correlation length $\xi$, as sketched in
Fig. 1b. In some views, including Adam-Gibbs and Random First Order Theory [8,
9, 10], $\xi$ characterizes a growing order associated with a thermodynamic
phase transition. Breaking this order costs an energy $E_{a}\sim\xi^{a}$ where
$a$ is some exponent. Kinetic constraints model (KCMs) [11, 12] are other
theories for which thermodynamics properties are instead trivial. Popular KCMs
consider independent defects that cost an energy $U$, whose density decreases
as $\exp(-U/T)$ under cooling, leading to a growing distance between defects
$\xi$. For certain rules on the kinetic of defects, including the popular
‘East model’ [13], $E_{a}$ grows with $\xi$ and becomes much larger than the
cost $U$ of moving a single defect.

By contrast, free volume [14] or elastic [15, 16, 17] models assume that the
activation energy is not cooperative as sketched in Fig. 1b: it is governed by
the energy barrier of elementary rearrangements, independently of $\xi$. In
elastic models $E_{a}$ is controlled by the high-frequency shear modulus
indeed known to grow in fragile liquids [15] or alternatively by a local
elastic modulus [16, 17]. These views are consistent with the observation that
thermal relaxation occurs more often in softer [18] or more plastic regions
[19, 20]. Obviously, the activation energy may also have a contribution from
localized events and one associated with cooperative dynamics. In some
scenarios, cooperative effects become dominant at low temperatures [21].

Figure 1:  a Temperature dependence of the relaxation time $\tau$ estimated
from the decay of the overlap function in our poly-disperse (see Sec. III). It
is customary to define the onset temperature $T_{\rm on}$ as the temperature
below which the energy of inherent structures starts to decrease. $T_{\rm
on}\simeq 0.7$ [22] for the considered model. We denote by $\tau_{\rm
on}\equiv\tau(T=T_{\rm on})$ the relaxation time at this temperature.
According to mode-coupling theory, the relaxation time diverges as
$\tau\propto(T-T_{c})^{-\gamma}$ (dashed lines). Estimates of $T_{c}$ values
fall in the shaded grey range, determined in Appendix B. b Schematic free-
energy landscape for models whose change of activation energy is dominated by
the change of the energy of local barriers (top), and for models where this
change is dominated by the growth of cooperative effects over a distance
$\xi$, corresponding to a distance in phase space $d_{\xi}$ (bottom). c Energy
change along the minimum energy path associated with a local relaxation event
of our considered model, at $T=0.47$. $\Delta r_{\rm path}$ indicates the norm
of the displacement field along the path in the units of $\rho^{-1/3}$, with
$\rho$ the number density. d Trajectories of the particles with the largest
displacement as the system moves along the minimum energy path with energy
profile illustrated in c. The size of each particle has been scaled down by a
factor of four for visualization purposes.

Even though local and cooperative theories are based on differing physical
assumptions, they both find empirical support through observed correlations
between the quantities they relate. Henceforth, our contention is that to
distinguish these views, one must go beyond correlations and systematically
measure the energy barrier of the localized activation events driving
structural relaxation, illustrated in Fig. 1c and d.

Previous numerical studies related the relaxation dynamics to the exploration
of the potential energy landscape $E(\bf{r})$. At low temperatures, this
exploration proceeds via consecutive activated transitions, or excitations,
through which the system transits from a minimum of $E(\bf{r})$, called
inherent structures (IS), to an adjacent one. Earlier works [23, 24] focused
on rather high temperatures and clarified that irreversibility occurs through
a sequence of elementary excitations, also evidenced experimentally [25, 26,
27], corresponding to a displacement of a temperature-independent number of
particles. Other studies considered the few lowest-energy excitations at lower
temperatures to investigate the plastic and quantum properties of glassy
solids [28, 29, 30, 31, 32]. Such low-energy excitations are well-known to be
already active at temperatures much smaller than the glass transition. They
thus cannot by themselves relax the material into a liquid. Overall, in the
absence of a comprehensive description of how the distribution of energy
barriers of localized relaxation events depends on temperature, it is hard to
establish if the gradual increase in these energy barriers drives the
formation of glasses, as suggested by local theories, or conversely, if the
dynamics’ slowdown results from the increasing spatial organization of these
events, as posited by cooperative theories.

In this work, we introduce SEER, an algorithm to systematically uncover high-
energy excitations in numerical glasses. SEER allows us to measure the
distribution $N(E_{\rm loc})$ of local energy barriers $E_{\rm loc}$ at
various temperatures, around a given configuration. We use it in a system of
continuously polydisperse particles. Such systems are known to reproduce the
main facts of the glass transition, and can be efficiently thermalized [33],
allowing us to explore a wide range of time scales. A remarkable result is
that $N(E_{\rm loc})$ is essentially shifted to higher energy under cooling.
This shift corresponds to an increase in local barriers, which must contribute
to the rise in activation energy $E_{a}(T)$, alongside potential cooperative
effects. Remarkably, our observations show that the increase in activation
energy resulting from these local barriers accurately and fully accounts for
the growth of the relaxation time over seven orders of magnitude. Thus local
effects, rather than cooperative ones, control the dynamics near the glass
transition.

Figure 2:  Schematic illustration of SEER’s working principle. a The algorithm
suppresses a discovered excitation by transforming it into a saddle
bifurcation; b in real space, it corresponds to adding restoring forces
opposite to displacements along the suppressed excitations.

##  II SEER: Systematic Excitation ExtRaction

We seek to uncover the excitations around a given inherent structure, which is
a minimum of $E(\bf{r})$ obtained by quenching an equilibrated configuration,
determining the number of excitations with activation energy $E_{\rm loc}$,
$N(E_{\rm loc})$. A standard method [23, 34, 35] to uncover excitations
consists in performing a small temperature cycle of duration $t_{w}$ at a
temperature $T_{w}$ before descending again the energy, possibly toward a new
minimum. $t_{w}$ and $T_{w}$ are chosen so that typically none or a single
excitation is triggered. In cases where a few excitations are triggered
simultaneously, algorithms were developed to disentangle them if they are
spatially separated [32], which we systematically use below. Unfortunately,
such a method detects a given excitation with a probability $\propto
t_{w}e^{-E_{\rm loc}/T_{w}}$: the first few lowest-energy excitations are
found extremely often, already at low $T_{w}$. At such temperatures, higher
energy ones are essentially never obtained.

To overcome this limitation and access higher-energy excitation, we have
developed a novel algorithm, SEER, whose detailed implementation is specified
in Appendix C. To access high-energy excitations, SEER modifies the energy
functional once an excitation is found. This modification is the minimal one
that ensures that this excitation cannot be triggered again. It allows us to
slowly increase the next cycling temperature $T_{w}$, so as to discover
sequentially higher energy excitations starting from the same configuration.

Specifically, consider an inherent structure ${\bf r^{*}}$. To lighten
notations, we consider ${\bf r^{*}}$ as specifying the reference position for
all particles, such that ${\bf r}$ denotes the displacement field from that
configuration. Once we discover an excitation via a temperature cycle as
discussed above, i.e. a new minimum ${\bf r}_{e}\neq{\bf 0}$ of the potential
energy, we modify the energy landscape along a single direction in phase
space. Specifically, we add a term $\frac{1}{2}k_{e}({\bf r}\cdot{\bf
r}_{e})^{2}$ to the Hamiltonian, which penalizes motion in the direction ${\bf
r}_{e}$ of the discovered excitation. We fix $k_{e}$ to the minimum value
ensuring that $E({\bf r})+\frac{1}{2}k_{e}({\bf r}\cdot{\bf r}_{e})^{2}$ has
not a minimum near ${\bf r}_{e}$, as illustrated in Fig. 2a,b. We repeat this
procedure for each discovered excitation, such that future thermal cycles
search excitations of the SEER energy:

$E_{{\mbox{\tiny{SEER}}}}({\bf r})=E({\bf
r})+\sum_{e=1}^{N_{e}}\frac{1}{2}k_{e}({\bf r}_{e}\cdot{\bf r})^{2},$ (1)

where $N_{e}$ is the number of excitations already in our catalogue. We add to
our catalogue the excitations of $E_{{\mbox{\tiny{SEER}}}}({\bf r})$ that are
also excitations of $E$, and stop the run when $N_{e}$ reaches a prescribed
value. We discuss the stability of SEER with respect to its control parameters
in Appendix E.

Once an excitation is detected, we use the nudge-elastic-band (NEB) method
[36] to determine the minimal energy path connecting the two inherent
structures, as illustrated in Fig. 1c,d. Obviously, this procedure is
performed using the real energy $E({\bf r})$, in the absence of added
potential. It yields both the barrier height $E_{\rm loc}$ and the number of
intermediate maxima.

In principle, SEER’s modification to the Hamiltonian may hinder the discovery
of new excitations in the vicinity of previously discovered ones. SEER
alleviates this effect by choosing a minimal perturbation $k_{e}$ to the
energy. More importantly, in a given run of SEER the order of apparition of
excitations is stochastic. To estimate the number of excitations $N(E_{\rm
loc})$ more accurately, we perform many runs, each delivering a catalogue of
$N_{e}$ excitations. We then merge these into a more complete library, as
detailed in Appendix D. By definition, $N(E_{\rm loc})$ counts the excitations
obtained when the number of merged catalogues diverges. Practically, $N(E_{\rm
loc})$ converges below some $E^{\rm conv}$ that depends on the number of
merged catalogues. Below, we consider 50 catalogues. It will correspond to a
number of excitations $N_{e}$ of the order of a hundred below $E^{\rm conv}$,
implying that SEERS only constraints a tiny fraction of the degrees of
freedom, $N_{e}/dN<1\%$.

##  III Choice of Model

We consider the three-dimensional soft repulsive particles of Ref. [22]. This
system belongs to the class of modern numerical models whose polydispersity
allows us to use the efficient ‘swap’ algorithm [37, 38, 39, 40] to reach
thermal equilibrium on a range of temperatures comparable to experiments. Such
models are receiving a considerable attention [33, 41, 42, 43, 44, 45, 46,
47], as they capture the hallmarks of the glass transition. The size of the
particles $\sigma$ is distributed as $p(\sigma)\propto\sigma^{-3}$ with
$\sigma\in[\sigma_{\rm min},2.2\sigma_{\rm min}]$. Similar systems capture
both the usual dynamical [39] and thermodynamical [48] properties of
experimental liquids. The interaction is given by the purely repulsive
potential

$V(r_{ij})=\epsilon\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{10}+\sum_{l=0}^{3}c_{2l}\left(\frac{r_{ij}}{\sigma_{ij}}\right)^{2l}\right]$
(2)

for $r_{ij}<x_{c}=1.4$ and vanishes at larger distances, with $c_{2l}$ chosen
to enforce continuity at $x_{c}$ up to $3$ derivatives, and with an
interaction range $\sigma_{ij}$ function of $\sigma_{i}$ and $\sigma_{j}$. We
study systems of $N=2000$ particles of mass $m$ at number density $\rho=0.58$
in cubic simulation boxes with periodic boundary conditions. We express mass
in units of $m$, temperature in units of $\epsilon$, lengths in units
$\rho^{-1/3}$, and time in units of $\sqrt{m\sigma_{\rm min}^{2}/\epsilon}$.

We then investigate the relaxation dynamics of equilibrium configurations via
molecular dynamics simulations (MD), as detailed in Appendix. We extract the
relaxation time $\tau$ from the decay of the overlap function

$Q(t)=\frac{1}{N}\Biggl{\langle}\sum_{i}^{N}e^{\left(\frac{\Delta
r_{i}(t)}{w}\right)^{2}}\Biggr{\rangle}$ (3)

where $\Delta r_{i}(t)$ is the displacement of particle $i$ during a time
interval $t$, and $w=0.3\rho^{-1/3}$. We define the relaxation time as
$Q(\tau)=1/e$.

Using MD from equilibrated configurations, we can access the dynamics on about
6 decades. However, previous studies [43, 44] have found in related models
that the shape of the correlation function does not significantly depend on
the temperature [49], a property referred to as time-temperature superposition
(TTS) [50]. We have verified in Appendix F that this property holds in our
model and use it to infer $\tau$ even when the correlation function does not
fully decay in our observation window, leading to the nine decades of
observation of Fig.1a. As has been reported in other highly-polydisperse [51],
in more conventional model systems [52] or in experimental data [53], the
model seems to display Arrhenius-like behavior (corresponding to straight
lines in the logarithmic representation of Fig.1a) above the onset temperature
$T_{\rm on}\simeq 0.7$ as well as at the lowest temperatures, with a super-
Arrhenius behavior in between. Our analysis below (using reheating methods)
will show that the activation energy is, in fact, not constant in any of the
temperatures probed.

##  IV Results of SEER

Figure 3:  a The number $N(E_{\rm loc})$ of distinct excitations with
activation energy $E_{\rm loc}$ saturates with the number of considered
catalogues, starting from small energies (Appendix E). $N(E_{\rm loc})$
increases by less than $5\%$ below an energy threshold (dashed line) when
adding $10$ more catalogues. The figure refers to $T=0.4$. b The discovered
excitations (circles) are simple (truly elementary) excitations with a single
energy maximum (squares) along the minimum energy path or, instead, present
multiple maxima. c Below the energy convergence threshold, a large fraction of
the discovered excitations are elementary. d The cumulative number of
elementary excitations shifts to higher energies as the temperature decreases.
At each $T$, we have at least $70$ elementary excitations with activation
energy below the convergence threshold. e The collapse of the cumulative
number of excitations on a master curve allows us to associate an energy shift
$\Delta E(T)$ to each temperature. (inset) Temperature dependence of the
energy shift $\Delta E(T)$, using $T_{c}=0.5$ as a reference. f Temperature
dependence of the number of particles $Np$ involved in an excitation, averaged
over all excitations with energy below the convergence threshold and on the
elementary one with a single maximum. Excitations become increasingly
localized as the temperature decreases. In d-f, data are averaged over $10$
configurations with $N=2000$. In e,f $N_{e}=50$.

###  IV.1 Temperature shifts the excitation density of state $N(E_{\rm loc})$

Fig. 3a illustrates $N(E_{\rm loc})$ for $T=0.4$ extracted from SEER, as a
function of the number of catalogs used to build the library. After 50 runs
$N(E_{\rm loc})$ approximately converges below some energy $E^{\rm conv}$. We
uncover about 200 excitations with $E_{\rm loc}<E^{\rm conv}$ (dashed line).
Many of these excitations have a single energy maximum along their reaction
path, as illustrated in Fig. 3b. The density $N(E_{\rm loc})$ of such single
excitations is compared with the total density in Fig. 3c. Considering all or
only single excitations lead to similar results, in what follows we focus on
the latter.

To study the temperature dependence of excitations, we consider the cumulative
number of elementary excitations $F_{T}(E_{\rm loc})$ in Fig. 3d. The figure
reveals a depletion of low-energy excitations under cooling. Our key
observation, illustrated in Fig. 3e, is that different curves $F_{T}(E_{\rm
loc})$ approximately collapse by simply shifting the energy by some value
$\Delta E(T)$. We show the temperature dependence for the energy shift $\Delta
E(T)$ in the inset of Fig. 3e. The energy shift approaches a constant at
temperatures close to the crossover one. It is consistent with the usual
notion that above the crossover temperature $T_{\rm on}\simeq 0.7$, the
properties of the inherent structures become temperature independent. As the
temperature decreases, $\Delta E(T)$ develops an approximate linear
temperature dependence.

Thus, the density of excitations evolves under cooling as if all excitations
were increasing their energy by some amount $\Delta E(T)$. Making the
assumption that the activation energy $E_{a}(T)$ is controlled by local
barriers then implies that it must also shift by $\Delta E(T)$ under cooling.
In other words, for $T_{1},T_{2}$ at low temperatures we must have:

$E_{a}(T_{1})-E_{a}(T_{2})=\Delta E(T_{1})-\Delta E(T_{2})$ (4)

This is our main prediction, which we can test below without fitting
parameters.

###  IV.2 Excitations become more localized under cooling

To show that individual excitations themselves do not become more and more
cooperative under cooling, we estimate the number of particles they involve
from their participation ratio. This number is $Np$, where $p=\frac{(\sum{\bf
u}_{i}^{2})^{2}}{N\sum{\bf u}_{i}^{4}}$ and $\bf u$ the exitation displacement
field. Fig. 3f illustrates the temperature dependence of $Np$ averaged over
all excitations below the energy convergence threshold. Excitations involve a
few particles and become increasingly localized at low temperatures.

The small number of particles involved in each excitation suggests that the
properties of the latter do not depend on the system size. Indeed, previous
studies [24] demonstrated that the statistical features of the energy
landscape of very small systems ($N<100$) already represent the thermodynamic
limit.

##  V Local barriers control relaxation

Figure 4: a Relaxation time as a function of the inverse temperature, as in
Fig.1a (circles). Solid squares indicate the values $\tau(T\to T^{\prime})$ of
the relaxation time of systems equilibrated at temperature $T$ when the
temperature is instantaneously changed to the value $T^{\prime}$. At each $T$,
lines are fit to $\tau(T\to T^{\prime})=t_{0}e^{E_{a}(T)/T^{\prime}}$ on both
reheating and equilibrated data so as to extract $t_{0}$ and $E_{a}(T)$.
Squares and circles overlap when $T=T^{\prime}$. The value of $t_{0}$ has no
clear $T$ dependence and, at each $T$, it slightly depends on the fitting
range. We adopt a conservative approach to account for this variability by
fixing $t_{0}$ to an average value and using three standard deviations as an
estimation of its error: $t_{0}\simeq(2.7\pm 1.3)10^{-5}\tau_{\rm on}$, where
$\tau_{\rm on}$ is the relaxation time at the onset temperature. At $T=0.5$,
we find $E_{a}(T=0.5)\simeq 9.4$ (other values of $E_{a}(T)$ can be extracted
from panel b). Interestingly, the fragility of this liquid is apparent, as the
solid lines (corresponding to constant activation energy) increase
systematically less fast than the $\tau(T)$ under cooling. Blue squares
indicate our predictions $\tau_{\rm pred}(T)$ over almost five orders of
magnitudes which follows from Eq.4 and from our measurement of the shift in
excitations’ activation energy: $\tau_{\rm pred}(T)=t_{0}e^{(E_{a}(0.5)+\Delta
E(T)-\Delta E(0.5))/T}$. b Change in activation energy $E_{a}(T)-E_{a}(0.5)$
(circles) inferred from panel a. Errors are dominated by the uncertainty in
$t_{0}$. The change in activation energy is quantitatively predicted, without
fitting parameters, by the excitations’ activation energy shift, $\Delta
E(T)-\Delta E(0.5)$ (squares). This result indicates that the change of
activation energy is governed by local barriers, instead of cooperative
effects on a growing length $\xi$. Finally, if the activation energy is
extracted from panel a using only the reheated data, one obtains a consistent
estimation of $E_{a}(T)-E_{a}(0.5)$ (diamonds).

The prediction Eq. 4 can be written alternatively as $E_{a}(T)=E_{0}+\Delta
E(T)$, which implies for the relaxation time grows:

$\tau(T)=t_{0}e^{\frac{E_{0}+\Delta E(T)}{T}}$ (5)

In this expression, two parameters enter, $t_{0}$ and $E_{0}$. The latter can
be readily removed to test our hypothesis by considering the change of
activation energy with temperature as in Eq. 4. Yet, this quantity is
interesting to estimate. The microscopic time $t_{0}$ depends on a diversity
of factors, including the entropy of individual barriers [54], and no formal
expression for it is available for the under-damped dynamics we use. It could
depend on temperature. Yet if the change of dynamics is dominated by the
energy change of local barriers, the evolution of $t_{0}$ can be neglected.
Below we test such a model, where $t_{0}$ is fixed.

There is no established procedure to measure the activation energy $E_{a}$ and
$t_{0}$. Previous works, see e.g. [55], focused on the apparent activation
energy $E_{\rm app}=\frac{dlog\tau}{dT}=E_{a}(T)-T\frac{dE_{a}}{dT}$.
Unfortunately, it overestimates the actual activation energy unless the latter
is constant [56]. Regarding this point, note that while an Arrhenius-like
behavior is often observed at low temperatures, it does not indicate a
constant activation energy as often assumed, but rather an activation energy
that depends linearly on temperature. The apparent activation energy cannot
distinguish these two scenarios.

We overcome these difficulties by introducing a novel empirical approach to
determine $t_{0}$ and $E_{a}(T)$. We consider equilibrated configurations at
some temperature $T$ and suddenly change the temperature to some value
$T^{\prime}$. If $T^{\prime}$ and $T$ do not greatly differ 111If $T^{\prime}$
is significantly larger than $T$, relaxation can occur via the nucleation of a
hot liquid phase invading the material [49, 79], we expect the instantaneous
relaxation time $\tau(T\to T^{\prime})$ to follow an Arrhenius behavior in
$T^{\prime}$ with a constant activation energy $E_{a}(T)$: $\tau(T\to
T^{\prime})=t_{0}e^{E_{a}(T)/T^{\prime}}$. Fig. 4a tests this expression
successfully as both $T$ and $T^{\prime}$ are varied. $t_{0}$ has no clear
temperature dependence, and all curves are well-fitted by a single $t_{0}$
value. From these measurements, we extract $E_{a}(T)$.

From these results, the hypothesis that the change of the activation energy
under cooling is governed by the change of local barriers energy, formalized
by Eq. 4, can now be tested. Fig. 4b compares these two changes, respectively
$E_{a}(T)-E_{a}(T=0.5)$ and $\Delta E(T)-\Delta E(T=0.5)$, taking $T=0.5$ as a
reference temperature. Remarkably, the agreement is quantitative. Equivalently
stated, this result implies that knowing the relaxation time at some
intermediate temperature, e.g. $T=0.5$, we can predict the dynamics at lower
temperatures and thus the fragility of the system, as illustrated in Fig. 4a.

These results do not imply that jumping over a single barrier is generally
sufficient to relax the system into a new state of similar energy. Excitations
tend to be in their lower-energy state, thus performing a single jump most
often increases energy, and can be followed by the reverse jump. Denote by
$n_{0}$ the number of local barriers that must act sequentially in the same
location to form an irreversible event, i.e. to find a new state of energy
similar to the original one. $n_{0}$ characterizes the cooperativity of the
dynamics. In theories such as RFOT or the East model, the growth of the
activation energy $E_{a}$ under cooling is due to a growth of cooperativity
$n_{0}$. Such an effect would contribute to the activation energy, in addition
of the individual growth of local barriers we computed. This contribution is
undetectable, thus the dependence $n_{0}$ on temperature can be neglected.

##  VI Conclusion

We have introduced SEER, a novel method to analyze the energy landscape of
glasses, which systematically measures the elementary rearrangements around
any given energy minimum. Our central result is that the density of state of
these excitations is shifted toward higher energy under cooling. Using a novel
method based on rapid temperature changes to analyze the dynamics and extract
the activation energy, we showed that this shift precisely predicts the liquid
slowing down. This finding reveals that up to milliseconds time scales for
which such liquids can be simulated, fragility is controlled by the change of
energy of individual excitations, which shrink in size under cooling. It is
not due to the fact that the dynamics becomes more and more cooperative, as is
often proposed.

At the practical level, the reheating methodology we introduced as well as
SEER have potential to be of broad applicability. Concerning the former,
developing sufficiently fast experimental temperature jumps [58] would be of
high interest to measure directly the activation energy in super-cooled
liquids. Concerning SEER, it is more directly usable in numerical models where
the energy landscape of inherent structures is indicative of the dynamics. It
excludes hard sphere interactions or interactions with very narrow wells [59,
60, 61], where the free energy landscape must be analyzed instead.
Implementing SEER would then require more involved methods, such as using
effective potentials [62, 63, 64, 65, 66].

We focused on modern numerical liquids made of poly-disperse particles, which
are receiving an immense attention [33, 41, 42, 43, 44, 45, 46, 47] as they
can be simulated on time scales closer to experimental conditions. It could be
argued that these systems are specific, with small particles moving faster
than larger ones [67]. Yet, just as for multicomponent metallic glasses [68,
69, 70] or silica where some atoms move much faster than others, these models
present the usual hallmarks of glassy dynamics [51, 52]. In specific systems,
however (e.g. such as those displaying a liquid-liquid phase transition [71]
or hexatic ordering [10]), cooperative effects may influence the excitation
energy, which will be interesting to investigate further.

Conceptually, our results indicate that a dynamical correlation length $\xi$
can grow while not affecting the relationship between relaxation time and
local barriers. Such a scenario is captured in simplified descriptions of
interacting excitations called ‘elasto-plastic models’, which are cellular
automata initially introduced to study the plasticity of amorphous materials
[72]. Their relaxation time is controlled, as in our case, by the density of
excitations around visited configurations $N(E_{\rm loc})$ [73, 74], which
also presents a gap or depletion similar to our observations. These models
display a growing correlation length under cooling [74]. It can be explained
in terms of slow avalanches of activated processes [75], whose statistics
matches the longest available MD simulations of dynamical heterogeneities
[76]. Thus, a description of the glass transition based on local barriers
appears quantitatively consistent with both the slowing down of the dynamics
presented here, and with the spatial correlations of the dynamics studied in
these works.

Finally, an important question for the future is what makes local barriers
grow under cooling. Do they follow a high-frequency elastic modulus [15], or
instead local (instead of global) elasticity [17], perhaps affected by the
presence of locally-favored structures [77, 78]? Alternatively, are they
controlled by the varying geometry of elementary rearrangements under cooling,
which appears to be shaped by the presence of a dynamical transition near
$T_{c}$ [34]? The systematic extraction of excitations SEER allows may give a
new handle to address this question in the future.

###### Acknowledgements.

We thank E. Lerner for providing some of the equilibrium configurations we
investigated. We thank the Simons collaboration as well as L. Berthier, G.
Biroli, C. Brito, C. Gavazzoni, E. Lerner, M. Muller, M. Popovic, M. Ozawa and
A. Tahaei for discussions. M.P.C. discloses support for the research of this
work from Singapore Ministry of Education [MOE-T2EP50221-0016]. M.W
acknowledges support from the Simons Foundation Grant (No. 454953 Matthieu
Wyart) and from the SNSF under Grant No. 200021-165509.

##  Appendix A Extremely low-temperatures

In Fig. 3e, we evaluated the energy shift at small temperatures $T=0.35$ and
$0.30$, whose large relaxation time is inaccessible to numerical simulations.
Figure 5b shows that the predicted shift in activation energy keeps its linear
temperature dependence also in this regime, as also reported in some
experimental data [53].

Figure 5: a Predicted relaxation time and b change in activation/gap energy up
to low-temperature values whose relaxation time is not accessible to numerical
simulations.

##  Appendix B Mode Coupling Critical Temperature

According to mode-coupling theory, the functional from $(T-T_{c})^{-\gamma}$
should describe the relaxation time for $\tau>\tau_{\rm on}$. We estimate
$T_{c}$ as a useful reference by examining the behavior of the relaxation time
for selected putative values of $T_{c}$ in Fig. 6. For $T_{c}>0.54$, MCT only
describes a small increase in relaxation time, up to $\tau\leq 10\tau_{\rm
on}$. For $T_{c}<0.51$, MCT underestimates the relaxation time in an
increasingly large temperature window. We thus estimate $T_{c}\in[0.51:0.54]$.

Figure 6: Power-law divergence of the relaxation time for selected putative
values of the critical temperature $T_{c}$. Data are shifted vertically for
clarity. Dashed lines mark the value of the relaxation time at the onset
temperature, $\tau_{\rm on}$.

##  Appendix C SEER: Systematic Excitation ExtRaction

SEER constructs a catalogue of excitations associated with an inherent
structure ${\bf r}=0$ of a system governed by the energy functional $E$. To
identify the excitations with the smallest activation energy, SEER searches
excitations via thermal cycles performed at a temperature $T_{w}$, which is
initially small $T_{w}\simeq 0.05$, and gradually increased by $\Delta
T_{w}\simeq 0.005$. The velocities ${\bf v(t=0)}$ of the initial configuration
are initially drawn from the equilibrium distribution at $T_{w}$. This is the
only source of randomness in the algorithm. When $T_{w}$ is changed, particle
velocities are set by rescaling these initial velocities.

We detail SEER’s working principle below:

  1. 1.

We set the temperature to $T_{w}$ by rescaling the initial velocities ${\bf
v}$.

  2. 2.

We perform a thermal cycle evolving the system up to a maximum time $t_{w}$,
while controlling the temperature via velocity rescaling. Thermal cycles at
temperature $T_{w}$ and duration $t_{w}$ uncover excitations with a typical
maximal activation energy $\propto T_{w}\log(t_{w}/t_{0})$. As such, if the
thermal cycle reaches the maximum allocated time without uncovering any
excitation, we expand the activation energy observation window increasing
$T_{w}$ by $\Delta T_{w}$ by rescaling the initial velocities set in 1.
Afterwards, we repeat this step. On the contrary, if the thermal cycle detects
an excitation, we perform the following steps.

  3. 2.

We determine by the dichotomy method the smallest temperature in the interval
$[T_{w}-\Delta T_{w}:T_{w}]$ at which a thermal cycle detects an excitation.
Temperature is always controlled by rescaling the velocities set in 1. By
lowering $T_{w}$, we increase the probability we resolve an elementary
transition.

  4. 3.

If the excitation is complex, we decompose it into elementary excitations
${\bf r}_{e}$ as in Ref. [32].

  5. 4.

If the excitations of $E_{{\mbox{\tiny{SEER}}}}$ are also excitations of $E$:

     * •

study the excitation with the nudge-elastic-band and add it to the catalogue.

     * •

for each excitation ${\bf r}_{e}$, we find the minimum $k_{e}$ that makes
$E+\frac{1}{2}k_{e}({\bf r}\cdot{\bf r}_{e})$ loose this excitation:
descending this modified energy from ${\bf r}_{e}$ goes back to ${\bf r}=0$.
We determine $k_{e}$ with $5\%$ accuracy through the dichotomy method.

     * •

We modify $E_{{\mbox{\tiny{SEER}}}}$ by adding the excitation suppressing
terms $\frac{1}{2}k_{e}({\bf r}\cdot{\bf r}_{e})$.

  6. 5.

We repeat the algorithm from 1 to search for a novel excitation.

##  Appendix D Removal of duplicate excitations

We associate many SEER’s generated catalogues with the same initial
configuration by changing the initial particle velocities. Merging these
catalogues yields our final catalogue. When merging catalogues, we remove
repeated excitations. The distribution of the distance $d_{ij}=\lVert{\bf
r}_{e}^{i}-{\bf r}_{e}^{j}\rVert$ between discovered excitations illustrated
in Fig.7 reveals that two excitations are coincident, $d_{ij}\simeq 0$, or
well separated. This finding allows us to robustly ensure our final catalogue
only contains distinct excitations by using a threshold on the distance.

Figure 7: Probability distribution of the distance $d_{ij}=\lVert{\bf
r}_{e}^{i}-{\bf r}_{e}^{j}\rVert$ between discovered excitations. In our final
catalogue, all excitations are at a distance larger than the indicated
threshold.

##  Appendix E Stability with respect to $N_{E}$ and $t_{w}$

Figure 8: a Activation energy dependence of the number of discovered
excitations at $T=0.50$ using 50 catalogues, each containing $N_{E}$
excitations. Convergence shifts towards higher activation energies as $N_{E}$
increases. b In the energy convergence range, the number of discovered
excitations $N(E)$ does not sensibly depend on the duration $t_{w}$ of the
thermal cycles.

We estimated in Fig. 4e and f excitations’ energy shift by analysing the
distribution of the energy barriers associated with the inherent structures of
equilibrium configurations, at each temperature. The data reported in the
figure are (i) obtained by quenching equilibrium configurations to a small
temperature $T=0.3$, and then quenching to zero temperature at a rate
$r=10^{-3}$; (ii) averaged over ten independent samples; and (iii) restricted
to elementary excitations with a single maximum. Fig. 9 illustrates the
robustness of our result by comparing different approaches to estimate the
energy shift $\Delta E$ of each configuration.

Concerning (i), we avoid instantaneous quenches as it leads to a few outlier
configurations with atypical $\Delta E$ at high $T$, as apparent in Fig. 9.
Direct inspection indicated that these configurations present at least one
excitation in its high-energy state, very close to a saddle transition (i.e.
with a small barrier). This situation can lead to instabilities in SEER, if
these excitations relax to low their low-energy state while other excitations
are studied. Eliminating such configurations (not shown), or annealing the
system, eliminates this problem and leads to consistent results independently
of the annealing rate, as shown in Fig. 9.

Concerning (ii), from the typical variance $\sigma^{2}$ of the ten data points
at each temperature, we find our typical error in the estimation of $\Delta E$
to be $\Delta\Delta E\simeq\sigma/\sqrt{10}\simeq 0.1$.

Concerning (iii), the analysis of all excitations (including complex ones,
with multiple maxima as shown in open triangles), rather than the elementary
ones, leads to analogous estimations of the gap energy at low temperatures,
and slightly smaller estimations at the highest temperatures. This observation
is plausibly related to the fact that, rigorously speaking, the density of
states of all excitations is ill-defined. Indeed, in a large system, pairs of
elementary excitations could be combined to generate a very large number of
complex excitations. SEER seeks to avoid these fortuitous combinations (that
can always occur by chance) by parsing complex excitations into elementary
ones. Parsing is simple when two excitations are far away but difficult if
they overlap. Overlaps will be more frequent at high temperatures, where
excitations are more extended.

Figure 9: Comparison between alternative estimations of the gap energy.

##  Appendix F Time temperature superposition

Figure 10: The overlap correlation functions collapse on a stretched
exponential function with stretching exponent $\beta\simeq 0.52$ (full line)
when time is rescaled by the relaxation time. We use this property to estimate
the relaxation time at the lowest temperatures, at which the overlap function
does not fully decay within our simulation time.

## References

  * Shortland  _et al._ [2007] A. Shortland, N. Rogers, and K. Eremin, Trace element discriminants between egyptian and mesopotamian late bronze age glasses, Journal of Archaeological Science 34, 781 (2007).
  * Fleming [1999] S. J. Fleming,  _Roman glass: reflections on cultural change_ (UPenn Museum of Archaeology, 1999).
  * Anderson [1995] D. Anderson, Through the glass lightly, Science 267, 1618 (1995).
  * Debenedetti and Stillinger [2001] P. Debenedetti and F. Stillinger, Supercooled liquids and the glass transition, Nature 410, 259 (2001).
  * Angell [1985] C. Angell, Strong and fragile liquids, Relaxations in complex systems , 3 (1985).
  * Kob  _et al._ [1997] W. Kob, C. Donati, S. J. Plimpton, P. H. Poole, and S. C. Glotzer, Dynamical heterogeneities in a supercooled lennard-jones liquid, Physical review letters 79, 2827 (1997).
  * Dalle-Ferrier  _et al._ [2007] C. Dalle-Ferrier, C. Thibierge, C. Alba-Simionesco, L. Berthier, G. Biroli, J.-P. Bouchaud, F. Ladieu, D. L’Hôte, and G. Tarjus, Spatial correlations in the dynamics of glassforming liquids: Experimental determination of their temperature dependence, Phys. Rev. E 76, 041510 (2007).
  * Lubchenko and Wolynes [2001] V. Lubchenko and P. Wolynes, Intrinsic Quantum Excitations of Low Temperature Glasses, Phys. Rev. Lett. 87, 195901 (2001).
  * Biroli and Bouchaud [2012] G. Biroli and J.-P. Bouchaud, The random first-order transition theory of glasses: a critical assessment, Structural Glasses and Supercooled Liquids: Theory, Experiment, and Applications , 31 (2012).
  * Tanaka  _et al._ [2010] H. Tanaka, T. Kawasaki, H. Shintani, and K. Watanabe, Critical-like behaviour of glass-forming liquids, Nat. Mater. 9, 324 (2010).
  * Keys  _et al._ [2011] A. S. Keys, L. O. Hedges, J. P. Garrahan, S. C. Glotzer, and D. Chandler, Excitations are localized and relaxation is hierarchical in glass-forming liquids, Physical Review X 1, 021013 (2011).
  * Ritort and Sollich [2003] F. Ritort and P. Sollich, Glassy dynamics of kinetically constrained models, Advances in physics 52, 219 (2003).
  * Garrahan and Chandler [2002] J. P. Garrahan and D. Chandler, Geometrical explanation and scaling of dynamical heterogeneities in glass forming systems, Physical review letters 89, 035704 (2002).
  * Turnbull and Cohen [1961] D. Turnbull and M. H. Cohen, Free-volume model of the amorphous phase: glass transition, The Journal of chemical physics 34, 120 (1961).
  * Dyre [2006] J. Dyre, Colloquium: The glass transition and elastic models of glass-forming liquids, Rev. Mod. Phys. 78, 953 (2006).
  * Rainone  _et al._ [2020] C. Rainone, E. Bouchbinder, and E. Lerner, Pinching a glass reveals key properties of its soft spots, Proceedings of the National Academy of Sciences 117, 5228 (2020).
  * Kapteijns  _et al._ [2021] G. Kapteijns, D. Richard, E. Bouchbinder, T. B. Schrøder, J. C. Dyre, and E. Lerner, Does mesoscopic elasticity control viscous slowing down in glassforming liquids?, The Journal of Chemical Physics 155, 74502 (2021).
  * Widmer-Cooper  _et al._ [2008] A. Widmer-Cooper, H. Perry, P. Harrowell, and D. Reichman, Irreversible reorganization in a supercooled liquid originates from localized soft modes, Nat. Phys. 4, 711 (2008).
  * Li  _et al._ [2022] Y.-W. Li, Y. Yao, and M. P. Ciamarra, Local Plastic Response and Slow Heterogeneous Dynamics of Supercooled Liquids, Physical Review Letters 128, 258001 (2022).
  * Lerbinger  _et al._ [2022] M. Lerbinger, A. Barbot, D. Vandembroucq, and S. Patinet, Relevance of Shear Transformations in the Relaxation of Supercooled Liquids, Physical Review Letters 129, 195501 (2022), arXiv:2109.12639 .
  * Mirigian and Schweizer [2013] S. Mirigian and K. S. Schweizer, Unified theory of activated relaxation in liquids over 14 decades in time, Journal of Physical Chemistry Letters 4, 3648 (2013).
  * Lerner [2019] E. Lerner, Mechanical properties of simple computer glasses, J. Non-Cryst. Solids 522, 119570 (2019).
  * Doliwa and Heuer [2003] B. Doliwa and A. Heuer, What does the potential energy landscape tell us about the dynamics of supercooled liquids and glasses?, Physical Review Letters 91, 235501 (2003).
  * Heuer [2008] A. Heuer, Exploring the potential energy landscape of glass-forming systems: from inherent structures via metabasins to macroscopic transport, J. Phys. Condens. Matter 20, 373101 (2008).
  * Simmons  _et al._ [2012] D. S. Simmons, M. T. Cicerone, Q. Zhong, M. Tyagi, and J. F. Douglas, Generalized localization model of relaxation in glass-forming liquids, Soft Matter 8, 11455 (2012).
  * Cicerone  _et al._ [2014] M. T. Cicerone, Q. Zhong, and M. Tyagi, Picosecond Dynamic Heterogeneity, Hopping, and Johari-Goldstein Relaxation in Glass-Forming Liquids, Physical Review Letters 113, 117801 (2014).
  * Cicerone  _et al._ [2023] M. T. Cicerone, K. Badilla-Nunez, J. Zahn, J. P. Stoppelman, and J. G. McDaniel, Excitation populations provide a thermodynamic order parameter for liquids (2023), arXiv:2201.12593 [cond-mat.soft] .
  * Ji  _et al._ [2019] W. Ji, M. Popović, T. de Geus, E. Lerner, and M. Wyart, Theory for the density of interacting quasilocalized modes in amorphous solids, Phys. Rev. E 99, 023003 (2019).
  * Schober  _et al._ [1993] H. Schober, C. Oligschleger, and B. Laird, Low-frequency vibrations and relaxations in glasses, J. Non. Cryst. Solids 156-158, 965 (1993).
  * Wang  _et al._ [2019] L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel, and E. Flenner, Low-frequency vibrational modes of stable glasses, Nat. Commun. 10, 26 (2019).
  * Khomenko  _et al._ [2020] D. Khomenko, C. Scalliet, L. Berthier, D. Reichman, and F. Zamponi, Depletion of two-level systems in ultrastable computer-generated glasses, Phys. Rev. Lett. 124, 225901 (2020).
  * Ji  _et al._ [2020] W. Ji, T. de Geus, M. Popović, E. Agoritsas, and M. Wyart, Thermal origin of quasilocalized excitations in glasses, Phys. Rev. E 102, 062110 (2020).
  * Ninarello  _et al._ [2017a] A. Ninarello, L. Berthier, and D. Coslovich, Models and algorithms for the next generation of glass transition studies, Phys. Rev. X 7, 021039 (2017a).
  * Ji  _et al._ [2022] W. Ji, T. W. de Geus, E. Agoritsas, and M. Wyart, Mean-field description for the architecture of low-energy excitations in glasses, Physical Review E 105, 044601 (2022).
  * Khomenko  _et al._ [2021] D. Khomenko, D. Reichman, and F. Zamponi, Relationship between two-level systems and quasilocalized normal modes in glasses, Phys. Rev. Materials 5, 055602 (2021).
  * Zarkevich and Johnson [2015] N. A. Zarkevich and D. D. Johnson, Nudged-elastic band method with two climbing images: Finding transition states in complex energy landscapes, The Journal of Chemical Physics 142, 24106 (2015).
  * Briano and Glandt [1984] J. Briano and E. Glandt, Statistical thermodynamics of polydisperse fluids, J. Chem. Phys. 80, 3336 (1984).
  * Gutiérrez  _et al._ [2015] R. Gutiérrez, S. Karmakar, Y. Pollack, and I. Procaccia, The static lengthscale characterizing the glass transition at lower temperatures, Europhys. Lett. 111, 56009 (2015).
  * Ninarello  _et al._ [2017b] A. Ninarello, L. Berthier, and D. Coslovich, Models and Algorithms for the Next Generation of Glass Transition Studies, Phys. Rev. X 7, 021039 (2017b).
  * Brito  _et al._ [2018] C. Brito, E. Lerner, and M. Wyart, Theory for Swap Acceleration near the Glass and Jamming Transitions for Continuously Polydisperse Particles, Phys. Rev. X 8, 031050 (2018).
  * Nishikawa  _et al._ [2022] Y. Nishikawa, M. Ozawa, A. Ikeda, P. Chaudhuri, and L. Berthier, Relaxation dynamics in the energy landscape of glass-forming liquids, Physical Review X 12, 10.1103/physrevx.12.021001 (2022).
  * Berthier and Reichman [2023] L. Berthier and D. R. Reichman, Modern computational studies of the glass transition, Nature Reviews Physics 5, 102 (2023).
  * Guiselin  _et al._ [2022] B. Guiselin, C. Scalliet, and L. Berthier, Microscopic origin of excess wings in relaxation spectra of supercooled liquids, Nature Physics 2022 , 1 (2022).
  * Scalliet  _et al._ [2022] C. Scalliet, B. Guiselin, and L. Berthier, Thirty milliseconds in the life of a supercooled liquid, Physical Review X 12, 041028 (2022).
  * Berthier  _et al._ [2019a] L. Berthier, P. Charbonneau, A. Ninarello, M. Ozawa, and S. Yaida, Zero-temperature glass transition in two dimensions, Nature communications 10, 1 (2019a).
  * Scalliet  _et al._ [2019] C. Scalliet, L. Berthier, and F. Zamponi, Nature of excitations and defects in structural glasses, Nat. Commun. 10, 5102 (2019).
  * Ciarella  _et al._ [2023] S. Ciarella, D. Khomenko, L. Berthier, F. C. Mocanu, D. R. Reichman, C. Scalliet, and F. Zamponi, Finding defects in glasses through machine learning, Nature Communications 14, 10.1038/s41467-023-39948-7 (2023).
  * Berthier  _et al._ [2019b] L. Berthier, M. Ozawa, and C. Scalliet, Configurational entropy of glass-forming liquids, The Journal of chemical physics 150, 160902 (2019b).
  * Berthier and Ediger [2020] L. Berthier and M. D. Ediger, How to ’measure’ a structural relaxation time that is too long to be measured?, The Journal of Chemical Physics 153, 44501 (2020).
  * Van Gurp and Palmen [1998] M. Van Gurp and J. Palmen, Time-temperature superposition for polymeric blends, Rheol. Bull 67, 5 (1998).
  * Berthier  _et al._ [2017] L. Berthier, P. Charbonneau, D. Coslovich, A. Ninarello, M. Ozawa, and S. Yaida, Configurational entropy measurements in extremely supercooled liquids that break the glass ceiling, Proc. Natl. Acad. Sci. U.S.A. 114, 11356 (2017).
  * Das and Sastry [2022] P. Das and S. Sastry, Crossover in dynamics in the Kob-Andersen binary mixture glass-forming liquid, Journal of Non-Crystalline Solids: X 14, 100098 (2022), arXiv:2204.04411v1 .
  * Mallamace  _et al._ [2010] F. Mallamace, C. Branca, C. Corsaro, N. Leone, J. Spooren, S.-H. Chen, and H. E. Stanley, Transport properties of glass-forming liquids suggest that dynamic crossover temperature is as important as the glass transition temperature, Proceedings of the National Academy of Sciences 107, 22457 (2010).
  * Pollak  _et al._ [1998] E. Pollak, H. Grabert, and P. Hänggi, Theory of activated rate processes for arbitrary frequency dependent friction: Solution of the turnover problem, The Journal of Chemical Physics 91, 4073 (1998).
  * Coslovich  _et al._ [2018] D. Coslovich, M. Ozawa, and W. Kob, Dynamic and thermodynamic crossover scenarios in the Kob-Andersen mixture: Insights from multi-CPU and multi-GPU simulations, The European Physical Journal E 41, 62 (2018).
  * Struik [1997] L. Struik, The apparent activation energy for mechanical and dielectric relaxation in glass-forming (polymeric) liquids: a misconception?, Polymer 38, 733 (1997).
  * Note [1] If $T^{\prime}$ is significantly larger than $T$, relaxation can occur via the nucleation of a hot liquid phase invading the material [49, 79].
  * Hecksher  _et al._ [2019] T. Hecksher, N. B. Olsen, and J. C. Dyre, Fast contribution to the activation energy of a glass-forming liquid, Proceedings of the National Academy of Sciences 116, 16736 (2019).
  * Eckert and Bartsch [2002] T. Eckert and E. Bartsch, Re-entrant glass transition in a colloid-polymer mixture with depletion attractions, Physical Review Letters 89, 125701 (2002).
  * Pham  _et al._ [2002] K. N. Pham, A. M. Puertas, J. Bergenholtz, S. U. Egelhaaf, A. Moussaïd, P. N. Pusey, A. B. Schofield, M. E. Cates, M. Fuchs, and W. C. K. Poon, Multiple glassy states in a simple model system, Science 296, 104 (2002).
  * Sciortino [2002] F. Sciortino, Disorderd materials: One liquid, two glasses, Nature Materials 1, 145 (2002).
  * Brito and Wyart [2009] C. Brito and M. Wyart, Geometric interpretation of previtrification in hard sphere liquids, The Journal of chemical physics 131, 149 (2009).
  * Degiuli  _et al._ [2015] E. Degiuli, E. Lerner, and M. Wyart, Theory of the jamming transition at finite temperature, J. Chem. Phys. 142, 164503 (2015).
  * Ghosh  _et al._ [2010] A. Ghosh, R. Mari, V. Chikkadi, P. Schall, J. Kurchan, and D. Bonn, Density of states of colloidal glasses and supercooled liquids, Soft matter 6, 3082 (2010).
  * Altieri [2018] A. Altieri, Higher-order corrections to the effective potential close to the jamming transition in the perceptron model, Physical Review E 97, 012103 (2018).
  * Altieri  _et al._ [2016] A. Altieri, S. Franz, and G. Parisi, The jamming transition in high dimension: an analytical study of the tap equations and the effective thermodynamic potential, Journal of Statistical Mechanics: Theory and Experiment 2016, 093301 (2016).
  * Pihlajamaa  _et al._ [2023] I. Pihlajamaa, C. C. L. Laudicina, and L. M. C. Janssen, Influence of polydispersity on the relaxation mechanisms of glassy liquids, Physical Review Research 5, 033120 (2023).
  * Meyer  _et al._ [2004] A. Meyer, J. Horbach, W. Kob, F. Kargl, and H. Schober, Channel Formation and Intermediate Range Order in Sodium Silicate Melts and Glasses, Physical Review Letters 93, 027801 (2004).
  * Voigtmann and Horbach [2009] T. Voigtmann and J. Horbach, Double Transition Scenario for Anomalous Diffusion in Glass-Forming Mixtures, Physical Review Letters 103, 205901 (2009).
  * Bartsch  _et al._ [2010] A. Bartsch, K. Rätzke, A. Meyer, and F. Faupel, Dynamic Arrest in Multicomponent Glass-Forming Alloys, Physical Review Letters 104, 195901 (2010).
  * Sastry and Austen Angell [2003] S. Sastry and C. Austen Angell, Liquid–liquid phase transition in supercooled silicon, Nature materials 2, 739 (2003).
  * Nicolas  _et al._ [2018] A. Nicolas, E. Ferrero, K. Martens, and J.-L. Barrat, Deformation and flow of amorphous solids: Insights from elastoplastic models, Rev. Mod. Phys. 90, 045006 (2018).
  * Popović  _et al._ [2021] M. Popović, T. de Geus, W. Ji, and M. Wyart, Thermally activated flow in models of amorphous solids, Phys. Rev. E 104, 025010 (2021).
  * Ozawa and Biroli [2022] M. Ozawa and G. Biroli, Elasticity, facilitation and dynamic heterogeneity in glass-forming liquids, arXiv preprint arXiv:2209.08861 (2022).
  * Tahaei  _et al._ [2023] A. Tahaei, G. Biroli, M. Ozawa, M. Popović, and M. Wyart, Scaling description of dynamical heterogeneity and avalanches of relaxation in glass-forming liquids, arXiv preprint arXiv:2305.00219 (2023).
  * Guiselin  _et al._ [2021] B. Guiselin, C. Scalliet, and L. Berthier, Microscopic origin of excess wings in relaxation spectra of deeply supercooled liquids, arXiv preprint 2103.01569 10.48550/arXiv.2103.01569 (2021).
  * Coslovich and Pastore [2007] D. Coslovich and G. Pastore, Understanding fragility in supercooled lennard-jones mixtures. i. locally preferred structures, The Journal of chemical physics 127, 124504 (2007).
  * Royall and Williams [2015] C. P. Royall and S. R. Williams, The role of local structure in dynamical arrest, Physics Reports 560, 1 (2015).
  * Mehri  _et al._ [2022] S. Mehri, L. Costigliola, and J. C. Dyre, Single-parameter aging in the weakly nonlinear limit, Thermo 2, 160 (2022).


# CONNECTION BETWEEN MID-INFRARED EMISSION PROPERTIES AND NARROW LINE REGION
OUTFLOWS IN TYPE 1 ACTIVE GALACTIC NUCLEI

Kai Zhang11affiliation: Key Laboratory for Research in Galaxies and Cosmology,
The University of Sciences and Technology of China, Chinese Academy of
Sciences, Hefei, Anhui 230026, China 2 2affiliationmark:  , Ting-Gui
Wang11affiliation: Key Laboratory for Research in Galaxies and Cosmology, The
University of Sciences and Technology of China, Chinese Academy of Sciences,
Hefei, Anhui 230026, China , Lin Yan33affiliation: Spitzer Science Center,
California Institute of Technology, 1200 E. California Blvd., Pasadena CA
91125, USA , and Xiao-Bo Dong11affiliation: Key Laboratory for Research in
Galaxies and Cosmology, The University of Sciences and Technology of China,
Chinese Academy of Sciences, Hefei, Anhui 230026, China<EMAIL_ADDRESS>
###### Abstract

The location of warm dust producing the Mid-infrared (MIR) emission in Type 1
Active Galactic Nuclei (AGNs) is complex and not yet fully known. We explore
this problem by studying how the MIR covering factor ($CF_{\rm MIR}$ = $L_{\rm
MIR}/L\mathrm{{}_{bol}}$) correlates with the fundamental parameters of AGN
accretion process (such as $L\mathrm{{}_{bol}}$, black hole mass
$M_{\mathrm{BH}}$, and Eddington ratio $L/L\mathrm{{}_{Edd}}$) and the
properties of narrow emission lines (as represented by [O III] $\lambda
5007$), using large data sets derived from the Sloan Digital Sky Spectroscopic
Survey (SDSS) and the Wide Infrared Sky Survey (WISE). Firstly we find that
the luminosity of the [O III] wing component ($L_{\rm wing}$) correlates more
tightly with the continuum luminosity ($\lambda L_{\lambda}(5100)$) than the
luminosity of the line core component ($L_{\rm core}$) does, which is in line
with our previous conclusion that the wing component, generally blueshifted,
originates from the polar outflows in the inner narrow-line region (NLR). We
then find that the MIR $CF$ shows the strongest correlation with $L_{\rm
wing}/L\mathrm{{}_{bol}}$ rather than with $L_{\rm core}/L\mathrm{{}_{bol}}$
or the above fundamental AGN parameters, and the correlation becomes stronger
as the infrared wavelength increases. We also confirm the anti-correlations of
$CF_{\rm MIR}$ with $L\mathrm{{}_{bol}}$ and $M_{\mathrm{BH}}$, and the lack
of dependence of $CF_{\rm MIR}$ on the Eddington ratio. These results suggest
that a large fraction of the warm dust producing MIR emission in AGNs is
likely embedded in polar outflows in the NLR instead of in the torus.

###### Subject headings:

galaxies: Seyfert— galaxies: kinematics and dynamics— infrared: galaxies—
quasars: emission lines

22affiliationtext: Key Laboratory for Research in Galaxies and Cosmology,
Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan
Road, Shanghai 200030, China

##  1\. Introduction

The mid-infrared (MIR) emission of radio quiet active galactic nuclei (AGNs)
may come from the dusty torus and/or narrow line region (NLR), as well as star
formation regions in the host galaxies. High-resolution imaging in the near-
infrared (NIR) provides direct evidence for the torus emission in a few nearby
AGNs (Bock et al. 2000; Radomski et al. 2003). At longer wavelengths,
observations showed extended MIR emission regions, coincident with the inner
NLR, in addition to the unresolved core, presumably the optically thick torus
(Tresch-Fienberg et al. 1987; Braatz et al. 1993; Cameron et al. 1993; Packham
et al. 2005; Hönig et al. 2010). Besides, MIR spectroscopic analysis revealed
an emission component, dominated by PAH emission, from star formation regions,
and its contribution increases with increasing wavelength and decreasing AGN
luminosity (e.g. Deo et al. 2009).

Recent high-resolution observations with the $MIDI$ instrument on $VLT$
demonstrated that in NGC1068, NGC 424, the Circinus galaxy and probably also
in NGC 3783, most of the MIR emission lies along the polar direction,
orthogonal to the major-axis of the NIR emission region (Raban et al. 2009;
Hönig et al. 2010, 2012; Tristram et al. 2012). Hönig et al. (2012) proposed
that the polar dust is entrained in radiatively driven dusty winds, originated
from the outerpart of the accretion disk (e.g. Königl & Kartje 1994; Keating
et al. 2012), or evaporated from the dusty torus (Dorodnitsyn et al. 2008).
Thus the strength of MIR emission may be closely connected to the properties
of the gas outflow. The outflow in the NLR is manifested by the blue-shifted
and relatively broad component in the high-ionization narrow emission lines
(e.g. Zhang et al. 2011). The NLR in nearby galaxies has been successfully
modeled by biconical outflows on scales of a few to tens of parsecs from the
nuclei in nearby Seyfert galaxies (Heckman et al. 1981; Das et al. 2005;
Crenshaw et al. 1999, 2010). It is well known that the NLR contains dust. It
would be interesting to explore if the MIR properties are linked to those of
narrow emission lines.

The Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) surveyed
the whole sky at 3.4, 4.6, 12 and 22 $\mu$m at sensitivities about 100-1000
better than previous all sky surveys. Taking the advantage of the WISE and the
Sloan Digital Sky Survey (SDSS) data releases, it is possible to construct a
large sample of AGNs with both measurements in MIR and optical to explore the
connection between the MIR properties and the AGN fundamental parameters or
NLR properties. In what follows, we adopt a $\Lambda$CDM cosmology with
$H_{\rm 0}$ = 70 km s-1 Mpc-1, $\Omega_{\rm m}$ = 0.3, and
$\Omega_{\rm\Lambda}$ = 0.7.

##  2\. The Sample and the Measurement

###  2.1. The Sample and Spectral Fit

From the spectroscopic data set of the SDSS Data Release 4 (Adelman-McCarthy
et al. 2006), we select 4178 Seyfert 1 galaxies and quasars (namely type 1
AGNs) that have reliable measurements of emission line fluxes and the
continuum suffers little host galaxy contamination. The sample is described in
detail in Dong et al. (2011). As described in Dong et al. (2008), we fit
simultaneously the continuum and emission lines in the spectral regimes near
H$\beta$+[O III] and H$\alpha$+[S II] with a model consisting of several
components: a power-law representing the featureless AGN continuum, the Fe II
templates, and a number of gaussians for emission lines. Each of the narrow
emission lines is fitted with a single Gaussian, except for the [O III]
$\lambda\lambda$4959, 5007 doublet lines, which are modeled with two
Gaussians. One accounts for the line core and the other for a possible blue
wing. The details of the fitting algorithm and spectral parameters can be
found in Dong et al. (2011) and Zhang et al. (2011).

###  2.2. Mid-Infared Data

We cross-correlate our SDSS Type 1 AGN catalog with the WISE All Sky Data
Release (Wright et al. 2010) catalog using a match radius of 3 arcseconds to
construct the mid-infrared detected type 1 AGN sample. Nearly all sources in
our sample are detected by WISE (4174 out of 4178). We refer to the WISE Vega
magnitude as $W1$, $W2$, $W3$ and $W4$ for the 3.4, 4.6, 12, and 22$\mu$m
bands, respectively. We reject the sources that have signal to noise ratio
less than 3 (S$/$N$<3$) in any WISE band. Radio loud AGNs are rejected
according to the catalog by Lu et al. (2010) to avoid contamination to [O III]
flux and continuum flux. The final sample has 3336 sources, and we call it
’extended sample’ hereafter. We further make a redshift cut of z$<$0.3 to
minimize potential cosmological evolution and uncertainty in the
$k$-correction that might complicate the correlation analysis. This ’base
sample’ consists of 608 sources. The magnitude in each band is converted to
$\nu F_{\nu}$ according to the published WISE zero points, and the
$k$-correction is made by interpolating the Spectral Energy Distribution (SED)
at rest-frame wavelengths of 3.4, 4.6, 12, and 17$\mu$m. The $k$-corrected
luminosities in each band are marked as $L1$, $L2$, $L3$ and $L4$. For every
source we connect the adjacent bands on the $\log\lambda F_{\lambda}$ versus
$\log\lambda$ diagram to approximate the MIR SED, and integrate the SED from
3.4 to 17$\mu$m to obtain the total MIR luminosity, denoted as $L_{\rm WISE}$.

##  3\. Distinction between the [O III] Core and Wing Component

Zhang et al. (2011) found that the equivalent width (EW) of the core component
of [O III] correlates negatively with $\lambda L_{\lambda}(5100)$ (Baldwin
effect) while the wing component displays no or a weak anti-Baldwin effect.
The different behaviors suggest that the two components come from distinct gas
components regulated by different physical processes. The gaseous kinematics
of the [O III] core traces the motion of stars in the gravitational potential
well of bulges (e.g. Nelson et al. 2000; Wang & Lu 2001; Komossa et al. 2008),
indicating that the [O III] core is emitted by the photo-ionized ISM in the
galaxy bulge. This is supported by imaging spectroscopic observations of
nearby Seyfert galaxies, which show that the gaseous and stellar kinematics
are coupled (e.g. Dumas et al. 2008). On the other hand, the wing is much
broader than the core but narrower than broad lines, indicating that its
emitting region is much closer to the nuclei than the core component.
Furthermore, the wings are usually blueshifted up to 1000 km s-1, suggesting
they are from an outflow component. High spatial resolution spectroscopic
observations of nearby Seyfert galaxies revealed individual high velocity
knots up to a velocity of 3000 km s-1 on scales from parsecs to up a hundred
parsecs (e.g. Crenshaw et al 2010). All these observations suggest that the
wing component is formed in the outflow at a few to a few hundred parsecs from
the nuclei.

To examine the difference between the core and the wing in detail, we plot the
luminosities of the [O III] core and wing against $\lambda L_{\lambda}(5100)$
in Figure 1 using the extended sample. Obviously, both the core and the wing
luminosities correlate with $\lambda L_{\lambda}(5100)$. The best fit of a
log-linear model yields a slope of 0.86$\pm$0.01 for the core component, and
1.14$\pm$0.01 for the wing component. These slopes are consistent with the
analysis of Zhang et al (2011) on the Baldwin effect. However, the plot
reveals an interesting feature that the correlation is much tighter for the
wing than for the core. Qualitatively, the $RMS$ scatters to the best fit of
the log-linear model are 0.29 dex for the wing, and 0.42 dex for the core. The
difference is significant according to $F-$test with a chance probability
$P<10^{-26}$. The $L_{\rm wing}$ correlates more tightly with the the broad
H$\alpha$ luminosity, which used for bolometric luminosity estimation in next
section, than the $L_{\rm core}$ does for the basic sample. Note that the
decomposition into the wing and the core, and other measurements introduce an
uncertainty of about 0.15 dex in the EWs of the wing and the core.111For every
source, we use best-fit 2-Gaussian model of the [O III] $\lambda 5007$ line
plus the measurement errors of the SDSS spectrum to generate a new mock
spectrum, then fit the mock [O III] line with two Gaussians. Statistically in
the extended sample, the standard deviations of the difference between the
output and input fluxes are 0.12 and 0.15 dex for the core and wing
components, respectively.

The tighter correlation between the wing and the continuum luminosity
indicates that outflows do not change dramatically from object to object,
while the larger scatter in the core component suggests that the variations of
ISM states like density, density profile, ionization state etc., are
pronounced. Physically, the ISM gas in the bulge is regulated by several
physical processes, such as gas accretion, star-formation, stellar wind and
AGN feedback. The competition between different processes may lead to very
diverse physical conditions of the gas among AGNs, which introduces a scatter
in the correlation. Interestingly, Zhang et al. (2011) found that the EW of
the [O III] core correlates strongest with the outflow velocity of the whole
[O III] profile (the outflow velocities being negative due to the outflow,
hence the correlation), suggesting the interplay of the outflow and ISM of the
galaxy bulge. This can be interpreted most directly in terms of AGN feedback:
high velocity outflows blow away the ISM in the host bulge. Also, continuum
variations may introduce a larger scatter in the core component than in the
wing component because the strength of the core component depends on the long-
term variability while the strength of the wing component follows more closely
the continuum variations (Zhang et al. 2013). The effect of the variability
depends on the light curves over time scales of $10^{3}$ years.

##  4\. Correlation between $CF_{\rm MIR}$ and AGN properties

We use the ratios of mid-infrared to bolometric luminosities to quantify the
MIR $CFs$, i.e. the relative contribution of MIR to bolometric luminosity. The
bolometric luminosity is estimated from the broad line H$\alpha$ luminosity
using a bolometric correction factor of 170, which is combining
$L\mathrm{{}_{bol}}$=12.17$\times$$\lambda L_{\lambda}(5100)$ in Richards et
al. (2006) and $\lambda L_{\lambda}(5100)$=14$\times L_{\rm H\alpha}$ in Stern
& Laor (2012).222The bolometric correction may depend on the accretion rate
and luminosity (e.g. Vasudevan & Fabian 2007; Lusso et al. 2012). This is
similar to the $130_{\div 2.4}^{\times 2.4}$ adopted in Stern & Laor (2012).
We do not use $\lambda L_{\lambda}(5100)$ directly as an estimate for the
bolometric luminosity, because it may be subject to host contamination and the
optical-UV slope gets bluer in higher luminosity AGN (See Section 4.1 in Zhang
et al. (2013) for details). The median bolometric luminosity of our sample is
$10^{45}$ erg s-1. In Figure 2 we plot the histograms of mid-infrared
luminosity to bolometric luminosity for the individual wavelengths (Panels a
to d) and for the total MIR luminosity (Panel e). All histograms show a log-
normal distribution. For $CF_{\rm WISE}$, the distribution is centered at log
$CF_{\rm WISE}$ = -0.37 (43%) and has a dispersion of 0.243 dex. To first
order, this means a covering factor of 43%, and is consistent with the type 2
fraction of $\sim$50% at this luminosity in literature (e.g. Hao et al. 2005;
Lawrence & Elvis 2010). However, one needs to caution that the type 2 fraction
may not be equal to the covering factor of warm dust. Anisotropy in primary
continuum and IR emission, presence of large scale obscuring material and hot
dust that emits at wavelength shorter than 3.4$\mu$m lead to deviation of one
from another (Calderone et al. 2012; Roseboom et al. 2013; Ma & Wang 2013).

The $CF_{\rm MIR}$ spreads over a dynamic range (the range it spans) of more
than 1.8 dex at the 17$\mu$m band and the dynamic range increases with
increasing wavelength. This trend is clearly seen in Figure 2 in panel (a) to
(d) respectively. The $\sigma$’s of $CF_{\rm MIR}$ increase from 0.200 dex at
3.4$\mu$m to 0.298 dex at 17$\mu$m. One possible reason of increasing
dispersion with wavelength is that star-formation contributes a large fraction
of emission at longer wavelengths. As discussed at the beginning of Section 5,
star-formation would contribute less than 10% at 17$\mu$m, thus introducing a
dispersion less than 0.04 dex. This alone is not enough to account for the
0.22 dex difference between dispersions of $CF1$ and $CF4$. The intrinsic
variance after subtracting the measurement errors in IR luminosity still show
a difference of 0.20 dex. There are two possibilities for this difference: 1)
The 17$\mu$m emission is from the dust torus but the fraction of warm dust in
the torus varies drastically from object to object due to either different
degree of clumpiness or different optical depth of the torus. 2) There is an
additional warm dust component, and the contribution from the additional
component varies. As discussed in Section 5.1, the torus explanation is
inconsistent with our findings, and we argue below that the latter is
plausible.

###  4.1. Correlations between $CF_{\rm MIR}$ and fundamental AGN properties

First, we analyze how $CF_{\rm MIR}$ correlate with redshift and the AGN
fundamental parameters, such as $L\mathrm{{}_{bol}}$, $M_{\mathrm{BH}}$ and
$L/L\mathrm{{}_{Edd}}$ using the basic sample. The results are summarized in
column 1, 2, 3, 4 of Table 1. We calculate the black hole masses based on the
FWHM of broad H$\beta$ and $\lambda L_{\lambda}(5100)$ using the formalism
presented in Wang et al. (2009, their Eqn. 11). The Eddington ratio
($L/L\mathrm{{}_{Edd}}$) is the ratios between the bolometric and Eddington
luminosities:
$L/L\mathrm{{}_{Edd}}=\frac{L\mathrm{{}_{bol}}}{L\mathrm{{}_{Edd}}}$,
$L\mathrm{{}_{Edd}}=1.26\times 10^{38}$ ($M_{\mathrm{BH}}$/$M_{\odot}$) erg
s-1. We use Spearman correlation coefficient ($\rho$) as a non-parametric
measure of the statistical dependence between two quantities.

We find that the $CF_{\rm MIR}$ correlate negatively with $L\mathrm{{}_{bol}}$
and $M_{\mathrm{BH}}$, but very weakly with $L/L\mathrm{{}_{Edd}}$. We get
$CF1\propto L\mathrm{{}_{bol}}^{-0.21}$, $CF2\propto
L\mathrm{{}_{bol}}^{-0.21}$, $CF3\propto L\mathrm{{}_{bol}}^{-0.31}$,
$CF4\propto L\mathrm{{}_{bol}}^{-0.35}$, and $CF_{\rm WISE}/\propto
L\mathrm{{}_{bol}}^{-0.26}$ respectively as shown in Figure 3. The typical
error of the slopes is 0.03. These results are fully consistent with previous
works (Maiolino et al. 2007; Treister et al. 2008; Calderone et al. 2012; Ma &
Wang 2013; Roseboom et al. 2013). We plot the best regression fits in
Calderone et al. (2012) and Ma & Wang (2013) on Panel (e) to illustrate the
consistency. Calderone et al. (2012) integrated the MIR luminosity at
2.1-13.7$\mu$m, and Ma & Wang (2013) used 3-10$\mu$m luminosity. To make a
fair comparison, we add 0.184 dex to the $CF_{\rm MIR}$ in Ma & Wang (2013) to
account for the shorter wavelength range. The 0.184 dex is the mean ratio of
3.4-17$\mu$m luminosity to 3.4-11.6$\mu$m luminosity in our sample. We do not
make any correction to the regression fit in Calderone et al. (2012) because
the gain in 2-3.4$\mu$m approximately equals to the loss in 13.7-17$\mu$m
according to the mean SED of all quasars in Richards et al. (2006). We can see
that the fitting lines from literatures and our best regression fit are
consistent with each other at this luminosity range. This further supports the
bolometric correction we adopt is reasonable.

The slopes of $CF_{\rm MIR}$ versus $M_{\mathrm{BH}}$ relations are -0.22,
-0.23, -0.35, -0.42, -0.30 for L1, L2, L3, L4, and $L_{\rm WISE}$
respectively, consistent with previous works (e.g. Mor et al. 2011; Ma & Wang
2013). The lack of dependence of $CF_{\rm MIR}$ on $L/L\mathrm{{}_{Edd}}$ is
reported by many authors (Cao et al. 2005; Mor et al. 2011; Ma & Wang 2013).
Considering the limited range of $L/L\mathrm{{}_{Edd}}$(1.5 dex) for our
sample, we refrain from over-interpreting this result.

###  4.2. Correlations between $CF_{\rm MIR}$ and narrow emission line
properties

Next, we explore whether the MIR emission is linked to the properties of
optical narrow emission lines as represented by [O III] $\lambda 5007$. We
present the results of correlation analysis between $CF_{\rm MIR}$ and the [O
III] total, core, wing luminosity to $L\mathrm{{}_{bol}}$ ratios using the
basic sample. The results are summarized in column 5, 6, 7 of Table 1. The
$CF_{\rm MIR}$ in all 4 bands correlate significantly with $L_{\rm
wing}$/$L\mathrm{{}_{bol}}$. Among them, the strongest one is between $CF4$
and $L_{\rm wing}$/$L\mathrm{{}_{bol}}$ with $\rho=0.52$. The correlation
weakens as wavelength decreases: with $\rho=$0.47, 0.36, 0.32 between $L_{\rm
wing}$/$L\mathrm{{}_{bol}}$ and $CF3$, $CF2$ and $CF1$, respectively. This is
in parallel with the increasing dynamic range of the $CFs$ with wavelengths.
Note that the bolometric correction has $RMS$ scatters of order 0.38 dex. So
it is likely that the weakening of the correlations is a combination of
decreasing dynamic range and relatively large uncertainty in the bolometric
correction. The physical reason of increasing dynamic range with increasing
wavelength is discussed in Section 5.3. We plot $CF4$ versus $L_{\rm
wing}$/$L\mathrm{{}_{bol}}$ in Figure 4. A log linear fit yields a slope of
1.19$\pm$0.04 with an $RMS$ scatter of 0.308 dex. This scatter is
significantly larger than the combinations of measurement errors: 0.11 dex in
L4/$L_{\rm H\alpha}$ and 0.15 dex in $L_{\rm wing}$/$L_{\rm H\alpha}$,
indicating that there is a considerable intrinsic scatter in this relation. No
correlation is found between outflow velocity and any IR properties in our
sample.

##  5\. Discussion

To facilitate discussions on the implication of the correlation between
$L_{\rm wing}$/$L\mathrm{{}_{bol}}$ and L4/$L\mathrm{{}_{bol}}$, we note that
the star-formation contamination is likely small in the MIR band for these
luminous type 1 AGNs in our sample. For most of our sources, the infrared
luminosity at 17$\mu$m is in the range of $10^{43}$ to $10^{46}$ erg s-1
($10^{9.5}$ to $10^{12.5}L_{\bigodot}$). According to Fig 10 of Deo et al.
(2009), starburst contribution is likely less than 10% at $17\mu$m and a few
percent in the L1 band. This is more likely to be true because we have removed
any source with substantial stellar contamination in optical. In the following
sections, we will discuss the location of the warm dust and its relation with
torus in light of the above results.

###  5.1. Classical Steady Torus?

Judging from studies of AGN at subarcsecond resolution in the MIR, the
majority of MIR emission originates from the central few pc region (Jaffe et
al. 2004; Tristram et al. 2007). In NGC 1068, about 75% of the MIR emission is
from the central unresolved core (Bock et al. 2000). This was believed to be
the torus, and a variety of torus models were constructed to match the NIR–MIR
SED (e.g. Fritz et al. 2006; Hönig et al. 2006; Schartmann et al. 2008;
Nenkova et al. 2008a, b). If the dust dominates the MIR emission, as the torus
acts to confine the ionization cone, there should be an anti-correlation
between the torus subtending angle and the covering factor of NLR; we thus
should obtain a negative correlation between $CF4$ and $L_{\rm
core}$/$L\mathrm{{}_{bol}}$. This is however in contradiction to the
observation.

###  5.2. Dusty NLR clouds?

NLR is theoretically expected to contain dust (Netzer & Laor 1993; Dopita et
al. 2002). Observationally, extended polar MIR emission is found around the
compact cores of several AGNs (See Section 1). Groves et al. (2006) modeled
the narrow line spectrum of NGC 4151 with dusty photo-ionization models, and
found that 10–50% of its 25$\mu$m emission is from the NLR dust. The imaging
observations also indicate that $>$27% of MIR emission emerges from a region
coincident with the NLR in NGC 4151 (Radomski et al. 2003). The warm dust
extends from the unresolved core to about 100pc. Mor et al. (2009) fitted the
MIR spectra of PG QSOs with 3 components: a black body, the clumpy torus
emission (Nenkova et al. 2008a,b) and the NLR dust. They derived that above
18$\mu$m, the NLR dust contributes as much as 40% of the emission. Schweitzer
et al. (2008) demonstrated the the NLR dust could produce the silicate
emission features observed in AGNs. As the cold gas and the dust are coupled
in the NLR, this would naturally explain the correlation between $CF4$ and
$L_{\rm wing}$/$L\mathrm{{}_{bol}}$ and the correlation weakens at shorter
infrared wavebands, qualitatively.

However, the NLR dust models do not consider the dust in the outflow component
of NLR, so they can not explain why $CF4$ is much more strongly correlated
with $L_{\rm wing}$/$L\mathrm{{}_{bol}}$ than with $L_{\rm
core}$/$L\mathrm{{}_{bol}}$. Furthermore, most of the infrared emission from 5
to 18$\mu$m is attributed to the torus emission, which is not consistent with
more recent higher resolution imaging observations in several nearby Seyfert
galaxies (see §1). In Mor et al. (2009), the distance of NLR dust is estimated
to be $R_{NLR}=295\times L_{46}^{0.5}$ pc, where $L_{46}$ is the bolometric
luminosity in unit of $10^{46}$ erg s-1. In high luminosity sources, we lack
observational evidence that the outflow may reach that far from the nuclei.
Although detailed model is needed to predict the infrared spectrum, it is
clear that any optical thin dust in the outflow region will be hotter than
that in the more distant classical NLR, and will contribute significantly to
emission below 15 $\mu$m. Alternatively, the dust clouds in outflows are
optically thick to infrared light, radiative transfer effect would result in
much colder infrared emission. We therefore do not think that dusty NLR clouds
on scales of $\gtrsim$100 pc are responsible for our observed correlation of
$L_{\rm wing}/L\mathrm{{}_{bol}}$ with $CF_{\rm MIR}$.

###  5.3. Dusty Outflows?

It was pointed out by Krolik and Begelman (1988) that a steady doughnut torus
suffers the problem that there is not enough vertical motions capable of
sustaining the H/R$\sim$1 geometry. One solution to this problem is assuming
that the dust is entrained in the magnetically driven accretion disk wind
(Emmering, Blandford & Shlosman 1992; Elitzur & Shlosman 2006), and the outer
part of the outflow acts as a clumpy torus. Alternatively, the IR radiation
pressure on dust may also offer enough pressure to support a thick torus, even
launches an outflow (Krolik 2007; Dorodnitsyn et al. 2011, 2012). In this type
of dynamical models, the outflow and the torus are naturally linked. A
connection between the outflow and the dust emission is expected. The
detection of dominant MIR emission in polar direction led Hönig et al. (2012)
to propose that the polar dust is entrained in a radiatively-driven dusty
winds. This dust component is an inward extent of NLR dust observed in many
nearby AGNs but lies only a few pc away from the nuclei. This is where the
outflow originates and accelerated. In this paradigm, the outflow gas and the
warm dust share similar covering factor, and would lead to our finding of a
positive correlation between $CF4$ and $L_{\rm wing}$/$L\mathrm{{}_{bol}}$
naturally. The contribution of the polar warm dust to the mid-infrared
emission increases with wavelength, so the correlation between $CF_{\rm MIR}$
and $L_{\rm wing}$/$L\mathrm{{}_{bol}}$ strengthens with wavelength. By
comparison, the static classical NLR traced by the [O III] core component lies
further from the nuclei. The dust embedded in the static NLR emit dominantly
at wavelength longer than 20$\mu$m, thus only connects to the warm dust
loosely. If the correlation between the warm dust and the outflow strength is
really driven by the dust in the polar winds, an implication is that the pc
scale MIR elongation in NGC 1068, the Circinus galaxy and NGC 424 may not be
special cases, but a common phenomenon for Seyfert galaxies. Most of the warm
dust may reside in the outflow instead of in the torus. Besides, the outflow
may have a significant effect on the dust re-distribution by transporting the
dust from the central region to larger scales.

We thank the referee for helpful suggestions that improve the paper
significantly. This work is supported by Chinese NSF grants NSFC 11233002,
11073019, and a Chinese Universities Scientific Fund (CUSF) for USTC. This
publication makes use of data products from the Wide-field Infrared Survey
Explorer, which is a joint project of the University of California, Los
Angeles, and the Jet Propulsion Laboratory/California Institute of Technology,
funded by the National Aeronautics and Space Administration. This study makes
use of data from the SDSS (see:
http://www.sdss.org/collaboration/credits.html).

## References

  * Adelman-McCarthy et al. (2006) Adelman-McCarthy, J. K., Agüeros, M. A., Allam, S. S., et al. 2006, ApJS, 162, 38
  * Bock et al. (2000) Bock, J. J., Neugebauer, G., Matthews, K., et al. 2000, AJ, 120, 2904
  * Braatz et al. (1993) Braatz, J. A., Wilson, A. S., Gezari, D. Y., Varosi, F., & Beichman, C. A. 1993, ApJ, 409, L5
  * Calderone et al. (2012) Calderone, G., Sbarrato, T., & Ghisellini, G. 2012, MNRAS, 425, L41
  * Cao (2005) Cao, X. 2005, ApJ, 619, 86
  * Cameron et al. (1993) Cameron, M., Storey, J. W. V., Rotaciuc, V., et al. 1993, ApJ, 419, 136
  * Crenshaw et al. (1999) Crenshaw, D. M., Kraemer, S. B., Boggess, A., et al. 1999, ApJ, 516, 750
  * Crenshaw et al. (2010) Crenshaw, D. M., Schmitt, H. R., Kraemer, S. B., Mushotzky, R. F., & Dunn, J. P. 2010, ApJ, 708, 419
  * Das et al. (2005) Das, V., Crenshaw, D. M., Hutchings, J. B., et al. 2005, AJ, 130, 945
  * Deo et al. (2009) Deo, R. P., Richards, G. T., Crenshaw, D. M., & Kraemer, S. B. 2009, ApJ, 705, 14
  * Dong et al. (2008) Dong, X., Wang, T., Wang, J., et al. 2008, MNRAS, 383, 581
  * Dong et al. (2011) Dong, X.-B., Wang, J.-G., Ho, L. C., et al. 2011, ApJ, 736, 86
  * Dopita et al. (2002) Dopita, M. A., Groves, B. A., Sutherland, R. S., Binette, L., & Cecil, G. 2002, ApJ, 572, 753
  * Dorodnitsyn et al. (2008) Dorodnitsyn, A., Kallman, T., & Proga, D. 2008, ApJ, 687, 97
  * Dorodnitsyn et al. (2011) Dorodnitsyn, A., Bisnovatyi-Kogan, G. S., & Kallman, T. 2011, ApJ, 741, 29
  * Dorodnitsyn & Kallman (2012) Dorodnitsyn, A., & Kallman, T. 2012, ApJ, 761, 70
  * Dumas et al. (2008) Dumas, G., Mundell, C. G., Emsellem, E., & Nagar, N. 2008, Astronomische Nachrichten, 329, 908
  * Elvis et al. (1994) Elvis, M., Wilkes, B. J., McDowell, J. C., et al. 1994, ApJS, 95, 1
  * Elitzur & Shlosman (2006) Elitzur, M., & Shlosman, I. 2006, ApJ, 648, L101
  * Emmering et al. (1992) Emmering, R. T., Blandford, R. D., & Shlosman, I. 1992, ApJ, 385, 460
  * Fritz et al. (2006) Fritz, J., Franceschini, A., & Hatziminaoglou, E. 2006, MNRAS, 366, 767
  * Groves et al. (2006) Groves, B., Dopita, M., & Sutherland, R. 2006, A&A, 458, 405
  * Hao et al. (2005) Hao, L., Strauss, M. A., Fan, X., et al. 2005, AJ, 129, 1795
  * Heckman et al. (1981) Heckman, T. M., Miley, G. K., van Breugel, W. J. M., & Butcher, H. R. 1981, ApJ, 247, 403
  * Hönig et al. (2006) Hönig, S. F., Beckert, T., Ohnaka, K., & Weigelt, G. 2006, A&A, 452, 459
  * Hönig et al. (2010) Hönig, S. F., Kishimoto, M., Gandhi, P., et al. 2010, A&A, 515, A23
  * Hönig et al. (2012) Hönig, S. F., Kishimoto, M., Antonucci, R., et al. 2012, arXiv:1206.4307
  * Jaffe et al. (2004) Jaffe, W., Meisenheimer, K., Röttgering, H. J. A., et al. 2004, Nature, 429, 47
  * Keating et al. (2012) Keating, S. K., Everett, J. E., Gallagher, S. C., & Deo, R. P. 2012, ApJ, 749, 32
  * Komossa et al. (2008) Komossa, S., Xu, D., Zhou, H., Storchi-Bergmann, T., & Binette, L. 2008, ApJ, 680, 926
  * Konigl& Kartje (1994) Konigl, A., & Kartje, J. F. 1994, ApJ, 434, 446
  * Krolik & Begelman (1988) Krolik, J. H., & Begelman, M. C. 1988, ApJ, 329, 702
  * Krolik (2007) Krolik, J. H. 2007, ApJ, 661, 52
  * Lawrence & Elvis (2010) Lawrence, A., & Elvis, M. 2010, ApJ, 714, 561
  * Lu et al. (2010) Lu, Y., Wang, T.-G., Dong, X.-B., & Zhou, H.-Y. 2010, MNRAS, 404, 1761
  * Lusso et al. (2012) Lusso, E., Comastri, A., Simmons, B. D., et al. 2012, MNRAS, 425, 623
  * Ma & Wang (2013) Ma, X.-C., & Wang, T.-G. 2012, arXiv:1212.4245
  * Maiolino et al. (2007) Maiolino, R., Shemmer, O., Imanishi, M., et al. 2007, A&A, 468, 979
  * Mor et al. (2009) Mor, R., Netzer, H., & Elitzur, M. 2009, ApJ, 705, 298
  * Mor & Trakhtenbrot (2011) Mor, R., & Trakhtenbrot, B. 2011, ApJ, 737, L36
  * Nelson (2000) Nelson, C. H. 2000, ApJ, 544, L91
  * Nenkova et al. (2008) Nenkova, M., Sirocky, M. M., Ivezić, Ž., & Elitzur, M. 2008, ApJ, 685, 147
  * Nenkova et al. (2008) Nenkova, M., Sirocky, M. M., Nikutta, R., Ivezić, Ž., & Elitzur, M. 2008, ApJ, 685, 160
  * Netzer & Laor (1993) Netzer, H., & Laor, A. 1993, ApJ, 404, L51
  * Packham et al. (2005) Packham, C., Radomski, J. T., Roche, P. F., et al. 2005, ApJ, 618, L17
  * Raban et al. (2009) Raban, D., Jaffe, W., Röttgering, H., Meisenheimer, K., & Tristram, K. R. W. 2009, MNRAS, 394, 1325
  * Radomski et al. (2003) Radomski, J. T., Piña, R. K., Packham, C., et al. 2003, ApJ, 587, 117
  * Richards et al. (2006) Richards, G. T., Lacy, M., Storrie-Lombardi, L. J., et al. 2006, ApJS, 166, 470
  * Roseboom et al. (2013) Roseboom, I. G., Lawrence, A., Elvis, M., et al. 2013, MNRAS, 429, 1494
  * Schartmann et al. (2008) Schartmann, M., Meisenheimer, K., Camenzind, M., et al. 2008, A&A, 482, 67
  * Schmitt et al. (2003) Schmitt, H. R., Donley, J. L., Antonucci, R. R. J., Hutchings, J. B., & Kinney, A. L. 2003, ApJS, 148, 327
  * Schweitzer et al. (2008) Schweitzer, M., Groves, B., Netzer, H., et al. 2008, ApJ, 679, 101
  * Stern & Laor (2012) Stern, J., & Laor, A. 2012, MNRAS, 423, 600
  * Treister et al. (2008) Treister, E., Krolik, J. H., & Dullemond, C. 2008, ApJ, 679, 140
  * Tresch-Fienberg et al. (1987) Tresch-Fienberg, R., Fazio, G. G., Gezari, D. Y., et al. 1987, ApJ, 312, 542
  * Tristram et al. (2007) Tristram, K. R. W., Meisenheimer, K., Jaffe, W., et al. 2007, A&A, 474, 837
  * Tristram et al. (2012) Tristram, K. R. W., Schartmann, M., Burtscher, L., et al. 2012, arXiv:1204.2823
  * Vasudevan & Fabian (2007) Vasudevan, R. V., & Fabian, A. C. 2007, MNRAS, 381, 1235
  * Wang & Lu (2001) Wang, T., & Lu, Y. 2001, A&A, 377, 52
  * Wang et al. (2009) Wang, J.-G., Dong, X.-B., Wang, T.-G., et al. 2009, ApJ, 707, 1334
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
  * Zhang et al. (2011) Zhang, K., Dong, X.-B., Wang, T.-G., & Gaskell, C. M. 2011, ApJ, 737, 71
  * Zhang et al. (2013) Zhang, K., Wang, T.-G., Gaskell, C. M., & Dong, X.-B. 2013, ApJ, 762, 51

Table 1Correlation Tables between the IR covering factors and the optical
parameters aafootnotemark:

| z | $L\mathrm{{}_{bol}}$ | $M_{\mathrm{BH}}$ | $L/L\mathrm{{}_{Edd}}$ |  $L_{[O\,{\tiny III}]}$/$L\mathrm{{}_{bol}}$ |  $L_{\rm core}$/$L\mathrm{{}_{bol}}$ |  $L_{\rm wing}$/$L\mathrm{{}_{bol}}$
---|---|---|---|---|---|---|---
$CF_{\rm WISE}$ | 0.06(1e-01) | -0.40(3e-25) | -0.38(2e-22) | -0.10(7e-03) | 0.28(1e-12) | 0.09(2e-02) | 0.45(1e-32)
$CF1$ | 0.08(2e-02) | -0.37(5e-22) | -0.32(8e-17) | -0.12(1e-03) | 0.10(1e-02) | -0.07(6e-02) | 0.32(3e-16)
$CF2$ | 0.09(1e-02) | -0.35(4e-20) | -0.33(2e-17) | -0.09(1e-02) | 0.14(5e-04) | -0.04(2e-01) | 0.36(8e-20)
$CF3$ | 0.05(1e-01) | -0.40(4e-25) | -0.38(5e-23) | -0.10(1e-02) | 0.31(2e-15) | 0.12(1e-03) | 0.47(4e-35)
$CF4$ | 0.02(5e-01) | -0.42(2e-27) | -0.40(3e-25) | -0.10(6e-03) | 0.43(9e-29) | 0.25(9e-11) |  0.52(1e-43)

a For each entry, we list the Spearman rank correlation coefficient ($\rho$)
and the probability of the null hypothesis that the correlation is not present
($P_{\mathrm{null}}$) in parenthesis, for the 608 objects in the basic
sample.$\rho\geq$0.50 are marked in Boldface.

Figure 1.— Left: The [O III] core luminosity versus $\lambda
L_{\lambda}(5100)$ using the extended sample (3336 sources). The solid line
shows the best fitting line, and the cyan contours show the density of the
distribution. The levels of the contours are 15, 30, 45, 60 and 75
respectively. Right: The [O III] wing luminosity against $\lambda
L_{\lambda}(5100)$ using the extended sample. The levels of the contours are
15, 30, 60 and 80.  Figure 2.— Distributions of MIR covering factors: $CF_{\rm
MIR}$ = $L_{\rm MIR}/L\mathrm{{}_{bol}}$ using the basic sample. The solid
line is the gaussian fit to the distribution, the centroid and $\sigma$ of the
gaussian are shown in the left-up corner. Panel (a) to (d) are $CF_{\rm MIR}$
distributions at each band, and Panel (e) is the distribution of $CF_{\rm
WISE}=L_{\rm WISE}$/$L\mathrm{{}_{bol}}$.  Figure 3.— $CF_{\rm MIR}$ against
$L\mathrm{{}_{bol}}$ of the basic sample. The blue line is the best regression
fit, and the red line is the mean value in each $L\mathrm{{}_{bol}}$ bin with
binsize of 0.3 dex starting from $L\mathrm{{}_{bol}}=10^{43.8}erg$ $s^{-1}$.
The slope of the regression fit is shown in the left-down corner. Panel (a) to
(d) are $CF_{\rm MIR}$ versus $L\mathrm{{}_{bol}}$ at each band, and Panel (e)
is the $CF_{\rm WISE}$ versus the bolometric luminosity. The green and purple
lines are the results from Ma & Wang (2013) and Calderone et al. (2012).
Figure 4.— $CF4$ versus $L_{\rm wing}/L\mathrm{{}_{bol}}$ using the basic
sample (608 sources). The correlation coefficient is shown in the left-up
corner. The solid line is the best regression fit.