compat=1.1.0 /tikzfeynman/every boson@@/.style= /tikz/draw=none,
/tikz/decoration=name=none, /tikz/postaction= /tikz/draw, /tikz/decoration=
zigzag,segment length=1mm, amplitude=0.6mm, , /tikz/decorate=true,

# Velocity decorrelation functions of high-energy cosmic rays propagating in
magnetic fields

Olivier Deligny  Laboratoire de Physique des 2 Infinis Irène Joliot-Curie
(IJCLab)
CNRS/IN2P3, Université Paris-Saclay, Orsay, France

(August 31, 2024)

###### Abstract

Diffusion tensor coefficients play a central role in describing cosmic-ray
transport in various astrophysical environments permeated with magnetic
fields, which are usually modeled as a fluctuating field on top of a mean
field. In this article, a formal derivation of these coefficients is presented
by means of the calculation of velocity decorrelation functions of particles.
It relies mainly on expanding the 2-pt correlation function of the
(fluctuating) magnetic field experienced by the particles between two
successive times in the form of an infinite Dyson series and retaining a class
of terms that converge to a physical solution. Subsequently, the velocity
decorrelation functions, themselves expressed as Dyson series, are deduced
from an iteration procedure that improves on the partial summation scheme. The
results are shown to provide approximate solutions compared to those obtained
by Monte-Carlo simulations as long as the Larmor radius of the particles is
larger than at least one tenth of the largest scale of the turbulence.

##  I Introduction

Cosmic rays compose less than one particle out of ten million in the
interstellar gas. Still, their average energy density is similar to that of
the gas. A small proportion of particles has therefore appropriated a
substantial part of the available energy. The population of GeV cosmic rays is
the most numerous and has important influence on the gas dynamics, heating and
cooling the inter-stellar medium, launching large-scale outflows and modifying
the phase structure of the background gas [1, 2, 3, 4]. Modeling their
transport calls for highly non-linear processes, as they excite magnetic
structures that control themselves the transport. When their streaming through
the background plasma is faster than the Alfvén speed, cosmic rays transfer
part of their energy and momentum to the gas by exciting instabilities [5, 6,
7, 8]. The streaming is possibly enhanced by the reduction of the Alfvén waves
damped by interstellar dust grains [9, 10] or by magnetrohydrodynamic
turbulence cascade [11]. Recent studies have probed the key role of those
cosmic rays on driving strong galactic winds or on matching the observed
ionization structure of the circum-galactic medium (the diffuse gas made of
the majority of galactic baryons outside of the galactic disk) [2, 3, 12], or
on star and galaxy formation [4].

The energy density of trans-GeV cosmic rays does not dominate any longer that
of the gas. Alfvén waves are then considered as an “extrinsic turbulence”
excited by turbulent motions in the thermal gas [13, 14]. The study presented
in this article focuses on the transport of cosmic rays in this type of
extrinsic turbulence. These environments constitute collisionless turbulent
plasmas in which particles propagate or are accelerated. The confinement and
transport of high-energy cosmic rays are governed by their scattering off the
fluctuating magnetic fields, which act as an effective source of collisions
[15]. Considering these collisions as a relaxation process, they tend to bring
the average velocity distribution of the particles to its isotropic mean [16].
Under these conditions, the flux of particles can be related to their gradient
of density by means of a diffusion tensor $D_{ij}$, which can be expressed in
terms of the magnetic field unit vector $\mathbf{b}$, the diffusion
coefficients parallel and perpendicular to the mean field $D_{\parallel}$ and
$D_{\perp}$, and the anti-symmetric diffusion coefficient $D_{A}$ describing
the particle drifts as [17]

$D_{ij}=D_{\perp}\delta_{ij}+\left(D_{\parallel}-D_{\perp}\right)b_{i}b_{j}+D_{A}\epsilon_{ijk}b_{k}.$
(1)

The determination of the coefficients proves to be a difficult task, as it
requires dealing with a highly non-linear problem in several respects. Robust
estimates rely therefore on numerical simulations exploring wide ranges of
particle rigidities and turbulence levels [18, 19, 20, 21, 22, 23, 24, 25, 26,
27, 28, 29, 30, 31, 32, 33, 34, 35]. They show that the various approximations
proposed in the literature fail to reproduce quantitatively and qualitatively
the numerical results in various ranges of turbulence levels [19, 36].
Estimations based on the classical scattering theory [37], on the quasi-linear
theory [38, 39], or on the intuitive ansatz using exponential decreases in the
velocity decorrelation functions proposed in [40] do not capture the memory
effects uncovered in simulations and yield to diffusion coefficients largely
underestimated or overestimated. In this paper, a formal derivation based on
first principles is explored. The range of application pertains to particles
with a Larmor radius larger than at least one tenth of the largest scale of
the turbulence. Consequently, the solution can be considered as an approximate
one to describe the transport of cosmic rays in the high-rigidity regime
(quasi-ballistic propagation) and in the rigidity range marking the transition
between gyroresonant scattering (i.e. Larmor radius of particles entering in
resonance with wavelengths of the turbulence) and quasi-ballistic propagation.

The derivation builds on the theoretical one presented in [25]. In that work,
an estimation of the average velocity of the particles propagating in the
turbulence as a function of time, expressed as a Dyson series, is achieved by
using a white-noise model for the 2-pt function of the magnetic field
experienced between two successive times. Such a modeling proves to be
accurate in the high rigidity regime, in which the values of the magnetic
field experienced by the particles decorrelate on time scales much smaller
than that of the scattering. Under these conditions, the summation of the
Dyson series gives rise to an exponential decay of the average velocity
characteristic of a Markovian process. The diffusion coefficients are then
inferred from velocity correlation function, $\langle v_{0i}v_{j}(t)\rangle$,
through a time integration [41],

$D_{ij}(t)=\int_{0}^{t}\mathrm{d}t^{\prime}\leavevmode\nobreak\ \langle
v_{0i}v_{j}(t^{\prime})\rangle,$ (2)

in the limit that $t\rightarrow\infty$. Here, $v_{0i}\equiv v_{i}(t=0)$ and
$\langle\cdot\rangle$ stands for the average quantities, taken over several
space and time correlation scales of the turbulent field. Since cosmic rays
are high-energy relativistic particles, the norm of the velocity is identified
to $c$ for convenience.

These results were extended to a range of rigidities gyroresonant with the
power spectrum of turbulence by allowing a short-time memory for the time
evolution of the particle velocities through a red-noise approximation for the
2-pt function of the magnetic field experienced between two successive times
[42]. The extension was limited, however, to the case of pure turbulence. In
this paper, the same kind of formalism is applied to any turbulence level. In
Section II, the relevant Dyson series and diagrammatic representations are
introduced to infer the time evolution of the particle velocities. The
modeling, beyond white- or red-noise approximations, of the 2-pt function of
the magnetic field experienced between two successive times is addressed in
Section III. The general diagrammatic technique to carry out the partial-
summation of the Dyson series for the velocity decorrelation functions is
presented in Section IV. The formalism is then applied to the case of parallel
diffusion coefficient in Section V, perpendicular diffusion coefficient in
Section VI, and anti-symmetric diffusion coefficient in Section VII. Based on
Monte-Carlo results (see Appendix A for details about the Monte-Carlo
generator), the limitations of the results obtained are underlined. General
conclusions are drawn in Section VIII. In addition, analytical approximations
obtained in the red-noise approximation for the 2-pt function of the magnetic
field experienced between two successive times are given in Appendix C.

##  II Transport of cosmic rays in magnetic fields

We are interested in determining the moments of $v_{i}(t)$ to derive a
workable expression for Eqn. 2. Adopting the convention of implicit summation
on repeated indices throughout the paper, the velocity of each test-particle
is governed by the Lorentz-Newton equation of motion,

$\dot{v}_{i}(t)=\delta\Omega\leavevmode\nobreak\ \epsilon_{ijk}v_{j}(t)\delta
b_{k}(t)+\Omega_{0}\leavevmode\nobreak\ \epsilon_{ijk}v_{j}(t)b_{0k}(t).$ (3)

Here, $\delta\Omega=c^{2}Z|e|\delta B/E$ is the gyrofrequency related to the
turbulence with root mean square $\delta B$ for each component, $\Omega_{0}$
is the one related to the mean field oriented, to fix the ideas, such that
$\mathbf{B}_{0}=B_{0}\mathbf{u}_{z}$, $Z|e|$ and $E$ the electric charge and
the energy of the particle respectively, and $\delta b_{k}(t)\equiv\delta
b_{k}(\mathbf{x}(t))$ the $k$-th component of the turbulence (expressed in
units of $\delta B$) at the spatial coordinate $\mathbf{x}(t)$, which
corresponds to the position of the test-particle at time $t$. A formal
solution for $\langle v_{i}(t)\rangle$ can be obtained by expressing the
solution of Eqn. 3 as an infinite number of Dyson series, each combining terms
in powers of $\delta\mathbf{b}$ coupled to terms in powers of
$\mathbf{B}_{0}$. Dealing with such an infinite number of Dyson series is
however hardly manageable. To circumvent this difficulty, we use the auxiliary
variable introduced in [25], $w_{i}(t)=R_{ij}(\Omega_{0}t)v_{j}(t)$, with
$\hat{R}(\Omega_{0}t)$ the rotation matrix of angle $\Omega_{0}t$ around
$\mathbf{u}_{z}$. The equation of motion for $\mathbf{w}$ is then

$R^{-1}_{ij}(\Omega_{0}t)\dot{w}_{j}(t)+\dot{R}^{-1}_{ij}(\Omega_{0}t)w_{j}(t)=\\\
\delta\Omega\leavevmode\nobreak\
\epsilon_{ijk}R^{-1}_{j\ell}(\Omega_{0}t)w_{\ell}(t)\delta
b_{k}(t)+\Omega_{0}\leavevmode\nobreak\
\epsilon_{ijk}R^{-1}_{j\ell}(\Omega_{0}t)w_{\ell}(t)b_{0k},$ (4)

which, taking advantage of $b_{0k}=\delta_{kz}$ throughout the paper, reads as

$\dot{w}_{i}(t)=\delta\Omega\leavevmode\nobreak\
R_{ij}(\Omega_{0}t)\epsilon_{jk\ell}R^{-1}_{km}(\Omega_{0}t)w_{m}(t)\delta
b_{\ell}(t)\\\ -R_{ij}(\Omega_{0}t)\dot{R}^{-1}_{jm}(\Omega_{0}t)w_{m}(t)\\\
+\Omega_{0}\leavevmode\nobreak\
R_{ij}(\Omega_{0}t)\epsilon_{jkz}R^{-1}_{km}(\Omega_{0}t)w_{m}(t).$ (5)

Noting that $\dot{R}^{-1}_{ij}(\Omega_{0}t)=\Omega_{0}\epsilon_{i\ell
z}R^{-1}_{\ell j}(\Omega_{0}t)$, the two terms of the second line of Eqn. 5
cancel so that

$\dot{w}_{i}(t)=\delta\Omega\leavevmode\nobreak\
R_{ij}(\Omega_{0}t)\epsilon_{jk\ell}R^{-1}_{km}(\Omega_{0}t)w_{m}(t)\delta
b_{\ell}(t).$ (6)

Hence, the equation of motion for $\mathbf{w}$ is similar to that of
$\mathbf{v}$ in a pure turbulence except for the action of the rotation
matrices: for each infinitesimal time step, $\mathbf{w}$ is rotated by
$-\Omega_{0}t$ around $\mathbf{u}_{z}$ prior to undergoing the impact of
$\delta\mathbf{b}$, and is rotated by $+\Omega_{0}t$ afterwards.

The formal solution for the average $\langle w_{i_{0}}(t)\rangle$ can be
expressed as a single Dyson series:

$\langle
w_{i_{0}}(t)\rangle=w_{0i_{0}}+\sum_{p=1}^{\infty}\delta\Omega^{p}\leavevmode\nobreak\
\epsilon_{k_{1}m_{1}n_{1}}\epsilon_{k_{2}m_{2}n_{2}}\dots\epsilon_{k_{p}m_{p}n_{p}}w_{0i_{p}}\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\dots\int_{0}^{t_{p-1}}\mathrm{d}t_{p}\\\
R_{i_{0}k_{1}}(\Omega_{0}t_{1})R_{i_{1}k_{2}}(\Omega_{0}t_{2})\cdots
R_{i_{p-1}k_{p}}(\Omega_{0}t_{p})R^{-1}_{m_{1}i_{1}}(\Omega_{0}t_{1})R^{-1}_{m_{2}i_{2}}(\Omega_{0}t_{2})\cdots
R^{-1}_{m_{p}i_{p}}(t_{p})\langle\delta b_{n_{1}}(t_{1})\dots\delta
b_{n_{p}}(t_{p})\rangle,$ (7)

using the shortcut notation $w_{0i_{0}}\equiv w_{i_{0}}(t=0)$. In the right
hand side of this expression, the expectation value $\langle\delta
b_{n_{1}}(t_{1})\dots\delta b_{n_{p}}(t_{p})\rangle$ can be related, in the
Gaussian regime, to all permutations of products of contraction of pairs by
using the Wick theorem,

$\langle\delta b_{n_{1}}(t_{1})\dots\delta
b_{n_{p}}(t_{p})\rangle=\sum_{\mathrm{pairings}}\prod_{j<\ell}\langle\delta
b_{n_{j}}(t_{j})\delta b_{n_{\ell}}(t_{\ell})\rangle,$ (8)

where the notation $\sum_{\mathrm{pairings}}\prod_{j<\ell}$ stands for the
$(2n_{p}-1)!!$ possible permutations of pairs with $t_{j}<t_{\ell}$. Without
loss of generality, we consider in this study the case of a 3D isotropic
turbulence without helicity. The 2-pt function is then dependent on the time
difference only,

$\langle\delta b_{n_{j}}(t_{j})\delta
b_{n_{\ell}}(t_{\ell})\rangle=\frac{\varphi(t_{j}-t_{\ell})}{3}\delta_{n_{j}n_{\ell}},$
(9)

with $\varphi(t)$ a function that describes the correlation of the turbulence
experienced by a test-particle along its path at two different times. An
expression for $\varphi(t)$, inferred from a formal derivation going beyond
the white-noise approximation or the red-noise one (see Appendix C), will be
presented in Section III.

On inserting Eqn. 8 and Eqn. 9 into Eqn. 7, the Dyson series reads as

$\langle
w_{i_{0}}(t)\rangle=w_{0i_{0}}+\sum_{p=1}^{\infty}\left(\frac{\delta\Omega^{2}}{3}\right)^{p}w_{0i_{2p}}\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\dots\int_{0}^{t_{2p-1}}\mathrm{d}t_{2p}\sum_{\mathrm{pairings}}\prod_{j<\ell}\\\
\left(R_{i_{j-1}k_{j}}(\Omega_{0}t_{j})R_{i_{\ell-1}k_{j}}(\Omega_{0}t_{\ell})R^{-1}_{m_{j}i_{j}}(\Omega_{0}t_{j})R^{-1}_{m_{j}i_{\ell}}(\Omega_{0}t_{\ell})-R_{i_{j-1}k_{j}}(\Omega_{0}t_{j})R_{i_{\ell-1}k_{\ell}}(\Omega_{0}t_{\ell})R^{-1}_{k_{\ell}i_{j}}(\Omega_{0}t_{j})R^{-1}_{k_{j}i_{\ell}}(\Omega_{0}t_{\ell})\right)\varphi(t_{j}-t_{\ell}),$
(10)

which is the relevant equation to determine the time evolution of the
auxiliary variable $\mathbf{w}(t)$ and subsequently of the particle velocity
$\mathbf{v}(t)=\hat{R}^{-1}(\Omega_{0}t)\mathbf{w}(t)$. The various terms of
the expansion can be conveniently represented using diagrammatic rules [43,
13]. In the following, we denote the “mass propagator” function $w(t)$,
defined such that $w_{i}(t)=w(t)\hat{u}_{i}$, by a double line, while a single
line stands for the corresponding “free propagator” corresponding to
$w^{(0)}(t)=1$. On the other hand, considering a contraction of a pair as an
“interaction” in which a free propagator is inserted, $\langle\delta
b_{i_{1}}(t_{j_{1}})\delta b_{i_{2}}(t_{j_{2}})\rangle=\langle\delta
b_{i_{1}}(t_{j_{1}})u^{(0)}(t)\delta b_{i_{2}}(t_{j_{2}})\rangle$, a curved
dotted line connecting two “vertices” then stands for a time-ordered
integration over an average product of two stochastic fields. In this manner,
the first term of the summation ($p=1$) is generically represented as

$t,i_{0}$$t_{1}$$t_{2}$$0,i_{2}$$=\leavevmode\nobreak\
\displaystyle\frac{\delta\Omega^{2}}{3}\displaystyle\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\leavevmode\nobreak\
\bigg{(}R_{i_{0}k_{1}}(\Omega_{0}t_{1})R_{i_{1}k_{1}}(\Omega_{0}t_{2})R^{-1}_{m_{1}i_{1}}(\Omega_{0}t_{1})R^{-1}_{m_{1}i_{2}}(\Omega_{0}t_{2})$$-\displaystyle
R_{i_{0}k_{1}}(\Omega_{0}t_{1})R_{i_{1}k_{2}}(\Omega_{0}t_{2})R^{-1}_{k_{2}i_{1}}(\Omega_{0}t_{1})R^{-1}_{k_{1}i_{2}}(\Omega_{0}t_{2})\bigg{)}\varphi(t_{1}-t_{2}),$
(11)

which, using the expression of the rotation matrices, reads explicitly as

$t,z$$t_{1}$$t_{2}$$0,z$$=\leavevmode\nobreak\
\displaystyle-2\frac{\delta\Omega^{2}}{3}\displaystyle\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\leavevmode\nobreak\
\cos{\Omega_{0}(t_{1}-t_{2})}\varphi(t_{1}-t_{2}),$ (12)

in the case $i_{0}=i_{2}=z$. For $i_{0}=i_{2}=x$, on the other hand, the
contraction of the rotation matrices leads to

$t,x$$t_{1}$$t_{2}$$0,x$$=\leavevmode\nobreak\
\displaystyle-\frac{\delta\Omega^{2}}{3}\displaystyle\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\leavevmode\nobreak\
(1+\cos{\Omega_{0}(t_{1}-t_{2})})\varphi(t_{1}-t_{2}).$ (13)

It will prove useful to express diagrams in the Laplace reciprocal space. In
that case, indices referring to initial and final times $0$ and $t$ are
removed. The generic change of variables $t=x+x_{1}+\cdots+x_{p}$,
$t_{1}=x_{1}+\cdots+x_{p}$, $t_{2}=x_{2}+\cdots+x_{p}$,$\cdots$, $t_{p}=x_{p}$
allows for sending all integration boundaries between 0 and $+\infty$ for the
$x_{i}$ variables. In this manner, the diagram of Eqn. 12 for instance becomes

$z$$t_{1}$$t_{2}$$z$$=\leavevmode\nobreak\
\displaystyle-2\frac{\delta\Omega^{2}}{3}\displaystyle\int_{0}^{\infty}\mathrm{d}x\leavevmode\nobreak\
\mathrm{e}^{-sx}\int_{0}^{\infty}\mathrm{d}x_{2}\leavevmode\nobreak\
\mathrm{e}^{-sx_{2}}\int_{0}^{\infty}\mathrm{d}x_{1}\leavevmode\nobreak\
\mathrm{e}^{-sx_{1}}\cos{\Omega_{0}(x_{1})}\varphi(x_{1})$$=\leavevmode\nobreak\
\displaystyle-2\frac{\delta\Omega^{2}}{3}\mathcal{L}^{2}[1](s)\mathcal{L}[\varphi(x)\cos{\Omega_{0}x}](s),$
(14)

with $\mathcal{L}[f(x)](s)$ the Laplace transform of $f(x)$ expressed as a
function of the variable $s$. For $p=2$, in addition to the unconnected two-
loop contribution, there are nested and crossed diagrams. In the case
$i_{0}=i_{2}=z$ for instance, the contribution of the nested diagram is

$t,z$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$0,z$$=\displaystyle\leavevmode\nobreak\
4\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\displaystyle\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\int_{0}^{t_{2}}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4}\leavevmode\nobreak\
\cos{\Omega_{0}\left(\frac{t_{2}-t_{3}}{2}\right)}\cos{\Omega_{0}\left(t_{1}-\frac{t_{2}-t_{3}}{2}-t_{4}\right)}\varphi(t_{1}-t_{4})\varphi(t_{2}-t_{3}),$
(15)

while that of the crossed diagram is

$t,z$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$0,z$$\,\,=\displaystyle\leavevmode\nobreak\
2\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\displaystyle\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\int_{0}^{t_{2}}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4}\leavevmode\nobreak\
\cos{\Omega_{0}\left(t_{1}-t_{2}-t_{3}+t_{4}\right)}\varphi(t_{1}-t_{3})\varphi(t_{2}-t_{4}),$
(16)

and so on and so forth for higher values of $p$ and other values of $i_{2p}$
and $i_{0}$ (see Appendix B for $p=2$ contributions with other configurations
of $i_{0}$ and $i_{2p}$). In these conditions, the mass propagator results
from the infinite summation

$t,j$$0,i$$t,j$$0,i$$t,j$$0,i$$t,j$$0,i$$t,j$$0,i$$t,j$$0,i$ (17)

which can be symbolically summarized by introducing a blob that stands for the
sum over all the possible connected diagrams as

$t,j$$0,i$$t,j$$0,i$$t,j$$0,i$ (18)

##  III Partial-summation approximation of the 2-pt function $\langle\delta
b_{n_{j}}(t_{j})\delta b_{n_{\ell}}(t_{\ell})\rangle$

One key ingredient in the calculation of Eqn. 18 is the 2-pt function of the
fluctuating magnetic field $\langle\delta b_{n_{j}}(t_{j})\delta
b_{n_{\ell}}(t_{\ell})\rangle$ experienced by a test particle. Based on Eqn.
9, we focus in this section on deriving an expression for the function
$\varphi(t)$.

The fluctuating magnetic field, considered as a Gaussian random field with
zero mean and $\langle\delta B_{i}(\mathbf{x})\delta
B^{\star}_{j}(\mathbf{x^{\prime}})\rangle=\delta B^{2}\delta_{ij}$, is denoted
as $\delta\mathbf{B}(\mathbf{k})$ in the reciprocal Fourier space. For an
homogeneous turbulence, the 2-pt correlation function between two components
of $\delta\mathbf{B}(\mathbf{x})$ is invariant under spatial translations.
This implies that two Fourier components of the field are uncorrelated at
different wavenumber vectors:

$\langle\delta B_{i}(\mathbf{k})\delta
B^{\star}_{j}(\mathbf{k^{\prime}})\rangle=P_{ij}(\mathbf{k})\delta(\mathbf{k}-\mathbf{k^{\prime}}),$
(19)

where $P_{ij}$ is the spectral tensor defined as the Fourier transform of the
2-pt correlation function. Combining the solenoidal nature of the field, its
isotropy and its absence of helicity, $P_{ij}$ takes the form [44]

$P_{ij}(\mathbf{k})=\frac{\mathcal{E}(k)}{4\pi
k^{2}}\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right),$ (20)

with $\mathcal{E}(k)$ the kinetic energy spectrum of the turbulence, which,
for a Kolmogorov turbulence, is defined between $k_{\mathrm{min}}$ and
$k_{\mathrm{max}}$ as

$\mathcal{E}(k)=\frac{(2\pi)^{2/3}\delta
B^{2}}{3\left(L_{\mathrm{max}}^{2/3}-L_{\mathrm{min}}^{2/3}\right)}k^{-5/3}.$
(21)

The minimum wavenumber vector $\mathbf{k}_{\mathrm{min}}$ is related to the
distance $L_{\mathrm{max}}$ over which the correlation function is non-zero
(size of the largest “eddies”), while the maximum one
$\mathbf{k}_{\mathrm{max}}$ is related to the scale $L_{\mathrm{min}}$ at
which the dissipation rate of the turbulence overcomes the energy cascade
rate.

With these ingredients, the 2-pt function of the fluctuating magnetic field
$\langle\delta b_{n_{j}}(t_{j})\delta b_{n_{\ell}}(t_{\ell})\rangle$
experienced by a test particle can be expressed as

$\langle\delta b_{i}(t)\delta
b_{j}(0)\rangle\simeq\int_{\mathbf{k}_{\mathrm{min}}}^{\mathbf{k}_{\mathrm{max}}}\mathrm{d}\mathbf{k}\frac{\mathcal{E}(k)}{4\pi
k^{2}}\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right)\langle\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}(t)}\rangle,$
(22)

where the Corrsin approximation has been used [45]. Different approximations
have been proposed to estimate the factor $\langle
e^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}(t)}\rangle$ [38, 39, 19]. In this
study, to evaluate it, we start from the formal expansion

$\langle\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}(t)}\rangle=\sum_{n\geq
0}\frac{\mathrm{i}^{n}}{n!}\int_{0}^{t}\mathrm{d}t_{1}\leavevmode\nobreak\
\dots\int_{0}^{t}\mathrm{d}t_{n}\langle(\mathbf{k}\cdot\mathbf{v}(t_{1}))\dots(\mathbf{k}\cdot\mathbf{v}(t_{n}))\rangle,$
(23)

where the substitution $\mathbf{x}(t)=\int_{0}^{t}\leavevmode\nobreak\
dt^{\prime}\leavevmode\nobreak\ \mathbf{v}(t^{\prime})$ has been used. Next,
the $n$-point correlation function entering into the integrand expression can
be substituted for the sum of all possible contraction of pairs (Wick
theorem). Because the process that draws at random the wavenumber vectors of
the turbulence is independent from that governing the velocity decorrelation
of the test particles, each pair can be approximated as

$\langle(\mathbf{k}\cdot\mathbf{v}(t_{1}))(\mathbf{k}\cdot\mathbf{v}(t_{2}))\rangle\simeq(kc)^{2}\langle\cos{\mathbf{\hat{k}}\cdot\mathbf{\hat{v}}(t_{1})}\cos{\mathbf{\hat{k}}\cdot\mathbf{\hat{v}}(t_{2})}\rangle.$
(24)

The arguments of the cosines are the pitch angles between the particle
velocities and the wavenumber vectors of the turbulence. The decorrelation of
the pitch angle is then assumed to decay exponentially,

$\langle(\mathbf{k}\cdot\mathbf{v}(t_{1}))(\mathbf{k}\cdot\mathbf{v}(t_{2}))\rangle\simeq(kc)^{2}\mathrm{e}^{-(t_{1}-t_{2})/\xi(k)},$
(25)

expression that requires the introduction of a correlation time scale $\xi$.
Guided by Monte-Carlo simulations that show a longer falloff timescale for the
2-pt function of the fluctuating magnetic field for large $k$ compared to
small $k$, a dependency in $(kc)^{-1}$ turns out to reproduce the main
features of $\varphi(t)$ for a reduced rigidity (Larmor radius conventionally
related to the turbulence only and expressed in units of the largest eddy
scale $L_{\mathrm{max}}$) $\rho=1$. An additional dependency in $\rho$ is
introduced through $\xi(k,\rho)=A\rho^{B}/(kc)$; $A\simeq 1$ and $B\simeq 0.5$
are found to provide a good compromise to cover the gyroresonant and high-
rigidity regimes. Some algebra then leads to

$\langle\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}(t)}\rangle\simeq\sum_{p\geq
0}(-(kc)^{2})^{p}\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\leavevmode\nobreak\
\dots\int_{0}^{t_{2p-1}}\mathrm{d}t_{2p}\\\
\sum_{\mathrm{pairings}}\prod_{\mathrm{pairs}\leavevmode\nobreak\
i<j}\mathrm{e}^{-(t_{i}-t_{j})/\xi(k)}.$ (26)

To evaluate the right hand side, only pairs with $j=i+1$ are retained. Under
this approximation, which corresponds to summing unconnected diagrams [43],
our estimate of
$\langle\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}(t)}\rangle$, denoted
with a subscript $0$, can be written in a compact non-linear manner:

$\langle\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}(t)}\rangle_{0}\simeq
1-(kc)^{2}\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\leavevmode\nobreak\
\mathrm{e}^{-\frac{t_{1}-t_{2}}{\xi(k)}}\langle\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}(t-t_{1})}\rangle_{0}.$
(27)

In the Laplace reciprocal space, the equation is then linear in
$\mathcal{L}[\langle\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}(t)}\rangle_{0}](s)$,
which reads as

$\mathcal{L}[\langle\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}(t)}\rangle_{0}](s)=\frac{1+s\xi(k)}{(1+s\xi(k))s+(kc)^{2}\xi(k)},$
(28)

so that
$\langle\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}(t)}\rangle_{0}$ can be
inferred from a numerical inverse Laplace transformation. Subsequently, we
infer the partial-summation approximation for the 2-pt function of the
fluctuating magnetic field experienced by a test particle as

$\langle\delta b_{i}(t)\delta
b_{j}(0)\rangle\simeq\frac{2\delta_{ij}}{3}\frac{(2\pi)^{2/3}\delta
B^{2}}{3\left(L_{\mathrm{max}}^{2/3}-L_{\mathrm{min}}^{2/3}\right)}\\\
\int_{k_{\mathrm{min}}}^{k_{\mathrm{max}}}\mathrm{d}k\leavevmode\nobreak\
k^{-5/3}\mathcal{L}^{-1}\bigg{[}\frac{1+s\xi(k)}{(1+s\xi(k))s+(kc)^{2}\xi(k)}\bigg{]}(t),$
(29)

expression from which the function $\varphi(t)$ is deduced by identification
with Eqn. 9.

Figure 1: 2-pt correlation function of the magnetic field experienced by a
test-particle as a function of the gyroperiod time scale $\delta\Omega t$ for
five different rigidities. The dotted lines are from Monte-Carlo simulations.

As benchmark values typical of those relevant for the propagation of high-
energy cosmic rays in the Galaxy, we use hereafter $\delta
B=1\leavevmode\nobreak\ \upmu$G, $L_{\mathrm{max}}/L_{\mathrm{min}}=100$ and
$L_{\mathrm{max}}=100\leavevmode\nobreak\ $pc. Examples of $\varphi(t)$
obtained from Eqn. 29 are shown as a function of the gyroperiod scale
$\delta\Omega\leavevmode\nobreak\ t$ as continuous lines in Fig. 1 for
different reduced rigidities, while the corresponding results from Monte-Carlo
simulations, taken from [46], are displayed as the dotted lines. Results are
displayed on a logarithmic scale for $\delta\Omega\leavevmode\nobreak\ t$ to
appreciate fine similarities and differences between the model and the Monte-
Carlo solution. As the rigidity increases, the time scale of correlation is
observed to decrease significantly in terms of $\delta\Omega t$. Overall, the
main features uncovered by the Monte-Carlo calculation are qualitatively well
reproduced by the model, especially at high rigidities (in particular the
oscillations around 0 that are absent in the quasi-linear theory or the red-
noise approximation). However, some differences are observed quantitatively.
The descent of the modeled $\varphi(t)$ with $\delta\Omega\leavevmode\nobreak\
t$ is delayed for $0.1\lesssim\rho\lesssim 1$ while, more importantly, it
becomes too fast for $\rho\lesssim 0.1$ and unable to capture the slow
decrease at large $\delta\Omega t$. Consequently, the use of the modeled
$\varphi(t)$ function is limited to the rigidity range such that $\rho\gtrsim
0.1$.

##  IV Summation scheme

To carry out a summation of the Dyson series (Eqn. 10 and Eqn. 18), we resort
to a two-step iteration procedure. In the first iteration, we calculate the
propagator that would be accurate in the case of a 2-pt function $\varphi(t)$
approximated by a Dirac function. As mentioned in the introduction, this
approximation proves to be accurate in the high-rigidity regime in pure
turbulence. It corresponds to the summation of the class of unconnected
diagrams and represents the simplest partial summation scheme of the Dyson
series [43]. However, this scheme leads to a non-physical solution in the case
of a non-zero mean field and in the gyroresonant regime. A more physical
solution is then obtained by including classes of nested and crossed diagrams
in the partial summation of the Dyson series, in which the first iteration of
the propagator is inserted into each ordered time. We now detail this scheme
below.

The propagator obtained by retaining only unconnected diagrams is denoted as
$\langle w_{i}(t)\rangle_{0}$. For the sake of clarity, we distinguish its
symbol from that of the final propagator by marking it with zigzags rather
than a double line. This propagator is the solution to the equation

$t,i$$t_{1}$$t_{2}$$0,j$$t,i$$0,j$$t,i$$0,j$ (30)

In view of the rules detailed in Section II governing the calculation of
diagrams, the second term of the equation reads as

$t,i$$t_{1}$$t_{2}$$0,j$$=\displaystyle\frac{\delta\Omega^{2}}{3}\displaystyle\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\leavevmode\nobreak\
\left(R^{-1}_{ij}(2\Omega_{0}(t_{1}-t_{2}))-(1+2\cos\Omega_{0}(t_{1}-t_{2}))R^{-1}_{ij}(\Omega_{0}(t_{1}-t_{2})\right)\varphi(t_{1}-t_{2})\langle
w_{i}(t-t_{1})\rangle_{0},$ (31)

which gives rise to a linear term for $\mathcal{L}[\langle
w_{i}(x)\rangle_{0}](s)$ in the Laplace reciprocal space:

$i$$t_{1}$$t_{2}$$j$$=\displaystyle\frac{\delta\Omega^{2}}{3}\mathcal{L}[\langle
w_{i}(x)\rangle_{0}](s)\mathcal{L}[1](s)\mathcal{L}\left[\left(R^{-1}_{ij}(2\Omega_{0}x)-(1+2\cos\Omega_{0}x)R^{-1}_{ij}(\Omega_{0}x)\right)\varphi(x)\right](s).$
(32)

The propagator $\langle w_{i}(t)\rangle_{0}$ is therefore obtained through a
numerical inverse of Laplace transform. Throughout this study, the Stehfest
algorithm is used, with Stehfest number $N=20$ [47].

The second summation scheme accounts, in addition to unconnected diagrams, for
contributions from nested and crossed diagrams:

$t,i$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$0,j$$t,i$$t_{1}$$t_{2}$$0,j$$t,i$$0,j$$t,i$$0,j$$\simeq$
(33)

which can be approximated by substituting internal double lines for zigzags:

$t,i$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$0,j$$t,i$$t_{1}$$t_{2}$$0,j$$t,i$$0,j$$t,i$$0,j$$\simeq$
(34)

In this manner, both the nested and crossed contributions give rise to linear
terms in $\mathcal{L}[\langle w_{i}(x)\rangle](s)$:

$i$$t_{1}$$t_{2}$$j$$=\leavevmode\nobreak\
\displaystyle\frac{\delta\Omega^{2}}{3}\mathcal{L}[\langle
w_{i}(x)\rangle](s)\mathcal{L}[1](s)\mathcal{L}\left[\left[R^{-1}_{ij}(2\Omega_{0}x)-(1+2\cos\Omega_{0}x)R^{-1}_{ij}(\Omega_{0}x)\right]\varphi(x)\langle
w_{i}(x)\rangle_{0}\right](s),$ (35)
$i$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$j$$\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ =\leavevmode\nobreak\
\displaystyle\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\mathcal{L}[\langle
w_{i}(x)\rangle](s)\mathcal{L}[1](s)\mathcal{L}\bigg{[}\big{[}1+R^{-1}_{ij}(\Omega_{0}(x_{1}+2x_{2}+x_{3}))(1+2\cos\Omega_{0}(x_{1}-x_{3}))$$-R^{-1}_{ij}(\Omega_{0}(2x_{1}+2x_{2}))-R^{-1}_{ij}(\Omega_{0}(2x_{2}+2x_{3}))\big{]}\varphi(x_{1}+x_{2})\varphi(x_{2}+x_{3})\langle
w_{i}(x_{1})\rangle_{0}\langle w_{i}(x_{2})\rangle_{0}\langle
w_{i}(x_{3})\rangle_{0}\bigg{]}(s).$ (36)

Note that the solution obtained for $\langle w_{i}(t)\rangle$ can be used as a
starting solution for iterating further based on Eqn. 34.

##  V Parallel diffusion

The velocity decorrelation function relevant for the parallel diffusion
corresponds to $i=j=z$ in equations of Section IV. Denoting for convenience as
$\hat{W}_{i}(s)$ the Laplace transform function $\mathcal{L}[\langle
w_{i}(t)\rangle](s)$, the first iterated propagator is inferred from Eqn. 30
in the Laplace space that reads as

$\hat{W}_{0z}(s)=\frac{1}{s}-\frac{2\delta\Omega^{2}}{3}\frac{\hat{W}_{0z}(s)}{s}\mathcal{L}\left[\varphi(x)\cos\Omega_{0}x\right](s),$
(37)

while the iterated propagator is inferred from Eqn. 34 as

$\hat{W}_{z}(s)=\frac{1}{s}-\frac{2\delta\Omega^{2}}{3}\frac{\hat{W}_{z}(s)}{s}\mathcal{L}\left[\varphi(x)\langle
w_{z}(x)\rangle_{0}\cos\Omega_{0}x\right](s)\\\
+2\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\frac{\hat{W}_{z}(s)}{s}\mathcal{L}\bigg{[}\varphi(x_{1}+x2)\varphi(x_{2}+x3)\\\
\langle w_{z}(x_{1})\rangle_{0}\langle w_{z}(x_{2})\rangle_{0}\langle
w_{z}(x_{3})\rangle_{0}\cos\Omega_{0}(x_{1}-x_{3})\bigg{]}(s).$ (38)

Figure 2: Parallel velocity decorrelation function of particles with reduced
rigidity $\rho=0.1$ for different values of $B_{0}$ ($\delta
B=1\leavevmode\nobreak\ \upmu$G), as a function of the gyroperiod time scale
$\delta\Omega t$. Dashed lines are from Monte-Carlo simulations.

The resulting velocity decorrelation functions $\langle
v_{0\parallel}v_{\parallel}(t)\rangle=\langle w_{0z}w_{z}(t)\rangle$ are shown
in Fig. 2 for $\rho=0.1$ and different values of $B_{0}$ (and $\delta
B=1\leavevmode\nobreak\ \upmu$G). For reference, results from Monte-Carlo
experiments are shown as the dashed lines. Overall, the main features of the
functions inferred from the simulations, namely the modulations on top of an
approximately exponential envelope that is decreasing slower with time for
increasing $B_{0}$ values, are captured by the calculation. In this rigidity
regime, the resonance between the Larmor radius of the particles with
wavelengths of the turbulence is the source of the modulations related to the
total angular frequency $\delta\Omega+\Omega_{0}$. They reflect memory effects
originating from large-scale wavenumber field lines around which particles
spiral while undergoing the imprint of a random walk caused by smaller
wavenumber vectors. As the intensity $B_{0}$ increases, the particles tend to
remain bound to the lines of the mean field for longer, and the decay takes
longer. Beyond similarities between the simulation and calculation results,
however, quantitative differences are observed in Fig. 2. The most notable one
concerns the global rate of falloff of the decorrelation functions that is
predicted to be too rapid for $B_{0}\lesssim 5\upmu$G by the calculation
compared to the Monte-Carlo simulations. The increase of the “decay time”
describing the approximately exponential envelope is indeed too slow for small
values of $B_{0}$, as clearly observed for $B_{0}=1\upmu$G, before to cross
the right range and to get too fast for $B_{0}\gtrsim 5\upmu$G. In other
words, the dependence of the decorrelation functions in $B_{0}$ is non-
linearly under-(over)estimated for $B_{0}\lesssim(\gtrsim)5\upmu$G. The
underestimation for small $B_{0}$, already visible for $B_{0}=0$, is
attributed at this stage on the one hand to the overestimation of $\varphi(t)$
observed in Fig. 1 for $\rho=0.1$, and on the other hand to an artefact due to
the partial summation of the Dyson series.

Figure 3: Same as Fig. 2 for $\rho=1$.

At higher rigidity ($\rho=1$), the Larmor radius of the particles is always
larger than the eddy sizes and scatterings can be considered independent one
from another. The process is Markovian and the decorrelation gets exponential
in the pure turbulence case, as demonstrated in [25] based on a white-noise
approximation to describe the 2-pt correlation function $\varphi(t)$. By
increasing $B_{0}$, the decay of the decorrelations gets slower. Similarly to
the case $\rho=0.1$, the calculation is observed to underestimate the “decay
time”; yet the values of $B_{0}$ leading to the underestimation span a much
wider range. Because the overestimation of $\varphi(t)$ observed in Fig. 1 for
$\rho=1$ is rather small, the velocity decorrelation function coincides quite
well with the Monte-Carlo one for $B_{0}=0$. Consequently, the differences for
$B_{0}>0$ stem predominantly from some incompleteness in the partial summation
of the Dyson series.

Figure 4: Parallel diffusion coefficient as a function of the reduced rigidity
for different values of $B_{0}$ ($\delta B=1\leavevmode\nobreak\ \upmu$G). The
Dashed line is from Monte-Carlo simulations in the pure turbulent case
($B_{0}=0)$.

The dependence in rigidity of the parallel diffusion coefficient
$D_{\parallel}$ as obtained from Eqn. 1 is shown in Fig. 4 for different
values of $B_{0}$, expressed in units of $cL_{\mathrm{max}}$. The dashed line
displays, for reference, the results obtained from Monte-Carlo simulations in
the case of pure turbulence. Despite the aforementioned differences between
the simulations and the calculation, $D_{\parallel}$ is observed to be
reproduced within a factor 2. In particular, the calculation slightly deviates
from the expected scalings in $\rho^{1/3}$ in the gyroresonant regime an in
$\rho^{2}$ in the quasi-ballistic regime [48, 49, 50]. As $B_{0}$ is
increasing, from the analysis of the decorrelation functions presented above,
$D_{\parallel}$ is expected to be more and more underestimated; yet the
calculation presented here provides genuine qualitative scalings that are also
reliable quantitatively around $\rho=0.1$.

##  VI Perpendicular diffusion

The velocity decorrelation function relevant for the perpendicular diffusion
corresponds to $i=j=x$ or $i=j=y$ in equations of Section IV. The calculation
proceeds the same way as for the parallel diffusion:

$\hat{W}_{0x}(s)=\frac{1}{s}-\frac{\delta\Omega^{2}}{3}\frac{\hat{W}_{0x}(s)}{s}\mathcal{L}\left[\varphi(x)\left(1+\cos{\Omega_{0}x}\right)\right](s),$
(39)

and

$\hat{W}_{x}(s)=\frac{1}{s}-\frac{\delta\Omega^{2}}{3}\frac{\hat{W}_{x}(s)}{s}\mathcal{L}\left[\varphi(x)\langle
w_{x}(x)\rangle_{0}(1+\cos\Omega_{0}x)\right](s)\\\
+2\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\frac{\hat{W}_{x}(s)}{s}\mathcal{L}\bigg{[}\varphi(x_{1}+x2)\varphi(x_{2}+x3)\\\
\langle w_{x}(x_{1})\rangle_{0}\langle w_{x}(x_{2})\rangle_{0}\langle
w_{x}(x_{3})\rangle_{0}\cos^{2}{\Omega_{0}\left(\frac{x_{1}}{2}+x_{2}+\frac{x_{3}}{2}\right)}\bigg{]}(s).$
(40)

The relevant decorrelation function is $\langle
v_{0\perp}v_{\perp}(t)\rangle=\langle
w_{0x}w_{x}(t)\rangle\cos\Omega_{0}t=\langle
w_{0y}w_{y}(t)\rangle\cos\Omega_{0}t$.

Figure 5: Perpendicular velocity decorrelation function of particles with
reduced rigidity $\rho=0.1$ for $\delta B=1\leavevmode\nobreak\ \upmu$G and
$B_{0}=1\upmu$G, as a function of the gyroperiod time scale $\delta\Omega t$.
Dashed line is from Monte-Carlo simulations.

A first illustration of the calculation is given in Fig. 5, where the
perpendicular decorrelation function is shown for $\rho=0.1$ and
$B_{0}=1\leavevmode\nobreak\ \upmu$G. The dashed line is from Monte-Carlo
simulations. As in the case of parallel diffusion in Section V, the envelope
of the decay is slightly underestimated by the calculation. However, the
modulations features, which have been shown from various simulations to be
responsible for the perpendicular sub-diffusive regime at early times, are
well reproduced. It is to be noted that none of the approximations proposed in
the literature could predict both the fast decaying envelope and the positions
of minimum and second maximum inherited from the modulations on top of the
decay.

Figure 6: Same as Fig. 5 for $\rho=1$ and $B_{0}=3\upmu$G. Figure 7:
Perpendicular diffusion coefficient as a function of the reduced rigidity for
different values of $B_{0}$ ($\delta B=1\leavevmode\nobreak\ \upmu$G). The
Dashed line is from Monte-Carlo simulations in the pure turbulent case
($B_{0}=0)$. Figure 8: Ratio between perpendicular and parallel diffusion
coefficient as a function of the reduced rigidity for different values of
$B_{0}$ ($\delta B=1\leavevmode\nobreak\ \upmu$G).

A second illustration is given in Fig. 6, for $\rho=0.1$ and
$B_{0}=1\leavevmode\nobreak\ \upmu$G. At such a high rigidity, the decay time
scale is longer than the gyroperiod time scale $\Omega_{0}t$, hence the
numerous oscillations. The results from the calculation or from the
simulations are almost indistinguishable.

From these two illustrations, we observe that the successes based on Eqn. 40
follow those presented in Section V in the case of parallel diffusion: the
calculation is able to reproduce the main features uncovered by numerical
simulations. However, the same limitations apply to the scaling of the results
with $B_{0}$, as a function of $\rho$.

In the same manner as in the parallel transport, the dependence in rigidity of
the perpendicular diffusion coefficient $D_{\perp}$ is shown in Fig. 7 for
different values of $B_{0}$, expressed in units of $cL_{\mathrm{max}}$. Due to
the oscillatory behavior of the velocity decorrelation functions, the
additional time integration should smooth out differences between the
simulations and the calculations. As $B_{0}$ is increasing, $D_{\perp}$ is
decreasing, as expected (for $B_{0}\rightarrow\infty$, a particle would be
spiraling around $\mathbf{B}_{0}$ at a fixed radius). More revealing is the
rigidity dependence of $D_{\perp}$, observed to rise more slowly than that of
$D_{\parallel}$ in the gyroresonant regime and, unlike $D_{\parallel}$, to
decrease in the high-rigidity regime. These dependencies are more clearly
highlighted in Fig. 8, where the ratio $D_{\perp}/D_{\parallel}$ is displayed.
The rise at low rigidities is in agreement with that revealed in Monte-Carlo
studies in which the turbulence dynamical range well covers the rigidities of
interest [33]. Furthermore, the decrease at high rigidities is also in
agreement with these simulations, as is the shift of the transition region
towards higher rigidities as $B_{0}$ is increasing.

##  VII Anti-symmetric diffusion

Figure 9: Anti-symmetric velocity decorrelation function of particles with
reduced rigidity $\rho=0.1$ for $\delta B=1\leavevmode\nobreak\ \upmu$G and
$B_{0}=1\upmu$G, as a function of the gyroperiod time scale $\delta\Omega t$.
Dashed line is from Monte-Carlo simulations.

Finally, we provide the velocity decorrelation function relevant for the anti-
symmetric diffusion by illustrating the case $j=y,i=x$:

$\hat{W}_{0x}(s)=\frac{1}{s}-\frac{\delta\Omega^{2}}{3}\frac{\hat{W}_{0x}(s)}{s}\mathcal{L}\left[\varphi(x)\sin{\Omega_{0}x}\right](s),$
(41)

and

$\hat{W}_{x}(s)=\frac{1}{s}-\frac{\delta\Omega^{2}}{3}\frac{\hat{W}_{x}(s)}{s}\mathcal{L}\left[\varphi(x)\langle
w_{x}(x)\rangle_{0}\sin{\Omega_{0}x}\right](s)\\\
+2\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\frac{\hat{W}_{x}(s)}{s}\mathcal{L}\bigg{[}\varphi(x_{1}+x2)\varphi(x_{2}+x3)\langle
w_{x}(x_{1})\rangle_{0}\\\ \langle w_{x}(x_{2})\rangle_{0}\langle
w_{x}(x_{3})\rangle_{0}\left(1-\sin{\Omega_{0}\left(x_{1}+2x_{2}+x_{3}\right)}\right)\bigg{]}(s).$
(42)

The relevant decorrelation function is $\langle
v_{0\perp}v_{\perp}(t)\rangle_{\mathrm{antisym}}=\langle
w_{0x}w_{y}(t)\rangle\sin\Omega_{0}t=-\langle
w_{0y}w_{x}(t)\rangle\sin\Omega_{0}t$. Similarly with the perpendicular
velocity decorrelation case of Section VI, the envelope of the decay is
slightly underestimated by the calculation compared to the Monte-Carlo results
in the case $\rho=0.1$ and $B_{0}=1\leavevmode\nobreak\ \upmu$G.

Finally, the dependence in rigidity of the anti-symmetric diffusion
coefficient $D_{\mathrm{a.-s.}}$ is observed to be linear, while that in
$B_{0}$ is observed to be inversely proportional. Overall, the shorthand
expression

$\frac{D_{\mathrm{a.-s.}}}{cL_{\mathrm{max}}\rho}\simeq\frac{\delta B}{B_{0}}$
(43)

is enough to capture the dependencies in the range of $B_{0}$ studied,
$0.1\leq\delta B/B_{0}\leq 1$.

##  VIII Conclusion

Velocity decorrelation functions of high-energy cosmic rays propagating in
magnetic fields have been obtained from the Dyson series governing the motion
of the particles. The absence of small parameter of expansion forbids any
perturbation theory to hold, and hence any truncation of the Dyson series
[51]. The partial-summation scheme developed in this study is shown to provide
an approximate solution that captures the main features uncovered by numerical
simulations in a range of rigidities covering the transition between the
gyroresonant and the quasi-ballistic regimes on the one hand, and the quasi-
ballistic regime itself on the other hand. Keeping in mind the limitations
underlined in terms of underestimation of the diffusion coefficients, the
calculation based on Eqn. 10 approximated by Eqn. 34 provides a rapid tool for
deriving approximate solutions without having to resort to heavy numerical
simulation campaigns.

The calculation was illustrated in the case of 3D isotropic turbulence
following a Kolmogorov power spectrum. This turbulence model is particularly
useful as a reference, as it has been widely used in the literature as a case
study for Monte Carlo simulations of test particles. However, the
computational techniques presented in this article are not limited to this
particular turbulence. Several astrophysical contexts require a 1D “slab”
turbulence model, a 2D isotropic one, 3D anisotropic ones with or without
helicity [52, 53, 54]. All these turbulence are described by different
spectral tensors $P_{ij}$ compared to the one used in Eqn. 20 that would give
specific expressions for the 2-pt function of the magnetic field compared to
Eqn. 9. A comprehensive study of these turbulence models is planned in a
future study.

A key ingredient of the approximate solution relies on the modeling of the
2-pt correlation function of the turbulence experienced by the particles
between two successive times. The approximation proposed in this study is
found to reproduce the main features uncovered, here again, by numerical
simulations as long as the reduced rigidity is larger than about one tenth of
the largest scale of the turbulence. Progress in the modeling of this 2-pt
function is needed to extend to lower rigidity the range of validity of the
type of calculation presented in this study.

###### Acknowledgements.

I gratefully thank Patrick Reichherzer and Leander Schlegel for providing me
the Monte-Carlo simulations presented in Fig. 1, and acknowledge funding from
ANR via the grant MultI-messenger probe of Cosmic Ray Origins (MICRO),
ANR-20-CE92-0052. This work was also made possible with the support of the
Institut Pascal at Université Paris-Saclay during the Paris-Saclay
Astroparticle Symposium 2022, with the support of the P2IO Laboratory of
Excellence (program “Investissements d’avenir” ANR-11-IDEX-0003-01 Paris-
Saclay and ANR-10-LABX-0038), the P2I axis of the Graduate School Physics of
Université Paris-Saclay, as well as IJCLab, CEA, IPhT, APPEC, the IN2P3 master
projet UCMN and EuCAPT ANR-11-IDEX-0003-01 Paris-Saclay and ANR-10-LABX-0038).
Finally, I thank Haris Lyberis for his help on an earlier development of this
work.

##  Appendix A Monte-Carlo generator

To serve as a reference for testing the model, a Monte-Carlo estimation of the
velocity correlation function is used. Following a scheme similar to that
widely used in the literature, a large number of particle trajectories in
given magnetic field configurations is simulated by solving numerically Eqn.
3. The numerical integration is performed using the standard Runge-Kutta
integrator.

The regular field is oriented along the $z$ axis,
$\mathbf{B}_{0}=B_{0}\mathbf{u}_{z}$, with $B_{0}$ constant in a given
configuration. On top of this constant field, a turbulence $\delta\mathbf{B}$
is added. To simulate numerically an isotropic and spatially homogeneous
turbulent field, an algorithm similar to that in [44, 18] is used. The recipe
consists in summing a large number $N_{\mathrm{m}}$ of plane waves
($N_{\mathrm{m}}=250$ in this study) with corresponding wave vector
$\mathbf{k}_{n}$, the direction, phase $\phi_{n}$ and polarisation of which
are chosen randomly:

$\delta\mathbf{B}(\mathbf{x})=\sqrt{2}\sum_{n=1}^{N_{\mathrm{m}}}\sum_{\alpha=1}^{2}\mathcal{E}_{n}(k_{n})\leavevmode\nobreak\
\mathbf{\hat{\xi}}_{n}^{\alpha}\leavevmode\nobreak\
\cos{(\mathbf{k}_{n}\cdot\mathbf{x}+\phi_{n}^{\alpha})}.$ (44)

To ensure the condition $\nabla\cdot\delta\mathbf{B}=0$, the two orthogonal
polarisation vectors $\mathbf{\hat{\xi}}_{n}^{\alpha}$ are oriented in the
plane perpendicular to the directions of the wave vectors. The wave number
distribution is built from a constant logarithmic spacing between
$k_{\mathrm{min}}$ and $k_{\mathrm{max}}$. The wave amplitudes satisfy
$\mathcal{E}_{n}^{2}(k_{n})=\mathcal{E}_{0}\delta
B^{2}k_{n}^{-5/3}(k_{n}-k_{n-1})$, where $\mathcal{E}_{0}$ is a normalisation
factor such that $\sum_{n}\mathcal{E}_{n}^{2}(k_{n})=\delta B^{2}$. In this
manner, the turbulence satisfies
$\langle\delta\mathbf{B}(\mathbf{x})\rangle=0$ and
$\langle\delta\mathbf{B}_{i}(\mathbf{x})\delta\mathbf{B}_{j}(\mathbf{x}^{\prime})\rangle=\delta
B^{2}\delta_{ij}$. The dynamic range of the turbulence explored here is
$L_{\mathrm{max}}/L_{\mathrm{min}}=100$.

##  Appendix B Perpendicular and anti-symmetric contributions for $p=2$

For completeness and clarity about the calculation rules of the diagrams,
explicit expressions of the perpendicular and anti-symmetric contributions are
provided in this Appendix in the case $p=2$.

Starting from the general expression of the crossed diagrams,

$t,i_{0}$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$0,i_{4}$$=\displaystyle\leavevmode\nobreak\
\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\displaystyle\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\int_{0}^{t_{2}}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4}\leavevmode\nobreak\
\varphi(t_{1}-t_{3})\varphi(t_{2}-t_{4})$$\times\displaystyle\leavevmode\nobreak\
\left(R_{i_{0}k_{1}}(\Omega_{0}t_{1})R_{i_{2}k_{1}}(\Omega_{0}t_{3})R^{-1}_{m_{1}i_{1}}(\Omega_{0}t_{1})R^{-1}_{m_{1}i_{3}}(\Omega_{0}t_{3})-R_{i_{0}k_{1}}(\Omega_{0}t_{1})R_{i_{2}k_{3}}(\Omega_{0}t_{3})R^{-1}_{k_{3}i_{1}}(\Omega_{0}t_{1})R^{-1}_{k_{1}i_{3}}(\Omega_{0}t_{3})\right)$$\times\displaystyle\leavevmode\nobreak\
\left(R_{i_{1}k_{2}}(\Omega_{0}t_{1})R_{i_{3}k_{2}}(\Omega_{0}t_{3})R^{-1}_{m_{2}i_{2}}(\Omega_{0}t_{1})R^{-1}_{m_{2}i_{4}}(\Omega_{0}t_{3})-R_{i_{1}k_{2}}(\Omega_{0}t_{1})R_{i_{3}k_{4}}(\Omega_{0}t_{3})R^{-1}_{k_{4}i_{2}}(\Omega_{0}t_{1})R^{-1}_{k_{2}i_{4}}(\Omega_{0}t_{3})\right),$
(45)

and using rotation-matrix properties such as
$R_{ij}(x_{1})R^{-1}_{jk}(x_{2})=R_{ik}(x_{1}-x_{2})$,
$R_{ij}(x)=R_{ji}(-x)=R^{-1}_{ij}(-x)$ and
$R_{ii}(x)=R^{-1}_{ii}(x)=1+2\cos\Omega_{0}x$, we obtain the expressions of
the perpendicular and anti-symmetric diagrams:

$t,y$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$0,y$$=$$t,x$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$0,x$$=\displaystyle\leavevmode\nobreak\
\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\displaystyle\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\int_{0}^{t_{2}}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4}\leavevmode\nobreak\
\varphi(t_{1}-t_{3})\varphi(t_{2}-t_{4})$$\times\displaystyle\leavevmode\nobreak\
\bigg{(}1-\cos{2\Omega_{0}(t_{1}-t_{3})}-\cos{2\Omega_{0}(t_{2}-t_{4})}+\cos{\Omega_{0}(t_{1}+t_{2}-(t_{3}+t_{4}))}[1+2\cos{\Omega_{0}(t_{2}+t_{3}-(t_{1}+t_{4}))}]\bigg{)},$
(46) $t,x$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$0,y$$=\leavevmode\nobreak\
-$$t,y$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$0,x$$=\displaystyle\leavevmode\nobreak\
\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\displaystyle\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\int_{0}^{t_{2}}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4}\leavevmode\nobreak\
\varphi(t_{1}-t_{3})\varphi(t_{2}-t_{4})$$\times\displaystyle\leavevmode\nobreak\
\bigg{(}-\sin{2\Omega_{0}(t_{1}-t_{3})}-\sin{2\Omega_{0}(t_{2}-t_{4})}+\sin{\Omega_{0}(t_{1}+t_{2}-(t_{3}+t_{4}))}[1+2\cos{\Omega_{0}(t_{2}+t_{3}-(t_{1}+t_{4}))}]\bigg{)}.$
(47)

In the Laplace reciprocal space, these diagrams read in terms of 3D numerical
integrations as

$y$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$y$$=$$x$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$x$$=\displaystyle\leavevmode\nobreak\
\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\mathcal{L}^{2}[1](s)\displaystyle\iiint\mathrm{d}x_{1}\mathrm{d}x_{2}\mathrm{d}x_{3}\leavevmode\nobreak\
\mathrm{e}^{-s(x_{1}+x_{2}+x_{3})}\varphi(x_{1}+x_{2})\varphi(x_{2}+x_{3})$$\times\displaystyle\leavevmode\nobreak\
\left[1+\left(1-2\cos{\Omega_{0}(x_{1}-x_{3})}+2\cos{\Omega_{0}x_{3}}\right)\cos{\Omega_{0}(x_{1}+2x_{2}+x_{3})}\right].$
(48) $x$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$y$$=\leavevmode\nobreak\
-$$y$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$x$$=\displaystyle\leavevmode\nobreak\
\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\mathcal{L}^{2}[1](s)\displaystyle\iiint\mathrm{d}x_{1}\mathrm{d}x_{2}\mathrm{d}x_{3}\leavevmode\nobreak\
\mathrm{e}^{-s(x_{1}+x_{2}+x_{3})}\varphi(x_{1}+x_{2})\varphi(x_{2}+x_{3})$$\times\displaystyle\leavevmode\nobreak\
\left[\left(1-2\cos{\Omega_{0}(x_{1}-x_{3})}+2\cos{\Omega_{0}x_{3}}\right)\sin{\Omega_{0}(x_{1}+2x_{2}+x_{3})}\right].$
(49)

The calculation of the nested diagrams proceeds the same manner. The general
expression reads as

$t,z$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$0,z$$=\displaystyle\leavevmode\nobreak\
\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\displaystyle\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\int_{0}^{t_{2}}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4}\leavevmode\nobreak\
\varphi(t_{1}-t_{4})\varphi(t_{2}-t_{3})$$\times\displaystyle\leavevmode\nobreak\
\left(R_{i_{0}k_{1}}(\Omega_{0}t_{1})R_{i_{3}k_{1}}(\Omega_{0}t_{4})R^{-1}_{m_{1}i_{1}}(\Omega_{0}t_{1})R^{-1}_{m_{1}i_{4}}(\Omega_{0}t_{4})-R_{i_{0}k_{1}}(\Omega_{0}t_{1})R_{i_{3}k_{4}}(\Omega_{0}t_{4})R^{-1}_{k_{4}i_{1}}(\Omega_{0}t_{1})R^{-1}_{k_{1}i_{4}}(\Omega_{0}t_{4})\right)$$\times\displaystyle\leavevmode\nobreak\
\left(R_{i_{1}k_{2}}(\Omega_{0}t_{2})R_{i_{2}k_{2}}(\Omega_{0}t_{3})R^{-1}_{m_{2}i_{2}}(\Omega_{0}t_{2})R^{-1}_{m_{2}i_{3}}(\Omega_{0}t_{3})-R_{i_{1}k_{2}}(\Omega_{0}t_{2})R_{i_{2}k_{3}}(\Omega_{0}t_{3})R^{-1}_{k_{3}i_{2}}(\Omega_{0}t_{2})R^{-1}_{k_{2}i_{3}}(\Omega_{0}t_{3})\right),$
(50)

which leads to

$t,y$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$0,y$$=$$t,x$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$0,x$$=\displaystyle\leavevmode\nobreak\
\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\displaystyle\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\int_{0}^{t_{2}}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4}\leavevmode\nobreak\
\varphi(t_{1}-t_{4})\varphi(t_{2}-t_{3})$$\times\displaystyle\leavevmode\nobreak\
\bigg{(}\cos{2\Omega_{0}(t_{1}+t_{3}-(t_{2}+t_{4}))}-\cos{\Omega_{0}(2(t_{1}-t_{4})+t_{3}-t_{2})}[1+2\cos{\Omega_{0}(t_{2}-t_{3})}]$$-\cos{\Omega_{0}(t_{1}-t_{4})}[1+2\cos{\Omega_{0}(2(t_{3}-t_{2})+t_{1}-t_{4})}]+\cos{\Omega_{0}(t_{1}-t_{4})}[1+2\cos{\Omega_{0}(t_{1}+t_{3}-(t_{2}+t_{4}))}][1+2\cos{\Omega_{0}(t_{2}-t_{3})}]\bigg{)},$
(51)

and to

$t,x$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$0,y$$=\leavevmode\nobreak\
-$$t,y$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$0,x$$=\displaystyle\leavevmode\nobreak\
\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\displaystyle\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\int_{0}^{t_{2}}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4}\leavevmode\nobreak\
\varphi(t_{1}-t_{4})\varphi(t_{2}-t_{3})$$\times\displaystyle\leavevmode\nobreak\
\bigg{(}\sin{2\Omega_{0}(t_{1}+t_{3}-(t_{2}+t_{4}))}-\sin{\Omega_{0}(2(t_{1}-t_{4})+t_{3}-t_{2})}[1+2\cos{\Omega_{0}(t_{2}-t_{3})}]$$-\sin{\Omega_{0}(t_{1}-t_{4})}[1+2\cos{\Omega_{0}(2(t_{3}-t_{2})+t_{1}-t_{4})}]+\sin{\Omega_{0}(t_{1}-t_{4})}[1+2\cos{\Omega_{0}(t_{1}+t_{3}-(t_{2}+t_{4}))}][1+2\cos{\Omega_{0}(t_{2}-t_{3})}]\bigg{)}.$
(52)

In the Laplace space, the diagrams read as

$y$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$y$$=$$x$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$x$$=\displaystyle\leavevmode\nobreak\
\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\mathcal{L}^{2}[1](s)\displaystyle\iiint\mathrm{d}x_{1}\mathrm{d}x_{2}\mathrm{d}x_{3}\leavevmode\nobreak\
\mathrm{e}^{-s(x_{1}+x_{2}+x_{3})}\varphi(x_{1}+x_{2}+x_{3})\varphi(x_{2})$$\times\displaystyle\leavevmode\nobreak\
\left(1+\cos{\Omega_{0}x_{2}}+\cos{\Omega_{0}(x_{1}+x_{3})}+\cos{\Omega_{0}(x_{1}+2x_{2}+x_{3})}\right),$
(53) $x$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$y$$=\leavevmode\nobreak\
-$$y$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$x$$=\displaystyle\leavevmode\nobreak\
\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\mathcal{L}^{2}[1](s)\displaystyle\iiint\mathrm{d}x_{1}\mathrm{d}x_{2}\mathrm{d}x_{3}\leavevmode\nobreak\
\mathrm{e}^{-s(x_{1}+x_{2}+x_{3})}\varphi(x_{1}+x_{2}+x_{3})\varphi(x_{2})$$\times\displaystyle\leavevmode\nobreak\
\left(\sin{\Omega_{0}x_{2}}+2\cos{\Omega_{0}x_{2}}\sin{\Omega_{0}(x_{1}+x_{2}+x_{3}})\right).$
(54)

##  Appendix C Red-noise approximation

Figure 10: Parallel velocity decorrelation function of particles with reduced
rigidity $\rho=1$ for $\delta B=1\leavevmode\nobreak\ \upmu$G and different
values of $B_{0}$ as a function of the gyroperiod time scale $\delta\Omega t$,
as obtained from the red-noise approximation. Dashed lines are from Monte-
Carlo simulations. Figure 11: Same as Fig. 10 for $\rho=10$, $\delta
B=1\leavevmode\nobreak\ \upmu$G and $B_{0}=0$.

We present in this Appendix results in the context of modeling $\varphi(t)$ as
a red-noise process with parameter $\tau$,

$\varphi(t)=\exp{(-t/\tau)}.$ (55)

Thanks to the property of the exponential function
$\varphi(t_{1}+t_{2})=\varphi(t_{1})\varphi(t_{2})$, analytical results (in
the Laplace reciprocal space) can be obtained for $\hat{W}_{0i}(s)$. We
restrict the range of application of the results presented in this Appendix to
the quasi-ballistic regime (that is, $\rho\geq 0.5$) for which an heuristic
expression for $\tau$ is

$\tau(\rho)\simeq\frac{L_{\mathrm{max}}}{16\leavevmode\nobreak\ \rho c}.$ (56)

Treating the case of the parallel diffusion, the contribution of the
unconnected diagrams to the zigzag propagator reads as

$z$$t_{1}$$t_{2}$$z$$=\displaystyle\leavevmode\nobreak\
-2\frac{\delta\Omega^{2}}{3}\frac{\hat{W}_{0z}(s)}{s}\frac{\tau(1+s\tau)}{(1+s\tau)^{2}+(\Omega_{0}\tau)^{2}}.$
(57)

To explore the impact that an improved zigzag propagator might have, we also
consider additional classes of diagrams to approximate the mass propagator
that would be obtained by substituting the thin line between $t_{2}$ and
$t_{1}$ for a zigzag (summation of all nested diagrams). To do so, we add to
Eqn. 30 the following contributions:

$z$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$z$$=\leavevmode\nobreak\ \displaystyle
2\tau^{3}\frac{\delta\Omega^{2}}{3}\frac{\hat{W}_{0z}(s)}{s}\frac{8+12(s\tau)^{3}+2(s\tau)^{4}-11(\Omega_{0}\tau)^{2}-(\Omega_{0}\tau)^{4}-12s\tau(-2+(\Omega_{0}\tau)^{2})+s^{2}(26\tau^{2}-3\Omega_{0}^{2}\tau^{4})}{(2+s\tau)((1+s\tau)^{2}+(\Omega_{0}\tau)^{2})^{2}((2+s\tau)^{2}+(\Omega_{0}\tau)^{2})},$
(58)

$z$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$t_{5}$$t_{6}$$z$$=\leavevmode\nobreak\
\displaystyle
2\tau^{5}\frac{\delta\Omega^{2}}{3}\frac{\hat{W}_{0z}(s)}{s}\frac{1}{(3+s\tau)((1+s\tau)^{2}+(\Omega_{0}\tau)^{2})^{2}((2+s\tau)^{2}+(\Omega_{0}\tau)^{2})^{2}((3+s\tau)^{2}+(\Omega_{0}\tau)^{2})}$$\times\displaystyle\bigg{[}2(1+s\tau)^{2}(2+s\tau)^{2}(3+s\tau)^{2}+(\Omega_{0}\tau)^{2}(70+178s\tau+139(s\tau)^{2}$$\displaystyle+44(s\tau)^{3}+5(s\tau)^{4}+27(\Omega_{0}\tau)^{2}+20s\tau^{3}\Omega_{0}^{2}+4s^{2}\tau^{4}\Omega_{0}^{2}+(\Omega_{0}\tau)^{4})\bigg{]}.$
(59)

The expressions for the iterated propagator require only 1D integrations:

$z$$t_{1}$$t_{2}$$z$$=\leavevmode\nobreak\
\displaystyle-2\frac{\delta\Omega^{2}}{3}\frac{\hat{W}_{0z}(s)}{s}\mathcal{L}[\mathrm{e}^{-x/\tau}\langle
w_{0z}(x)\rangle_{0}\cos{\Omega_{0}x}](s),$ (60)
$z$$t_{1}$$t_{2}$$t_{3}$$t_{4}$$z$$=\leavevmode\nobreak\ \displaystyle
2\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\frac{\hat{W}_{0z}(s)}{s}\mathcal{L}[\mathrm{e}^{-2x/\tau}\langle
w_{0z}(x)\rangle_{0}](s)$$\leavevmode\nobreak\
\displaystyle\times\left(\mathcal{L}^{2}[\mathrm{e}^{-x/\tau}\langle
w_{0z}(x)\rangle_{0}\cos{\Omega_{0}x}](s)+\mathcal{L}^{2}[\mathrm{e}^{-x/\tau}\langle
w_{0z}(x)\rangle_{0}\sin{\Omega_{0}x}](s)\right).$ (61)

Results are displayed in Fig. 10 for $\rho=1$ and different values of $B_{0}$
(the case $B_{0}=100\leavevmode\nobreak\ \upmu$G is not shown as the
calculation yields to nonphysical results with regions such that $\langle
v_{0\parallel}v_{\parallel}(t)\rangle/c^{2}>1$). Also shown in Fig. 11 is the
decorrelation function obtained for $\rho=10$: the (exponential) decay is
observed to be $\sim 100$ times that obtained for $\rho=1$, yielding a
diffusion coefficient scaling effectively as $\rho^{2}$, as expected. We note
that including or not the classes of diagrams depicted in Eqn. 58 and Eqn. 59
changes the results very little; this reinforces the soundness of Eqn. 30 to
approximate the mass propagator during the zeroth iteration. Overall, the red-
noise approximation for the 2-pt function of the turbulence experienced by
particles between two successive times is therefore competitive in the high-
rigidity regime; yet the range applicability for $B_{0}$ is reduced.

## References

  * Pfrommer  _et al._ [2017] C. Pfrommer, R. Pakmor, C. M. Simpson, and V. Springel, Astrophys. J. Lett. 847, L13 (2017), arXiv:1709.05343 [astro-ph.GA] .
  * Butsky and Quinn [2018] I. S. Butsky and T. R. Quinn, The Astrophysical Journal 868, 108 (2018).
  * Ji  _et al._ [2020] S. Ji, T. Chan, C. B. Hummels, P. F. Hopkins, J. Stern, D. Keres, E. Quataert, C.-A. Faucher-Giguère, and N. Murray, Monthly Notices of the Royal Astronomical Society 496, 4221–4238 (2020).
  * Hopkins  _et al._ [2022] P. F. Hopkins, I. S. Butsky, G. V. Panopoulou, S. Ji, E. Quataert, C.-A. Faucher-Giguère, and D. Kereš, Mon. Not. Roy. Astron. Soc. 516, 3470 (2022), arXiv:2109.09762 [astro-ph.HE] .
  * Kulsrud and Pearce [1969] R. Kulsrud and W. P. Pearce, Astrophys. J. 156, 445 (1969).
  * Mannheim and Schlickeiser [1994] K. Mannheim and R. Schlickeiser, Astron. Astrophys. 286, 983 (1994), arXiv:astro-ph/9402042 .
  * Ensslin  _et al._ [2007] T. A. Ensslin, C. Pfrommer, V. Springel, and M. Jubelgas, Astron. Astrophys. 473, 41 (2007), arXiv:astro-ph/0603484 .
  * Guo and OH [2008] F. Guo and S. P. OH, Mon. Not. Roy. Astron. Soc. 384, 251 (2008), arXiv:0706.1274 [astro-ph] .
  * Squire  _et al._ [2021] J. Squire, P. F. Hopkins, E. Quataert, and P. Kempski, Monthly Notices of the Royal Astronomical Society 502, 2630–2644 (2021).
  * Ji  _et al._ [2022] S. Ji, J. Squire, and P. F. Hopkins, Monthly Notices of the Royal Astronomical Society 513, 282–295 (2022).
  * Farmer and Goldreich [2004] A. J. Farmer and P. Goldreich, Astrophys. J. 604, 671 (2004), arXiv:astro-ph/0311400 [astro-ph] .
  * Quataert  _et al._ [2021] E. Quataert, Y.-F. Jiang, and T. A. Thompson, Monthly Notices of the Royal Astronomical Society 510, 920–945 (2021).
  * Frisch [1966] U. Frisch, Annales d’Astrophysique 29, 645 (1966).
  * Skilling [1971] J. Skilling, Astrophys. J. 170, 265 (1971).
  * Jokipii [1972] J. R. Jokipii, Astrophys. J. 172, 319 (1972).
  * Bhatnagar  _et al._ [1954] P. L. Bhatnagar, E. P. Gross, and M. Krook, Phys. Rev. 94, 511 (1954).
  * Jones [1990] F. C. Jones, Astrophys. J. 361, 162 (1990).
  * Giacalone and Jokipii [1999] J. Giacalone and J. R. Jokipii, Astrophys. J. 520, 204 (1999).
  * Casse  _et al._ [2001] F. Casse, M. Lemoine, and G. Pelletier, Phys. Rev. D 65, 023002 (2001).
  * Candia and Roulet [2004] J. Candia and E. Roulet, JCAP 10, 007, arXiv:astro-ph/0408054 .
  * Parizot [2004] E. Parizot, Nucl. Phys. B Proc. Suppl. 136, 169 (2004), arXiv:astro-ph/0409191 .
  * De Marco  _et al._ [2007] D. De Marco, P. Blasi, and T. Stanev, in  _30th International Cosmic Ray Conference_ , Vol. 2 (2007) pp. 195–198, arXiv:0705.3184 [astro-ph] .
  * Globus  _et al._ [2008] N. Globus, D. Allard, and E. Parizot, Astron. Astrophys. 479, 97 (2008), arXiv:0709.1541 [astro-ph] .
  * Hauff  _et al._ [2010] T. Hauff, F. Jenko, A. Shalchi, and R. Schlickeiser, Astrophys. J. 711, 997 (2010).
  * Plotnikov  _et al._ [2011] I. Plotnikov, G. Pelletier, and M. Lemoine, Astron. Astrophys. 532, A68 (2011), arXiv:1105.0618 [astro-ph.HE] .
  * Harari  _et al._ [2014] D. Harari, S. Mollerach, and E. Roulet, Phys. Rev. D 89, 123001 (2014), arXiv:1312.1366 [astro-ph.HE] .
  * Fatuzzo and Melia [2014] M. Fatuzzo and F. Melia, Astrophys. J. 784, 131 (2014), arXiv:1402.5469 [astro-ph.HE] .
  * Sonsrettee  _et al._ [2015] W. Sonsrettee, P. Subedi, D. Ruffolo, W. H. Matthaeus, A. P. Snodin, P. Wongpan, and P. Chuychai, Astrophys. J. 798, 59 (2015).
  * Snodin  _et al._ [2016] A. P. Snodin, A. Shukurov, G. R. Sarson, P. J. Bushby, and L. F. S. Rodrigues, Mon. Not. Roy. Astron. Soc. 457, 3975 (2016), arXiv:1509.03766 [astro-ph.HE] .
  * Subedi  _et al._ [2017] P. Subedi _et al._ , Astrophys. J. 837, 140 (2017), arXiv:1612.09507 [physics.space-ph] .
  * Giacinti  _et al._ [2018] G. Giacinti, M. Kachelriess, and D. V. Semikoz, JCAP 07, 051, arXiv:1710.08205 [astro-ph.HE] .
  * Reichherzer  _et al._ [2020] P. Reichherzer, J. Becker Tjus, E. G. Zweibel, L. Merten, and M. J. Pueschel, Mon. Not. Roy. Astron. Soc. 498, 5051 (2020), arXiv:1910.07528 [astro-ph.HE] .
  * Dundovic  _et al._ [2020] A. Dundovic, O. Pezzi, P. Blasi, C. Evoli, and W. H. Matthaeus, Phys. Rev. D 102, 103016 (2020), arXiv:2007.09142 [astro-ph.HE] .
  * Reichherzer  _et al._ [2022] P. Reichherzer, L. Merten, J. Dörner, J. Becker Tjus, M. J. Pueschel, and E. G. Zweibel, Appl. Sciences 4, 15 (2022), arXiv:2104.13093 [astro-ph.HE] .
  * Reichherzer  _et al._ [2023] P. Reichherzer, A. F. A. Bott, R. J. Ewart, G. Gregori, P. Kempski, M. W. Kunz, and A. A. Schekochihin, (2023), arXiv:2311.01497 [astro-ph.HE] .
  * Fraschetti and Giacalone [2012] F. Fraschetti and J. Giacalone, Astrophys. J. 755, 114 (2012), arXiv:1206.6494 [astro-ph.HE] .
  * Forman and Gleeson [1975] M. A. Forman and L. J. Gleeson, Astrophys. Space Sci. 32, 77 (1975).
  * Jokipii [1966] J. R. Jokipii, Astrophys. J. 146, 480 (1966).
  * Jokipii [1973] J. R. Jokipii, Astrophys. J. 183, 1029 (1973).
  * Bieber and Matthaeus [1997] J. W. Bieber and W. H. Matthaeus, Astrophys. J. 485, 655 (1997).
  * Kubo [1957] R. Kubo, J. Phys. Soc. Jap. 12, 570 (1957).
  * Deligny [2021] O. Deligny, Astrophys. J. 920, 87 (2021), arXiv:2107.04391 [astro-ph.HE] .
  * Bourret [1962] R. C. Bourret, Il Nuovo Cimento XXVI, 3833 (1962).
  * Batchelor [1970] G. K. Batchelor,  _The theory of homogeneous turbulence_ (Cambridge University Press, 1970).
  * Corrsin [1959] S. Corrsin,  _Advances in Geophysics and Atmospheric Diffusion_ , Vol. 6 (New York: Academic, 1959).
  * Deligny  _et al._ [2023] O. Deligny, P. Reichherzer, and L. Schlegel, PoS ICRC2023, 509 (2023).
  * Stehfest [1970] H. Stehfest, Commun. ACM 13, 4749 (1970).
  * Wentzel [1974] D. G. Wentzel, Ann. Rev. Astron. Astrophys. 12, 71 (1974).
  * Bhattacharjee and Sigl [2000] P. Bhattacharjee and G. Sigl, Phys. Rept. 327, 109 (2000), arXiv:astro-ph/9811011 .
  * Aloisio and Berezinsky [2004] R. Aloisio and V. Berezinsky, Astrophys. J. 612, 900 (2004), arXiv:astro-ph/0403095 .
  * Kraichnan [1961] R. H. Kraichnan, Journal of Mathematical Physics 2, 124 (1961).
  * Goldreich and Sridhar [1995] P. Goldreich and S. Sridhar, Astrophys. J. 438, 763 (1995).
  * Cho and Lazarian [2002] J. Cho and A. Lazarian, Phys. Rev. Lett. 88, 245001 (2002).
  * Cho and Lazarian [2003] J. Cho and A. Lazarian, Monthly Notices of the Royal Astronomical Society 345, 325 (2003), arXiv:astro-ph/0301062 [astro-ph] .


# Synthetic gauge fields enable high-order topology on Brillouin real
projective plane

Jinbing Hu  Department of Physics, The University of Hong Kong, Pokfulam, Hong
Kong, China College of Optical-Electrical Information and Computer
Engineering, University of Shanghai for Science and Technology, Shanghai
200093, China Songlin Zhuang  College of Optical-Electrical Information and
Computer Engineering, University of Shanghai for Science and Technology,
Shanghai 200093, China Yi Yang<EMAIL_ADDRESS>Department of Physics, The
University of Hong Kong, Pokfulam, Hong Kong, China

###### Abstract

The topology of the Brillouin zone, foundational in topological physics, is
always assumed to be a torus. We theoretically report the construction of
Brillouin real projective plane ($\mathrm{RP}^{2}$) and the appearance of
quadrupole insulating phase, which are enabled by momentum-space nonsymmorphic
symmetries stemming from $\mathbb{Z}_{2}$ synthetic gauge fields. We show that
the momentum-space nonsymmorphic symmetries quantize bulk polarization and
Wannier-sector polarization nonlocally across different momenta, resulting in
quantized corner charges and an isotropic binary bulk quadrupole phase
diagram, where the phase transition is triggered by a bulk energy gap closing.
Under open boundary conditions, the nontrivial bulk quadrupole phase manifests
either trivial or nontrivial edge polarization, resulting from the violation
of momentum-space nonsymmorphic symmetries under lattice termination. We
present a concrete design for the $\mathrm{RP}^{2}$ quadrupole insulator based
on acoustic resonator arrays and discuss its feasibility in optics, mechanics,
and electrical circuits. Our results show that deforming the Brillouin
manifold creates opportunities for realizing high-order band topology.

The topology of Brillouin zone (BZ) is widely known as a torus, around which
the state evolution encodes the topology of solids and artificial materials.
However, broadly speaking, a torus is only one closed compact manifold and can
thus be transformed into other manifolds; this transition should lead to rich
topological consequences. To this end, the Brillouin torus was recently
transformed into a Klein bottle under a momentum-space non-symmorphic
($\mathbf{k}\\!-\\!{\mathrm{NS}}$) symmetry that performs glide reflection in
momentum space Chen  _et al._ (2022). In contrast to the conventional
nonsymmorphic symmetries that act on real space and feature symmetry-invariant
momenta (e.g. Refs. Wieder  _et al._ (2018); Lu  _et al._ (2016); Armitage
_et al._ (2018); Zhao and Schnyder (2016); Wieder  _et al._ (2016); Liu  _et
al._ (2014); Fang and Fu (2015)), $\mathbf{k}\\!-\\!{\mathrm{NS}}$ symmetries
act on momentum space and symmetry invariant momenta may not be identified
Chen  _et al._ (2022). The key to creating such
$\mathbf{k}\\!-\\!{\mathrm{NS}}$ symmetries is synthetic gauge fields
Aidelsburger  _et al._ (2018), which can projectively modify point-group
symmetries into space-group ones Zhao  _et al._ (2020); Chen  _et al._ (2022,
2023). Beyond the torus, only the Brillouin Klein bottle was synthesized under
a single $\mathbf{k}\\!-\\!{\mathrm{NS}}$ symmetry so far Chen  _et al._
(2022), leading to interesting nonlocal first-order topology.

The electric polarization of crystals Resta (1992); King-Smith and Vanderbilt
(1993); Vanderbilt (2018) was recently extended to higher electric multipole
moments, giving rise to topological multipole insulators characterized by
high-order band topology Benalcazar  _et al._ (2017a, b); Schindler  _et al._
(2018). Although first proposed in tight-binding models, topological
quadrupole insulators have been observed in various physical platforms (e.g.
Refs. Peterson  _et al._ (2018); Mittal  _et al._ (2019); Qi  _et al._ (2020);
Ni  _et al._ (2020); Xue  _et al._ (2020); Serra-Garcia  _et al._ (2018);
Imhof  _et al._ (2018)). In addition, the topological quadrupole phase was
also observed in nonsymmorphic systems, where quantized quadrupole moment is
protected by real-space nonsymmorphic reflections Lin  _et al._ (2020); Zhang
_et al._ (2020); Liu  _et al._ (2023); Mao  _et al._ (2022); Liu  _et al._
(2022). However, high-order band topology on manifolds beyond the Brillouin
torus has not been investigated yet.

In this work, we theoretically introduce a quadrupole topological insulator on
a two-dimensional Brillouin real projective plane ($\mathrm{RP}^{2}$) via
synthetic gauge fields. Two $\mathbf{k}\\!-\\!{\mathrm{NS}}$ symmetries,
emerged from $Z_{2}$ gauge fields and the associated checkerboard flux
pattern, partition the original BZ torus into several non-unique
$\mathrm{RP}^{2}$ manifolds, and any of them can be equivalently used for
topology diagnosis. The symmetries of $\mathbf{k}\\!-\\!{\mathrm{NS}}$ lead to
the quantization of bulk and Wannier-sector polarization nonlocally across
different momenta. An isotropic binary bulk phase diagram featured by bulk
energy gap closing is obtained, in contrast to the quadrupole insulator
associated with Wannier band closing Benalcazar  _et al._ (2017a); Peterson
_et al._ (2018); Mittal  _et al._ (2019); Qi  _et al._ (2020); Serra-Garcia
_et al._ (2018); Imhof  _et al._ (2018). For a nontrivial bulk quadrupole
protected by $\mathbf{k}\\!-\\!{\mathrm{NS}}$ symmetries, the edge
polarizations can still be trivial or non-trivial since lattice termination
unavoidably violates the $\mathbf{k}\\!-\\!{\mathrm{NS}}$ symmetry, indicating
a breakdown of the bulk-boundary correspondence. Finally, we present a
concrete design of the $\mathrm{RP}^{2}$ quadrupole insulator using acoustics
and discuss other possible realizations in optics, mechanics, and circuits.

Figure 1: Real projective plane in momentum space. a. Folding a real
projective plane (left) from a rectangle (right), whose opposite sides are
glued in an anti-parallel manner (denoted by color arrows). b. The
$\mathbf{k}\\!-\\!{\mathrm{NS}}$ mirror symmetries $M_{x}$ and $M_{y}$ divide
the first Brillouin zone (BZ) $[-\pi,\pi)\times[-\pi,\pi)$ into four quadrants
I-IV. c. The $\mathbf{k}\\!-\\!{\mathrm{NS}}$ inversion $P_{\pi}$ further
divides each quadrant into diagonal and anti-diagonal partners (e.g. green and
blue colors in the I and III quadrants. The red dots indicate two momenta
linked by $P_{\pi}$. d. $M_{x}$, $M_{y}$ and $P_{\pi}$ jointly divide the
first BZ into sixteen plaquettes. A reduced BZ can be defined using any four
plaquettes of distinct colors, e.g. the square at the zone center surrounded
by red and violet arrows. Other examples include the squares in black solid
and dashed lines.

Brillouin real projective plane. $\mathrm{RP}^{2}$ is a two-dimensional
surface generalizing the Möbius strip. It is constructed from a square by
identifying two pairs of opposite edges with a half twist (Fig. 1a).
Mathematically, the construction of $\mathrm{RP}^{2}$ can be represented as a
unit square ($[0,1]\times[0,1]$) with opposite edges following the equivalence
relations (see Sec. I of SM_ )

$\displaystyle(x,0)$ $\displaystyle\sim(1-x,1)\leavevmode\nobreak\
\leavevmode\nobreak\ \mathrm{for}\leavevmode\nobreak\ \leavevmode\nobreak\
0\leq x\leq 1,$ (1a) $\displaystyle(0,y)$
$\displaystyle\sim(1,1-y)\leavevmode\nobreak\ \leavevmode\nobreak\
\mathrm{for}\leavevmode\nobreak\ \leavevmode\nobreak\ 0\leq y\leq 1.$ (1b)

We aim to realize these conditions to realize $\mathrm{RP}^{2}$ in momentum
space. The $\mathbf{k}\\!-\\!{\mathrm{NS}}$ reflections $M_{x}$ and $M_{y}$
meet the need, as they constrain the wave vector $\textbf{k}=(k_{x},k_{y})$
via

$\displaystyle M_{x}:(k_{x},k_{y})$
$\displaystyle\rightarrow(-k_{x},k_{y}+\pi),$ (2a) $\displaystyle
M_{y}:(k_{x},k_{y})$ $\displaystyle\rightarrow(k_{x}+\pi,-k_{y}),$ (2b)

which realizes Eq. 1 for $(k_{x},k_{y})\in[-\pi/2,\pi/2]\times[-\pi/2,\pi/2]$
(see Sec. II of SM_ ). In addition, the simultaneous presence of $M_{x}$ and
$M_{y}$ naturally defines a $\mathbf{k}\\!-\\!{\mathrm{NS}}$ inversion
$P_{\pi}\equiv M_{x}M_{y}$, which confines the wavevectors as

$\displaystyle P_{\pi}:(k_{x},k_{y})\rightarrow(\pi-k_{x},\pi-k_{y}).$ (3)

$M_{x}$, $M_{y}$, and $P_{\pi}$ further partition the Brillouin torus.
Specifically, $M_{x}$ and $M_{y}$ divides the first BZ
($[-\pi,\pi]\times[-\pi,\pi]$) into four quadrants (Fig. 1b). Further, each
quadrant is divided into two pairs by $P_{\pi}$: the diagonal and off-diagonal
pairs, and each pair is linked by $P_{\pi}$ (Fig. 1c). Therefore, $M_{x}$,
$M_{y}$, and $P_{\pi}$ jointly divide the first BZ is divided into sixteen
plaquettes (Fig. 1d) (see Sec. II of SM_ ). A reduced BZ can be defined on any
four plaquettes of distinct colors, e.g. the four plaquettes
$(k_{x},k_{y})\in[-\pi/2,\pi/2]\times[-\pi/2,\pi/2]$ at the zone center in
Fig. 1d. Notably, the boundaries of the reduced BZ connect in an anti-parallel
manner via $M_{x}$ (red arrows) and $M_{y}$ (violet arrows), thereby realizing
an $\mathrm{RP}^{2}$ in momentum space. There are three remarkable features of
the $\mathbf{k}\\!-\\!{\mathrm{NS}}$ symmetry group. First, analogous to their
real-space counterparts, the $\mathbf{k}\\!-\\!{\mathrm{NS}}$ symmetries
exhibit rich Abelian and non-Abelian algebra depending on the
(anti-)commutation between $M_{x}$ and $M_{y}$ Wieder  _et al._ (2018); Lu
_et al._ (2016); Zhang  _et al._ (2020); Wang  _et al._ (2016); Yang  _et al._
(2022). Second, the choice of the reduced $\mathrm{RP}^{2}$ BZ is non-unique
(e.g. squares in black solid and dashed lines in Fig. 1d, as long as four
plaquettes in distinct colors are enclosed). Third, the reduced
$\mathrm{RP}^{2}$ BZ remains a closed manifold, enabling adiabatic calculation
of topological invariants.

Figure 2: Lattice model featuring $\mathbf{k}\\!-\\!{\mathrm{NS}}$ reflections
$M_{x}$ and $M_{y}$. a-b. Threading $\pi$ flux projectively modify the
conventional reflections $m_{i}\leavevmode\nobreak\ (i=x,y)$ into
$\mathbf{k}\\!-\\!{\mathrm{NS}}$ ones $M_{i}\leavevmode\nobreak\ (i=x,y)$,
which anticommute with the translational symmetries $L_{j}\leavevmode\nobreak\
(j=x,y)$ along the other direction. c. Model with checkerboard $\pi$ flux
enables anti-commutative $M_{x}$ and $M_{y}$. Solid and dashed lines indicate
real positive and negative hoppings. d. Iso-energy contour plot at $E=0.5$ for
$\gamma_{x}=0.75$ and $\gamma_{y}=0.68$. e. The middle two bands in
$(\gamma_{x},\gamma_{y})$ space; red circle indicates gap closing. f. Bulk
spectrum as a function of $\gamma_{y}$ for $\gamma_{x}=0$ with red dots
showing gap closing at $\gamma_{y}=\pm\sqrt{2}$.

Tight-binding implementation. Consider the square lattice of four sites (Fig.
2a), the presence of $\pi$ flux enables an anti-commutative relation between
reflection $M_{x}$ and the translation operator $L_{y}$ along $y$, i.e.,
$M_{x}L_{y}M_{x}^{-1}L_{y}^{-1}=-1$. The right-hand negative sign can be
reformulated as a half translation along $k_{y}$, i.e., $-1=e^{iG_{y}a/2}$
(where $a$ is the lattice constant and $G_{y}$ the reciprocal lattice vector
along $y$) Chen  _et al._ (2022). Therefore, besides reversing $k_{x}$,
$M_{x}$ also contains a half translation along $k_{y}$, fulfilling Eq. (2a).
Similarly, the relation Eq. (2b) can be implemented by adding the $\pi$ flux
via a different gauge connection in Fig. 2b. Jointly, we arrange the negative
hopping along both $x$ and $y$ directions (Fig. 2a-b) to realize a
checkerboard flux pattern (Fig. 2c) featuring $\mathbf{k}\\!-\\!{\mathrm{NS}}$
$M_{x}$ and $M_{y}$. The associated spinless Bloch Hamiltonian is

$\displaystyle H(k_{x},k_{y})=\left[\begin{array}[]{cccc}0&a_{+}&b_{+}&0\\\
a_{+}^{*}&0&0&b_{-}\\\ b_{+}^{*}&0&0&a_{-}\\\ 0&b_{-}^{*}&a_{-}^{*}&0\\\
\end{array}\right],$ (8)

with $a_{\pm}=\gamma_{x}\pm\lambda_{x}e^{-ik_{x}}$,
$b_{\pm}=\pm\gamma_{y}+\lambda_{y}e^{-ik_{y}}$, and $*$ denoting complex
conjugate. Here, $\gamma_{x}$ and $\gamma_{y}$ ($\lambda_{x}$ and
$\lambda_{y}$) represent the intra-cell (inter-cell) hopping amplitudes along
$x$ and $y$, respectively. Without loss of generality, we set unity lattice
constants and $\lambda_{x}=\lambda_{y}=1$.

Under the link arrangements in Fig. 2c, the $\mathbf{k}\\!-\\!{\mathrm{NS}}$
symmetries satisfying Eqs. (2) and (3) have the form:
$M_{x}=\sigma_{3}\otimes\tau_{1}$, $M_{y}=\sigma_{1}\otimes\tau_{0}$, and
$P_{\pi}=i\sigma_{2}\otimes\tau_{1}$, together generating a non-Abelian group
$D_{8}$. Here, $\sigma$’s and $\tau$’s are Pauli matrices acting on sites
along $y$ and $x$, respectively. An example band structure validates the BZ
partition in Fig. 1d, where the band in the reduced $\mathrm{RP}^{2}$ BZ can
be extended to generate that of the full BZ (see Sec. III of SM_ ). Such
shrinking of the BZ results from the gauge transformation induced by the
negative coupling, which doubles the lattice constant along $x$ and $y$ (see
Sec. IV of SM_ ).

The model is an insulator at zero energy unless the bulk gap closing condition

$\displaystyle\sqrt{\gamma_{x}^{2}+\gamma_{y}^{2}}=\sqrt{\lambda_{x}^{2}+\lambda_{y}^{2}}$
(9)

is satisfied, as shown in Fig. 2e and its surface cut at $\lambda_{x}=0$ in
Fig. 2f [see Sec. III for a detailed derivation of Eq. (9)]. The isotropic
feature of Eq. (9) implies that the bulk-energy gap closing condition is non-
separable in the intra- and inter-cell hopping, which is helpful for
understanding the quadrupole phase below.

Figure 3:  Edge properties. a-d. Edge Wannier spectra $\\{\nu_{x}\\}$ and edge
polarization $p_{x}(R_{y})$ are calculated under the parameters
$(\gamma_{x},\gamma_{y})=(1.25,1.25)$ (a), $(0.25,1.25)$ (b), $(1.25,0.25)$
(c), and $(0.25,0.25)$ (d). e-f. Edge spectrum under $\gamma_{x}=0.25$,
$\gamma_{y}=0.5$ (e) where densities of edge (left) and bulk (right) bands
both exhibit momentum-space glide reflection.

Dipole moment quantized by $\mathbf{k}\\!-\\!{\mathrm{NS}}$ symmetries.
Although containing momentum half translation,
$\mathbf{k}\\!-\\!{\mathrm{NS}}$ $M_{x}$, $M_{y}$, and $P_{\pi}$ quantize
electric polarization (see Sec. V of SM_  for detailed derivation):

$\displaystyle\\{p_{x}^{j}(k_{y})\\}\xlongequal[]{M_{x}}\\{-p_{x}^{j}(k_{y}+\pi)\\}\leavevmode\nobreak\
\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ \textbf{1}.$ (10a)
$\displaystyle\\{p_{y}^{j}(k_{x})\\}\xlongequal[]{M_{y}}\\{-p_{y}^{j}(k_{x}+\pi)\\}\leavevmode\nobreak\
\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ \textbf{1}.$ (10b)
$\displaystyle\\{p_{x}^{j}(k_{y})\\}\xlongequal[]{P_{\pi}}\\{-p_{x}^{j}(\pi-
k_{y})\\}\leavevmode\nobreak\ \leavevmode\nobreak\
\mathrm{mod}\leavevmode\nobreak\ \textbf{1}.$ (10c)
$\displaystyle\\{p_{y}^{j}(k_{x})\\}\xlongequal[]{P_{\pi}}\\{-p_{y}^{j}(\pi-
k_{x})\\}\leavevmode\nobreak\ \leavevmode\nobreak\
\mathrm{mod}\leavevmode\nobreak\ \textbf{1}.$ (10d)

The quantization properties above are distinct from those of their non-
projective counterparts. For example, in Eq. (10a), $M_{x}$ quantizes only the
total polarization $p_{x}$ as 0 or 1/2 but not for each $k_{y}$, in contrast
to the conventional point-group reflection that quantizes polarization for
each $k_{y}$ Alexandradinata  _et al._ (2014); Benalcazar  _et al._ (2017b).
This feature stems from the fact that no wavevector k is invariant under
$M_{x}$, as indicated by the momentum half translation on the right hand of
Eq. (10a). In addition, $\mathbf{k}\\!-\\!{\mathrm{NS}}$ $P_{\pi}$ enforces a
vanished bulk polarization as the Berry phases of the two filled bands come in
pair $(-\nu,\nu)$.

The system supports quantized edge polarization. Figs. 3a-d show the Wannier
bands (top) and edge polarization $p_{x}^{\rm e}$ (bottom) of the $y$-ribbon
structure (i.e., periodic in $x$ and open in $y$) for
$(\gamma_{x},\gamma_{y})=(0.75\pm 0.5,0.75\pm 0.5)$, where the edge
polarization $p_{i}^{\rm e}$ is nontrivial for
$|\gamma_{i}|<|\lambda_{i}|\leavevmode\nobreak\ \leavevmode\nobreak\ (i=x,y)$.
The momentum glide reflection can be validated by the wave function. We
compute the edge spectrum of the $y$-ribbon for $\gamma_{x}=0.25$ and
$\gamma_{y}=0.5$ (Fig. 3e), where both the charge density of the edge bands
(red in Fig. 3e and Fig. 3f left) and the bulk bands (gray in Fig. 3e and Fig.
3f right) exhibit a half translation along $k_{x}$, as enforced by the $M_{y}$
symmetry.

Figure 4: Quadrupole insulator on RP2 from $\mathbf{k}\\!-\\!{\mathrm{NS}}$
symmetries. a-b. Gapped Wannier bands $\nu_{x}$ and $\nu_{y}$ showing vanished
bulk polarizations. c-d. Nontrivial Wannier-sector polarizations. e.
Appearance of four zero corner modes under full open boundaries for
parameters: $\gamma_{x}=0.25$, $\gamma_{y}=0.2$, and $N_{x}=N_{y}=25$ (number
of the unit cell along $x$ and $y$). f. Charge density of corner states (red
in Fig. 4e). g. Isotropic binary bulk phase diagram showing the quadrupole
$\mathbf{p}^{\nu}=(1/2,1/2)$ and the trivial phase. Their boundary (red
circle) is associated with a bulk gap closing. h. Edge polarization
$\mathbf{p}^{\mathrm{e}}$. Only the central square (orange) features
simultaneously nontrivial edge polarization
$(p^{\mathrm{e}}_{x},p^{\mathrm{e}}_{y})=(1/2,1/2)$ and corner modes. The four
segments (light gray regions), although featuring nontrivial bulk quadrupole
moments (see g), have nontrivial edge polarization only along a single
direction, and corner modes are absent therein. The boundaries between the
square and the segments are associated with edge energy gap closing (blue and
green lines; see Sec. VII of SM_ ).

Quadrupole moment quantized by $\mathbf{k}\\!-\\!{\mathrm{NS}}$ symmetries.
$M_{x}$ and $M_{y}$ further quantize the quadrupole moment nonlocally by
restricting the Wanner-sector polarization $p_{x}^{\nu_{y}}$ and
$p_{y}^{\nu_{x}}$ as (Sec. V of SM_ )

$\displaystyle\\{p_{x}^{\nu_{y}}(k_{y})\\}$
$\displaystyle\xlongequal[]{M_{x}}\\{-p_{x}^{\nu_{y}}(k_{y}+\pi)\\}\leavevmode\nobreak\
\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ \textbf{1},$ (11)
$\displaystyle\\{p_{y}^{\nu_{x}}(k_{x})\\}$
$\displaystyle\xlongequal[]{M_{y}}\\{-p_{y}^{\nu_{x}}(k_{x}+\pi)\\}\leavevmode\nobreak\
\leavevmode\nobreak\ \mathrm{mod}\leavevmode\nobreak\ \textbf{1}.$ (12)

Whereas $p_{x}^{\nu_{y}}$ and $p_{y}^{\nu_{x}}$ can be either 0 or 1/2, they
must additionally satisfy $p_{x}^{\nu_{y}}=p_{y}^{\nu_{x}}$ because the phase
transition occurs under an isotropic bulk band closing [see Eq. (9) and Fig.
2e]. This feature differentiates the present model from the originally
reported quadrupole model Benalcazar  _et al._ (2017a, b), where the phase
transition accompanies a Wannier band closing and the Wanner-sector
polarizations are independent along different directions. Shown in Figs. 4a-b
are the gapped Wannier bands symmetric around $\nu_{x,y}=0.5$ in the
$\mathrm{RP}^{2}$ reduced BZ, confirming vanished bulk polarization
$p_{x}=p_{y}=0$. Such vanished bulk polarization is enforced by
$\mathbf{k}\\!-\\!{\mathrm{NS}}$ symmetry $P_{\pi}$, accordant with the fact
that besides being symmetric around $k_{i}=0\leavevmode\nobreak\ (i=x,y)$, the
Wannier band is also symmetric around $k_{i}=\pm\pi/2\leavevmode\nobreak\
(i=x,y)$ (see Sec. VI of SM_ ). Figs. 4c and d show the Wannier-sector
polarizations $p_{x}^{\nu_{y}}=p_{y}^{\nu_{x}}=0.5$ for the filled Wannier
bands (red in Fig. 4a-b), confirming a nontrivial quadrupole moment.
Accordingly, the energy spectrum under full open boundaries (Fig. 4e) contains
four degenerate in-gap corner modes at zero energy, and the associated corner
charges are quantized as $\pm 1/2$ (Fig. 4f).

Fig. 4g shows the binary bulk phase diagram, where the nontrivial quadrupole
phase region resides inside the red circle and trivial elsewhere. The bulk
quadrupole phase transition is triggered by a bulk energy gap closing rather
than by Wannier gap closing, indicating a feature of the intrinsic high-order
phase Geier  _et al._ (2018); Queiroz  _et al._ (2019). In contrast, in the
edge polarization phase diagram (Fig. 4h), only the central square (orange)
shows nontrivial edge polarization along both directions (consequently, zero-
energy corner modes under full open boundaries), while the four segments
(gray) exhibit nontrivial edge polarization along one and only one direction.
The square and the segments are separated by a gap closing (vertical blue and
horizontal green lines) in the edge spectrum (see Sec. VII of SM_ ), a feature
of the extrinsic high-order phase where corner modes can be removed without
closing the bulk gap Geier  _et al._ (2018). Therefore, the bulk and edge
phase diagrams exhibit features of the intrinsic and extrinsic high-order
phases, respectively; such a violation of bulk-boundary correspondence relates
to the nonsymmorphic nature of the $\mathbf{k}\\!-\\!{\mathrm{NS}}$
symmetries, which can no longer be obeyed without periodicity.

Figure 5: Experimental designs. a. Acoustic links for positive (+) and
negative (-) hoppings via V-shaped and cross-shaped connections between two
identical resonators, respectively. b. Acoustic resonator array using the link
building blocks, where blue color denotes negative coupling. c. Bulk spectrum.
The bands of interest are colored in blue, between which a large gap
$\approx$700\text{\,}\mathrm{Hz}$$ (light gray) appears. d. Eigenfrequencies
of an open-boundary structure, exhibiting four in-gap modes. e. Pressure
distribution of one in-gap corner mode. f. Other feasible platforms: optical
ring resonators (top), mechanical metamaterials (middle), and circuits
(bottom).

Implementation of Brillouin real projective plane. We propose an experimental
design to realize the $\mathrm{RP}^{2}$ BZ and the associated high-order
topology based on acoustic resonator arrays Qi  _et al._ (2020); Xue  _et al._
(2022); Li  _et al._ (2022); Xue  _et al._ (2020); Ni  _et al._ (2020). The
cavity resonator emulating the site of Fig. 2c is designed to work at the
$p_{z}$ mode where both positive and negative couplings between nearest
neighboring resonators are possible Qi  _et al._ (2020); these are achieved by
linking two identical resonators in the cross and V-shaped manners,
respectively (Fig. 5a). In particular, a V-shaped (cross) connected bi-
resonator provides the negative (positive) coupling because the in-phase
coupled mode is at a lower (higher) frequency, while the out-of-phase coupled
mode is at a higher (lower) frequency (see Sec. VIII of SM_ ). The coupling
strength is strongly related to the tube cross-section and tube position with
respect to the middle of the resonator. Without loss of generality, we fix the
tube thickness and scan the tube width and position, from which distinct
coupling strengths can be designed. To get a wide band gap, a large coupling
contrast is utilized between intra-cell and inter-cell hoppings that are
$47\text{\,}\mathrm{Hz}$ and $482.2\text{\,}\mathrm{Hz}$, respectively.

Based on the link building blocks, we design the model (Fig. 2c) by
constructing an acoustic resonator array featuring the checkerboard $\pi$ flux
configuration in Fig. 5b, where the blue color indicates negative coupling.
The bulk energy spectrum is present in Fig. 5c, where the bands of interest
(blue) are separated by a large gap $\approx$700\text{\,}\mathrm{Hz}$$ (light
gray). We computed the eigenfrequency spectrum for an array of $5\times 6$
cells under the open boundary condition. The results are shown in Fig. 5d,
which reveals the presence of four in-gap corner modes, as verified by their
corresponding sound pressure distributions located at the corners of the array
(Fig. 5e). The same $\mathrm{RP}^{2}$-BZ model can be realized with coupled
optical ring resonators (Fig. 5f top), mechanical metamaterials (Fig. 5f
middle), and circuits (Fig. 5f bottom); we provide more detailed discussions
on these platforms in Sec. IX of SM_ .

In summary, we report the construction of a Brillouin real projective plane by
means of momentum-space nonsymmorphic ($\mathbf{k}\\!-\\!{\mathrm{NS}}$)
symmetries that are enabled by $\mathbb{Z}_{2}$ synthetic gauge fields.
Theoretical derivation reveals that $\mathbf{k}\\!-\\!{\mathrm{NS}}$
symmetries quantize bulk polarization and Wannier-sector polarization
nonlocally across different momenta, giving rise to a high-order quadrupole
phase with a phase transition of bulk energy gap closing. At last, we show
that the Brillouin real projective plane and its quadrupole phase can be
viably realized on many engineered physical platforms.

Note added. During the completion of this work about high-order topology on
Brillouin $\mathrm{RP}^{2}$, we noticed three relevant recent preprints Wang
_et al._ (2023); Zhu  _et al._ (2023); Tao  _et al._ (2023), which used
momentum-space nonsymmorphic symmetries for first-order and high-order
topology in two- and three-dimensional Brillouin Klein bottles.

## References

  * Chen  _et al._ (2022) Z. Chen, S. A. Yang, and Y. Zhao, Nature Communications 13, 2215 (2022).
  * Wieder  _et al._ (2018) B. J. Wieder, B. Bradlyn, Z. Wang, J. Cano, Y. Kim, H.-S. D. Kim, A. M. Rappe, C. Kane, and B. A. Bernevig, Science 361, 246 (2018).
  * Lu  _et al._ (2016) L. Lu, C. Fang, L. Fu, S. G. Johnson, J. D. Joannopoulos, and M. Soljačić, Nature Physics 12, 337 (2016).
  * Armitage  _et al._ (2018) N. Armitage, E. Mele, and A. Vishwanath, Reviews of Modern Physics 90, 015001 (2018).
  * Zhao and Schnyder (2016) Y. Zhao and A. P. Schnyder, Physical Review B 94, 195109 (2016).
  * Wieder  _et al._ (2016) B. J. Wieder, Y. Kim, A. Rappe, and C. Kane, Physical review letters 116, 186402 (2016).
  * Liu  _et al._ (2014) C.-X. Liu, R.-X. Zhang, and B. K. VanLeeuwen, Physical Review B 90, 085304 (2014).
  * Fang and Fu (2015) C. Fang and L. Fu, Physical Review B 91, 161105 (2015).
  * Aidelsburger  _et al._ (2018) M. Aidelsburger, S. Nascimbene, and N. Goldman, Comptes Rendus Physique 19, 394 (2018).
  * Zhao  _et al._ (2020) Y. Zhao, Y.-X. Huang, and S. A. Yang, Physical Review B 102, 161117 (2020).
  * Chen  _et al._ (2023) Z. Chen, Z. Zhang, S. A. Yang, and Y. Zhao, Nature Communications 14, 743 (2023).
  * Resta (1992) R. Resta, Ferroelectrics 136, 51 (1992).
  * King-Smith and Vanderbilt (1993) R. King-Smith and D. Vanderbilt, Physical Review B 47, 1651 (1993).
  * Vanderbilt (2018) D. Vanderbilt,  _Berry phases in electronic structure theory: electric polarization, orbital magnetization and topological insulators_ (Cambridge University Press, 2018).
  * Benalcazar  _et al._ (2017a) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Science 357, 61 (2017a).
  * Benalcazar  _et al._ (2017b) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Physical Review B 96, 245115 (2017b).
  * Schindler  _et al._ (2018) F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. Parkin, B. A. Bernevig, and T. Neupert, Science advances 4, eaat0346 (2018).
  * Peterson  _et al._ (2018) C. W. Peterson, W. A. Benalcazar, T. L. Hughes, and G. Bahl, Nature 555, 346 (2018).
  * Mittal  _et al._ (2019) S. Mittal, V. V. Orre, G. Zhu, M. A. Gorlach, A. Poddubny, and M. Hafezi, Nature Photonics 13, 692 (2019).
  * Qi  _et al._ (2020) Y. Qi, C. Qiu, M. Xiao, H. He, M. Ke, and Z. Liu, Physical Review Letters 124, 206601 (2020).
  * Ni  _et al._ (2020) X. Ni, M. Li, M. Weiner, A. Alù, and A. B. Khanikaev, Nature communications 11, 2108 (2020).
  * Xue  _et al._ (2020) H. Xue, Y. Ge, H.-X. Sun, Q. Wang, D. Jia, Y.-J. Guan, S.-Q. Yuan, Y. Chong, and B. Zhang, Nature communications 11, 2442 (2020).
  * Serra-Garcia  _et al._ (2018) M. Serra-Garcia, V. Peri, R. Süsstrunk, O. R. Bilal, T. Larsen, L. G. Villanueva, and S. D. Huber, Nature 555, 342 (2018).
  * Imhof  _et al._ (2018) S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, _et al._ , Nature Physics 14, 925 (2018).
  * Lin  _et al._ (2020) Z.-K. Lin, H.-X. Wang, Z. Xiong, M.-H. Lu, and J.-H. Jiang, Physical Review B 102, 035105 (2020).
  * Zhang  _et al._ (2020) X. Zhang, Z.-K. Lin, H.-X. Wang, Z. Xiong, Y. Tian, M.-H. Lu, Y.-F. Chen, and J.-H. Jiang, Nature communications 11, 65 (2020).
  * Liu  _et al._ (2023) Y. Liu, B. Lei, P. Yu, L. Zhong, K. Yu, and Y. Wu, Mechanical Systems and Signal Processing 198, 110433 (2023).
  * Mao  _et al._ (2022) N. Mao, H. Wang, Y. Dai, B. Huang, and C. Niu, npj Computational Materials 8, 154 (2022).
  * Liu  _et al._ (2022) Z. Liu, G. Wei, and J.-J. Xiao, Physical Review B 106, 085415 (2022).
  * (30) See Supplementary Materials for notes on (i) two-dimensional manifolds, (ii) momentum-space nonsymmorphic symmetries, (iii) hopping strength and bulk gap closing, (iv) gauge transformation, (v) quantization of Wannier-sector polarization on RP2, (vi) phase diagram of bulk quadrupole and edge polarization, and (vii) designs of $\mathrm{RP}^{2}$ BZ on engineered platforms.
  * Wang  _et al._ (2016) Z. Wang, A. Alexandradinata, R. J. Cava, and B. A. Bernevig, Nature 532, 189 (2016).
  * Yang  _et al._ (2022) Y. Yang, H. C. Po, V. Liu, J. D. Joannopoulos, L. Fu, and M. Soljačić, Physical Review B 106, L161108 (2022).
  * Alexandradinata  _et al._ (2014) A. Alexandradinata, X. Dai, and B. A. Bernevig, Physical Review B 89, 155114 (2014).
  * Geier  _et al._ (2018) M. Geier, L. Trifunovic, M. Hoskam, and P. W. Brouwer, Physical Review B 97, 205135 (2018).
  * Queiroz  _et al._ (2019) R. Queiroz, I. C. Fulga, N. Avraham, H. Beidenkopf, and J. Cano, Physical Review Letters 123, 266802 (2019).
  * Xue  _et al._ (2022) H. Xue, Z. Wang, Y.-X. Huang, Z. Cheng, L. Yu, Y. Foo, Y. Zhao, S. A. Yang, and B. Zhang, Physical Review Letters 128, 116802 (2022).
  * Li  _et al._ (2022) T. Li, J. Du, Q. Zhang, Y. Li, X. Fan, F. Zhang, and C. Qiu, Physical Review Letters 128, 116803 (2022).
  * Wang  _et al._ (2023) Y. Wang, C. Zhang, Z. Y. Chen, B. Liang, Y. X. Zhao, and J. Cheng, arXiv preprint arXiv:2305.07174 10.48550/arXiv.2305.07174 (2023).
  * Zhu  _et al._ (2023) Z. Zhu, L. Yang, J. Wu, Y. Meng, X. Xi, B. Yan, J. Chen, J. Lu, X. Huang, W. Deng, C. Shang, P. P. Shum, Y. Yang, H. Chen, G.-G. Liu, Z. Liu, and Z. Gao, arXiv preprint arXiv:2305.08450 10.48550/arXiv.2305.08450 (2023).
  * Tao  _et al._ (2023) Y.-L. Tao, M. Yan, M. Peng, Q. Wei, Z. Cui, S. A. Yang, G. Chen, and Y. Xu, arXiv preprint arXiv:2305.09174 10.48550/arXiv.2305.09174 (2023).