# Bringing together two paradigms of non-equilibrium:
Driven dynamics of aging systems

Diego Tapias<EMAIL_ADDRESS>Institute for
Theoretical Physics, Georg-August-Universität Göttingen, 37077 Göttingen,
Germany  Charles Marteau  Department of Mathematics, King’s College London,
London WC2R 2LS, UK  Department of Physics and Astronomy, University of
British Columbia, Vancouver V6T 1Z1, Canada  Fabián Aguirre-López  Department
of Mathematics, King’s College London, London WC2R 2LS, UK  Chair of
Econophysics and Complex Systems, École polytechnique, 91128 Palaiseau Cedex,
France  LadHyX UMR CNRS 7646, École polytechnique, 91128 Palaiseau Cedex,
France  Peter Sollich<EMAIL_ADDRESS>Institute for
Theoretical Physics, Georg-August-Universität Göttingen, 37077 Göttingen,
Germany  Department of Mathematics, King’s College London, London WC2R 2LS, UK

###### Abstract

There are two fundamental paradigms for non-equilibrium dynamics: on the one
hand, aging towards an equilibrium state that cannot be reached on reasonable
timescales; on the other, external driving that can lead to non-equilibrium
steady states. We explore how these two mechanisms interact, by studying the
behaviour of trap models, which are paradigmatic descriptions of slow glassy
dynamics, when driven by trajectory bias towards high or low activity. To
diagnose whether the driven systems continue to age, we establish a framework
for mapping the biased dynamics to a Markovian time evolution with time-
dependent transition rates. We find that the original aging dynamics reacts in
two qualitatively distinct ways to the driving: it can be destroyed by driving
of any nonzero strength (_fragile_ aging), whereby the dynamics either reaches
an active steady state or effectively freezes; or it can persist within a
finite range of driving strengths around the undriven case (_robust_ aging).
This classification into fragile and robust aging could form the basis for
distinguishing different universality classes of aging dynamics.

In non-equilibrium statistical physics there are two fundamental paradigms for
bringing a system out of equilibrium. The first consists in slowing the
dynamics down to such an extent that the system fails to reach equilibrium on
practical timescales and instead ages, with its properties depending on the
time elapsed since preparation. This phenomenon is found in glasses [1, 2, 3,
4, 5, 6, 7, 8, 9] after abrupt decreases in temperature or volume, beyond the
regime of the supercooled (metastable equilibrium) liquid, but is in fact much
more widespread, occurring in many amorphous or disordered materials [10, 11,
12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49],
athermal systems such as crumpled sheets [50, 51], active matter [52, 53, 54],
cellular protein condensates [55, 56], etc.

The second fundamental paradigm is that of a system brought out of equilibrium
by external 111Active matter systems can also be seen as within this class,
with the driving then arising from e.g. self-propulsion forces acting on the
constituent particles; indeed recent research has revealed some intriguing
connections between external driving by shear and driving by activity [99,
100, 101, 102, 53, 103, 104, 105, 106, 107]. driving such as mechanical stress
that leads to shear flow, or a magnetic field that creates a current around a
ring. For steady driving a non-equilibrium steady state is then often reached,
which balances energy input and dissipation. Of particular interest within
this class of non-equilibrium states have recently been systems driven by
trajectory bias, see e.g. [58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71]. In this approach, the ensemble of equilibrium trajectories is biased by
conditioning on the value of a non-equilibrium observable such as a time-
averaged current. An essentially equivalent (by analogy with equilibrium
statistical mechanics) trajectory ensemble can then be constructed by
application of an appropriate biasing field that is conjugate to the desired
observable. Either way, the bias favours trajectories that in equilibrium
would be rare fluctuations, thus driving the system out of equilibrium. One
attraction of such ensembles, known variously as _biased_ , _tilted_ or
$s$–ensemble [63, 72, 73, 67], is that they are “minimally biased” in the
sense that they maximize the entropy in trajectory space [74] subject to the
desired constraint.

So far, the above two paradigms have largely been considered separately 222One
recent exception comes from studies of the plastic activity under driving
oscillatory strain, which exhibits slow power laws suggestive of aging [108]..
In particular, studies of driving by trajectory bias have focussed on cases
where non-equilibrium steady states are reached, and the dynamical phase
transitions between those as the bias is varied [62, 63, 64, 65, 66]. Here we
bring the two paradigms together and study the driven dynamics of aging
systems. Our results show the existence of two qualitatively different classes
of aging dynamics: fragile aging that is destroyed by any bias, leading either
to an active steady state or an inactive frozen state, and robust aging that
persists even in the presence of bias. As regards our methodology, the
$s$–ensemble can generally be studied by analysing the properties of a biased
master operator [76, 58, 63, 64, 65]; as long as this has a finite spectral
gap, the driving leads to a steady state [64]. To analyse the driven dynamics
of aging systems, we therefore construct an appropriately generalized approach
that is applicable also to systems with a vanishing spectral gap.

We consider in this study trap models (TMs), which are paradigmatic
descriptions of slow, glassy dynamics. They represent such dynamics as a
sequence of jumps between local energy minima or “traps” in configuration
space, and have played a key role in understanding the physical properties of
glasses, see e.g. [77, 12, 78, 79, 80, 81]. What makes them attractive for our
purposes is that – in the absence of driving – they exhibit well understood
aging behaviour [12, 82, 83, 44, 84, 85]. The key signature of this is that
correlation functions become explicitly dependent on two times rather than
just time differences. Equivalently, one can say that relaxation timescales
grow (in principle without bound) with the age of the system, measured from
the time of its preparation in the glassy phase.

Mathematically, TMs are defined as continuous time, discrete state Markov
chains with random energy landscapes [22]. The disorder in the energies (or
trap depths) generates a broad distribution of relaxation times characteristic
of systems with glassy dynamics [9, 80]. This translates into a vanishing
spectral gap (in the thermodynamic limit) of the master operator, and in fact
the spectral density of the slowest modes generically exhibits a power-law
singularity below the glass transition temperature [86, 87, 88, 89].

Figure 1: Dynamical $(s,T)$ phase diagrams for (a) Bouchaud and (b)
Barrat–Mézard trap models; (c) dynamical free energy and (d) mean
time–averaged activity as a function of the bias parameter $s$ for
Barrat–Mézard model, for different temperatures $T$ as shown. In (a) and (b)
the solid line on the $T$-axis above $T_{\rm{g}}$ indicates a first–order
phase transition. Below $T_{\rm{g}}$ the transition is continuous and for BM
it is shifted away from $s=0$ to negative $s$ for $T<1/2$, indicating a robust
aging phase. The circles in (c,d) correspond to the circle on the dynamical
phase transition line in (b).

In this work, we investigate the effects of driving on the aging dynamics of
the two most widely studied mean–field TMs as defined by Bouchaud [12] and
Barrat and Mézard (BM) [90]. The driving is chosen to bias trajectories
towards high or low activity, defined as the number of jumps from one trap to
another. A pictorial representation of our main result is given in Fig. 1a &
b, where we show the phase diagrams in the plane of bias parameter $s$ and
temperature $T$. These diagrams indicate that, as analysed in detail below,
the BM trap model has a robust aging phase at low $T$ while aging in the
Bouchaud model is always fragile, thus revealing two qualitatively distinct
patterns in the response of aging to driving.

We outline briefly the commonalities and differences between the TMs we
consider, before introducing the driving by activity bias. Both models models
assume a disordered landscape of traps, with depths $E>0$ drawn from an
exponential distribution $\rho(E)={\rm{e}}^{-E}$ whose width sets our units of
energy. This leads to the emergence of a glass transition at temperature
$T_{\rm{g}}=1$; above this, there is an well–defined equilibrium phase with an
associated Boltzmann distribution. Below $T_{\rm{g}}$ there is no steady–state
solution for the master equation governing the dynamics (eq. (3) below),
resulting in slow, glassy relaxation and aging. The two models differ in the
choice of transition rates between configurations (or traps): for the Bouchaud
TM they are of Arrhenius form while the BM model posits Glauber rates:

$\displaystyle W_{\rm{Bouchaud}}(E\to E^{\prime})$
$\displaystyle={\rm{e}}^{-E/T}$ (1) $\displaystyle W_{\rm{BM}}(E\to
E^{\prime})$ $\displaystyle=1/[1+{\rm{e}}^{-(E^{\prime}-E)/T}]$ (2)

Assuming mean field connectivity where any trap can be reached from any other,
the (undriven) dynamics of a TM in the thermodynamic limit of a large number
traps is then described by a master equation for the probability $P(E,\tau)$
of being in a trap of depth $E$ at time $\tau$:

$\displaystyle\partial_{\tau}P(E,\tau)$ $\displaystyle=-r(E)P(E,\tau)+$
$\displaystyle\rho(E)\int_{0}^{\infty}dE^{\prime}\,W(E^{\prime}\to
E)P(E^{\prime},\tau)$ (3)

with $r(E)=\int dE^{\prime}\,W(E\to E^{\prime})\rho(E^{\prime})$ the escape
rate from a trap of depth $E$. This is simply ${\rm{e}}^{-E/T}$ for Bouchaud,
and $\propto{\rm{e}}^{-E}$ for sufficiently deep traps in BM [91].

To implement the driving we bias trajectories of the above two TMs according
to their activity [64, 65, 67] $\mathcal{A}$, defined as the total number of
jumps between traps along a trajectory of duration $t$; we write
$\bar{a}=\mathcal{A}/t$ for the time-averaged activity. Introducing a
conjugate bias field $s$, the driven trajectory ensemble is then defined by
reweighting each trajectory with
${\rm{e}}^{-s\mathcal{A}}={\rm{e}}^{-st\bar{a}}$, with positive $s$ favouring
smaller $\mathcal{A}$. The appropriate normalization constant is the dynamical
partition function $Z(s,t)=\langle{\rm{e}}^{-st\bar{a}}\rangle_{0}$, the
average being over the trajectory distribution of the undriven dyamics
($s=0$). Associated with this is a dynamical free energy

$\displaystyle\psi(s)=\lim_{t\to\infty}\frac{1}{t}\ln\langle{\rm{e}}^{-st\bar{a}}\rangle_{0}$
(4)

that can alternatively be seen as the scaled cumulant generating function of
$\bar{a}$. By direct analogy with equilibrium thermodynamics, derivatives
w.r.t. the field $s$ give averages over the biased ensemble:

$\displaystyle-\psi^{\prime}(s)=-\frac{1}{t}\frac{\partial}{\partial s}\ln
Z(s,t)=\frac{\langle\bar{a}\,{\rm{e}}^{-st\bar{a}}\rangle_{0}}{Z(s,t)}=\langle\bar{a}\rangle_{s}$
(5)

where we have left the limit $t\to\infty$ implicit and the last equality
defines the mean time-averaged activity in the driven system. Higher-order
cumulants can be obtained similarly, as well as – via a Legendre transform and
under appropriate assumptions – the rate function for the large deviations of
$\bar{a}$ [65, 67].

The dynamical free energy $\psi(s)$ can be calculated as the largest
eigenvalue of a biased master operator $\mathbb{W}(s)$ [76, 58]. The latter
retains the diagonal elements of the original master operator, i.e. the escape
rates $r(E)$, but changes the off–diagonal elements, i.e. the transition
rates, into

$\displaystyle W(E\to E^{\prime},s)=W(E\to E^{\prime}){\rm{e}}^{-s}$ (6)

This comes from the fact that the weight of any trajectory in the driven
ensemble is multiplied by ${\rm{e}}^{-s}$ for each jump. For the Bouchaud TM,
the resulting eigenvalue problem for $\psi(s)$ can be reduced to a single
self-consistent equation [91]. Determining the mean time-averaged activity
from (5) we find that $\langle\bar{a}\rangle_{s}$ is positive for $s<0$ at all
$T$, while it vanishes for $s>0$ (where also $\psi(s)=0$). In the first regime
we therefore have an active steady state while for $s>0$ the system is
inactive, see Fig. 1a. In these broad features the phase diagram resembles
previous results for kinetically constrained models, with a dynamical phase
transition at $s=0$. There, however, the transition is always first order [62,
63]. For the Bouchaud model, this scenario applies only above the glass
transition, $T>1$, while for lower temperatures we find that the transition is
second order. This in fact makes sense: at $s=0$ for $T<1$ the system ages and
its mean time-averaged activity $\langle\bar{a}\rangle_{0}$ vanishes for
$t\to\infty$. The driven activities $\langle\bar{a}\rangle_{s}$ for $s<0$ then
extrapolate to this zero value as $s\to 0$, making the transition continuous;
we find explicitly the scaling $\langle\bar{a}\rangle_{s}\sim|s|^{\beta-1}$
(with $\beta=1/T$ as usual) [91]. The point $(s=0,T=1)$, i.e. the glass
transition of the undriven system, thus becomes a tricritical point in the
$(s,T)$-phase diagram, at which the order of the transition changes from first
to second (we refer to Refs. [69, 70] for other recent examples of tricritical
behavior in the context of dynamical phase transitions.)

For the BM model we find for $T>1/2$ a qualitatively similar dynamical phase
diagram [91], with a transition from an active to an inactive phase at $s=0$
that is first order for $T>1$ and second order for $1/2<T<1$. This agreement
between the two models for $T>1/2$ is physically reasonable because in this
regime the slow dynamics in the BM model can be shown to arise from activation
across effective energy barriers [92, 93, 94], mirroring the explicitly
activated dynamics (1) in the Bouchaud TM. For lower temperatures, $T<1/2$,
the BM dynamical phase diagram becomes qualitatively different, with the phase
boundary moving towards negative $s$ (Fig. 1b). Using an appropriate scaling
assumption for the eigenvector [91] we can determine its location in closed
form as

$\displaystyle s^{*}=\ln(\sin(\pi T))$ (7)

The absence of a transition at $s=0$ already hints at the fact that low
temperature aging in the BM model is robust to driving, a point we will
substantiate below. The shift of the phase boundary away from $s=0$ also
reflects a change in character of (undriven) BM aging: for $T<1/2$ the slow
dynamics is governed by entropic barriers, i.e. by the scarcity of ever deeper
traps [90, 92, 93].

While the above free energy calculations can distinguish active and inactive
dynamical phases, we need more precise tools to study whether the driven
systems still show aging. In undriven TMs, the aging regime exhibits a
vanishing time-averaged activity, i.e.
$\langle\bar{a}\rangle=\langle\mathcal{A}\rangle/t\to 0$ for $t\to\infty$, but
as the system keeps evolving its total activity $\langle\mathcal{A}\rangle$
still diverges (sublinearly) with trajectory length $t$. We therefore now ask
whether this remains the case in the inactive ($\langle\bar{a}\rangle\to 0$)
regions of our phase diagrams. This would be straightforward for the undriven
dynamics: one writes $\mathcal{A}=\int_{0}^{t}d\tau\,a_{\tau}$ where
$a_{\tau}$ is the time-dependent activity or rate of jumps. Using that the
mean of the latter is the mean escape rate

$\langle a_{\tau}\rangle=\int dE\,P(E,\tau)r(E)$ (8)

one then only needs to solve the Master equation (3) to determine $\langle
a_{\tau}\rangle$. This reasoning does not work, however, for the auxiliary
master operator (6) as this does not represent a physical, probability-
conserving dynamics. Consequently, we develop a formalism to represent the
driven dynamics as an equivalent “auxiliary” Markov process, extending
existing work [64, 65, 72, 67] to aging systems. The appropriate auxiliary
master operator (see also remarks for diffusion processes in [72, 73]) has
transition rates [91]

$\displaystyle W^{\rm{aux}}(E\to E^{\prime},\tau)=W(E\to
E^{\prime})\,{\rm{e}}^{-s}\,\frac{q(E^{\prime},\tau)}{q(E,\tau)}$ (9)

with $q(E,\tau)$ the solution of the backward master equation for
$\mathbb{W}(s)$,

$\displaystyle-{\partial_{\tau}}q(E,\tau)$ $\displaystyle=-r(E)q(E,\tau)$ (10)
$\displaystyle+{\rm{e}}^{-s}\int dE^{\prime}q(E^{\prime},\tau)W(E\to
E^{\prime})\rho(E^{\prime})$

and final time boundary condition $q(E,t)=1$. The corresponding auxiliary
escape rates are given by

$\displaystyle r^{\rm aux}_{\tau}(E)=r(E)-\partial_{\tau}\ln(q(E,\tau))$ (11)

Within this probability-conserving Markov dynamics we can now compute the
desired time-dependent activity as the mean escape rate:

$\displaystyle\langle a_{\tau}\rangle_{s}$ $\displaystyle=\int
dE\,P_{s}(E,\tau)[r(E)-\partial_{\tau}\ln(q(E,\tau)]$ (12)

It then remains to find the distribution of trap depths in the driven systems,
$P_{s}(E,\tau)$, that enters into this average. This has the form

$\displaystyle P_{s}(E,\tau)=\frac{q(E,\tau)P_{\rm{f}}(E,\tau)}{Z(s,t)}$ (13)

where $P_{\rm{f}}(E,\tau)$ is the solution of the forward master equation for
$\mathbb{W}(s)$. In the simple case where $\mathbb{W}(s)$ has a spectral gap
[64], $q(E,\tau)\propto{\rm{e}}^{\psi(s)(t-\tau)}$ while $P_{s}(E,\tau)$
reaches a steady state for large $\tau$ 333Both statements apply away from the
ends of a long trajectory, i.e. once $\tau$ and $t-\tau$ are significantly
larger than the inverse spectral gap. and one retrieves the known result
$r^{\rm aux}(E)=r(E)+\psi(s)$ [59, 60] and a time-independent nonzero activity

$\displaystyle\langle a_{\tau}\rangle_{s}\to\langle\bar{a}\rangle=\langle
r\rangle_{s}+\psi(s)$ (14)

This approach of the activity to a steady state plateau can be seen for the
most negative $s$-values in Fig. 2, consistent with an active phase. In the
inactive phase, we observe in the BM model with $T<1/2$ (Fig. 2b) that the
activity $\langle a_{\tau}\rangle_{s}$ decays like in the undriven case $s=0$,
as $\sim 1/\tau$ 444This can be confirmed analytically at $T=0$, for any bias
$s$ [91].: the system continues to age in the presence of driving, its total
activity $\langle\mathcal{A}\rangle_{s}$ growing as $\ln t$. This is what we
mean by robust aging, and it is confirmed by Fig. 3b: the distribution of trap
depths continues to shift (logarithmically) to larger values with increasing
$\tau$, just as in undriven aging.

The effect of driving on aging in the BM model is rather different for
$T>1/2$: while for $s=0$ one has $\langle
a_{\tau}\rangle_{0}\sim\tau^{1-\beta}$ [92] and so a divergent total activity
$\langle\mathcal{A}\rangle_{0}\sim t^{2-\beta}$ demonstrating aging, any $s>0$
leads to a faster activity decay $\langle
a_{\tau}\rangle_{s}\sim\tau^{-(\alpha+1)}$ with $\alpha>0$ and thus a finite
total activity $\langle\mathcal{A}\rangle_{s}$. This means the inactive phase
is frozen here, with the system coming to rest after a finite number of jumps.
Correspondingly, $P_{s}(E,\tau)$ in Fig. 3a shows a rather rapid transfer of
probability with increasing $\tau$ into traps that are sufficiently deep for
the system not to make further jumps in the remainder of the trajectory.
Overall, the aging in the undriven system ($s=0$) is fragile: any positive
drive freezes the system, while any negative $s$ produces an active steady
state. The Bouchaud model behaves similarly throughout its glass phase, with
an aging activity decay $\langle a_{\tau}\rangle_{0}\sim\tau^{T-1}$ in the
undriven case $s=0$ [12] that changes to $\langle
a_{\tau}\rangle_{s}\sim\tau^{-(T+1)}$ [91] and hence a frozen, finite total
activity phase, for any $s>0$.

Figure 2: Mean time–dependent activity for two temperatures and different
driving parameters $s$ as shown, for BM trap model. Dashed lines: power law
scaling for the undriven case ($s=0$).  Figure 3: Time evolution of trap depth
distribution for BM trap model with fixed driving parameter $s=0.25$ at
temperatures (a) above and (b) below $1/2$. Note in (b) the continual shift of
the distribution to larger trap depths indicating robust aging.

Summarizing, we have shown that the interplay of aging and driving by
trajectory bias leads to significant new insights into aging dynamics. In the
trap models we studied, which in the absence of driving exhibit aging below
their glass transition temperature, this aging dynamics can respond in two
fundamentally different ways to driving. In the BM model for $T>1/2$ and in
the Bouchaud model, aging is fragile: any drive towards high activity ($s<0$)
produces an active steady state, while driving towards low activity leads to
frozen dynamics arresting after a finite number of jumps. For $T<1/2$ in the
BM model, on the other hand, aging is robust, remaining qualitatively
unchanged by moderate driving towards either higher or lower activity. We
demonstrated these differences using a general analysis framework for
Markovian systems driven by trajectory bias, which is complete in the sense
that it can be deployed whether or not the underlying undriven system exhibits
aging. What is fascinating in particular is that the response to driving
detects, without additional physical input via e.g. specifically designed
coarse-graining [94, 97, 85], the qualitative differences between aging
dynamics governed by energetic barriers – which we found to be fragile – and
by entropic effects. The resulting robust aging behavior that we found in the
low–temperature BM model may have broader implications for the physics of
supercooled liquids, bearing in mind a recent proposal that entropic
relaxation is dominant in this regime [98].

More generally, the response of aging systems to driving that we have studied
here, by bringing together the two key non-equilibrium paradigms, provides new
analysis tools for non-equilibrium dynamics and could eventually lead to a
classification of aging into distinct universality classes.

Acknowledgements: FAL acknowledges financial support through a scholarship
from Conacyt (Mexico).

## References

  * Kob and Barrat [1997] W. Kob and J. L. Barrat, Aging effects in a Lennard-Jones glass, Phys. Rev. Lett. 78, 4581 (1997).
  * Leheny and Nagel [1998] R. L. Leheny and S. R. Nagel, Frequency-domain study of physical aging in a simple liquid, Phys. Rev. B 57, 5154 (1998).
  * Kob  _et al._ [2000] W. Kob, J.-L. Barrat, F. Sciortino, and P. Tartaglia, Aging in a simple glass former, J. Phys. Condens. Matter 12, 6385 (2000).
  * Sciortino and Tartaglia [2001] F. Sciortino and P. Tartaglia, Extension of the fluctuation-dissipation theorem to the physical aging of a model glass-forming liquid, Phys. Rev. Lett. 86, 107 (2001).
  * Martelozzo  _et al._ [2002] V. C. Martelozzo, A. B. Schofield, W. C. K. Poon, and P. N. Pusey, Structural aging of crystals of hard-sphere colloids, Phys. Rev. E 66, 021408 (2002).
  * Viasnoff  _et al._ [2003] V. Viasnoff, S. Jurine, and F. Lequeux, How are colloidal suspensions that age rejuvenated by strain application?, Faraday Discuss. 123, 253 (2003).
  * Martinez  _et al._ [2010] V. A. Martinez, G. Bryant, and W. Van Megen, Aging dynamics of colloidal hard sphere glasses, J. Chem. Phys. 133, 114906 (2010).
  * Warren and Rottler [2013] M. Warren and J. Rottler, Quench, Equilibration, and Subaging in Structural Glasses, Phys. Rev. Lett. 110, 10.1103/PhysRevLett.110.025501 (2013).
  * Arceri  _et al._ [2021] F. Arceri, F. P. Landes, L. Berthier, and G. Biroli, A Statistical Mechanics Perspective on Glasses and Aging, in  _Encyclopedia of Complexity and Systems Science_, edited by R. A. Meyers (Springer Berlin Heidelberg, Berlin, Heidelberg, 2021) pp. 1–68.
  * Struik [1978] L. C. E. Struik,  _Physical Aging in Amorphous Polymers and Other Materials_ (Elsevier, Houston, 1978).
  * Sibani and Hoffmann [1989] P. Sibani and K. H. Hoffmann, Hierarchical models for aging and relaxation of spin-glasses, Phys. Rev. Lett. 63, 2853 (1989).
  * Bouchaud [1992] J. P. Bouchaud, Weak ergodicity breaking and aging in disordered-systems, J Phys Fr. 2, 1705 (1992).
  * Cugliandolo and Kurchan [1993] L. F. Cugliandolo and J. Kurchan, Analytical solution of the off-equilibrium dynamics of a long- range spin-glass model, Phys. Rev. Lett. 71, 173 (1993).
  * Baldassarri  _et al._ [1995] A. Baldassarri, L. F. Cugliandolo, J. Kurchan, and G. Parisi, On the out-of-equilibrium order-parameter in long-range spin- glasses, J. Phys. A 28, 1831 (1995).
  * Barrat and Mézard [1995a] A. Barrat and M. Mézard, Phase-space diffusion and low-temperature aging, J Phys Fr. 5, 941 (1995a).
  * Bouchaud and Dean [1995] J. P. Bouchaud and D. S. Dean, Aging on Parisi’s tree, J Phys Fr. 5, 265 (1995).
  * Franz and Hertz [1995] S. Franz and J. Hertz, Glassy transition and aging in a model without disorder, Phys. Rev. Lett. 74, 2114 (1995).
  * Hodge [1995] I. M. Hodge, Physical aging in polymer glasses, Science 267, 1945 (1995).
  * Ritort [1995] F. Ritort, Glassiness in a model without energy barriers, Phys. Rev. Lett. 75, 1190 (1995).
  * Bouchaud  _et al._ [1996] J. P. Bouchaud, L. Cugliandolo, J. Kurchan, and M. Mézard, Mode-coupling approximations, glass theory and disordered- systems, Phys. A 226, 243 (1996).
  * Cugliandolo  _et al._ [1996] L. F. Cugliandolo, J. Kurchan, and P. Le Doussal, Large time out-of-equilibrium dynamics of a manifold in a random potential, Phys. Rev. Lett. 76, 2390 (1996).
  * Monthus and Bouchaud [1996] C. Monthus and J.-P. Bouchaud, Models of traps and glass phenomenology, J. Phys. Math. Gen. 29, 3847 (1996).
  * Rieger [1996] H. Rieger, Aging in disordered systems, Phys. A 224, 267 (1996).
  * Cugliandolo  _et al._ [1997] L. F. Cugliandolo, D. S. Dean, and J. Kurchan, Fluctuation-dissipation theorems and entropy production in relaxational systems, Phys. Rev. Lett. 79, 2168 (1997).
  * Vincent  _et al._ [1997] E. Vincent, J. Hammann, M. Ocio, J.-P. Bouchaud, and L. F. Cugliandolo, Slow dynamics and aging in spin-glasses, in  _Complex Behaviour of Glassy Systems: Proceedings of the XIV Sitges Conference Sitges, Barcelona, Spain_, Vol. 492 (Springer Berlin Heidelberg, 1997) pp. 184–219, arxiv:cond-mat/9607224 .
  * Bouchaud  _et al._ [1998] J. P. Bouchaud, L. F. Cugliandolo, J. Kurchan, and M. Mézard, Out of equilibrium dynamics in spin-glasses and other glassy systems, in  _Spin Glas. Random Fields_ , edited by A. P. Young (World Scientific, Singapore, 1998) pp. 161–223.
  * Berthier  _et al._ [2000] L. Berthier, J. L. Barrat, and J. Kurchan, A two-time-scale, two-temperature scenario for nonlinear rheology, Phys. Rev. E 61, 5464 (2000).
  * Cloître  _et al._ [2000] M. Cloître, R. Borrega, and L. Leibler, Rheological aging and rejuvenation in microgel pastes, Phys. Rev. Lett. 85, 4819 (2000).
  * Fielding  _et al._ [2000] S. M. Fielding, P. Sollich, and M. E. Cates, Aging and rheology in soft materials, J. Rheol. 44, 323 (2000).
  * Rinn  _et al._ [2000] B. Rinn, P. Maass, and J. P. Bouchaud, Multiple scaling regimes in simple aging models, Phys. Rev. Lett. 84, 5403 (2000).
  * Kurchan [2001] J. Kurchan, Rheology, and how to stop aging, in  _Jamming Rheol. Constrained Dyn. Microsc. Macrosc. Scales_ (Taylor and Francis, London and New York, 2001) pp. 72–79.
  * Ramos and Cipelletti [2001] L. Ramos and L. Cipelletti, Ultraslow dynamics and stress relaxation in the aging of a soft glassy system, Phys. Rev. Lett. 87, 245503 (2001).
  * Ben Arous  _et al._ [2002] G. Ben Arous, A. Bovier, and V. Gayrard, Aging in the random energy model, Phys. Rev. Lett. 88, 087201 (2002).
  * Cipelletti and Ramos [2002] L. Cipelletti and L. Ramos, Slow dynamics in glasses, gels and foams, Curr. Opin. Colloid Interface Sci. 7, 228 (2002).
  * Castillo  _et al._ [2002] H. E. Castillo, C. Chamon, L. F. Cugliandolo, and M. P. Kennett, Heterogeneous aging in spin glasses, Phys. Rev. Lett. 88, 237201 (2002).
  * Barkai and Cheng [2003] E. Barkai and Y. C. Cheng, Aging continuous time random walks, J. Chem. Phys. 118, 6167 (2003).
  * Ben Arous  _et al._ [2003] G. Ben Arous, A. Bovier, and V. Gayrard, Glauber dynamics of the random energy model II. Aging below the critical temperature, Comm. Math. Phys. 236, 1 (2003).
  * Fielding and Sollich [2003] S. M. Fielding and P. Sollich, Equivalence of driven and aging fluctuation-dissipation relations in the trap model, Phys. Rev. E 67, 011101 (2003).
  * Montanari and Ricci-Tersenghi [2003] A. Montanari and F. Ricci-Tersenghi, Microscopic description of aging dynamics: Fluctuation- dissipation relations, effective temperature, and heterogeneities, Phys. Rev. Lett. 90, 017203 (2003).
  * Cohen-Addad  _et al._ [2004] S. Cohen-Addad, R. Höhler, and Y. Khidas, Origin of the slow linear viscoelastic response of aqueous foams, Phys. Rev. Lett. 93, 028302 (2004).
  * Ritort [2004] F. Ritort, Spontaneous relaxation in generalized oscillator models with glassy dynamics, J. Phys. Chem. B 108, 6893 (2004).
  * Mayer  _et al._ [2006] P. Mayer, S. Leonard, L. Berthier, J. P. Garrahan, and P. Sollich, Activated aging dynamics and negative fluctuation-dissipation ratios, Phys. Rev. Lett. 96, 030602 (2006).
  * Fielding  _et al._ [2009] S. M. Fielding, M. E. Cates, and P. Sollich, Shear banding, aging and noise dynamics in soft glassy materials, Soft Matter 5, 2378 (2009).
  * Berthier and Biroli [2011] L. Berthier and G. Biroli, Theoretical perspective on the glass transition and amorphous materials, Rev. Mod. Phys. 83, 587 (2011).
  * Sollich  _et al._ [2017] P. Sollich, J. Olivier, and D. Bresch, Aging and linear response in the Hébraud–Lequeux model for amorphous rheology, J. Phys. A: Math. Theor. 50, 165002 (2017).
  * Joshi and Petekidis [2018] Y. M. Joshi and G. Petekidis, Yield stress fluids and ageing, Rheol. Acta 57, 521 (2018).
  * Lubchenko and Wolynes [2018] V. Lubchenko and P. G. Wolynes, Aging, Jamming, and the Limits of Stability of Amorphous Solids, J. Phys. Chem. B 122, 3280 (2018).
  * Parley  _et al._ [2020] J. T. Parley, S. M. Fielding, and P. Sollich, Aging in a mean field elastoplastic model of amorphous solids, Phys. Fluids 32, 127104 (2020).
  * Parley  _et al._ [2022a] J. T. Parley, R. Mandal, and P. Sollich, Mean-field description of aging linear response in athermal amorphous solids, Phys. Rev. Materials 6, 065601 (2022a).
  * Lahini  _et al._ [2017] Y. Lahini, O. Gottesman, A. Amir, and S. M. Rubinstein, Nonmonotonic Aging and Memory Retention in Disordered Mechanical Systems, Phys. Rev. Lett. 118, 085501 (2017).
  * Shohat  _et al._ [2023] D. Shohat, Y. Friedman, and Y. Lahini, Logarithmic aging via instability cascades in disordered systems, Nat. Phys. 10.1038/s41567-023-02220-2 (2023).
  * Janssen  _et al._ [2017] L. M. C. Janssen, A. Kaiser, and H. Löwen, Aging and rejuvenation of active matter under topological constraints, Sci. Rep. 7, 10.1038/s41598-017-05569-6 (2017).
  * Mandal and Sollich [2020] R. Mandal and P. Sollich, Multiple Types of Aging in Active Glasses, Phys. Rev. Lett. 125, 218001 (2020).
  * Janzen and Janssen [2021] G. Janzen and L. M. C. Janssen, Aging in thermal active glasses, ArXiv210505705 Cond-Mat (2021), arxiv:2105.05705 [cond-mat] .
  * Jawerth  _et al._ [2020] L. Jawerth, E. Fischer-Friedrich, S. Saha, J. Wang, T. Franzmann, X. Zhang, J. Sachweh, M. Ruer, M. Ijavi, S. Saha, J. Mahamid, A. A. Hyman, and F. Jülicher, Protein condensates as aging Maxwell fluids, Science 370, 1317 (2020).
  * Takaki  _et al._ [2023] R. Takaki, L. Jawerth, M. Popović, and F. Jülicher, Theory of Rheology and Aging of Protein Condensates, PRX Life 1, 013006 (2023).
  * Note [1] Active matter systems can also be seen as within this class, with the driving then arising from e.g. self-propulsion forces acting on the constituent particles; indeed recent research has revealed some intriguing connections between external driving by shear and driving by activity [99, 100, 101, 102, 53, 103, 104, 105, 106, 107].
  * Lebowitz and Spohn [1999] J. L. Lebowitz and H. Spohn, A gallavotti–cohen-type symmetry in the large deviation functional for stochastic dynamics, Journal of Statistical Physics 95, 333 (1999).
  * Evans [2004a] R. M. L. Evans, Rules for transition rates in nonequilibrium steady states, Phys. Rev. Lett. 92, 150601 (2004a).
  * Evans [2004b] R. M. L. Evans, Detailed balance has a counterpart in non-equilibrium steady states, Journal of Physics A: Mathematical and General 38, 293 (2004b).
  * Lecomte  _et al._ [2007] V. Lecomte, C. Appert-Rolland, and F. Van Wijland, Thermodynamic formalism for systems with markov dynamics, Journal of statistical physics 127, 51 (2007).
  * Garrahan  _et al._ [2007] J. P. Garrahan, R. L. Jack, V. Lecomte, E. Pitard, K. van Duijvendijk, and F. van Wijland, Dynamical first-order phase transition in kinetically constrained models of glasses, Physical review letters 98, 195702 (2007).
  * Garrahan  _et al._ [2009] J. P. Garrahan, R. L. Jack, V. Lecomte, E. Pitard, K. van Duijvendijk, and F. van Wijland, First-order dynamical phase transition in models of glasses: an approach based on ensembles of histories, Journal of Physics A: Mathematical and Theoretical 42, 075007 (2009).
  * Jack and Sollich [2010] R. L. Jack and P. Sollich, Large deviations and ensembles of trajectories in stochastic models, Progress of Theoretical Physics Supplement 184, 304 (2010).
  * Chetrite and Touchette [2013] R. Chetrite and H. Touchette, Nonequilibrium microcanonical and canonical ensembles and their equivalence, Physical review letters 111, 120601 (2013).
  * Jack and Sollich [2015] R. L. Jack and P. Sollich, Effective interactions and large deviations in stochastic processes, The European Physical Journal Special Topics 224, 2351 (2015).
  * Jack [2020] R. L. Jack, Ergodicity and large deviations in physical systems with stochastic dynamics, The European Physical Journal B 93, 1 (2020).
  * Fodor  _et al._ [2022] É. Fodor, R. L. Jack, and M. E. Cates, Irreversibility and biased ensembles in active matter: Insights from stochastic thermodynamics, Annual Review of Condensed Matter Physics 13, 215 (2022).
  * Agranov  _et al._ [2022] T. Agranov, M. E. Cates, and R. L. Jack, Entropy production and its large deviations in an active lattice gas, Journal of Statistical Mechanics: Theory and Experiment 2022, 123201 (2022).
  * Agranov  _et al._ [2023] T. Agranov, M. E. Cates, and R. L. Jack, Tricritical behavior in dynamical phase transitions, Physical Review Letters 131, 017102 (2023).
  * Carugno  _et al._ [2023] G. Carugno, P. Vivo, and F. Coghi, Delocalization-localization dynamical phase transition of random walks on graphs, Phys. Rev. E 107, 024126 (2023).
  * Chetrite and Touchette [2015a] R. Chetrite and H. Touchette, Nonequilibrium markov processes conditioned on large deviations, in  _Annales Henri Poincaré_ , Vol. 16 (Springer, 2015) pp. 2005–2057.
  * Chetrite and Touchette [2015b] R. Chetrite and H. Touchette, Variational and optimal control representations of conditioned and driven processes, Journal of Statistical Mechanics: Theory and Experiment 2015, P12001 (2015b).
  * Pressé  _et al._ [2013] S. Pressé, K. Ghosh, J. Lee, and K. A. Dill, Principles of maximum entropy and maximum caliber in statistical physics, Rev. Mod. Phys. 85, 1115 (2013).
  * Note [2] One recent exception comes from studies of the plastic activity under driving oscillatory strain, which exhibits slow power laws suggestive of aging [108].
  * Derrida and Lebowitz [1998] B. Derrida and J. L. Lebowitz, Exact large deviation function in the asymmetric exclusion process, Physical review letters 80, 209 (1998).
  * Dyre [1987] J. C. Dyre, Master-equation approach to the glass transition, Phys. Rev. Lett. 58, 792 (1987).
  * Denny  _et al._ [2003] R. A. Denny, D. R. Reichman, and J.-P. Bouchaud, Trap models and slow dynamics in supercooled liquids, Physical review letters 90, 025503 (2003).
  * Heuer [2008] A. Heuer, Exploring the potential energy landscape of glass-forming systems: from inherent structures via metabasins to macroscopic transport, Journal of Physics: Condensed Matter 20, 373101 (2008).
  * Scalliet  _et al._ [2021] C. Scalliet, B. Guiselin, and L. Berthier, Excess wings and asymmetric relaxation spectra in a facilitated trap model, The Journal of Chemical Physics 155 (2021).
  * Ridout  _et al._ [2023] S. A. Ridout, I. Tah, and A. J. Liu, Building a” trap model” of glassy dynamics from a local structural predictor of rearrangements, arXiv preprint arXiv:2308.07913 (2023).
  * Godreche and Luck [2001] C. Godreche and J. Luck, Statistics of the occupation time of renewal processes, Journal of Statistical Physics 104, 489 (2001).
  * Arous  _et al._ [2002] G. B. Arous, A. Bovier, and V. Gayrard, Aging in the random energy model, Physical review letters 88, 087201 (2002).
  * Bertin [2013] E. Bertin, Ageing in the trap model as a relaxation further away from equilibrium, Journal of Physics A: Mathematical and Theoretical 46, 095004 (2013).
  * Baity-Jesi  _et al._ [2018] M. Baity-Jesi, G. Biroli, and C. Cammarota, Activated aging dynamics and effective trap model description in the random energy model, Journal of Statistical Mechanics: Theory and Experiment 2018, 013301 (2018).
  * Margiotta  _et al._ [2018] R. G. Margiotta, R. Kühn, and P. Sollich, Spectral properties of the trap model on sparse networks, Journal of Physics A: Mathematical and Theoretical (2018).
  * Margiotta  _et al._ [2019] R. G. Margiotta, R. Kühn, and P. Sollich, Glassy dynamics on networks: local spectra and return probabilities, Journal of Statistical Mechanics: Theory and Experiment 2019, 093304 (2019).
  * Tapias  _et al._ [2020] D. Tapias, E. Paprotzki, and P. Sollich, From entropic to energetic barriers in glassy dynamics: the barrat–mézard trap model on sparse networks, Journal of Statistical Mechanics: Theory and Experiment 2020, 093302 (2020).
  * Tapias and Sollich [2022] D. Tapias and P. Sollich, Localization properties of the sparse barrat-mézard trap model, Physical Review E 105, 054109 (2022).
  * Barrat and Mézard [1995b] A. Barrat and M. Mézard, Phase space diffusion and low temperature aging, Journal de Physique I 5, 941 (1995b).
  * [91] See supplemental material for the development of the formalism and its application to trap models (including numerical details), which also includes Refs. [109-111].
  * Bertin [2003] E. M. Bertin, Cross-over from entropic to thermal dynamics in glassy models, Journal of physics A: mathematical and general 36, 10683 (2003).
  * Sollich [2006] P. Sollich, Trap models with slowly decorrelating observables, Journal of Physics A: Mathematical and General 39, 2573 (2006).
  * Cammarota and Marinari [2015] C. Cammarota and E. Marinari, Spontaneous energy-barrier formation in entropy-driven glassy dynamics, Physical Review E 92, 010301 (2015).
  * Note [3] Both statements apply away from the ends of a long trajectory, i.e. once $\tau$ and $t-\tau$ are significantly larger than the inverse spectral gap.
  * Note [4] This can be confirmed analytically at $T=0$, for any bias $s$ [91].
  * Cammarota and Marinari [2018] C. Cammarota and E. Marinari, Numerical evidences of universal trap-like aging dynamics, Journal of Statistical Mechanics: Theory and Experiment 2018, 043303 (2018).
  * Baity-Jesi  _et al._ [2021] M. Baity-Jesi, G. Biroli, and D. R. Reichman, Revisiting the concept of activation in supercooled liquids, The European Physical Journal E 44, 1 (2021).
  * Henkes  _et al._ [2011] S. Henkes, Y. Fily, and M. C. Marchetti, Active jamming: Self-propelled soft particles at high density, Phys. Rev. E 84, 10.1103/PhysRevE.84.040301 (2011).
  * Nandi and Gov [2018] S. K. Nandi and N. S. Gov, Effective temperature of active fluids and sheared soft glassy materials, Eur. Phys. J. E 41, 117 (2018).
  * Liao and Xu [2018] Q. Liao and N. Xu, Criticality of the zero-temperature jamming transition probed by self-propelled particles, Soft Matter 14, 853 (2018).
  * Mo  _et al._ [2020] R. Mo, Q. Liao, and N. Xu, Rheological similarities between dense self-propelled and sheared particulate systems, Soft Matter 16, 3642 (2020).
  * Mandal and Sollich [2021a] R. Mandal and P. Sollich, Shear-induced orientational ordering in an active glass former, Proc. Natl. Acad. Sci. U. S. A. 118, e2101964118 (2021a).
  * Mandal and Sollich [2021b] R. Mandal and P. Sollich, How to study a persistent active glassy system, J. Phys.: Condens. Matter 33, 184001 (2021b).
  * Morse  _et al._ [2021] P. K. Morse, S. Roy, E. Agoritsas, E. Stanifer, E. I. Corwin, and M. L. Manning, A direct link between active matter and sheared granular systems, Proc. Natl. Acad. Sci. U. S. A. 118, e2019909118 (2021).
  * Villarroel and Düring [2021] C. Villarroel and G. Düring, Critical yielding rheology: From externally deformed glasses to active systems, Soft Matter 17, 9944 (2021).
  * Keta  _et al._ [2023] Y.-E. Keta, R. Mandal, P. Sollich, R. L. Jack, and L. Berthier, Intermittent relaxation and avalanches in extremely persistent active matter, Soft Matter 19, 3871 (2023).
  * Parley  _et al._ [2022b] J. T. Parley, S. Sastry, and P. Sollich, Mean-Field Theory of Yielding under Oscillatory Shear, Phys. Rev. Lett. 128, 198001 (2022b).
  * Note [5] In writing equation (S74), we have assumed that $N$ is still large. A more precise calculation that does not neglect terms of $O(1/N)$ shows that the exact eigenvalue equation for a finite number of configurations is $\frac{1}{N}\sum@\slimits@_{i=1}^{N}\frac{\mathrm{e}^{-E_{i}/T}}{\psi(s)+c\mathrm{e}^{-E_{i}/T}}=\mathrm{e}^{s}$ with $c=(N-1+\mathrm{e}^{-s})/N$. As can be seen, this correction does not affect the conclusions for large $N$ as $c\to 1$.
  * De Bacco  _et al._ [2016] C. De Bacco, A. Guggiola, R. Kühn, and P. Paga, Rare events statistics of random walks on networks: localisation and other dynamical phase transitions, Journal of Physics A: Mathematical and Theoretical 49, 184003 (2016).
  * et al. [2018] V. B. S. et al., Arpack.jl: Julia wrapper for the arpack library designed to solve large-scale eigenvalue problems, https://github.com/JuliaLinearAlgebra/Arpack.jl (2018).

Bringing together two paradigms of non-equilibrium: Driven dynamics of aging
systems
Supplementary Material

##  .1 Biased ensembles of trajectories

In this section we present the formalism of biased ensembles of trajectories
applied to a general Markov chain, with the aim of deriving, and providing
intuition for, the main equations presented in the text. We follow closely the
setup in Ref. [64]. We present the formalism for a discrete configuration
space first. The appropriate modifications for the thermodynamic limit with
its continuous trap depth distributions is given in a later section.

Consider a Markov chain on a set of configurations labelled by $\mathcal{C}$.
The system makes stochastic transitions between configurations with transition
rates $W(\mathcal{C^{\prime}}\to\mathcal{C})$. The probability
$P(\mathcal{C},\tau)$ of finding the system in configuration $\mathcal{C}$ at
time $\tau$ then evolves according to the master equation

$\frac{\partial}{\partial\tau}P(\mathcal{C},\tau)=-r(\mathcal{C})P(\mathcal{C},\tau)+\sum_{\mathcal{C^{\prime}}\neq\mathcal{C}}W(\mathcal{C^{\prime}}\to\mathcal{C})P(\mathcal{C^{\prime}},\tau)$
(S1)

with
$r(\mathcal{C})=\sum_{\mathcal{C^{\prime}}\neq\mathcal{C}}W(\mathcal{C}\to\mathcal{C^{\prime}})$
the escape rate from configuration $\mathcal{C}$. We define the vector
$|P(\tau)\rangle=\sum_{\mathcal{C}}p(\mathcal{C},\tau)|\mathcal{C}\rangle$ (we
use bra-ket notation here; $|\mathcal{C}\rangle$ can simply be viewed as the
Euclidean unit vector in direction $\mathcal{C}$) and introduce the compact
notation $\partial_{\tau}|P(\tau)\rangle=\mathbb{W}|P(\tau)\rangle$ where the
elements of the master operator $\mathbb{W}$ are
$\langle\mathcal{C}|\mathbb{W}|\mathcal{C^{\prime}}\rangle=W(\mathcal{C^{\prime}}\to\mathcal{C})(1-\delta_{\mathcal{C},\mathcal{C^{\prime}}})-\delta_{\mathcal{C},\mathcal{C^{\prime}}}r(\mathcal{C})$.

We consider current-like trajectory observables of the form

$\displaystyle\mathcal{A}=\sum_{{\rm{jumps}}\,\mathcal{C}\to\mathcal{C}^{\prime}}\alpha(\mathcal{C},\mathcal{C}^{\prime})$
(S2)

where $\alpha(\mathcal{C},\mathcal{C}^{\prime})=1$ (for
$\mathcal{C}\neq\mathcal{C}^{\prime}$) gives the activity but other choices
can be used to capture currents of particles or entropy etc. We are interested
in trajectory ensembles that force the average $\langle\mathcal{A}\rangle$ to
be far from its average value $\langle\mathcal{A}\rangle_{0}$. This is most
conveniently done not via a hard constraint but with a bias in the probability
assigned to each trajectory $\pi$:

$P[\pi,s]=\frac{P[\pi,0]\exp\left(-s\mathcal{A}\right)}{Z(s,t)}$ (S3)

After discretizing the total trajectory length $t$ into small steps $\Delta
t$, a trajectory $\pi$ is represented as a sequence of configurations
$\\{\mathcal{C}_{\tau}\\}$ for $\tau=0,\Delta t,2\Delta t,\ldots,t$. In terms
of the original master operator $\mathbb{W}$, the biased probability of a
trajectory can then be written as

$\displaystyle P[\pi,s]$ $\displaystyle=\frac{\left[\prod_{\tau=0}^{t-\Delta
t}\langle\mathcal{C}_{\tau+\Delta t}|(\mathbbold{1}+\Delta
t\mathbb{W})|\mathcal{C}_{\tau}\rangle\mathrm{e}^{-s\alpha(\mathcal{C}_{\tau},\mathcal{C}_{\tau+\Delta
t})}\right]}{Z(s,t)}\times p(\mathcal{C}_{0},0)$ (S4)

with a generic initial distribution
$|P(0)\rangle=\sum_{\mathcal{C}}p(\mathcal{C},0)|\mathcal{C}\rangle$. Here we
use the convention $\alpha(\mathcal{C},\mathcal{C}^{\prime})=0$ for
$\mathcal{C}=\mathcal{C}^{\prime}$, so that the reweighting factor in eq. (S4)
contributes only when the configuration changes, as defined in (S2). Because
of this, the factors in the product can be seen as matrix elements of a biased
master operator defined as (see also [61, 63])

$\displaystyle\langle\mathcal{C}|\mathbb{W}(s)|\mathcal{C^{\prime}}\rangle$
$\displaystyle=\langle\mathcal{C}|\mathbb{W}|\mathcal{C^{\prime}}\rangle\mathrm{e}^{-s\alpha(\mathcal{C^{\prime}},\mathcal{C})}=\left\\{\begin{array}[]{ll}W(\mathcal{C^{\prime}}\to\mathcal{C})\mathrm{e}^{-s\alpha(\mathcal{C}^{\prime},\mathcal{C})},&\mathcal{C}\neq\mathcal{C^{\prime}},\\\
-r(\mathcal{C}),&\mathcal{C}=\mathcal{C^{\prime}}.\end{array}\right.$ (S7)

Summing equation (S4) over all possible trajectories gives the normalization
factor or partition function

$Z(s,t)=\langle e|(\mathbbold{1}+\Delta t\mathbb{W}(s))^{t/\Delta
t}|P(0)\rangle\to\langle e|\mathrm{e}^{\mathbb{W}(s)t}|P(0)\rangle$ (S8)

where $\langle e|=\sum_{\mathcal{C}}\langle\mathcal{C}|$ is a projection
state.

If the original stochastic model has a finite state space and is irreducible,
one sees from the spectral decomposition of the biased master operator

$\displaystyle\mathbb{W}(s)=\sum_{\alpha}|V_{\alpha}\rangle\lambda_{\alpha}\langle
U_{\alpha}|$ (S9)

that the dynamical free energy (as introduced in eq. (4)) has a long–time
limit [58, 61, 64, 67] and is explicitly given by

$\psi(s)=\lim_{t\to\infty}\frac{1}{t}\ln
Z(s,t)=\max_{\alpha}\lambda_{\alpha}.$ (S10)

Thus the dynamical free energy is the largest eigenvalue of $\mathbb{W}(s)$,
which dominates the exponential in (S8) for large $t$. If $\mathbb{W}(s)$ has
a spectral gap between between its largest and second largest eigenvalue, then
the biased dynamics approaches a steady state away from the ends ($\tau$ close
to $0$ or $t$) of a long enough trajectory. The extent of the temporal
boundary regions is given by the inverse spectral gap.

The systems that we study in the manuscript are interesting because they do
not necessarily exhibit a spectral gap. As a matter of fact, in the aging
regime the associated relaxation timescales grow without bound (see Refs. [86,
88] for the spectral distribution of the relaxation rates in Bouchaud and BM
models, respectively). Then even the unbiased master operator $\mathbb{W}$ has
eigenvalues arbitrarily close to the largest eigenvalue of zero, i.e. no gap
in the spectrum. In order to analyse these systems and in particular the time-
dependence along their trajectories, one needs to introduce a time–dependent
_auxiliary master operator_ ; this construction is described below.

##  .2 Construction of time-dependent auxiliary Markov operator

The need for this construction arises from the fact that the biased path
probability (S4) is not in the usual form for stochastic Markov dynamics. This
is because the corresponding biased master operator does not in general
conserve probability, i.e. $\langle e|\mathbb{W}(s)\neq 0$. However a
multiplicative (“Doob”) transform of $\mathbb{W}(s)$ can be used to bring it
back into the standard form, even for systems without a spectral gap that can
exhibit aging. This was shown by Chetrite and Touchette [65, 72] in rather
abstract terms for diffusion processes; we make it concrete here for Markov
jump processes on discrete state spaces.

One starts by introducing multiplicative factors $q(\mathcal{C},\tau)$ into
the path probability (S4):

$\displaystyle
P[\pi,s]=\frac{q(\mathcal{C}_{0},0)p(\mathcal{C}_{0},0)}{Z(s,t)}\left[\prod_{\tau=0}^{t-\Delta
t}q(\mathcal{C}_{\tau+\Delta t},\tau+\Delta t)\langle\mathcal{C}_{\tau+\Delta
t}|(\mathbbold{1}+\Delta t\mathbb{W}(s))|\mathcal{C}_{\tau}\rangle
q^{-1}(\mathcal{C}_{\tau},\tau)\right]$ (S11)

These factors all cancel except for the one at the final time, which therefore
has to be $q(\mathcal{C},t)=1$. The other factors should be chosen so that for
every timestep one has a probability-conserving transition matrix, which
requires

$\sum_{\mathcal{C^{\prime}}}q(\mathcal{C^{\prime}},\tau+\Delta
t)\langle\mathcal{C^{\prime}}|(\mathbbold{1}+\Delta
t\mathbb{W}(s))|\mathcal{C}\rangle q^{-1}(\mathcal{C},\tau)=1$ (S12)

Given the reweighting factors $q(\mathcal{C^{\prime}},\tau+\Delta t)$ at some
time, the factors $q(\mathcal{C},\tau)$ at the previous timestep are therefore
determined. If one collects the $q$-factors in a time-dependent vector
$\langle q(\tau)|=\sum_{\mathcal{C}}\langle\mathcal{C}|q(\mathcal{C},\tau)$,
the relation (S12) reads

$\langle q(\tau)|=\langle q(\tau+\Delta t)|(\mathbbold{1}+\Delta
t\mathbb{W}(s))$ (S13)

or in the continuous time limit $\Delta t\to 0$

$-\partial_{\tau}\langle q(\tau)|=\langle q(\tau)|\mathbb{W}(s)\,,$ (S14)

which is the _backward_ master equation for the biased master operator.
Together with the final-time boundary condition $\langle q(t)|=\langle e|$ one
gets the formal solution

$\langle q(\tau)|=\langle e|\mathrm{e}^{\mathbb{W}(s)(t-\tau)}\,.$ (S15)

Once $\langle q(\tau)|$ has been found, the auxiliary master operator
$\mathbb{W}^{\rm aux}(\tau)$ can be determined. For this, we use that the
contribution from each timestep in (S11) must equal
$\langle\mathcal{C}_{\tau+\Delta t}|(\mathbbold{1}+\Delta t\mathbb{W}^{\rm
aux}(\tau))|\mathcal{C}_{\tau}\rangle$. Dividing both sides of the resulting
equality by $\Delta t$, and simplifying the notation by replacing
$\mathcal{C}_{\tau+\Delta t}\to\mathcal{C^{\prime}}$ and
$\mathcal{C}_{\tau}\to\mathcal{C}$ we get

$\displaystyle\frac{q(\mathcal{C},\tau+\Delta
t)-q(\mathcal{C},\tau)}{q(\mathcal{C},\tau)\Delta
t}\delta_{\mathcal{C},\mathcal{C^{\prime}}}+\frac{q(\mathcal{C^{\prime}},\tau+\Delta
t)}{q(\mathcal{C},\tau)}\langle\mathcal{C^{\prime}}|\mathbb{W}(s)|\mathcal{C}\rangle=\langle\mathcal{C^{\prime}}|\mathbb{W}^{\rm
aux}(\tau)|\mathcal{C}\rangle$ (S16)

For $\Delta t\to 0$ this simplifies (using (S14)) to

$\displaystyle\langle\mathcal{C^{\prime}}|\mathbb{W}^{\rm
aux}(\tau)|\mathcal{C}\rangle$
$\displaystyle=\frac{q(\mathcal{C^{\prime}},\tau)}{q(\mathcal{C},\tau)}\langle\mathcal{C^{\prime}}|\mathbb{W}(s)|\mathcal{C}\rangle-\frac{\delta_{\mathcal{C},\mathcal{C^{\prime}}}}{q(\mathcal{C},\tau)}\sum_{\tilde{\mathcal{C}}}q(\tilde{\mathcal{C}},\tau)\langle\tilde{\mathcal{C}}|\mathbb{W}(s)|\mathcal{C}\rangle$
(S17)

Therefore the time-dependent auxiliary master operator $\mathbb{W}^{\rm
aux}(\tau)$ is obtained by reweighting all transition rates of the unbiased
master operator by an additional multiplicative correction, i.e.

$W^{\rm
aux}(\mathcal{C}\to\mathcal{C^{\prime}},\tau)=W(\mathcal{C}\to\mathcal{C^{\prime}})\mathrm{e}^{-s\alpha(\mathcal{C},\mathcal{C}^{\prime})}\frac{q(\mathcal{C^{\prime}},\tau)}{q(\mathcal{C},\tau)}$
(S18)

for all $\mathcal{C}\neq\mathcal{C^{\prime}}$. Equation (S18) corresponds to
equation (9) in the manuscript.

Before proceeding we comment briefly on our notation for the “backwards
factor” $q$. From equation (S15), it is clear that the dependence of $q$ on
$\tau$ is via the time difference $t-\tau$. We generally leave the
$t$-dependence implicit and write either $\langle q(\tau)|$ or
$q(\mathcal{C},\tau)$ as used above. For concrete calculations below in the
thermodynamic limit, where we will need to perform Laplace transforms with
respect to the variable $t-\tau$, we will use the notation
$q(E,\tau)=Q(E,t-\tau)$ where appropriate.

Having established the auxiliary transition rates in (S18), it now remains to
find the auxiliary escape rates, i.e. the (negative) diagonal terms of the
auxiliary master operator. These follow by combining equations (S14) and (S17)
as @@ I have added minus signs in front of the matrix elements - can you check
please? @@

$r_{\tau}^{\rm aux}(\mathcal{C})=-\langle\mathcal{C}|\mathbb{W}^{\rm
aux}(\tau)|\mathcal{C}\rangle=-\langle\mathcal{C}|\mathbb{W}(s)|\mathcal{C}\rangle-\partial_{\tau}\ln
q(\mathcal{C},\tau)$ (S19)

With the above results in hand, we proceed to show the relation between the
mean auxiliary escape rate and the time-dependent activity (equation (12)).
The argument is in fact the same as for standard Markov jump processes, but we
include it here for completeness. Let us write the mean contribution to
$\mathcal{A}$ from jumps between time $\tau=k\Delta t$ und $\tau+\Delta t$ as

$\displaystyle\langle\alpha(\mathcal{C}_{k},\mathcal{C}_{k+1})\rangle$
$\displaystyle=\sum_{\mathcal{C}_{k+1},\mathcal{C}_{k}}P_{s}(\mathcal{C}_{k},\mathcal{C}_{k+1})\alpha(\mathcal{C}_{k},\mathcal{C}_{k+1})$
(S20)

Here the weighting is with the joint probability
$P_{s}(\mathcal{C}_{k},\mathcal{C}_{k+1})$ in the biased process of the two
configurations at the beginning and end of the time interval,
$\mathcal{C}_{k}$ and $\mathcal{C}_{k+1}$. This probability is
$P_{s}(\mathcal{C}_{k},\mathcal{C}_{k+1})=P_{s}(\mathcal{C}_{k},\tau)\,\Delta
t\,W^{\rm{aux}}(\mathcal{C}_{k}\to\mathcal{C}_{k+1},\tau)$ (provided the two
configurations are distinct; if they are identical then there is no
contribution to (S20) anyway). It is for this step that the probability-
conserving nature of $\mathbb{W}^{\rm aux}(\tau)$ is crucial; had we used
$\mathbb{W}(s)$, the sum over the later states in the trajectory would have
given an additional nontrivial factor, which in fact is proportional to
$q(\mathcal{C}_{k},\tau)$.

Inserting the expression for $P_{s}(\mathcal{C}_{k},\mathcal{C}_{k+1})$, the
mean contribution to our trajectory observable $\mathcal{A}$ between timesteps
$k$ and $k+1$ is

$\displaystyle\langle\alpha(\mathcal{C}_{k},\mathcal{C}_{k+1})\rangle$
$\displaystyle=\sum_{\mathcal{C}_{k+1},\mathcal{C}_{k}}P_{s}(\mathcal{C}_{k},\tau)W^{\rm{aux}}(\mathcal{C}_{k}\to\mathcal{C}_{k+1},\tau)\alpha(\mathcal{C}_{k},\mathcal{C}_{k+1})\,dt$
(S21)

For the activity, where
$\alpha(\mathcal{C},\mathcal{C}^{\prime})=1-\delta_{\mathcal{C},\mathcal{C}^{\prime}}$,
the sum over $\mathcal{C}_{k+1}$ just gives the escape rate. Summing also over
all timesteps up to $K=t/\Delta t$ to get the total activity $\mathcal{A}$ one
has

$\displaystyle\langle\mathcal{A}\rangle_{s}$
$\displaystyle=\sum_{k=0}^{K-1}\sum_{\mathcal{C}_{k}}P_{s}(\mathcal{C}_{k},\tau)r^{\rm{aux}}_{\tau}(\mathcal{C}_{k})dt$
(S22)

In the continuous time limit this reads

$\displaystyle\langle\mathcal{A}\rangle_{s}=\int_{0}^{t}d\tau\,\langle
r^{\rm{aux}}_{\tau}\rangle_{s}\equiv\int_{0}^{t}d\tau\,\langle
a_{\tau}\rangle_{s}$ (S23)

where the last equivalence is just the definition of the time-dependent
activity $a_{\tau}$. Overall, one sees that the mean time-dependent activity
$\langle a_{\tau}\rangle_{s}$ is just $\langle
r^{\rm{aux}}_{\tau}\rangle_{s}$, the average escape rate (calculated from the
auxiliary master operator) in the biased process. This demonstrates equation
(12) in the manuscript.

For concrete calculations based on the above formalism we need an expression
for the probability (again in the driven process) of being in configuration
$\mathcal{C}$ at time $\tau$, $P_{s}(\mathcal{C},\tau)$. If we sum the
trajectory probability (S11) over the configurations at all other times we get

$\displaystyle P_{s}(\mathcal{C},\tau)$ $\displaystyle=\frac{\langle
e|\mathrm{e}^{\mathbb{W}(s)(t-\tau)}|\mathcal{C}\rangle\langle\mathcal{C}|\mathrm{e}^{\mathbb{W}(s)\tau}|P(0)\rangle}{Z(s,t)}$
(S24)

The first factor in the numerator is $q(\mathcal{C},\tau)$; we define the
second one as $P_{\rm f}(\mathcal{C},\tau)$. To obtain eq. (13) in the main
text we just need to show that this solves the forward master equation for
$\mathbb{W}(s)$. But this follows directly because

$P_{\rm f}(\mathcal{C},\tau)=\langle\mathcal{C}|P_{\rm f}(\tau)\ \
\mbox{with}\ \ |P_{\rm
f}(\tau)\rangle=\mathrm{e}^{\mathbb{W}(s)\tau}|P(0)\rangle$ (S25)

and because

${\partial_{\tau}}|P_{\rm f}(\tau)\rangle=\mathbb{W}(s)|P_{\rm
f}(\tau)\rangle$ (S26)

Concerning the physical interpretation, one sees that $P_{\rm
f}(\mathcal{C},\tau)$ contains the effect of the bias propagated forwards from
the past. Note that even though $P_{\rm f}(\mathcal{C},\tau)$ is calculated
from $\mathbb{W}(s)$ and the initial distribution $P(\mathcal{C},0)$ exactly
in the same way as for conventional Markov dynamics, it is not a normalized
probability distribution because $\mathbb{W}(s)$ does not conserve
probability.

The factor $q(\mathcal{C},\tau)$ can now be interpreted similarly as the
effect of the bias acting in the future, propagated backwards to the current
time. Equation (S24) shows that the actual probability of being in
configuration $\mathcal{C}$ is simply a product of these past and future
factors, suitably normalized by $Z(s,t)$. From equations (S8) and (S24) one
also deduces that the normalization factor is nothing but

$\displaystyle Z(s,t)$ $\displaystyle=\sum_{\mathcal{C}}P_{\rm
f}(\mathcal{C},t)$ (S27)

Having set out the general construction for driven aging systems, we describe
briefly how it simplifies in the conventional scenario where the driven
(biased) dynamics reaches a steady state because the biased master operator
has a gap in its spectrum below the leading eigenvalue $\psi(s)$. In such a
scenario, the exponential of $\mathbb{W}(s)$ is dominated by the largest
eigenvalue in (S9) at times longer than the inverse spectral gap. This largest
eigenvalue is $\psi(s)$, so that

$\displaystyle Z(s,t)$ $\displaystyle=\langle
e|V_{0}\rangle\mathrm{e}^{\psi(s)t}\langle U_{0}|P(0)\rangle,$ (S28)
$\displaystyle\langle q(\tau)|$ $\displaystyle=\langle
e|V_{0}\rangle\mathrm{e}^{\psi(s)(t-\tau)}\langle U_{0}|,$ (S29)
$\displaystyle|P_{\rm f}(\tau)\rangle$
$\displaystyle=|V_{0}\rangle\mathrm{e}^{\psi(s)\tau}\langle U_{0}|P(0)\rangle$
(S30)

with $U_{0}$ and $V_{0}$ the top left and right eigenvectors of the biased
master operator $\mathbb{W}(s)$. Bearing in mind the eigenvector normalization
$\langle U_{0}|V_{0}\rangle=1$, one thus retrieves the known [63, 64]
$\tau$-independent steady state distribution $P_{s}(\mathcal{C},\tau)=\langle
U_{0}|\mathcal{C}\rangle\langle\mathcal{C}|V_{0}\rangle$. One also finds that
the escape rate shift $-\partial_{\tau}\ln q(\mathcal{C},\tau)=\psi(s)$ is
just the dynamical free energy,

$\displaystyle r^{\rm aux}(\mathcal{C})=r(\mathcal{C})+\psi(s)$ (S31)

which is in agreement with previous results [60, 59] (see also [64]) and leads
to (14) in the main text.

Summarizing, in this section we have extended existing results for the
representation of driven $s$–ensemble dynamics in terms of a physical
(“auxiliary”) master operator to the case where the $\tau$-dependence of the
latter cannot be ignored, for instance because of the existence of aging
effects. The key identities are (S18) and (S19) for the elements of the
auxiliary master operator, the equations (S14) and (S26) governing the
backward and forward factors, the resulting expression for the
$\tau$-dependent configuration probabilities (S24) and finally, for the
important case of driving by activity bias, the expression (S23) for the total
and time-dependent activities.

We note finally that the formalism applies without changes also if one adds to
the biasing observable a time-integrated contribution of the form

$\Delta t\sum_{k}b(\mathcal{C}_{k})=\int_{0}^{t}d\tau\,b(\mathcal{C}_{\tau})$
(S32)

The biased master operator then merely acquires an additional contribution of
$-s\,b(\mathcal{C})$ to its diagonal element
$\langle\mathcal{C}|\mathbb{W}(s)|\mathcal{C}\rangle$, with a corresponding
additional term $+s\,b(\mathcal{C})$ on the r.h.s. of (S31).

##  .3 Thermodynamic limit and top eigenvector systems

In this subsection we show how to rewrite the master equation (S1) in the
thermodynamic limit of an infinite number of traps identified by their depth
$E$, which in the limit becomes a continuous variable; the result is eq. (3)
in the main text. We also derive the thermodynamic limit version of the
forward (eq. (S26)) and backward (eq. (S14)) master equations. We additionally
write the discrete and continuous version of the top eigenvalue system for the
left and right eigenvectors.

As in the main text we assume that each configuration (or trap) $\mathcal{C}$
has an energy (trap depth) $E_{\mathcal{C}}$. The transition rates are taken
to be functions of the energies only and we normalize them so that the escape
rates in a system of $N$ traps with mean-field connectivity are $O(1)$, hence

$\displaystyle
W(\mathcal{C^{\prime}}\to\mathcal{C})=\frac{1}{N}W(E_{\mathcal{C^{\prime}}}\to
E_{\mathcal{C}})$ (S33)

where the difference between the Bouchand and BM models lies only in the
function $W(E\to E^{\prime})$. We introduce the probability for the system to
be in a trap with depth $E$ at time $\tau$ as
$P(E,\tau)=\sum_{\mathcal{C}}P(\mathcal{C},\tau)\delta(E-E_{\mathcal{C}})$ and
rewrite equation (S1) as an equation for $P(E,t)$. This only requires
multiplying both sides by $\delta(E-E_{\mathcal{C}})$ and summing over all
configurations:

$\displaystyle{\partial_{\tau}}P(E,\tau)$
$\displaystyle=-\frac{1}{N}\sum_{\mathcal{C},\mathcal{C^{\prime}}\neq\mathcal{C}}W(E_{\mathcal{C}}\to
E_{\mathcal{C^{\prime}}})P(\mathcal{C},\tau)\delta(E-E_{\mathcal{C}})+\frac{1}{N}\sum_{\mathcal{C},\mathcal{C^{\prime}}\neq\mathcal{C}}W(E_{\mathcal{C^{\prime}}}\to
E_{\mathcal{C}})P(\mathcal{C^{\prime}},\tau)\delta(E-E_{\mathcal{C}})$ (S34)

The constraint on the sums can safely removed without changing the equality as
the terms with $\mathcal{C}=\mathcal{C^{\prime}}$ cancel. In the transition
rate factor of the last term we can replace $E_{\mathcal{C}}$ by $E$ because
of the Dirac delta; the remaining sum over $\mathcal{C}$ then gives the
density of states
$\rho(E)=\frac{1}{N}\sum_{\mathcal{C}}\delta(E-E_{\mathcal{C}})$:

$\displaystyle\frac{\partial}{\partial\tau}P(E,\tau)$
$\displaystyle=-\frac{1}{N}\sum_{\mathcal{C^{\prime}}}W(E\to
E_{\mathcal{C^{\prime}}})P(E,\tau)+\sum_{\mathcal{C^{\prime}}}W(E_{\mathcal{C^{\prime}}}\to
E)P(\mathcal{C^{\prime}},\tau)\rho(E)$ (S35)

Finally, we introduce a factor of $1=\int
dE^{\prime}\,\delta(E^{\prime}-E_{\mathcal{C^{\prime}}})$ on the r.h.s. of
equation (S35) and perform the sum over $\mathcal{C^{\prime}}$. The result is

$\displaystyle{\partial_{\tau}}P(E,\tau)$
$\displaystyle=-r(E)P(E,\tau)+\rho(E)\int
dE^{\prime}P(E^{\prime},\tau)W(E^{\prime}\to E)$ (S36)

which corresponds to equation (3) in the main text, with the escape rate

$r(E)=\frac{1}{N}\sum_{\mathcal{C^{\prime}}}W(E\to
E_{\mathcal{C^{\prime}}})=\int dE^{\prime}\,W(E\to
E^{\prime})\rho(E^{\prime})$ (S37)

Note that the escape rate for finite $N$ is
$r_{\mathcal{C}}=1/N\sum_{\mathcal{C^{\prime}}\neq\mathcal{C}}W(E_{\mathcal{C}}\to
E_{\mathcal{C^{\prime}}})$. Because of the restriction
$\mathcal{C^{\prime}}\neq\mathcal{C}$ this differs from $r(E)$ by a term of
$O(1/N)$, but this becomes irrelevant in the thermodynamic limit $N\to\infty$.
Such corrections are ignored by default in all derivations that follow.

The calculation for the continuous limit of the forward master equation (eq.
(S26)) is essentially the same as for the unbiased case. One only has to
consider the biased master operator as defined in (S7) instead of the
conventional one. The result is, for $N\to\infty$,

$\displaystyle{\partial_{\tau}}P_{\rm f}(E,\tau)$ $\displaystyle=-r(E)P_{\rm
f}(E,\tau)+{\rm{e}}^{-s}\rho(E)\int dE^{\prime}P_{\rm
f}(E^{\prime},\tau)W(E^{\prime}\to E)$ (S38)

For the backward equation (eq. (S14)) we show the calculation in more detail
because some care is required in the definition of the thermodynamic limit of
the backward factor:

$\displaystyle
q(E,\tau)\rho(E)=\frac{1}{N}\sum_{\mathcal{C}}q(\mathcal{C},\tau)\delta(E-E_{\mathcal{C}})$
(S39)

The factor $\rho(E)$ on the l.h.s. is needed because in the unbiased case the
backward factor should reduce to a unit constant,
$q(E,\tau)=q(\mathcal{C},\tau)=1$, corresponding to the conservation of
probability. From (S15) the same has to hold generally, even in the presence
of a trajectory bias, at the final time $\tau=t$. We now start from equation
(S14), insert the decomposition $\langle
q(\tau)|=\sum_{\mathcal{C}}q(\mathcal{C},\tau)\langle C|$ and multiply by
$|\mathcal{C^{\prime}}\rangle$ from the right to project onto configuration
$\mathcal{C^{\prime}}$:

$\displaystyle-\sum_{\mathcal{C}}{\partial_{\tau}}q(\mathcal{C},\tau)$
$\displaystyle\langle\mathcal{C}|\mathcal{C^{\prime}}\rangle=\sum_{\mathcal{C}}q(\mathcal{C},\tau)\langle\mathcal{C}|\mathbb{W}(s)|\mathcal{C^{\prime}}\rangle$
(S40) $\displaystyle-{\partial_{\tau}}q(\mathcal{C^{\prime}},\tau)$
$\displaystyle=\sum_{\mathcal{C}}q(\mathcal{C},\tau)\big{[}W(\mathcal{C^{\prime}}\to\mathcal{C}){\rm{e}}^{-s}(1-\delta_{\mathcal{C},\mathcal{C^{\prime}}})-r(\mathcal{C^{\prime}})\delta_{\mathcal{C},\mathcal{C^{\prime}}}\big{]}$
(S41)
$\displaystyle=\sum_{\mathcal{C},\mathcal{C}\neq\mathcal{C^{\prime}}}q(\mathcal{C},\tau)W(\mathcal{C^{\prime}}\to\mathcal{C}){\rm{e}}^{-s}-q(\mathcal{C^{\prime}},\tau)\sum_{\mathcal{C},\mathcal{C}\neq\mathcal{C^{\prime}}}W(\mathcal{C^{\prime}}\to\mathcal{C})$
(S42)

where in the last line we have used the definition (for finite $N$) of the
escape rate $r(\mathcal{C^{\prime}})$. Now express the transition rates in
terms of energies (eq. (S33)), multiply both sides of the equation by
$\delta(E-E_{\mathcal{C^{\prime}}})/N$ and sum over all configurations
$\mathcal{C^{\prime}}$. This yields

$\displaystyle-\rho(E){\partial_{\tau}}q(E,\tau)=\frac{{\rm{e}}^{-s}}{N^{2}}\sum_{\mathcal{C}\neq\mathcal{C^{\prime}}}q(\mathcal{C},\tau)W(E_{\mathcal{C^{\prime}}}\to
E_{\mathcal{C}})\delta(E-E_{\mathcal{C^{\prime}}})-\frac{1}{N^{2}}\sum_{\mathcal{C}\neq\mathcal{C^{\prime}}}q(\mathcal{C^{\prime}},\tau)W(E_{\mathcal{C^{\prime}}}\to
E_{\mathcal{C}})\delta(E-E_{\mathcal{C^{\prime}}})$ (S43)

The constraint $\mathcal{C}\neq\mathcal{C^{\prime}}$ only gives $O(1/N)$
corrections so can be dropped; performing the sum over $\mathcal{C^{\prime}}$
then yields

$\displaystyle-\rho(E){\partial_{\tau}}q(E,\tau)={\rm{e}}^{-s}\frac{\rho(E)}{N}\sum_{\mathcal{C}}q(\mathcal{C},\tau)W(E\to
E_{\mathcal{C}})-\frac{1}{N}q(E,\tau)\rho(E)\sum_{\mathcal{C}}W(E\to
E_{\mathcal{C}})$ (S44)

Finally we cancel the common factor of $\rho(E)$ and again switch to an
integration over trap depths on the r.h.s. by inserting $1=\int
dE^{\prime}\delta(E^{\prime}-E_{\mathcal{C}})$ and performing the remaining
sum over $\mathcal{C}$:

$\displaystyle-{\partial_{\tau}}q(E,\tau)$ $\displaystyle={\rm{e}}^{-s}\int
dE^{\prime}q(E^{\prime},\tau)W(E\to E^{\prime})\rho(E^{\prime})-q(E,\tau)\int
dE^{\prime}\,W(E\to E^{\prime})\rho(E^{\prime})$ (S45)

As the final integral is just $r(E)$, this is equation (Bringing together two
paradigms of non-equilibrium: Driven dynamics of aging systems) in the main
text.

In the explicit calculations below we will often work with Laplace transforms.
As $q(E,\tau)$ depends on the trajectory length $t$ and the running time
$\tau$ only via the time difference $\Delta t=t-\tau$ (cf. eq. (S15)), it will
then be helpful to work with the quantity $Q(E,\Delta t)=q(E,\tau)$, for which
equation (S45) becomes

$\displaystyle\frac{\partial}{\partial\Delta t}Q(E,\Delta t)$
$\displaystyle=-r(E)Q(E,\Delta t)+{\rm{e}}^{-s}\int
dE^{\prime}Q(E^{\prime},\Delta t)W(E\to E^{\prime})\rho(E^{\prime})$ (S46)

We finally write down, for later reference, the equations for the top left and
top right eigenvectors of the biased master operator, in the discrete case
(finite $N$) and in the continuous ($N\to\infty$) limit. These eigenvectors
correspond directly to solutions of the backward and forward master equations,
respectively, with the time derivatives replaced by the top eigenvalue, which
according to eq. (S10) is the dynamical free energy $\psi(s)$. For the top
right eigenvector we use the notation $V_{0}(\mathcal{C})$ in the discrete
case. The eigenvalue problem then reads

$\displaystyle-
V_{0}(\mathcal{C})r(\mathcal{C})+{\rm{e}}^{-s}\sum_{\mathcal{C^{\prime}}\neq\mathcal{C}}V_{0}(\mathcal{C^{\prime}})W(\mathcal{C^{\prime}}\to\mathcal{C})=\psi(s)V_{0}(\mathcal{C})$
(S47)

and becomes in the continuous limit (using that
$V_{0}(E)=\sum_{\mathcal{C}}\delta(E-E_{\mathcal{C}})V_{0}(\mathcal{C})$)

$\displaystyle-V_{0}(E)r(E)$ $\displaystyle+{\rm{e}}^{-s}\rho(E)\int
dE^{\prime}W(E^{\prime}\to E)V_{0}(E^{\prime})=\psi(s)V_{0}(E)$ (S48)

Notice that the left hand side is exactly the same as the right hand side of
equation (S38) with $P_{\rm f}$ replaced by $V_{0}$, as it should be.

For the left eigenvector $U_{0}(\mathcal{C})$ we have the following condition
for a discrete set of states

$\displaystyle-
U_{0}(\mathcal{C^{\prime}})r(\mathcal{C^{\prime}})+\sum_{\mathcal{C}\neq\mathcal{C^{\prime}}}U_{0}(\mathcal{C})W(\mathcal{C^{\prime}}\to\mathcal{C}){\rm{e}}^{-s}=\psi(s)U_{0}(\mathcal{C^{\prime}})$
(S49)

The continuous limit reads (with
$U_{0}(E)\rho(E)=1/N\sum_{\mathcal{C}}\delta(E-E_{\mathcal{C}})U_{0}(\mathcal{C})$)

$\displaystyle-U_{0}(E)r(E)$ $\displaystyle+{\rm{e}^{-s}}\int
dE^{\prime}\rho(E^{\prime})W(E\to
E^{\prime})U_{0}(E^{\prime})=\psi(s)U_{0}(E)$ (S50)

with the l.h.s. again matching eq. (S45). Note that because the activity is
invariant under time-reversal and because the undriven dynamics obeys detailed
balance, i.e. $e^{\beta E}W(E\to E^{\prime})=e^{\beta
E^{\prime}}W(E^{\prime}\to E)$, the left and right eigenvector problems can be
related to each other and one finds in particular that $V_{0}(E)\propto
U_{0}(E)\rho(E)e^{\beta E}$.

##  .4 Bouchaud trap model

In this section we provide calculational details for the activity–biased
dynamics in the Bouchaud trap model. We start for completeness with the
calculation of the mean time–averaged activity in the unbiased case, which
will act as a baseline for the results in the presence of driving by
trajectory bias.

We note that as the Bouchaud trap model is a renewal process, some of our
results can be related to those from the classical reference [82]. We will
comment on the specific instances of this connection below. The link is made
via the distribution of interarrival times, i.e. times between jumps, which is
given by

$\displaystyle\rho(\tau)=\int_{1}^{\infty}d\tilde{\tau}\,\frac{1}{\tilde{\tau}}\exp(-\tau/\tilde{\tau})\tilde{\rho}(\tilde{\tau})$
(S51)

in terms of $\tilde{\rho}(\tilde{\tau})=T\tilde{\tau}^{-(T+1)}$, the mean
trapping time distribution (which in turn follows from
$\rho(E)=\mathrm{e}^{-E}$ and $\tilde{\tau}=1/r(E)=\mathrm{e}^{E/T}$).
Evaluating the integral in (S51) gives

$\displaystyle\rho(\tau)=T\tau^{-(T+1)}[T\Gamma(T)-\Gamma(1+T,\tau)]$ (S52)

with $\Gamma(\cdot)$ and $\Gamma(\cdot,\cdot)$ the complete and (upper)
incomplete Gamma function, respectively. The latter can be neglected
asymptotically and a comparison with equation (1.6) from reference [82] then
reveals that the mapping between the notation there and here is $\theta\mapsto
T$, $\tau_{0}\mapsto(T\Gamma(T))^{1/T}$.

###  .4.1 Zero bias

Our starting point is the master equation (3) with the appropriate transition
rates (eq. (1)) and escape rates $r(E)=W(E\to E^{\prime})=\mathrm{e}^{-E/T}$:

$\displaystyle{\partial_{\tau}}P(E,\tau)$
$\displaystyle=-{\rm{e}}^{-E/T}P(E,\tau)+\rho(E)\int_{0}^{\infty}dE^{\prime}\,{\rm{e}}^{-E^{\prime}/T}P(E^{\prime},\tau)$
(S53)

We consider throughout a uniform distribution across traps as initial
condition, corresponding to $P(E,0)=\rho(E)={\rm{e}}^{-E}$. Solving the
dynamical equation in the Laplace domain [22], with $z$ conjugate to $\tau$, a
straightforward calculation yields

$\displaystyle\hat{P}(E,z)$
$\displaystyle=\frac{{\rm{e}}^{-E}}{(z+{\rm{e}}^{-E/T})f(z)}$ (S54)

with (all integrations over $E$ run from 0 to $\infty$ in what follows)

$f(z)=\int
dE\,\frac{z\,{\rm{e}}^{-E}}{z+{\rm{e}}^{-E/T}}={}_{2}F_{1}(1,T,1+T,-1/z)$
(S55)

which as shown can be expressed in terms of the Gauss hypergeometric function
${}_{2}F_{1}$. The mean time–dependent activity, which is equal to the mean
escape rate by eq. (S23), is

$\langle a_{\tau}\rangle_{0}=\int dE\,P(E,\tau)\mathrm{e}^{-E/T}$ (S56)

and has Laplace transform

$\langle\hat{a}(z)\rangle_{0}=\int
dE\,\hat{P}(E,z)\mathrm{e}^{-E/T}=\frac{1}{f(z)}\int
dE\,\frac{{\rm{e}}^{-E-E/T}}{z+{\rm{e}}^{-E/T}}=\frac{1-f(z)}{f(z)}$ (S57)

While the above expressions are pleasingly explicit, the inverse Laplace
transform cannot be obtained in closed form. We therefore focus on the
physically most interesting limit of long times, $z\to 0$. With the variable
change $x=\mathrm{e}^{-E}$ we obtain

$f(z)=\int_{0}^{1}dx\,\frac{z}{z+x^{\beta}}$ (S58)

where $\beta=1/T$ as usual. For temperatures below the glass transition,
$T<1$, the substitution $y=x/z^{T}$ transforms this to

$f(z)=z^{T}\int_{0}^{z^{-T}}\\!\\!\frac{dy}{1+y^{\beta}}$ (S59)

For $z\to 0$ the upper integration boundary diverges and the integral
evaluates to $\pi T/\sin(\pi T)$. Combining with the straightforward expansion
for $T>1$, one has overall for small $z$

$f(z)\approx\left\\{\begin{array}[]{lcl}z/(1-\beta)&\mbox{for}&T>1\\\
z^{T}\frac{\pi T}{\sin(\pi T)}&\mbox{for}&T<1\end{array}\right.$ (S60)

Inserting into (S57) and transforming back to the time domain shows that for
long times the activity $\langle a_{\tau}\rangle_{0}$ approaches the constant
value $1-\beta$ when $T>1$, indicating a conventional equilibrium state. For
$T<1$, on the other hand,

$\displaystyle\langle a_{\tau}\rangle_{0}$ $\displaystyle\approx\frac{\sin(\pi
T)}{\pi T\Gamma(T)}\tau^{T-1}$ (S61)

This prediction compares well against results from numerical solutions of the
master equation as shown in Fig. S1. The mean total activity can be obtained
as the integral of the last expression (cf. eq. (S23)), with the result

$\displaystyle\langle\mathcal{A}\rangle_{0}\approx\frac{\sin(\pi T)}{\pi
T^{2}\Gamma(T)}t^{T}$ (S62)

which corresponds exactly to equation (3.6) in reference [82].

Figure S1: Evolution of the mean time–dependent activity with time for
different biases at $T=0.4$ for the Bouchaud TM. Solid lines are from
numerical solutions, obtained as the time–domain version of equation (S93).
The analytical solution for the unbiased case (eq. (S61)) is shown by the
green dashed line, while the red and blue dashed lines display the predictions
(S95) for positive biases. Total trajectory length (for the biased cases)
$t=10^{6}$.

###  .4.2 Bias towards high activity

Our starting point is the right eigenvalue equation (S48) with the appropriate
transition and escape rates (see eq. (1))

$\displaystyle\psi(s)V_{0}(E)$
$\displaystyle=-{\rm{e}}^{-E/T}V_{0}(E)+{\rm{e}}^{-s}\rho(E)\int
dE^{\prime}{\rm{e}}^{-E^{\prime}/T}V_{0}(E^{\prime})$ (S63)

which may be rearranged to read

$\displaystyle V_{0}(E)$ $\displaystyle=\frac{{\rm{e}}^{-s}\rho(E)\int
dE^{\prime}{\rm{e}}^{-E^{\prime}/T}V_{0}(E^{\prime})}{\psi(s)+{\rm{e}}^{-E/T}}\propto\frac{\rho(E)}{\psi(s)+{\rm{e}}^{-E/T}}$
(S64)

Inserting back into eq. (S63) gives the condition

$\int
dE\,\rho(E)\frac{\mathrm{e}^{-E/T}}{\psi(s)+\mathrm{e}^{-E/T}}=\mathrm{e}^{s}$
(S65)

For the case of an exponential density of states that we consider throughout,
$\rho(E)=\mathrm{e}^{-E}$ for $E>0$, the function on the l.h.s. becomes simply
$1-f(\psi(s))$ with $f$ defined in (S55), giving

$f(\psi(s))=1-\mathrm{e}^{s}$ (S66)

The function $f(z)$ increases monotonically with $z$ from $f(0)=0$ to
$f(z\to\infty)=1$; the condition (S66) therefore has exactly one positive
solution $\psi(s)$ for all $s<0$. The other eigenvalues are all negative as
one can see from the finite-$N$ form of the eigenvalue condition given in
(S74) below. As the discussion there shows, the second eigenvalue approaches
zero for $N\to\infty$ so that the spectrum has a gap, of size $\psi(s)$.

Figure S2: (Left) Dynamical free energy in the Bouchaud trap model against
bias parameter $s$, for two different temperatures. (Right) Mean time–averaged
activity. The kink at $s=0$ shows a first order transition from an active
($s<0$) to an inactive dynamical phase. The free energy is obtained as the
numerical solution of condition (S66). The activity curves are obtained by
first–order perturbation theory as explained in the text.

To summarize, for bias towards high activity ($s<0$) the dynamical free energy
of the Bouchaud model is obtained by solving equation (S65) with an
exponential density of states, or explicitly (S66). The resulting free energy
is plotted as a function of $s$ for two different temperatures in Fig. S2. In
the plot we have anticipated that the dynamical free energy goes to zero for
$s>0$, as we will show in the next subsection. The derivative
$-\psi^{\prime}(s)=\langle\bar{a}\rangle_{s}$ is nonzero for $s<0$ so that we
always obtain a nonzero time–averaged activity: stationarity in the Bouchaud
model is restored by bias to high activity for all temperatures, even those
($T<1$) where the unbiased model is in its glassy phase and exhibits aging.

Differences between the regime below and above the glass transition become
visible, however, if we consider the behaviour of the dynamical free energy
around $s=0$. For small negative $s$ we can expand the exponential in equation
(S65) as ${\rm{e}}^{s}=1-|s|$, giving the condition $|s|=f(\psi(s))$. Using
then the small-$z$ behaviour (S60) of $f(z)$ one finds that around the phase
boundary $s=0$, the dynamical free energy on the active side ($s<0$) behaves
as

$\displaystyle\psi(s)$ $\displaystyle\approx(1-\beta)|s|\,,\qquad T>1$ (S67)
$\displaystyle\psi(s)$ $\displaystyle\approx b_{0}|s|^{\beta}\,,\qquad T<1$
(S68)

with $b_{0}=\left({\sin(\pi T)}/{\pi T}\right)^{1/T}$. Equations (S67) and
(S68) tell us that the dynamical phase transition is first order at high
temperatures, but continuous in the temperature range where the unbiased
system is glassy. Correspondingly, the time–averaged activity
$\langle{a}\rangle_{s}=\psi^{\prime}(s)$ has a jump at the transition for
$T>1$, while it approaches zero continuously for $T<1$ (cf. Fig. S2(right)).

Let us connect the above result with those of reference [82], focussing on the
glassy regime $T<1$. One can exploit the Gärtner–Ellis theorem to compute the
rate function $I(a)$ of the time–averaged activity $\bar{a}=\mathcal{A}/t$ in
the long–time limit, see e.g. [67], as

$\displaystyle I(\bar{a})=\sup_{s}[-s\bar{a}-\psi(s)]$ (S69)

Using now equation (S68) (together with $\psi(s)=0$ for $s>0$, see next
subsection), one finds directly that the scaling of the rate function around
the minimum $a=0$ must be (for $\bar{a}>0$)

$\displaystyle I(\bar{a})\sim\bar{a}^{\frac{1}{1-T}}$ (S70)

Thus to exponential order for large $t$, the distribution of the time–averaged
activity (at zero bias) is

$\displaystyle
P(\bar{a})\asymp\mathrm{e}^{-tI(\bar{a})}\sim\mathrm{e}^{-ct\bar{a}^{1/(1-T)}}$
(S71)

with $\asymp$ meaning an equality up to subexponential prefactors and $c$ some
$T$-dependent constant. In terms of the total activity, eq. (S71) translates
into

$\displaystyle
P(\mathcal{A})\asymp\mathrm{e}^{-ct^{-T/(1-T)}\mathcal{A}^{1/(1-T)}}=\mathrm{e}^{-c(\mathcal{A}/t^{T})^{1/(1-T)}}$
(S72)

This scaling of the activity distribution corresponds exactly to that found
for “large” activities in reference [82], by a saddle-point approximation in
equation (3.8). Interestingly, the scaling regime considered there was that of
activities of order of the average (S62), i.e. $\mathcal{A}\sim t^{T}$, while
in the large deviation rate function $I(\bar{a})$ one considers
$\mathcal{A}\sim t$. What the comparison above shows, then, is that the result
for the probability distribution of $\mathcal{A}$ in that regime (for large
$\mathcal{A}/t^{T}$) agrees with that of our calculation for small $|s|$, i.e.
small $\bar{a}=\mathcal{A}/t$. In other words, the two approaches match for
$t^{T}\ll\mathcal{A}\ll t$ where the relevant scaling regimes meet. Note that
our numerical results for the dynamical free energy for large $|s|$, where the
asymptotic expression (S68) no longer applies, have no analogue in reference
[82] as they relate to larger activity values with $a=\mathcal{A}/t$ of order
one.

Figure S3: Steady-state trap depth distributions (eq. (S73)) in the biased
Bouchaud trap model for different values of $s$. Left: $T=2$, with
corresponding equilibrium distribution shown as a solid line. Right: $T=1/2$,
showing that the distribution develops a maximum around increasingly deep
traps as $|s|$ decreases.

We can gain further insight into the behaviour for bias towards high activity
by studying the steady state distribution of trap depths. As explained after
(S30), the steady state probability for a given configuration is
$P_{s}(\mathcal{C})=\langle
U_{0}|\mathcal{C}\rangle\langle\mathcal{C}|V_{0}\rangle$, and correspondingly
in the thermodynamic limit $P_{s}(E)=U_{0}(E)V_{0}(E)$. We have already
determined the right eigenvector $V_{0}(E)$; the left eigenvector $U_{0}(E)$
can be found similarly starting from equation (S50) or, as discussed below
that equation, one can use directly that $U_{0}(E)\propto
V_{0}(E)/[\rho(E)e^{\beta E}]$. Either way one finds for the steady state
distribution

$P_{s}(E)\propto\frac{\mathrm{e}^{-(1-\beta)E}}{(1+\psi(s)\mathrm{e}^{\beta
E})^{2}}$ (S73)

and some example plots are shown in Fig. S3. While negative bias leads to
active phases for all $T$ as we saw above, the behaviour of $P_{s}(E)$ as
$s\to 0$ shows the signature of the glass transition: for $T>1$, $P_{s}(E)$
smoothly approaches the Boltzmann equilibrium distribution of the unbiased
system; in the glass regime $T<1$, on the other hand, $P_{s}(E)$ develops a
maximum that shifts to arbitrarily large trap depths as $s\to 0$. This makes
sense physically: without bias, the system ages into deeper traps over time.
Any high activity bias stops this aging, but it does so in increasingly deep
regions of the landscape as the bias gets weaker.

###  .4.3 Bias towards low activity

We next look at the effect of biasing towards lower values than typical
activity using a bias parameter $s>0$. That the dynamical free energy in this
regime tends to zero for $N\to\infty$ can be seen easily starting from the
analogue of the condition (S65) for a finite number $N$ of traps, which reads

$\displaystyle\frac{1}{N}\sum_{\mathcal{C}}\frac{\mathrm{e}^{-E_{\mathcal{C}}/T}}{\psi(s)+\mathrm{e}^{-E_{\mathcal{C}}/T}}=\mathrm{e}^{s}$
(S74)

Let us assume that the trap depths are arranged in increasing order:
$0<E_{1}<\ldots<E_{N}$. It is clear that the left hand side of equation (S74)
as a function of $\psi$ has poles at $\\{-\mathrm{e}^{-E_{\mathcal{C}}/T}\\}$,
and therefore the largest eigenvalue $\psi(s)$ must lie between
$-\mathrm{e}^{-\beta E_{N}}$ and 0 555In writing equation (S74), we have
assumed that $N$ is still large. A more precise calculation that does not
neglect terms of $O(1/N)$ shows that the exact eigenvalue equation for a
finite number of configurations is
$\frac{1}{N}\sum_{i=1}^{N}\frac{\mathrm{e}^{-E_{i}/T}}{\psi(s)+c\mathrm{e}^{-E_{i}/T}}=\mathrm{e}^{s}$
with $c=(N-1+\mathrm{e}^{-s})/N$. As can be seen, this correction does not
affect the conclusions for large $N$ as $c\to 1$.. Because $E_{N}$ diverges
for $N\to\infty$, then $\psi(s)\to 0$ in the limit. The implication is that
$s>0$ always produces an inactive dynamical phase, where the time–averaged
activity $\langle{\bar{a}}\rangle_{s}=-\psi^{\prime}(s)$ also goes to zero.
Note that by extension of the above argument, the second largest eigenvalue of
$\mathbb{W}(s)$ lies between $-\mathrm{e}^{-\beta E_{N-1}}$ and
$-\mathrm{e}^{-\beta E_{N}}$. These two bounds also vanish for large $N$ and
so there is no gap spectral gap in the thermodynamic limit. Thus to understand
the dynamics in the inactive biased phase, we need to use the general
formalism introduced in an earlier section, with the time-dependent weight
factors $q(E,\tau)$.

To look at the time-dependent trap depth distribution $P_{s}(E,\tau)$
(equation (S24)) for a given $s$ and total trajectory length $t$ we need three
ingredients: the forward probability $P_{\rm{f}}(E,\tau)$, the backward factor
$q(E,\tau)=Q(E,t-\tau)$ and the normalizing partition function $Z(s,t)$. If
one Laplace transforms these with respect to the relevant time variables
$\tau$, $t-\tau$ and $t$, respectively, explicit expressions can be found. As
before we choose an initial state where the system finds itself in a randomly
chosen trap, corresponding to $P_{\rm{f}}(E,0)=\rho(E)=\mathrm{e}^{-E}$. The
initial condition for the backward factor is, according to the discussion
after equation (S11), $q(E,t)=Q(E,0)=1$.

The forward (eq. (S38)) and backward (eq. (S46)) equations in Laplace space
are

$\displaystyle z\hat{P}_{\rm{f}}(E,z)-P_{\rm{f}}(E,0)$
$\displaystyle=-{\rm{e}}^{-E/T}\hat{P}_{\rm{f}}(E,z)+{\rm{e}}^{-s}\rho(E)\int
dE^{\prime}\,{\rm{e}}^{-E^{\prime}/T}\hat{P}_{\rm{f}}(E^{\prime},z)$ (S75)
$\displaystyle z\hat{Q}(E,z)-Q(E,0)$
$\displaystyle=-{\rm{e}}^{-E/T}\hat{Q}(E,z)+{\rm{e}}^{-s}{\rm{e}}^{-E/T}\int
dE^{\prime}\,\rho(E^{\prime})\hat{Q}(E^{\prime},z)$ (S76)

with corresponding solutions

$\displaystyle\hat{P}_{\rm{f}}(E,z)$
$\displaystyle=\frac{{\rm{e}}^{-E}}{(z+{\rm{e}}^{-E/T})\tilde{f}(z)}$ (S77)
$\displaystyle\hat{Q}(E,z)$
$\displaystyle=\frac{1}{z+{\rm{e}}^{-E/T}}\left(1+\frac{\mathrm{e}^{-s}f(z)}{\tilde{f}(z)}\frac{\mathrm{e}^{-E/T}}{z}\right)$
(S78)

with the abbreviation

$\tilde{f}(z)=1-{\rm{e}}^{-s}(1-f(z)).$ (S79)

In the unbiased case $s=0$ one has $\tilde{f}(z)=f(z)$ and the solution for
$\hat{P}_{\rm{f}}$ agrees with the one for $\hat{P}$ in (S54) as it must,
while similarly $\hat{Q}(E,z)$ simplifies to $1/z$, i.e. $Q(E,\Delta t)=1$ as
expected.

For the case with bias towards low activity of interest here, $s>0$, the key
difference is that $\tilde{f}(z)$ now has a finite limit $1-\mathrm{e}^{-s}$
for $z\to 0$. Using this in (S77) and transforming back to the time domain
then gives directly for long times

$\displaystyle P_{\rm{f}}(E,\tau)$
$\displaystyle\approx\frac{\mathrm{e}^{-r(E)\tau}\mathrm{e}^{-E}}{1-\mathrm{e}^{-s}}$
(S80)

with $r(E)=\mathrm{e}^{-E/T}$ as before. For $\hat{Q}$, if we assume that the
relevant traps are those with escape rate $r(E)$ of order $z$, the second term
in brackets in (S78) can be neglected for $z\to 0$ because $f(z)\to 0$, giving

$\displaystyle q(E,\tau)$ $\displaystyle\approx\mathrm{e}^{-r(E)(t-\tau)}$
(S81)

The partition function is found from equation (S27), which in the
thermodynamic limit reads

$\displaystyle Z(s,t)=\int_{0}^{\infty}dE\,P_{\rm{f}}(E,t)$ (S82)

The Laplace transform w.r.t. total trajectory length $t$ of this is, from
(S77),

$\hat{Z}(s,z)=\frac{f(z)}{z\tilde{f}(z)}$ (S83)

Using (S60) one then sees that for long times

$\displaystyle Z(s,t)$
$\displaystyle\approx\frac{T\Gamma(T)}{1-\mathrm{e}^{-s}}t^{-T}$ (S84)

(This result applies also for $T>1$, where the singular term from (S60) for
$T<1$ is no longer the leading small $z$ contribution but still gives the
dominant large $t$ behaviour.) The complete distribution according to equation
(13) is then

$\displaystyle P_{s}(E,\tau)$
$\displaystyle\approx\frac{t^{T}}{T\Gamma(T)}\mathrm{e}^{-E}\mathrm{e}^{-\mathrm{e}^{-E/T}t}$
(S85)

Something remarkable happens here: the forward factor $P_{\rm{f}}(E,\tau)$ and
the backward factor $q(E,\tau)$ both depend on the time $\tau$ along the
trajectory, but this dependence cancels from the overall trap depth
distribution. The physical interpretation is that after a (large enough,
$\tau\gg 1$) transient time the trap depth distribution no longer changes for
the remainder of the trajectory. Intriguingly, this distribution is
independent of the specific value of the bias $s$, as long as it is positive:
the $s$ dependence of $P_{\rm f}(E,\tau)$ and $Z(s,t)$ cancels from
$P_{s}(E,\tau)$.

Further intuition can be obtained by noting that $P_{s}(E,\tau)$ in (S85) is
peaked at values of $E$ around $T\ln t$ because it can be written as a
function of $\delta E=E-T\ln t$ alone. For values of $E$ around the peak one
has $r(E)t=\mathrm{e}^{-E/T}t\sim 1$, meaning that the system will only make a
number of jumps of $O(1)$ within the entire duration of the trajectory. The
implication is that low activity-trajectories are realized by the system
rapidly descending into traps that are so deep that it remains essentially
stuck for the rest of the time, the dynamics then being “frozen”.

The character of the low activity-dynamics can also be seen from numerical
solutions for the time evolution of $P_{s}(E,\tau)$ as shown in Fig. S4. Even
at $\tau=1$ this distribution has a contribution proportional to the limiting
form (S85), because the bias prefers trajectories that already start in deep
traps. As $\tau$ grows, the weight of the peak at $E$ of order unity is
shifted into the large $\tau$-peak as trajectories starting from shallow traps
rapidly descend.

Figure S4: Trap depth distribution $P_{s}(E,\tau)$ in the Bouchaud trap model,
at temperature $T=0.6$, bias parameter $s=0.25$, total trajectory length
$t=10^{4}$, evaluated at different running times $\tau$. Also shown is the
prediction from equation (S85).

We can substantiate the above picture by looking in more detail at the
activity along biased trajectories. Considering first the mean long
time–averaged activity, one finds from equations (S84) and (5) that

$\langle\bar{a}\rangle_{s}=-\frac{1}{t}\partial_{s}\ln
Z(s,t)=\frac{1}{t}\frac{1}{\mathrm{e}^{s}-1}$ (S86)

Correspondingly the total activity is

$\langle\mathcal{A}\rangle_{s}=t\langle\bar{a}\rangle_{s}=\frac{1}{\mathrm{e}^{s}-1}$
(S87)

consistent with our expectation that the system will only make a finite number
of jumps along the entire trajectory. As an aside we note that this number
diverges as $s\to 0$ as it must, to connect with the extensive ($\sim t$)
number of jumps that we found in the active phase ($s<0$, see for instance
Fig. S2).

To get more detailed information we next consider the mean time–dependent
activity. This can be evaluated using equation (12). Using the approximations
(S81) and (S85), for which we assumed deep traps (with $e^{-E/T}\sim z$) to be
dominant, would actually yield a vanishing result. This tells us that we need
a more accurate treatment that includes the contribution from the shallow
traps. To achieve this, we first rewrite (12) as

$\langle a_{\tau}\rangle_{s}=\int dE\,\left[r(E)+\frac{\partial_{\Delta
t}Q(E,\Delta t)}{Q(E,\Delta t)}\right]\frac{Q(E,\Delta
t)P_{\rm{f}}(E,\tau)}{Z(s,t)}=\int dE\,\frac{S(E,\Delta
t)P_{\rm{f}}(E,\tau)}{Z(s,t)}$ (S88)

where

$S(E,\Delta t)=r(E)Q(E,\Delta t)+\partial_{\Delta t}Q(E,\Delta t)$ (S89)

Now from (S78), the Laplace transform of this function is

$\hat{S}(E,z)=\left(1+\frac{\mathrm{e}^{-s}f(z)}{\tilde{f}(z)}\frac{\mathrm{e}^{-E/T}}{z}\right)-Q(E,\Delta
t=0)=\mathrm{e}^{-s}\mathrm{e}^{-E/T}\hat{\sigma}(z)$ (S90)

with

$\hat{\sigma}(z)=\frac{f(z)}{z\tilde{f}(z)}$ (S91)

Carrying out the integration over $E$ in (S88) we then obtain for the mean
time-dependent activity

$\langle a_{\tau}\rangle_{s}=\frac{\sigma(\Delta t)\pi(\tau)}{Z(s,t)}$ (S92)

where

$\hat{\pi}(z)=\mathrm{e}^{-s}\int
dE\,\mathrm{e}^{-E/T}\hat{P}_{\rm{f}}(E,z)=\mathrm{e}^{-s}\,\frac{1-f(z)}{\tilde{f}(z)}$
(S93)

and $\sigma(\Delta t)$ and $\pi(\tau)$ are the inverse Laplace transforms of
$\hat{\sigma}(z)$ and $\hat{\pi}(z)$, respectively. For $\sigma(\Delta t)$, we
see by comparing (S83) and (S91) that in fact $\sigma(\Delta t)=Z(s,\Delta t)$
and so from (S84) $\sigma(\Delta t)/Z(s,t)=(\Delta t/t)^{-T}$, which is unity
if we are not near the end of the trajectory, i.e. as long as $\tau\ll t$ so
that $\Delta t\approx t$. Moving on to $\pi(\tau)$, for long times $f(z)$ is
small and we can expand

$\hat{\pi}(z)=\mathrm{e}^{-s}\frac{1-f(z)}{1-\mathrm{e}^{-s}+\mathrm{e}^{-s}f(z)}\approx\frac{1}{\mathrm{e}^{s}-1}-f(z)\,\frac{\mathrm{e}^{s}}{(\mathrm{e}^{s}-1)^{2}}$
(S94)

The singular term from (S60) then determines the long-time behaviour (for
$1\ll\tau\ll t$), giving for the mean time–dependent activity the final result

$\langle
a_{\tau}\rangle_{s}=\pi(\tau)=\frac{T^{2}\Gamma(T)\mathrm{e}^{s}}{(\mathrm{e}^{s}-1)^{2}}\tau^{-T-1}$
(S95)

The power law exponent here is $-(T+1)$ as stated in the main text. The
prediction (S95) is tested against numerics in Fig. S1, showing good agreement
including the prefactor.

##  .5 Barrat–Mézard trap model

In this section, we consider activity–biased dynamics in the Barrat–Mézard
trap model. Because of the more complicated transition rates (eq. (2))
obtaining analytical results is more difficult than for the Bouchaud case.
Here, we focus on the regime that is qualitatively different from the Bouchaud
trap model, namely temperatures $T<1/2$ (cf. the phase diagrams in Fig. 1 in
the main text).

We start as in the previous section with the calculation of the mean
time–dependent activity in the unbiased case. Then, we set out the arguments
that allow us to obtain the phase boundary of equation (7). Finally, we
calculate the mean time–dependent activity in the presence of a bias, at zero
temperature, and discuss the validity of the resulting scaling across the
whole robust aging phase (as shown e.g. in Fig. 2(b) in the main text).

###  .5.1 Zero bias

The master equation (3) with the appropriate transition rates (eq. (2)) and
corresponding escape rate $r(E)$ reads

$\displaystyle\frac{\partial}{\partial\tau}P(E,\tau)$
$\displaystyle=-r(E)P(E,\tau)+\rho(E)\int
dE^{\prime}\frac{P(E^{\prime},\tau)}{1+{\rm{e}}^{-(E-E^{\prime})/T}}$ (S96)

with initial condition $P(E,0)=\rho(E)={\rm{e}}^{-E}$. The escape rate is

$r(E)=\int dE^{\prime}\,W(E\to E^{\prime})\rho(E^{\prime})=\int
dE^{\prime}\frac{\mathrm{e}^{-E^{\prime}}}{1+{\rm{e}}^{-(E^{\prime}-E)/T}}=\int
dE^{\prime}\frac{\mathrm{e}^{-E^{\prime}}\mathrm{e}^{-E/T}}{\mathrm{e}^{-E/T}+{\rm{e}}^{-E^{\prime}/T}}=f(\mathrm{e}^{-E/T})$
(S97)

with $f(z)$ as defined in (S55). In the long time limit, typical trap depths
are large and we can focus on the leading behaviour (S60) of $f(z)$ for small
arguments, giving

$\displaystyle r(E)\approx\frac{\pi T}{\sin(\pi T)}{\rm{e}}^{-E}$ (S98)

This is proportional to the fraction $\mathrm{e}^{-E}$ of traps of depth
greater than that of the initial trap, which is consistent with the slow
dynamics being governed by entropic barriers. The inverse of $r(E)$ yields the
corresponding entropic timescale [93] and it is then natural, following Ref.
[92], to introduce the scaling variable $u^{-1}=\tau{\rm{e}}^{-E}$. The
behavior for large and small $u$ thus corresponds to deep and shallow traps,
respectively. It turns out that in terms of $u$ an asymptotic solution of
equation (S96) can be found (see details in [92, 93]). When converted back to
the physical variables this solution reads

$\displaystyle P(E,\tau)$ $\displaystyle\sim\tau{\rm{e}}^{-E}\,,\quad$
$\displaystyle\mathrm{e}^{E}/\tau\gg 1$ (S99) $\displaystyle P(E,\tau)$
$\displaystyle\sim\tau^{1-1/T}\mathrm{e}^{-E}{\rm{e}}^{E/T}\,,\quad$
$\displaystyle\mathrm{e}^{E}/\tau\ll 1$ (S100)

We comment briefly on the regime $1/2<T<1$ where activation effects dominate
[92, 94]. One then finds, in analogy with the situation in the Bouchaud model,
that the mean time–dependent activity is dominated by the shallow traps, i.e.
by the part of the distribution corresponding to equation (S100). With this we
can compute the scaling of the mean time–dependent activity as

$\displaystyle\langle a_{\tau}\rangle_{0}$ $\displaystyle=\int
dE\,P(E,\tau)r(E)\sim\tau^{1-1/T}$ (S101)

This was the result used in Fig. 2(a).

For the regime $T<1/2$ that we mostly focus on in this section one finds, on
the other hand, that the mean time–dependent activity is governed by typical
traps with $\mathrm{e}^{E}\sim\tau$ [92, 93]. This results in $\langle
a_{\tau}\rangle_{0}\sim r(E)\sim 1/\tau$, with an exponent that no longer
depends on temperature $T$.

###  .5.2 Dynamical phase boundary for $T<1/2$

Here we give the arguments to derive equation (7), which describes the phase
boundary of the continuous dynamical phase transition at $T<1/2$ (cf. phase
diagram in Fig. 1(b)). We tackle the problem by approaching the boundary from
the active side, i.e. $s<s^{*}$.

Our starting point is the left eigenvector equation (S50), with the transition
rates (2). We write it as follows

$\displaystyle-r(E)U_{0}(E)+\mathrm{e}^{-s}\int
dE^{\prime}\frac{\rho(E^{\prime})U_{0}(E^{\prime})}{1+{\rm{e}}^{-(E^{\prime}-E)/T}}=\psi(s)U_{0}(E)$
(S102)

with $r(E)\approx\frac{\pi T}{\sin(\pi T)}{\rm{e}}^{-E}$ (eq. (S98)) and
$\rho(E)=\mathrm{e}^{-E}$. With the change of variables $x={\rm{e}}^{-E}$
equation (S102) becomes

$\displaystyle\mathrm{e}^{-s}\int_{0}^{1}dx^{\prime}\frac{U_{0}(x^{\prime})}{1+(x^{\prime}/x)^{\beta}}=U_{0}(x)(\psi(s)+c_{T}x)$
(S103)

with $c_{T}=\frac{\pi T}{\sin(\pi T)}$ and $\beta=1/T$. The form of the r.h.s.
of (S103) suggests a solution with the scaling $U_{0}(x)\propto g(x/\psi(s))$
with $g$ a master curve to be determined; numerical results supporting this
scaling form are shown in Fig. S5. Inserting the assumed scaling and changing
variables once more to $y=x/\psi(s)$ we get

$\displaystyle\mathrm{e}^{-s}\int_{0}^{1/\psi(s)}dy^{\prime}\frac{g(y^{\prime})}{1+(y^{\prime}/y)^{\beta}}=g(y)(1+c_{T}y)$
(S104)

and at the phase boundary the upper integration boundary can be taken to
infinity. Figure S5 suggests that the resulting $g(y)$ exhibits a crossover
between two different power laws at $y=1$. The exponent we find numerically
for small $y$ is consistent with the value $\beta$. Indeed, evaluating the two
sides of (S104) to leading order for small $y$ gives

$\displaystyle
y^{\beta}\mathrm{e}^{-s}\int_{0}^{\infty}dy^{\prime}\frac{g(y^{\prime})}{y^{\prime}{}^{\beta}}=g(y)$
(S105)

and hence $g(y)\sim y^{\beta}$. Of course the remaining integral on the l.h.s.
has to be finite for this conclusion to hold: we find below that it is.

For large $y$ the exponent is nontrivial and we next establish the equations
required to determine it. We rewrite eq. (S104) for $y\gg 1$ as

$\displaystyle\mathrm{e}^{-s}\int_{0}^{\infty}dy^{\prime}\frac{g(y^{\prime})}{1+(y^{\prime}/y)^{\beta}}\approx
g(y)c_{T}y$ (S106)

and insert the power law scaling $g(y)\sim y^{\alpha}$, with the exponent
$\alpha$ to be determined, leading to

$\displaystyle
c_{T}\mathrm{e}^{s}=\int_{0}^{\infty}dy^{\prime}\frac{y^{\prime\alpha}}{1+y^{\prime\beta}}$
(S107)

This first condition relates $\alpha$ and $s$, but does not yet determine the
position $s^{*}$ of the phase boundary.

Figure S5: Normalized left (top) eigenvector components $U_{0}(E)$ as a
function of the reduced variable $y=x/\psi(s)=\mathrm{e}^{-E}/\psi(s)$
obtained from the numerical solution of equation (S102) for different values
of $s$ close to the phase boundary $s^{*}$. The resulting dynamical free
energies are shown in the legend. Temperature: $T=0.2$. Collapse of the curves
into a single (master) curve indicates correctness of the assumption
$U_{0}(x)\propto g(x/\psi(s))$. Dashed line: $g(y)\sim y^{\beta}$. Dotted
line: $g(y)\sim y^{\beta/2-1}$.

For the second equation we use a physical argument, based on an analogy with
the unbiased case [94]. There, the dynamics changes at $T=1/2$ from being
dominated by activation across effective energy barriers to slowing down by
entropic barriers. This change is related to the fact that for $T>1/2$, the
system is more likely to move towards shallower traps in each jump, and so
eventually to the top of the landscape; while for $T<1/2$ it moves on average
towards deeper traps. At $T=1/2$, the probabilities of moving to shallower and
to deeper traps,

$P_{\uparrow}(E)=\frac{1}{r(E)}\int_{0}^{E}dE^{\prime}\,W(E\to
E^{\prime})\rho(E^{\prime})\quad\mbox{and}\quad
P_{\downarrow}(E)=\frac{1}{r(E)}\int_{E}^{\infty}dE^{\prime}\,W(E\to
E^{\prime})\rho(E^{\prime})$ (S108)

are exactly equal (for deep traps, $E\gg 1$) as one can check by explicit
calculation using $\rho(E^{\prime})=\mathrm{e}^{-E^{\prime}}$. Based on the
picture that at the dynamical phase boundary for $T<1/2$ the driving is just
strong enough to stop the system from aging, i.e. to allow it to return to
shallow traps, we now posit that the analogous criterion holds there with the
appropriate weighting by the left eigenvector $U_{0}(E^{\prime})$, so that

$P_{\uparrow}(E)\propto\int_{0}^{E}dE^{\prime}\,W(E\to
E^{\prime})\rho(E^{\prime})U_{0}(E^{\prime})\quad\mbox{and}\quad
P_{\downarrow}(E)\propto\int_{E}^{\infty}dE^{\prime}\,W(E\to
E^{\prime})\rho(E^{\prime})U_{0}(E^{\prime})$ (S109)

should again be equal. In terms of the scaling function $g(y)$ this equality
reads, by analogy with the l.h.s. of (S104),

$\displaystyle\int_{y}^{\infty}dy^{\prime}\frac{g(y^{\prime})}{1+(y^{\prime}/y)^{\beta}}=\int_{0}^{y}dy^{\prime}\frac{g(y^{\prime})}{1+(y^{\prime}/y)^{\beta}}$
(S110)

or in terms of $w=y^{\prime}/y$,

$\displaystyle\int_{1}^{\infty}dw\,\frac{g(wy)}{1+w^{\beta}}=\int_{0}^{1}dw\,\frac{g(wy)}{1+w^{\beta}}$
(S111)

Finally we take the limit $y\to\infty$ because at the dynamical phase
transition, where $\psi(s)\to 0$, all traps effectively have large
$y=\mathrm{e}^{-E}/\psi(s)$. With the assumed power law scaling
$g(wy)\sim(wy)^{\alpha}$ for large $y$ our condition then becomes

$\displaystyle\int_{1}^{\infty}dw\,\frac{w^{\alpha}}{1+w^{\beta}}=\int_{0}^{1}dw\,\frac{w^{\alpha}}{1+w^{\beta}}$
(S112)

giving as the non–trivial solution for the exponent $\alpha$

$\displaystyle\alpha=\frac{\beta}{2}-1$ (S113)

This prediction for the exponent determining the large-$y$ behaviour of $g(y)$
is fully consistent with our scaling collapse of the numerical eigenproblem
solutions in Fig. S5. It also ensures that the integrand in (S105) scales as
$y^{\prime}{}^{-\beta/2-1}$ for large $y^{\prime}$ and so the integral
converges. Finally, inserting the exponent value (S113) into equation (S107)
and solving for $s$ one finds expression (7) in the manuscript for the
location of the dynamical phase boundary.

###  .5.3 Zero temperature

In the unbiased model, we saw above that the decay of the time-dependent
activity at $T<1/2$ has a $T$-independent power law; in fact the physics of
the slow dynamics here is driven by an entropic mechanism [92, 93] that exists
down to zero temperature. We therefore now focus on this limit and solve for
the scaling of the mean time–averaged activity at $T=0$, allowing for any
value of the bias (including the undriven case $s=0$). We find that the
activity continues to scale as $1/\tau$, consistent with our numerical
observation that this scaling applies for any $T<1/2$ in the robust aging
phase.

At zero temperature, the transition rates (eq. (2)) reduce to Heaviside
functions ($\Theta(x)=0$ for $x<0$ and $\Theta(x)=1$ for $x\geq 0$) and the
escape rate becomes the (complement of the) cumulative distribution of the
density of states, namely
$r(E)=\int_{0}^{\infty}dE^{\prime}\Theta(E^{\prime}-E)\rho(E^{\prime})=\int_{E}^{\infty}dE^{\prime}\rho(E^{\prime})$.
In what follows, we will consider a general distribution $\rho(E)$ because for
zero temperature and zero bias it is known that the results are in fact
independent of $\rho(E)$ [90].

We write for this scenario the forward (eq. (S38)) and backward (eq. (S46))
equations as

$\displaystyle\frac{\partial}{\partial\tau}P_{\rm{f}}(E,\tau)$
$\displaystyle=-r(E)P_{\rm{f}}(E,\tau)+\mathrm{e}^{-s}\rho(E)\int_{0}^{E}dE^{\prime}P_{\rm{f}}(E^{\prime},\tau)$
(S114) $\displaystyle\frac{\partial}{\partial\Delta t}Q(E,\Delta t)$
$\displaystyle=-r(E)Q(E,\Delta
t)+\mathrm{e}^{-s}\int_{E}^{\infty}dE^{\prime}\rho(E^{\prime})Q(E^{\prime},\Delta
t)$ (S115)

with $\Delta t=t-\tau$ as before. It will be convenient to transform those
equations from trap depths $E$ to the corresponding escape rates $r(E)$. We
take into account that $P_{\rm{f}}(E,\tau)$ transforms as a probability
distribution, i.e.
$P_{\rm{f}}(E,\tau)=P_{\rm{f}}(r,\tau)\left|\frac{dr}{dE}\right|=P_{\rm{f}}(r,\tau)\rho(E)$;
whereas $Q(E,t-\tau)$ transforms as a scalar function. Hence the transformed
equations are

$\displaystyle\frac{\partial}{\partial\tau}P_{\rm{f}}(r,\tau)$
$\displaystyle=-rP_{\rm{f}}(r,\tau)+\mathrm{e}^{-s}\int_{r}^{1}dr^{\prime}P_{\rm{f}}(r^{\prime},\tau)$
(S116) $\displaystyle\frac{\partial}{\partial\Delta t}Q(r,\Delta t)$
$\displaystyle=-rQ(r,\Delta
t)+\mathrm{e}^{-s}\int_{0}^{r}dr^{\prime}\,Q(r^{\prime},\Delta t)$ (S117)

Now we Laplace transform the equations, using the initial conditions
$P_{\rm{f}}(r,0)=\frac{P_{\rm{f}}(E,0)}{\rho(E)}=1$ and $Q(r,0)=1$, to yield

$\displaystyle\hat{P}_{\rm{f}}(r,z){(z+r)}$
$\displaystyle=1+\mathrm{e}^{-s}\int_{r}^{1}dr^{\prime}\hat{P}_{\rm{f}}(r^{\prime},z)$
(S118) $\displaystyle\hat{Q}(r,z)(z+r)$
$\displaystyle=1+\mathrm{e}^{-s}\int_{0}^{r}dr^{\prime}\,\hat{Q}(r^{\prime},z)$
(S119)

These equations have the explicit solutions

$\displaystyle\hat{P}_{\rm{f}}(r,z)$
$\displaystyle=\frac{(z+1)^{b}}{(z+r)^{b+1}}$ (S120)
$\displaystyle\hat{Q}(r,z)$ $\displaystyle=\frac{(z+r)^{b-1}}{z^{b}}$ (S121)

with $b={\rm{e}}^{-s}$. The partition function in Laplace space follows as

$\displaystyle\hat{Z}(s,z)=\int_{0}^{1}dr\,\hat{P}_{\rm{f}}(r,z)=\frac{(z+1)^{b}}{bz^{b}}-\frac{1}{b}$
(S122)

We now focus as usual on the long-time limit $z\ll 1$, where we expect the
dynamics to concentrate around deep traps $r\ll 1$. In $\hat{P}_{\rm{f}}(r,z)$
the singularity with the largest real part is at $z=-r$. The numerator
$(z+1)^{b}$ is unity there to leading order so we can approximate
$P_{\rm{f}}(r,z)\approx(z+r)^{-b-1}$. For $\hat{Q}(z)$ the dominant
singularity is at $z=0$ so that to leading order $\hat{Q}(r,z)=r^{b-1}z^{-b}$.
Using the analogous argument for $\hat{Z}(s,z)$ and carrying out the inverse
Laplace transforms we obtain

$\displaystyle P_{\rm{f}}(r,\tau)$
$\displaystyle\approx\frac{{\rm{e}}^{-r\tau}\tau^{b}}{\Gamma(1+b)}$ (S123)
$\displaystyle q(r,\tau)$
$\displaystyle\approx\frac{(t-\tau)^{b-1}r^{b-1}}{\Gamma(b)}$ (S124)
$\displaystyle Z(s,t)$ $\displaystyle\approx\frac{t^{b-1}}{\Gamma(1+b)}$
(S125)

and hence

$\displaystyle
P_{s}(r,\tau)\approx\frac{{\rm{e}}^{-r\tau}r^{b-1}\tau^{b}(t-\tau)^{b-1}}{\Gamma(b)t^{b-1}}$
(S126)

The auxiliary escape rate can now be calculated from eq. (12) and becomes
$r+(b-1)/(t-\tau)$. The correction term is negligibly small as long as we are
far from the end of the trajectory, i.e. $\tau\ll t$. We therefore have

$\displaystyle\langle a_{\tau}\rangle_{s}$
$\displaystyle\approx\int_{0}^{1}dr\,P_{s}(r,\tau)r$ (S127)
$\displaystyle\approx\frac{(t-\tau)^{b-1}b}{\tau t^{b-1}}$ (S128)
$\displaystyle\approx\frac{\mathrm{e}^{-s}}{\tau}$ (S129)

assuming again $\tau\ll t$. In Figure S6 we test this prediction against
numerical solutions of the forward and backward equations, and see that both
the asymptotic scaling and its prefactor are captured correctly. The mean
total activity is from (S125)

$\langle\mathcal{A}\rangle_{s}=-\partial_{s}\ln Z(s,t)=\mathrm{e}^{-s}\ln t$
(S130)

to leading order for large $t$. This is consistent with (S129); it also shows
that the system does not freeze, i.e. that the mean number of jumps diverges
with $t$. The fact that the system continues to age can also be seen from the
explicit form of the time-dependent distribution across traps (again for
$\tau\ll t$)

$P_{s}(r,\tau)\propto
P_{\rm{f}}(r,\tau)q(r,t-\tau)\propto\mathrm{e}^{-r\tau}(r\tau)^{b-1}$ (S131)

This shows simple aging scaling, i.e. a dependence only on $r(E)\tau$, so that
the distribution shifts to deeper and deeper traps as $\tau$ increases.

In the main text, in Fig. 2(b), we show that the scaling from eq. (S129) is
consistent also with numerical results for finite temperatures in the robust
aging phase. This supports our conjecture that the scaling $\langle
a_{\tau}\rangle\sim\tau^{-1}$ is independent of temperature throughout this
phase, with both temperature and bias only affecting the prefactor as for the
case $T=0$.

Figure S6: Mean time–dependent activity (eq. (12)) for $T=0$ obtained via the
numerical integration of equations (S114) and (S115) (solid lines) compared
against the scaling prediction (S129) (dashed lines), for different biases $s$
as shown. Total trajectory length $t=10^{5}$.

##  .6 First order perturbation theory

In this section, we explain how we obtain the mean time–averaged activity as
shown in e.g. Fig. 1(d) (see also Fig. S2(right)) , without having to take
numerical derivatives of the dynamical free energy $\psi(s)$. We work for
simplicity with a finite number $N$ of traps $\mathcal{C}$ and exploit the
fact that $\mathbb{W}(s)$ can be symmetrized by a similarity transformation
involving the diagonal matrix associated with the formal equilibrium
(Boltzmann) distribution for $s=0$, i.e.
${\mathbb{P}}_{\rm{eq}}\propto{\rm{diag}}({\rm{e}}^{\beta E_{\mathcal{C}}})$
(see also the comments after (S50)). This equilibrium solution always exists
for finite $N$, though in the glassy regime the time to reach it diverges with
$N$.

Explicitly, the similarity transformation reads

$\displaystyle\tilde{\mathbb{W}}(s)={\mathbb{P}}_{\rm{eq}}^{-1/2}\mathbb{W}(s){\mathbb{P}}_{\rm{eq}}^{1/2}$
(S132)

or in components

$\displaystyle\tilde{W}_{\mathcal{C}\mathcal{C^{\prime}}}(s)={\rm{e}}^{-\beta
E_{\mathcal{C}}/2}W_{\mathcal{C}\mathcal{C^{\prime}}}(s){\rm{e}}^{\beta
E_{\mathcal{C^{\prime}}}/2}\,.$ (S133)

Denoting by $\tilde{V}_{0}$ the normalised eigenvector corresponding to the
top eigenvalue of $\tilde{\mathbb{W}}(s)$, we have

$\displaystyle\psi(s)=\langle\tilde{V}_{0}|\tilde{\mathbb{W}}(s)|\tilde{V}_{0}\rangle$
(S134)

One can now deploy standard first–order perturbation theory (as used for
instance in Ref. [110]; the key feature is that the $s$-dependence of
$\tilde{V}_{0}$ does not contribute) to express the first derivative of the
free energy as

$\displaystyle\psi^{\prime}(s)$
$\displaystyle=\langle\tilde{V}_{0}|\tilde{\mathbb{W}}^{\prime}(s)|\tilde{V}_{0}\rangle$
(S135)
$\displaystyle=\sum_{\mathcal{C},\mathcal{C^{\prime}}}\tilde{V}_{0}(\mathcal{C})\tilde{W}^{\prime}_{\mathcal{C}\mathcal{C^{\prime}}}(s)\tilde{V}_{0}(\mathcal{C^{\prime}})$
(S136)
$\displaystyle=\sum_{\mathcal{C},\mathcal{C^{\prime}}}\tilde{V}_{0}(\mathcal{C}){\rm{e}}^{-\beta
E_{\mathcal{C}}/2}W^{\prime}_{\mathcal{C}\mathcal{C^{\prime}}}(s){\rm{e}}^{\beta
E_{\mathcal{C^{\prime}}}/2}\tilde{V}_{0}(\mathcal{C^{\prime}})$ (S137)

Now we use that the eigenvectors of the original biased master operator can be
written in terms of the symmetrized eigenvector as $\langle
U_{0}|=\langle\tilde{V}_{0}|{\mathbb{P}}_{\rm{eq}}^{-1/2}$ and
$|V_{0}\rangle={\mathbb{P}}_{\rm{eq}}^{1/2}|\tilde{V}_{0}\rangle$ (see
discussion after (S50) and Ref. [86] for an illustration in the case of the
unbiased Bouchad trap model). Additionally, we use the expression of the
master operator as given by equation (S7). This results in

$\displaystyle-\psi^{\prime}(s)=\sum_{\mathcal{C},\mathcal{C^{\prime}},\mathcal{C}\neq\mathcal{C^{\prime}}}U_{0}(\mathcal{C})W(\mathcal{C^{\prime}}\to\mathcal{C}){\rm{e}}^{-s}V_{0}(\mathcal{C^{\prime}})$
(S138)

and in the thermodynamic limit

$\displaystyle-\psi^{\prime}(s)$ $\displaystyle=\int
dE\,dE^{\prime}\,U_{0}(E)V_{0}(E^{\prime})W(E^{\prime}\to
E){\rm{e}}^{-s}\rho(E)$ (S139)

One can easily verify that the r.h.s. of equation (S139) corresponds to the
mean time–averaged activity by using the equation for the top right
eigenvector (eq. (S48)) to replace the integral over $E^{\prime}$. Then
equation (S139) transforms to

$\displaystyle-\psi^{\prime}(s)$ $\displaystyle=\int
dE\,U_{0}(E)V_{0}(E)(\psi(s)+r(E))$ (S140)

and one can recognize this as equation (14) of the manuscript.

For Barrat–Mézard transition rates, eq. (S139) becomes

$\displaystyle-\psi^{\prime}(s)=\int
dE\,dE^{\prime}\,\frac{U_{0}(E)V_{0}(E^{\prime})\rho(E){\rm{e}}^{-s}}{1+{\rm{e}}^{-\beta(E-E^{\prime})}}=\langle\bar{a}\rangle_{s}$
(S141)

and this expression was used to generate the data for Fig. 1(d) in the main
text.

##  .7 Numerical details

Here we provide a brief summary of how we generated the main figures of the
manuscript. The dynamical free energy in Fig. 1(c) for the active regime is
obtained as the top eigenvalue of a biased master operator $\mathbb{W}(s)$ of
size $N=2^{13}$. An instance of the master operator with elements given by eq.
(S7) can be generated after randomly sampling $N$ energies from to $\rho(E)$
and associate to each configuration one of these energies. To use first–order
perturbation theory for the determination of the mean time–averaged activity
we obtain the top right eigenvector using the Arnoldi iteration (as
implemented in the Arpack package [111]) and then exploit detailed balance to
obtain the top left–eigenvector (see discussion after eq. (S50)).

The numerical calculation of the mean time–dependent activity (see e.g. Fig.
2) proceeds in two steps. First, we find the biased probability distribution
at a given time (eq. (13), see e.g. Fig. 3 and Fig. S4) and then we compute
the mean auxiliary escape rate according to eq. (12). To obtain
$P_{s}(E,\tau)$ we start by discretizing the energy range $E\in[0,15]$ with a
grid of step size $\Delta E=0.02$, we also fix the total time (or trajectory
length) $t=10^{4}$ and consider a discrete time step $\Delta\tau=5\times
10^{-2}$. Then we discretize the forward and backward master equations for
$\mathbb{W}(s)$ (eq. (S38) and (Bringing together two paradigms of non-
equilibrium: Driven dynamics of aging systems), respectively), replacing the
time derivatives by finite differences. Finally, we use the initial conditions
$P_{\rm f}(E,0)={\rm{e}}^{-E}$ and $q(E,t)=1$ and then combine the numerical
solutions according to equation (13). To evaluate the mean time–dependent
activity (eq. (12)) it then remains to estimate the derivative of the
logarithm of the backwards factor. This is done on–the–fly during the
numerical integration, using again finite differences.