††footnotetext: *Equal contributions.

# Microscopic origin of excess wings in relaxation spectra of supercooled
liquids

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

###### Abstract

Glass formation is encountered in diverse materials. Experiments have revealed
that dynamic relaxation spectra of supercooled liquids generically become
asymmetric near the glass transition temperature, $T_{g}$, where an extended
power law emerges at high frequencies. The microscopic origin of this “wing”
remains unknown, and was so far inaccessible to simulations. Here, we develop
a novel computational approach and study the equilibrium dynamics of model
supercooled liquids near $T_{g}$. We demonstrate the emergence of a power law
wing in numerical spectra, which originates from relaxation at rare, localised
regions over broadly-distributed timescales. We rationalise the asymmetric
shape of relaxation spectra by constructing an empirical model associating
heterogeneous activated dynamics with dynamic facilitation, which are the two
minimal physical ingredients revealed by our simulations. Our work offers a
glimpse of the molecular motion responsible for glass formation at relevant
experimental conditions.

The formation of amorphous solids results from the rapid growth of the
structural relaxation time $\tau_{\alpha}$ of the supercooled liquid Berthier
and Biroli (2011). Molecular motion occurs on a timescale of about $10^{-10}$
s at the onset temperature of glassy behaviour but takes about $100$ s at the
experimental glass transition temperature $T_{g}$ Schmidtke  _et al._ (2012).
Over the last decades, dielectric, mechanical and light scattering experiments
kept developing to probe molecular motion over a broader frequency range with
increased accuracy Lunkenheimer  _et al._ (2000); Körber  _et al._ (2020);
Schmidtke  _et al._ (2013); Gainaru  _et al._ (2009); Flämig  _et al._ (2020);
Hecksher  _et al._ (2017). This progress reveals that the temperature
evolution of $\tau_{\alpha}$ is just the tip of the iceberg, as relaxation
spectra $\chi^{\prime\prime}(\omega)$ measured near $T_{g}$ exhibit relaxation
processes taking place over an extremely large frequency window Schneider  _et
al._ (2000); Lunkenheimer  _et al._ (2002); Dixon  _et al._ (1990); Leheny and
Nagel (1998). The overall shift of relaxation spectra is accompanied by an
equivalent broadening of about 12 decades, which is the other side of the same
coin. A microscopic explanation of these slow dynamics is at the heart of
glass transition research Berthier and Biroli (2011).

High-temperature spectra reflect near exponential relaxation in the picosecond
range, but low-$T$ spectra broaden into a two-step process with a stretched
exponential relaxation at low frequency $\omega\approx 1/\tau_{\alpha}$ and a
microscopic peak remaining at the picosecond timescale. In 1990, Nagel and
coworkers Menon  _et al._ (1992); Dixon  _et al._ (1990); Menon and Nagel
(1995); Leheny and Nagel (1998) showed that for a number of molecular liquids
the structural relaxation peak extends much further at high frequencies
$\omega\tau_{\alpha}\gg 1$ and transforms into a power law,
$\chi^{\prime\prime}(\omega)\sim\omega^{-\sigma}$, with a small exponent
$\sigma(T)\in\left[0.2,0.4\right]$ decreasing with temperature Menon and Nagel
(1995). Using logarithmic scales, this resembles a “wing” in “excess” of the
$\alpha$-peak. At $T_{g}$, the wing extends over the mHz-MHz range with an
amplitude about 100 times smaller than the $\alpha$-peak. A universal scaling
comprising the excess wing was proposed Menon and Nagel (1995), which can be
altered by additional microscopic processes Wu (1991); Ngai and Paluch (2004).
While this universality is debated Schneider  _et al._ (2000); Lunkenheimer
_et al._ (2002), the presence of an excess contribution often taking the form
of a wing is not Blochowicz  _et al._ (2003); Gainaru  _et al._ (2009).

Elucidating the nature of molecular motions responsible for the small signal
in these excess wings appears daunting. Yet, experiments managed to
characterise its heterogeneous nature Bauer  _et al._ (2013); Duvvuri and
Richert (2003) and aging properties Lunkenheimer  _et al._ (2005). So far,
computer simulations were unable to access the required range of equilibration
temperatures and timescales to even address the question. Physical
interpretations and empirical models have been proposed to explain the shape
of relaxation spectra. Some of them couple slow translational motion to an
“additional” degree of freedom (e.g., rotational) Diezemann  _et al._ (1999);
Domschke  _et al._ (2011). Others invoke spatially heterogeneous dynamics to
construct a broad distribution of timescales of static Sethna  _et al._
(1991); Stevenson and Wolynes (2010); Viot  _et al._ (2000); Chamberlin
(1999); Dyre and Schrøder (2000) or kinetic Berthier and Garrahan (2005)
origin. The winged asymmetric shape then requires specific physics, such as
geometric frustration Viot  _et al._ (2000), lengthscale-dependent dynamics
Chamberlin (1999), or dynamic facilitation Berthier and Garrahan (2005). With
specific choices, these approaches yield relaxation spectra comprising excess
wings, but direct microscopic investigations testing the underlying hypotheses
are still lacking.

Here, we show that computer simulations can now directly observe excess wings
and assess their microscopic origin. We take advantage of the recent swap
Monte Carlo algorithm Ninarello  _et al._ (2017) to efficiently produce
equilibrated configurations of a supercooled liquid with $\tau_{\alpha}\approx
100$ s. We observe their physical relaxation dynamics over ten decades in
time, up to $20$ ms. We are thus able to probe for the first time the
temperature and time regimes where excess wings are observed in experiments.
We report the emergence of a power law (a wing) in numerical spectra with the
same characteristics as in experiments. We demonstrate that it is caused by a
sparse population of localised regions, whose relaxation times are power law
distributed. These relaxed regions then coarsen by dynamic facilitation. We
construct an empirical model to illustrate how heterogeneous dynamics and
dynamic facilitation generically lead to asymmetric, winged relaxation
spectra.

We study size-polydisperse mixtures of $N$ soft repulsive spheres in two and
three dimensions, as described in the Methods section. These models are
representative computational glass-formers Berthier  _et al._ (2019a, 2017).
We use the swap Monte Carlo algorithm designed in Ref. Berthier  _et al._
(2019b) to generate $n_{s}\in\left[200,450\right]$ independent equilibrium
configurations at temperatures $T$ down to the extrapolated experimental glass
transition temperature $T_{g}$. Each equilibrium configuration is then taken
as the initial condition of a multi-CPU molecular dynamics (MD) simulation
(without swap). The $n_{s}$ independent simulations run for up to a simulation
time $t_{\textrm{max}}=1.5\times 10^{7}$ in $3d$ (one week on 2 CPUs for
$N=1200$). We push a few $2d$ simulations to unprecedentedly long times, up to
$t_{\textrm{max},2d}=6\times 10^{8}$, representing a computational time of
several months. By using the relaxation time at the onset of glassy dynamics
to relate numerical and experimental timescales, our longest simulations
translate into a physical time of about $20$ ms for systems having an
equilibrium relaxation time $\tau_{\alpha}\approx 10^{2}$ s. This strategy is
key to observe excess wings, which would otherwise be buried underneath the
structural relaxation in conventional approaches Yu  _et al._ (2017). The $2d$
and $3d$ models behave similarly, so we present quantitative results for the
$3d$ model ($N=1200$) in Figs. 1, 3 and illustrate the relaxation process in
Fig. 2 with $2d$ snapshots ($N=10000$), which are easier to visualize.
Quantitative results for the $2d$ model are provided in the Supplementary
Information (SI).

We investigate the spatio-temporal evolution of the relaxation dynamics using
averaged and particle-resolved dynamic observables. In $3d$, we measure the
self-intermediate scattering function $F_{s}(t)$, averaged over the $n_{s}$
independent runs. We define the relaxation time $\tau_{\alpha}$ by
$F_{s}(\tau_{\alpha})=e^{-1}$. In $2d$, collective long-ranged fluctuations
affect the measurement of $F_{s}(t)$. We instead focus on observables which
are blind to these fluctuations Illing  _et al._ (2017) and define
$\tau_{\alpha}$ via the bond-orientational correlation function Flenner and
Szamel (2015). In both two and three dimensions, we investigate the relaxation
process at the particle scale via the bond-breaking correlation $C_{B}^{i}(t)$
which quantifies the fraction of nearest neighbours lost by particle $i$ after
time $t$. Starting from $C_{B}^{i}(t=0)=1$, it decreases as rearrangements
take place close to particle $i$, and reaches zero when its local environment
is completely renewed. Precise definitions of the correlation functions are
provided in the Methods section.

To connect with experimental results obtained in the frequency domain, we
compute a dynamic susceptibility $\chi^{\prime\prime}(\omega)$ from a
distribution of relaxation times $G(\log\tau)$ Blochowicz  _et al._ (2003);
Berthier and Garrahan (2005)

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

where the distribution $G$ is related to the derivative of a time correlation
function, $G(\log t)\approx-\mathrm{d}F_{s}(t)/\mathrm{d}\log t$ in $3d$. We
use the bond-breaking correlation function instead of $F_{s}$ in $2d$. We
discuss the numerical evaluation of $\chi^{\prime\prime}$ in the Methods
section, whereas the discussion on the statistical noise and the comparison to
direct Fourier transforms are in the Supplementary Information, Section I and
Figure 2.

Figure 1:  Emergence of excess wings in a three-dimensional glass-former near
the glass transition temperature. (a) Self-intermediate scattering function
$F_{s}(t)$ at various temperatures. (b) Relaxation time $\tau_{\alpha}$
rescaled by its value $\tau_{o}$ at the onset temperature. Symbols are
directly measured data (squares), or obtained using time temperature
superposition (TTS, circles). A conservative Arrhenius extrapolation locates
$T_{g}=0.056$, where $\tau_{\alpha}(T_{g})/\tau_{o}=10^{12}$ (dashed line).
(c) Relaxation spectra for the same temperatures as in panel (a). The dashed
lines represent the estimated $\alpha$-peaks. Close to $T_{g}$, the spectra
lie above the $\alpha-$peak and display a power law signal with an exponent
$\sigma\approx 0.38$ (full line), in quantitative agreement with excess wings
observed experimentally.

We start by presenting equilibrium measurements of $F_{s}(t)$ in $3d$ in Fig.
1(a), concentrating on the unexplored low-$T$ regime below the mode-coupling
crossover $T_{\rm mct}\approx 0.095$. The latter is determined by a power law
fit of $\tau_{\alpha}(T)$ in the range $\tau_{\alpha}/\tau_{o}<10^{3}$, where
$\tau_{o}\approx 3$ is the value of $\tau_{\alpha}$ at the onset temperature
$T_{o}\approx 0.20$ Ninarello  _et al._ (2017). At all temperatures, the
correlations display a fast initial decay near $t\approx\tau_{o}$, due to fast
dynamical processes. At larger times, we observe a much slower decay to zero.
As $T$ decreases, the relaxation time grows and eventually exits the
numerically accessible time window. At the lowest investigated temperatures
near $T_{g}$, correlations appear almost constant over more than 7 decades in
time, suggesting near-complete dynamic arrest. We recall that thanks to the
swap algorithm, all measurements reflect genuine equilibrium dynamics, even
when $\tau_{\alpha}$ is larger than the simulated time by many orders of
magnitude.

Our strategy allows us to directly observe the $\alpha$-relaxation when
$\tau_{\alpha}<t_{\textrm{max}}$, equivalently $\tau_{\alpha}/\tau_{o}\lesssim
5\times 10^{6}$ down to $T=0.0755$, see Figs. 1(a,b). In this regime, the
relaxation is well-described by a stretched exponential
$F_{0}e^{-(t/\tau_{\alpha})^{\beta}}$ with an almost constant stretching
exponent $\beta\approx 0.56$, the amplitude $F_{0}$ modestly changing with
temperature. We use this time temperature superposition (TTS) property to
estimate $\tau_{\alpha}$ for $0.07\leq T\leq 0.0755$, where the decorrelation
of $F_{s}(t)$ is sufficient Berthier and Ediger (2020), and obtain
$\tau_{\alpha}$ over roughly 2 additional decades, see Fig. 1(b). We finally
use an Arrhenius law to extrapolate $\tau_{\alpha}$ over 4 more decades to get
a safe lower bound for the experimental glass temperature $T_{g}\approx
0.056$, defined by $\tau_{\alpha}(T_{g})/\tau_{o}=10^{12}$ Ninarello  _et al._
(2017), see the Methods section for details.

The corresponding relaxation spectra are shown in Fig. 1(c) for the $3d$
model. They all display a peak at high frequency $\omega\approx 1/\tau_{o}$,
corresponding to the short-time decay of $F_{s}(t)$. A low-frequency peak near
$\omega\approx 1/\tau_{\alpha}$ is also visible. As $T$ decreases, this
$\alpha$-peak shifts to lower frequencies and eventually exits the accessible
frequency window. When the $\alpha$-peak is not directly measured, we
extrapolate its shape by inserting the above stretched exponential form for
$F_{s}(t)$ into Eq. (1). We use $\beta=0.56$, $\tau_{\alpha}$ given by the
Arrhenius extrapolation, and a constant $F_{0}$. The tiny temperature
dependence of $F_{0}$ is immaterial on the logarithmic scale of Fig. 1(c). The
resulting $\alpha$-peaks are shown in Fig. 1(c) with dashed lines that
smoothly merge into the measured data at the highest temperatures, validating
our procedure.

As $T$ decreases, the measured susceptibility and the $\alpha$-peak deviate
increasingly from one another, the data being systematically in excess of the
$\alpha$-peak. Since the Arrhenius extrapolation underestimates
$\tau_{\alpha}$, this excess is (at worst) slightly underestimated and cannot
be accounted by a vertical shift which would require unphysical values of
$F_{0}$ and $\beta$. At the lowest $T$, where the $\alpha$-peak no longer
interferes with the measurements, the spectra are well described by a power
law $\chi^{\prime\prime}(\omega)\sim\omega^{-\sigma}$ at low frequencies, with
an exponent $\sigma\approx 0.38$ slightly decreasing with $T$, and an
amplitude about 100 times smaller than the $\alpha$-peak. The relaxation
spectra of the $2d$ model in Supplementary Figure 3 exhibit similar features
with an exponent $\sigma_{2d}\approx 0.45$, which is quite close to the one
found in $3d$. In our simulations, the measured spectra do not exhibit a
secondary peak separated from the $\alpha$-relaxation, and cannot be
interpreted using an additive $\beta$-process Yu  _et al._ (2017). Therefore,
close to $T_{g}$, the numerical spectra follow a power law over a similar
frequency range, with a similar exponent and a similar amplitude as the excess
wings obtained experimentally, suggesting that simulated glass-formers display
excess wings resembling observations in molecular liquids.

Figure 2: Visualisation of spatially heterogeneous and facilitated dynamics.
Relaxation in the $2d$ system at $T_{2d}=0.09$ with
$\tau_{\alpha}/\tau_{o}=10^{8}$. Frames are logarithmically spaced between
$t=2\times 10^{-3}\tau_{\alpha}$ (top left) and $t=0.6\tau_{\alpha}$ (bottom
right) from left to right and top to bottom. Particles are coloured according
to $C_{B}^{i}(t)$ from blue [immobile, $C_{B}^{i}(t)=1$] to red [relaxed,
$C_{B}^{i}(t)=0$]. The linear size of the simulation box is 100.

We take advantage of the atomistic resolution offered by simulations to
explore the microscopic origin of excess wings and provide a physical
interpretation of the spectral shapes.We illustrate the relaxation dynamics
with $2d$ snapshots, which are easier to render and interpret. We confirm that
the same mechanisms are observed in $3d$. In Fig. 2 we show $2d$ snapshots
illustrating how structural relaxation proceeds at a temperature $T_{2d}=0.09$
(we estimate $T_{g,2d}\approx 0.07$) for which $\tau_{\alpha}/\tau_{o}\approx
10^{8}$, corresponding to around 10 ms in physical time. This temperature is
the lowest for which the $\alpha$-relaxation can be observed in the numerical
window, and is considerably lower than the mode-coupling crossover near
$T_{{\rm mct},2d}\approx 0.12$. Images are shown at logarithmically-spaced
times $t$ in the range $t/\tau_{\alpha}\in\left[10^{-3},1\right]$. Particles
are coloured according to $C_{B}^{i}(t)$: red particles have relaxed, blue
ones have not. We present in Fig. S3 the relaxation spectrum measured at this
temperature.

For $t\ll\tau_{\alpha}$, relaxation starts at a sparse population of localised
regions which emerge independently throughout the sample over broadly
distributed times. This conclusion holds over a large range of temperatures
down to $T_{g}$ in both $d=2,3$. As time increases, newly relaxed regions
continue to appear, but a second mechanism becomes apparent in Fig. 2 as
regions that have relaxed in one frame typically appear larger in the next.
This growth of relaxed regions in Fig. 2 is the signature of dynamic
facilitation Chandler and Garrahan (2010). More precisely, we observe that
from one frame to the next, relaxation events keep accumulating at similar
locations, which results in mobile particles undergoing multiple relaxations
and mobility propagating to nearby particles. Also, the slowest regions are
typically “invaded” at $t\gg\tau_{\alpha}$ from their faster boundaries.
Dynamic facilitation has been identified before at high temperatures above the
mode-coupling crossover Chandler and Garrahan (2010); Keys  _et al._ (2011);
Vogel and Glotzer (2004). Our investigations show that it becomes a central
physical mechanism for structural relaxation near $T_{g}$.

We concentrate on the early times where power law spectra are observed.
Visualisation suggests that clusters of relaxed particles appear at sparse
locations. We now establish that these early relaxation events are responsible
for the excess wing. To this end, we define mobile ($C_{B}^{i}<0.55$) and
immobile ($C_{B}^{i}\geq 0.55$) particles; the threshold value near 0.5 is
determined requiring self-consistency with alternative mobility definitions
based on displacements. We identify connected clusters of mobile particles by
performing a nearest neighbour analysis (details in the Methods section), and
investigate the statistical properties of relaxed clusters. In particular, we
find that the excess wing regime at $t/\tau_{\alpha}\ll 1$ is dominated by the
appearance of new clusters, whereas the growth of existing clusters dominates
at later times. We report in Fig. 3 the distribution $\Pi(\tau)$ of waiting
times $\tau$ for the appearance of new clusters in $3d$. For $T\leq 0.07$, we
cannot measure the entire distribution, which is thus determined up to an
uninteresting prefactor. The corresponding $2d$ results are shown in
Supplementary Figure 4.

Figure 3: Microscopic origin of excess wings. Waiting-time distribution of
newly relaxing clusters in $3d$ from $T_{\rm mct}$ (right) to $T_{g}$ (left),
with the same color code as Fig. 1(a). Approaching $T_{g}$, the distributions
develop a power law tail at $\tau\ll\tau_{\alpha}$, with an exponent $0.38$
that directly accounts for the excess wings in the spectra of Fig. 1(c).

At the highest investigated temperature, near $T_{\rm mct}$, the distribution
$\Pi(\log_{10}\tau)$ in Fig. 3 is already very broad, with clusters appearing
as early as $10^{-4}\tau_{\alpha}$. The distribution peaks near
$0.1\tau_{\alpha}$, when dynamic facilitation starts to dominate, and has a
cutoff around $10\tau_{\alpha}$. As $T$ decreases below the mode-coupling
crossover, a power law tail emerges at $\tau\ll\tau_{\alpha}$. For $T\leq
0.07$, the power law extends over at least 6 decades, with a nearly constant
exponent $\Pi(\log_{10}\tau)\sim\tau^{0.38}$ for the $3d$ model. The
relaxation of localised clusters at early times is extremely broadly
distributed, presumably stemming from an equally broad distribution of
activation energies.

Remarkably, if we plug the measured distribution of waiting times in Eq. (1),
a power law $\Pi(\log_{10}\tau)\sim\tau^{0.38}$ directly translates into a
power law $\chi^{\prime\prime}(\omega)\sim\omega^{-0.38}$ in the spectra,
which is thus valid for $\omega\tau_{\alpha}\gg 1$. The agreement with the
data in Fig. 1(c) is therefore quantitative. A similar agreement is found in
$2d$ with the exponent $\sigma_{2d}=0.45$, see Supplementary Figure 4. This
analysis demonstrates that the high-frequency power law in
$\chi^{\prime\prime}(\omega)$ stems from the relaxation of a sparse population
of clusters characterised by a broad distribution of relaxation times.

This microscopic view of the power law wing alone does not explain why it
appears in excess of the $\alpha$-peak observed at larger times when dynamic
facilitation sets in. To explain this point, we construct an empirical model
based on our numerical observations. We first imagine that the liquid can be
decomposed into independent domains characterised by a local relaxation time,
see Fig. 4(a). This heterogeneous viewpoint is mathematically captured by trap
models Dyre (1987); Bouchaud (1992). To introduce dynamic facilitation as the
second key ingredient, we construct a facilitated trap model, assuming that a
given local relaxation event may now affect the state of the other traps, see
Fig. 4(b). To provide a qualitative, generic description of relaxation
spectra, we analyse the simplest version of such a model and assume, in a
mean-field spirit, that dynamic facilitation equally affects all traps. A more
local version was designed in Refs. Rehwald  _et al._ (2010); Rehwald and
Heuer (2012) for different purposes.

Figure 4: Facilitated trap model generically predicts asymmetric winged
relaxation spectra. (a) The liquid is modeled as a collection of traps with
energies $E$, distributed according to $\rho(E)$. (b) Relaxation is thermally
activated and affects the energy of the other traps by a random amount,
proportional to $\Delta$. (c) Relaxation spectra $\chi^{\prime\prime}(\omega)$
in the absence ($\Delta=0$) and presence ($\Delta=0.05$) of dynamic
facilitation at $\alpha=1.1$ and $T=0.629$. Dynamic facilitation compresses
(arrow) the low-frequency part of the underlying spectrum, giving rise to a
sharper $\alpha$-peak well-fitted by the spectrum of a stretched exponential
(dashed). The high-frequency part of the underlying spectrum, unaffected by
facilitation, is well described by a power law $\omega^{-0.2}$ (line). The
excess wing thus corresponds to the beginning of the relaxation process.

We consider $N$ traps with energy levels $E>0$ drawn from a distribution
$\rho(E)$, and assume activated dynamics. The energy $E$ of a trap is renewed
after a Poisson-distributed timescale of mean $\langle\tau(E)\rangle=e^{E/T}$.
Since deep traps take much longer to relax than shallow ones, the system is
dynamically heterogeneous. Following Ref. Arkhipov and Baessler (1994), we use
$\rho(E)\propto e^{-E^{\alpha}}$, with $\alpha\in\left[1,2\right]$ to smoothly
interpolate between the much-studied Gaussian Dyre (1987); Rehwald  _et al._
(2010) and exponential Bouchaud (1992) distributions. Dynamics at temperature
$T$ leads to the equilibrium energy distribution
$P_{\text{eq}}(T,E)\propto\rho(E)e^{E/T}$. Whenever a trap relaxes, the energy
of all other traps is shifted by a random amount uniformly distributed in the
interval $\left[-\frac{\Delta}{\sqrt{N}},\frac{\Delta}{\sqrt{N}}\right]$,
using a Metropolis filter to leave the equilibrium distribution
$P_{\text{eq}}$ unchanged. This coupling between traps mimics dynamic
facilitation Rehwald  _et al._ (2010). The relaxation spectra
$\chi^{\prime\prime}(\omega)$ is computed either analytically ($\Delta=0$), or
by simulating the facilitated model ($\Delta>0$).

The model is specified by two parameters $(\alpha,\Delta)$, for which
equilibrium dynamics can be studied at any temperature $T$. We have
systematically investigated this parameter space, and find spectra with
quantitative differences but generic features Scalliet  _et al._ (2021). In
Fig. 4(c), we select $(\alpha=1.1,\Delta=0.05)$ at $T=0.629$ for aesthetic
reasons, as this produces a spectrum qualitatively resembling experimental and
numerical ones close to $T_{g}$. Fitting the $\alpha$-peak to the frequency
representation of a stretched exponential reveals an excess wing at high
frequencies. However, in the absence of dynamic facilitation ($\Delta=0$) one
obtains the blue spectrum, with the same high-frequency behaviour, but which
extends much further at low frequencies. Indeed, without facilitation each
trap relaxes independently, and the equilibrium distribution $P_{\rm eq}$
determines the dynamic spectrum, which is broad and relatively symmetric. In
the presence of facilitation, $\Delta>0$, shallow traps still relax
independently and are essentially unaffected. Crucially, deep traps now
receive small kicks whenever a shallow trap relaxes, and their energies slowly
diffuse towards the most probable value. This accelerates their relaxation,
which eventually affects the tail of the relaxation time distribution. As a
result, dynamic facilitation “compresses” the low-frequency part of the
underlying spectrum (blue), as hinted in Ref. Xia and Wolynes (2001), and
highlighted by the arrow in Fig. 4(c). We thus interpret the winged,
asymmetric spectrum as a broad underlying distribution of relaxation
timescales (well described by a power law at early times) compressed by
dynamic facilitation at long times. Ironically, in our picture, the
$\alpha$-peak itself is in “excess” of a much broader underlying time
distribution with a high-frequency power law shape. In this view, the excess
wing forms an integral part of the structural relaxation.

Our study frontally attacks a central question regarding the relaxation
dynamics of supercooled liquids near the experimental glass transition and
paves the way for many more studies of a totally unexplored territory now made
accessible to modern computer studies. Enlarging further the family of
available computer glass-formers would also help filling the gap with the more
complex molecular systems studied experimentally.

Acknowledgments– We thank G. Biroli, M. Ediger and J. Kurchan for discussions,
and S. Nagel for detailed explanations about experiments. Some simulations
were performed at MESO@LR-Platform at the University of Montpellier. This work
was supported by a grant from the Simons Foundation (#454933, L.B.), the
European Research Council under the EU’s Horizon 2020 Program, Grant No.
740269 (C.S.), a Herchel Smith Postdoctoral Research Fellowship (C.S.), a
Ramon Jenkins Research Fellowship from Sidney Sussex College, Cambridge (C.S.)
and Capital Fund Management – Fondation pour la Recherche (B.G.).

## Methods

### Glass-forming computer models

We study a non-additive, continuously polydisperse mixture of spherical
particles of equal mass $m$ in two and three dimensions ($d=2,3$) Ninarello
_et al._ (2017). Two particles $i$ and $j$, at a distance $r_{ij}$ from one
another interact via the repulsive potential

$v(r_{ij})=\epsilon\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}+c_{0}+c_{2}\left(\frac{r_{ij}}{\sigma_{ij}}\right)^{2}+c_{4}\left(\frac{r_{ij}}{\sigma_{ij}}\right)^{4},$
(2)

if $r_{ij}/\sigma_{ij}<x_{c}=1.25$. The constants
$c_{0}=-28\epsilon/x_{c}^{12},c_{2}=48\epsilon/x_{c}^{14},c_{4}=-21\epsilon/x_{c}^{16}$
ensure continuity of the potential and its first two derivatives at the cutoff
$x_{c}$. The particles’ diameters $\sigma_{i}$ are distributed from
$\mathcal{P}(\sigma)=\mathcal{A}/\sigma^{3}$ with $\mathcal{A}$ a
normalisation constant, $\sigma_{\mathrm{max}}/\sigma_{\mathrm{min}}=2.219$.
We use the average diameter $\overline{\sigma}$ as unit length, $\epsilon$ as
unit energy (the Boltzmann constant is set to unity) and
$\sqrt{m\overline{\sigma}^{2}/\epsilon}$ as unit time. In these units,
$\sigma_{\mathrm{min}}=0.73$ and $\sigma_{\mathrm{max}}=1.62$. We employ a
non-additive cross-diameter rule
$\sigma_{ij}=0.5(\sigma_{i}+\sigma_{j})(1-0.2|\sigma_{i}-\sigma_{j}|)$ to
avoid fractionation and crystallization at low temperature Ninarello  _et al._
(2017). We simulate the glass-forming model at number density of particles
$\rho=N/L^{d}=1$ in a cubic/square box of linear size $L$ using periodic
boundary conditions. We consider various system sizes: $N=1200,10000$ in $3d$
and $N=2000,10000$ in $2d$.

### Preparation of equilibrated configurations

The model glass-forming liquid is efficiently simulated at equilibrium with
the swap Monte Carlo algorithm. We employ the hybrid swap Monte
Carlo/Molecular Dynamics algorithm implemented in the LAMMPS package ($2d/3d$)
or homemade code ($3d$), with optimal parameters, as described in Ref.
Berthier  _et al._ (2019b). We prepare $n_{s}\in\left[200,450\right]$
independent equilibrated configurations at temperatures down to the
experimental glass transition temperature.

### Molecular dynamics simulations

The equilibrium configurations generated by the swap algorithm are used as
initial conditions for standard molecular dynamics (MD) simulations with
integration time step equal to $0.01$. In $3d$, we run conventional MD (NVE)
simulations and NVT simulations in $2d$ using a Nosé-Hoover thermostat. The
simulations are either run using a homemade MD code or with the LAMMPS
package, which allows us to run multi-CPU simulations and perform extremely
long runs for relatively large systems (e.g., two months on 24 CPUs for Fig.
2).

### Relating experimental and numerical timescales

We measure the relaxation time at the onset of glassy dynamics as reference
time, and use this value to translate numerical timescales into experimental
ones. In experiments, many supercooled liquids have $\tau_{o}\approx 10^{-10}$
s. We measure $\tau_{o}\approx 3$ in $2d$ and $3d$ simulations. In $3d$, the
longest simulation time is $t_{\textrm{max}}=1.5\times 10^{7}=5\times
10^{6}\tau_{o}$. We therefore simulate the equilibrium relaxation at $T_{g}$
over $0.5$ ms. In $2d$, we ran monthslong simulations to reach
$t_{\textrm{max}}^{2d}=6\times 10^{8}=2\times 10^{8}\tau_{o}$. Our numerical
approach therefore allows us to observe the equilibrium dynamics over $20$ ms
at $T_{g}$, which is a giant leap forward in equilibrium simulations of
supercooled liquids.

### Average dynamic observables

In $3d$, we monitor the relaxation dynamics via the self-intermediate
scattering function

$F_{s}(t)=\left\langle\frac{1}{N}\sum_{i=1}^{N}\cos\left[\bm{q}\cdot\delta\bm{r}_{i}(t)\right]\right\rangle_{\bm{q},n_{s}},$
(3)

where $\delta\bm{r}_{i}(t)$ is the displacement of particle $i$ over time $t$.
The brackets indicate the ensemble average over $n_{s}$ independent runs along
with an angular average over wavevectors with $|\bm{q}|=6.9$ (first peak in
the total structure factor).

In $2d$, collective long-ranged fluctuations give rise to a spurious
contribution to the displacements of particles Illing  _et al._ (2017) which
affects the measurement of $F_{s}(t)$ and makes it ill-suited to capture the
glassy slowdown. We instead study the dynamics through the evolution of the
local environment of particles, instead of their displacements. We define a
bond-orientational correlation function $C_{\Psi}(t)$ Flenner and Szamel
(2015). We introduce the six-fold bond-orientational order parameter of
particle $i$

$\Psi_{i}(t)=\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}e^{\mathrm{i}6\theta_{ij}(t)},$
(4)

where $n_{i}$ is the number of neighbours of $i$ at time $t$. Neighbours are
particles $j$ with $r_{ij}<1.45$ (first minimum in the radial distribution
function). Alternative definitions of neighbours, e.g., via Voronoi
tessellation or solid-angle based method van Meel  _et al._ (2012), lead to
the same quantitative results. Here $\theta_{ij}(t)$ is the angle between the
$x$-axis and the axis connecting $i$ and $j$ at time $t$, without loss of
generality thanks to rotational invariance. The bond-orientational correlation
function is defined as

$C_{\Psi}(t)=\left\langle\frac{\sum_{i}\Psi_{i}(t)\left[\Psi_{i}(0)\right]^{*}}{\sum_{i}|\Psi_{i}(0)|^{2}}\right\rangle_{n_{s}},$
(5)

where the brackets denote the ensemble average over $n_{s}$ independent runs,
and the star is the conjugate complex. In $2d$, we define the relaxation time
via $C_{\Psi}(\tau_{\alpha})=e^{-1}$.

### Mobility at the single-particle level

When analysing the mobility at the single-particle level, we first need a
criterion to distinguish between mobile and immobile particles. In $3d$, we
have considered several mobility definitions which all give quantitatively
similar results. The first mobility definition is based on displacements. To
remove fast dynamical processes, we use the conjugate-gradient method and find
the inherent structure (IS) of a configuration at time $t$,
$\\{\bm{r}_{i}^{IS}(t)\\}$. Particle $i$ is defined as mobile at time $t$ if
$|\bm{r}_{i}^{IS}(t)-\bm{r}_{i}^{IS}(0)|>0.8$ Schrøder  _et al._ (2000). This
cutoff is between the first minimum and the second maximum of the self part of
the van Hove function
$G_{s}(r,t)=\langle\delta(r-|\bm{r}_{i}^{IS}(t)|)\rangle_{i,n_{s}}$ in the
time regime where $F_{s}(t)$ is almost constant. This first mobility
definition is however not convenient in $2d$ because of the collective long-
ranged fluctuations which affect the translational dynamics.

A second mobility definition is based on changes in the particle’s local
environment. At time $t=0$, we find the number $n_{i}$, and identity of
particle $i$’s neighbours, defined as particles $j$ with
$r_{ij}/\sigma_{ij}<1.485$ in $3d$ (1.3 in $2d$), corresponding to the first
minimum in the rescaled pair correlation function $g(r_{ij}/\sigma_{ij})$. We
define the bond-breaking correlation as the fraction of remaining neighbours
at time $t$

$C_{B}^{i}(t)=\frac{n_{i}(t|0)}{n_{i}},$ (6)

where $n_{i}(t|0)$ is the number of particles neighbour of $i$ at $t=0$ and
still neighbour at $t$. To avoid short time oscillations in $C_{B}^{i}$ caused
by particles frequently exiting/entering the shell defining neighbours, we use
a slightly larger cutoff to define neighbours at $t>0$, namely
$r_{ij}/\sigma_{ij}<1.7$ (in $d=2,3$). We compute the bond-breaking
correlation function

$C_{B}(t)=\left\langle\frac{1}{N}\sum_{i=1}^{N}C_{B}^{i}(t)\right\rangle_{n_{s}},$
(7)

averaged over $n_{s}$ independent runs.

A particle is defined as mobile at $t$ if $C_{B}^{i}(t)<0.55$, i.e., if it has
lost half of its initial neighbours. The cutoff value ensures that the set of
particles identified as mobile in this way significantly overlap with that
identified via the displacement criterion. We then introduce clusters of
mobile particles. Two particles $i$ and $j$ mobile at time $t$ belong to the
same cluster if $r_{ij}<1.5$ in $3d$ and 1.4 in $2d$, close to the first
minimum of $g(r)$.

### Relevant temperature scales

We determine three temperature scales relevant to the glassy slowdown: the
onset temperature of glassy dynamics $T_{o}$, the mode-coupling crossover
temperature $T_{\rm mct}$ below which conventional MD simulations cannot reach
equilibrium, and the extrapolated experimental glass transition temperature
$T_{g}$. In $3d$, $T_{o}=0.2,T_{\rm mct}=0.095,T_{g}=0.056$. In $2d$,
$T_{o,2d}=0.2,T_{\rm mct,2d}=0.12,T_{g,2d}=0.07$. We fit the high-temperature
$\tau_{\alpha}$ data to an Arrhenius law, and identify the onset $T_{o}$ as
the temperature below which $\tau_{\alpha}$ is super-Arrhenius. We note
$\tau_{o}=\tau_{\alpha}(T_{o})$. The mode-coupling crossover temperature
$T_{\rm mct}$ is obtained by fitting the data with a power law
$\tau_{\alpha}(T)\propto(T-T_{\rm mct})^{-\gamma}$ in the regime
$0\leq\log_{10}(\tau_{\alpha}/\tau_{o})\leq 3$ Götze (2008), with $\gamma=2.7$
and $2.5$ in $d=2,3$ respectively. Given that
$\log_{10}(\tau_{\alpha}/\tau_{o})\approx 4$ at $T_{\rm mct}$, this
temperature delimits the regime $T>T_{\rm mct}$ where MD alone can reach
equilibrium, from the regime $T<T_{\rm mct}$ where the swap algorithm is
needed to perform equilibrium simulations. The experimental glass transition
temperature $T_{g}$ is defined by
$\log_{10}(\tau_{\alpha}(T_{g})/\tau_{o})=12$. In $3d$, the longest simulation
time is $t_{\textrm{max}}=1.5\times 10^{7}=5\times 10^{6}\tau_{o}$, so we can
directly access $\log_{10}(\tau_{\alpha}/\tau_{o})\lesssim 7$. We thus need to
extrapolate our data over 5 decades to locate $T_{g}$. We increase the
accuracy of the extrapolation by using time-temperature superposition (TTS),
which is well-obeyed in our model Berthier and Ediger (2020). In the
temperature regime where correlation functions reach $e^{-1}$, the second step
of the relaxation is well-fitted by a stretched exponential
$F_{0}e^{-(t/\tau_{\alpha})^{\beta}}$. The stretching exponent $\beta\simeq
0.56$ in $3d$ (in $2d$, $\beta\simeq 0.6$ for $C_{\Psi}$ and $\beta\simeq
0.67$ for $C_{B}$) is almost temperature-independent, and the amplitude
$F_{0}$ slightly increases with decreasing temperature. Fixing $\beta$, we
estimate $\tau_{\alpha}$ at temperatures where decorrelation is sufficient to
perform accurate TTS, extending our measurements over $\sim 2$ decades. We
extrapolate $\tau_{\alpha}$ over the 4 remaining decades using an Arrhenius
fit $\tau_{\alpha}(T)\propto e^{E_{A}/T}$ with $E_{A}=2.67$ in $3d$ ($2.97$ in
$2d$), and locate $T_{g}$. Importantly, the Arrhenius extrapolation is a safe
choice as it at worst underestimates relaxation times.

### Computation of relaxation spectra

The computation of relaxation spectra $\chi^{\prime\prime}(\omega)$ first
requires to differentiate the correlation function with respect to the
logarithm of time. We use a first-order finite difference approximation.
Namely, if configurations are stored at logarithmically-spaced times
$\\{t_{k}\\}_{k=1\dots n}$, we have for $k>1$

$\frac{\mathrm{d}F_{s}(t_{k})}{\mathrm{d}\log
t}=\frac{F_{s}(t_{k})-F_{s}(t_{k-1})}{\log(t_{k})-\log(t_{k-1})}.$ (8)

The integral in Eq. (1) is then evaluated by

$\chi^{\prime\prime}(\omega)=-\sum_{k=2}^{n}\frac{\mathrm{d}F_{s}(t_{k})}{\mathrm{d}\log
t}\frac{\omega t_{k}}{1+(\omega
t_{k})^{2}}\log\left(\frac{t_{k}}{t_{k-1}}\right).$ (9)

We use the bond-breaking correlation function $C_{B}$ instead of $F_{s}$ in
$2d$. In the Supplementary Information, we discuss errors which arise from
computing the spectrum when $F_{s}$ does not decay to zero. We also discuss
issues related to statistical noise and the comparison to direct Fourier
transforms.

### Trap model

We consider traps with energy levels $E>0$ drawn from the exponential power
distribution

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

and take $E_{0}=1$ in the following. We assume that dynamics at temperature
$T$ is thermally activated. The energy $E$ of a trap is renewed after a
Poisson-distributed timescale of mean $\langle\tau(E)\rangle=e^{E/T}$. The
equilibrium energy distribution at temperature $T$ is

$P_{\text{eq}}(T,E)=\frac{\rho(E)e^{E/T}}{Z(T)}~{}~{}\text{where}~{}~{}Z(T)=\int_{0}^{\infty}\mathrm{d}E\rho(E)e^{E/T}.$
(11)

We monitor relaxation dynamics by computing the average persistence function
$p(t)$. In the absence of dynamic facilitation, the persistence can be
directly computed

$p(t)=\int_{0}^{\infty}\mathrm{d}EP_{\text{eq}}(T,E)e^{-t/\langle\tau(E)\rangle}.$
(12)

In the absence of dynamic facilitation, the average persistence is evaluated
using Mathematica (NIntegrate, working precision 30). We then calculate the
relaxation spectrum $\chi^{\prime\prime}(\omega)$ by following the procedure
described previously, replacing $F_{s}(t)$ with the persistence $p(t)$. We
compute the persistence $p(t)$ over a time interval large enough to observe
full decorrelation, $[10^{-10},10^{70}]$ for $\alpha=1.1$, $T=0.629$, and
minimise errors in the relaxation spectrum.

### Simulations of the facilitated trap model

We consider a system composed of $N$ traps. We initialise the simulation with
an equilibrium condition by sampling the traps’ energies directly from the
equilibrium distribution $P_{\text{eq}}(T,E)$. Since the cumulative
probability distribution of energies $\mathcal{C}_{\text{eq}}$ cannot be
computed explicitly, we use Mathematica to evaluate it, and to numerically
construct the reciprocal function $\mathcal{E}={{\mathcal{C}}_{\rm eq}}^{-1}$.
For each of the $N$ traps, we generate $X$ uniformly distributed in $[0,1]$,
and assign it an energy $E=\mathcal{E}(X)$. This procedure generates an
initial condition in equilibrium. Each trap is assigned a renewal time
exponentially distributed, with mean $e^{E/T}$. We initialise the persistence
$p_{i}(t=0)$ of all traps to one.

The dynamics proceeds as follows. First, we identify the trap $i_{o}$ with the
smallest renewal time $\tau_{\text{min}}$, which will relax first. We update
all other traps by subtracting $\tau_{\text{min}}$ to their renewal time
$\tau_{i}$. When the trap $i_{o}$ relaxes, its persistence is set to zero,
$p_{i_{o}}=0$ and we give it a new energy value sampled from $\rho(E)$, and a
new renewal time, as described above.

This relaxation event then affects all other traps. We attempt to displace
their energy by a random amount $\delta E$ (different for each trap) uniformly
distributed in
$\left[-\frac{\Delta}{\sqrt{N}},\frac{\Delta}{\sqrt{N}}\right]$: $E\rightarrow
E^{\prime}=E+\delta E$. The scaling with $N$ ensures that the resulting
dynamics is independent on $N$. We then accept or reject this attempt in order
to leave the equilibrium probability distribution $P_{\text{eq}}$ unchanged.
To this end, we introduce an effective potential $V=-T\log P_{\text{eq}}$, and
compute the change in effective potential $\delta
V=T(E^{\prime\alpha}-E^{\alpha})-\delta E$. We then use the Metropolis filter:
if $\delta V<0$, the change in energy is accepted, otherwise, it is accepted
with probability $\exp(-\delta V/T)$. When accepted, we pick a new renewal
time exponentially distributed with average $e^{E^{\prime}/T}$. When the move
is completed, we again determine which of the traps is the next one to relax,
and proceed as before.

We measure the average persistence
$p(t)=\left\langle\sum_{i}p_{i}(t)/N\right\rangle$, where the brackets
indicate average over independent runs, and where the sum runs over all traps.
We simulate the dynamics of the model until the total persistence is equal to
zero.

## References

  * Berthier and Biroli (2011) L. Berthier and G. Biroli, Reviews of modern physics 83, 587 (2011).
  * Schmidtke  _et al._ (2012) B. Schmidtke, N. Petzold, R. Kahlau, M. Hofmann, and E. Rössler, Physical Review E 86, 041507 (2012).
  * Lunkenheimer  _et al._ (2000) P. Lunkenheimer, U. Schneider, R. Brand, and A. Loid, Contemporary Physics 41, 15 (2000).
  * Körber  _et al._ (2020) T. Körber, R. Stäglich, C. Gainaru, R. Böhmer, and E. A. Rössler, The Journal of Chemical Physics 153, 124510 (2020).
  * Schmidtke  _et al._ (2013) B. Schmidtke, N. Petzold, R. Kahlau, and E. A. Rössler, The Journal of Chemical Physics 139, 084504 (2013).
  * Gainaru  _et al._ (2009) C. Gainaru, R. Kahlau, E. A. Rössler, and R. Böhmer, The Journal of chemical physics 131, 184510 (2009).
  * Flämig  _et al._ (2020) M. Flämig, M. Hofmann, N. Fatkullin, and E. Rössler, The Journal of Physical Chemistry B 124, 1557 (2020).
  * Hecksher  _et al._ (2017) T. Hecksher, D. H. Torchinsky, C. Klieber, J. A. Johnson, J. C. Dyre, and K. A. Nelson, Proceedings of the National Academy of Sciences 114, 8710 (2017).
  * Schneider  _et al._ (2000) U. Schneider, R. Brand, P. Lunkenheimer, and A. Loidl, Phys. Rev. Lett. 84, 5560 (2000).
  * Lunkenheimer  _et al._ (2002) P. Lunkenheimer, R. Wehn, T. Riegger, and A. Loidl, Journal of Non-Crystalline Solids 307-310, 336 (2002).
  * Dixon  _et al._ (1990) P. K. Dixon, L. Wu, S. R. Nagel, B. D. Williams, and J. P. Carini, Phys. Rev. Lett. 65, 1108 (1990).
  * Leheny and Nagel (1998) R. L. Leheny and S. R. Nagel, Journal of Non-Crystalline Solids 235-237, 278 (1998).
  * Menon  _et al._ (1992) N. Menon, K. P. O’Brien, P. K. Dixon, L. Wu, S. R. Nagel, B. D. Williams, and J. P. Carini, Journal of Non-Crystalline Solids 141, 61 (1992).
  * Menon and Nagel (1995) N. Menon and S. R. Nagel, Phys. Rev. Lett. 74, 1230 (1995).
  * Wu (1991) L. Wu, Phys. Rev. B 43, 9906 (1991).
  * Ngai and Paluch (2004) K. L. Ngai and M. Paluch, The Journal of Chemical Physics 120, 857 (2004).
  * Blochowicz  _et al._ (2003) T. Blochowicz, C. Tschirwitz, S. Benkhof, and E. A. Rössler, The Journal of Chemical Physics 118, 7544 (2003).
  * Bauer  _et al._ (2013) T. Bauer, P. Lunkenheimer, S. Kastner, and A. Loidl, Phys. Rev. Lett. 110, 107603 (2013).
  * Duvvuri and Richert (2003) K. Duvvuri and R. Richert, The Journal of Chemical Physics 118, 1356 (2003).
  * Lunkenheimer  _et al._ (2005) P. Lunkenheimer, R. Wehn, U. Schneider, and A. Loidl, Phys. Rev. Lett. 95, 055702 (2005).
  * Diezemann  _et al._ (1999) G. Diezemann, U. Mohanty, and I. Oppenheim, Phys. Rev. E 59, 2067 (1999).
  * Domschke  _et al._ (2011) M. Domschke, M. Marsilius, T. Blochowicz, and T. Voigtmann, Phys. Rev. E 84, 031506 (2011).
  * Sethna  _et al._ (1991) J. P. Sethna, J. D. Shore, and M. Huang, Phys. Rev. B 44, 4943 (1991).
  * Stevenson and Wolynes (2010) J. D. Stevenson and P. G. Wolynes, Nature physics 6, 62 (2010).
  * Viot  _et al._ (2000) P. Viot, G. Tarjus, and D. Kivelson, The Journal of Chemical Physics 112, 10368 (2000).
  * Chamberlin (1999) R. V. Chamberlin, Phys. Rev. Lett. 82, 2520 (1999).
  * Dyre and Schrøder (2000) J. C. Dyre and T. B. Schrøder, Rev. Mod. Phys. 72, 873 (2000).
  * Berthier and Garrahan (2005) L. Berthier and J. P. Garrahan, The Journal of Physical Chemistry B 109, 3578 (2005).
  * Ninarello  _et al._ (2017) A. Ninarello, L. Berthier, and D. Coslovich, Phys. Rev. X 7, 021039 (2017).
  * Berthier  _et al._ (2019a) L. Berthier, P. Charbonneau, A. Ninarello, M. Ozawa, and S. Yaida, Nature communications 10, 1 (2019a).
  * Berthier  _et al._ (2017) L. Berthier, P. Charbonneau, D. Coslovich, A. Ninarello, M. Ozawa, and S. Yaida, Proceedings of the National Academy of Sciences 114, 11356 (2017).
  * Berthier  _et al._ (2019b) L. Berthier, E. Flenner, C. J. Fullerton, C. Scalliet, and M. Singh, Journal of Statistical Mechanics: Theory and Experiment 2019, 064004 (2019b).
  * Yu  _et al._ (2017) H.-B. Yu, R. Richert, and K. Samwer, Science advances 3, e1701577 (2017).
  * Illing  _et al._ (2017) B. Illing, S. Fritschi, H. Kaiser, C. L. Klix, G. Maret, and P. Keim, Proceedings of the National Academy of Sciences 114, 1856 (2017).
  * Flenner and Szamel (2015) E. Flenner and G. Szamel, Nature communications 6, 1 (2015).
  * Berthier and Ediger (2020) L. Berthier and M. D. Ediger, The Journal of Chemical Physics 153, 044501 (2020).
  * Chandler and Garrahan (2010) D. Chandler and J. P. Garrahan, Annual Review of Physical Chemistry 61, 191 (2010).
  * Keys  _et al._ (2011) A. S. Keys, L. O. Hedges, J. P. Garrahan, S. C. Glotzer, and D. Chandler, Phys. Rev. X 1, 021013 (2011).
  * Vogel and Glotzer (2004) M. Vogel and S. C. Glotzer, Phys. Rev. Lett. 92, 255901 (2004).
  * Dyre (1987) J. C. Dyre, Phys. Rev. Lett. 58, 792 (1987).
  * Bouchaud (1992) J.-P. Bouchaud, Journal de Physique I 2, 1705 (1992).
  * Rehwald  _et al._ (2010) C. Rehwald, O. Rubner, and A. Heuer, Phys. Rev. Lett. 105, 117801 (2010).
  * Rehwald and Heuer (2012) C. Rehwald and A. Heuer, Phys. Rev. E 86, 051504 (2012).
  * Arkhipov and Baessler (1994) V. I. Arkhipov and H. Baessler, The Journal of Physical Chemistry 98, 662 (1994).
  * Scalliet  _et al._ (2021) C. Scalliet, B. Guiselin, and L. Berthier, The Journal of Chemical Physics 155, 064505 (2021).
  * Xia and Wolynes (2001) X. Xia and P. G. Wolynes, Phys. Rev. Lett. 86, 5526 (2001).
  * van Meel  _et al._ (2012) J. A. van Meel, L. Filion, C. Valeriani, and D. Frenkel, The Journal of chemical physics 136, 234107 (2012).
  * Schrøder  _et al._ (2000) T. B. Schrøder, S. Sastry, J. C. Dyre, and S. C. Glotzer, The Journal of Chemical Physics 112, 9834 (2000).
  * Götze (2008) W. Götze,  _Complex dynamics of glass-forming liquids: A mode-coupling theory_ , Vol. 143 (OUP Oxford, 2008).


# FedDBL: Communication and Data Efficient Federated Deep-Broad Learning for
Histopathological Tissue Classification

Tianpeng Deng†, Yanqi Huang†, Zhenwei Shi†, Jiatai Lin, Qi Dou, Zaiyi Liu,
Xiao-jing Guo, C. L. Philip Chen, Chu Han This work was supported by Key-Area
Research and Development Program of Guangdong Province (No. 2021B0101420006),
Regional Innovation and Development Joint Fund of National Natural Science
Foundation of China (No. U22A20345), National Science Foundation for Young
Scientists of China (No. 62102103), Natural Science Foundation for
Distinguished Young Scholars of Guangdong Province (No. 2023B1515020043).
(Equal contribution: Tianpeng Deng, Yanqi Huang and Guoqiang
Han.)(Corresponding author: Xiao-jing Guo, C. L. Philip Chen and Chu
Han.)Tianpeng Deng, Jiatai Lin, Guoqiang Han and C. L. P. Chen are with the
School of Computer Science and Engineering, South China University of
Technology, Guangzhou, Guangdong, 510006, China.Yanqi Huang, Zhenwei Shi,
Zaiyi Liu and Chu Han are with the Department of Radiology, Guangdong
Provincial People’s Hospital (Guangdong Academy of Medical Sciences), Southern
Medical University, Guangzhou 510080, China and Guangdong Provincial Key
Laboratory of Artificial Intelligence in Medical Image Analysis and
Application, Guangzhou 510080, China.Qi Dou is with the Department of Computer
Science and Engineering, The Chinese University of Hong Kong, Hong Kong,
999077, China.Xiao-jing Guo is with the Department of Breast Pathology and
Lab, Tianjin Medical University Cancer Institute and Hospital, National
Clinical Research Center of Cancer, Key Laboratory of Cancer Prevention and
Therapy, Tianjin, Tianjin’s Clinical Research Center for Cancer, Key
Laboratory of Breast Cancer Prevention and Therapy, Tianjin Medical
University, Ministry of Education, Tianjin 300060, China.

###### Abstract

Histopathological tissue classification is a fundamental task in computational
pathology. Deep learning-based models have achieved superior performance but
centralized training with data centralization suffers from the privacy leakage
problem. Federated learning (FL) can safeguard privacy by keeping training
samples locally, but existing FL-based frameworks require a large number of
well-annotated training samples and numerous rounds of communication which
hinder their practicability in the real-world clinical scenario. In this
paper, we propose a universal and lightweight federated learning framework,
named Federated Deep-Broad Learning (FedDBL), to achieve superior
classification performance with limited training samples and only one-round
communication. By simply associating a pre-trained deep learning feature
extractor, a fast and lightweight broad learning inference system and a
classical federated aggregation approach, FedDBL can dramatically reduce data
dependency and improve communication efficiency. Five-fold cross-validation
demonstrates that FedDBL greatly outperforms the competitors with only one-
round communication and limited training samples, while it even achieves
comparable performance with the ones under multiple-round communications.
Furthermore, due to the lightweight design and one-round communication, FedDBL
reduces the communication burden from 4.6GB to only 276.5KB per client using
the ResNet-50 backbone at 50-round training. Since no data or deep model
sharing across different clients, the privacy issue is well-solved and the
model security is guaranteed with no model inversion attack risk. Code is
available at https://github.com/tianpeng-deng/FedDBL.

###### Index Terms:

Federated Learning, Data and Communication Efficiency, Deep-Broad Learning,
Histopathological Tissue Classification

##  I Introduction

Tissue classification [gurcan2009histopathological, FUCHS2011515], also known
as tissue phenotyping, aims to use computer algorithms to automatically
recognize different tissue types in the Whole Slide Images (WSIs). It is one
of the fundamental tasks in computational pathology [srinidhi2021deep,
wang2019weakly] which can parse the landscape of tumor microenvironment for
precise predictions of cancer diagnosis [bulten2022artificial], prognosis
[fu2020pan, pages2018international] and treatment response
[vanguri2022multimodal]. With the advancements of deep learning algorithms and
the growing number of open data [kather2019predicting, zhao2020artificial,
yang2022two], this problem has been well studied with outstanding
classification performance [hatami2021deep]. In clinical practice, however, it
still faces ethical, regulatory and legal obstacles where centralized data
collection may lead to privacy leakage, especially the RAW data.

Federated Learning (FL) [yang2019federated, le2021federated] framework
provides a promising solution to protect user privacy by only sharing the
intermediate results or the model parameters instead of the raw data, which
has been widely studied in medical image analysis [pati2022federated,
sheller2020federated]. But only very few attempts [saldanha2022swarm,
shen2022tmi, ke2021isbi] have been made in computational pathology and the
research progress still lags behind other medical image modalities
[rauniyar2022federated] due to the following two obstacles.

The first one is the data dependency problem. Since most of the existing FL
frameworks are constructed based on deep learning models. They are data-hungry
and commonly require a large amount of well-annotated samples. However,
labeling histopathological images is time-consuming, expertise-dependent and
expensive [greenwald2022whole, pati2021reducing]. When without enough training
samples, existing models may not achieve favorable performance. Another
obstacle is the communication overhead. The training procedure of traditional
FL models needs multiple cloud-client iterations to achieve global
convergence. However, deep learning models are with tens of millions of
parameters, which greatly increases the communication burden when with
multiple communication rounds. Although some recent works use self-supervised
learning SSL methods with unlabelled data to reduce the demand of labelled
data [yan2023label, kassem2022federated], they still need multiple
communication rounds to train a stable domain-specific pre-trained model. Lack
of training samples may further amplify the communication burden because deep
learning models commonly require more iterations to converge when with limited
training samples. Moreover, frequent communications may increase the chance of
being attacked, such as man-in-the-middle attacks [wang2020man].

Therefore, it is urged to construct a data-efficient and communication-
efficient FL model for histopathological tissue classification. In this paper,
we proposed a simple and effective solution for histopathological tissue
classification, which considers not only the data sharing problem, but also
the data dependency, communication efficiency, model robustness and model
inversion attack. Our proposed model Federated Deep-Broad Learning, FedDBL in
short, contains three integrated components, including a common federated
learning framework, a pre-trained deep learning (DL) backbone and a broad
learning (BL) inference system [BLS, gong2022research]. The federated learning
framework serves for decentralized training to avoid data sharing across
different medical centers or institutions. The pre-trained DL backbone can
provide stable and robust deep features when there are not enough training
labels even with domain-irrelevant pre-trained DL backbone. It can also
effectively avoid the model inversion attack since no back-propagation is
calculated for gradients [zhu2019deep]. The BL system is a lightweight
classifier with good approximation capability which can greatly shorten the
transmission time and overcome the data dependency problem. Fig. 1
comprehensively demonstrates the strengths of FedDBL compared with the
centralized learning and the conventional federated learning ways.

Figure 1: Overall comparison among centralized training, traditional DL-based
FL and our proposed FedDBL paradigms. (a) Centralized learning gathers data
from all the clients which cannot protect the patient’s privacy. (b)
Traditional DL-based FL preserves privacy by transmitting the model parameters
to the central server without sharing the raw data. However, the communication
overload highly depends on the model size and the number of communication
rounds. (c) Our proposed FedDBL not only protects privacy, but also
dramatically saves the communication burden through a super lightweight
trainable broad learning system.

Extensive experiments with five-fold cross-validation are conducted to
demonstrate the superiority of FedDBL in several aspects, including data
dependency, communication efficiency, flexibility and the practicability of
the model encryption. When with enough training data, FedDBL can mostly
outperform conventional FL strategies and achieve comparable or even better
classification performance compared with centralized learning strategy. When
reducing the training samples in the data dependency experiment, FedDBL still
maintains a high-level performance and greatly outperforms both centralized
learning and conventional FL frameworks, even with only 1% training samples.
FedDBL is also flexible to any deep learning architectures to support data-
and communication-efficient histopathological tissue classification. Another
spotlight of FedDBL is communication efficiency. Compared with the
conventional FL frameworks, FedDBL’s one-round training manner reduces the
upload workload from 4.609GB to 276.5KB (over 17,000 times faster) with
ResNet-50 backbone compared to traditional 50-round iterative training. Thanks
to the tiny model size, FedDBL is also computationally efficient in model
encryption which can further upgrade the privacy protection level. The main
contributions of this paper can be summarized as follows:

  * •

We propose a novel federated learning approach (FedDBL) for histopathological
tissue classification to preserve patients’ privacy.

  * •

To the best of our knowledge, FedDBL is the first study that considers
communication efficiency and data efficiency simultaneously which reduces the
communication overhead of each client by around 17,000 times on ResNet-50 with
extremely limited training samples (only 1%).

  * •

FedDBL is a simple, effective and easy-to-use algorithm that associates three
classical modules, including a robust pre-trained deep learning feature
extractor, a fast broad learning inference system and a simple federated
learning framework. It is highly extendable that allows to replace any module
with a more advanced one.

  * •

Extensive experiments demonstrate that FedDBL drastically relieves the
dependence on training data and reduces the communication overhead while
maintaining outstanding classification performance, which promotes its
clinical practicability.

##  II Related Works

###  II-A Histopathological Tissue Classification

High-resolution WSIs offer a wide range of tissue phenotypes where the pixel-
level annotation is time-consuming and requires a great deal of biomedical
knowledge ([srinidhi2021deep]), making patch-level histopathological tissue
classification an alternate solution for automated analysis in computer-aided
tumor diagnosis ([kather2019predicting, xue2021selective,
abdeltawab2021pyramidal]).

Due to the rapid development of computer vision, the most popular natural
image classification models can be transferred into histopathological tissue
phenotyping. However, it still suffers from the data dependency problem with a
huge annotation burden ([ayyad2021role]). Thus, various approaches have been
proposed to reduce the annotation effort. [han2022multi] proposed a multi-
layer pseudo-supervision approach with a progressive dropout attention
mechanism to convert patch-level labels into pseudo-pixel-level labels. An
extra classification gate mechanism was presented which reduced the false-
positive rate for non-predominant category classification and improved the
segmentation performance in return. [xue2021selective] utilized a generative
adversarial network (GAN) to generate pseudo samples to expand the training
data. [dolezal2022uncertainty] cropped WSIs into tiles for training the
uncertainty quantization model and solved the problem of domain shift in
external validation data. In order to get rid of lacking image annotations,
[wang2022transformer] employed unsupervised contrastive learning to obtain a
robust initialized model with moderate feature representation of the
histopathological feature space, with no annotation burden. Our previous study
([lin2022pdbl]) introduced pyramidal deep-broad learning (PDBL) as a pluggable
module for any CNN backbone to further improve histopathological tissue
classification performance.

Besides that, another unexplored challenge is the patient privacy issue. Only
a few attempts ([saldanha2022swarm, saldanha2022direct]) have been made in
federated learning for computational pathology, which will be discussed in the
following subsection. And to the best of our knowledge, we are the first study
to consider privacy protection in histopathological tissue classification.

###  II-B Federated Learning

####  II-B1 Federated Learning in Medical Image Analysis

Because of the ethical issue, federated learning (FL) has been widely adopted
in medical applications to preserve the patients’ privacy ([pati2022federated,
warnat2021swarm, sheller2020federated]). In medical imaging, FL has witnessed
a boost in interest ([kaissis2020secure]), such as MRI reconstruction
([guo2021multi, li2020multi]), CT lesion segmentation ([yang2021federated])
and etc. In the COVID-19 pandemic, COVID-19-related applications with data
from different medical centers or even from different countries become the
most urgent demand in the real-world clinical scenario while FL greatly
advances the diagnostic performance ([bai2021advancing]). [dayan2021federated]
used 20 institutes’ data across the global for predicting the future oxygen
requirements of symptomatic patients suffering from COVID-19.
[dou2021federated] proposed a federated model to detect COVID-19 lung
abnormalities with good generalization capability on unseen multinational
datasets.

####  II-B2 Federated Learning in Computational Histopathology

In histopathological images, a swarm learning architecture with blockchain
protocols has been proposed to predict the mutational status
([saldanha2022swarm]). However, compared with other medical imaging
modalities, there are few studies ([saldanha2022direct]) that adopt federated
learning in histopathological images for the following reasons. First, the
digitalization of pathology is unpopular. Pathological diagnosis still relies
on observing specimens under a microscope. Second, image annotation is also an
obstacle for the histopathological image process since only pathologists are
capable to label WSIs which greatly increases the difficulties of acquiring
well-annotated data. Third, due to the gigapixel resolution of WSIs, the size
of the deep learning model is generally large, which increases the
communication burden in networking.

There are technical solutions in FL to the high communication overhead
problem, such as compressing the model size ([reisizadeh2020fedpaq,
jhunjhunwala2021adaptive]). [reisizadeh2020fedpaq] proposed FedPAQ to reduce
the interactive overhead of FL by compressing the model with lower bit-
precision and [jhunjhunwala2021adaptive] proposed an adaptive quantization
strategy to achieve communication efficiency.

However, the underlying assumption of existing studies is that there should be
enough samples for model training where they may not be able to take into
account both communication efficiency and limited data issue
([kamp2021federated, zhang2023two]). In this study, we fully consider the
specialty of histopathological images, the difficulties of data labeling and
the communication efficiency in the real-world clinical scenario, which has
never been discussed in decentralized computational pathology.

##  III Methodology

In this section, we introduce our framework Federated Deep-Broad Learning
(FedDBL). This framework is designed for privacy-preserving tissue
classification with limited training samples and extremely low communication
overhead. In the following subsections, we first describe the intuitive
thinking and problem setting in Section III-A. The overall framework and the
methodology of FedDBL are shown in Section III-B. Finally, we demonstrate the
implementation details in Section III-C.

###  III-A Problem Setting

As a classical upstream task in computational pathology, existing tissue
classification approaches have achieved outstanding performance under an ideal
condition with enough training samples by centralized learning. However, they
might face the following obstacles in the real-world clinical scenario.

Annotation burden: Collecting enough well-labeled training samples is
expensive and time-consuming because it requires labelers with medical
background.

Privacy preservation: The raw data should not be shared across different
medical institutions (or clients) to preserve the patient’s privacy.
Transmitting raw data may break the principle of medical ethics.

Communication cost: The communication overhead has always been a challenge in
federated learning models affected by many compound factors, such as the model
size, the communication rounds, the model convergence speed, the network
bandwidth and etc.

To resolve the aforementioned challenges, we propose a simple and effective
FL-based framework, demonstrated in Fig. 2. First, we abandon conventional
end-to-end training manner since limited training samples may harm the
robustness of the deep learning model and decrease the convergence speed.
Therefore, we separate feature extraction and inference for local training in
each client. A pre-trained deep feature extractor (CNN backbone) is introduced
to avoid the feature extractor being affected by the training sample bias from
different clients in order to guarantee the robustness of extracted features.
Then an independent broad learning inference system [BLS, lin2022pdbl] serves
for fast inference. Finally, we apply a classical weighted averaging as in
FedAvg [mcmahan2017communication], to fuse the broad learning inference
systems from all the clients.

Figure 2: The overall architecture of FedDBL with three modules, deep feature
extraction module, broad inference module and federated decentralized module.
(a) Deep feature extraction module serves for extracting multi-scale deep-
broad features from low level to high level by a pre-trained DL backbone.
Features of all the patches are stored in a local feature bank. (b) Broad
inference module is introduced for fast inference by a broad learning system.
(c) Federated decentralized module applies a classical federated average
approach to aggregate the broad learning weights from different clients.

###  III-B FedDBL Architecture and Formulation

As shown in Fig. 2, FedDBL consists of three modules, deep feature extraction
module (DL-module), broad inference module (BL-module) and federated
decentralized module (Fed-module). DL-module together with BL-module serves
for local training on the client side. Fed-module is executed on the server
side. Algorithm 1 provides the details of the entire FedDBL pipeline.

Let
$\mathcal{D}_{1},\mathcal{D}_{2},\cdots,\mathcal{D}_{k},\cdots,\mathcal{D}_{K}$
denote the local training sets from $K$ clients with the dataset size of
$n_{k}$ for each client $\mathcal{D}_{k}$. The total number of training
samples is denoted as $N=\sum_{k=1}^{K}{n_{k}}$. For each sample $X$ with
ground truth $Y$ in $\mathcal{D}$, DL-module with pre-trained parameters
$\Theta$ extracts the features and stores them in the local feature bank
$\mathbf{B}$. Then BL-module calculates the weights $W_{client}$ of broad
learning system. By the federated aggregation approach, we can obtain the
global weight $W_{global}$. The workflows of the server and the clients are
demonstrated in Algorithm 1 and Algorithm 2, respectively.

Input : A set of $K$ clients

Output : A global model $W_{global}$

1 Prepare pre-trained DL backbone parameters $\Theta$

2Initialize BL system setting

3for  _each client $k$ _in parallel__ do

4
$W_{client}^{k}\leftarrow$ClientExecution$\left(\Theta,\mathcal{D}_{k}\right)$

5 end for

6$W_{global}\leftarrow$ Fed-module$(W_{client}^{1},\cdots,W_{client}^{K})$

return $W_{global}$

Algorithm 1 FedDBL framework (Server Execution)

Input : Pre-trained DL backbone $\Theta$, training set $\mathcal{D}$ with $n$
training samples

Output : Deep-broad learning model $W_{client}$.

1

/* DL-module */

2 for  _training sample $X$ in $\mathcal{D}$_ do

3 for  _$s$ -th stage in $\Theta$_ do

$f_{X}^{s}\leftarrow\Theta^{s}\left(X\right)$ // Feature extraction

$\mathbf{e}_{X}^{s}=\frac{1}{H_{X}^{s}\times
W_{X}^{s}}\sum_{i=1}^{H_{X}^{s}}\sum_{j=1}^{W_{X}^{s}}f_{X}^{s}(i,j)$ //
Adaptive global average pooling

4

5 end for

$\mathbf{b}_{X}=\mathbf{e}_{X}^{1}\parallel\mathbf{e}_{X}^{2}\parallel\cdots\parallel\mathbf{e}_{X}^{S}$
// Concatenation

6

7 end for

8Obtain $\left\\{\mathbf{b}_{i}|i=1,2,\cdots,n\right\\}$

9

$\mathbf{B}\leftarrow\sigma\left(\left[\mathbf{b_{1}}^{\mathbf{T}},\mathbf{b_{2}}^{\mathbf{T}},\dots,\mathbf{b}_{n}^{\mathbf{T}}\right]\right)^{\mathbf{T}}$
// Normalization transformation

10

/* BL-module */

11 Initialize BL system setting defined by central server

12

$\mathbf{B}^{+}\leftarrow\lim_{\lambda\to
0}\left(\mathbf{B}\mathbf{B}^{\mathbf{T}}+\lambda
E\right)^{-1}\mathbf{B}^{\mathbf{T}}$ // Solve Pseudo-inverse

13

$W_{client}\leftarrow\mathbf{B}^{+}Y$ // Calculate BL model weight

14

15return $W_{client}$

Algorithm 2 FedDBL framework (Client Execution)

####  III-B1 Deep Feature Extraction Module

A large number of samples and repeated backpropagation are required in
standard DL training to achieve a good feature representation ability. When
suffering from the insufficient data problem, the model training procedure
might be unstable which leads to poor feature representation and model
overfitting. Our previous study [lin2022pdbl] reveals that directly adopting a
stable pre-trained model for feature extraction is more favorable to the model
performance than training the model with limited samples, even the pre-trained
model was trained by an irrelevant image domain (ImageNet111https://image-
net.org/). Inspired by this idea, we use a pre-trained CNN model with no
further training to extract the deep features. Notice that, the selection of
the pre-trained models is flexible, and can be from any image domain. We have
conducted an experiment to justify the flexibility in Section IV. Another
advantage of using pre-trained models is to avoid model inverse attacks since
the training samples are all unseen. To enrich the feature representation, we
extract multi-stage features from low-level to high-level, details as below.

As illustrated in DL-module of Algorithm 2, each client $k$
$(k\in\left[1,\cdots,K\right])$ downloads the pre-trained DL backbone as
feature extractor $\Theta$ and extracts multi-stage deep features
$\mathbf{b}_{X}$ of training sample $X$ locally (we omit $k$ for simplicity),
where $\mathbf{b}_{X}$ consists of multiple stages’ features
$\Theta^{s}(X)(s\in\left[1,\cdots,S\right])$. The features of the entire
dataset $\mathcal{D}_{k}$ are stored in the local feature bank $\mathbf{B}$.
Then the local feature bank $\mathbf{B}$ will be passed to broad inference
module. Since neither the training data nor the feature bank is shared across
different clients, there is no privacy leakage risk for the RAW data in deep
feature extraction module.

####  III-B2 Broad Inference Module

With the local feature bank $\mathbf{B}$, each client $k$ can conduct a local
BL system [BLS] through BL-module (Algorithm 2) for fast inference. By solving
the Eq. (1) optimization problem, an optimal BL model $W_{client}$ can be
obtained rapidly through the pseudo-inverse method (Eq. (2)).

$W_{client}=\underset{W_{init}}{{\arg\min}}\left\|\mathbf{B}W_{init}-Y\right\|_{2}^{2}+\gamma\left\|W_{init}\right\|_{2}^{2}$
(1) $W_{client}=\mathbf{B}^{+}Y=\lim_{\lambda\to
0}\left(\mathbf{B}\mathbf{B}^{\mathbf{T}}+\lambda
E\right)^{-1}\mathbf{B}^{\mathbf{T}}Y$ (2)

where $Y$ represents the ground-truth label matrix, $\mathbf{B}$ is local
feature bank in the form of matrix. $W_{init}$ is the initialized broad
learning weights. $E$ is the identity matrix, $\lambda$ is a constant
parameter and $\gamma$ is the regularization parameter. The pseudo-inverse
method of solving BL model considerably reduces the computational burden while
achieving high communication efficiency. For the inference process, the
predicted results can be calculated by $Y_{test}=\mathbf{B}_{test}W_{client}$
after extracting test samples’ deep features with the largest probabilistic
value.

Thanks to the lightweight broad learning model $W_{client}$, the communication
efficiency is drastically improved compared with the conventional DL-based FL
frameworks.

####  III-B3 Federated Decentralized Module

In this module, we conduct a federated learning framework for decentralized
learning. Given the broad learning model $W_{client}^{k}$ of each client $k$,
we first upload the models from all the clients to the central server. And
then general federated aggregation methods can be applied to aggregate them.
Here, we use the most common weighted averaging way for model aggregation as
adopted in FedAvg [mcmahan2017communication], FedProx [li2020federated] and
FedPAQ [reisizadeh2020fedpaq].

$W_{global}=\sum_{k=1}^{K}\frac{n_{k}}{N}W_{client}^{k}$ (3)

where $W_{global}$ is the global model from the server, $n_{k}$ is the number
of training samples in client $k$ and $N$ is the total number of training
samples. A larger training dataset will contribute more to the global model.
Since we only share the broad learning model for once, the communication
efficiency and the patient’s privacy are guaranteed.

###  III-C Implementation Details

All of our experiments are implemented in Pytorch on a workstation with an
NVIDIA RTX 3090 and the i9-11900K CPU with 16 cores. We use the cross-entropy
loss for the baseline centralized training with batch size $20$. The SGD
optimizer is set as follows: the learning rate is $1e^{-3}$, the weight decay
is $1e^{-4}$ and the momentum is 0.9. The patches are $224\times 224$ under
$20\times$ WSIs. Different client numbers are used depending on the datasets.

We adopt three well-known federated aggregation methods, FedAvg
[mcmahan2017communication], FedProx [li2020federated] and FedPAQ
[reisizadeh2020fedpaq], for comparison. And the centralized model is trained
as the baseline. FedProx has the parameter $\mu$ to adjust the effect of the
proximal term on the loss function. Here we set $\mu$ as $1$ which has better
performance.

##  IV Experiments

In this section, we present the details of the datasets and conduct various
experiments to demonstrate the performance and efficiency of the proposed
_FedDBL_. Section IV-A shows two open datasets and the experimental settings
in the federated learning framework. In Section LABEL:sub:one-round, we
compare _FedDBL_ with centralized learning baselines, conventional federated
learning baselines and one-round federated learning baselines. The
effectiveness is comprehensively discussed in Section LABEL:sub:multi-round.
We use Matthews Correlation Coefficient (MCC), Accuracy and F1-score as the
evaluation metrics in all the experiments.

###  IV-A Datasets and Experimental Settings

TABLE I: Statistics of MC-CRC. #1 denotes TCGA, #2 denotes Kather, #3 denotes Guangdong Provincial People’s Hospital and #4 denotes Yunnan Cancer Hospital. |  A DI |  B ACK |  D EB |  L YM |  M UC |  M US |  N ORM |  S TR |  T UM
---|---|---|---|---|---|---|---|---|---