# Using Kernel-Based Statistical Distance to Study the Dynamics of Charged
Particle Beams in Particle-Based Simulation Codes

Chad E. Mitchell<EMAIL_ADDRESS>Robert D. Ryne  Lawrence Berkeley
National Laboratory, Berkeley, CA 94720, USA  Kilean Hwang  Facility for Rare
Isotope Beams, Michigan State University, East Lansing, MI 48824, USA

###### Abstract

Measures of discrepancy between probability distributions (statistical
distance) are widely used in the fields of artificial intelligence and machine
learning. We describe how certain measures of statistical distance can be
implemented as numerical diagnostics for simulations involving charged-
particle beams. Related measures of statistical dependence are also described.
The resulting diagnostics provide sensitive measures of dynamical processes
important for beams in nonlinear or high-intensity systems, which are
otherwise difficult to characterize. The focus is on kernel-based methods such
as Maximum Mean Discrepancy, which have a well-developed mathematical
foundation and reasonable computational complexity. Several benchmark problems
and examples involving intense beams are discussed. While the focus is on
charged-particle beams, these methods may also be applied to other many-body
systems such as plasmas or gravitational systems.

###### pacs:

##  I Introduction

When modeling the dynamics of charged particle beams, the following question
often arises. Given two ensembles of simulation particles, how similar are
they? In particular, when do two particle ensembles represent the same
underlying phase space density? This question is central to validating the
random sampling algorithm used for initial beam generation, to comparing
particle-based output across multiple simulation codes, to matching the
particle beam successfully into a periodic transport system, and for studying
the long-time phase space evolution of beams in circular and multi-pass
systems.

Two beams are typically compared using their first and second moments,
followed by qualitative examination of the beam phase space. For beams at high
intensity, in the presence of collective instabilities, or in the presence of
highly nonlinear transport, the details of the distribution (including higher-
order moments and nonlinear correlations) become increasingly important. Such
systems may exhibit filamentation of the beam phase space, developing
structure on finer and finer scales, in some cases relaxing to a quasi-
stationary state. It is then valuable to have quantitative particle-based
diagnostics that can characterize the nonlinear dynamical processes of regular
or chaotic mixing and collisionless relaxation [1, 2, 3, 4, 5, 6].

These problems can best be addressed by implementing a two-sample measure of
statistical distance with well-understood mathematical properties, which can
be used to compare particle populations. Statistical distances are widely used
in machine learning (ML), information theory, statistics, probability theory,
and data mining. Unfortunately, such quantities tend to have high
computational complexity. For beam physics applications, the complexity must
scale well with the number of particles $n$ (as $O(n^{2})$ or better) and with
the phase space dimension $d$ (up to at least $d=6$).

We describe the kernel-based statistical distance known as Maximum Mean
Discrepancy (MMD) [7], which has recently had a major impact in the ML
community. Broadly speaking, kernel methods allow nonlinear problems involving
higher-order statistics to be treated using linear methods, by embedding the
set of probability distributions into a reproducing kernel Hilbert space
(RKHS) [9, 8]. We do not discuss the RKHS formalism, although many results are
most naturally viewed in this context. The use of MMD as a statistical
distance leads naturally to a measure of statistical dependence or
correlation, known here as the Hilbert Schmidt Correlation (HSCor) [10]. These
two diagnostics can provide powerful quantitative tools to study the beam
dynamics processes mentioned above.

The structure of this paper is as follows. In Section II we review several
concepts of statistical distance. Section III describes the properties of
Maximum Mean Discrepancy and its implementation in particle-based tracking
codes. Section IV describes the properties of the Hilbert Schmidt Correlation
and its corresponding implementation. In Section V, we demonstrate how these
tools can be used for beam dynamics applications, using several benchmark
problems involving symplectic maps. In Section VI, we apply these tools to
examples involving high-intensity beams with self-consistent space charge.
Section VII contains a summary and conclusions. There are two Appendices.

##  II Statistical Distance

Consider a beam described by its kinetic distribution function $f$ on an
appropriate phase space of dimension $d$ with coordinates
$X=(q_{1},p_{1},\ldots,q_{d/2},p_{d/2})$, normalized so that

$\int f(X)dX=1.$ (1)

Then $f$ defines a probability density on the phase space. By Liouville’s
theorem, the condition (1) is preserved during the beam evolution, and we
often use $f_{t}$ to denote the phase space density at time $t$.

We wish to define a distance $\rho$ between pairs of probability densities $f$
and $g$, such that $\rho$ satisfies the following natural conditions:

$\displaystyle{\rm i})\text{ non-negativity: }\rho(f,g)\geq 0,$
$\displaystyle{\rm ii})\text{ symmetry: }\rho(f,g)=\rho(g,f),$
$\displaystyle{\rm iii})\text{ the triangle inequality:
}\rho(f,g)\leq\rho(f,h)+\rho(h,g)$ $\displaystyle\text{ for any probability
density $h$},$ $\displaystyle{\rm iv})\text{ identity: }\rho(f,g)=0\text{ if
and only if $f=g$ }$ $\displaystyle\text{(except on a possible set of zero
probability)}.$

That is, $\rho$ should define a metric on the set of all probability densities
[11].

The distance should also capture the concept of relaxation of beams, so that
$f_{t}$ relaxes to $f_{eq}$ in the “coarse-grained” sense if and only if
$\rho(f_{t},f_{eq})\rightarrow 0$ as $t\rightarrow\infty$. The concept of
relaxation in the “coarse-grained” sense is well-formalized by the
probabilistic concept of weak convergence [12, 13]: we say $f_{t}\Rightarrow
f_{eq}$ if and only if:

$\lim_{t\rightarrow\infty}\int f_{t}(X)\phi(X)dX=\int f_{eq}(X)\phi(X)dX$ (3)

for each bounded, continuous function $\phi$ on the phase space. Informally,
(3) states that the ensemble average of each well-behaved observable $\phi$
must approach the ensemble average of $\phi$ over $f_{eq}$ as
$t\rightarrow\infty$. This is a natural concept of convergence to use when
characterizing the long-time behavior of beams.

The description of the beam as a dynamically evolving probability density is
valid provided that the beam consists of a single particle species, without
charge loss. In the presence of charge loss, one is often interested in the
behavior of the part of the beam that lies within a bounded subregion of the
phase space. The tools described in this paper can also be adapted to treat
this case, using a technique to be described in Section VI.

In the following subsections, we briefly describe several indicators of
statistical distance that are used in ML and pattern recognition, information
science, probability, and statistics. For a general survey, see for example
[14].

###  II.1 Kullback-Liebler Divergence

The most widely-used statistical distance is the Kullback-Liebler (KL)
divergence $D_{KL}$, given by:

$D_{KL}(f||g)=\int f(X)\log\frac{f(X)}{g(X)}dX,$ (4)

which originated in information theory as a measure of the relative entropy of
one probability density with respect to another [15, 16]. It has the property
that $D_{KL}(f||g)\geq 0$, with equality if and only if $f=g$. (It satisfies
metric conditions i and iv above.) However, the integral in (4) is not defined
for all pairs of probability densities $f$ and $g$, and $D_{KL}$ fails to
satisfy conditions ii and iii, so it is not a metric in the above sense.

Note that $D_{KL}$ is invariant under any symplectic time-evolution map
$\mathcal{M}_{t}$, since:

$\displaystyle D_{KL}(f_{t}||g_{t})$ $\displaystyle=\int
f_{t}(X)\log\frac{f_{t}(X)}{g_{t}(X)}dX$ (5) $\displaystyle=\int
f_{0}(\mathcal{M}_{t}^{-1}(X))\log\frac{f_{0}(\mathcal{M}_{t}^{-1}(X))}{g_{0}(\mathcal{M}_{t}^{-1}(X))}dX$
$\displaystyle=\int
f_{0}(X^{\prime})\log\frac{f_{0}(X^{\prime})}{g_{0}(X^{\prime})}dX^{\prime}=D_{KL}(f_{0}||g_{0}),$

where we used the fact that $\mathcal{M}_{t}$ has Jacobian determinant 1. In
particular, if a density $f_{eq}$ is invariant under $\mathcal{M}_{t}$, in the
sense that

$f_{eq}(\mathcal{M}_{t}^{-1}(X))=f_{eq}(X),$ (6)

then $D_{KL}(f_{t}||f_{eq})$ is independent of $t$. Thus, $D_{KL}$ cannot
capture weak convergence of the form $f_{t}\Rightarrow f_{eq}$ to an invariant
density $f_{eq}$ under a symplectic time-evolution.

The KL divergence has been used in kinetic simulations [17], and it is well-
motivated by statistical mechanics considerations. However, typical algorithms
for computing $D_{KL}$ require binning the (possibly highly-filamented)
distribution functions $f$ and $g$, which becomes increasingly problematic in
a phase space of dimension $>2$. Although gridless two-sample estimation
algorithms also exist [18], the rate of convergence with particle number can
be arbitrarily slow, with a convergence rate that varies with the
distribution. This makes the quantity $D_{KL}$ difficult to estimate reliably
from samples, especially in spaces of high dimension.

###  II.2 Wasserstein Metric

The $p$-Wasserstein distance $(p=1,2,\ldots)$ is defined by:

$W_{p}(f,g)=\left(\min_{h}\int|X-Y|^{p}h(X,Y)dXdY\right)^{1/p}$ (7)

where $|X-Y|$ denotes the Euclidean distance between points $X$ and $Y$. Here
the minimum is taken [19] over all joint probability densities $h$ with
marginal densities $f$ and $g$, so that:

$f(X)=\int h(X,Y)dY,\quad g(Y)=\int h(X,Y)dX.$ (8)

The distances $W_{p}$ originated in the theory of optimal transport [22, 20,
21], where the case $p=1$ is also known as the Kantorovich-Rubinstein metric
or the Earth Mover’s Distance (EMD).

Note that (7) is guaranteed to be finite when the densities $f$ and $g$ both
have finite moments of order $p$. On the set of all such densities, $W_{p}$ is
known to satisfy all the metric conditions (i-iv). It is also known that
$W_{p}$ correctly captures the concept of weak convergence, in the sense that
$W_{p}(f_{t},f_{eq})\rightarrow 0$ (as $t\rightarrow\infty$) if and only if
$f_{t}\Rightarrow f_{eq}$ and the $p$th moments of $f_{t}$ converge to those
of $f_{eq}$.

Due to its desirable geometric properties, the Wasserstein distance has been
applied in ML to tasks such as shape matching, image retrieval, graphics, and
to the statistical analysis of detector events in high-energy colliders [23,
24, 25]. However, the estimation of $W_{p}$ from sample data requires solving
a linear optimization problem with a computational complexity of $O(n^{3}\log
n)$, where $n$ is the number of samples [26, 27]. Furthermore, the sample
estimate converges to the population value as $O(n^{-1/d})$ for $d>2$. This
makes $W_{p}$ challenging to use for practical beam dynamics simulations,
which typically require $n\geq 10^{5}$ and $d\geq 4$. Since there are few
cases in which (7) can be determined in closed form, algorithms for computing
$W_{p}$ are also difficult to benchmark.

###  II.3 Maximum Mean Discrepancy

A kernel $k$ is a symmetric, real-valued function defined on pairs of phase
space points that is positive definite, in the sense that:

$\sum_{i=1}^{N}\sum_{j=1}^{N}c_{i}c_{j}k(X_{i},X_{j})\geq 0,$ (9)

for each $N=1,2,\ldots$, each finite set of points $X_{1},\ldots,X_{N}$ and
real numbers $c_{1},\ldots,c_{N}$. Typical examples are provided in Appendix
A.

To each kernel $k$ is associated a Hilbert space (RKHS) consisting of real-
valued functions on the phase space. The Maximum Mean Discrepancy (MMD)
between two probability densities $f$ and $g$ is then defined by [7]:

$\gamma_{k}(f,g)=\max_{\phi}\left|\int f(X)\phi(X)dX-\int
g(X)\phi(X)dX\right|$ (10)

where the maximum is taken over all functions $\phi$ in the RKHS with
$||\phi||_{k}\leq 1$, $||\cdot||_{k}$ denotes the Hilbert space norm, and
$|\cdot|$ denotes the absolute value.

When $k$ is bounded, the quantity in (10) is defined for all probability
densities $f$ and $g$, and $\gamma_{k}$ satisfies the metric conditions
(i-iii). Additional restrictions on $k$ are used to ensure that $\gamma_{k}$
satisfies condition (iv), and that $\gamma_{k}$ captures weak convergence, in
the sense previously described. The class of kernels satisfying these
restrictions has been extensively studied [28, 29, 30, 31, 32], and it
includes most of the kernels widely used in ML, including those described in
Appendix A.

Due to its well-developed mathematical foundation, its connection to other ML
kernel methods such as support vector machines [8], its applicability to
domains more general than Euclidean space, and its relative ease of
estimation, the distance $\gamma_{k}$ has become a powerful tool in
statistical two-sample (homogeneity) testing for ML applications. Estimation
of (10) from sample data can be achieved using $O(n^{2})$ operations, and the
sample estimate converges to the population value as $O(n^{-1/2})$,
independently of the dimension $d$. Furthermore, approximations also exist
that can be computed with complexity $O(n)$, making $\gamma_{k}$ a practical
quantity for beam dynamics applications.

##  III Properties of Maximum Mean Discrepancy

Given a kernel $k$, the maximum appearing in (10) can be evaluated exactly by
using the properties of its corresponding RKHS. As a result, the MMD between
two probability densities $f$ and $g$ can be expressed using the explicit
integral formula [28]:

$\gamma_{k}(f,g)=\left(\iint
k(X,X^{\prime})\Delta(X)\Delta(X^{\prime})dXdX^{\prime}\right)^{1/2},$ (11)

where $\Delta=f-g$. For certain choices of $k$, the quantity (11) coincides
with other well-known indicators of statistical distance. For example, in the
special case that $k(X,X^{\prime})=|X|+|X^{\prime}|-|X-X^{\prime}|$, (11)
appears in the statistics literature as the energy distance [34, 33].

###  III.1 Basic properties

Given any kernel $k$, one may construct a corresponding kernel $k_{N}$ by:

$k_{N}(X,X^{\prime})=\frac{k(X,X^{\prime})}{\sqrt{k(X,X)k(X^{\prime},X^{\prime})}}.$
(12)

The condition that $k$ be positive definite (9) then implies that $k_{N}$ is
positive definite with $|k_{N}|\leq 1$. For simplicity, we will assume that
all kernels are so normalized. It then follows from (11) that the distance
$\gamma_{k}$ is dimensionless with:

$0\leq\gamma_{k}\leq 2.$ (13)

It is natural to choose a kernel that reflects the underlying properties of
the domain, so we often consider kernels that are translation-invariant, in
the sense that $k(X,X^{\prime})=k(X+\delta X,X^{\prime}+\delta X)$ for any
phase space displacement $\delta X$. A continuous, translation-invariant
kernel can be written in terms of its Fourier components as:

$k(X,X^{\prime})=\int e^{i(X-X^{\prime})\cdot\omega}\Lambda(\omega)d\omega.$
(14)

When $k$ is normalized, $\Lambda$ is a probability density [11] on the space
of frequencies $\omega=(\omega_{1},\ldots,\omega_{d})$. Using (14) in (11),
one finds that:

$\gamma_{k}(f,g)=\left(\int\left|\phi_{f}(\omega)-\phi_{g}(\omega)\right|^{2}\Lambda(\omega)d\omega\right)^{1/2},$
(15)

where $\phi_{f}$ and $\phi_{g}$ denote the Fourier transforms of the densities
$f$ and $g$, normalized so that:

$\phi_{f}(\omega)=\int e^{iX\cdot\omega}f(X)dX.$ (16)

When the probability density $\Lambda$ is also an integrable function that is
strictly positive everywhere, it is possible to prove that (15) satisfies the
metric conditions (i-iv) and correctly reflects the weak convergence of
probability distributions, as previously described.

One additional property of $\gamma_{k}$ is also useful. If
$\\{e_{l}:l=1,2,\ldots\\}$ denotes an orthonormal basis for the RKHS
associated with the kernel $k$, then we may define a complete sequence of beam
“moments” $m_{l}$ by:

$m_{l}(f)=\int e_{l}(X)f(X)dX\quad(l=1,2,\ldots).$ (17)

It follows from (11) that the MMD between two distributions $f$ and $g$ may be
written in terms of these moments as:

$\gamma_{k}(f,g)=\sqrt{\sum_{l=1}^{\infty}\left|m_{l}(f)-m_{l}(g)\right|^{2}}.$
(18)

In particular, for every $l=1,2,\ldots$ we have:

$\left|m_{l}(f)-m_{l}(g)\right|\leq\gamma_{k}(f,g).$ (19)

Thus when $\gamma_{k}$ is small, all of the moments $m_{l}$ of the
distributions $f$ and $g$ must nearly coincide. (An example is provided in
Appendix A.)

###  III.2 Sample estimate

A direct estimate of (11) from particle (sample) data is given by [7]:

$\displaystyle\gamma_{k}^{2}(f,g)$
$\displaystyle=\frac{1}{m^{2}}\sum_{i,j=1}^{m}k(X_{i},X_{j})-\frac{2}{mn}\sum_{i,j=1}^{m,n}k(X_{i},Y_{j})$
$\displaystyle+\frac{1}{n^{2}}\sum_{i,j=1}^{n}k(Y_{i},Y_{j}),$ (20)

where the $m$ particle phase space coordinates $\\{X_{j}\\}_{j=1}^{m}$ are
sampled from the distribution $f$, and the $n$ particle phase space
coordinates $\\{Y_{j}\\}_{j=1}^{n}$ are sampled from the distribution $g$.
Note that we allow $m\neq n$.

An alternative grouping of the sum (20) that yields superior numerical
performance in practice is given by:

$\gamma_{k}^{2}(f,g)=\sum_{i,j=1}^{m+n}c_{i}c_{j}k(\hat{X}_{i},\hat{X}_{j}),$
(21)

where

$c_{j}=\begin{cases}\frac{1}{m},&1\leq j\leq m\\\ -\frac{1}{n},&m+1\leq j\leq
m+n\end{cases}$ (22)

and $\\{\hat{X}_{j}\\}_{j=1}^{m+n}$ contains the phase space coordinates
sampled from the distribution $f$, followed by the phase space coordinates
sampled from $g$, so that:

$\hat{X}_{j}=\begin{cases}X_{j},&1\leq j\leq m\\\ Y_{j-m},&m+1\leq j\leq
m+n\end{cases}.$ (23)

In the form (21), it is clear from (9) that the estimate satisfies
$\gamma_{k}^{2}\geq 0$, and that the computational complexity is
$O((m+n)^{2})$.

In the special case that the kernel $k$ is translation-invariant, the
complexity can be reduced by using the spectral representation (15) to
approximate $\gamma_{k}$ as the sum [35, 36]:

$\gamma_{k}^{2}(f,g)=\frac{1}{L}\sum_{l=1}^{L}\left|\sum_{j=1}^{m+n}c_{j}e^{i\omega_{l}\cdot\hat{X}_{j}}\right|^{2}$
(24)

where the $\\{\omega_{l}\\}_{l=1}^{L}$ denote $L$ frequency vectors that are
randomly sampled from the spectral probablity density $\Lambda$, and the
quantities $\\{c_{j}\\}_{j=1}^{m+n}$ and $\\{\hat{X}_{j}\\}_{j=1}^{m+n}$ are
given by (22-23). The complexity of (24) is $O(L(m+n))$, where in practice it
is sufficient to use $L\ll m+n$. This results in significant speed-up when
translation-invariant kernels are used to compute $\gamma_{k}$.

###  III.3 Statistical error

It is shown in [26, 27] that the sample estimate (20) converges to the
population value (11) as $O(m^{-1/2}+n^{-1/2})$ when $m,n\rightarrow\infty$,
independently of the dimension of the underlying space. A similar result [35,
36] may be obtained for the estimate (24).

In practice, the sample estimate in (20-21) or (24) is used to test the
hypothesis that two distributions are distinct, $f\neq g$. To set a confidence
threshold for this test, one must know the statistical distribution of (20)
under the null hypothesis that $f=g$. This problem has been treated in detail
[7, 33]. For our purposes, it is enough to know that the rms value of
$\gamma_{k}$ under the hypothesis that $f=g$ is given by taking the expected
value of (20), which is given by:

$\gamma_{k}^{\rm
noise}=\operatorname{E}[\gamma_{k}^{2}]^{1/2}=\sqrt{\frac{m+n}{mn}}\left(1-||f||_{k}^{2}\right)^{1/2},$
(25)

where the notation $||\cdot||_{k}$ denotes

$||f||_{k}^{2}=\iint k(X,X^{\prime})f(X)f(X^{\prime})dXdX^{\prime}.$ (26)

An identical result is obtained by using (24). Thus (25) represents a
statistical noise level for $\gamma_{k}$ that is associated with the use of a
finite number of particles, and when $m=n$, we see that $\gamma_{k}^{\rm
noise}\propto 1/\sqrt{n}$. The probability $P$ that $\gamma_{k}$ exceeds a
threshold value $\tau>0$ is then bounded above by Markov’s inequality, which
states that for any $\tau>0$:

$P(\gamma_{k}>\tau)\leq\left(\frac{\gamma_{k}^{\rm noise}}{\tau}\right)^{2}.$
(27)

It follows that (given the null hypothesis) large deviations above the noise
floor $\gamma_{k}^{\rm noise}$ are rare. A more detailed investigation [7]
reveals that the bound in (27) is loose–that is, the probability on the left-
hand side in (27) is much smaller and decays more quickly with $\tau$ than
(27) alone would suggest.

###  III.4 Numerical implementation

The expressions in (21) and (24) are straightforward to implement in a
parallel particle-based beam dynamics simulation code. The resulting numerical
diagnostic, which we denote by MMD, may be used to compare an evolving
particle distribution $f_{t}$ with itself after successive $t$-intervals
$\Delta t$ (e.g., lattice periods), or to compare the evolving particle
distribution against a fixed reference distribution (usually an initial or
predicted final distribution). The frequency samples
$\\{\omega_{l}\\}_{l=1}^{L}$ in (24) may be generated once and stored at
initialization of the simulation, or the samples may be drawn independently at
each evaluation of the MMD. The latter is the approach favored in the
literature.

An algorithm can also be implemented to estimate $\gamma_{k}$ using the
representation given in (18). When the basis functions
$\\{e_{l}\\}_{l=1}^{\infty}$ are known, it is straightforward to estimate the
moments (17) from particle samples. However, the number of basis functions
required to obtain convergence of the sum (18) grows rapidly with the phase
space dimension $d$. For a Gaussian kernel with $d\geq 2$, we found that this
algorithm was outperformed by the spectral algorithm (24) for all examples
tested.

A remark about the choice of kernel $k$ is in order. For simplicity, we use a
Gaussian kernel (67) along each unbounded phase space coordinate. The kernel
width $\sigma$ is chosen to coincide with a typical rms beam size. In some
cases, it is natural to consider periodic domains (for example, if one models
a longitudinal beam slice with periodic boundary conditions). Along a phase
space coordinate that is naturally periodic, we use a Poisson kernel (71) with
the appropriate periodicity. For dynamical problems, it is important that the
kernel remain fixed throughout the simulation.

##  IV Hilbert Schmidt Correlation

Recall that two random variables $X$ and $Y$ described by a joint probability
density $P_{XY}$ are said to be independent when $P_{XY}=P_{X}P_{Y}$, where
$P_{X}$ and $P_{Y}$ are the marginal densities given by:

$\displaystyle P_{X}(X)$ $\displaystyle=\int P_{XY}(X,Y)dY,$ (28a)
$\displaystyle P_{Y}(Y)$ $\displaystyle=\int P_{XY}(X,Y)dX.$ (28b)

Given a metric $\rho$ on the set of such joint probability densities, a
natural measure of deviation from independence between $X$ and $Y$ is given by
the distance $\rho(P_{XY},P_{X}P_{Y})$. This motivates the following
definitions.

###  IV.1 Definition and properties

For simplicity, we assume that $X$ and $Y$ take their values in the same space
$\mathbb{R}^{d}$, on which a kernel $k$ is defined. We define a new kernel
$\kappa$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ by:

$\kappa((X,Y),(X^{\prime},Y^{\prime}))=k(X,X^{\prime})k(Y,Y^{\prime}).$ (29)

It follows that $\kappa$ is symmetric and positive definite (9).

The Hilbert Schmidt correlation $\mathcal{R}_{k}$ between $X$ and $Y$ is then
defined by:

$\mathcal{R}_{k}^{2}(X,Y)=\frac{\gamma_{\kappa}^{2}(P_{XY},P_{X}P_{Y})}{\gamma_{\kappa}(P_{XX},P_{X}P_{X})\gamma_{\kappa}(P_{YY},P_{Y}P_{Y})}.$
(30)

The quantity in the numerator is known as the Hilbert Schmidt Independence
Criterion (HSIC) [10, 37, 38]. The normalizing factor in the denominator
appears in [33], and is designed to ensure that $\mathcal{R}_{k}\leq 1$. The
joint densities $P_{XX}$ and $P_{YY}$ represent the limiting case of perfect
correlation when $X=Y$, namely:

$\displaystyle P_{XX}(X,Y)$ $\displaystyle=P_{X}(X)\delta(Y-X),$ (31)
$\displaystyle P_{YY}(X,Y)$ $\displaystyle=P_{Y}(Y)\delta(X-Y).$

It may be shown that the quantity in (30) satisfies:

$\displaystyle 0\leq\mathcal{R}_{k}\leq 1,$ (32)
$\displaystyle\mathcal{R}_{k}=0\text{ if and only if }X\text{ and }Y\text{ are
independent},$ $\displaystyle\mathcal{R}_{k}=1\text{ if }X\text{ and }Y\text{
are identical}.$

It follows that $\mathcal{R}_{k}$ provides a natural measure of (possibly
nonlinear) correlation between $X$ and $Y$. The special case when
$k(X,X^{\prime})=|X|+|X^{\prime}|-|X-X^{\prime}|$ is known in the statistics
literature as the distance correlation (dCor) [33, 39, 34, 40].

When the kernel $k$ is translation-invariant, we may use the representation
(15) to write the numerator of (30) as:

$\displaystyle\gamma_{\kappa}^{2}(P_{XY},P_{X}P_{Y})=$
$\displaystyle\int\left|\phi_{XY}(\omega,\omega^{\prime})-\phi_{X}(\omega)\phi_{Y}(\omega^{\prime})\right|^{2}\Lambda(\omega)\Lambda(\omega^{\prime})d\omega
d\omega^{\prime},$ (33)

where $\Lambda$ is the spectral density of $k$ defined in (14), and

$\displaystyle\phi_{XY}(\omega,\omega^{\prime})=\int e^{i\left(\omega\cdot
X+\omega^{\prime}\cdot Y\right)}P_{XY}(X,Y)dXdY,$ (34a)
$\displaystyle\phi_{X}(\omega)=\phi_{XY}(\omega,0),\quad\quad\phi_{Y}(\omega^{\prime})=\phi_{XY}(0,\omega^{\prime}).$
(34b)

Corresponding expressions for the factors in the denominator of (30) are
obtained from (33) by taking $Y\mapsto X$ or $X\mapsto Y$ as appropriate, and
using (31) to write:

$\displaystyle\phi_{XX}(\omega,\omega^{\prime})$ $\displaystyle=\int
e^{i\left(\omega\cdot X+\omega^{\prime}\cdot Y\right)}P_{XX}(X,Y)dXdY$
$\displaystyle=\phi_{X}(\omega+\omega^{\prime}),$ (35a)
$\displaystyle\phi_{YY}(\omega,\omega^{\prime})$ $\displaystyle=\int
e^{i\left(\omega\cdot X+\omega^{\prime}\cdot Y\right)}P_{YY}(X,Y)dXdY$
$\displaystyle=\phi_{Y}(\omega+\omega^{\prime}).$ (35b)

###  IV.2 Sample estimate

An estimate of the numerator of (30) from particle (sample) data is given by:

$\displaystyle\gamma_{\kappa}^{2}(P_{XY},P_{X}P_{Y})$
$\displaystyle=\frac{1}{m^{2}}\sum_{i,j=1}^{m}k(X_{i},X_{j})k(Y_{i},Y_{j})$
$\displaystyle+\frac{1}{m^{4}}\sum_{i,j,q,r=1}^{m}k(X_{i},X_{j})k(Y_{q},Y_{r})$
$\displaystyle-\frac{2}{m^{3}}\sum_{i,j,q=1}^{m}k(X_{i},X_{j})k(Y_{i},Y_{q})$
(36)

where the $m$ pairs $\\{(X_{j},Y_{j})\\}_{j=1}^{m}$ are sampled from $P_{XY}$.
The corresponding expressions appearing in the denominator of (30) are
obtained from (36) by taking $Y_{j}\mapsto X_{j}$ and $X_{j}\mapsto Y_{j}$,
respectively. It can be shown that (36) can be computed with $O(m^{2})$
complexity by introducing $O(m)$ storage.

In the special case that the kernel $k$ is translation invariant, this
complexity can be further reduced by working in frequency space using (33). An
efficient estimate is given by:

$\displaystyle\gamma_{\kappa}^{2}$
$\displaystyle(P_{XY},P_{X}P_{Y})=\frac{1}{L}\sum_{k=1}^{L}\left|\frac{1}{m}\sum_{j=1}^{m}e^{i\omega_{2k-1}\cdot
X_{j}}e^{-i\omega_{2k}\cdot Y_{j}}\right.$ (37)
$\displaystyle\left.-\left(\frac{1}{m}\sum_{j=1}^{m}e^{i\omega_{2k-1}\cdot
X_{j}}\right)\left(\frac{1}{m}\sum_{j=1}^{m}e^{-i\omega_{2k}\cdot
Y_{j}}\right)\right|^{2},$

where $\\{\omega_{k}\\}_{k=1}^{2L}$ denote $2L$ frequency vectors that are
randomly sampled from the spectral probability density $\Lambda$ associated
with $k$. Note that the complexity of (37) is given by $O(mL)$. (This is a
variant of the random Fourier features estimate appearing in [41, 36].)

###  IV.3 Statistical error

In practice, the sample estimate in (36) or (37) is used to test the
hypothesis that two random variables $X$ and $Y$ are independent. To set a
confidence threshold for this test, one needs to know the statistical
distribution of (36) under the null hypothesis that $P_{XY}=P_{X}P_{Y}$. This
problem has been treated in detail [10, 37, 38, 33]. For our purposes, it is
enough to know the rms value of $\mathcal{R}_{k}$ under the hypothesis that
$X$ and $Y$ are independent. This is given by taking the expected value of
(36), which yields:

$\displaystyle\operatorname{E}[\gamma_{k}^{2}(P_{XY},P_{X}P_{Y})]$
$\displaystyle\quad\quad=\frac{(m-1)}{m^{2}}\left(1-||P_{X}||_{k}^{2}\right)\left(1-||P_{Y}||_{k}^{2}\right),$
(38)

where $||\cdot||_{k}$ has the same meaning as in (26). Thus we have, to
leading order in $1/m$:

$\mathcal{R}_{k}^{\rm
noise}=\frac{1}{\sqrt{m}}\frac{\left(1-||P_{X}||_{k}^{2}\right)^{1/2}\left(1-||P_{Y}||_{k}^{2}\right)^{1/2}}{\gamma_{\kappa}(P_{XX},P_{X}P_{X})\gamma_{\kappa}(P_{YY},P_{Y}P_{Y})}.$
(39)

Note that (39) is fully determined by the marginal densities $P_{X}$ and
$P_{Y}$ through (31).

Given the null hypothesis, an inequality corresponding to (27) holds after
$\gamma_{k}$ is replaced by $\mathcal{R}_{k}$, indicating that large
deviations above the noise floor value $\mathcal{R}_{k}^{\rm noise}$ are rare.

###  IV.4 Numerical Implementation

In practice, the random variables $X$ and $Y$ described above may represent
two phase space coordinates within a single beam (e.g., $X=z$ and $Y=p_{z}$)
or two $d$-tuples of phase space coordinates (e.g., $X=(x,p_{x})$ and
$Y=(y,p_{y})$). In this case, computation of $\mathcal{R}_{k}(X,Y)$ using (36)
or (37) returns a measure of correlation between phase space coordinates
within the beam.

Alternatively, let $X$ denote the vector of phase space coordinates for a
particle within the beam at initial time (or lattice location) $t=0$, and let
$Y$ denote the vector of phase space coordinates for the same particle at a
later time $t$. Then $\mathcal{R}_{k}(X,Y)$ measures the correlation of a
particle’s coordinates at time $t$ with the particle’s initial coordinates,
and this quantity will be denoted $\mathcal{R}_{k}(t)$ for simplicity. The
dynamical evolution of this quantity is intimately related to mixing. (See
Appendix B.)

Numerical evaluation of $\mathcal{R}_{k}(t)$ requires that each particle be
assigned a unique index $j$, so that one may construct the particle pairs
$\\{(X_{j},Y_{j})\\}_{j=1}^{m}$ at each desired evaluation time $t$. In
particular, the arrays containing the initial and final phase space
coordinates of particle $j$ must be stored on the same processor, which
requires appropriate bookkeeping and possible communication.

##  V Applications to Idealized and Exactly-Solvable Models

In this section, we illustrate several applications of the above tools to
dynamical problems involving beams: 1) to compare two beams for benchmarking
and quantifying beam mismatch, 2) to detect nonlinear phase space correlations
and coupling between phase planes, 3) to verify numerically that a beam that
is generated from a well-matched distribution remains stationary, 4) to study
the relaxation of a non-equilibrium beam to a stationary state, and 5) to
measure the rate of chaotic mixing and decay of correlations within a beam
during its evolution.

To aid in benchmarking the diagnostics MMD $(\gamma_{k}$) and HSCor
$(\mathcal{R}_{k})$, idealized distributions and exactly-solvable models
involving maps are used. The next section will discuss realistic examples
involving high-intensity beams.

###  V.1 Distribution comparison and mismatch

Let $f$ and $f^{\prime}$ denote two centered Gaussian distributions with
covariance matrices $\Sigma$ and $\Sigma^{\prime}$, respectively. To determine
the distance between these distributions, we consider an arbitrary Gaussian
kernel $k$ of the form:

$k(X,X^{\prime})=\exp\left(-\frac{1}{2}(X-X^{\prime})^{T}S(X-X^{\prime})\right),$
(40)

where $S$ is any symmetric, positive definite matrix. Such a matrix may always
be decomposed as $S=A^{T}A$ for some matrix $A$. The MMD between $f$ and
$f^{\prime}$ may be obtained using (11) as:

$\displaystyle\gamma_{k}^{2}(f,f^{\prime})$
$\displaystyle=\det\left(I+2\Sigma_{N}\right)^{-1/2}+\det(I+2\Sigma_{N}^{\prime})^{-1/2}$
$\displaystyle-2\det(I+\Sigma_{N}+\Sigma_{N}^{\prime})^{-1/2},$ (41)

where $I$ is the identity matrix, and the normalized covariance matrices are:

$\Sigma_{N}=A\Sigma
A^{T},\quad\quad\Sigma_{N}^{\prime}=A\Sigma^{\prime}A^{T}.$ (42)

This yields a large class of examples for benchmarking the computation of MMD
in any dimension.

As a special case of (41), consider two Gaussian beams $f$ and $f^{\prime}$
described on a 2D phase space $(x,p_{x})$ with identical emittance $\epsilon$,
with covariance matrices:

$\Sigma=\epsilon\begin{pmatrix}\beta&-\alpha\\\
-\alpha&\gamma\end{pmatrix},\quad\quad\Sigma^{\prime}=\epsilon\begin{pmatrix}\beta^{\prime}&-\alpha^{\prime}\\\
-\alpha^{\prime}&\gamma^{\prime}\end{pmatrix}.$ (43)

Taking $S=\Sigma^{-1}$ in (40) and computing the MMD between $f$ and
$f^{\prime}$ using (41) gives:

$\gamma_{k}^{2}(f,f^{\prime})=\frac{1}{3}-\frac{1}{\sqrt{5+4\zeta}}.$ (44)

Here the result is expressed in terms of the linear beam mismatch parameter
$\zeta$, given by:

$\zeta=\frac{1}{2}\left(\beta\gamma^{\prime}-2\alpha\alpha^{\prime}+\gamma\beta^{\prime}\right),\quad\zeta\geq
1.$ (45)

The same result is obtained by taking $S=(\Sigma^{\prime})^{-1}$ in (40). Note
that (44) vanishes when $\zeta=1$, and increases monotonically with increasing
mismatch $\zeta$. Thus, when the MMD can be expressed in terms of the linear
mismatch, the result behaves as expected.

In addition to detecting differences based on the second beam moments, the MMD
detects differences in the details of two distributions. For a 4D example
relevant to beams, consider a K-V distribution $f_{KV}$ and a (4D) Gaussian
distribution $f_{G}$ with the same second moments. Using (15) gives the exact
result:

$\gamma_{k}(f_{KV},f_{G})=\sqrt{\frac{1}{9}-\frac{1}{2e}+\frac{I_{2}(4)+I_{3}(4)}{2e^{4}}},$
(46)

where $I_{n}$ denotes the modified Bessel function of order $n$. This
corresponds to the numerical value $\gamma_{k}(f_{KV},f_{G})\approx 0.12863$.
For comparison, a numerical computation of $\gamma_{k}$ using (24) from two
sampled beams with $m=n=10^{5}$ particles and $L=10^{4}$ frequency samples
gives an estimated value of 0.128.

###  V.2 Detecting phase space correlations

For a Gaussian distribution of any dimension, one may detect linear
correlations among the various degrees of freedom by using (41) in (30). As an
example, consider a Gaussian distribution on a 2D phase space $(q,p)$ with the
covariance matrix:

$\Sigma=\begin{pmatrix}1&r\\\ r&1\\\ \end{pmatrix},\quad f(q,p)=\frac{1}{2\pi
a}e^{-(q^{2}+p^{2}-2rqp)/2a^{2}},$ (47)

where $a=\sqrt{1-r^{2}}$. Using a Gaussian kernel of unit width, the
correlation $\mathcal{R}_{k}$ between the variables $q$ and $p$ is given by:

$\mathcal{R}_{k}^{2}=\frac{g(r)}{g(1)},\quad
g(r)=\frac{1}{3}+\frac{1}{\sqrt{9-4r^{2}}}-\frac{2}{\sqrt{9-r^{2}}}.$ (48)

This result is expressed in terms of the standard linear correlation
coefficient $r$. Note that $\mathcal{R}_{k}$ increases monotonically from 0 to
1 as $|r|$ increases from 0 to 1.

The quantity $\mathcal{R}_{k}$ is useful for detecting nonlinear correlations,
even when the exact structure of the correlation is unknown. Consider the case
of a Gaussian beam with a quadratic correlation in the longitudinal phase
space $(z,\delta)$:

$f(z,\delta)=\frac{1}{2\pi\sigma\tau}e^{-z^{2}/{2\sigma^{2}}}e^{-(\delta+az^{2})^{2}/2\tau^{2}}.$
(49)

The linear correlation between $z$ and $\delta$ in (49) vanishes, since one
may verify that $\langle{z\delta\rangle}=0$. Using the dimensionless variables
$\bar{z}=z/\sigma$ and $\bar{\delta}=\delta/(a\sigma^{2})$, and choosing a
Gaussian kernel of width 1, one may evaluate the correlation $\mathcal{R}_{k}$
between $\bar{z}$ and $\bar{\delta}$. The result is shown in Fig. 1 as a
function of the dimensionless parameter $\bar{\tau}=\tau/(a\sigma^{2})$. The
result is independent of the sign of $a$. We see that $\mathcal{R}_{k}$
becomes small as the quadratic coefficient $a$ becomes small or as the
uncorrrelated energy spread $\tau$ becomes large, as one might expect.

Figure 1: Hilbert-Schmidt correlation $\mathcal{R}_{k}$ between variables $z$
and $\delta$ in the longitudinal phase space of a bunch with a quadratic
correlation (49). (Inset) Sampled particles ($10^{4}$) for the case
$\tau/a\sigma^{2}=1$, yielding the computed value $\mathcal{R}_{k}\approx
0.375$.

###  V.3 Testing for stationarity

Often one must characterize the degree to which a given particle distribution
remains stationary over many periods of dynamical evolution. This can be
quantified by computing the distance $\gamma_{k}(f_{t},f_{0})$ between the
initial distribution $f_{0}$ and the distribution $f_{t}$ after $t$ periods.

For example, consider the 2D nonlinear symplectic map given by [42]:

$\begin{pmatrix}q^{f}\\\
p^{f}\end{pmatrix}=\begin{pmatrix}\cos\phi&\sin\phi\\\
-\sin\phi&\cos\phi\end{pmatrix}\begin{pmatrix}q\\\
p\end{pmatrix},\quad\phi=\psi+\frac{\alpha}{2}\left(q^{2}+p^{2}\right),$ (50)

where $\psi>0$ and $\alpha>0$ are constants. This may be viewed as a simple
model of a betatron phase advance in a single plane that increases linearly
with the action $J=(q^{2}+p^{2})/2$.

A Gaussian distribution of the form:

$f(q,p)=\frac{1}{2\pi\epsilon_{0}}e^{-(q^{2}+p^{2})/2\epsilon_{0}}$ (51)

is an explicit function of the action $J$, and is therefore invariant under
the map (50). Sampling $n=10^{4}$ particles from (51) and tracking them under
iterates of the map (50), we compare the distribution at each iteration with
the initial distribution. The result is shown in Fig. 2.

Figure 2: Dynamics of a 2D beam with $n=10^{4}$ particles sampled from a
matched Gaussian distribution (51) evolving under iteration of the map (50).
The quantity $\gamma_{k}(f_{t},f_{0})$ is shown as a function of the iteration
number $t$ for 4 distinct random seeds, showing that the distribution remains
stationary.

The MMD distance $\gamma_{k}(f_{t},f_{0})$ to the initial distribution is
nonzero after the first iteration, but remains near $10^{-2}$ over the time
interval observed. Due to the finite number of particles, the value of
$\gamma_{k}(f_{t},f_{0})$ experiences statistical fluctuations around the
predicted rms value given by (25) (black curve).

###  V.4 Relaxation to a stationary state

If the distribution (51) is given an initial centroid offset $q\mapsto
q+q_{0}$ with $q_{0}>0$, the beam will filament and converge weakly to a
stationary equilibrium of the form [43, 42]:

$f_{eq}(q,p)=\frac{1}{2\pi\epsilon_{0}}e^{-(q^{2}+p^{2}+q_{0}^{2})/2\epsilon_{0}}I_{0}\left(\frac{q_{0}}{\epsilon_{0}}\sqrt{q^{2}+p^{2}}\right).$
(52)

Using a Gaussian kernel of width $\sigma$ in expression (15), one may solve
exactly for the time evolution of the MMD distance to equilibrium:

$\gamma_{k}^{2}(f_{t},f_{eq})=2s^{2}\sum_{n=1}^{\infty}\nu_{n}e^{-\nu_{n}(1+2d^{2}s^{2}\tau_{n}^{2})}I_{n}(\nu_{n}),$
(53)

where $I_{n}$ denotes the modified Bessel function of order $n$ and

$\nu_{n}=\frac{d^{2}}{1+d^{2}s^{2}(1+\tau_{n}^{2})},$ (54)

is given in terms of the dimensionless parameters:

$\tau_{n}=nt\alpha\epsilon_{0},\quad
d^{2}=\frac{q_{0}^{2}}{2\epsilon_{0}},\quad s=\frac{\sigma}{q_{0}}.$ (55)

In particular, we see that for large $t$, (53) converges to zero as
$\gamma_{k}\sim O(1/t^{2})$.

Fig. 3 illustrates the result obtained from tracking $10^{5}$ particles for
the case $\epsilon_{0}=0.01$, $q_{0}=0.5$, $\psi=0.3$, $\alpha=0.1$. We use a
Gaussian kernel with $\sigma=1$. The inset shows the difference from the
prediction. By $t=500$, the distribution has converged to equilibrium to
within the resolution set by the particle noise.

Figure 3: Evolution of a distribution with $n=10^{5}$ particles sampled from
(51) with initial offset $q\mapsto q+q_{0}$ under iteration of the map (50).
The distance of the distribution to the predicted equilibrium (52) after $t$
iterates is shown. (Red curve) Analytical prediction (53). (Black, inset)
Absolute error in the numerical result obtained using (21). (Blue, inset)
Absolute error in the numerical result obtained using (24). (Red, dashed)
Prediction (25) of numerical noise using the distribution $f_{eq}$.

This shows that MMD provides a diagnostic capable of measuring the dynamical
relaxation of a beam to a stationary state.

###  V.5 Mixing and decay of correlations

Although the term “mixing” is sometimes used to refer to any process involving
filamentation and relaxation of the beam in phase space, we distinguish
between regular mixing (which is characteristic of nonlinear integrable
systems) and mixing in the ergodic sense (which is characteristic of systems
with widespread chaos). Here, we refer to the latter meaning of the term, as
it is formalized in ergodic theory [44]. (See Appendix B.)

A simple illustration of chaotic mixing behavior is given by the Arnold cat
map, which is the 2D symplectic map given by:

$\begin{pmatrix}q^{f}\\\ p^{f}\end{pmatrix}=\begin{pmatrix}2&1\\\
1&1\end{pmatrix}\begin{pmatrix}q\\\ p\end{pmatrix}\mod 2\pi,$ (56)

where we assume that $q$ and $p$ each have period $2\pi$. The uniform density:

$f(q,p)=\frac{1}{(2\pi)^{2}},\quad q,p\in[0,2\pi)$ (57)

is invariant under (56) since the map is area-preserving.

The Hilbert Schmidt correlation provides a quantitative measure of mixing and
the resulting decay of correlations over time. To illustrate this, we compute
the value $\mathcal{R}_{k}(X(t),X(0))$, where $X(0)$ is the random variable
denoting a particle’s initial phase space coordinates, as sampled from (57),
and $X(t)$ denotes the particle’s phase space coordinates after $t$ iterations
of the map (56). For simplicity, we denote this quantity by
$\mathcal{R}_{k}(t)$. In a periodic domain, it is appropriate to use a kernel
that reflects the underlying domain periodicity. Using a Poisson kernel with
parameter $\sigma$ in each dimension (Appendix A), we may compute the value of
$\mathcal{R}_{k}(t)$ explicitly. The value after $t$ iterations is given
exactly by the sum:

$\mathcal{R}_{k}(t)=\left(\frac{1-\sigma^{2}}{2\sigma}\right)\left[\sum_{\begin{subarray}{c}n_{1},n_{2}=-\infty\\\
n_{1}\neq 0,n_{2}\neq
0\end{subarray}}^{\infty}\sigma^{C(n_{1},n_{2},t)}\right]^{1/2},$ (58)

where the exponent is $(t=1,2,3,\ldots)$:

$\displaystyle C(n_{1},$ $\displaystyle n_{2},t)=|n_{1}|+|n_{2}|$
$\displaystyle+\left|F_{2t+1}n_{1}+F_{2t}n_{2}\right|+|F_{2t}n_{1}+F_{2t-1}n_{2}|$
(59)

given in terms of the usual Fibonacci sequence:

$F_{1}=F_{2}=1,\quad F_{t}=F_{t-1}+F_{t-2}.$ (60)

Fig. 4 shows the decay of the quantity $\mathcal{R}_{k}(t)$ as a function of
the iteration number $t$ for the case $\sigma=1/2$. Because the prediction
(58) is only defined for nonnegative integer values of $t$, the red curve
shown connects these values using smooth interpolation. Note that mixing for
the map (56) is very rapid. After only 4-5 iterations, the correlations
between the initial and final phase space coordinates are at or below the
level expected due to numerical particle noise, as given by (39) and indicated
by the black line. The inset shows that no correlation is visible by eye
between the initial and final coordinates.

Figure 4: Dynamics of a beam with $n=10^{5}$ particles sampled from the
density (57) evolving under iteration of the map (56). The quantity
$\mathcal{R}_{k}(t)$ is shown as a function of the iteration number $t$,
illustrating the decay of correlations due to mixing. The red curve shows the
prediction (58), the black points are the results of simulation, and the black
curve denotes the expected rms value due to noise (39). (Inset) Plot of
initial $q$ vs. final $q$ after 10 iterations, showing no visible correlation.

##  VI Application to High-Intensity Beams

###  VI.1 Beam in a constant focusing channel

As our first example with space charge, we consider a long (unbunched) intense
beam in a constant focusing channel that is initialized in a stationary
thermal equilibrium, so the 4D beam distribution takes the form:

$f_{0}(x,p_{x},y,p_{y})\propto e^{-H(x,p_{x},y,p_{y})/kT},$ (61)

where $H$ denotes the self-consistent Hamiltonian:

$H=\frac{1}{2}(p_{x}^{2}+p_{y}^{2})+\frac{1}{2}\Omega^{2}(x^{2}+y^{2})+\frac{q\phi(x,y,s)}{\beta^{2}\gamma^{3}mc^{2}},$
(62)

and $\phi$ is a solution of the 2D Poisson equation:

$\nabla_{\perp}^{2}\phi=-\frac{\lambda}{\epsilon_{0}}\int
f_{0}(x,p_{x},y,p_{y})dp_{x}dp_{y}.$ (63)

Here $\lambda$ denotes the charge per unit length, and $\epsilon_{0}$ denotes
the vacuum permittivity.

For simulation, we consider a proton beam with a kinetic energy of 200 MeV and
a beam current of 20 A in an external focusing of strength $\Omega=0.628$
$m^{-1}$ (corresponding to a 2.7 T solenoid field). The temperature $kT$ is
chosen to yield the initial emittance $\epsilon_{x,rms}=\epsilon_{y,rms}=1.24$
$\mu$m. The tune depression due to space charge is then given by
$\nu/\nu_{0}\approx 0.55$.

Fig. 5 illustrates the MMD distance between the initial distribution and the
distribution at time $t$, for four distinct random seeds. To compute (24), a
Gaussian kernel was used. The kernel width parameter $\sigma$ associated with
each phase space dimension was matched to the corresponding rms width of the
distribution (61). For each random seed, self-consistent tracking with
transverse (2D) space charge using $n=10^{6}$ particles was performed using
the symplectic gridless spectral solver described in [45]. Notice that the
distribution remains stationary to within the level expected due to particle
noise (black line). Compare this dynamical behavior to that shown in Fig. 2.

Figure 5: Test of stationarity for an unbunched (4D) beam with $n=10^{6}$
particles sampled from a thermal equilibrium distribution (61) propagating in
a linear constant focusing channel. The quantity $\gamma_{k}(f_{t},f_{0})$ is
shown as a function of propagation distance $t$ for 4 distinct random seeds.

As a second example, consider a proton beam with the same energy and emittance
as above. However, instead of the stationary distribution (61), we use the
initial distribution [6]:

$f_{0}(x,p_{x},y,p_{y})=\frac{1}{\pi^{2}r_{m}^{2}p_{m}^{2}}\Theta(r_{m}-r)\Theta(p_{m}-p)$
(64)

where $r=\sqrt{x^{2}+y^{2}}$, $p=\sqrt{p_{x}^{2}+p_{y}^{2}}$, and $\Theta$
denotes the unit step function. The distribution (64) is not stationary, but
to minimize fluctuations of the rms beam size, we match the beam in an rms
sense by setting $p_{m}^{2}=\Omega^{2}r_{m}^{2}-K_{pv}$, where $K_{pv}$
denotes the generalized beam perveance. In the absence of precise knowledge of
the final equilibrium state, we compare the particle distribution at
successive time intervals separated by $\Delta t=L$, where $L=6.20$ m denotes
the period of linearized envelope oscillations about the equilibrium beam
size. The result is shown in Fig. 6. On the time scale shown, the beam appears
to undergo relaxation toward a final distribution containing a small but
visible low-density halo (inset). This fast relaxation appears to be a
property of the collisionless Vlasov-Poisson system, and it is to be
distinguished from slow relaxation due to collisional effects [4], which are
not included in the numerical model.

Figure 6: Collisionless relaxation of an unbunched (4D) beam with $n=10^{6}$
particles sampled from the distribution (64) in a linear constant focusing
channel (62), where $t$ denotes the number of linearized envelope periods. The
MMD between the distribution on successive periods $\gamma_{k}(f_{t},f_{t-1})$
decays to the level of noise as the beam relaxes to a stationary state.
(Inset) The final particle distribution (one quadrant).

The system described in this section cannot be mixing due to the existence of
invariants of motion (for example, the angular momentum). Indeed, the computed
quantity $\mathcal{R}_{k}(t)$ converges rapidly to a nonzero value. This
behavior will be examined using a more complex example in the following
section.

###  VI.2 Beam in a periodic FODO lattice

As a simple application to an intense beam in a periodic focusing structure,
we consider a lattice consisting of a single FODO cell (Fig. 7).

Figure 7: Matched (K-V) beam envelopes for a 10 A proton beam at 200 MeV in
the FODO cell used in Section VI.b. Red rectangle - focusing quad. Blue
rectangle - defocusing quad.

We consider a proton beam with a kinetic energy of 200 MeV and a beam current
of 10 A with an initial emittance of $\epsilon_{x,rms}=\epsilon_{y,rms}=1$
$\mu$m (unnormalized). The initial distribution is Gaussian of the form:

$f_{0}(x,p_{x},y,p_{y})=\frac{1}{(2\pi)^{2}\epsilon_{x,rms}\epsilon_{y,rms}}e^{-\frac{1}{2}X^{T}\Sigma^{-1}X},$
(65)

where $X=(x,p_{x},y,p_{y})$ and $\Sigma$ denotes the covariance matrix.
Although the second beam moments are chosen so that the beam is matched in an
rms-sense, the distribution is not matched in detail. Fig. 7 shows the
evolution of the matched beam envelopes over a single period. The zero-current
phase advance is $60.1^{\circ}$ per cell, while the 10 A phase advance is
$25.9^{\circ}$ per cell, so that space charge plays a significant role.

Fig. 8 illustrates the MMD distance between the initial distribution and the
distribution after $t$ periods (blue), together with the MMD distance between
successive periods (red). To compute (24), a Gaussian kernel was used whose
width along each dimension matches the initial rms beam size in that
dimension. In each case, self-consistent tracking with transverse (2D) space
charge using $n=10^{6}$ particles was performed using the algorithm described
in [45]. The distribution relaxes quickly over the first 100 periods, but
fluctuations above the noise level persist on a much longer time scale.

Figure 8: Dynamics of an unbunched (4D) beam with $n=10^{6}$ particles sampled
from the distribution (65) in the FODO channel shown in Fig. 7. The MMD
between the distribution on successive lattice periods
$\gamma_{k}(f_{t},f_{t-1})$ decays to near (but remains slightly above) the
level of noise (black line). The MMD to the initial distribution
$\gamma_{k}(f_{t},f_{0})$ is largely unchanged after the first 100 periods.

Fig. 9 shows the correlation $\mathcal{R}_{k}(t)$ between the distribution at
time $t$ and the distribution at $t=0$. Note that $\mathcal{R}_{k}(t)$
converges to a fixed nonzero value within just a few periods, and then remains
constant. This indicates that particle coordinates remain correlated with
their initial values indefinitely, and that the dynamics is not mixing. This
generally suggests the existence of invariants of motion in some regions of
the phase space. A 2D plot of $y(0)$ versus $y(1000)$ is also shown. Note that
the correlations are not easily visible. In fact, if $\mathcal{R}_{k}$ is
computed using only the initial and final values of $y$, neglecting all other
coordinates, then the corresponding value is 0.08. This shows that
$\mathcal{R}_{k}$ can quantify correlations present in higher dimensions that
are not easily visualized by viewing 2D projections.

Figure 9: Dynamics of an unbunched (4D) beam with $n=10^{6}$ particles sampled
from the distribution (65) in the FODO channel shown in Fig. 7. The
correlation $\mathcal{R}_{k}(t)$ between the distribution after $t$ periods
and the initial distribution at $t=0$ is shown. Correlations appear to persist
indefinitely. (Inset) Plot of initial $y$ vs. final $y$ after 1K periods,
showing that the correlations are not easily visible in low-dimensional
projections.

###  VI.3 Treatment of beam loss

In the presence of beam loss, it may be of interest to study the dynamics of
the beam on a bounded subregion $E$ of the phase space (e.g., defined by the
vacuum chamber or by the dynamic aperture). For example, one may study the
relaxation of the distribution defined by those particles that remain
indefinitely within the region $E$. In this case, the beam distribution
function $f_{t}$ satisfies:

$\int_{E}f_{t}(X)dX=\chi_{t},\quad 0\leq\chi_{t}\leq 1,$ (66)

where $\chi_{t}$ denotes the (possibly time-dependent) fraction of beam
particles that lie within the region $E$.

The formalism of the previous section may be modified to treat this case as
follows.

  * •

In the calculation of (20-24) and (36-37), sum only over those particles that
lie within the desired region $E$. This is equivalent to modifying the kernel
$k$ by setting $k(X,X^{\prime})=0$ if $X\notin E$ or $X^{\prime}\notin E$.

  * •

In the calculation of (20-24) and (36-37), replace the weight coefficients
$1/m$ by particle weights $w_{j}$ with $\sum_{j=1}^{m}w_{j}=\chi_{t}$, and
similarly for $1/n$.

One may verify that many of the desired mathematical properties of the MMD and
HSICor still hold for these modified statistics.

##  VII Conclusions

Modeling beams in the presence of high intensity space charge, collective
instabilities, or strong nonlinear focusing can result in dynamical processes
that are difficult to characterize using typical numerical diagnostics based
on second-order moments. We have introduced two numerical diagnostics
originating in the ML literature with highly desirable mathematical properties
that are straightforward to implement in parallel particle-based simulation
codes. The first is a measure of statistical distance known as the Maximum
Mean Discrepancy ($\gamma_{k}$), which serves as a measure of similarity
between two particle ensembles. The second is a measure of statistical
dependence or correlation between random variables known as the Hilbert
Schmidt Correlation ($\mathcal{R}_{k}$). These quantities are useful for a
variety of applications, including: matching a beam into a periodic transport
system, numerical benchmarking, detecting possibly nonlinear phase space
correlations, characterizing relaxation to a (quasi-)stationary state, and
characterizing mixing or decay of correlations within the beam.

It is important to note that the quantitative results obtained will depend on
the choice of the kernel $k$. On one hand, this kernel-dependence may be
viewed as a disadvantage of the diagnostics described here. On the other hand,
one may view the kernel as a natural way to parameterize a large family of
possible diagnostics, all of which correctly capture the same underlying
physical processes. (This is a consequence of the mathematical properties
described in Sections III-IV.) In numerical experiments the observed dynamical
evolution of $\gamma_{k}$ or $\mathcal{R}_{k}$ was largely independent of the
choice of kernel, although this remains an active area of investigation. An
alternative and parameter-free statistical distance with similar mathematical
properties is the Wasserstein distance $W_{p}$ (Section II.B). In the future,
the authors plan to investigate the feasibility of using efficient
approximations to $W_{p}$ [46] for beam dynamics applications.

The diagnostics described here are well-suited to applications involving large
simulation ensembles. For example, quantities involving $\gamma_{k}$ or
$\mathcal{R}_{k}$ may be used as objectives for accelerator design
optimization or for training machine learning models that require detailed
information about the beam distribution function. This raises the possibility
of tailoring the final beam phase space density by using large-scale automated
machine tuning.

Finally, although we have focused on the case of charged particle beams, it is
clear that these tools can be applied without change to kinetic simulations of
other many-body systems such as plasmas or gravitational systems, which rely
on the tracking of large particle ensembles.

##  VIII Acknowledgments

This work was supported by the Director, Office of Science of the U.S.
Department of Energy under Contract No. DE-AC02-05CH11231, and made use of
computer resources at the National Energy Research Scientific Computing
Center. The authors acknowledge support from the U.S. DOE Early Career
Research Program under the Office of High Energy Physics.

## Appendix A: Commonly-used kernels

This Appendix lists several of the kernels $k$ most commonly used for
discrimination testing and independence testing in ML. These kernels are all
translation-invariant and normalized as in (12), so that $k(X,X)=1$ for all
$X$. All the kernels listed here have the property that the quantity
$\gamma_{k}$ in (11) satisfies the metric conditions (i-iv) and captures weak
convergence, as described in Section II.

Gaussian kernel $(\sigma>0)$ of dimension $d$:

$\displaystyle k_{\rm
Gaussian}(X,X^{\prime})=e^{-|X-X^{\prime}|^{2}/2\sigma^{2}},$ (67a)
$\displaystyle\Lambda_{\rm
Gaussian}(\omega)=\frac{e^{-\sigma^{2}|\omega|^{2}/2}}{(2\pi\sigma^{2})^{d/2}}.$
(67b)

Laplace kernel $(\sigma>0)$ of dimension $d$:

$\displaystyle k_{\rm Laplace}(X,X^{\prime})=e^{-|X-X^{\prime}|/\sigma},$
(68a) $\displaystyle\Lambda_{\rm
Laplace}(\omega)=\frac{\sigma^{d}}{\pi^{(d+1)/2}}\frac{\Gamma\left(\frac{d+1}{2}\right)}{(1+\sigma^{2}\omega^{2})^{(d+1)/2}},$
(68b)

where $\Gamma$ denotes the gamma function.

Matérn kernel ($\nu>0$, $\sigma>0$) of dimension $d$:

$\displaystyle k_{\rm
Matern}(X,X^{\prime})=\frac{2^{1-\nu}}{\Gamma(\nu)}\zeta^{\nu}K_{\nu}\left(\zeta\right),$
(69a) $\displaystyle\Lambda_{\rm
Matern}(\omega)=\frac{2^{s}\sigma^{d}\nu^{\nu}\Gamma(s)}{(2\pi)^{d/2}\Gamma(\nu)}\left({2\nu}+\sigma^{2}\left|\omega\right|^{2}\right)^{-s},$
(69b)

where $K_{\nu}$ is the modified Bessel function of order $\nu$, and we
abbreviate:

$\zeta=\frac{\sqrt{2\nu}}{\sigma}|X-X^{\prime}|,\quad\quad s=\nu+\frac{d}{2}.$
(70)

Poisson kernel ($0<\sigma<1$) of dimension 1:

$\displaystyle k_{\rm
Poisson}(X,X^{\prime})=\frac{(1-\sigma)^{2}}{1-2\sigma\cos(X-X^{\prime})+\sigma^{2}},$
(71a) $\displaystyle\Lambda_{\rm
Poisson}(\omega)=\sum_{n=-\infty}^{\infty}\left(\frac{1-\sigma}{1+\sigma}\right)\sigma^{|n|}\delta(\omega-n).$
(71b)

Another class of kernels $k_{\rm Wendland}$ (of dimension $d$) is constructed
by using a polynomial with compact support in the separation distance
$r=|X-X^{\prime}|/\sigma$, where $\sigma>0$. See [47, 48] for a detailed
description of these. An example for $d=1$ is given by:

$\displaystyle k_{\rm Wend}(X,X^{\prime})=(1-r)^{3}(1+3r),\quad\quad r<1,$
(72a) $\displaystyle\Lambda_{\rm
Wend}(\omega)=\frac{24\sigma}{\pi}\left\\{\frac{2+\cos\omega\sigma}{(\omega\sigma)^{4}}-\frac{3\sin\omega\sigma}{(\omega\sigma)^{5}}\right\\}.$
(72b)

When working with a Gaussian kernel, an orthonormal basis for the RKHS is
given by $\\{e_{n}\\}_{n=0}^{\infty}$ where [49]:

$e_{n}(X)=\sqrt{\frac{1}{\sigma^{2n}n!}}e^{-X^{2}/2\sigma^{2}}X^{n}.$ (73)

When working with the Poisson kernel, a (complex) orthonormal basis for the
RKHS is given by $\\{e_{n}\\}_{n=-\infty}^{\infty}$ with:

$e_{n}(X)=\sqrt{\frac{1-\sigma}{1+\sigma}}\sigma^{|n|/2}e^{inX}.$ (74)

A corresponding real basis is easily constructed.

We can construct a kernel $k$ of higher dimension using kernels $k^{(j)}$
$(j=1,\ldots,d)$ of lower dimension. For example, if we write
$X=(X_{1},\ldots,X_{d})$ and
$X^{\prime}=(X_{1}^{\prime},\ldots,X_{d}^{\prime})$, then

$k(X,X^{\prime})=\prod_{j=1}^{d}k^{(j)}(X_{j},X_{j}^{\prime}).$ (75)

If the kernels $k^{(j)}$ are translation-invariant with spectral densities
$\Lambda^{(j)}$, then setting $\omega=(\omega_{1},\ldots,\omega_{d})$ we have:

$\Lambda(\omega)=\prod_{j=1}^{d}\Lambda^{(j)}(\omega_{j}).$ (76)

If the RKHS associated with the kernel $k^{(j)}$ has basis $e_{n}^{(j)}$, then
the RKHS associated with $k$ has basis

$e_{n}(X)=\prod_{j=1}^{d}e^{(j)}_{n_{j}}(X_{j}),$ (77)

where $n=(n_{1},\ldots,n_{d})$ ranges over all possible indices.

## Appendix B: Definition of mixing

Let $\mathcal{M}$ be a map, and let $f$ denote a probability density that is
invariant under the map, in the sense that:

$f(\mathcal{M}^{-1}(X))=f(X).$ (78)

Then the map $\mathcal{M}$ is mixing with respect to the invariant density $f$
if for any two sets $A$ and $B$ [50]:

$\lim_{t\rightarrow\infty}P(\mathcal{M}^{-t}(A)\cap B)=P(A)P(B),$ (79)

where

$P(A)=\int_{A}f(X)dX$ (80)

denotes the probability that a point lies in $A$. Informally, for any sets $A$
and $B$, the sequence of sets $\mathcal{M}^{-t}(A)$ becomes asymptotically
independent of $B$ as $t\rightarrow\infty$ [44].

There are many equivalent formulations of the condition (79). For our
purposes, it is enough to know that a map $\mathcal{M}$ is mixing with respect
to a density $f$ if and only if $\mathcal{R}_{k}(t)\rightarrow 0$ as
$t\rightarrow\infty$, where $\mathcal{R}_{k}(t)$ denotes the correlation
between the particle coordinates sampled from $f$ at $t=0$ and the
corresponding particle coordinates at later time $t$. (See Section IV.) This
holds for all kernels $k$ satisfying the desirable properties described in
Section II.C.

## References

  * [1] H. Kandrup el al, “Chaotic collisionless evolution in galaxies and charged-particle beams,” Ann. N. Y. Acad. Sci. 1045, 12 (2005).
  * [2] C. Bohn, “Chaotic dynamics in charged-particle beams: possible analogs of galactic evolution,” Ann. N. Y. Acad. Sci. 1045, 34 (2005).
  * [3] D. Stratakis et al, “Experimental and numerical study of phase mixing of an intense beam,” Phys. Rev. ST Accel. Beams 12, 064201 (2009).
  * [4] Y. Levin et al, “Nonequilibrium statistical mechanics of systems with long-range interactions,” Phys. Rep. 535, 1-60 (2014).
  * [5] B. B. Kadomtsev and O. P. Pogutse, “Collisionless relaxation in systems with Coulomb interactions,” Phys. Rev. Lett. 25, 1155 (1970).
  * [6] Y. Levin, R. Pakter, and T. Teles, “Collisionless Relaxation in Non-Neutral Plasmas,” Phys. Rev. Lett. 100, 040604 (2008).
  * [7] A. Gretton et al, “A Kernel Two-Sample Test”, Journal of Machine Learning Research 13, 723-773 (2012).
  * [8] T. Hofmann, B. Scholkopf, and A. J. Smola, “Kernel Methods in Machine Learning,” The Annals of Statistics 36, 1171-1220 (2008).
  * [9] V. Paulsen and M. Raghupathi, An Introduction to the Theory of Reproducing Kernel Hilbert Spaces, Cambridge University Press, 2016.
  * [10] A. Gretton et al, “Measuring Statistical Dependence with Hilbert-Schmidt Norms”, Technical report no. 140, Max Planck Institute for Biological Cybernetics (2005), http://www.gatsby.ucl.ac.uk/~gretton/papers/GreHerSmoBouSch05a.pdf.
  * [11] By a “probability density” we mean, more generally, a probability measure defined on the $\sigma$-algebra of Borel subsets of the phase space $\mathbb{R}^{d}$. Such a measure may not possess a density in the strict sense of a function in $L^{1}(\mathbb{R}^{d})$. For example, our “probability density” here is allowed to contain Dirac delta functions.
  * [12] R. M. Dudley, Real Analysis and Probability, Cambridge Studies in Advanced Mathematics 74, Cambridge University Press, 2002.
  * [13] P. Billingsley, Convergence of Probability Measures, 2nd ed, Wiley, NY, 1999.
  * [14] M. Markatou, D. Karlis, and Y. Ding, “Distance-Based Statistical Inference,” Annu. Rev. Stat. Appl. 8, 301-27 (2021).
  * [15] S. Kullback and R. Leibler, “On Information and Sufficiency,” Annals of Mathematical Statistics 22, 79-86 (1951).
  * [16] S. Kullback, Information Theory and Statistics, Dover Publications, Inc., Mineola, NY (1968).
  * [17] C. Granero-Belinchon, S. Roux, and N. Garnier, “Kullback-Leibler divergence measure of intermittency: Application to turbulence,” Phys. Rev. E 97, 013107 (2018).
  * [18] Q. Wang et al, “Divergence Estimation for Multidimensional Densities via $k$-Nearest-Neighbor Distances”, IEEE Transactions on Information Theory 55, 2392 (2009).
  * [19] Strictly speaking, the infimum of this quantity must be taken over all “joint” probability measures defined on the Borel subsets of $\mathbb{R}^{d}\times\mathbb{R}^{d}$ whose marginals have densities $f$ and $g$, respectively. See comment [11].
  * [20] L. V. Kantorovich, “On the Translocation of Masses”, Dokl. Akad. Nauk USSR 37, pp. 227-229 (1942). [English translation: Journal of Mathematical Sciences 133, 1381-1382 (2006)].
  * [21] L. N. Wasserstein, “Markov Processes Over Denumerable Products of Spaces Describing Large Systems of Automata,” Probl. Inform. Transmission 5, pp. 47-52 (1969).
  * [22] C. Villani, Topics in Optimal Transportation, American Mathematical Society, 2003.
  * [23] P. Komiske, E. Metodiev, and J. Thaler, “Metric Space of Collider Events,” Phys. Rev. Lett. 123, 041801, 2019.
  * [24] M. Erdmann el, “Generating and Refining Particle Detector Simulations Using the Wasserstein Distance in Adversarial Networks,” Computing and Software for Big Science 2:4 (2018).
  * [25] T. Cai, J. Cheng, and N. Craig, “Linearized Optimal Transport for Collider Events,” Phys. Rev. D 102, 116019 (2020).
  * [26] B. Sriperumbudur et al, “On the Empirical Estimation of Integral Probability Metrics”, Electronic Journal of Statistics 6, 1550-1599 (2012).
  * [27] B. Sriperumbudur et al, “Non-parametric Estimation of Integral Probability Metrics”, in Proc. 2010 IEEE International Symposium on Information Theory, 13-18 June 2010.
  * [28] B. Sriperumbudur et al, “Hilbert Space Embeddings and Metrics on Probability Measures”, Journal of Machine Learning Research 11, 1517-1561 (2010).
  * [29] B. Sriperumbudur, K. Fukumizu, and G. Lanckriet, “Universality, Characteristic Kernels and RKHS Embeddings of Measures,” Journal of Machine Learning Research 12, 2389-2410 (2011).
  * [30] B. Sriperumbudur, “On the Optimal Estimation of Probability Measures in the Weak and Strong Topologies,” Bernoulli 22, 1839-1893 (2016).
  * [31] C. J. Simon-Gabriel and B. Schölkopf, “Kernel Distribution Embeddings: Universal Kernels, Characteristic Kernels and Kernel Metrics on Distributions,” Journal of Machine Learning Research 19, 1-29 (2018).
  * [32] C. J. Simon-Gabriel et al, “Metrizing Weak Convergence with Maximum Mean Discrepancies” (2021), https://arxiv.org/abs/2006.09268.
  * [33] D. Sejdinovic et al, “Equivalence of Distance-Based and RKHS-Based Statistics in Hypothesis Testing”, The Annals of Statistics 41, 2263-2291 (2013).
  * [34] G. J. Székely and M. L. Rizzo, “Energy statistics: A class of statistics based on distances,” Journal of Statistical Planning and Inference 143, 1249-1272 (2013).
  * [35] J. Zhao and D. Meng, “FastMMD: Ensemble of circular discrepancy for efficient two-sample test”, Neural Computation 27, 1345-1372 (2015).
  * [36] A. Rahimi and B. Recht, “Random Features for Large-Scale Kernel Machines,” Adv. Neural Inf. Process. Syst. 20 (2007).
  * [37] A. Gretton et al, “Measuring Statistical Dependence with Hilbert-Schmidt Norms”, Lecture Notes in Computer Science, pp. 63-77 (2005).
  * [38] A. Gretton et al, “A Kernel Statistical Test of Independence”, in Proc. of the 20th International Conference on Neural Information Processing Systems, p. 585-592 (2007).
  * [39] G. Székely, G. Rizzo, and M. Bakirov, “Measuring and testing dependence by correlation of distances”, Ann. Stat. 35, 2769-2794 (2007).
  * [40] R. Lyons, “Distance Covariance in Metric Spaces”, The Annals of Probability 41, 3284-3305 (2013).
  * [41] Q. Zhang et al, “Large-scale kernel methods for independence testing”, Stat. Comput. 28, 113-130 (2018).
  * [42] C. Mitchell, “Weak Convergence to Equilibrium of Statistical Ensembles in Integrable Hamiltonian Systems”, Journal of Mathematical Physics 60, 052702 (2019).
  * [43] R. E. Meller et al, “Decoherence of Kicked Beams,” Technical Note SSC-N-360 (1987), https://lss.fnal.gov/archive/other/ssc/ssc-n-360.pdf
  * [44] P. Walters, An Introduction to Ergodic Theory, Springer-Verlag, New York, 1982.
  * [45] J. Qiang, “Symplectic multiparticle tracking model for self-consistent space-charge simulation,” Phys. Rev. Accel. Beams 20, 014203 (2017).
  * [46] S. Nietert, Z. Goldfeld, K. Kato, “Smooth $p$-Wasserstein distance: structure, empirical approximation, and statistical applications,” Proceedings of the 38th International Conference on Machine Learning, PMRL 139, 8172 (2021).
  * [47] H. Wendland, “Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree,” Adv. Comput. Math. 4, 389-396 (1995).
  * [48] A. Chernih and S. Hubbert, “Closed form representations and properties of the generalised Wendland functions,” J. Approx. Theory 177, 17-33 (2014).
  * [49] H. Minh, “Some Properties of Gaussian Reproducing Kernel Hilbert Spaces and Their Implications for Function Approximation and Learning Theory”, Constructive Approximation 32, 307-338 (2010).
  * [50] This is the definition of strong mixing. It is assumed that $A$ and $B$ are “nice” (Borel) subsets of the phase space $\mathbb{R}^{d}$, and the map $\mathcal{M}$ is Borel-measurable. This holds, for example, if $\mathcal{M}$ is continuous.