# First-principles Fermi acceleration in magnetized turbulence

Martin Lemoine  Institut d’Astrophysique de Paris,
CNRS – Sorbonne Université,
98 bis boulevard Arago, F-75014 Paris, France

###### Abstract

This work provides a concrete implementation of E. Fermi’s model of particle
acceleration in magnetohydrodynamic (MHD) turbulence, connecting the rate of
energization to the gradients of the velocity of magnetic field lines, which
it characterizes within a multifractal picture of turbulence intermittency. It
then derives a transport equation in momentum space for the distribution
function. This description is shown to be substantiated by a large-scale
numerical simulation of strong MHD turbulence. The present, general framework
can be used to model particle acceleration in a variety of environments.

_Introduction–_ Particle energization through scatterings off inhomogeneous,
random moving structures is a universal process [1, 2], which has stirred
considerable interest in various branches of physics: primarily astrophysics,
with applications ranging from solar flares [3] to more remote phenomena
involving plasmas in extreme conditions, e.g. [4], but also statistical plasma
physics [5] or high energy density physics [6]. Remarkably, the two papers of
E. Fermi [1, 2] represent the first concrete scenarios for the origin of non-
thermal particles in the Universe. While the literature has placed significant
emphasis on acceleration at shock fronts, numerical experiments have
demonstrated that stochastic acceleration can be efficient [7, 8, 9], notably
so at large turbulent Alfvén velocity [10, 11, 12, 13, 14], in the sense that
it produces extended, hard powerlaw distributions of suprathermal particles.
Besides, the stochastic Fermi process assuredly plays a role in the vicinity
of shock fronts [15, 16, 17], just as it seemingly controls part of the
energization in reconnection environments [18, 19].

While the overall scenario is commonly pictured as originally formulated by E.
Fermi – a sequence of discrete, point-like interactions between a particle and
infinitely massive, perfectly conducting plasma clouds – its implementation in
a realistic turbulent context has remained a challenge [20, 21, 22], to the
extent that phenomenological applications rely on a Fokker-Planck model
parameterized by a momentum diffusion coefficient.

The present Letter proposes a novel approach to this problem and formulates a
transport equation to describe the evolution of the distribution function in
momentum space. It is first shown that particle momenta obey a continuous-time
random walk (CTRW), whose random force scales as the gradients of the velocity
of magnetic field lines, coarse-grained on a scale comparable to the particle
gyroradius $r_{\rm g}\equiv pc/eB$ ($p$ momentum, $B=|\bm{B}|$ with $\bm{B}$
magnetic field). A key observation is that those gradients are subject to
intermittency on small length scales. Hence, the random forces are neither
Gaussian, nor white noise in time and consequently, the random walk deviates
from Brownian motion, just as the transport equation, which is derived here
from known properties of CTRW, differs from Fokker-Planck. This equation is
characterized by the statistics of velocity gradients, which are captured via
a multifractal description of turbulence intermittency. This framework is
eventually shown to reproduce the time- and momentum-dependent Green functions
obtained by tracking a large number of test particles in a large-scale MHD
simulation. The present formalism thus provides a successful implementation of
stochastic Fermi acceleration in realistic, collisionless MHD turbulence.

_A (continuous-time) random walk picture–_ To evaluate energy gains/losses in
the original Fermi model, it proves convenient to boost to the scattering
center frame where the motional electric field $\bm{E}$ vanishes. The
generalization of that model to a continuous random flow similarly tracks the
particle momentum in the instantaneous (here, non-inertial) frame
$\mathcal{R}_{\not{E}}$ in which $\bm{E}$ vanishes [23, 24], which, in ideal
MHD, moves at velocity $\bm{v_{E}}=c\,\bm{E}\times\bm{B}/B^{2}$. In that
frame, momentum gains or losses scale in direct proportion to the (lab frame)
spatio-temporal gradients of the velocity field $\bm{v_{E}}$, as expressed by
$\Gamma_{\rm acc}$, $\Gamma_{\parallel}$ and $\Gamma_{\perp}$ below. In
detail, the momentum $p$ of particles with gyroradius $r_{\rm g}\ll\ell_{\rm
c}$ ($\ell_{\rm c}$ coherence length of the turbulence) evolves as

$\dot{p}\,=\,p\left\\{\Gamma_{\rm
acc}+\Gamma_{\parallel}+\Gamma_{\perp}\right\\}\,,$ (1)

with $\Gamma_{\rm acc}=-v^{-1}\,\mu\,\bm{b}\cdot\partial_{t}\bm{v_{E}}$ ($v$
particle velocity; $\mu=\bm{p}\cdot\bm{b}/p$ pitch-angle cosine with respect
to the magnetic field direction $\bm{b}=\bm{B}/B$);
$\Gamma_{\parallel}=-\,{\mu}^{2}\,\bm{b}\cdot\left(\bm{b}\cdot\bm{\nabla}\right)\bm{v_{E}}$
and
$\Gamma_{\perp}=-\,\left(1-{\mu}^{2}\right)\left[\bm{\nabla}\cdot\bm{v_{E}}-\bm{b}\cdot\left(\bm{b}\cdot\bm{\nabla}\right)\bm{v_{E}}\right]/2$.
For simplicity, the present work focuses on the sub-relativistic limit
$v_{E}\ll c$. Equation (1) – more precisely, its relativistic limit – has been
shown to account for the bulk of acceleration in numerical simulations of
collisionless turbulence [25], putting the present model on solid footing. It
generalizes the two contributions originally identified by E. Fermi:
interactions with moving magnetic mirrors are captured by $\Gamma_{\perp}$,
while orbits in dynamic, curved magnetic field lines are represented by
$\Gamma_{\parallel}$; the remaining term $\Gamma_{\rm acc}$ describes the
effective gravity force associated to acceleration/deceleration of the field
lines; of second order in $v_{E}/c$, it is subdominant in the sub-relativistic
limit, unless $v\ll v_{\rm A}$.

All quantities in Eq. (1) are understood to be coarse-grained on wavenumber
scale $k_{\rm g}\sim r_{\rm g}^{-1}$ (length scales $l_{\rm g}\sim
2\pi\,r_{\rm g}$), where $r_{\rm g}=pc/eB$ denotes the particle gyroradius in
the lab frame. This procedure filters out wavenumbers $k>k_{\rm g}$, whose
contribution averages out over a gyro-orbit, to retain the larger scales that
shape the velocity structures responsible for acceleration (in accord with the
original Fermi picture).

Henceforth, Eq. (1) is simplified to the symbolic $\dot{p}=p\,\Gamma_{l_{\rm
g}}$, $\Gamma_{l_{\rm g}}$ representing an aggregate (random) force exerted by
electromagnetic fields, coarse-grained on scale $l_{\rm g}$; order of unity
factors related to $\mu$ are thus omitted; we also consider relativistic
particles ($v\sim c$) to ease the discussion. For technical details concerning
the modeling of this random process, see [Supp. Mat. A]. Separate now
fluctuations from the mean,
$\Gamma_{l}=\langle\Gamma_{l}\rangle+\delta\Gamma_{l}$, the average carrying
over the statistical realizations of the turbulent flow:
$\langle\Gamma_{l}\rangle$ characterizes systematic heating, while the random
$\delta\Gamma_{l}$ represents the diffusive part. If $\delta\Gamma_{l}$ were
Gaussian distributed, and its time correlation function that of white noise,
the process would describe Brownian motion, in one-to-one correspondence with
a Fokker-Planck equation for the distribution function [26]. As anticipated
above, however, those random forces are neither Gaussian in amplitude, nor
white noise in time: at small scales, they develop large powerlaw tails as a
result of intermittency, while at large scales, the coherence time of the
random force $\gtrsim l_{\rm g}/c$ cannot be regarded as infinitesimal.

To obtain the transport equation, we first observe that the process
$\dot{p}=p\,\Gamma_{l_{\rm g}}$ can be described as a CTRW: unlike discretized
Brownian motion, which operates at a fixed and uniform time step, the random
walk is here defined by the joint probability
$\phi\left(p|p^{\prime};\,t-t^{\prime}\right)$ to jump from $p^{\prime}$ to
$p$ in time $\Delta t=t-t^{\prime}$, with both $\Delta p=p-p^{\prime}$ and
$\Delta t$ regarded as random variables. Expectations are $\Delta t\sim l_{\rm
g}(p^{\prime})/c$ – thus, a function of $p^{\prime}$ – and $\Delta\ln
p\sim\Gamma_{l_{\rm g}}\Delta t$. We will assume $\Delta t$ to be
exponentially distributed with mean waiting time $t_{p}\equiv l_{\rm
g}(p^{\prime})/c$ and $\Delta\ln\,p$ distributed as $\Gamma_{l_{\rm g}}l_{\rm
g}/c$ [$l_{\rm g}=l_{\rm g}(p^{\prime})$], see [Supp. Mat. A] for methodology.
The random walk is then entirely defined by the statistics of the velocity
gradients $\Gamma_{l_{\rm g}}$.

_Statistics of momentum jumps–_ In turbulence theories, such statistics are
conveniently described within a multifractal analysis [27, 28, 29], which
ascribes to each position $\bm{x}$ a local scaling exponent $h(\bm{x})$ for
gradients on coarse-graining scale $l$, viz.

$\Gamma_{l}(\bm{x})\,\sim\,\Gamma_{\ell_{\rm c}}(\bm{x})\,\left(l/\ell_{\rm
c}\right)^{h(\bm{x})}\,,$ (2)

and which describes the set of locations $\bm{x}$ with index $h(\bm{x})$ as a
fractal of dimension $d(h)$. The statistics of $\Gamma_{l}$ are thus entirely
captured by the probability distribution function (p.d.f.)
${\textsf{p}}_{\Gamma_{\ell_{\rm c}}}$ of $\Gamma_{\ell_{\rm c}}$ and by the
spectrum $d(h)$, since the probability of being at $\bm{x}$ in a set with
exponent $h$ on scale $l$ evolves as the volume filling fraction $l^{D-d(h)}$
($D$ number of spatial dimensions). The gradient $\Gamma_{\ell_{\rm
c}}(\bm{x})$ on the coherence scale $\ell_{\rm c}$ is naturally modeled as a
Gaussian variable with standard deviation $\sigma_{\rm c}\sim v_{\rm
A}/\ell_{\rm c}$, where $v_{\rm A}$ denotes the Alfvén velocity of the
turbulent component. The spectrum $d(h)$ can take different forms, the
simplest being log-normal [30], modern descriptions of the statistics of
Elsässer fields in MHD turbulence rather relying on log-Poisson models [31,
32, 33, 34]. We use the former log-normal form, as it provides a simple and
satisfactory description of the statistics of the gradients of $\bm{v_{E}}$;
see [Supp. Mat. B], which includes Refs. [35, 36, 37, 38]. We thus derive the
p.d.f. ${\textsf{p}}_{\Gamma_{l}}$ of $\Gamma_{l}$ as (using $D=3$) [29]

$\displaystyle{\textsf{p}}_{\Gamma_{l}}$ $\displaystyle\,\sim\,\int{\rm
d}\Gamma_{\ell_{\rm c}}\,{\textsf{p}}_{\Gamma_{\ell_{\rm c}}}\,\int{\rm
d}h\,l^{3-d(h)}\,\delta\left[\Gamma_{l}-\Gamma_{\ell_{\rm
c}}\,\left(l/\ell_{\rm c}\right)^{h}\right]\,.$

This p.d.f. remains to be properly normalized. In this formulation, the
gradient statistics ${\textsf{p}}_{\Gamma_{l}}$ on all scales reduce to a
function of the main quantities $v_{\rm A}$, $\ell_{\rm c}$ and the few
parameters characterizing $d(h)$, which themselves depend on the physical
properties of the turbulence. This offers a first-principles connection
between the fundamental statistics of turbulence intermittency and the physics
of particle energization.

_The transport equation–_ The CTRW is exactly equivalent to the following
kinetic equation for the volume averaged distribution $n_{p}(t)=4\pi
p^{2}f(p,\,t)$, where $f(p,\,t)$ represents the angle-averaged distribution
function [39, 40]:

$\displaystyle\partial_{t}\,n_{p}(t)\,=\,\int_{0}^{+\infty}{\rm
d}p^{\prime}\int_{0}^{t}{\rm d}t^{\prime}\,\biggl{[}$
$\displaystyle\psi\left(p|p^{\prime};\,t-t^{\prime}\right)n_{p^{\prime}}(t^{\prime})$
$\displaystyle-\psi\left(p^{\prime}|p;\,t-t^{\prime}\right)n_{p}(t^{\prime})\biggr{]}\,.$
(4)

The kernel $\psi\left(p|p^{\prime};\,t-t^{\prime}\right)$ differs from the
CTRW jump distribution probability
$\phi\left(p|p^{\prime};\,t-t^{\prime}\right)$ introduced earlier, yet the two
are related as follows. Denoting with a tilde symbol the Laplace transform in
time, and $\nu$ the Laplace variable conjugate to $t-t^{\prime}$,

$\tilde{\psi}\left(p|p^{\prime};\,\nu\right)\,=\,\frac{\nu\,\tilde{\phi}\left(p|p^{\prime};\,\nu\right)}{1-\tilde{\phi}_{p^{\prime}}\left(\nu\right)}\,,$
(5)

with the short-hand notation
$\tilde{\phi}_{p^{\prime}}\left(\nu\right)\equiv\int_{0}^{+\infty}{\rm
d}p\,\tilde{\phi}\left(p|p^{\prime};\,\nu\right)$, the subscript
${}_{p^{\prime}}$ emphasizing the dependence on $p^{\prime}$. As discussed
above, we characterize the CTRW with a joint probability distribution of the
form:

$\phi\left(p|p^{\prime};\,t-t^{\prime}\right)\,=\,\varphi\left(p|p^{\prime}\right)\,\frac{e^{-(t-t^{\prime})/t_{p^{\prime}}}}{t_{p^{\prime}}}\,,$
(6)

recalling that $t_{p^{\prime}}\equiv l_{\rm g}(p^{\prime})/c$. Then,

$\tilde{\psi}\left(p|p^{\prime};\,\nu\right)\,=\,\frac{\varphi\left(p|p^{\prime}\right)}{t_{p^{\prime}}}\,,$
(7)

in which case the transport equation becomes local in time [40],

$\partial_{t}\,n_{p}(t)\,=\,\int_{0}^{+\infty}{\rm
d}p^{\prime}\,\left[\frac{\varphi\left(p|p^{\prime}\right)}{t_{p^{\prime}}}n_{p^{\prime}}(t)-\frac{\varphi\left(p^{\prime}|p\right)}{t_{p}}n_{p}(t)\right]\,.$
(8)

It takes the form of a master equation for a Markov process, balancing gains
and losses at respective rates $t_{p^{\prime}}$ and $t_{p}$. As
$\varphi\left(p|p^{\prime}\right)$ represents the p.d.f. to jump to $p$ from
$p^{\prime}$ in any given amount of time, it is normalized through
$\int_{0}^{+\infty}{\rm d}p\,\varphi\left(p|p^{\prime}\right)=1$. Recalling
that $\Delta\ln p$ is distributed as $\Gamma_{l_{\rm g}}l_{\rm g}/c$, the
p.d.f. $\varphi(p|p^{\prime})$ derives from ${\textsf{p}}_{\Gamma_{l}}$
through

$\varphi\left(p|p^{\prime}\right)=\frac{1}{p\,t_{p^{\prime}}}{\textsf{p}}_{\Gamma_{l_{\rm
g}}}\,,$ (9)

at $\Gamma_{l_{\rm g}}=\ln(p/p^{\prime})/t_{p^{\prime}}$, with $l_{\rm
g}=l_{\rm g}(p^{\prime})$.

_Test against numerical experiments–_

Figure 1: Statistics $|\Gamma_{l}|\,{\textsf{p}}_{|\Gamma_{l}|}$ of the
absolute values of the gradients $\Gamma_{l}$ ($\Gamma_{l}$ expressed in units
of $c/\ell_{\rm c}$). Symbols: values recorded in a $1\,024^{3}$ MHD
simulation at various coarse-graining scales $l$, as indicated; squares (resp.
diamonds): gradients measured along (resp. perpendicular to) the mean magnetic
field direction as coarse-grained on scale $l$. The p.d.f. reveals powerlaw
tails at large values of $|\Gamma_{l}|$ on small scales. Solid red line:
adjustment of a multi-fractal lognormal model; solid blue line: a broken
powerlaw approximation; dotted line: shot noise level associated to the finite
sampling variance.

The above model is now tested on a direct numerical simulation of driven
incompressible MHD turbulence ($1\,024^{3}$ with $1\,024$ time steps) [41,
42]. This simulation has no guide field, and its units have been set to obtain
an Alfvén velocity $v_{\rm A}=0.4\,c$; this offers a reasonable compromise
between the limits of applicability of this (sub-relativistic) simulation and
the value of $v_{\rm A}$ needed to observe acceleration within the duration of
the simulation ($2.8\,\ell_{\rm c}/c$); see [Supp. Mat. B] for the methods
used to extract numerical data from this simulation.

The statistics of the absolute values of the velocity gradients $\Gamma_{l}$,
parallel and perpendicular to the mean magnetic field and coarse-grained on
scale $l$, are shown in Fig. 1 for various values of $l$. The red solid line
represents the adjustment obtained for a log-normal spectrum
$d(h)=3-\left(h-h_{\textsc{mf}}\right)^{2}/\left(2\sigma_{\textsc{mf}}^{2}\right)$,
parameterized by $h_{\textsc{mf}}=-0.2$ and $\sigma_{\textsc{mf}}=0.9$, see
[Supp. Mat. B]. As the random walk has been simplified to one aggregate force
term $\Gamma_{l}$, we have chosen to tune those parameters to provide a fair
reconstruction of the ensemble of parallel and perpendicular gradients, rather
than fitting one or the other. The solid blue line shows an adjustment of
${\textsf{p}}_{|\Gamma_{l}|}$ by a broken powerlaw approximation, which has
the advantage of speeding up the numerical integration of the kinetic
transport equation. It takes the form

${\textsf{p}}_{\Gamma_{l}}\,\sim\,\left[1+\left(\frac{\Gamma_{l}}{\sigma_{\textsc{bp}}(l)}\right)^{k_{0}(l)/k_{1}}\right]^{-k_{1}}\,,$
(10)

with $\sigma_{\textsc{bp}}(l)$ and $k_{0}(l)$ $l-$dependent quantities;
$\sigma_{\textsc{bp}}(l)$ characterizes the width of the distribution before
it turns over into the powerlaw behavior with index $\simeq-k_{0}(l)$, while
$k_{1}=3$ ensures the smoothness of that transition from core to wing, see
[Supp. Mat. B] for details and methodology. Both models use a width
$\sigma\simeq 0.3\,v_{\rm A}/\ell_{\rm c}$, with $\sigma=\sigma_{\rm c}$ for
$\Gamma_{\ell_{\rm c}}$ [resp. $\sigma=\sigma_{\textsc{bp}}(\ell_{\rm c})$]
for the multifractal (resp. broken powerlaw) model.

From ${\textsf{p}}_{\Gamma_{l}}$, we derive $\varphi\left(p|p^{\prime}\right)$
using Eq. (9) then integrate the transport equation Eq. (8) to compare the
theoretical spectra $n_{p}(t)$ with experimental ones obtained by tracking a
large number of particles in the MHD simulation. We remark here that, first,
those two models for ${\textsf{p}}_{\Gamma_{l}}$ provide slightly different
fits to the measured gradient statistics, therefore the comparison of the
corresponding $n_{p}(t)$ provides a first glance at how the choice of
parameters affects those spectra. Second, ${\textsf{p}}_{\Gamma_{l}}$ extends
to values of $\Gamma_{l}$ beyond the range where it can be measured in the
simulation, allowing in particular for unbounded gains in momentum over a
finite timescale $l_{\rm g}/c$; to regularize this, we multiply
$\varphi\left(p|p^{\prime}\right)$ by a cut-off $\exp\left[-\left(\ln p-\ln
p^{\prime}\right)^{2}\right]$, which bounds the maximum gain to the order of
unity on that timescale. This theoretical maximum is well motivated here,
since $v_{\rm A}$ is not small compared to $c$ [24]; different choices are
possible, but the overall influence does not exceed that associated to the
uncertainty affecting ${\textsf{p}}_{\Gamma_{l}}$. Third, the model predicts
the evolution of momentum in the frame $\mathcal{R}_{\not{E}}$, hence the
comparison to the experimental $n_{p}(t)$ requires a boost to the simulation
frame, which slightly broadens the particle distribution. To minimize this
effect, we inject particles at a given momentum $p_{0}$ in the
$\mathcal{R}_{\not{E}}$ frame at time $t=0$, integrate the equation then boost
the theoretical spectrum to the simulation frame. We thus effectively measure
a Green function convoluted with this boost. Finally, the numerical
distributions of the velocity gradients have small yet non-zero mean values,
implying a slight advection drift toward increasing momenta; it has been taken
into account in adjusting ${\textsf{p}}_{\Gamma_{l}}$ to the data. Further
details on these procedures are provided in [Supp. Mat.B].

Figure 2: Particle spectra obtained at the final time $t=2.8\ell_{\rm c}/c$
for various initial momenta as indicated in units of $p_{\rm c}$; $p_{\rm c}$
is such that $l_{\rm g}(p_{\rm c})=\ell_{\rm c}$. The red and blue lines
correspond to the two models shown in Fig. 1.  Figure 3: Same as Fig. 2, now
showing the particle spectra at various times, for the smallest initial
injection momentum, as indicated.

Those theoretical models (red and blue curves, following the conventions of
Fig. 1) are compared in Figs. 2 and 3 to the momentum- and time-dependent
Green functions obtained by tracking a large number of test particles in the
simulation cube. The test particles have been propagated in the time-dependent
snapshots of the simulation, i.e. the time evolution of the electromagnetic
fields has been properly taken into account. In Fig. 2, the spectra are
plotted vs $p/p_{0}$ for increasing values of the initial momentum $p_{0}$,
from $p_{0}/p_{\rm c}=0.067$ to $p_{0}/p_{\rm c}=2.1$, thus covering a dynamic
range of $1.5$ decade; $p_{\rm c}$ denotes the momentum such that $l_{\rm
g}(p_{\rm c})=\ell_{\rm c}$, i.e. the coarse-graining becomes comparable to
the coherence scale. At $p_{0}>p_{\rm c}$, the spectrum narrows down: the
acceleration decreases because the particle then sees the turbulence as a
collection of incoherent cells of small extent relative to its gyroradius. In
Fig. 3, the spectra are plotted for $p_{0}/p_{\rm c}=0.067$ up to the maximal
integration time of the simulation, $t=2.8\,\ell_{\rm c}/c$. This rather
satisfactory comparison between theoretical and experimental Green functions
supports the present picture.

_Discussion and perspectives–_ As in the original Fermi picture, the
efficiency of acceleration scales with the turbulent Alfvén velocity $v_{\rm
A}$ 111In the presence of a mean field $B_{0}$ (corresponding Alfvén velocity
$v_{\rm A\,0}$), $v_{\rm A}=(\delta B/B_{0})v_{\rm A\,0}$ ($\delta B$
turbulent component), therefore $\delta B/B_{0}$ impacts acceleration at fixed
$v_{\rm A\,0}$., which controls the magnitude of the random force through its
influence on $\Gamma_{\ell_{\rm c}}$. At lower $v_{\rm A}$, one observes a
softer spectrum at a given time, but as time passes, the spectrum becomes
harder; two different values of $v_{\rm A}$ eventually give comparable spectra
at times rescaled by $1/v_{\rm A}^{2}$. Figures 2 and 3 indicate that powerlaw
spectra are generic, as observed in kinetic simulations [10, 11]. Those
powerlaws find their origin in the large excursions that particles momenta can
undergo in sparse regions of intense gradients [24, 44], here captured by the
extended wings of $\varphi\left(p^{\prime}|p\right)$. In the absence of those,
a Kramers-Moyal expansion would indeed reduce Eq. (8) to the Fokker-Planck
form with diffusion coefficient $D_{pp}=\int{\rm
d}p^{\prime}\,\left(p^{\prime}-p\right)^{2}\,\varphi\left(p^{\prime}|p\right)/t_{p}$
and advection coefficient $A_{p}=\int{\rm
d}p^{\prime}\,\left(p^{\prime}-p\right)\,\varphi\left(p^{\prime}|p\right)/t_{p}$.

The present scenario naturally accounts for the recent observation that some
particles can see their energy increase exponentially fast in localized
regions, up to an energy gain of a few [9], by virtue of
$\dot{p}=p\,\Gamma_{l}$ with $\Gamma_{l}$ varying on scales of extent $l$; see
also [24, 45]. As scattering is here dominated by magnetic mirrors, which
become intermittent on small scales, one also anticipates anomalous spatial
transport. Inspection of numerical data confirms that some particles preserve
their pitch-angles over long distances $\gtrsim\ell_{\rm c}$, while others
suffer strong deflection on short distances. This might account for spatial
superdiffusion events observed in recent simulations [46] and trapping in
others [8].

The present approach differs from that of Ref. [13], which extracts from
numerical experiments empirical momentum-dependent functions
$A_{p}=\langle\Delta p\rangle/\Delta t$ and $D_{pp}=\langle\Delta
p^{2}\rangle/\Delta t$. It also differs from Ref. [47], which models transport
in momentum space through a fractional Fokker-Planck equation, describing the
random walk as a Lévy process. The present $\varphi(p^{\prime}|p)$ is more
akin to a truncated Lévy flight, which contains flat tails but well-defined
high-order moments. This implies that, at asymptotic times, the random walk
will behave as some Brownian motion, due to central-limit convergence;
however, the extended wings of $\varphi\left(p^{\prime}|p\right)$ slow down
this convergence quite appreciably, as expressed by the Berry-Esséen theorem
[48], hence, for practical matters, the present non-Fokker-Planck form remains
required. This transport equation remains amenable to extensions: for
instance, wave-particle resonant interactions, if effective [49, 50], could be
included by adding an extra diffusion term; similarly for particle heating in
small-scale non-ideal electric fields, which ensures injection into the Fermi
process in kinetic simulations [10, 11, 12, 13, 14]; finally, including
standard radiative and escape terms would allow to model the emerging spectra
from astrophysical sources.

In summary, the present work has established a novel, general framework for
implementing stochastic Fermi acceleration in a realistic, collisionless
turbulent bath, opening a connection between the statistics of the
intermittent gradients of the velocity of magnetic field lines and the rate of
energization. Improved insight on turbulence intermittency could ultimately
allow a first-principles calculation of particle spectra. The applicability of
this formalism thus extends beyond that of the numerical simulation of
incompressible MHD turbulence against which it was successfully tested.

_Acknowledgments:_ This work has been supported by the Sorbonne Université
DIWINE Emergence-2019 program and by the ANR (UnRIP project, Grant No.
ANR-20-CE30-0030). The possibility to use the resources of the JH Turbulence
Database (JHTDB), which is developed as an open resource by the Johns Hopkins
University, under the sponsorship of the National Science Foundation, is
gratefully acknowledged.

## References

  * Fermi [1949] E. Fermi, On the Origin of the cosmic radiation, Phys. Rev. 75, 1169 (1949).
  * Fermi [1954] E. Fermi, Galactic magnetic fields and the origin of cosmic radiation., Astrophys. J. 119, 1 (1954).
  * Petrosian [2012] V. Petrosian, Stochastic Acceleration by Turbulence, Sp. Sc. Rev. 173, 535 (2012).
  * Tramacere  _et al._ [2011] A. Tramacere, E. Massaro, and A. M. Taylor, Stochastic Acceleration and the Evolution of Spectral Distributions in Synchro-Self-Compton Sources: A Self-consistent Modeling of Blazars’ Flares, Astrophys. J. 739, 66 (2011).
  * Zaslavskiĭ and Chirikov [1972] G. M. Zaslavskiĭ and B. V. Chirikov, Reviews of Topical Problems: Stochastic Instability of Non-Linear Oscillations, Soviet Physics Uspekhi 14, 549 (1972).
  * Beyer  _et al._ [2018] K. A. Beyer, B. Reville, A. F. A. Bott, H.-S. Park, S. Sarkar, and G. Gregori, Analytical estimates of proton acceleration in laser-produced turbulent plasmas, Journal of Plasma Physics 84, 905840608 (2018).
  * Dmitruk  _et al._ [2003] P. Dmitruk, W. H. Matthaeus, N. Seenu, and M. R. Brown, Test Particle Acceleration in Three-dimensional Magnetohydrodynamic Turbulence, Astrophys. J. 597, L81 (2003).
  * Trotta  _et al._ [2020] D. Trotta, L. Franci, D. Burgess, and P. Hellinger, Fast Acceleration of Transrelativistic Electrons in Astrophysical Turbulence, Astrophys. J. 894, 136 (2020).
  * Pezzi  _et al._ [2022] O. Pezzi, P. Blasi, and W. H. Matthaeus, Relativistic Particle Transport and Acceleration in Structured Plasma Turbulence, Astrophys. J. 928, 25 (2022).
  * Zhdankin  _et al._ [2017] V. Zhdankin, G. R. Werner, D. A. Uzdensky, and M. C. Begelman, Kinetic turbulence in relativistic plasma: From thermal bath to nonthermal continuum, Phys. Rev. Lett. 118, 055103 (2017).
  * Comisso and Sironi [2018] L. Comisso and L. Sironi, Particle acceleration in relativistic plasma turbulence, Phys. Rev. Lett. 121, 255101 (2018).
  * Comisso and Sironi [2019] L. Comisso and L. Sironi, The interplay of magnetically dominated turbulence and magnetic reconnection in producing nonthermal particles, Astrophys. J. 886, 122 (2019).
  * Wong  _et al._ [2020] K. Wong, V. Zhdankin, D. A. Uzdensky, G. R. Werner, and M. C. Begelman, First-principles demonstration of diffusive-advective particle acceleration in kinetic simulations of relativistic plasma turbulence, Astrophys. J. 893, L7 (2020).
  * Zhdankin [2021] V. Zhdankin, Particle Energization in Relativistic Plasma Turbulence: Solenoidal versus Compressive Driving, Astrophys. J. 922, 172 (2021).
  * Matsumoto  _et al._ [2017] Y. Matsumoto, T. Amano, T. N. Kato, and M. Hoshino, Electron Surfing and Drift Accelerations in a Weibel-Dominated High-Mach-Number Shock, Phys. Rev. Lett. 119, 105101 (2017).
  * Trotta  _et al._ [2021] D. Trotta, F. Valentini, D. Burgess, and S. Servidio, Phase space transport in the interaction between shocks and plasma turbulence, Proceedings of the National Academy of Science 118, e2026764118 (2021).
  * Nakanotani  _et al._ [2022] M. Nakanotani, G. P. Zank, and L. L. Zhao, Turbulence-dominated Shock Waves: 2D Hybrid Kinetic Simulations, Astrophys. J. 926, 109 (2022).
  * Drake  _et al._ [2006] J. F. Drake, M. Swisdak, H. Che, and M. A. Shay, Electron acceleration from contracting magnetic islands during reconnection, Nature 443, 553 (2006).
  * Guo  _et al._ [2019] F. Guo, X. Li, W. Daughton, P. Kilian, H. Li, Y.-H. Liu, W. Yan, and D. Ma, Determining the Dominant Acceleration Mechanism during Relativistic Magnetic Reconnection in Large-scale Systems, Astrophys. J. 879, L23 (2019).
  * Bykov and Toptygin [1993] A. M. Bykov and I. Toptygin, Particle kinetics in highly turbulent plasmas (renormalization and self-consistent field methods), Physics Uspekhi 36, 1020 (1993).
  * Bouchet  _et al._ [2004] F. Bouchet, F. Cecconi, and A. Vulpiani, Minimal stochastic model for Fermi’s acceleration, Phys. Rev. Lett. 92, 040601 (2004).
  * Burby  _et al._ [2013] J. W. Burby, A. I. Zhmoginov, and H. Qin, Hamiltonian mechanics of stochastic acceleration, Phys. Rev. Lett. 111, 195001 (2013).
  * Lemoine [2019] M. Lemoine, Generalized Fermi acceleration, Phys. Rev. D 99, 083006 (2019).
  * Lemoine [2021] M. Lemoine, Particle acceleration in strong MHD turbulence, Phys. Rev. D 104, 063020 (2021).
  * Bresci  _et al._ [2022] V. Bresci, M. Lemoine, L. Gremillet, L. Comisso, L. Sironi, and C. Demidem, Nonresonant particle acceleration in strong turbulence: Comparison to kinetic and MHD simulations, Phys. Rev. D 106, 023028 (2022).
  * Risken [1989] H. Risken,  _The Fokker-Planck Equation. Methods of Solution and Applications_ (Springer, Berlin, New-York, 1989).
  * Frisch and Parisi [1985] U. Frisch and G. Parisi, On the singularity structure of fully developed turbulence, in  _Turbulence and Predictability in Geophysical Fluid Dynamics_ , edited by M. Ghil, R. Benzi, and G. Parisi (North Holland, Amsterdam, 1985) pp. 84–87.
  * Benzi  _et al._ [1984] R. Benzi, G. Paladin, A. Vulpiani, and G. Parisi, On the multifractal nature of fully developed turbulence and chaotic systems, Journal of Physics A Mathematical General 17, 3521 (1984).
  * Benzi  _et al._ [1991] R. Benzi, L. Biferale, G. Paladin, A. Vulpiani, and M. Vergassola, Multifractality in the statistics of the velocity gradients in turbulence, Phys. Rev. Lett. 67, 2299 (1991).
  * Dubrulle [2019] B. Dubrulle, Beyond Kolmogorov cascades, Journal of Fluid Mechanics 867, P1 (2019).
  * Biskamp and Müller [2000] D. Biskamp and W.-C. Müller, Scaling properties of three-dimensional isotropic magnetohydrodynamic turbulence, Physics of Plasmas 7, 4889 (2000).
  * Müller  _et al._ [2003] W.-C. Müller, D. Biskamp, and R. Grappin, Statistical anisotropy of magnetohydrodynamic turbulence, Phys. Rev. E 67, 066302 (2003).
  * Padoan  _et al._ [2004] P. Padoan, R. Jimenez, Å. Nordlund, and S. Boldyrev, Structure Function Scaling in Compressible Super-Alfvénic MHD Turbulence, Phys. Rev. Lett. 92, 191102 (2004).
  * Chandran  _et al._ [2015] B. D. G. Chandran, A. A. Schekochihin, and A. Mallet, Intermittency and Alignment in Strong RMHD Turbulence, Astrophys. J. 807, 39 (2015).
  * She and Leveque [1994] Z.-S. She and E. Leveque, Universal scaling laws in fully developed turbulence, Phys. Rev. Lett. 72, 336 (1994).
  * Dubrulle [1994] B. Dubrulle, Intermittency in fully developed turbulence: Log-Poisson statistics and generalized scale covariance, Phys. Rev. Lett. 73, 959 (1994).
  * Yang  _et al._ [2019] Y. Yang, M. Wan, W. H. Matthaeus, Y. Shi, T. N. Parashar, Q. Lu, and S. Chen, Role of magnetic field curvature in magnetohydrodynamic turbulence, Phys. Plasmas 26, 072306 (2019).
  * Yuen and Lazarian [2020] K. H. Yuen and A. Lazarian, Curvature of magnetic field lines in compressible magnetized turbulence: Statistics, magnetization predictions, gradient curvature, modes, and self-gravitating media, Astrophys. J. 898, 66 (2020).
  * Klafter and Silbey [1980] J. Klafter and R. Silbey, Derivation of the Continuous-Time Random-Walk Equation, Phys. Rev. Lett. 44, 55 (1980).
  * Berkowitz  _et al._ [2006] B. Berkowitz, A. Cortis, M. Dentz, and H. Scher, Modeling non-Fickian transport in geological formations as a continuous time random walk, Reviews of Geophysics 44, RG2003 (2006).
  * Li  _et al._ [2008] Y. Li, E. Perlman, M. Wan, Y. Yang, C. Meneveau, R. Burns, S. Chen, A. Szalay, and G. Eyink, A public turbulence database cluster and applications to study Lagrangian evolution of velocity increments in turbulence, J. Turbul. 9, N31 (2008).
  * Eyink  _et al._ [2013] G. Eyink, E. Vishniac, C. Lalescu, H. Aluie, K. Kanov, K. Bürger, R. Burns, C. Meneveau, and A. Szalay, Flux-freezing breakdown in high-conductivity magnetohydrodynamic turbulence, Nature (London) 497, 466 (2013).
  * Note [1] In the presence of a mean field $B_{0}$ (corresponding Alfvén velocity $v_{\rm A\tmspace+{.1667em}0}$), $v_{\rm A}=(\delta B/B_{0})v_{\rm A\tmspace+{.1667em}0}$ ($\delta B$ turbulent component), therefore $\delta B/B_{0}$ impacts acceleration at fixed $v_{\rm A\tmspace+{.1667em}0}$.
  * Lemoine and Malkov [2020] M. Lemoine and M. A. Malkov, Power-law spectra from stochastic acceleration, Mon. Not. R. Astron. Soc. 499, 4972 (2020).
  * Brunetti and Lazarian [2016] G. Brunetti and A. Lazarian, Stochastic reacceleration of relativistic electrons by turbulent reconnection: A mechanism for cluster-scale radio emission?, Mon. Not. R. Astron. Soc. 458, 2584 (2016).
  * Maiti  _et al._ [2022] S. Maiti, K. Makwana, H. Zhang, and H. Yan, Cosmic-ray Transport in Magnetohydrodynamic Turbulence, Astrophys. J. 926, 94 (2022).
  * Isliker  _et al._ [2017] H. Isliker, L. Vlahos, and D. Constantinescu, Fractional transport in strongly turbulent plasmas, Phys. Rev. Lett. 119, 045101 (2017).
  * Shlesinger [1995] M. F. Shlesinger, Comment on “Stochastic Process with Ultraslow Convergence to a Gaussian: The Truncated Lévy Flight”, Phys. Rev. Lett. 74, 4959 (1995).
  * Chandran [2000] B. D. G. Chandran, Scattering of energetic particles by anisotropic magnetohydrodynamic turbulence with a Goldreich-Sridhar power spectrum, Phys. Rev. Lett. 85, 4656 (2000).
  * Yan and Lazarian [2002] H. Yan and A. Lazarian, Scattering of cosmic rays by magnetohydrodynamic interstellar turbulence, Phys. Rev. Lett. 89, 281102 (2002).

Supplemental Material to “First-principles Fermi acceleration in magnetized
turbulence”

##  Appendix A Methods – stochastic model

The random walk is written: $\dot{p}=p\,\Gamma_{l_{\rm g}}$, where
$\Gamma_{l_{\rm g}}$ is expressed as a gradient of the velocity field
$\bm{v_{E}}$ coarse-grained on scale $l$ [Eq. (1)]. If $\Gamma_{l_{\rm g}}$
were Gaussian distributed, with a vanishing coherence time, meaning a p.d.f.
characterized by $\langle\Gamma_{l_{\rm g}}\rangle$ (drift) and
$\left\langle\delta\Gamma_{l_{\rm g}}(t)\delta\Gamma_{l_{\rm
g}}(t^{\prime})\right\rangle\propto\delta(t-t^{\prime})$
($\delta\Gamma_{l_{\rm g}}\equiv\Gamma_{l_{\rm g}}-\langle\Gamma_{l_{\rm
g}}\rangle$) [26], the random walk would be described as Brownian motion and
the transport equation would take a Fokker-Planck form. This is not the case
here, since the random process describes the influence of velocity structures
whose scales $l$ enter the range $l_{\rm g}<l<\ell_{\rm c}$, and whose
amplitudes are not Gaussian distributed. One must thus consider the full
transition kernel $\phi\left(p|p^{\prime};\,t-t^{\prime}\right)$ that defines
the probability to transit from momentum $p^{\prime}$ to $p$, from time
$t^{\prime}$ to $t$. The stochastic process then takes the form of a
continuous-time random walk: the time step between two jumps becomes random
distributed and its distribution depends on the location in momentum space.
This description fits well the present picture, because the coarse-graining
scale $l_{\rm g}(p^{\prime})$ depends on momentum. In order to simplify the
description, and to bring the transport equation to a tractable form, two
assumptions are made: (1) structures on the coarse-graining scale $l_{\rm g}$
dominate the acceleration process; (2) the p.d.f. of the time interval $\Delta
t$ between two momentum jumps is exponentially distributed, meaning ${\rm
p.d.f.}(\Delta t)=\exp(-\Delta t/\tau)/\tau$, with $\tau=t_{p^{\prime}}\equiv
l_{\rm g}/c$ by application of (1) (for particles moving at $v\sim c$).

Assumption (2) is weak, in the sense that a Monte Carlo implementation of the
above stochastic process indicates that other p.d.f. for the transition time,
with equal mean waiting time, lead to nearly identical particle spectra. The
choice of an exponential p.d.f. simplifies however the derivation of the
transport equation.

Assumption (1) rests on the premise that the stochastic force increases in
amplitude as one goes to smaller scales, and so does intermittency. This also
agrees with the discussion of the comparative influence of various scales $l$
on energization in [24]. In order to test this assumption, the transport
equation has been integrated to obtain particle spectra for different choices
of the dominant scale $l$. Namely, the mean waiting time $\tau$ between two
jumps has been set to $\tau={\rm max}\left(l/c,\,l_{\rm g}/c\right)$, while
$\Delta\ln p$ has been assumed to scale as $\Gamma_{c\,\tau}\,\tau$ as before.
Setting $l=l_{\rm g}$ for all $p^{\prime}$, one recovers $\tau=t_{p^{\prime}}$
and the transport equation discussed in the main text; by contrast, fixing $l$
to some constant value mimics a situation in which structures on scale $l$
dominate the physics of acceleration. This exercise reveals that the spectra
become steeper for $l\gtrsim 0.1\ell_{\rm c}$, indicating that the largest
scales provide a sub-dominant contribution to acceleration. This falls in line
with the above expectations: as intermittency becomes weak on large scales,
the development of powerlaw tails is suppressed. For the values of $p$
considered here in the comparison to numerical data, i.e. $l_{\rm
g}(p_{0})/\ell_{\rm c}\gtrsim 0.07$, this observation supports assumption (1).

For momenta $p\geq p_{\rm c}$, such that $l_{\rm g}\geq l_{\rm c}$, the mean
waiting time must be kept constant at $\ell_{\rm c}/c$, as this properly
describes the coherence time of the random force that is exerted on the
particle. On such timescale, the particles suffers energy gain or loss of the
order of $|\Delta\ln p|\simeq(v_{\rm A}/c)\ell_{\rm c}/l_{\rm
g}\sim\Gamma_{\ell_{\rm c}}\,\ell_{\rm c}^{2}/l_{\rm g}c$, the additional
factor $\ell_{\rm c}/l_{\rm g}$ accounting for the fact that those particles
decouple from the turbulence due to their large magnetic rigidity.

Finally, the stochastic process can be generalized to sub-relativistic
particles ($v\ll c$) through the direct substitution $t_{p}\rightarrow l_{\rm
g}/v$ for the mean waiting time, in which case $\Delta\ln p$ is distributed as
$\Gamma_{l_{\rm g}}l_{\rm g}/v$. The definition of $\Gamma_{\rm acc}$
indicates that it scales as $v_{\rm A}/v$ relatively to the other two; hence,
at $v\ll v_{\rm A}$, this term should dominate the energization. As the
intermittent statistics of $\Gamma_{\rm acc}$ remain comparable to those of
$\Gamma_{\parallel}$ and $\Gamma_{\perp}$, this should not affect the particle
spectra significantly.

##  Appendix B Methods – data extraction and manipulation out of a $1024^{3}$
MHD numerical simulation

The numerical data is extracted from a direct numerical simulation of
incompressible MHD on a periodic grid containing $1\,024^{3}$ cells, sampled
on $1\,024$ time steps. This simulation, which is made available for public
use from the Johns Hopkins University Turbulence Database, available from:
http://turbulence.pha.jhu.edu/Forced_MHD_turbulence.aspx. [42], is visco-
resistive, with magnetic Reynolds number at the Taylor scale $\simeq 140$, and
magnetic Prandtl number unity. The $1\,024$ cube snapshots have been archived
after the simulation has reached a statistically stationary state; the
properties of the turbulence evolve little over those additional $1\,024$ time
steps. The turbulence is driven through an external random force acting on
velocity fluctuations on a stirring scale equal to $0.5$ in units of the size
of the simulation box. The corresponding correlation length for the magnetic
field, as defined in the database, is $\simeq 0.1$ in those units. In the
present work, we adopt $\ell_{\rm c}=0.14$ for reasons detailed below.

There is no guide field, hence the Alfvén velocity $v_{\rm A}$ denotes that of
turbulent fluctuations. The units of the simulation are fixed to obtain
$v_{\rm A}=0.41\,c$; accordingly, the rms velocity (as measured without
coarse-graining) is $\langle\delta u^{2}\rangle^{1/2}\simeq 0.4\,c$. This
choice of units allows to remain within the limits of applicability of the
simulation, which is not special-relativistic, and at the same time to offer
the possibility to observe particle acceleration on the finite timescale of
the simulation ($2.8\,\ell_{\rm c}/c$). A fraction $0.1\,$% of the volume
contains velocities in excess of $0.8\,c$, and $0.03\,$% of the volume in
excess of $0.88\,$c, therefore this compromise is reasonable. See also further
below.

###  B.1 Statistics of velocity gradients

The statistics of velocity jumps are obtained through Monte Carlo sampling of
the simulation cube. For a given coarse-graining scale $l$ (expressed in units
of the cube size), $N_{l}\simeq 1/{l}^{3}$ points are drawn at random through
the cube; $N_{l}$ is thus chosen so as to obtain a maximum number of points
without oversampling the volume. $N_{l}$ decreases rapidly with increasing
coarse-graining scale $l$, up to $N_{l}\sim\mathcal{O}(300)$ on the coherence
scale $l\simeq\ell_{\rm c}$. Correspondingly, the sample variance $\sim
1/N_{l}$ becomes non-negligible on those scales.

Using the numerical tools provided by the database, the values of the plasma
velocity ${\bm{v_{\rm p}}}_{l}$, the magnetic field $\bm{B}_{l}$ and their
spatial derivatives $\partial_{i}{{v_{\rm p}}_{j}}_{l}$,
$\partial_{i}{B_{j}}_{l}$ are coarse-grained on scale $l$. Those are used
compute the velocity field ${\bm{v_{E}}}_{l}$ through

${\bm{v_{E}}}_{l}={\bm{v_{\rm p}}}_{l}-\left({\bm{v_{\rm
p}}}_{l}\cdot\bm{b}_{l}\right)\,\bm{b}_{l}\,,$ (B1)

with $\bm{b}_{l}=\bm{B}_{l}/\left|\bm{B}_{l}\right|$. This amounts to
calculating the electric field using ideal Ohm’s law $\bm{E}_{l}=-{\bm{v_{\rm
p}}}_{l}\times\bm{B}_{l}/c$ then
${\bm{v_{E}}}_{l}=c\,\bm{E}_{l}\times\bm{B}_{l}/B_{l}^{2}$. Resistive
corrections to ideal Ohm’s law can be neglected on the scales of interest. The
spatial derivatives of ${\bm{v_{E}}}_{l}$ are similarly derived from those of
${\bm{v_{\rm p}}}_{l}$ and $\bm{B}_{l}$.

Those field values and their derivatives can be used to reconstruct the
quantities entering the force terms, in particular $\Gamma_{\parallel}$ and
$\Gamma_{\perp}$: keeping track of the scale $l$ over which they are defined
through a subscript,

$\displaystyle{\Gamma_{\parallel}}_{l}$
$\displaystyle\,=\,-\bm{b}_{l}\cdot\left(\bm{b}_{l}\cdot\bm{\nabla}\right)\,{\bm{v_{E}}}_{l}\,,$
$\displaystyle{\Gamma_{\perp}}_{l}$
$\displaystyle\,=\,-\frac{1}{2}\left[\bm{\nabla}\cdot{\bm{v_{E}}}_{l}-\bm{b}_{l}\cdot\left(\bm{b}_{l}\cdot\bm{\nabla}\right)\,{\bm{v_{E}}}_{l}\right]$
(B2)

###  B.2 Multifractal model

Those gradients ${\Gamma_{\parallel}}_{l}$ and ${\Gamma_{\perp}}_{l}$ are
random variables that become increasingly non-Gaussian at small scales
$l\ll\ell_{\rm c}$. Such statistics are conveniently discussed within the
multifractal picture of turbulence intermittency [27], which ascribes a local
scaling exponent $h(\bm{x})$ to a position $\bm{x}$, and which describes the
set of locations $\bm{x}$ with index $h(\bm{x})$ as a fractal of dimension
$d(h)$. The spectrum $d(h)$ is a continuous function of $h$.

Intermittency can also be captured through the structure functions, e.g.
$S_{n}(l)=\langle|\Gamma_{l}|^{n}\rangle$, and their scaling exponents
$S_{n}(l)\propto l^{\zeta_{n}}$. The scaling exponent $\zeta_{n}$ and the
spectrum $d(h)$ are however in one-to-one correspondence via the Legendre
transform $\zeta_{n}={\rm min}_{h}\left[h\,n+D-d(h)\right]$ [27] ($D$ number
of spatial dimensions). Therefore, using one or the other is a matter of
choice. We rely here on $d(h)$ to match the statistics of $\Gamma_{l}$ on
scale $l$, in order to characterize the probability distribution function
(p.d.f.) of momentum jumps.

This spectrum $d(h)$ is generally described by a log-normal form [30],

$d(h)=D_{0}-\frac{\left(h-h_{\textsc{mf}}\right)^{2}}{2\sigma_{\textsc{mf}}^{2}}\,,$
(B3)

or a log-Poisson model [35, 36, 30]

$d(h)=D_{0}-\frac{(h-h_{0})}{\ln\beta}\left[1-\ln\left(\frac{h-h_{0}}{(D-D_{0})|\ln\beta|}\right)\right]\,.$
(B4)

In these formulations, $D_{0}$ represent the dimension of the most
intermittent structures, while $h_{\textsc{mf}}$, $\sigma_{\textsc{mf}}$,
$h_{0}$ and $\beta$ are parameters. Currently, those parameters are adjusted
to simulation data, but ultimately, one may hope to gain fundamental insight
on their values.

Theories of turbulence intermittency commonly discuss the statistics of the
Elsässer fields $\bm{\delta{z_{\pm}}_{l}}=\bm{\delta v_{l}}\pm\bm{\delta
b_{l}}$, expressed in terms of the plasma velocity fluctuation $\bm{\delta
v_{l}}={\bm{v_{\rm p}}}_{l}\left(\bm{x}+\bm{l}\right)-{\bm{v_{\rm
p}}}_{l}\left(\bm{x}\right)$ coarse-grained on scale $l$ and the magnetic
field fluctuation (normalized to Alfvén velocity units) $\bm{\delta
b_{l}}=\left[\bm{B}_{l}\left(\bm{x}+\bm{l}\right)-\bm{B}_{l}\left(\bm{x}\right)\right]/\sqrt{4\pi\rho}$
($\rho$ plasma mass density). The quantities $\delta
z_{\pm}\equiv\bm{\delta{z_{\pm}}_{l}}\cdot\bm{l}/|\bm{l}|$ are well described
by log-Poisson models with e.g., $D_{0}=2$, $h_{0}=\frac{1}{9}$ and
$\beta\simeq 0.69$ [31, 32, 33, 34]. Such a model reproduces satisfactorily
the p.d.f. of $\delta z_{\pm}$ in the present MHD simulation, as explicitly
verified.

The measured statistics of ${\Gamma_{\parallel}}_{l}$ and
${\Gamma_{\perp}}_{l}$ differ from those of $\delta z_{\pm}/|\bm{l}|$ (the
additional $1/|\bm{l}|$ is introduced to convert the original Elsässer
quantity into a gradient). The p.d.f. of ${\Gamma_{\parallel}}_{l}$ and
${\Gamma_{\perp}}_{l}$ are systematically more extended and they display
harder powerlaws at large values of the arguments. This difference can be
attributed to the following: ${\Gamma_{\parallel}}_{l}$ and
${\Gamma_{\parallel}}_{l}$ are defined in terms of the velocity field
$\bm{v_{E}}$, whereas $\delta z_{\pm}$ is related to the plasma velocity
$\bm{v_{\rm p}}$; however, gradients of those velocity fields differ by
quantities related to the curvature of magnetic field lines [24], while the
statistics of the field line curvature displays hard powerlaw tails with slope
close to $-2.5$ [37, 38]. In the following, we use a simple log-normal model,
which turns out to provide a satisfactory description of their p.d.f.; more
work is warranted to relate the statistics of ${\Gamma_{\parallel}}_{l}$ and
${\Gamma_{\perp}}_{l}$ with those of $\delta z_{\pm}$.

###  B.3 Multi-fractal fitting procedure

The multi-fractal model for ${\Gamma_{\parallel}}_{l}$ and
${\Gamma_{\perp}}_{l}$ is constructed as described in the text. On the
coherence scale $\ell_{\rm c}$ of the turbulence, Gaussian statistics are well
motivated to model $p_{\Gamma_{\ell_{\rm c}}}$ while on smaller length scales,
the statistics of velocity gradients become increasingly non-Gaussian. We have
set $\ell_{\rm c}=0.14$ in units of the cube size, rather than $0.1$ as it is
defined in the database, because non-Gaussian wings remain visible on the
latter scale and start to disappear on the former. This rescaling exerts no
significant influence on the comparison of the theoretical Green functions to
the data.

The multi-fractal model corresponds to the log-normal spectrum
$d(h)=3-\left(h-h_{\textsc{mf}}\right)^{2}/\left(2\sigma_{\textsc{mf}}^{2}\right)$,
characterized by $\sigma_{\textsc{mf}}=0.9$ and $h_{\textsc{mf}}=-0.2$.

The distributions $p_{\Gamma_{l}}$ (with $\Gamma_{l}$ denoting generically
${\Gamma_{\parallel}}_{l}$ or ${\Gamma_{\perp}}_{l}$) are slightly asymmetric
around zero. This asymmetry appears to depend on scale, it differs between the
parallel and perpendicular direction and it is not always well resolved out of
the noise of the numerical data. Yet, it implies a non-vanishing mean value
for $\Gamma_{l}$ and therefore net advection in momentum space. This mean
value can be measured by averaging over the simulation cube; it is discussed
in Ref. [24]. Net advection can also be measured in the spectra obtained by
particle tracking in the MHD simulation (discussed next). To model this
asymmetry while retaining a minimal number of parameters, we proceed as
follows. We first break the distribution $p_{\Gamma_{l}}$ over positive and
negative values of $\Gamma_{l}$. Both parts are characterized by the same
parameters except the width $\sigma_{\rm c}$ characterizing the fluctuations
on the integral scale, which takes value $\sigma^{+}_{\rm c}\simeq
0.10\,c/\ell_{\rm c}$ for positive values of $\Gamma_{l}$, and
$\sigma^{-}_{\rm c}\simeq 0.07\,c/\ell_{\rm c}$ for negative values of
$\Gamma_{l}$.

The distribution shown in the main text plots the distribution for the
absolute values $|\Gamma_{l}|$, thus summing those two partial distributions.
Overall, the distributions of $\Gamma_{l}$ on all scales are characterized by
four parameters, whose uncertainty is of the order of $0.1$ for
$\sigma_{\textsc{MF}}$ and $h_{\textsc{MF}}$, and $0.05\,c/\ell_{\rm c}$ for
what regards $\sigma^{\pm}_{\rm c}$. This uncertainty exceeds slightly the
difference between $\sigma^{+}_{\rm c}$ and $\sigma^{-}_{\rm c}$: their
difference has been adjusted through its influence over the advection rate, as
measured from the final particle spectra.

A drawback of the above multifractal model is that the integral that defines
$p_{\Gamma_{l}}$ cannot be expressed analytically. To speed up computations,
it proves useful to introduce an approximation of the $\Gamma_{l}$ statistics,
modeled as the broken powerlaw described in the main text. The parameters used
to adjust the gradient statistics are:
$\sigma_{\textsc{bp}}(l)=\sigma_{\textsc{bp}}(\ell_{\rm c})\left(l/\ell_{\rm
c}\right)^{-0.3}$, $k_{0}(l)=2.8+2\left(l/\ell_{\rm c}\right)^{2}$ and
$k_{1}=3$. The dependence of $k_{0}$ on $l$ is to guarantee that, as
$l\rightarrow\ell_{\rm c}$, where the distribution should become approximately
Gaussian, the powerlaw wing does not play a role anymore.

For the same reasons as above, we break this distribution into two sub-
distributions, one for positive values of $\Gamma_{l}$, one for negative
values. They are characterized by the same parameters, except
$\sigma_{\textsc{bp}}(\ell_{\rm c})$, which breaks into
$\sigma^{+}_{\textsc{bp}}(\ell_{\rm c})=0.17\,c/\ell_{\rm c}$ and
$\sigma^{-}_{\textsc{bp}}(\ell_{\rm c})=0.13\,c/\ell_{\rm c}$. The
uncertainties on those parameters are estimated to be $\simeq
0.05\,c/\ell_{\rm c}$ for $\sigma_{\textsc{bp}}(\ell_{\rm c})$, and $\simeq
0.2$ for $k_{0}(l)$.

Finally, it should be stressed that better fits could be obtained by allowing
for more degrees of freedom, but the choice has been made here to retain a
maximum of simplicity in view of the uncertainty characterizing the numerical
data and the number of simplifications made in building up the theoretical
model.

###  B.4 Particle tracking

To reconstruct the Green functions describing the evolution of the momentum
distribution function in time, a large number of particles are tracked in the
simulation (typically $24\,000$ per value of the initial momentum). Each test
particle is initialized at time $t=0$ at a random position, with a random
orientation in velocity space and with a given momentum $p_{0}$ in the
$\mathcal{R}_{\not{E}}$ frame in which the electric field vanishes.

The particle position and momentum is advanced using a Boris pusher. At each
time step of duration $0.1\,r_{\rm g}(p_{0})/c$, $r_{\rm g}(p_{0})$ denoting
the initial gyroradius $p_{0}c/eB$, the pusher queries the database to
retrieve the values of $\bm{v_{\rm p}}$ and $\bm{B}$ at the position of the
particle and at the corresponding time. The time evolution of the
electromagnetic fields is thus properly taken into account. Those values are
obtained using high-order Lagrangian interpolation. The field values are not
coarse-grained here, hence the code advances the particles using the true
time-dependent Lorentz force. The electric field is reconstructed using ideal
Ohm’s law $\bm{E}=-\bm{v_{\rm p}}\times\bm{B}/c$. Finally, the Green functions
are reconstructed by calculating the momentum distribution function of
particles, as a function of time $t$ and as a function of initial momentum
$p_{0}$.

In order to compare the predictions of the theoretical transport equation with
the numerical data, the theoretical spectra, defined in
$\mathcal{R}_{\not{E}}$, are boosted to the simulation frame. Noting that the
boost multiplies the energies in $\mathcal{R}_{\not{E}}$ by a factor
$(1+\mu\,v_{E})/\sqrt{1-v_{E}^{2}/c^{2}}$, with $\mu$ the cosine of the angle
between the particle direction and the velocity field $\bm{v_{E}}$ at that
location ($v_{E}=|\bm{v}_{E}|$), we draw a large number of values $\mu$
uniformly distributed in $[-1,1]$ and we extract a large number of velocity
values $v_{E}$ at random points in the simulation to construct that kernel.
The MHD simulation is not special-relativistic, hence values of $v_{E}$ in
excess of $0.87c$ (corresponding to Lorentz factor $2$) are discarded; this
concerns however a tiny fraction of the numerical grid ($\lesssim 0.03\,\%$).

We consider values of $p_{0}$ such that the wavenumber scale $k\sim r_{\rm
g}^{-1}$ associated to the initial gyroradius $r_{\rm g}(p_{0})$ falls within
the inertial range of the turbulence cascade, above the dissipative range.
This corresponds, in practice to values $10^{-2}\,\ell_{\rm c}\lesssim r_{\rm
g}(p_{0})\lesssim 0.15\,\ell_{\rm c}$. In terms of coarse-graining scales,
which are used in the theoretical computation, this corresponds to
$0.06\lesssim l_{\rm g}/\ell_{\rm c}\lesssim 1$.

###  B.5 Sensitivity to the fitting procedure

Figure B1: Particle spectra $p\,n_{p}(t)$ obtained by tracking a large number
of particles in the MHD simulation (symbols), with theoretical spectra
overlaid. Solid blue line: spectrum for the broken powerlaw model with
parameters adopted in the main text; dashed (resp. dash-dotted) red lines:
spectra obtained by decreasing (resp. increasing) $\sigma_{\textsc{bp}}$ by
$0.05$; long-dashed (resp. short-dash-dotted) blue lines: spectra obtained for
$k_{0}(l)=2.6+2\left(l/\ell_{\rm c}\right)^{2}$ [resp.
$k_{0}(l)=3.0+2\left(l/\ell_{\rm c}\right)^{2}$]; densely dashed magenta line:
changing the cut-off function $\exp\left[-\left(\ln p-\ln
p^{\prime}\right)^{2}/a\right]$ from $a=1$ (main text) to $a=0.3$; densely
dotted magenta line: same, but with $a=3$. See text for details.

The spectra $n_{p}(t)$ that are derived by integrating the kinetic equation
depend directly on the model of velocity gradient statistics. The two models
plotted in red (multi-fractal) and blue (broken powerlaw) offer a first glance
at this influence. Additionally, we plot in Fig. B1 a number of spectra
obtained by varying the parameters of the broken powerlaw approximation within
the quoted uncertainties. The solid blue corresponds to the spectrum discussed
in the main text, with the choice of parameters indicated above. The dashed
and dash-dotted red lines show the spectra obtained by respectively,
decreasing or increasing the values of $\sigma_{\textsc{bp}}(\ell_{\rm c})$ by
$0.05\,c/\ell_{\rm c}$. The long-dashed and short-dash-dotted blue lines show
the spectra obtained by changing the $2.8$ exponent of $k_{0}$ to,
respectively, $2.6$ and $3.0$. Finally, the densely dashed and densely dotted
magenta lines show the effect of changing the cut-off function
$\exp\left[-\left(\ln p-\ln p^{\prime}\right)^{2}/a\right]$ from $a=1$ (main
text) to respectively, $a=0.3$ and $a=3$. The model for $a=0.3$ also assumes:
$k_{0}(l)=2.6+2(l/\ell_{\rm c})^{2}$,
$\sigma_{\textsc{bp}}^{-}=0.15\,c/\ell_{\rm c}$,
$\sigma_{\textsc{bp}}^{+}=0.20\,c/\ell_{\rm c}$, and
$\sigma_{\textsc{bp}}\propto l^{-0.2}$; the model for $a=3.$ assumes:
$k_{0}(l)=3.0+2(\l/\ell_{\rm c})^{2}$, other parameters unchanged.


In this paper, we propose an algorithm, Epochal Difficult Captions, to supplement the training of any model for the Automated Audio Captioning task. Epochal Difficult Captions is an elegant evolution to the keyword estimation task that previous work have used to train the encoder of the AAC model. Epochal Difficult Captions modifies the target captions based on a curriculum and a difficulty level determined as a function of current epoch. Epochal Difficult Captions can be used with any model architecture and is a lightweight function that does not increase training time. We test our results on three systems and show that using Epochal Difficult Captions consistently improves performance.

§ INTRODUCTION

The Automated Audio Captioning task is a cross-modal task where the model is trained to generate a descriptive sentence for an input audio. This task has many practical real world use cases, such as closed captioning for people who are deaf or hard of hearing. Automated Audio Captioning describes the sound events that are happening in both the foreground and background of the audio, as opposed to Automated Speech Recognition, which predicts the low time frame phonemes and utterances.

The Automated Audio Captioning task requires the model to be able to capture both high time frame and low time frame events that are occurring throughout the audio. As such, previous work have often used pretrained models such as PANNs that have already been trained on audio tagging in order to aid the model in generating factually correct caption. In a similar vein, several authors [1, 2] have used keyword estimation, a supplementary and intermediate objective where the model has to predict one or a few keywords from the encoder to help the decoder to generate better captions. The keyword estimation task has shown to have potential in helping the model perform better. However, the keyword estimation task requires heavy preprocessing to create the audio-keywords pairs for training. Furthermore, while the different authors have the same motivation behind keyword estimation task, their implementations and audio-keywords pairs are often different and subject to their own interpretation. This makes it hard to compare results and methods across different papers. We find that there is an elegant solution to this. Epochal Difficult Captions modify the training targets directly during the training loop. Based on a curriculum, the training captions start off with a higher keywords to sentence ratio during the earlier epochs in the the training run. This is done by removing stopwords directly from the batched training targets. As the training progresses, the keywords to sentence ratio will be gradually reduced until the sentence reverts to its original unmodified sentence.

We test Epochal Difficult Captions using two different publicly available systems [3, 4] from previous work. Both systems are of the typical encoder-decoder structure. The encoder consists of either a single [4] pretrained CNN or both [3] a pretrained CNN and transformer encoder. The decoder used is a 2 layer transformer decoder. The extracted audio features from the encoder is passed to the decoder for cross attention. Both systems are trained using cross entropy, but [4] has the additional option of using reinforcement learning to optimize the metrics directly. We experiment using both cross entropy, and reinforcement learning where available. We find that using Epochal Difficult Captions consistently improves performance across these settings with no noticeable extra compute time and power needed. Epochal Difficult Captions is used only during training and hence does not affect any part of the inference process.

The rest of this paper follows this structure. In section <ref>, we introduce related work pertaining to Automated Audio Captioning datasets, popular model architectures and training objectives. We also discuss the keyword estimation task, a popular intermediate task in Automated Audio Captioning. In Section <ref>, we explain our approach, Epochal Difficult Captions for Curriculum Learning. This section shows how we derive our difficulty curve and levels, and how stopwords are removed as an function of the difficulty level. In Section <ref>, we specify our experimental settings, details, and process to evaluate the effectiveness of Epochal Difficult Captions. In Section <ref>, we examine the performance of the models trained with Epochal Difficult Captions for Curriculum Learning. Finally, we give our concluding remarks in Section <ref>.

§ RELATED WORK

§.§ Audio Captioning Datasets

The Clotho Dataset v2.1 [5] consists of 6974 audio samples of duration from 15 to 30 seconds. Each audio clip has 5 corresponding captions annotated by humans. Each caption has a length of 8 to 20 words, and the captions describe the events in that audio. The Clotho Dataset v1 was augmented by including more training samples and a validation split to form Clotho v2.1. The Clotho Dataset was created to address the critiques of the larger Audiocaps dataset [6], which was found to contain biases in the annotation process. There are also other niche datasets like the Hospital Scene and Car Scene dataset [7] which are Mandarin annotated.

In this work, we focus only on the Clotho Dataset v2.1, henceforth referred to as the Clotho Dataset.

§.§ Model Architectures and Objectives

Most work use an encoder-decoder model architecture. The encoder is used to encode audio information to produce audio features, which is then passed to the decoder for cross attention and to generate the caption. Most architectures converge to a similar theme - using a pretrained convolutional neural network such as PANNs for encoding, and a transformer decoder for decoding. However, the training process often differs from the pretraining stage and an intermediate or supplementary objective.

For instance, the AT-CNN [8], and the Audio Captioning Transformer [9] pretrains their model encoder on audio tagging as an pretraining task before finetuning their model directly on Automated Audio Captioning. Other authors also use self supervised supplementary objectives [3] or caption retrieval [10] to guide the model training. Self-Critical Sequence
Training [11], a reinforcement learning tactic used in Image Captioning, has also been used to optimize the model for Audio Captioning. Keyword Estimation, which we will go into detail in Section <ref>, has also been a popular choice [12, 2] of an intermediate objective.

In this work, we build on the work of [3] and [4] by using their model architectures that has been proven to achieve competitive results. Henceforth, we refer [4] and [3] as System 1 and System 2 respectively. Figure <ref> delineates their system architectures.

Detailed overview of System 1 and System 2. Both System 1 and 2 uses a pretrained CNN10 in the encoder. System 2 has an additional 2 layer transformer appended to the CNN10 in the encoder. For the decoder, both System 1 and 2 uses a transformer decoder. However, System 1 uses pretrained word2vec word embeddings, while System 2 trains the word embeddings from scratch.

§.§ Keyword Estimation

Several authors have applied the Keyword Estimation task as a intermediate task to aid the training of the model for Automated Audio Captioning. This is done by decomposing the Automated Audio Captioning task into two smaller tasks, Keyword Estimation and Caption Generation. The model encoder then performs Keyword Estimation to predict a set of keywords, which is then passed to the decoder for generation. The Keyword Estimation task helps to mitigate the problem of indeterminacy in word selection [12] hence reducing the solution space for the decoder to generate better captions.

Since the Clotho Dataset does not come with labels for keywords, authors often use their own interpretation to create labels for the keyword Estimation task. [12] uses a rule based keyword extraction method to produce keywords labels for training. Frequent word lemmas of nouns, verbs, adjectives, and adverbs as regarded as keywords. On the other hand, [2] extracts their keywords by using the Natural Language Toolkit[https://github.com/nltk/nltk] to discard `useless words' such as make, go, etc. Remaining verbs are used in their original forms, and the nouns are not changed.

Firstly, given keyword estimation requires some form of pre-processing subject to the interpretation of the researcher, it is difficult to compare the effectiveness of the keyword estimation task across different pieces of work. Secondly, decomposing the Automated Audio Captioning task into subtasks keyword estimation and caption generation requires the decoder to rely on the prediction of the keywords from the encoder. This might result in error propagation if the keywords are inaccurate. With these considerations, we propose the Epochal Difficult Captions for curriculum learning, which has the same motivation behind keyword estimation, but does not require the decomposition of Automated Audio Captioning into subtasks.

§ EPOCHAL DIFFICULT CAPTIONS FOR CURRICULUM LEARNING

Epochal Difficult Captions for Curriculum Learning is an algorithm to modify the target captions relative to the current stage of the training setup. The difficulty level, $D$, is a function of the current epoch. $D$ indicates how the curriculum modifies the target captions. In this paper, we use the stopwords curriculum. While we have another frequency-based curriculum, we omitted it due to space constraints.

§.§ Difficulty level as a function of epoch

The difficulty level in curriculum learning dictates the complexity of the captions and it should increase from a lower to a higher magnitude relative to the current epoch. We use an exponential function (Equation <ref>) followed by a floor and ceiling function (Equation <ref>) to generate a difficulty value $D$ for all non-negative epoch values \(D  \in  (0, 1)\).

\begin{equation} \label{eq:diff}
D^{\prime} = 1 - e^{( -\alpha * epoch )}
\end{equation}

\begin{equation} \label{eq:diff2}
0.05,           & \text{if } D < 0.05 \\
D^{\prime},     & \text{otherwise}
\end{cases}
\end{equation}

where the hyperparameter \(\alpha\) in Equation <ref> controls the rate of increase of difficulty. At $D = 1$, the original captions are not modified. The difficulty level should asymptotically tend to 1 nearing the max epoch. Therefore, \(\alpha\) requires the max epoch hyperparameter to be taken into consideration. In our experiments, we use \(\alpha = 0.20\) when max epoch is 30, \(\alpha = 0.10\) when max epoch is 60, and \(\alpha = 0.05\) when max epoch is 100.

axis y line*=left,
ylabel=total tokens per epoch,
legend pos=south east
only marks,
mark= *,
mark size=1pt]
table[x=difficulty, y=total_tokens_per_epoch]

axis y line*=right,
axis x line=none,
ylabel=average unique tokens per batch,
legend pos=south east
only marks,
mark= x,
mark size=2pt]
table[x=difficulty, y=avg_tokens_per_batch]
refstyle=plot_epoch_d_100total tokens per epoch
refstyle=plot_epoch_avgtokensperbatch_100average unique tokens per batch

100 epochs

axis y line*=left,
ylabel=total tokens per epoch,
legend pos=south east
only marks,
mark= *,
mark size=1pt]
table[x=difficulty, y=total_tokens_per_epoch]

axis y line*=right,
axis x line=none,
ylabel=average unique tokens per batch,
legend pos=south east
only marks,
mark= x,
mark size=2pt]
table[x=difficulty, y=avg_tokens_per_batch]
refstyle=plot_epoch_d_60total tokens per epoch
refstyle=plot_epoch_avgtokensperbatch_60average unique tokens per batch

60 epochs

axis y line*=left,
ylabel=total tokens per epoch,
legend pos=south east
only marks,
mark= *,
mark size=1pt]
table[x=difficulty, y=total_tokens_per_epoch]

axis y line*=right,
axis x line=none,
ylabel=average unique tokens per batch,
legend pos=south east
only marks,
mark= x,
mark size=2pt]
table[x=difficulty, y=avg_tokens_per_batch]
refstyle=plot_epoch_d_30total tokens per epoch
refstyle=plot_epoch_avgtokensperbatch_30average unique tokens per batch

30 epochs

Scatter plot of $D$ (x-axis) against total tokens per epoch (left y-axis) and average unique tokens per batch (right y-axis) for different epoch values 100 (a), 60(b), 30(c)

§.§ Stopwords Curriculum

In the stopwords curriculum, the target captions for each audio are modified by removing stopwords. Stopwords are common words that carry little useful information and can be obtained from the public nltk library. Some examples of stopwords are ‘for’, ‘do’, ‘its’, ‘yours’, ‘such’, ‘into’. Removal of these stopwords mimic the keyword estimation task described in Section <ref> by only including unique and less prevalent words in the caption.

Stopwords Curriculum

[1]
$y$ - the original target caption
$D$ - difficulty value (Equation <ref>, <ref>)
$SW$ - list of stopwords.
$y^\prime$ - modified caption.
$P(remove) \gets (1 - D)$
$y^\prime \gets []$
$word \in y$
$word \in stop\_words$
$s \gets$ sample uniform(0,1)
$s <  P(remove)$
$y^\prime \gets new\_sentence$
$y^\prime \gets y^\prime + word$
$y^\prime \gets y^\prime + word$

Algorithm <ref> describes the stopwords curriculum. Each stopword in the caption is an candidate for removal depending on the current difficulty level of the current epoch. The greater the difficulty level $D$ is, the lower the probability of removing a stopword from the caption.

In Figure <ref>, we plot the relationship between $D$, the total number of words used per epoch in training, and the average number of unique words per batch during training. At low difficulty levels, the total number of tokens per epoch increases rapidly with each epoch at low difficulty levels. At higher difficulty levels above 0.9, the total number of tokens fluctuates around the same range. Using such a token distribution based on a fixed curriculum allows us to gradually expose the model to progressively more unique words.

Model                                        BLEU1     BLEU2     BLEU3     BLEU4     ROUGEL     METEOR     CIDEr     SPICE     SPIDEr

System 1 - cross entropy   	0.559   	0.358   	0.237   	0.153   	0.169   	0.374   	0.382   	0.116   	0.249

System 1 - cross entropy + stopwords   	0.558   	0.362   	0.242   	0.159   	0.170   	0.375   	0.391   	0.115   	0.253

System 1 - scst   	0.641   	0.417   	0.277   	0.174   	0.182   	0.407   	0.432   	0.124   	0.278

System 1 - scst + stopwords	   0.642   	0.409   	0.272   	0.172   	0.182   	0.402   	0.444   	0.124   	0.284

System 2 - cross entropy     0.553     0.367     0.248     0.160     0.162     0.372     0.359     0.111      0.235

System 2 - cross entropy + stopwords        0.558   	0.376   	0.258	    0.172   	0.167   	0.376   	0.381   	0.115   	0.248

Comparison of performance of systems trained on Epochal Difficult Captions with the Stopwords curriculum against their counterpart

§ EXPERIMENTAL DETAILS

§.§ Data

We use Clotho dataset v2.1 for all our experiments as mentioned in Section <ref>. Since we are building on the previous systems of two authors [3, 4], we use the same settings to preprocess the waveforms into the spectrograms as inputs into the models. For both systems, we use 64 Mel-bands, sampling rate of 44100, FFT window length of 1024, and a hop size of 512.

§.§ Training and Evaluation

For our experiments on both System 1 and System 2, we try to replicate the original work without using Epochal Difficult Captions to act as a control experiment, then rerun the same experiments with Epochal Difficult Captions. Where possible, we use the original hyperparameters.

For System 1 [4], we use a batch size of 32 with no gradient accumulation for 30 epochs and 60 epoch with no early stopping. System 1 has two different training setups. The first setup first optimizes on a cross entropy loss for 30 epochs, then uses a reinforcement learning technique, Self Critical Sequence Training (SCST) [11], to train the model for 60 epochs. Learning rate is set to $1 \times 10^{-4}$ for the cross entropy loss and $5 \times 10^{-5}$ for the reinforcement training. SpecAugmentation [13] is applied to all log mel-spectrogram inputs as an data augmentation tactic. We do apply label smoothing. For inference, we perform beam search with beam size of 3 for decoding.

For System 2 [3], we use a batch size of 64 with a gradient accumulation steps of 4 for 100 epochs. Learning rate is set to $3 \times 10^{-4}$ and SpecAugmentation [13] is applied to all log mel-spectrogram inputs as an data augmentation tactic. We do not apply label smoothing. For inference, we perform beam search with beam size of 4 for decoding.

The COCO image captioning evaluation process [14] is used to evaluate the generated caption. There are a few metrics used. The BLEUn scores [15] measures of n-gram overlap between the generated and reference caption, while ROUGEL [16] bases its score off the longest common sequence. METEOR [17] uses a harmonic mean of unigram precision and recall to score the generated sequence. CIDEr [18] uses the average cosine similarity between the candidate sentence and the reference sentences to produce a score. Next, SPICE [19] uses semantic propositional information to score the generated captions. Finally, the SPIDEr metric is a combination of both SPICE and CIDEr. We use all the aforementioned metrics for evaluation. The SPIDEr metric is however the metric that researchers try to beat.

§ EXPERIMENTAL RESULTS AND ANALYSIS

We compare the effectiveness of using Epochal Difficult Captions for Curriculum Learning in 3 different settings. The results are shown in Table <ref>. From the table, we observe that the most meaningful metric, the SPIDEr score, increases across all 3 systems. For System 2 particularly, the effect of using Epochal Difficult Captions leads to an improvement of all metrics.

As described in Section <ref>, Epochal Difficult Captions allows the model to gradually transit from training on prominent keywords to the original unmodified caption. We believe this allows the model to focus on more important words such as nouns and action verbs in the earlier stage of the training, thus providing the model with a better ability to recognize objects and actions from the log mel spectrogram. As the training proceeds, filler words and syntactically important terms are consecutively added epoch-wise to enable the model to generate grammatically sound sentences.

§ CONCLUSION

This work introduces the use of Epochal Difficult Captions for Curriculum Learning. Epochal Difficult Captions is an elegant evolution to the keyword estimation problem. We have examined the difficulty curve for Curriculum Learning and specified the stopwords Curriculum. We have shown that Epochal Difficult Captions can be easily added into any Automated Audio Captioning system during the training stage. Using Epochal Difficult Captions for Curriculum Learning, performance across 3 different settings in 2 systems improves across multiple metrics and the SPIDEr score.

[1]
Yuma Koizumi, Ryo Masumura, Kyosuke Nishida, Masahiro Yasuda, and Shoichiro
“A Transformer-based Audio Captioning Model with Keyword
arXiv:2007.00222 [cs, eess, stat], Aug. 2020,
arXiv: 2007.00222.

[2]
Zhongjie Ye, Helin Wang, Dongchao Yang, and Yuexian Zou,
“Improving the performance of automated audio captioning via
integrating the acoustic and semantic information,” 2021.

[3]
Andrew Koh, Fuzhao Xue, and Eng Siong Chng,
“Automated audio captioning using transfer learning and
reconstruction latent space similarity regularization,” 2021.

[4]
Xinhao Mei, Qiushi Huang, Xubo Liu, Gengyun Chen, Jingqian Wu, Yusong Wu,
Jinzheng ZHAO, Shengchen Li, Tom Ko, H. Tang, Xi Shao, Mark D. Plumbley, and
Wenwu Wang,
“An encoder-decoder based audio captioning system with transfer and
reinforcement learning,”
in Proceedings of the 6th Detection and Classification of
Acoustic Scenes and Events 2021 Workshop (DCASE2021), Barcelona, Spain,
November 2021, pp. 206–210.

[5]
Konstantinos Drossos, Samuel Lipping, and Tuomas Virtanen,
“Clotho: An Audio Captioning Dataset,”
arXiv:1910.09387 [cs, eess], Oct. 2019,
arXiv: 1910.09387.

[6]
Chris Dongjoo Kim, Byeongchang Kim, Hyunmin Lee, and Gunhee Kim,
“AudioCaps: Generating Captions for Audios in The
in Proceedings of the 2019 Conference of the North
American Chapter of the Association for Computational Linguistics:
Human Language Technologies, Volume 1 (Long and Short Papers),
Minneapolis, Minnesota, June 2019, pp. 119–132, Association for
Computational Linguistics.

[7]
Xuenan Xu, Heinrich Dinkel, Mengyue Wu, and Kai Yu,
“Audio Caption in a Car Setting with a Sentence-Level
arXiv:1905.13448 [cs, eess], Oct. 2020,
arXiv: 1905.13448.

[8]
Xuenan Xu, Heinrich Dinkel, Mengyue Wu, Zeyu Xie, and Kai Yu,
“Investigating local and global information for automated audio
captioning with transfer learning,”
in ICASSP 2021 - 2021 IEEE International Conference on
Acoustics, Speech and Signal Processing (ICASSP), 2021, pp. 905–909.

[9]
Xinhao Mei, Xubo Liu, Qiushi Huang, Mark D. Plumbley, and Wenwu Wang,
“Audio captioning transformer,” 2021.

[10]
Tianyi Zhang, Varsha Kishore, Felix Wu, Kilian Q. Weinberger, and Yoav Artzi,
“Bertscore: Evaluating text generation with bert,” 2020.

[11]
Steven J. Rennie, Etienne Marcheret, Youssef Mroueh, Jarret Ross, and Vaibhava
“Self-critical sequence training for image captioning,” 2016.

[12]
Yuma Koizumi, Yasunori Ohishi, Daisuke Niizumi, Daiki Takeuchi, and Masahiro
“Audio Captioning using Pre-Trained Large-Scale Language
Model Guided by Audio-based Similar Caption Retrieval,”
arXiv:2012.07331 [cs, eess], Dec. 2020,
arXiv: 2012.07331.

[13]
Daniel S. Park, William Chan, Yu Zhang, Chung-Cheng Chiu, Barret Zoph, Ekin D.
Cubuk, and Quoc V. Le,
“Specaugment: A simple data augmentation method for automatic speech
Interspeech 2019, Sep 2019.

[14]
Tsung-Yi Lin, Michael Maire, Serge Belongie, Lubomir Bourdev, Ross Girshick,
James Hays, Pietro Perona, Deva Ramanan, C. Lawrence Zitnick, and Piotr
“Microsoft coco: Common objects in context,” 2015.

[15]
Kishore Papineni, Salim Roukos, Todd Ward, and Wei-Jing Zhu,
“Bleu: A method for automatic evaluation of machine translation,”
in Proceedings of the 40th Annual Meeting on Association for
Computational Linguistics, USA, 2002, ACL '02, p. 311–318, Association for
Computational Linguistics.

[16]
Chin-Yew Lin,
“ROUGE: A package for automatic evaluation of summaries,”
in Text Summarization Branches Out, Barcelona, Spain, July
2004, pp. 74–81, Association for Computational Linguistics.

[17]
Michael Denkowski and Alon Lavie,
“Meteor universal: Language specific translation evaluation for any
target language,”
in Proceedings of the EACL 2014 Workshop on Statistical Machine
Translation, 2014.

[18]
Ramakrishna Vedantam, C. Lawrence Zitnick, and Devi Parikh,
“Cider: Consensus-based image description evaluation,”
in 2015 IEEE Conference on Computer Vision and Pattern
Recognition (CVPR), 2015, pp. 4566–4575.

[19]
Peter Anderson, Basura Fernando, Mark Johnson, and Stephen Gould,
“Spice: Semantic propositional image caption evaluation,” 2016.