# Monte Carlo simulations of glass-forming liquids beyond Metropolis

Ludovic Berthier  Laboratoire Charles Coulomb (L2C), Université de Montpellier
& CNRS (UMR 5221), 34095 Montpellier, France  Gulliver, UMR CNRS 7083, ESPCI
Paris, PSL Research University, 75005 Paris, France  Federico Ghimenti
Laboratoire Matière et Systèmes Complexes (MSC), Université Paris Cité & CNRS
(UMR 7057), 75013 Paris, France  Frédéric van Wijland  Laboratoire Matière et
Systèmes Complexes (MSC), Université Paris Cité & CNRS (UMR 7057), 75013
Paris, France

###### Abstract

Monte Carlo simulations are widely employed to measure the physical properties
of glass-forming liquids in thermal equilibrium. Combined with local Monte
Carlo moves, the Metropolis algorithm can also be used to simulate the
relaxation dynamics, thus offering an efficient alternative to molecular
dynamics. Monte Carlo simulations are however more versatile, because
carefully designed Monte Carlo algorithms can more efficiently sample the
rugged free energy landscape characteristic of glassy systems. After a brief
overview of Monte Carlo studies of glass-formers, we define and implement a
series of Monte Carlo algorithms in a three-dimensional model of polydisperse
hard spheres. We show that the standard local Metropolis algorithm is the
slowest, and that implementing collective moves or breaking detailed balance
enhances the efficiency of the Monte Carlo sampling. We use time correlation
functions to provide a microscopic interpretation of these observations.
Seventy years after its invention, the Monte Carlo method remains the most
efficient and versatile tool to compute low-temperatures properties in
supercooled liquids.

##  I Introduction: Monte Carlo simulations of glass-formers

Understanding how glasses are formed upon slowly supercooling liquids and
predicting their physical properties remain broad and important research
topics [1, 2, 3]. From a fundamental viewpoint, one would like to obtain a
description of these emerging properties starting from the sole knowledge of
how microscopic degrees of freedom move and interact [4, 3, 2, 5, 6]. This
statistical mechanics program is still incomplete, because fluctuations are
typically hard to control analytically in disordered systems in finite
dimensions [6]. Computer simulations thus represent a central tool, since
finite dimensional models can be simulated with an exquisite control of the
microscopic interactions leading to glass formation [7].

Although Monte Carlo [8] and Molecular Dynamics [9] were invented just a few
years apart, computer studies of the glass transition in supercooled liquids
have been largely dominated by Molecular Dynamics [10]. This can be justified
by the fact that the glass transition is primarily a kinetic phenomenon where
the dramatic slowing down of the dynamics upon cooling eventually leads to the
emerging solidity of the amorphous solid [1]. Instead, the Monte Carlo
approach was devised to provide accurate numerical estimates of
configurational averages that involve an equilibrium sampling of the
configuration space [8, 11, 12, 13]. Thus, Monte Carlo simulations of particle
systems were used to study the static structure and thermodynamic properties
of liquids, thereby usefully accompanying the development of liquid state
theories [14]. In supercooled liquids, however, the structure evolves very
little as the system approaches the glass transition, at least at the level of
two-point density correlation functions [4, 5, 15, 16].

The dominance of molecular dynamics studies receded over the last two decades
for at least two reasons. First, it is now understood that there is more to
the structure and thermodynamics of supercooled liquids than two-body
functions [6, 15]. More complex correlation functions have been shown to be
relevant [17, 18, 19], while configurational entropy [20, 21, 22] and
constrained free energies [23, 24, 25] are also useful indicators of the glass
transition. These call for the development of efficient Monte Carlo techniques
to probe configurational averages in equilibrium conditions. Second, it has
been realized that details of the microscopic dynamics do not affect the
glassy dynamics at long times [26]. In particular, local Monte Carlo moves
lead to dynamic relaxation pathways that are equivalent to the ones explored
by the Molecular Dynamics [27, 28, 29]. Thus, local Monte Carlo dynamics has
become an efficient alternative to Molecular Dynamics to sample not only
static but also dynamic properties of supercooled liquids in the canonical
ensemble [7].

A direct consequence of the dynamic equivalence between local Monte Carlo and
Molecular Dynamics is that sampling the equilibrium distribution becomes very
difficult in both approaches when temperature decreases because structural
relaxation times, which are good indicators of sampling efficiency and ergodic
exploration, grow rapidly as the glass transition is approached. For the
reasons given above, algorithms that can speedup equilibration and equilibrium
sampling in deeply supercooled liquids are very much needed and sought after.
Within the Monte Carlo approach, the price to pay is that the dynamic
information is often lost as algorithms take unphysical shortcuts in
configuration space to achieve a faster sampling [11, 12, 13]. This is maybe
not as bad as it sounds, since samples that are efficiently equilibrated at
very low temperatures can still be used as initial conditions for Molecular
Dynamics or local Monte Carlo simulations [30].

Walking in the footsteps of the Metropolis algorithm, several directions can
already be followed to improve the sampling efficiency of Monte Carlo
simulations of supercooled liquids [31]. One could invent more complicated
particle moves, going beyond the small single particle translational move. A
cluster algorithm was designed for dense systems of hard disks that proposes
moves involving a large number of particles [32]. Another example is the swap
Monte Carlo algorithm which exchanges the positions of randomly chosen pairs
of particles. It was first proposed in the context of simulations of multi-
component solids [33], the swap Monte Carlo was later introduced in
supercooled liquids [34], and its efficiency has been improved and
demonstrated across a broad range of systems and spatial dimensions [35, 36,
37, 38, 39]. Another generic strategy within the Metropolis realm is to add
additional degrees of freedom to the system, such as temperature, to explore a
larger phase space. Parallel tempering is the most studied example of this
family [40], where a range of temperatures are simultaneously simulated using
replicas of the system, with infrequent temperature swaps proposed between
replicas [41, 42].

The versatility of Monte Carlo simulations becomes obvious when methods
reaching beyond the realm of the Metropolis acceptance rule are defined. As
noted in many textbooks, the Metropolis strategy that ensures detailed balance
is sufficient to reach thermal equilibrium but is not necessary, as global
balance (with proof of ergodicity) is sufficient [43]. In practice, defining
Monte Carlo algorithms that break detailed balance and ensure Boltzmann
sampling is not easy, but several such algorithms have been proposed recently
in different contexts [44, 45]. In particular, lifting Monte Carlo algorithms
introduce additional degrees of freedom to ensure global balance while
breaking detailed balance [46, 47, 48, 49, 50]. For dense particle systems,
the event-chain Monte Carlo algorithm was first proposed for hard disks [44,
51], and later generalized to soft potentials [52]. We have recently
introduced a nonequilibrium algorithm [53] where lifting is applied to both
translational and diameter degrees of freedom, thus generalizing both the
event-chain and swap Monte Carlo algorithms.

This brief overview of Monte Carlo approaches applied to supercooled liquids
suggests that more efficient algorithms are still crucially needed. Here, we
build on our previous work on two-dimensional hard disks [53] and assess
whether the proposed nonequilibrium algorithms remain efficient in three
dimensions. To this end, we define and implement a similar family of Monte
Carlo algorithms for a polydisperse model of hard spheres known to display
glassy dynamics [54]. We then compare their performances to reach and sample
equilibrium as well as for nonequilibrium compressions towards the jamming
transition. Overall, our results show that the performances of nonequilibrium
swap algorithms remain very good in three dimensions.

The manuscript is organized as follows. In Sec. II we define the studied model
of hard spheres. In Sec. III, we introduce the various Monte Carlo algorithms
we analyze. In Sec. IV we confirm that all algorithms properly sample the same
equilibrium Boltzmann distribution, by following the equation of state. In
Sec. V we compare the efficiency of all Monte Carlo algorithms. In Sec. VI we
study various time correlation functions to provide microscopic insights into
the performances of the algorithms. Finally, in Sec. VII we use these Monte
Carlo algorithms to prepare very dense jammed packings of spheres. We conclude
the paper in Sec. VIII.

##  II Polydisperse hard sphere model

We simulate a polydisperse mixture of $N=1000$ hard spheres in three
dimensions. The pair interaction is $V(r)=0$ when the distance between
particle pairs $(i,j)$ satisfies $r_{ij}>\sigma_{ij}$ with
$\sigma_{ij}=(\sigma_{i}+\sigma_{j})/2$, and $V(r)=\infty$ otherwise.

The diameters $\\{\sigma_{i},i=1\cdots N\\}$ of the particles are drawn from a
power law distribution [54, 55], $\pi(\sigma)\propto\sigma^{-3}$, with
boundaries $\sigma_{\rm min}$ and $\sigma_{\rm max}$ chosen so that the
polydispersity
$\Delta\equiv\frac{\sqrt{\overline{\sigma^{2}}-\overline{\sigma}^{2}}}{\overline{\sigma}}$
is $\Delta\approx 23\%$, with the overline denoting an average over
$\pi(\sigma)$. We study the dynamics of the system as a function of the
packing fraction $\phi$ defined as

$\phi=\frac{\pi N\overline{\sigma^{3}}}{6V},$ (1)

where $V$ is the volume of the system. Periodic boundary conditions are used.
In Secs. V and VI we use the $NVT$ ensemble, where the volume $V$ is fixed. We
move to the $NPT$ ensemble to perform compressions to measure the equation of
state in Sec. IV and to arrive to the jamming point in Sec. VII. In that case,
the volume of the simulation box is allowed to vary to maintain the pressure
$P$ constant.

We use reduced units to present our results: the packing fraction is non-
dimensional, as well as the compressibility factor $Z=P/(\rho k_{B}T)$, with
$\rho=N/V$ the number density and $k_{B}$ the Boltzmann constant. The average
particle diameter $\overline{\sigma}$ sets the unit length. For each algorithm
presented in Sec. III, we define an elementary timescale $t_{\rm move}$ such
that $\tau_{0}=Nt_{\rm move}$ typically represents a Monte Carlo sweep across
the system. Therefore we use $\tau_{0}$ as the time unit when presenting time-
dependent quantities and relaxation timescales.

##  III Monte Carlo algorithms: Metropolis and beyond

In this section we introduce the Monte Carlo algorithms implemented to
simulate the hard sphere model presented in Sec. II. Following [53], we
carefully introduce a unit of time $t_{\text{move}}$ for each algorithm, from
which a quantitative comparison of their performances can be achieved. In Fig.
1 we sketch the various types of Monte Carlo moves described in this section.

Figure 1: Elementary Monte Carlo moves. (a) In the Metropolis Monte Carlo
(MMC) algorithm, randomly chosen particles sequentially perform small random
displacements, which are accepted according to the Metropolis rule. (b) In the
Swap algorithm, randomly chosen pairs of particles exchange their positions
with an acceptance probability given by the Metropolis rule. (c) In the Event
Chain Monte Carlo (ECMC) , directed chains of particle displacements are
implemented. This violates detailed balance, but preserves the Boltzmann
distribution. (d) In the collective Swap (cSwap) algorithm, collective swap
moves are implemented. This violates detailed balance but preserves the
Boltzmann distribution.

###  III.1 Metropolis algorithms: Local moves (MMC) and Swap algorithm

In the simplest version of the Metropolis Monte Carlo [8] (MMC) algorithm,
during an elementary time $t_{\text{move}}$ a particle is randomly and
uniformly selected, and a small random displacement is proposed, see Fig.
1(a). The displacement vector is uniformly drawn from a cube of linear size
$\delta$ centered around the particle. The proposed displacement is accepted
if it does not generate any overlap with the neighboring particles, which is
the consequence of using the Metropolis rule for the hard core potential
$V(r)$ defined in Sec. II. Therefore, the MMC algorithm is local and satisfies
detailed balance, thus leading to an equilibrium sampling of the Boltzmann
distribution. For the simulations presented in this work we choose
$\delta=0.115\overline{\sigma}$, obtained by optimizing the relaxation time of
the system at intermediate packing fraction [27].

In the swap Monte Carlo algorithm [34, 35], which we simply call ‘Swap’, pairs
of particles can perform swap moves. A single swap move takes place during a
time interval $t_{\text{move}}$. During a swap move, a pair of particles is
randomly and uniformly selected, and their positions are exchanged if this
generates no overlap among the spheres, see Fig. 1(b). This acceptance rule is
again the consequence of the Metropolis rule for the hard core potential
$V(r)$ in the case of a proposed particle swap. The sketch in Fig. 1(b) shows
that it is equivalent to considering that the two swapped particles exchange
their diameters at fixed position, or exchange their positions at fixed
diameters. While the latter view is more natural, we adopt the former because
it allows a simpler definition of single-particle time correlation functions
[35], as studied in Sec. VI. In the Swap algorithm, we alternate randomly
between sets of $N$ swap moves and sets of $N$ MMC moves. A set of $N$ swap
moves is performed with probability $p_{\text{Swap}}$, while a set of $N$ MMC
moves is performed with probability $1-p_{\text{Swap}}$. Following previous
work [54, 35], we take $p_{\text{Swap}}=0.2$, which optimizes the relaxation
time for the decay of the self-intermediate scattering function, see Eq. (3),
and definitions below.

###  III.2 Beyond Metropolis: Event-chain (ECMC) and collective swap (cSwap)

Both the local Metropolis and the Swap algorithms in Sec. III.1 obey detailed
balance, which is a sufficient condition to achieve a proper sampling of the
Boltzmann distribution in the steady state. Both for translational moves and
diameter swaps, the design of algorithms that break detailed balance while
preserving the Boltzmann distribution in their stationary state is possible.
Such strategy goes beyond the seminal Metropolis approach [8]. We now present
two possible dynamics for the hard sphere model.

In the original ECMC algorithm [44], the configuration space of the particles’
positions $\text{\bf r}^{N}$ is lifted to an extended space $\\{\text{\bf
r}^{N},i,\text{\bf v}\\}$, which now includes the label $i$ of an active
particle and a direction of motion v. During a time $t_{\text{move}}$, the
active particle $i$ moves along the direction v, until it collides with a
neighbor $j$. When this happens, a lifting event takes place, and the activity
label changes from $i$ to $j$. In this way, after a time $nt_{\text{move}}$, a
directed chain of $n$ particles has been displaced along the direction v, see
Fig. 1(c). When the individual displacements of the particles involved in a
chain add up to a total length $\ell$, the activity label is uniformly
resampled among the $N$ particles, while the self-propulsion direction v is
uniformly resampled in the set $\\{\text{\bf e}_{x},\text{\bf e}_{y},\text{\bf
e}_{z}\\}$ in order to restart another chain. Contrary to Metropolis
algorithms, there is no need to satisfy any acceptance rule, which renders the
algorithm rejection-free. Also, while detailed balance is broken (chain moves
cannot be reverted), the Boltzmann distribution is the solution of the global
balance condition [44]. The efficiency of ECMC is roughly independent on the
chain length, as long as $\ell$ is of the order of a few particle diameters
and smaller than the linear size of the system. We choose
$\ell=10\overline{\sigma}$.

It is of course possible to combine ECMC with Swap to define the SwapECMC
algorithm [53]. To this end, we choose between an ECMC move which builds a
chain of length $\ell$ and $N$ swap moves with a probability
$p_{\text{SwapECMC}}$ and $(1-p_{\text{SwapECMC}})$, respectively. We take
$p_{\text{SwapECMC}}=0.2$.

In the cSwap algorithm, the lifting framework used in ECMC is implemented in
diameter space, rather than in position space [53]. A one dimensional array
storing the label of the particles in order of increasing diameters is created
at the beginning of a simulation. We define the operators
$\mathcal{L},\mathcal{R}$, that, for a given particle $i$, select respectively
the particle to the left or to the right of $i$ along the ordered array,
applying periodic boundary conditions. A configuration
$\mathcal{C}\equiv\\{\text{\bf r}^{N},\sigma^{N},i\\}$ of the system includes
all positions $\text{\bf r}^{N}$, diameters $\sigma^{N}$, and an activity
label $i$, which denotes an active particle. During a time $t_{\text{move}}$,
the largest possible diameter $\sigma_{j}$ (corresponding to particle $j$)
that can be reached by particle $i$ without generating overlaps is identified
and particle $i$ is inflated: $\sigma_{i}\leftarrow\sigma_{j}$. The $n_{s}$
non-active particles with a diameter comprised between $\sigma_{i}$ and
$\sigma_{j}$ (including $j$), perform a small deflation, moving one step to
the left in the ordered array of particle diameters:
$\sigma_{\mathcal{L}j}\leftarrow\sigma_{\mathcal{L}^{2}j},\ldots,\sigma_{\mathcal{L}^{n_{s}-1}j}\leftarrow\sigma_{\mathcal{L}^{ns}j}$.
Finally, the activity label switches from $i$ to $\mathcal{L}i$. To warrant
ergodicity, the activity label is resampled uniformly among all the particles
with probability $1-\frac{1}{N_{\text{cSwap}}}$. Overall, this results in a
collective swap move, as sketched in Fig. 1(d). The algorithm is rejection-
free, it does not invoke the Metropolis acceptance rule, and it breaks
detailed balance. In a simulation, a set of $N_{\text{cSwap}}$ moves is
performed with probability $p_{\text{cSwap}}$, while with probability
$1-p_{\text{cSwap}}$ we perform a set of $N$ local translational Metropolis
moves. We fix the probability $p_{\text{cSwap}}=0.2$, and we optimize the
number of Swap moves $N_{\text{cSwap}}$ at a density $\phi=0.624$, where
equilibrated samples are obtained with a modest effort using the cSwap
algorithm. We thus set $N_{\text{cSwap}}=500$.

We recently showed [53] that for the hard sphere interaction $V(r)$, the
stationary distribution of the cSwap algorithm, once marginalized over the
lifted degrees of freedom, is the Boltzmann distribution. By construction, the
particle size distribution is preserved by an elementary cSwap move. In
addition, the collective swap moves are performed with probability
$p_{\text{cSwap}}$, while elementary translational moves as in MMC are
performed with probability $(1-p_{\text{cSwap}})$.

Finally, ECMC and cSwap can be combined in the cSwapECMC algorithm [53]. In
this algorithm, we perform $N_{\text{cSwap}}$ cSwap moves with probability
$p_{\text{cSwapECMC}}$, and we otherwise displace a chain of length $\ell$
using the ECMC procedure with probability $1-p_{\text{cSwapECMC}}$. We set
$p_{\text{cSwapECMC}}=0.2$.

##  IV Equilibrium sampling: Equation of state

All the algorithms proposed in Sec. III can be shown to lead to a proper
sampling of the same Boltzmann distribution in the steady state. Thus,
equilibrium sampling is preserved even if the algorithms that break detailed
balance are, technically, evolving out of equilibrium. To numerically
demonstrate that this is the case, we perform slow compressions of the hard
sphere system using MMC, ECMC, Swap, and cSwap starting from a relatively
dilute fluid state at $\phi=0.53$. To slowly compress the system we introduce
Monte Carlo moves of the volume [13].

In a volume move, a volume change $\delta V$ is sampled uniformly in an
interval of width $\delta V_{\text{max}}$ centered around $0$. The new
configuration proposed for the system is then a box of volume $V+\delta V$,
with particle positions rescaled by a factor $s=\left(\frac{V+\delta
V}{V}\right)^{1/3}$. This configuration is accepted with a Metropolis
acceptance rate, tuned to ensure sampling of the $NPT$ ensemble. The new
configuration is accepted with probability $\text{e}^{-\beta P\delta V+N\log
s}$ if no overlap among the spheres are generated by the volume change.

Figure 2: Equation of state $Z=Z(\phi)$ obtained from four different
algorithms during slow compressions. The black dashed line is the Carnahan-
Starling formula $Z_{\rm CS}$ from Eq. (2).

Volume moves are introduced in a different way for each algorithm. For MMC,
Swap and cSwap, we replace with probability $1-\frac{1}{N}$ a single
translational move with a volume move. This choice is justified by the
necessity to prevent the algorithm from becoming periodic [13]. In ECMC, the
presence of a finite chain length ensures that the algorithm is unbiased [52].
The volume moves are thus implemented by replacing the displacement of a chain
of length $\ell$ with a volume move with probability $1/2$.

We use a slow compression rate, $\dot{P}=10^{-6}$, for all algorithms. This is
small enough that the system can easily reach the Boltzmann distribution in
the fluid state. However, when the packing fraction becomes large, glassy
dynamics emerges and for each algorithm there will exist a packing fraction
above which thermal equilibrium can no longer be achieved.

In Fig. 2 we present the evolution of the compressibility factor $Z$ with the
packing fraction $\phi$ for four algorithms. As expected, the algorithms all
produce the same equation of state $Z=Z(\phi)$ in the stationary fluid state.
The numerical results are in good agreement with the extension to polydisperse
systems of the empirical Carnahan-Starling formula [56]:

$Z_{\text{CS}}=\frac{1}{1-\phi}+\frac{3\phi}{(1-\phi)^{2}}\frac{\overline{\sigma}\overline{\sigma^{2}}}{\overline{\sigma^{3}}}+\frac{\phi^{2}(3-\phi)}{(1-\phi)^{3}}\frac{\overline{\sigma^{2}}^{3}}{\overline{\sigma^{3}}^{2}},$
(2)

which requires moments of the particle size distribution as input. The mutual
agreement between algorithms confirms that they all probe the same stationary
distribution.

In addition, for each algorithm we observe a packing fraction above which
clear deviations from the stationary equation of state are observed. This
happens near $\phi\approx 0.60$ for MMC and ECMC, and near $\phi\approx 0.64$
for Swap and cSwap. At these packing fractions, the system fails to relax
within the compression timescale and the measured equation of state differs
from the equilibrium result. Consistently with previous results [54], the
introduction of diameter moves considerably extends the stationary regime.
There are subtle physical differences between MMC and ECMC and Swap and cSwap
that go beyond the mere reduced convergence time, but further explorations are
required to reveal them, as we do in the next Section.

##  V Relative efficiency of Monte Carlo algorithms

Figure 3: Evolution of the relaxation time $\tau_{\alpha}$ for MMC, ECMC,
Swap, SwapECMC, cSwap, cSwapECMC algorithms as a function of the packing
fraction $\phi$.

To study the relative efficiency of the various Monte Carlo algorithms, we
need to define a relevant timescale that captures the rate at which the
configuration space is explored. To this end, we define the self-intermediate
scattering function $F_{s}$, defined as

$F_{s}(\text{\bf
k},t)\equiv\frac{1}{N}\sum_{i}\left\langle\text{e}^{-i\text{\bf
k}\cdot\Delta\text{\bf r}_{i}(t)}\right\rangle,$ (3)

with $\Delta\text{\bf r}_{i}(t)\equiv\text{\bf r}_{i}(t)-\text{\bf
r}_{i}(0)-\Delta\text{\bf r}_{\text{com}}(t)$, and $\Delta\text{\bf
r}_{\text{com}}$ is the displacement of the center of mass of the system over
the time interval $[0,t]$, and the brackets $\langle\ldots\rangle$ indicate an
ensemble average over the initial conditions and the dynamics. The modulus of
the wavevector k is chosen to correspond to the first peak of the static
structure factor, $\lvert\text{\bf
k}\rvert\approx\frac{2\pi}{\overline{\sigma}}$. From the time decay of
$F_{s}$, the relaxation time $\tau_{\alpha}$ is defined as $F_{s}(\text{\bf
k},t=\tau_{\alpha})\equiv\text{e}^{-1}$.

Figure 4: Relative efficiency of nonequilibrium ECMC versus equilibrium MMC.
The speedup due to collective chain moves decreases slowly with the density.

For a given packing fraction, we prepare equilibrated samples by performing
simulations in the $NVT$ ensemble using the cSwapECMC algorithm for a duration
of $\sim 10\tau_{\alpha}$. Using these equilibrated samples as initial
conditions, we then perform $NVT$ simulations using the six different
algorithms (MMC, ECMC, Swap, SwapECMC, cSwap, cSwapECMC) in the stationary
regime. We average $F_{s}(\text{\bf k},t)$ over 40 independent runs for ECMC
and MMC, and over 20 independent runs for the other algorithms. We also
average over all the wavevector of the same modulus. From these curves, we
extract $\tau_{\alpha}$ by means of an interpolation on a logarithmic time
axis. It is known that $\tau_{\alpha}$ is a good indicator of how fast
structural relaxation occurs at large $\phi$. Therefore the relative
efficiency of the Monte Carlo algorithms can be obtained by comparing how fast
the various algorithms decorrelate the structure of the system.

The relaxation times for the different algorithms as a function of the packing
fraction $\phi$ are shown in Fig. 3. Clearly, the curves split into two
families: one in which only translational moves are performed (MMC, ECMC) and
one in which diameter moves are combined with particle translations (Swap,
SwapECMC, cSwap, cSwapECMC). The latter is much more efficient that the
former.

Full relaxation of the system using translational moves only becomes hard when
$\phi\approx 0.61$. Although both are close, ECMC maintains an edge over MMC
across the whole range of packing fractions investigated in our work. However,
the gap between the two algorithms becomes smaller as we probe higher
densities. This is illustrated in Fig 4, where the ratio of the relaxation
times for the two algorithms is plotted as a function of $\phi$, reaching a
factor about 6 near $\phi=0.61$.

Figure 5: Efficiency of SwapECMC, cSwap and cSwapECMC relative to the
Metropolis Swap algorithm, as a function of the packing fraction $\phi$.
Moving beyond Metropolis systematically improves the performances of Monte
Carlo algorithms at large $\phi$.

When diameter swaps moves are introduced, it is possible to equilibrate
samples up to much higher densities near $\phi\approx 0.65$. The relative
efficiency of the different swap Monte Carlo algorithms is presented in Fig.
5. Comparing Swap to SwapECMC shows that replacing equilibrium local particle
translations with collective chain moves marginally improves the efficiency,
reducing the relaxation time by a factor $3$, which does not depend sensibly
on $\phi$. On the other hand, the efficiency of cSwap increases over that of
Swap with density, reaching a factor $6$ at $\phi=0.648$. Finally, cSwapECMC
combines the advantage of the nonequilibrium nature of both diameter exchanges
and translational moves, reaching a speedup of $20$ at $\phi=0.648$,
suggesting that the two speedups are nearly multiplicative.

The results obtained in this Section are very similar to our recent findings
for the two-dimensional hard-disk system [53].

##  VI Insights from microscopic dynamics

To gain physical insights into the relative performances of the various Monte
Carlo algorithms, we measure their respective time correlation functions and
compare their behaviors.

###  VI.1 Mean-squared displacement

Figure 6: The mean-squared displacement, Eq. (4), for MMC and ECMC at
$\phi=0.604$ are very similar, apart from a global time rescaling.

The mean-squared displacement $\Delta r^{2}(t)$ of the particles relative to
the center of mass is defined as

$\Delta r^{2}(t)\equiv\frac{1}{N}\sum_{i}\langle\left[\Delta\text{\bf
r}_{i}(t)\right]^{2}\rangle.$ (4)

The behavior of $\Delta r^{2}(t)$ for ECMC and MMC is compared at high packing
fraction in Fig. 6. Both curves show a very similar time dependence. They
display a short time diffusive behavior followed by a plateau at intermediate
times, reflecting particle caging. They eventually display diffusive behavior
again at long times. The plateau height is similar in both curves. These
curves seemingly only differ by a global rescaling of the time, implying that
both the short-time and long-time diffusion constants are accelerated by a
similar amount when event-chain translational moves are introduced. The
similarity between these two functions suggests that with the ECMC algorithm
the system follows dynamic relaxation pathways that are similar to the ones of
followed when resorting to MMC, up to the fact that these paths are explored
faster.

Figure 7: Mean-squared displacement, Eq. (4) at high density $\phi=0.672$ for
MMC, ECMC, Swap, cSwap and cSwapECMC. Swap moves increase the plateau height
corresponding to caged particle motion.

We also investigate the mean-squared displacement at a much larger packing
fraction, $\phi=0.672$, where equilibration was achieved using the cSwapECMC
algorithm which is the most efficient. We use 15 independent equilibrated
configurations obtained via cSwapECMC as initial conditions for trajectories
run with the other algorithms. The results, shown in Fig. 7, are therefore
representative of the stationary relaxation dynamics explored by the various
algorithms over the time window of the simulations [30].

At this very large density, MMC and ECMC are totally unable to relax the
system, and particles simply explore their cages over the entire time window,
thus producing a very clear plateau in the mean-squared displacements. Again,
ECMC reaches the plateau a bit faster than MMC, but the plateau values are
identical.

The situation is different for Swap, which closely follows the MMC data at
very short times but reaches a plateau value that is about three times as
large as that of MMC and ECMC. This suggests that particle swaps dramatically
enhance the amplitude of the in-cage motion performed by the particles, as
noted before [35]. Moving from Swap to SwapECMC produces a slightly faster
approach to a plateau value that is similar in both cases. Finally, moving
from SwapECMC to the most efficient cSwapECMC, we observe a similar time-
dependence at short times, but the plateau height is also slightly larger,
with a long-time dynamics that starts to enter the diffusive regime before all
other algorithms, which is consistent with the speedup observed by measuring
$\tau_{\alpha}$ in Section V.

###  VI.2 Non-Gaussian parameter

Figure 8: Non-Gaussian parameter $\alpha_{2}(t)$ for MMC and ECMC. for each
algorithm, the time has been rescaled by the corresponding relaxation time
$\tau_{\alpha}$.

The decreasing efficiency of ECMC relative to MMC at high packing fraction
suggests that the collective displacement of particle chains struggle to
produce displacements that are significantly different from the one produced
by MMC. To confirm this picture at the level of single particle displacements,
we measure the non-Gaussian parameter $\alpha_{2}(t)$, defined as

$\alpha_{2}(t)\equiv\frac{3\sum_{i}\langle\Delta\text{\bf
r}_{i}(t)^{4}\rangle}{5\sum_{i}\langle\Delta\text{\bf
r}_{i}(t)^{2}\rangle^{2}}-1.$ (5)

This quantity measures the deviation of the probability distribution function
of particle displacements from a Gaussian. It is a proxy to quantify the
heterogeneity of single-particle displacements in glassy systems [57, 58].

A plot of $\alpha_{2}(t)$ for MMC, ECMC, Swap and cSwap is shown in Fig. 8. In
this graph, the time is rescaled by the relaxation time $\tau_{\alpha}$ of
each algorithm. Unlike for Molecular Dynamics, the value of $\alpha_{2}(t)$ at
very short times is non zero since both the Metropolis translations and the
ECMC chain displacements are, by construction, non-Gaussians in this regime.
Moreover, the value of $\alpha_{2}(t)$ at short times is larger for ECMC than
MMC, as a reflection of the physical differences between the type of moves
performed by the two algorithms. At larger times, however the two functions
for ECMC and MMC, become virtually identical with a similar non-monotonic time
dependence and a peak of equal amplitude near $t/\tau_{\alpha}\sim 0.1$. This
suggests that the nature of single particle displacement at large times is
similar for ECMC and MMC algorithms. This is consistent with the picture
suggested by the mean-squared displacement of Sec. VI.1, and the displacements
fields obtained in two dimensions [53]. The behavior for the Swap and cSwap
algorithms is instead very different. At short times, the non-Gaussian
parameter has a similar value for both algorithms, since they both exploit
Metropolis Monte Carlo moves for the particle translations. The non-monotonic
dependence near the relaxation time is much weaker, with a much smaller
maximum. Upon rescaling time by $\tau_{\alpha}$, the curves for Swap and cSwap
are close, with cSwap being even smaller than Swap. The impressive change in
the behavior of the non-Gaussian parameter upon introducing Swap moves is
consistent with the analysis of the mean-squared displacement in Sec. VI.1,
and suggests that swap moves considerably decrease the dynamic heterogeneity
[35].

###  VI.3 Diameter dynamics

To understand the effect of collective swap moves introduced in cSwap, we
focus on the dynamics of the particle diameters. As explained in Sec. III,
swap Monte Carlo moves can be interpreted as providing a time dependence to
the particle diameters, promoting $\sigma_{i}$ from a constant to a
fluctuating degree of freedom [36].

Figure 9: Time series of the diameter of three tagged particles for (a) Swap
and (b) cSwap dynamics in the dilute fluid. The continuous time random walk
picture for diameter dynamics in Swap becomes directed as a result of the
collective swap moves in cSwap.

In Fig. 9 we show short time series of the diameter $\sigma_{i}(t)$ for three
randomly-selected tagged particles at a relatively low packing fraction,
$\phi=0.2$ for both Swap and cSwap algorithms. Since the particles diameters
do not change when translational moves are performed, in the time series only
swap moves are taken into account. The time $t$ in Fig. 9 thus counts the
number of cSwap or Swap moves performed. Under Swap dynamics, the diameters
perform a continuous time random walk. There are moments where the diameter is
constant either because particle $i$ is not selected, or because swap moves
are not accepted. When the diameter changes as a result of a successful swap,
its value can increase or decrease, so that at long times the diameter
performs a diffusive exploration of the particle size distribution
$\pi(\sigma)$.

The dynamics is very different under cSwap dynamics, where a form of directed
motion is apparent. Each time series now display a succession of continuous
decrease of the diameter, interrupted by sharp inflations. These upwards jumps
correspond to the situation where the tagged particle has inherited the
activity tag and a collective swap move is performed from this particle. In
between the sudden expansion steps, the diameter undergoes small deflations
and moments of stasis. The deflations take place when the particle participate
in a collective swap move, while being inactive. The moments of stasis happen
when the particle is neither active nor involved in a collective swap.

Figure 10: Evolution of $\langle n_{s}\rangle$, the average number of
particles involved in collective swap moves with the packing fraction.

Unlike ECMC, within cSwap there is no need to prescribe a length for the
collective moves. Instead, the number of particles $n_{s}$ participating to a
collective move is set by the maximal possible expansion of the active
particle. Clearly, this expansion must be a function of the packing fraction.
To analyze this, we measure $\langle n_{s}\rangle$, the average number of
particles that participate in collective swaps, see Fig. 10. While $\langle
n_{s}\rangle\sim N$ in the dilute limit $\phi\to 0$, this number decreases
when the packing fraction increases. We observe $\langle n_{s}\rangle\approx
190$ at $\phi=0.20$ and a much smaller number, $\langle n_{s}\rangle\approx
30$ near $\phi=0.6$ where MMC becomes unable to equilibrate. Near $\phi=0.65$
where cSwap also becomes very slow, the collective swap moves still involve a
significant number of particles, $\langle n_{s}\rangle\approx 20$. Recalling
that in addition collective swaps are never rejected, the enhanced efficiency
of cSwap over swap is easily accounted for by these observations.

Figure 11: Time relaxation of $C_{\sigma}(t)$ defined in Eq. (6) for Swap,
SwapECMC, cSwap and cSwapECMC algorithms. The dynamics of diameter
fluctuations follows the translational dynamics.

Finally, we track the diameter dynamics by measuring the time autocorrelation
function of diameter fluctuations [35], defined as

$C_{\sigma}(t)\equiv\left\langle\frac{\sum_{i}\delta\sigma_{i}(t)\delta\sigma_{i}(0)}{\sum_{i}\delta\sigma_{i}(0)^{2}}\right\rangle,$
(6)

with $\delta\sigma_{i}(t)\equiv\sigma_{i}(t)-\overline{\sigma}$. The time
evolution of $C_{\sigma}(t)$ at a fixed packing fraction $\phi=0.642$ is shown
in Fig. 11. We choose this packing fraction because a full decorrelation can
be still be observed for the four Monte Carlo algorithms containing swap
moves.

This figure demonstrates that the hierarchy of algorithmic efficiency studied
in Sec. V is respected when considering the relaxation dynamics of diameter
fluctuations. Interestingly, the decay of $C_{\sigma}(t)$ is accelerated for
SwapECMC with respect to Swap, even though the only difference between the two
algorithms lies in the position sector. This confirms the physical idea that
the acceleration provided by the swap moves comes from a dynamic interplay of
particle translations and diameter fluctuations [35], so that accelerating one
type of degrees of freedom also speeds up all others.

##  VII Monte Carlo algorithms to produce jammed sphere packings

All Monte Carlo algorithms can be used to prepare jammed configurations of
hard spheres [59]. Starting from an equilibrated configuration in the fluid at
$\phi_{\text{init}}$, we suddenly apply a very large pressure using $NPT$
simulations combined with MMC, Swap and cSwap algorithms, so that only volume
moves with $\delta V<0$ are accepted. The high applied pressure is $\beta
P\overline{\sigma}^{3}=10^{12}$, and the maximum step for box moves is
$V_{\text{max}}=5\times 10^{-5}V_{\text{init}}$, with $V_{\text{init}}$ the
volume of the initial equilibrium configuration at $\phi_{\text{init}}$. We
carry out compressions using the MMC, Swap and cSwap algorithm for a time
$t=4\times 10^{7}$. As the compression proceeds and the packing fraction
increases, most of the translational MMC moves and the box moves are rejected,
slowing down the compression process. To compensate for this effect, we
implement and adaptive step technique: we track the running average
${p_{\text{acc}}}$ of the acceptance rate for Metropolis translational moves.
When this running average drops below a threshold ${p}_{\text{th}}=0.2$, we
rescale the step size for proposed translational volume moves by a factor $2$.
The threshold for the running average of the acceptance rate is rescaled as
well, making this procedure iterative: $\delta\leftarrow\delta/2$, $\delta
V_{\text{max}}\leftarrow\delta V_{\text{max}}/2$,
${p}_{\text{th}}\leftarrow{p}_{\text{th}}/2$.

Figure 12: Jamming packing fraction achieved through high-pressure
compressions in the NPT ensemble, starting from configurations equilibrated at
$\phi_{\text{init}}$ using MMC, Swap and cSwap.

After the pressure quench, the packing fraction of the system slowly increases
until saturation to a packing fraction that corresponds to the jamming point,
$\phi_{\text{J}}$. By construction, $\phi_{\text{J}}$ may still depend on the
initial condition $\phi_{\rm init}$ and on all details of the Monte Carlo
algorithm employed to compress the system [60, 61, 62].

The evolution of $\phi_{\text{J}}$ with $\phi_{\text{init}}$ is shown in Fig.
12. For conventional MMC, $\phi_{\text{J}}$ exhibits a plateau at low
$\phi_{\text{init}}$, with a value $\phi_{\text{J,pl}}^{(\text{MMC})}\approx
0.66$. This value corresponds roughly to the jamming point defined by O’Hern
and coworkers [63] for the specific polydisperse model studied here. The
jamming packing fraction then increases when $\phi_{\text{init}}$ grows above
$\phi^{*(\text{MMC})}_{\text{init}}\approx 0.59$. The largest
$\phi_{\text{J}}$ reached with this protocol is
$\phi^{(\text{MMC})}_{\text{J}}\approx 0.6947$, which is consistent with
previous results [64].

Adding swap moves during the compression gives access to considerably denser
jammed packings [54, 65, 66]. The plateau height of the plateau at low
$\phi_{\rm init}$ becomes $\phi_{\text{J,pl}}^{(\text{Swap})}\approx 0.706$
and is slightly larger for cSwap, $\phi_{\text{J,pl}}^{(\text{cSwap})}\approx
0.709$. It is remarkable that a fast compression with elementary swap moves
from a dilute initial condition produces denser packings than a very painful
equilibration at large $\phi_{\rm init}$ followed by a compression involving
translational moves only [54].

Figure 13: Jammed packing at $\phi_{\text{J}}=0.7202\ldots$ obtained combining
compressions in the NPT ensemble with the cSwap algorithm.

Combining careful equilibration at larger $\phi_{\rm init}$ to compressions
involving swap moves produces even denser jammed configurations, reaching
$\phi_{\text{J}}^{\text{Swap}}\approx 0.719$. This value is slightly improved
when resorting to the cSwap algorithm during compressions,
$\phi_{\text{J}}^{\text{cSwap}}\approx 0.720$. This is significantly larger
than in previous efforts for the same model [64], although this value remains
smaller than the close packing density for a periodic crystal of monodisperse
spheres, $\phi\simeq 0.7404$. An example of a disordered very dense jammed
packing with $\phi_{\rm J}=0.7202$ is shown in Fig. 13.

The small but systematic dependence of $\phi_{\text{J}}$ on $\phi_{\rm init}$
in Fig. 12 for Swap and cSwap seems qualitatively different from the results
we have obtained in two dimensional hard disks where this dependence was
negligible [53]. As a result, for two-dimensional systems a rapid compression
with swap moves from a dilute or a very glassy initial conditions led to the
same jamming packing fraction [53, 67]. This is no longer the case in three
dimensions, where exploring dense equilibrium states in the fluid helps
producing denser jammed packings.

To confirm this idea, we compare the range of jamming packing fractions
obtained using our Monte Carlo methods to the packings obtained by gradient-
descent methods combining translational and diameter moves [65, 66]. We
obtained a jammed polydisperse configuration from the work of Hagh et al. [66]
which has a different particle size distribution and a polydispersity
$\Delta\approx 25\%$. The augmented gradient-descent protocol yields
$\phi_{\text{J}}^{gd}\approx 0.712$. Using this jammed configuration as a
starting point, we first decompress the system to equilibrate the system in
the dilute fluid and lose all memory of the initial jammed packing. When we
rapidly compress the system from a dilute fluid at $\phi_{\text{init}}=0.52$
we get $\phi_{\text{J,pl}}^{(\text{MMC})}\approx 0.66$, which corresponds to
the jamming point of O’Hern et al. This becomes
$\phi_{\text{J,pl}}^{(\text{Swap})}\approx 0.710$ when swap moves are
introduced in rapid compressions from dilute initial conditions, which is
roughly consistent with the gradient-descent method. We then carefully
equilibrate the system using $NVT$ simulations and the cSwapECMC algorithm at
$\phi_{\text{init}}=0.652$ and use these equilibrium configurations as initial
conditions for compressions using cSwap. The additional effort spent in
preparing a dense equilibrated fluid configuration allows us to reach a
jamming packing fraction $\phi_{\text{J}}^{\text{cSwap}}\approx 0.717$, which
is denser than the result obtained from gradient-descent methods alone.

##  VIII Conclusion

The design of efficient sampling algorithms in systems such as glass-formers
whose dynamics is controlled by proliferating dynamical bottlenecks has been a
challenge on its own since computing hardware has become part of the standard
toolbox of theoreticians. When it comes to sampling a set of degrees of
freedom according to a given distribution, algorithmic progress is in
principle limited by strict mathematical bounds, as discussed in [68], and it
is in general simply not known how close one is from the optimally fastest
algorithm. This algorithm chimera can be approached on the basis of physical
intuition, along with trial and error. But the rules of the game can also be
tweaked by introducing out-of-the-box degrees of freedom. This can be done
reversibly: this is what the swap algorithm achieves by turning particles into
objects of fluctuating sizes. Or this can be done irreversibly, as with event-
chain or lifting methods, which exploit the ballistic acceleration of
correlated clusters of particles. The present work has presented how a
combination of these methods can bring significant speed-up gains even in such
uncooperative systems as glass-formers.

However, much remains to be done, and to begin with, adapting our methods to
particles interacting by means of soft continuous potentials stands out as a
primary goal. Swap performances do not decrease with continuous potentials
[35], and some extensions of ECMC to continuous potentials were successful
[52, 69]. The most salient challenge will be to design a version of cSwap for
glass-formers with continuous potentials to extend the range of applicability
of irreversible Monte Carlo methods in the field of supercooled liquids. Yet
implementing collective swaps involving finite energy variations will require
care so as to preserve global balance. All these perspectives directly follow
from our work; they should help the development of efficient and versatile
sampling methods for systems with a complex energy landscape.

###### Acknowledgements.

We thank V. F. Hagh and E. I. Corwin for providing us one of the hard spheres
systems studied in this work. We also thank Y. Nishikawa for interesting
exchanges on cSwap. We acknowledge the financial support of the ANR THEMA
AAPG2020 grant.

## Author declarations

### Conflicts of interest

The authors have no conflicts of interest to disclose.

### Author contributions

Ludovic Berthier Conceptualization (equal), Supervision (equal),
writing/review & editing (equal). Federico Ghimenti: Software, writing/review
& editing (equal). Frédéric van Wijland: Conceptualization (equal),
Supervision (equal), writing/review & editing (equal).

## Data availability

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

## References

  * Ediger  _et al._ [1996] M. D. Ediger, C. A. Angell, and S. R. Nagel, Supercooled liquids and glasses, The journal of physical chemistry 100, 13200 (1996).
  * Cavagna [2009] A. Cavagna, Supercooled liquids for pedestrians, Physics Reports 476, 51 (2009).
  * Berthier and Ediger [2016] L. Berthier and M. D. Ediger, Facets of glass physics, Physics today 69, 40 (2016).
  * Götze [2009] W. Götze,  _Complex dynamics of glass-forming liquids_ , Vol. 143 (Oxford University Press on Demand, 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).
  * Parisi  _et al._ [2020] G. Parisi, P. Urbani, and F. Zamponi,  _Theory of simple glasses_ (Cambridge University Press, 2020).
  * Berthier and Reichman [2023] L. Berthier and D. R. Reichman, Modern computational studies of the glass transition, Nature Reviews Physics 5, 1 (2023).
  * Metropolis  _et al._ [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Equation of State Calculations by Fast Computing Machines, The Journal of Chemical Physics 21, 1087 (1953).
  * Alder and Wainwright [1959] B. J. Alder and T. E. Wainwright, Studies in Molecular Dynamics. I. General Method, The Journal of Chemical Physics 31, 459 (1959).
  * Barrat and Roux [1991] J.-L. Barrat and J.-N. Roux, Molecular dynamics simulations of supercooled liquids near the glass transition, Journal of non-crystalline solids 131, 255 (1991).
  * Newman and Barkema [1999] M. E. Newman and G. T. Barkema,  _Monte Carlo methods in statistical physics_ (Clarendon Press, 1999).
  * Landau and Binder [2021] D. Landau and K. Binder,  _A guide to Monte Carlo simulations in statistical physics_ (Cambridge university press, 2021).
  * Frenkel and Smit [2001] D. Frenkel and B. Smit,  _Understanding molecular simulation: from algorithms to applications_ , Vol. 1 (Elsevier, 2001).
  * Hansen and McDonald [2013] J.-P. Hansen and I. R. McDonald,  _Theory of simple liquids: with applications to soft matter_ (Academic press, 2013).
  * Royall and Williams [2015] C. P. Royall and S. R. Williams, The role of local structure in dynamical arrest, Physics Reports 560, 1 (2015).
  * Nussinov [2005] Z. Nussinov, The frustration-based approach of supercooled liquids and the glass transition: a review and critical assessment, Journal of Physics: Condensed Matter 17, R1143 (2005).
  * Bouchaud and Biroli [2004] J.-P. Bouchaud and G. Biroli, On the Adam-Gibbs-Kirkpatrick-Thirumalai-Wolynes scenario for the viscosity increase in glasses, The Journal of Chemical Physics 121, 7347 (2004).
  * Biroli  _et al._ [2008] G. Biroli, J.-P. Bouchaud, A. Cavagna, T. S. Grigera, and P. Verrocchio, Thermodynamic signature of growing amorphous order in glass-forming liquids, Nature Physics 4, 771 (2008).
  * Berthier and Kob [2012] L. Berthier and W. Kob, Static point-to-set correlations in glass-forming liquids, Phys. Rev. E 85, 011102 (2012).
  * Sciortino  _et al._ [1999] F. Sciortino, W. Kob, and P. Tartaglia, Inherent structure entropy of supercooled liquids, Phys. Rev. Lett. 83, 3214 (1999).
  * Sastry [2001] S. Sastry, The relationship between fragility, configurational entropy and the potential energy landscape of glass-forming liquids, Nature 409, 164 (2001).
  * Berthier  _et al._ [2019a] L. Berthier, M. Ozawa, and C. Scalliet, Configurational entropy of glass-forming liquids, The Journal of Chemical Physics 150, 160902 (2019a).
  * Franz and Parisi [1997] S. Franz and G. Parisi, Phase diagram of coupled glassy systems: A mean-field study, Phys. Rev. Lett. 79, 2486 (1997).
  * Berthier and Coslovich [2014] L. Berthier and D. Coslovich, Novel approach to numerical measurements of the configurational entropy in supercooled liquids, Proceedings of the National Academy of Sciences 111, 11668 (2014).
  * Guiselin  _et al._ [2022] B. Guiselin, L. Berthier, and G. Tarjus, Statistical mechanics of coupled supercooled liquids in finite dimensions, SciPost Physics 12, 091 (2022).
  * Gleim  _et al._ [1998] T. Gleim, W. Kob, and K. Binder, How does the relaxation of a supercooled liquid depend on its microscopic dynamics?, Phys. Rev. Lett. 81, 4404 (1998).
  * Berthier and Kob [2007] L. Berthier and W. Kob, The monte carlo dynamics of a binary lennard-jones glass-forming mixture, Journal of Physics: Condensed Matter 19, 205130 (2007).
  * Berthier [2007a] L. Berthier, Revisiting the slow dynamics of a silica melt using monte carlo simulations, Phys. Rev. E 76, 011507 (2007a).
  * Berthier [2007b] L. Berthier, Efficient measurement of linear susceptibilities in molecular simulations: Application to aging supercooled liquids, Phys. Rev. Lett. 98, 220601 (2007b).
  * Scalliet  _et al._ [2022] C. Scalliet, B. Guiselin, and L. Berthier, Thirty milliseconds in the life of a supercooled liquid, Phys. Rev. X 12, 041028 (2022).
  * Barrat and Berthier [2023] J.-L. Barrat and L. Berthier, Computer simulations of the glass transition and glassy materials, Comptes Rendus. Physique 24, 1 (2023).
  * Dress and Krauth [1995] C. Dress and W. Krauth, Cluster algorithm for hard spheres and related systems, Journal of Physics A: Mathematical and General 28, L597 (1995).
  * Kranendonk and Frenkel [1991] W. Kranendonk and D. Frenkel, Computer simulation of solid-liquid coexistence in binary hard sphere mixtures, Molecular physics 72, 679 (1991).
  * Grigera and Parisi [2001] T. S. Grigera and G. Parisi, Fast monte carlo algorithm for supercooled soft spheres, Phys. Rev. E 63, 045102 (2001).
  * Ninarello  _et al._ [2017] A. Ninarello, L. Berthier, and D. Coslovich, Models and algorithms for the next generation of glass transition studies, Phys. Rev. X 7, 021039 (2017).
  * Berthier  _et al._ [2019b] L. Berthier, E. Flenner, C. J. Fullerton, C. Scalliet, and M. Singh, Efficient swap algorithms for molecular dynamics simulations of equilibrium supercooled liquids, Journal of Statistical Mechanics: Theory and Experiment 2019, 064004 (2019b).
  * Berthier  _et al._ [2019c] L. Berthier, P. Charbonneau, and J. Kundu, Bypassing sluggishness: Swap algorithm and glassiness in high dimensions, Phys. Rev. E 99, 031301 (2019c).
  * Parmar  _et al._ [2020] A. D. S. Parmar, M. Ozawa, and L. Berthier, Ultrastable metallic glasses in silico, Phys. Rev. Lett. 125, 085505 (2020).
  * Alvarez-Donado  _et al._ [2023] R. Alvarez-Donado, S. Bonfanti, and M. Alava, Simulated multi-component metallic glasses akin to experiments, arXiv preprint arXiv:2309.05806 (2023).
  * Hukushima and Nemoto [1996] K. Hukushima and K. Nemoto, Exchange monte carlo method and application to spin glass simulations, Journal of the Physical Society of Japan 65, 1604 (1996).
  * Yamamoto and Kob [2000] R. Yamamoto and W. Kob, Replica-exchange molecular dynamics simulation for supercooled liquids, Phys. Rev. E 61, 5473 (2000).
  * Jung  _et al._ [2024] G. Jung, G. Biroli, and L. Berthier, Normalizing flows as an enhanced sampling method for atomistic supercooled liquids, arXiv preprint arXiv:2404.09914 (2024).
  * Krauth [2006] W. Krauth,  _Statistical mechanics: algorithms and computations_ , Vol. 13 (OUP Oxford, 2006).
  * Bernard  _et al._ [2009] E. P. Bernard, W. Krauth, and D. B. Wilson, Event-chain monte carlo algorithms for hard-sphere systems, Phys. Rev. E 80, 056704 (2009).
  * Turitsyn  _et al._ [2011] K. S. Turitsyn, M. Chertkov, and M. Vucelja, Irreversible monte carlo algorithms for efficient sampling, Physica D: Nonlinear Phenomena 240, 410 (2011).
  * Chen  _et al._ [1999] F. Chen, L. Lovász, and I. Pak, Lifting markov chains to speed up mixing, in  _Proc. 31st Ann. ACM Symp. Theory of Comp._ (1999) pp. 275–281.
  * Diaconis  _et al._ [2000] P. Diaconis, S. Holmes, and R. M. Neal, Analysis of a nonreversible markov chain sampler, Annals of Applied Probability , 726 (2000).
  * Sakai and Hukushima [2013] Y. Sakai and K. Hukushima, Dynamics of one-dimensional ising model without detailed balance condition, Journal of the Physical Society of Japan 82, 064003 (2013).
  * Sakai and Hukushima [2016] Y. Sakai and K. Hukushima, Eigenvalue analysis of an irreversible random walk with skew detailed balance conditions, Phys. Rev. E 93, 043318 (2016).
  * Vucelja [2016] M. Vucelja, Lifting—A nonreversible Markov chain Monte Carlo algorithm, American Journal of Physics 84, 958 (2016).
  * Bernard and Krauth [2011] E. P. Bernard and W. Krauth, Two-step melting in two dimensions: First-order liquid-hexatic transition, Phys. Rev. Lett. 107, 155704 (2011).
  * Michel  _et al._ [2014] M. Michel, S. C. Kapfer, and W. Krauth, Generalized event-chain Monte Carlo: Constructing rejection-free global-balance algorithms from infinitesimal steps, The Journal of Chemical Physics 140, 054116 (2014).
  * Ghimenti  _et al._ [2024] F. Ghimenti, L. Berthier, and F. van Wijland, Irreversible monte carlo algorithms for hard disk glasses: from event-chain to collective swaps, arXiv preprint arXiv:2402.06585 (2024).
  * Berthier  _et al._ [2016] L. Berthier, D. Coslovich, A. Ninarello, and M. Ozawa, Equilibrium sampling of hard spheres up to the jamming density and beyond, Phys. Rev. Lett. 116, 238002 (2016).
  * Berthier  _et al._ [2017] L. Berthier, P. Charbonneau, D. Coslovich, A. Ninarello, M. Ozawa, and S. Yaida, Configurational entropy measurements in extremely supercooled liquids that break the glass ceiling, Proceedings of the National Academy of Sciences 114, 11356 (2017).
  * Boublík [1970] T. Boublík, Hard-sphere equation of state, The Journal of chemical physics 53, 471 (1970).
  * Rahman [1964] A. Rahman, Correlations in the motion of atoms in liquid argon, Phys. Rev. 136, A405 (1964).
  * Kob  _et al._ [1997] W. Kob, C. Donati, S. J. Plimpton, P. H. Poole, and S. C. Glotzer, Dynamical heterogeneities in a supercooled lennard-jones liquid, Phys. Rev. Lett. 79, 2827 (1997).
  * Berthier and Witten [2009] L. Berthier and T. A. Witten, Glass transition of dense fluids of hard and compressible spheres, Phys. Rev. E 80, 021502 (2009).
  * Donev  _et al._ [2004] A. Donev, S. Torquato, F. H. Stillinger, and R. Connelly, Comment on “jamming at zero temperature and zero applied stress: The epitome of disorder”, Phys. Rev. E 70, 043301 (2004).
  * O’Hern  _et al._ [2004] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Reply to “comment on ‘jamming at zero temperature and zero applied stress: The epitome of disorder’ ”, Phys. Rev. E 70, 043302 (2004).
  * Chandler and Garrahan [2010] D. Chandler and J. P. Garrahan, Dynamics on the way to forming glass: Bubbles in space-time, Annual review of physical chemistry 61, 191 (2010).
  * O’Hern  _et al._ [2002] C. S. O’Hern, S. A. Langer, A. J. Liu, and S. R. Nagel, Random packings of frictionless particles, Phys. Rev. Lett. 88, 075507 (2002).
  * Ozawa  _et al._ [2017] M. Ozawa, L. Berthier, and D. Coslovich, Exploring the jamming transition over a wide range of critical densities, SciPost Physics 3, 027 (2017).
  * Brito  _et al._ [2018] C. Brito, E. Lerner, and M. Wyart, Theory for swap acceleration near the glass and jamming transitions for continuously polydisperse particles, Phys. Rev. X 8, 031050 (2018).
  * Hagh  _et al._ [2022] V. F. Hagh, S. R. Nagel, A. J. Liu, M. L. Manning, and E. I. Corwin, Transient learning degrees of freedom for introducing function in materials, Proceedings of the National Academy of Sciences 119, e2117622119 (2022).
  * Bolton-Lum  _et al._ [2024] V. Bolton-Lum, R. C. Dennis, P. Morse, and E. Corwin, The ideal glass and the ideal disk packing in two dimensions, arXiv preprint arXiv:2404.07492 (2024).
  * Krauth [2021] W. Krauth, Event-chain monte carlo: foundations, applications, and prospects, Frontiers in Physics 9, 663457 (2021).
  * Nishikawa  _et al._ [2023] Y. Nishikawa, W. Krauth, and A. C. Maggs, Liquid-hexatic transition for soft disks, Phys. Rev. E 108, 024103 (2023).


Identifying the multifractal set on which energy dissipates
in a turbulent Navier-Stokes fluid

John D. Gibbon

Department of Mathematics, Imperial College London, London SW7 2AZ, UK

###### Abstract

The rich multifractal properties of fluid turbulence illustrated by the work
of Parisi and Frisch are related explicitly to Leray’s weak solutions of the
three-dimensional Navier-Stokes equations. Directly from this correspondence
it is found that the set on which energy dissipates, $\mathbb{F}_{m}$, has a
range of dimensions $\mathfrak{D}_{m}=3/m$ ($1\leq m\leq\infty$), and a
corresponding range of sub-Kolmogorov dissipation inverse length scales
$L\eta_{m}^{-1}\leq Re^{3/(1+\mathfrak{D}_{m})}$ spanning $Re^{3/4}$ to
$Re^{3}$. Correspondingly, the multifractal model scaling parameter $h$, must
obey $h\geq h_{min}$ with $-{\scriptstyle\frac{2}{3}}\leq
h_{min}\leq{\scriptstyle\frac{1}{3}}$.

In memory of Charlie Doering (1956-2021)

##  1 The need for a mathematical telescope

Despite the cliché that “turbulence is the last great unsolved problem in
classical physics” [1], the challenge that it poses to physicists,
mathematicians and engineers remains [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
14, 15, 16, 17, 18, 19]. A natural correspondence is proposed that helps us
relate the rich multifractality of the statistical properties of homogeneous
and isotropic fluid turbulence [2, 3, 4] to the methods of mathematical
analysis which are used to study weak solutions of the incompressible Navier-
Stokes equations [20, 21, 22, 24, 23, 25, 26, 27, 28], including the Millenial
regularity problem [29]. Heretofore the two approaches have taken divergent
paths which has caused a lacuna in our understanding. Multifractality is
manifest in the wide sweep of vortical structures that appear in turbulent
flows, ranging from three-dimensional vortices at many scales, to semi-broken
filaments whose fractal dimension is less than unity : see Fig. 1 for
illustrative examples. In telescopic terms, the simultaneous existence of many
vortical structures each with its own fractal dimension requires the
mathematical equivalent of an adjustable focus, somewhat in the manner of a
zoom lens. Thus, an adjustable parameter with a wide range of values is
needed. The multifractal model (MFM) of Parisi and Frisch possesses this
latter feature through the allowed variation of its scaling parameter $h$ [2],
while the Navier-Stokes equations seemingly do not. The MFM was developed
after it had become clear that Kolmogorov’s 1941 theory needed to be modified
in order to take account of the effects of intermittency [3] : see also [30,
31, 32, 33, 34, 35].

Figure 1: (left) Iso-surface plot of the magnitude of the vorticity
$|\mbox{\boldmath$\omega$}|$ (two standard deviations above the mean) from a
$512^{3}$ pseudospectral direct numerical simulation of the unforced $3D$
Navier-Stokes equations at a representative time at an integral-scale Reynolds
number of 874. Plot courtesy of Nadia B. Padhan, IISc, Bangalore. The second
and third are filled-contour plots of the vorticity from $2048^{2}$
pseudospectral direct numerical simulations of the forced $2D$ Navier-Stokes
equations in a statistically steady state, courtesy of Kiran Kolluru, IISc,
Bangalore : (middle) forward-cascade-dominated turbulence
($Re_{\lambda}=1957$) ; (right) inverse-cascade-dominated turbulence
($Re_{\lambda}=1729$).

The question is this : can this telescopic property be derived from Leray’s
weak solution formulation of the three-dimensional incompressible Navier-
Stokes equations [20, 21, 22, 24, 23, 25, 26, 27, 28]? From the Navier-Stokes
side of the fence, no mathematical tools seemingly exist that would allow us
to perform rigorous analysis on fractal sets. Standard methods of Navier-
Stokes analysis involve integration over the full domain volume, which has the
unfortunate consequence of washing out the delicate nature of the set on which
energy dissipates [5, 6]. This paper circumvents these problems by showing how
the existence of a multifractal set on which energy dissipates is natural to
the problem. Results can be extracted from the Leray formalism which are then
compared to the multifractal theory of Parisi and Frisch [2].

In recent work, Dubrulle and Gibbon [36] addressed the issue of the
equivalence of the Navier-Stokes equations and the multifractal model. They
used the Paladin-Vulpiani111The argument used in [32, 33] to estimate this
scale is to equate the turnover time to the viscous diffusion time. It is also
based on a Reynolds number defined in terms of the energy dissipation rate.
[32, 33] inverse dissipation length scale

$L\eta^{-1}_{PV}\sim Re^{1/(1+h)}$ (1.1)

as an ad hoc bridge between the two theories, where $h$ is the multifractal
scaling parameter. In this paper we are able to go a step further by showing
that this scale, at least in an inequality form, appears naturally out of the
correspondence between the two theories without any ad hoc introduction.

##  2 Recasting Leray’s weak solution formulation of the Navier-Stokes
equations

###  2.1 A zoom lens

The incompressible Navier-Stokes equations are

$\left(\partial_{t}+\mbox{\boldmath$u$}\cdot\nabla\right)\mbox{\boldmath$u$}+\nabla
p=\nu\Delta\mbox{\boldmath$u$}+\mbox{\boldmath$f$}(\mbox{\boldmath$x$})$ (2.1)

subject to
$\mbox{div}\,\mbox{\boldmath$u$}=\mbox{div}\,\mbox{\boldmath$f$}=0$, on a
periodic domain $V_{d}=[0,\,L]^{d}$ for dimensions $d=3,\,2,\,1$. Results on
weak solutions can be expressed in many ways [20] but most appear in time
averaged form [21, 22, 24, 23, 25, 26, 27, 28]. In [27], it was shown that for
$1\leq m\leq\infty$, the Sobolev norms

$\|\nabla^{n}\mbox{\boldmath$u$}\|_{2m}=\left(\int_{V_{d}}|\nabla^{n}\mbox{\boldmath$u$}|^{2m}\,dV\right)^{1/2m}\,,$
(2.2)

under the Navier-Stokes invariance property
$\mbox{\boldmath$x$}^{\prime}=\lambda^{-1}\mbox{\boldmath$x$}$ ;
$t^{\prime}=\lambda^{-2}t$ and
$\mbox{\boldmath$u$}=\lambda^{-1}\mbox{\boldmath$u$}^{\prime}$, have the
scaling property

$\displaystyle\|\nabla^{n}\mbox{\boldmath$u$}\|_{2m}$ $\displaystyle=$
$\displaystyle\lambda^{-1/\alpha_{n,m,d}}\|\nabla^{{}^{\prime}n}\mbox{\boldmath$u$}^{\prime}\|_{2m}\,$
(2.3) $\displaystyle\alpha_{n,m,d}$ $\displaystyle=$
$\displaystyle\frac{2m}{2m(n+1)-d}\,.$ (2.4)

It has been shown in [27, 28] that for $n\geq 1$, the weighted time averages
of the dimensionless family

$F_{n,m,d}=\nu^{-1}L^{1/\alpha_{n,m,d}}\|\nabla^{n}\mbox{\boldmath$u$}\|_{2m}$
(2.5)

are finite such that222The proof of (2.6) in [27] depends heavily on the
result in [21]. The case $n=0$ is valid only for $m>3$.

$\left<F_{n,m,d}^{(4-d)\alpha_{n,m,d}}\right>_{T}<\infty\,.$ (2.6)

The brackets $\left<\cdot\right>_{T}$ are a time average up to time $T>0$. The
set of exponents $\alpha_{n,m,d}$ in (2.6) clearly play a significant role in
Navier-Stokes regularity properties333It was shown in [27] that to prove full
regularity in the $d=3$ case a factor of $2\alpha_{n,m,3}$ would be needed in
the exponent of (2.6). This result remains elusive.. A series of well known
weak solution results can be recovered by choosing various values of $n$ and
$m$, which all have the same level of regularity [27]. Increasing $m$
amplifies the larger features in the solution in proportion to the smaller
until $m=\infty$ is reached which corresponds to the maximum. Therefore,
adjusting the value of $m$ acts as a focus control on a telescope that allows
us to zoom into different features in a landscape. The subscript labelling on
$\alpha_{n,m,d}$ encodes the nature of the norm of which it is an exponent :
thus $(n,m,d)$ means “$n$ derivatives in $L^{2m}$ in integer-$d$ spatial
dimensions”. However, $\alpha_{n,m,d}$ can be rewritten as

$\displaystyle\alpha_{n,m,d}$ $\displaystyle=$
$\displaystyle\frac{2}{2(n+1)-\mathfrak{D}_{m}}\,,\qquad\mathfrak{D}_{m}=d/m\,,$
(2.7) $\displaystyle=$ $\displaystyle\alpha_{n,1,\mathfrak{D}_{m}}\,.$

The subscripts on $\alpha_{n,1,\mathfrak{D}_{m}}$ now suggest that the problem
can now be recast into “$n$ derivatives in $L^{2}$ on a set, designated as
$\mathbb{F}_{m}$, which has a range of dimensions $\mathfrak{D}_{m}=d/m$”. The
range of dimensions $\mathfrak{D}_{m}=d/m$ suggests that $\mathbb{F}_{m}$ is
multifractal in character, with a corresponding set of co-dimensions

$\mathfrak{C}_{m}=d\left(1-\frac{1}{m}\right)\,.$ (2.8)

Mimicking (2.6), the recasting of $\alpha_{n,m,d}$ into
$\alpha_{n,1,\mathfrak{D}_{m}}$ suggests that the family of time averages

$X_{n,\mathfrak{D}_{m}}(T)=\left<F_{n,1,\mathfrak{D}_{m}}^{(4-\mathfrak{D}_{m})\,\alpha_{n,1,\mathfrak{D}_{m}}}\right>_{T}\,,$
(2.9)

may be a key set of objects. The $X_{n,m,\mathfrak{D}_{m}}$ are neither
necessarily equal to the averages in (2.6) nor would it be an easy task to
prove they are bounded above because there exists no mechanism for performing
analysis on fractal or multifractal domains. Nevertheless, the case $n=1$ is
an exception in that it can be bounded because of the cancellation of the
factor of $(4-\mathfrak{D}_{m})$ in

$\left[(4-\mathfrak{D}_{m})\alpha_{n,1,\mathfrak{D}_{m}}\right]_{n=1}=2\,,\qquad
1\leq m\leq\infty\,,$ (2.10)

for every value of $0\leq\mathfrak{D}_{m}\leq d$. It is this case we consider
in the next subsection.

###  2.2 The three-dimensional case when $n=1$

The $n=1$ case of (2.9) for $d=3$ corresponds, for each $m$, to a time-
averaged energy dissipation rate $\varepsilon_{m}$ on the multifractal set
$\mathbb{F}_{m}$. Each $\varepsilon_{m}$ is defined as

$\varepsilon_{m}=\nu
L^{-\mathfrak{D}_{m}}\left<\int_{\mathbb{F}_{m}}|\nabla\mbox{\boldmath$u$}|^{2}dV\right>_{T}\,,\quad\mathfrak{D}_{m}=3/m\,.$
(2.11)

Noting that $\mathbb{F}_{m}\subset\mathbb{T}^{3}$, (2.11) can be re-written as

$\displaystyle\varepsilon_{m}$ $\displaystyle\leq$ $\displaystyle\nu
L^{3-\mathfrak{D}_{m}}L^{-3}\left<\int_{\mathbb{T}^{3}}|\nabla\mbox{\boldmath$u$}|^{2}dV\right>_{T}$
(2.12) $\displaystyle\leq$ $\displaystyle
c\,L^{-1-\mathfrak{D}_{m}}\nu^{3}Re^{3}\,.$

The $Re^{3}$ upper bound is derived for generic narrow-band forcing considered
by Doering and Foias [16] where the Reynolds number $Re$ is defined by

$Re=UL/\nu\qquad\mbox{with}\qquad
U^{2}=L^{-3}\left<\|\mbox{\boldmath$u$}\|^{2}_{2}\right>_{T}\,.$ (2.13)

(2.12) shows that $\varepsilon_{m}\nu^{-3}$ can be expressed in terms of a set
of inverse lengths $\eta_{m}^{-1}$ defined such that

$\varepsilon_{m}\nu^{-3}:=\eta_{m}^{-1-\mathfrak{D}_{m}}$ (2.14)

leading to

$L\eta_{m}^{-1}\leq c\,Re^{\frac{3}{1+\mathfrak{D}_{m}}}\,.$ (2.15)

(2.15) is a formula444It was common in the literature of that period to assign
a single value to the fractal dimension introduced from another source, as in
[31]. that appears in Kraichnan [31], but without any identification of the
nature of $\mathbb{F}_{m}$. The size of the upper bound on $L\eta_{m}^{-1}$
ranges from the standard inverse Kolmogorov bound $Re^{3/4}$ when $m=1$ to
$Re^{3}$ when $m=\infty$.

###  2.3 The two-dimensional case

What of the two-dimensional case? Advantage can be taken of the lack of a
vortex stretching term
($\mbox{\boldmath$\omega$}\cdot\nabla\mbox{\boldmath$u$}=0$) in the vorticity
formulation to find the enstrophy dissipation rate [16, 5]

$\nu
L^{-2}\left<\int_{\mathbb{T}^{2}}|\nabla\mbox{\boldmath$\omega$}|^{2}dV\right>_{T}\leq
c\,\nu^{3}L^{-6}Re^{3}\,.$ (2.16)

In a full $\mathbb{T}^{2}$ domain, (2.16) leads to an $Re^{1/2}$ bound on the
set of inverse Kraichnan lengths $L\eta_{m,ens}^{-1}$. The equivalent of
(2.11) for the enstrophy dissipation rate on $\mathbb{F}_{m}$ is

$\varepsilon_{m,ens}=\nu
L^{-\mathfrak{D}_{m}}\left<\int_{\mathbb{F}_{m}}|\nabla\mbox{\boldmath$\omega$}|^{2}dV\right>_{T}\,,$
(2.17)

with $\mathfrak{D}_{m}=2/m$. Then $\varepsilon_{m,ens}\nu^{-3}$ can be
expressed in terms of a set of inverse lengths $\eta_{m,ens}^{-1}$ such that
$\varepsilon_{m,ens}\nu^{-3}=\eta_{m,ens}^{-4-\mathfrak{D}_{m}}$, leading to

$L\eta_{m,ens}^{-1}\leq c\,Re^{\frac{3}{4+\mathfrak{D}_{m}}}\,,$ (2.18)

and thus a spread of $Re^{1/2}$ to $Re^{3/4}$. In this two-dimensional case
Alexakis and Doering [17] derived a tighter bound of $Re^{2}$ on the right
hand side of (2.16) when the forcing has a single wave-number or has a
constant flux. An amendation of the numerator in (2.18) gives a spread ranging
from $Re^{1/3}$ to $Re^{1/2}$, thus implying even less multifractality.

##  3 Correspondence with the multifractal model

Kolmogorov’s 1941 theory (K41) of three-dimensional, homogeneous, isotropic
turbulence [3] revolves around the $p$th-order velocity structure function
$S_{\\!p}$ defined as

$S_{\\!p}=\left<\left|\mbox{\boldmath$u$}(\mbox{\boldmath$x$}+\mbox{\boldmath$r$})-\mbox{\boldmath$u$}(\mbox{\boldmath$x$})\right|^{p}\right>_{stat.av.}\,.$
(3.1)

The scaling properties of $S_{\\!p}$ have their origin in the invariance
property enjoyed by the incompressible Euler equations under the
transformation $\mbox{\boldmath$x$}^{\prime}=\lambda^{-1}\mbox{\boldmath$x$}$
; $t^{\prime}=\lambda^{h-1}t$ and
$\mbox{\boldmath$u$}=\lambda^{h}\mbox{\boldmath$u$}^{\prime}$. The scaling
parameter $h$ can take any value, unlike the Navier-Stokes equations where it
must satisfy $h=-1$. The invariance property shows that $S_{\\!p}$ scales as
$S_{\\!p}\sim r^{hp}$. As explained in [3], the axioms of K41 lead to
$h={\scriptstyle\frac{1}{3}}$, which is consistent with the “four-fifths law”
$S_{3}=-{\scriptstyle\frac{4}{5}}\varepsilon r$. The problem with the result
$S_{\\!p}\sim r^{p/3}$ is that experimental data deviate from the linear in
$p$ scaling but instead lie on a concave curve $\zeta_{p}$ below the line
${\scriptstyle\frac{1}{3}}p$. The multifractal model of Parisi and Frisch [2,
3] adapted the K41 formalism by relaxing the requirement that
$h={\scriptstyle\frac{1}{3}}$ by insisting that the energy dissipation rate is
invariant in $h$ only in an average sense. They achieved this by introducing
the probability $P(h)\sim r^{3-D(h)}$ of observing a given scaling exponent
$h$ at the scale $r$. Each value of $h$ belongs to a given fractal set of
dimension $D(h)$. This procedure produces the result

$S_{\\!p}\sim
r^{\zeta_{p}}\quad\mbox{where}\quad\zeta_{p}=\min_{h}\left[hp+3-D(h)\right]\,.$
(3.2)

This must be constrained by the four-fifths law which requires $\zeta_{3}=1$,
thus leading to the inequalities for $D(h)$ and the co-dimension $C(h)$ (which
sum to the spatial dimension)

$D(h)\leq 3h+2\quad\Rightarrow\quad C(h)\geq 1-3h\,.$ (3.3)

Given that K41 and the MFM models revolve around the flow being homogeneous
and isotropic, any comparison with results from the Navier-Stokes equations
will likely only be valid when the flow has reached a fully turbulent state :
i.e. when $T$ is sufficiently large.

To make the two theories touch, it is natural to suggest a correspondence
between $\mathfrak{D}_{m}$ and $D(h)$ and $\mathfrak{C}_{m}$ and $C(h)$, which
leads to a simple inequality relating $m$ and $h$

$\frac{3}{m}\leq 2+3h\,,$ (3.4)

from which we find $h\geq h_{min}$ with
$h_{min}=m^{-1}-{\scriptstyle\frac{2}{3}}$. The two opposite limits $m\to 1$
and $m\to\infty$ give

$-{\scriptstyle\frac{2}{3}}\leq h_{min}\leq{\scriptstyle\frac{1}{3}}\,,$ (3.5)

which is precisely the range found in [36].

A more modern version of the multifractal model uses the theory of large
deviations to re-express $P(h)$ as $P(h)\sim r^{C(h)}$, where $C(h)$ is the
multifractal spectrum [3, 35, 6]. $C(h)$ therefore plays the role as the co-
dimension. In the theory of large deviations it is possible that $C(h)\geq 3$,
the domain dimension, corresponding to $D(h)<0$, but only with an
infinitesimal probability. This would correspond to $m<0$, which must be
excluded here. The exponent $3/(1+\mathfrak{D}_{m})$ in (2.15) is bounded
below by

$\frac{3}{1+\mathfrak{D}_{m}}\geq\frac{1}{1+h}\,.$ (3.6)

The left hand side of (3.6) is derived directly from the Navier-Stokes
equations, while the right hand side comes from the MFM and can be recognized
as the Paladin-Vulpiani inverse dissipation scale $\eta^{-1}_{PV}$ defined in
(1.1). Paladin and Vulpiani [32] treated (3.3) as an equality, which would
change the $\geq$ in (3.6) to equality. Thus we see that the Paladin-Vulpiani
scale appears naturally in this formulation.

Finally, given that there exists a continuum of sub-Kolmogorov scales from the
$n=1$ case alone, the role of higher derivatives ($n>1$) in the
$X_{n,\mathfrak{D}_{m}}$ in (2.9) is an open question and one that is relevant
to the regularity problem. Dubrulle and Gibbon [36] have shown that the lower
bound on $C(h)\geq 1-3h$ holds for all values of $n$, meaning that one can do
no better than the result given by the four-fifths law at $n=1$. In (3.5),
$h\geq-{\scriptstyle\frac{2}{3}}$ bounds $h$ away from the dangerous value of
$h=-1$, where real singularities can occur [37, 3, 35, 36]. Recent work in
[38, 39] has characterized small structures using values of $h$.

Acknowledgments : The author thanks the Isaac Newton Institute for
Mathematical Sciences, Cambridge for support and hospitality during the
programme Mathematical Aspects of Fluid Turbulence : where do we stand? in
2022, when work on this paper was undertaken. This work was supported by EPSRC
grant no EP/R014604/1&#34. Thanks are also due to Berengere Dubrulle (Saclay)
and Dario Vincenzi (Nice) for discussions. I am indebted to Nadia Padhan and
Kiran Kolluru of IISc Physics Department, Bangalore for the plots in Fig. 1.

## References

  * [1] K. Sreenivasan & G. Falkovich, Lessons from Hydrodynamic Turbulence, Physics Today, 59 (2006) 43–49.
  * [2] G. Parisi and U. Frisch, in Turbulence and Predictability in Geophysical Fluid Dynamics. Proc. Int. School of Physics “E. Fermi” 1983 Varenna, Italy, edited by M. Ghil, R. Benzi and G. Parisi, pp. 84–87 North-Holland, (Amsterdam 1985).
  * [3] U. Frisch, Turbulence : The Legacy of A. N. Kolmogorov, Cambridge University Press (Cambridge, 1995).
  * [4] C. Meneveau and K. R. Sreenivasan, The multifractal nature of turbulent energy dissipation, J. Fluid Mech. 224 (1991) 429–484.
  * [5] M. K. Verma, Energy transfer in fluid flows : Multiscale and Spectral Perspectives, Cambridge University Press (Cambridge, 2019).
  * [6] B. Dubrulle, Beyond Kolmogorov cascades, J. Fluid Mech. 867 (2019) 1–63.
  * [7] P. K. Yeung and K. Ravikumar, Advancing understanding of turbulence through extreme-scale computation : Intermittency and simulations at large problem sizes, Phys. Rev. Fluids 5 (2020) 110517.
  * [8] G. E. Elsinga, T. Ishihara and J. C. R. Hunt, Extreme dissipation and intermittency in turbulence at very high Reynolds numbers, Proc. R. Soc. A 476 (2020) 20200591.
  * [9] D. Buaria, A. Pumir, E. Bodenschatz, and P. K. Yeung, Extreme velocity gradients in turbulent flows, New J. Phys. 21 (2019) 043004.
  * [10] T. Ishihara, K. Morishita, M. Yokokawa, A. Uno, and Y. Kaneda, Energy spectrum in high-resolution direct numerical simulation of turbulence, Phys. Rev. Fluids 1 (2016) 082403(R).
  * [11] M. Buzzicotti, A. Bhatnagar, L. Biferale, A. Lanotte and S. S. Ray, Lagrangian statistics for Navier-Stokes turbulence under Fourier-mode reduction : fractal and homogeneous decimations, New J. Phys. 18 (2016) 113047.
  * [12] D. A. Donzis, J. D. Gibbon, A. Gupta, R. M. Kerr, R. Pandit and D. Vincenzi, Vorticity moments in four numerical simulations of the $3D$ Navier-Stokes equations, J. Fluid Mech. 732 (2013) 316–31.
  * [13] K. R. Sreenivasan and C. Meneveau, The Fractal Facets of Turbulence, J. Fluid Mech. 173 (1986) 357–386.
  * [14] S. Douady, Y. Couder and M. E. Brachet, Direct observation of the intermittency of intense vorticity filaments in turbulence, Phys. Rev. Lett. 67 (1991) 983–986.
  * [15] E. Aurell, U. Frisch, J. Lutsko and M. Vergassola, On the multifractal properties of the energy dissipation derived from turbulence data, J. Fluid Mech. 238 (1992) 467–486.
  * [16] C. R. Doering and C. Foias, Energy dissipation in body-forced turbulence, J. Fluid Mech. 467 (2002) 289–306.
  * [17] A. Alexakis and C. R. Doering, Energy and enstrophy dissipation in steady state 2D turbulence, Phys. Letts. A 359 (2006) 652–657.
  * [18] L. Chevillard, S. G. Roux, E. Leveque, N. Mordant, J. F. Pinton and A. Arneodo, Intermittency of velocity time increments in turbulence, Phys. Rev. Letts. 95 (2005) 064501.
  * [19] T. Ishihara, T. Gotoh and Y. Kaneda, Study of high-Reynolds number isotropic turbulence by direct numerical simulation, Annu. Rev. Fluid Mech. 41 (2009) 165–180.
  * [20] J. Leray, Sur le mouvement d’un liquide visqueux emplissant l’éspace, Acta Math. 63 (1934) 193–248.
  * [21] C. Foias, C. Guillopé and R. Temam, New a priori estimates for Navier-Stokes equations in dimension 3, Comm. Part. Diff. Equns. 6(3) (1981) 329–359.
  * [22] C. R. Doering and J. D. Gibbon, Applied Analysis of the Navier-Stokes Equations, Cambridge University Press (Cambridge, 1995).
  * [23] C. Foias, O. Manley, R. Rosa and R. Temam, Navier-Stokes Equations and Turbulence, Cambridge University Press (Cambridge, 2001).
  * [24] L. Escauriaza, G. Seregin and V. Sverak, On $L^{3,\infty}$-Solutions to the Navier-Stokes equations and backward uniqueness, Russ. Math. Surveys 58 (2003) 211–262.
  * [25] A. Cheskidov and R. Shvydkoy, A unified approach to regularity problems for the $3D$ Navier-Stokes and Euler equations : the use of Kolmogorov’s dissipation range, J. Math. Fluid Mech. 16 (2014) 263–273.
  * [26] J. C. Robinson, J. L. Rodrigo and W. Sadowski, The Three-dimensional Navier-Stokes Equations : Classical Theory ; Cambridge Studies in Advanced Mathematics, Cambridge University Press (Cambridge, 2016).
  * [27] J. D. Gibbon, Weak and Strong Solutions of the 3D Navier-Stokes Equations and Their Relation to a Chessboard of Convergent Inverse Length Scales, J. Nonlin. Sci. 29(1) (2019) 215–228.
  * [28] J. D. Gibbon, Intermittency, cascades and thin sets in three-dimensional Navier-Stokes turbulence, EPL 131 (2020) 64001.
  * [29] https://www.claymath.org/millennium-problems/navier%E2%80%93stokes-equation
  * [30] R. Benzi, G. Paladin, G. Parisi and A. Vulpiani, On the multifractal nature of fully developed turbulence and chaotic systems, J. Phys. A : Math. Gen. 17 (1984) 3521–3531.
  * [31] R. H. Kraichnan, Intermittency and attractor size in isotropic turbulence, Phys. Fluids 28 (1985) 10–11.
  * [32] G. Paladin and A. Vulpiani, Degrees of freedom of turbulence, Phys. Rev. A 35 (1987) 1971–1973.
  * [33] G. Paladin and A. Vulpiani, Anomalous scaling laws in multifractal objects, Phys. Rep. 156 4 (1987) 147–225.
  * [34] M. Nelkin, Multifractal scaling of velocity derivatives in turbulence, Phys. Rev. A 42 (1990) 7226–7229.
  * [35] G. L. Eyink (2008) Turbulence Theory – unpublished lecture notes chapters 3c,d ; The Johns Hopkins University, Baltimore, http://www.ams.jhu.edu/~eyink/Turbulence/notes/
  * [36] B. Dubrulle and J. D. Gibbon, A correspondence between the multifractal model of turbulence and the Navier–Stokes equations, Phil. Trans. R. Soc. A 380 (2022) 20210092.
  * [37] L. Caffarelli, R. Kohn and L. Nirenberg, Partial regularity of suitable weak solutions of the Navier-Stokes equations, Comm. Pure Appl. Math. 35 (1982) 771–831.
  * [38] F. Nguyen, J.-P. Laval, P. Kestener, A. Cheskidov, R. Shvydkoy and B. Dubrulle, Local estimates of Hölder exponents in turbulent vector fields, Phys. Rev. E 99 (2019) 053114.
  * [39] F. Nguyen, J.-P. Laval and B. Dubrulle, Characterizing most irregular small-scale structures in turbulence using local Hölder exponents, Phys. Rev. E 102 (2020) 063105.