whenever $A\cap B=\emptyset$}.$

######  Lemma 6.3.

Let $A$ be a bounded Borel subset of $\mathbb{R}^{d}$ with non empty interior
such that $|\partial A|=0$. Then

(6.10) $\lim_{\varepsilon\to
0}\varepsilon^{d}n_{\varepsilon}(A)=\lim_{\varepsilon\to
0}\varepsilon^{d}\tilde{n}_{\varepsilon}(A)=\gamma_{d}|A|.$

Furthermore, for every $\varepsilon>0$ let $\Sigma_{\varepsilon}\subset A$ be
an optimal subset for $n_{\varepsilon}(A)$ (resp. for
$\tilde{n}_{\varepsilon}(A)$). Then the associated empirical measure
$\mu_{\varepsilon}=\frac{1}{n_{\varepsilon}}\sum_{i=1}^{n_{\varepsilon}}\delta_{x_{i}}$
converges tightly to the uniform probability density on $A$ as $\varepsilon\to
0$.

###### Proof.

First we show (6.10) when $A=Q_{a}$ for $a>0$. Observe that
$n_{\varepsilon}(Q_{a})=S(Q_{\frac{a}{\varepsilon}})$ while
$S(Q_{k_{\varepsilon}})\leq S(Q_{\frac{a}{\varepsilon}})\leq
S(Q_{k_{\varepsilon}+1})$ being $k_{\varepsilon}$ the integer part of
$a/\varepsilon$. Then the equality $\lim_{\varepsilon\to
0}\varepsilon^{d}n_{\varepsilon}(Q_{a})=\gamma_{d}a^{d}$ follows since, by
(6.7), we have $S(Q_{k_{\varepsilon}})\sim\gamma_{d}\,k_{\varepsilon}^{-d}$ as
$\varepsilon\to 0$.

On the other hand, for $\varepsilon<\delta<a$, we have
$n_{\varepsilon}(Q_{a-\delta})\leq{\tilde{n}_{\varepsilon}}(Q_{a})\leq
n_{\varepsilon}(Q_{a})$. Thus by applying the previous convergence to $Q_{a}$
and $Q_{a-\delta}$ and then sending $\delta\to 0$, we deduce that
$\lim_{\varepsilon\to
0}\varepsilon^{d}\,{\tilde{n}_{\varepsilon}}(Q_{a})=\gamma_{d}a^{d}.$

Let now $A$ be a Borel set with non empty interior $\mathring{A}$ such that
$|\partial A|=0$. We consider the family of hypercubes
$\mathcal{A}\vcentcolon=\left\\{{Q(x,a)\;\colon
x\in\overline{A},a>0}\right\\}$ being $Q(x,a)=x+Q_{a}$. As $\overline{A}$ is
compact, for every $\delta>0$ we may find a finite subfamily $\\{Q_{i}:i\in
I\\}\subset\mathcal{A}$ such that $\overline{A}\subset\cup_{i\in I}Q_{i}$ and
$\sum_{i\in I}|Q_{i}|\leq|\overline{A}|+\delta$. By the subadditivity property
of the set function $n_{\varepsilon}$ and by using the first step, we get

$\limsup_{\varepsilon\to 0}\varepsilon^{d}n_{\varepsilon}(A)\leq\sum_{i\in
I}\limsup_{\varepsilon\to
0}\varepsilon^{d}n_{\varepsilon}(Q_{i})\leq\gamma_{d}\sum_{i\in
I}|Q_{i}|\leq\gamma_{d}(|\overline{A}|+\delta)=\gamma_{d}(|A|+\delta),$

where we used also that $|\partial A|=0$; by sending $\delta\to 0$ we get

$\limsup_{\varepsilon\to
0}\varepsilon^{d}n_{\varepsilon}(A)\leq\gamma_{d}|A|.$

On the other hand by Besicovitch’s Covering Theorem (generalized by Morse to
families of hypercubes, see for instance [21]), there exists a countable
family $\left\\{{Q_{n}}\right\\}\subset\mathcal{A}$ such that
$Q_{n}\subset\mathring{A}$, $Q_{m}\cap Q_{n}=\emptyset$ if $m\neq n$, with
$\left|\mathring{A}\setminus\cup_{n}Q_{n}\right|=0.$ By the superadditivity of
the set function ${\tilde{n}_{\varepsilon}}$ and by applying first step to
each $Q_{n}$, we deduce that

$\liminf_{\varepsilon\to
0}\varepsilon^{d}{\tilde{n}}_{\varepsilon}(A)\geq\sum_{n}\liminf_{\varepsilon\to
0}\varepsilon^{d}n_{\varepsilon}(Q_{n})\geq\gamma_{d}\sum_{i\in
I}|Q_{i}|\geq\gamma_{d}|\mathring{A}|\ .$

Recalling that $\tilde{n}_{\varepsilon}(A)\leq n_{\varepsilon}(A)$ and that
$\partial A$ was assumed to be of vanishing Lebesque measure, we deduce the
converse inequality $\liminf_{\varepsilon\to
0}\varepsilon^{d}n_{\varepsilon}(A)\geq\gamma_{d}|A|,$ whence (6.10).

Let us now prove the second assertion of 6.3. Let $\Sigma_{\varepsilon}$ be an
optimal subset of $A$ associated with $n_{\varepsilon}(A)$. Up to a
subsequence we may assume that $\mu_{\varepsilon}$ converges tightly to a
probability measure $\mu$ supported on the compact set $\overline{A}$. For
every $x_{0}\in\mathbb{R}^{d}$ and $\delta>0$, by applying (6.10), we find
that:

$\mu(Q_{\delta}(x_{0}))\leq\liminf_{\varepsilon\to
0}\frac{\sharp(\Sigma_{\varepsilon}\cap
Q_{\delta}(x_{0}))}{n_{\varepsilon}(A)}\leq\limsup_{\varepsilon\to
0}\frac{n_{\varepsilon}(Q_{\delta}(x_{0})\cap
A)}{n_{\varepsilon}(A)}=\frac{|Q_{\delta}(x_{0})\cap A|}{|A|}.$

Since $x_{0}$ and $\delta$ are arbitrary, we infer that $\mu$ is an absolutely
continuous measure with a density $\frac{d\mu}{dx}\leq\frac{1}{|A|}$ a.e. in
$\overline{A}$. Thus as $\mu(\overline{A})=1$ (and $|\overline{A}|=|A|$), we
conclude that $\mu$ is the uniform probability density on $A$ and it is the
unique cluster point of $\mu_{\varepsilon}$ as $\varepsilon\to 0$. In view of
the first equality in (6.10), the previous arguments work as well when
substituting $n_{\varepsilon}$ with $\tilde{n}_{\varepsilon}$. ∎

######  Remark 6.4.

In view of 6.3, the relation between $\gamma_{d}$ given in (6.7) and the best
density of spheres packing constant in $\mathbb{R}^{d}$ is now clear since the
maximal volume fraction of non overlapping spheres of diameter $\varepsilon$
which can be placed in a regular subset $A$ of $\mathbb{R}^{d}$, namely
$\tilde{n}_{\varepsilon}(A)\,\omega_{d}\left({\frac{\varepsilon}{2}}\right)^{d}$
(being $\omega_{d}$ the volume of the unit sphere), is asymptotically equal to
$\gamma_{d}\,\frac{\omega_{d}}{2^{d}}$. The exact value of this volume
fraction is well known for $d\leq 3$, where the optimal configuration can be
recovered from a periodic lattice (regular hexagonal lattice for $d=2$). We
refer to [10] for a survey in more dimensions of the celebrated sphere packing
problem.

As a consequence of 6.3 we have

######  Lemma 6.5.

Let $\ell$ be the repulsive potential as defined in (6.4) and assume that
there exists $(\rho_{N})\in\mathcal{P}(\overline{\Omega})$ such that
$\limsup_{N\to\infty}F_{N}(\rho_{N})<+\infty$. Then the parameter $\kappa$
defined in (6.2) satisfies the inequality
$\kappa\leq\gamma_{d}\left|{\Omega}\right|$. Conversely, if the previous
inequality is strict, then
$\limsup_{N\to\infty}\inf\widetilde{F_{N}}<+\infty$.

###### Proof.

Since $F_{N}$ is an indicator function, the finiteness of $F_{N}(\rho_{N})$
implies that $\inf F_{N}=-F_{N}^{*}(0)=0$ which in view of (6.6) is equivalent
to the existence of a subset $\Sigma\subset\overline{\Omega}$ such that
$\sharp(\Sigma)=N$ and $|x-y|\geq\varepsilon_{N}$ for all
$(x,y)\in\Sigma^{2}\setminus\Delta.$ Therefore it is necessary that $N\leq
n_{\varepsilon_{N}}(\overline{\Omega})$ for $N$ large. By applying (6.10) to
$A=\overline{\Omega}$ and to the sequence $\varepsilon_{N}$, we get
$n_{\varepsilon_{N}}(\overline{\Omega})\sim\varepsilon_{N}^{-d}\gamma_{d}|\Omega|$
as $N\to\infty$, thus concluding to the desired inequality since, with the
help of (6.2), we have

$\lim_{N\to\infty}\frac{N}{n_{\varepsilon_{N}}(\overline{\Omega})}=\frac{\kappa}{\gamma_{d}|\Omega|}\leq
1.$

Conversely, if $\frac{\kappa}{\gamma_{d}|\Omega|}<1$, then the inequality
$N\leq n_{\varepsilon_{N}}(\overline{\Omega})$ holds for large $N$, ensuring
the existence of a set $\Sigma\subset\overline{\Omega}$ such that
$\sharp(\Sigma)=N$ and $\widetilde{F_{N}}(\mu_{\Sigma})=0$. ∎

###### Proof of 6.1.

Since $F_{N}\leq\widetilde{F_{N}}$, it is enough to establish the
$\Gamma$-$\liminf$ inequality for $(F_{N})$ and the $\Gamma$-$\limsup$
inequality for $(\widetilde{F_{N}})$. We proceed in two steps:

Step 1 ($\Gamma$-$\liminf$ inequality). As $F_{N}$ is an indicator function,
proving the $\Gamma$-$\liminf$ inequality amounts to show that, for any
sequence $(\rho_{N})$ such that $F_{N}(\rho_{N})=0$ and
$\rho_{N}\rightharpoonup\rho$ in $\mathcal{P}(\overline{\Omega})$, it holds
$\rho\in\mathcal{K}$. Let $\varphi\in C^{0}_{+}(\overline{\Omega})$. Then, by
applying the Fenchel inequality, we have:

(6.11)
$\left\langle{v,\rho}\right\rangle=\lim_{N\to\infty}\left\langle{v,\rho_{N}}\right\rangle\leq\limsup_{N\to\infty}F_{N}^{*}(v).$

By selecting an optimal subset $\Sigma_{N}\subset\overline{\Omega}$ in (6.6),
we deduce that

$F_{N}^{*}(v)=\left\langle{v,\mu_{\Sigma_{N}}}\right\rangle\ =\
\int_{0}^{+\infty}\frac{\sharp\left(\left\\{{v>t}\right\\}\cap\Sigma_{N}\right)}{N}dt,$

where in the last equality we used the fact that
$\mu_{\Sigma_{N}}(\left\\{{v>t}\right\\})=\frac{\sharp\left(\left\\{{v>t}\right\\}\cap\Sigma_{N}\right)}{N}$.

Next we invoke the definition of the set function $n_{\varepsilon_{N}}$ to
infer that:

$F_{N}^{*}(v)\leq\int_{0}^{+\infty}\frac{n_{\varepsilon_{N}}\left(\left\\{{v>t}\right\\}\right)}{N}dt.$

By observing that, for a.e. $t\in[0,\sup v]$, the open subset
$\left\\{{v>t}\right\\}$ is non empty with a Lebesgue negligible boundary, we
apply (6.10) to pass to the limit $N\to\infty$ in the previous equality with
the help of dominated convergence: rCl lim sup_N →∞ F_N^*(v) & ≤ lim_N→∞
∫_0^+∞ nεN({v ¿ t })N dt
= ∫_0^supv γdκ —{ v ¿ t }— dt = γdκ ∫_Ωv(x) dx.

Therefore, from (6.11), we see that the inequality
$\left\langle{v,\rho}\right\rangle\leq\frac{\gamma_{d}}{\kappa}\int_{\Omega}v(x)dx$
holds for every non negative continuous test function $v$. It follows that
$\rho$ is an absolutely continuous probability on $\overline{\Omega}$ with a
density $u$ such that $u\leq\frac{\gamma_{d}}{\kappa}$ a.e. This implies that
$\rho\in\mathcal{K}$ since
$\frac{\gamma_{d}}{\kappa}=\frac{1}{\theta|\Omega|}.$

Step 2. $\Gamma$-$\limsup$ inequality We need to show that
$\mathcal{K}\subset\mathcal{K}_{\infty}$ where

$\mathcal{K}_{\infty}\vcentcolon=\left\\{{\rho\in\mathcal{P}(\overline{\Omega})\;\colon\exists\Sigma_{N}\subset\overline{\Omega},F_{N}(\mu_{\Sigma_{N}})=0,\mu_{\Sigma_{N}}\rightharpoonup\rho}\right\\}.$

Since $\mathcal{K}_{\infty}$ is a closed subset of
$\mathcal{P}(\overline{\Omega})$ equipped with the tight convergence, it is
enough to show that $\mathcal{K}_{0}\subset\mathcal{K}_{\infty}$ being
$\mathcal{K}_{0}$ a dense subset of $\mathcal{K}$, namely (see 6.6 below):

(6.12) $\mathcal{K}_{0}=\left\\{{\rho\in\mathcal{P}(\overline{\Omega})\
\;\colon\ \rho=\left(\sum_{i\in I}t_{i}\mathbbm{1}_{A_{i}}\right)dx\ ,\
t_{i}\leq u^{*}}\right\\},$

where $u^{*}:=\frac{1}{\theta|\Omega|}$ and $\\{A_{i},i\in I\\}$ is a finite
family of disjoint open subsets such that $A_{i}\subset\subset\Omega$ and
$|\partial A_{i}|=0$.

Let $\rho\in\mathcal{K}_{0}$ be in the form above and denote
$\omega=\cup_{i\in I}A_{i}$, where $I=\left\\{{1,\dotsc,K}\right\\}.$ By
setting $A_{0}=\Omega\setminus\overline{\omega}$, we obtain a partition
$\cup_{i=0}^{K}A_{i}$ of full measure in $\Omega$. Note that the condition
$\sum_{i\in I}t_{i}|A_{i}|=1$ with $t_{i}\leq u^{*}$ implies that $1\leq
u^{*}|\omega|$. Thus the volume ratio of $\omega$ satisfies
$\frac{|\omega|}{|\Omega|}\geq\theta$ (recall that $\theta<1$).

Let us fix a parameter $\eta<1$ that ultimately will be sent to $1$. Then, for
every $i\in I$, we set:

(6.13) $\delta_{i}\vcentcolon=\left(\frac{u^{*}}{\eta
t_{i}}\right)^{1/d},\quad\varepsilon_{i,N}=\delta_{i}\varepsilon_{N}\quad\text{(thus
$\varepsilon_{i,N}>\varepsilon_{n}$)}\,.$

According to (6.10), for any $i\in I$, there exists a set $\Sigma_{i,N}\subset
A_{i}$ such that:

$\bigcup_{x\in\Sigma_{i,N}}B(x,\varepsilon_{i,N})\subset
A_{i},\quad\sharp(\Sigma_{i,N})=\tilde{n}_{\varepsilon_{i,N}}(A_{i})\sim\frac{\gamma_{d}}{\varepsilon_{N}^{d}\delta_{i}^{d}}|A_{i}|,\quad\mu_{\Sigma_{i,N}}\rightharpoonup\frac{\mathbbm{1}_{A_{i}}}{|A_{i}|}.$

Then we set $\Sigma^{\prime}_{N}:=\cup_{i\in I}\Sigma_{i,N}$ and, taking into
account (6.2), (6.13) and the equality $\sum_{i\in I}t_{i}|A_{i}|=1$, we
obtain:

$\lim_{N\to\infty}\frac{\sharp(\Sigma^{\prime}_{N})}{N}=\eta,\quad\frac{\sum_{x\in\Sigma^{\prime}_{N}}\delta_{x}}{N}\rightharpoonup\eta\rho.$

Since $\eta<1$, we know that $\sharp(\Sigma^{\prime}_{N})<N$ for large $N$ and
we need completing $\Sigma^{\prime}_{N}$ with a subset $\Sigma_{0,N}\subset
A_{0}$ such that $\sharp(\Sigma_{0,N})=N-\sharp(\Sigma^{\prime}_{N})$ and
$\\{B(x,\varepsilon_{N})\;\colon x\in\Sigma_{0,N}\\}$ is a family of disjoints
balls in $A_{0}$. This requires that
$N-\sharp(\Sigma^{\prime}_{N})\leq\tilde{n}_{\varepsilon_{N}}(A_{0})$. This
condition is indeed satisfied for $\eta$ close to $1$ and large $N$ since, by
(6.10), it holds $\lim_{N\to\infty}\frac{\sharp(\Sigma_{0,N})}{N}=1-\eta$
while
$\lim_{N\to\infty}\frac{\tilde{n}_{\varepsilon_{N}}(A_{0})}{N}=u^{*}|A_{0}|$.
Possibly passing to a subsequence, we may assume that
$\mu_{\Sigma_{0,N}}\rightharpoonup\rho_{0,\eta}$ where
$\rho_{0,\eta}\in\mathcal{P}(\overline{A_{0}})$.

Summarizing we have obtained a set
$\Sigma_{N}=\Sigma^{\prime}_{N}\cup\Sigma_{0,N}$ such that
$\sharp(\Sigma_{N})=N$ and
$\mu_{\Sigma_{N}}=(1-\eta)\mu_{\Sigma^{\prime}_{N}}+\eta\mu_{\Sigma_{0,N}}$.
In addition $\mu_{\Sigma_{N}}$ is admissible since by construction
$\\{B(x,\varepsilon_{n})\;\colon x\in\Sigma_{N}\\}$ is a family of disjoint
balls contained in $\Omega$. Therefore the weak limit of $\mu_{\Sigma_{N}}$ as
$N\to\infty$ namely $\rho_{\eta}\vcentcolon=\eta\rho+(1-\eta)\rho_{0,\eta}$
belongs to $\mathcal{K}_{\infty}$. Since this is true for any $\eta<1$ close
to $1$, by sending $\eta\to 1$ we conclude that $\rho\in\mathcal{K}_{\infty}$.
∎

######  Lemma 6.6.

Assume that $\theta$ given by (6.9) satisfies $\theta<1$ and let
$\mathcal{K}_{0}$ be the subset of $\mathcal{K}$ defined in (6.12). Then
$\mathcal{K}_{0}$ is weakly* dense in $\mathcal{K}$.

###### Proof.

Let us consider first $\rho\in\mathcal{K}$ of the form $\rho=udx$ where $u\in
C(\Omega,[0,u^{*}])$ ($u^{*}=\frac{1}{\theta|\Omega|}$) and
$\operatorname{supp}u\subset\subset\Omega$. Given $\delta>0$, we select an
increasing sequence of real numbers $t_{0}=0<t_{1}<\dots t_{k-1}<t_{k}=u^{*}$
such that $t_{i}-t_{i-1}\leq\delta$ and $|\\{u=t_{i}\\}|=0$ for all $i\geq 1$;
then we set $\Omega_{i}=\left\\{{u>t_{i}}\right\\}$ for $i\geq 0$ and
$A_{i}=\Omega_{i-1}\setminus\overline{\Omega_{i}}$ for $1\leq i\leq k$.

We see that the $A_{i}$’s are non empty disjoint open subsets of $\Omega$ such
that $\cup_{i=1}^{k}A_{i}$ is of full Lebesgue measure in
$\Omega_{0}=\left\\{{u>0}\right\\}$ and we have $\Omega_{0}\Subset\Omega$ by
our assumption on the support of $u$. Therefore we obtain an element
$\rho_{\delta}=u_{\delta}dx\in\mathcal{K}_{0}$ by setting:

$u_{\delta}\vcentcolon=\sum_{i=1}^{k}u_{i}\mathbbm{1}_{A_{i}}\quad\text{where}\quad
u_{i}=\frac{1}{|A_{i}|}\int_{A_{i}}u(x)dx.$

As $|u_{\delta}-u|\leq\delta$, we clearly have
$\rho_{\delta}\rightharpoonup\rho$ hence $\rho\in\overline{\mathcal{K}_{0}}$.

By using the classical convolution approximation argument (which preserves the
integral constraint $\int_{\Omega}u(x)dx=1$ and the inequality $u\leq u^{*}$),
we can deduce that the same conclusion applies to all $\rho\in\mathcal{K}$
with compact support in $\Omega$. The last step consists in showing that any
element $\rho=u\ dx$ in $\mathcal{K}$ can be approximated by a sequence
$\rho_{n}=u_{n}dx$ in $\mathcal{K}$ such that
$\operatorname{supp}u_{n}\Subset\Omega$. Thanks to the assumption
$u^{*}|\Omega|>1$, the non negative function $u^{*}-u$ has a positive integral
over $\Omega$. Then we may take a non decreasing sequence of compact sets
$K_{n}$ such that $\Omega=\cup_{n}K_{n}$ and consider the following sequence:

$u_{n}=\mathbbm{1}_{K_{n}}(u+s_{n}(u^{*}-u))\quad\text{where}\quad
s_{n}=\frac{\int_{\Omega\setminus K_{n}}udx}{\int_{K_{n}}(u^{*}-u)dx}.$

As $\lim_{n}\int_{K_{n}}(u^{*}-u_{n})dx=\int_{\Omega}(u^{*}-u)dx>0$, we have
$s_{n}\to 0$ and therefore $u_{n}\leq u^{*}$ for large $n$. In addition $\int
u_{n}dx=1$ while $\operatorname{supp}u_{n}\Subset\Omega$ by construction. As
$u_{n}\to u$ in $L^{1}(\Omega)$, we infer that
$u\in\overline{\mathcal{K}_{0}}$. ∎

##  Appendix A Tools and notation from convex analysis

Let $X$ be a topological vector space, and $f\colon X\to(-\infty,+\infty]$.
The _lower semi-continuous envelope_ of $f$, denoted $\overline{f}$ is the
greatest lower semi-continuous function below $f$, _i.e._ ,

$\overline{f}(x)=\sup\left\\{{g(x)\;\colon g\leq f,g\ \text{ is
l.s.c.}}\right\\}=\inf_{x_{n}\to x}\\{\liminf_{n\to\infty}f(x_{n})\\}.$

The _convex hull_ of $f$, denoted by $\operatorname{cl}f$, is the largest
convex lower semi-continuous function below $f$. It may be defined as the
function whose epigraph is the closed convex hull of the epigraph of $f$ in
$X\times\overline{\mathbb{R}}$. Notice that $f$ is lower semi-continuous iff
$\overline{f}=f$ and that if $f$ is convex then
$\operatorname{cl}f=\overline{f}$.

The _Legendre-Fenchel conjugate_ of $f$, denoted by $f^{*}$, is defined on
$X^{*}$ as

$f^{*}(v)=\sup\left\\{{\left\langle{v,x}\right\rangle-f(x)\;\colon x\in
X}\right\\},$

while its biconjugate , denoted by $f^{**}$, is given by

$f^{**}(x)=\sup\left\\{{\left\langle{v,x}\right\rangle-f^{*}(v)\;\colon v\in
X^{*}}\right\\}.$

We list some well-known properties (see for instance [3, 27])

(A.1) $\displaystyle\begin{cases}(a)&\text{$f^{*}$ is convex and lower semi-
continuous;}\\\ (b)&\text{ if $f\leq g$, then $f^{*}\geq g^{*}$; }\\\
(c)&\text{if $\\{f^{*}<+\infty\\}$ is non empty, then
$f^{**}=\operatorname{cl}f$ ;}\\\ (d)&\text{$(\overline{f})^{*}=f^{*}$ and
$(\overline{f})^{**}=f^{**}$;}\\\ (e)&\text{if $f$ is convex, then
$\overline{f}=f^{**}$.}\\\ \end{cases}$

Let $f\colon X\to(-\infty,+\infty]$ be a convex l.s.c. function and $x_{0}\in
X$ such that $f(x_{0})<+\infty.$ Then, for every $z\in X$, the map
$t\mapsto\frac{1}{t}(f(x_{0}+tz)-f(x_{0}))$ is non decreasing on
$(0,+\infty)$. We define the recession function of $f$ as follows:

$f_{\infty}(z):=\sup_{t>0}\frac{f(x_{0}+tz)-f(x_{0})}{t}=\lim_{t\to\infty}\frac{f(x_{0}+tz)-f(x_{0})}{t}.$

The next result shows that the definition above does not depend of $x_{0}$.

######  Lemma A.1.

$f_{\infty}$ coincides with the support function of the subset ${\rm
dom}\,f^{*}:=\\{v\in X^{*}\;\colon f^{*}(v)<+\infty\\}$, that is:

(A.2) $f_{\infty}(z)\ =\ \sup\\{\left\langle{v,x}\right\rangle\;\colon
v\in{\rm dom}\,f^{*}\\}.$

Accordingly $f_{\infty}$ is convex l.s.c. and positively one-homogeneous on
$X$.

##  Appendix B Some basic results in $\Gamma$-convergence theory

For an extensive study we refer to the monographs [1, 5, 13]. Recall that a
sequence of functionals $(F_{n})$ defined on a topological space $X$ is said
to be $\Gamma-$ convergent to $F:X\mapsto(-\infty,+\infty]$ if the two
following conditions hold:

  1. i)

($\Gamma-\liminf$) If $x_{n}\to x$ in $X$, then
$\liminf_{n\to\infty}F_{n}(x_{n})\geq F(x)$.

  2. ii)

($\Gamma-\limsup$) For all $x\in X$, there is a recovering sequence $(x_{n})$
in $X$ such that $x_{n}\to x$ and $\limsup_{n\to\infty}F_{n}(x_{n})\leq F(x)$.

It is easy to show that if $F_{n}\xrightarrow{\Gamma}F$, then $F$ is l.s.c on
$X$ . Furhermore we have the following equivalence:

$F_{N}\xrightarrow{\Gamma}F\iff\overline{F_{n}}\xrightarrow{\Gamma}F.$

As a consequence, one can readily check that, if $(F_{n})$ is monotone non
decreasing, then we have $F_{N}\xrightarrow{\Gamma}F$ where
$F=\sup_{n}\overline{F_{n}}.$

In addition we will use the two following results:

######  Proposition B.1 ([13, Proposition 5.9]).

Let $X$ be a separable Banach space, $(G_{n})_{n\in\mathbb{N}}$ a sequence of
equi-Lipschitz functionals on $X$. Then:

$G_{n}\xrightarrow{\Gamma}G\iff G_{n}(x)\to G(x)\ ,\ \forall x\in X.$

######  Proposition B.2.

Let $X$ be a separable Banach space and $X^{*}$ its dual endowed with the
weak* topology. Let $(F_{n})_{n\in\mathbb{N}}$ be a sequence of convex
functionals on $X^{*}$ satisfying the equicoercivity condition

$\sup F_{n}(\rho_{n})<+\infty\ \Rightarrow\
\sup_{n}\left\lVert{\rho_{n}}\right\rVert_{X^{*}}<+\infty.$

Then the following equivalence hold:

$F_{N}\xrightarrow{\Gamma}F\quad\text{in $X^{*}$}\ \iff\
F_{N}^{*}\xrightarrow{\Gamma}F^{*}\quad\text{in $X$}.$

###### Proof.

See [1, Theorem 3.11] in the reflexive case and D. Azé [2] for a
generalization to the non-reflexive case . ∎

## References

  * [1] H. Attouch.  Variational convergence for functions and operators.  Applicable Mathematics Series. Pitman (Advanced Publishing Program), Boston, MA, 1984.
  * [2] D. Azé.  Convergence des variables duales dans des problèmes de transmission à travers des couches minces par des méthodes d’épi-convergence.  Ricerche Mat., 35(1):125–159, 1986.
  * [3] G. Bouchitté.  Convex Analysis and Duality.  Encyclopedia of Mathematical physics. Academic Press, 2006.
  * [4] Guy Bouchitté, Giuseppe Buttazzo, Thierry Champion, and Luigi De Pascale.  Relaxed multi-marginal costs and quantization effects.  In Annales de l’Institut Henri Poincaré C, Analyse non linéaire. Elsevier, 2020.
  * [5] Andrea Braides et al.  Gamma-convergence for Beginners, volume 22.  Clarendon Press, 2002.
  * [6] Jan Brandts, Sergey Korotov, and Michal Křížek.  Simplicial partitions with applications to the finite element method.  Springer Monographs in Mathematics. Springer, Cham, [2020] ©2020.
  * [7] Giuseppe Buttazzo, Thierry Champion, and Luigi De Pascale.  Continuity and estimates for multimarginal optimal transportation problems with singular costs.  Applied Mathematics & Optimization, 78(1):185–200, 2018.
  * [8] J. A. Carrillo, J. Mateu, M. G. Mora, L. Rondi, L. Scardia, and J. Verdera.  The equilibrium measure for an anisotropic nonlocal energy.  Calc. Var. Partial Differential Equations, 60(3):Paper No. 109, 28, 2021.
  * [9] Gustave Choquet.  Diamètre transfini et comparaison de diverses capacités.  Séminaire Brelot-Choquet-Deny. Théorie du potentiel, 3, 1958-1959.  talk:4.
  * [10] J. H. Conway and N. J. A. Sloane.  Sphere packings, lattices and groups, volume 290 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences].  Springer-Verlag, New York, third edition, 1999.  With additional contributions by E. Bannai, R. E. Borcherds, J. Leech, S. P. Norton, A. M. Odlyzko, R. A. Parker, L. Queen and B. B. Venkov.
  * [11] Codina Cotar, Gero Friesecke, and Brendan Pass.  Infinite-body optimal transport with coulomb cost.  Calculus of Variations and Partial Differential Equations, 54(1):717–742, 2015.
  * [12] Codina Cotar and Mircea Petrache.  Next-order asymptotic expansion for n-marginal optimal transport with coulomb and riesz costs.  Advances in Mathematics, 344:137–233, 2019.
  * [13] Gianni Dal Maso.  An introduction to $\Gamma$-convergence, volume 8.  Springer Science & Business Media, 2012.
  * [14] Persi Diaconis and David Freedman.  Finite exchangeable sequences.  The Annals of Probability, pages 745–764, 1980.
  * [15] Otto Frostman.  Potentiel d’équilibre et capacité des ensembles avec quelques applications à la théorie des functions,lund.  Doctoral thesis, 115 s., 1935.
  * [16] Thomas Leblé and Sylvia Serfaty.  Large deviation principle for empirical fields of log and Riesz gases.  Invent. Math., 210(3):645–757, 2017.
  * [17] Christian Licht and Gérard Michaille.  Global-local subadditive ergodic theorems and application to homogenization in elasticity.  Annales Mathématiques Blaise Pascal, 9(1):21–62, 2002.
  * [18] Elliott H. Lieb and Barry Simon.  The Thomas-Fermi theory of atoms, molecules and solids.  Advances in Math., 23(1):22–116, 1977.
  * [19] Simone Di Marino, Mathieu Lewin, and Luca Nenna.  Grand-canonical optimal transport.  arXiv.2201.06859, 2022.
  * [20] Bertrand Maury.  Congested transport and microscopic and macroscopic scales.  In European Congress of Mathematics, pages 427–442. Eur. Math. Soc., Zürich, 2018.
  * [21] Anthony P. Morse.  Perfect blankets.  Trans. Amer. Math. Soc., 61:418–442, 1947.
  * [22] Brendan Pass.  Optimal transportation with infinitely many marginals.  Journal of Functional Analysis, 264(4):947–963, 2013.
  * [23] Mircea Petrache and Sylvia Serfaty.  Next order asymptotics and renormalized energy for Riesz interactions.  J. Inst. Math. Jussieu, 16(3):501–569, 2017.
  * [24] Sylvia Serfaty.  Coulomb gases and Ginzburg-Landau vortices.  Zurich Lectures in Advanced Mathematics. European Mathematical Society (EMS), Zürich, 2015.
  * [25] Sylvia Serfaty.  Systems of points with Coulomb interactions.  In Proceedings of the International Congress of Mathematicians—Rio de Janeiro 2018. Vol. I. Plenary lectures, pages 935–977. World Sci. Publ., Hackensack, NJ, 2018.
  * [26] Monica Skoge, Aleksandar Donev, Frank H. Stillinger, and Salvatore Torquato.  Packing hyperspheres in high-dimensional Euclidean spaces.  Phys. Rev. E (3), 74(4):041127, 11, 2006.
  * [27] Constantin Zalinescu.  Convex analysis in general vector spaces.  World scientific, 2002.


aainstitutetext: Department of Physics and Astronomy, University of
California, Irvine, CA 92697, USAbbinstitutetext: Physics Division, Lawrence
Berkeley National Laboratory, Berkeley, CA 94720, USAccinstitutetext: Berkeley
Institute for Data Science, University of California, Berkeley, CA 94720,
USAddinstitutetext: Jagellonian University, Krakow, Poland

# Towards a Deep Learning Model for Hadronization

Aishik Ghosh b Xiangyang Ju b,c Benjamin Nachman d and Andrzej Siodmok

###### Abstract

Hadronization is a complex quantum process whereby quarks and gluons become
hadrons. The widely-used models of hadronization in event generators are based
on physically-inspired phenomenological models with many free parameters. We
propose an alternative approach whereby neural networks are used instead. Deep
generative models are highly flexible, differentiable, and compatible with
Graphical Processing Unit (GPUs). We make the first step towards a data-driven
machine learning-based hadronization model by replacing a component of the
hadronization model within the Herwig event generator (cluster model) with a
Generative Adversarial Network (GAN). We show that a GAN is capable of
reproducing the kinematic properties of cluster decays. Furthermore, we
integrate this model into Herwig to generate entire events that can be
compared with the output of the public Herwig simulator as well as with
$e^{+}e^{-}$ data.

##  1 Introduction

Simulations are essential tools for nearly all aspects of data analysis at
particle colliders (see e.g., Ref. Buckley:2011ms ). These simulations are
rooted in particle and nuclear physics and must model a large range in energy
scales. At the smallest distance scales, various forms of perturbation theory
offer accurate, first-principles descriptions of hard-scatter particle
reactions and collinear parton shower radiation. The conversion from quarks
and gluons to hadrons is performed using hadronization models. Such approaches
are physically inspired but are ultimately phenomenological models with many
parameters that must be fit to data. There are currently two main
hadronization models, each inspired by a different description of strong
dynamics in the low-energy region. The linear confining potential motivated
the string model Andersson:1983ia ; Sjostrand:1984ic  implemented in Pythia
Sjostrand:2007gs ; Sjostrand:2006za  and preconfinement Amati:1979fg ;
Bassetto:1979vy  inspired the cluster model Webber:1983if  in Herwig
Corcella:2000bw ; Bahr:2008pv ; Bellm:2015jjp ; Bellm:2019zci  and Sherpa
Gleisberg:2008ta ; Sherpa:2019gpd . In both models, there is an intermediate
object between quarks/gluons and hadrons. This intermediate object (string or
cluster) takes as input the kinematic and flavor information from quarks and
gluons and then has an approximately universal fragmentation into different
hadron species that carry some fraction of the object’s momentum.

While existing hadronization models have been used successfully in a large
number of phenomenological and experimental studies at the Large Hadron
Collider and beyond, there is also significant room for innovation. Existing
models are not flexible enough to describe all of the properties of
hadronization (see e.g. Ref. 2004.03540 ). Even so, these models still have a
large number of parameters that need to be fit to data, which are adjusted
(‘tuned’) using semi-automated programs like Professor Buckley:2009bj .
Existing tuning methods are not able to process high-dimensional observables
or simultaneously tune many parameters because they rely on relatively simple
surrogate models to approximate the dependence of the data on the model. A
number of recently proposed automated tuning approaches employ sophisticated
surrogate models Ilten:2016csi ; Andreassen:2019nnm ; Wang:2021gdl , but they
all still require approximating complex relationships in high dimensions and
therefore often are limited to relatively low-dimensional parameter spaces.

One natural alternative to the existing hadronization simulations is deep
generative modeling. Machine learning-based generators are highly flexible and
differentiable by construction, which can aid parameter tuning. Three standard
approaches to deep generative models include Generative Adversarial Networks
(GANs) Goodfellow:2014:GAN:2969033.2969125 ; Creswell2018 , (Variational)
Autoencoders (VAEs) kingma2014autoencoding ; Kingma2019 , and Normalizing
Flows (NFs) 10.5555/3045118.3045281 ; Kobyzev2020 . While first proposed in
high energy physics (HEP) to emulate an entire parton shower
deOliveira:2017pjk  or detector simulations Paganini:2017dwg ;
Paganini:2017hrr , deep generative models have now been proposed for many
aspects of HEP simulations including matrix element generation
Bendavid:2017zhk ; Butter:2019cae ; Alanazi:2020jod ; Bothmann:2020ywa ;
Gao:2020zvv ; Gao:2020vdv ; Butter:2021csz , parton showers deOliveira:2017pjk
; Monk:2018zsb ; Carrazza:2019cnt ; Kansal:2020svm ; Lai:2020byl ;
Kansal:2021cqp ; Orzari:2021suh ; Tsan:2021brw ; Touranakou:2022qrp , detector
simulation Paganini:2017hrr ; Paganini:2017dwg ; Vallecorsa:2019ked ;
SHiP:2019gcl ; Chekalina:2018hxi ; ATL-SOFT-PUB-2018-001 ; Carminati:2018khv ;
Vallecorsa:2018zco ; Musella:2018rdi ; Erdmann:2018kuh ; Deja:2019vcv ;
Derkach:2019qfk ; Erdmann:2018jxd ; Oliveira:DLPS2017 ; deOliveira:2017rwa ;
Hooberman:DLPS2017 ; Belayneh:2019vyx ; buhmann2020getting ; 2009.03796 ;
Maevskiy:2020ank ; deja2020endtoend ; 1816035 ; Rehm:2021zow ; Rehm:2021zoz ;
Rehm:2021qwm ; Khattak:2021ndw ; Anderlini:2021qpm ; Fanelli:2019qaq ;
Lu:2020npg ; Buhmann:2021caf ; Krause:2021ilc ; Hariri:2021clz ;
Buhmann:2021lxj ; ATLAS:2021pzo ; Krause:2021wez ; Bieringer:2022cbs , and
more (see Ref. Butter:2022rso ; Butter:2020tvl ; Feickert:2021ajf  for
reviews). Using neural networks for modeling non-pertrubative inputs has a
long history in the context of Parton Distribution Functions (PDFs) (Ref.
Forte:2002fg  through Ref. Ball:2021leu ). Similarly to hadronization models,
PDFs cannot be calculated using perturbation theory. In contrast to
hadronization, PDFs are modeled as deterministic functions that are evolved in
energy scale using perturbation theory Gribov:1972ri ; Dokshitzer:1977sg ;
Altarelli:1977zs .

On the path towards a fully flexible, data-optimized, machine learning-based
hadronization model, we demonstrate the first step by training a GAN to mimic
a component of the cluster hadronization implementation in Herwig. In
particular, we replace part of the cluster decayer inside Herwig with a GAN
using the Open Neural Network Exchange (ONNX) bai2019  interface to call the
neural network inside the C++ code. This GAN-based cluster decayer, HADML, is
trained on Herwig. Future work will add additional complexity (cluster to
cluster decays, color reconnection of clusters Gieseke:2012ft ;
Gieseke:2019wcl ; Bellm:2019wrh , etc.) and will ultimately lead to a model
that can be trained (tuned) on data. This ultimate model will benefit from
new, high-dimensional future measurements Arratia:2021otl  that will provide
the necessary constraining power for the flexible neural network approaches.

This paper is organized as follows. Section 2 briefly introduces details of
the Herwig Monte Carlo event generator and how we interface a GAN in the
hadronization stage. Then, Sec. 2 presents the first numerical results with
the HADML hadronization model. The paper ends with conclusions and outlook in
Sec. 4.

##  2 Methods

###  2.1 Dataset

The training data was created using the hadronization cluster model
Webber:1983if  . The cluster model is based on t’Hooft’s planar diagram theory
tHooft:1973alw : the dominant color structure of Quantum Chromodynamics (QCD)
diagrams in the perturbation expansion in $1/N_{\mathrm{c}}$ can be
represented in a planar form using color lines, which is commonly known as the
limit $N_{\mathrm{c}}\to\infty$. The resulting color topology in Monte Carlo
events with partons in the final-state color features open color lines after
the parton showers. Following a non-perturbative isotropic decay of any left
gluons in the parton jets to quark-antiquark pairs, the event finally consists
of color-connected partons in color triplet or anti-triplet states. These
parton pairs form color-singlet clusters. This is so-called color
preconfinement Amati:1979fg : the tendency of the partons generated in the
parton shower to be arranged in color singlet clusters (pre-hadrons) with
limited extension in both coordinate and momentum space. The principle of
color preconfinement states that the mass distribution of these clusters is
independent of the hard-scattering process and its center-of-mass energy. The
cluster mass spectrum is not only universal but also peaked at low masses;
therefore, most of the clusters decay into two hadrons and some just into one
hadron. However, there is a small fraction of clusters that are too heavy for
this to be a reasonable approach. Therefore, these heavy clusters are first
split into lighter clusters before they decay. Such decays of massive clusters
are beyond the scope of this publication, and we will consider it in future
work. Since the kinematics of a cluster decaying into a single hadron is
trivial, our training data set only includes cases of decay into two hadrons.
To further simplify the training data, we consider only decays into pairs of
$\pi^{0}$. Each decay in our data set was described with the following
information: the four-momentum of the cluster, the four-momenta of the two
hadrons together with their flavor (encoded as a Particle Data Group (PDG)
Tanabashi:2018oca  code), and the Pert flag. Pert = 1 means that hadrons that
contain a parton produced in the perturbative stage of the event remember the
direction of the parton in the rest frame of the cluster. To create the
training data, we used $e^{+}e^{-}$ collisions at $\sqrt{s}=91.2$ GeV
generated by Herwig version 7.2.1. The only modification to the default
generator settings was the change that the hadrons produced from cluster
decays were on the mass shell***This setting can be achieved by adding the
command: set ClusterDecayer:OnShell Yes in the input file..

###  2.2 GAN Model and Training

We trained a conditional GAN to simulate the cluster decays. In a GAN, there
is a Generator neural network (Generator for short) and a Discriminant neural
network (Discriminator for short). Inputs to the Generator are the cluster’s
four vectors ($E$, $p_{x}$, $p_{y}$, $p_{z}$), and $N$ features sampled from a
Gaussian distribution. The $N$ numbers are called noise. $N$ is a
hyperparameter and set to be 10. Outputs of the Generator are the polar angle,
$\phi$, and azimuthal angle, $\theta$, of the leading hadron’s momentum in the
spherical coordinate system in the cluster frame, in which the two hadrons are
created back-to-back. With the two angular variables, $\theta$ and $\phi$, and
the cluster’s four vector, we reconstruct the four vectors of the two outgoing
hadrons as a postprocessing step. Inputs to the Discriminator are just the two
angular variables coming from either the Generator, labeled as background, or
those from the Herwig, labeled as signal. The output of the Discriminator is a
score that is higher for events from the Herwig and lower for events from the
Generator. The Discriminator is trained to separate signal from background.
However, the Generator is trained to yield signal-like Discriminator score.

The GAN is based on multilayer perceptrons (MLPs). Both the Generator and the
Discriminator are composed of a two-layer perceptron. Each perceptron consists
of a sequence of Keras keras  modules: a fully connected (dense) network of a
hidden size of 256, a batch normalization layer, and a LeakyReLU activation
function leakyReLU . These parameters were not extensively optimized.

To help train a GAN, we preprocessed the training data. The incoming cluster’s
four vector is scaled so that their values are between -1 and 1; so are the
two angular variables ($\phi$ and $\theta$). In this way, all inputs and
outputs are within the same scale. Finally, we use the $\tanh$ activation
function as the last layer of the Generator. The Discriminator and the
Generator are trained separately and alternately by two independent Adam
optimizers adam , both with a learning rate of $10^{-4}$, for about 1000
epochs.

Figure 1: Generator loss and discriminator loss and progressive best
Wasserstein distance as a function of the training epochs for training a GAN
with events where two partons are with Pert = 0. Both Generator and
Discriminator loss are the binary-crossentropy loss, and the Discriminator
loss is divided by two for visualization purposes. The progressive Wasserstein
distance is gauged by the right side of the $y$ axis.

Figure 1 shows the evolution of the Discriminator loss, which is divided by
two for visualization purposes, the Generator loss, and the progressive best
total Wasserstein distances†††This is a common metric in machine learning that
quantifies the minimal ‘work’ required to transform one density into another,
where work, in this case, is defined as the integral of the density multiplied
by the distance moved.  Vaserstein ; Kantorovich  for training a GAN with
events where two partons are with Pert = 0. The total Wasserstein distance
summing over the distances of all variables, is calculated after training for
one epoch and only the smallest value is plotted. At the beginning of the
training (epoch $<70$), even though the Generator loss is going up, we see a
rapid drop in the Wasserstein distance until the Generator loss is beyond 0.8.
For more than 100 epochs, the Discriminator keeps outperforming the Generator
as seen by the increasing Generator loss and the decreasing Discriminator
loss. This situation is changed around epoch 200 and finally, the two networks
reach an equilibrium around epoch 250. Beyond epoch 600, we only see about
0.002 improvements in the Wasserstein distance. The best model for events with
partons of Pert = 0, is found at the epoch 849 with a total Wasserstein
distance of 0.0228. A similar analysis was performed when training events with
at least one parton with Pert = 1.

###  2.3 Integration into Herwig

Each part of Herwig is implemented as a C++ class that contains the
implementation of the Herwig physics models, inheriting from an abstract base
class in ThePEG Lonnblad:2006pt . The ClusterHandronizationHandler is the
class that controls the cluster hadronization model. Our ultimate goal will be
to replace the entire ClusterHandronizationHandler with its ML counterpart.
However, since in these studies, we concentrate on the decay of clusters into
two hadrons, it was sufficient to modify ClusterDecayer \- a helper class of
the ClusterHandronizationHandler that controls this process. The generative
model trained in Python using TensorFlow is converted into the ONNX format
bai2019  and integrated into the Herwig chain using the C++ API of ONNX
Runtime onnxruntime . The advent of the ONNX format makes it possible to train
a model in one software and hardware environment and then apply it in a
completely different environment. ONNX Runtime is well suited for running fast
neural network inference as part of a large C++ workflow, and by using it, we
avoid having to integrate and maintain TensorFlow AbaBar16Tensorflow  within
the Herwig framework.

All preprocessing and postprocessing steps performed for training are repeated
within Herwig for inference. The entire simulation chain, including the GAN,
is then run in Herwig in order to produce the final comparisons and results.

##  3 Results

Figure 2:  Pseudorapidity (left panels) and transverse momentum (right panels)
distribution of $\pi^{0}$ from decays of Pert=0 (upper panels) and Pert=1
(lower panels) clusters produced in $e^{+}e^{-}$ collisions at $\sqrt{s}=91.2$
GeV.

Figure 3:  Pseudorapidity (left panels) and transverse momentum (right panels)
distribution of $\pi^{0}$ from decays of Pert=0 (upper panels) and Pert=1
(lower panels) clusters produced in $e^{+}e^{-}$ collisions at $\sqrt{s}=192$
GeV.

Section 3.1 provides low-level results of individual cluster decays while Sec.
3.2 includes full event simulations and comparisons to $e^{+}e^{-}$ data.

###  3.1 Low-level Validation

Since the training data contained only clusters produced in $e^{+}e^{-}$
collisions at $\sqrt{s}=91.2$ GeV that decayed into $\pi^{0}$ pairs, we begin
by comparing the $\pi^{0}$ kinematic variables generated by HADML and Herwig
precisely in such decays. The data generated by Herwig, with which we compared
the results of HADML in this section, were not used for training. In Fig. 2 we
show the distribution of the pseudorapidity (left panels) and transverse
momentum distribution (right panels) of $\pi^{0}$ from the decays of the Pert
= 0 (upper panels) and Pert=1 (lower panels) clusters. As expected, we see
that the transverse momentum spectra of pions coming from clusters containing
“perturbative” quarks (Pert=1) are harder compared to those containing only
non-perturbative partons (Pert=0). However, the most important observation
from Fig. 2 is that Herwig 7 \+ HADML (labeled on figures as H7+HADML) matches
the pseudorapidity of the pions generated by Herwig 7 with the cluster model
(labeled as H7 on figures). Transverse momentum spectra that extend over
several orders of magnitude are also well approximated by H7+HADML. Taking a
closer look at these distributions, we see minor differences for low
transverse momenta in the case of clusters that have a memory of perturbative
quarks (bottom-left panel in Fig. 2). Such small differences are, of course,
acceptable, especially since the information about the four-momentum of
partons that make up the clusters were not used for training. Taking this
additional information into account in the training process will likely
eliminate these minor differences. However, this is beyond the scope of this
publication, and we will leave this problem for future work.

It is crucial that the hadronization model is universal, i.e., that it works
independently of the hard process or collision energy. As we described in the
Sec. 2.1 the cluster model has this property. To test whether HADML also is
universal, we decided to repeat the comparison made at the beginning of this
section, but this time generating events with collision energies twice as high
as those used in the training data. In Fig. 3 we show $\pi^{0}$ kinematic
variables generated by H7+HADML and Herwig 7 in $e^{+}e^{-}$ collisions at
$\sqrt{s}=192$ GeV. We can see that all distributions are described very
similarly by both models, which reassured us that the HADML model is also
universal.

The last thing we need to check before using HADML to simulate the decay of
all clusters into hadron pairs in Herwig is whether the model is able to
describe the kinematics of other hadrons than $\pi^{0}$. In Fig. 4 we present
the pseudorapidity (left panels) and transverse momentum (right panels)
distribution of $\pi^{\pm}$ and $\pi^{0}$ (first row), kaons (second row) and
lambdas (third row). We see that the distributions differ for the various
hadrons, but they are all described almost identically by both models. This
encouraged us to perform a comparison with experimental data in which the
kinematics of all hadrons‡‡‡Except for a small number of hadrons that come
from the decay of a cluster into a single hadron for which the kinematics is
trivial. in Herwig are generated by HADML model.

Figure 4:  Pseudorapidity (left panels) and transverse momentum (right panels)
distribution of $\pi^{\pm}$ and $\pi^{0}$ (first row), Kaons (second row) and
Lambdas (third row).

###  3.2 Full-event Validation

In this section, we generate full events using HADML integrated into Herwig
and compare the results also to data from LEP§§§Note that the data are for
illustration only - given that the GAN is trained on Herwig, we cannot expect
it to outperform Herwig. Tuning to data is a longer-term goal of this research
(see Sec. 4).. In particular, we consider an analysis from DELPHI with data
collected at $\sqrt{s}=91.2$ GeV DELPHI:1996sen  using
RIVET¶¶¶https://rivet.hepforge.org/analyses/DELPHI_1996_S3430090.
Buckley:2010ar . These events correspond to hadronic $Z$ boson decays with a
number of event shape and identified hadron spectra. These data have been used
for hadronization parameter tuning DELPHI:1996sen ; Reichelt:2017hts .

Figure 5:  Normalized, differential cross-sections of Thurst (top left),
Thrust major (top right), Thrust minor (lower left), and Sphericity (lower
right) for Herwig, Herwig with HADML, and for data from DELPHI at LEP. Error
bars on the predictions represent statistical uncertainties.

Figure 5 shows histograms of various event shapes. Thrust Farhi:1977sg ;
Brandt:1964sa  is the quintessential $e^{+}e^{-}$ event shape:

$\displaystyle
T=\max_{\vec{n}}\left(\frac{\sum|\vec{p}_{i}\cdot\vec{n}|}{\sum|\vec{p}_{i}|}\right)\,,$
(1)

where the sum runs over all final state particle three momenta. The direction
$\vec{n}$ that maximizes the argument of Eq. 1 is called the Thrust axis.
Thrust major is defined similarly to Eq. 1 but with $\vec{n}$ replaced with
vectors transverse to the Thrust axis and Thrust minor is the same, but with
an optimization only over directions perpendicular to both the Thurst and
Thurst major axes. The Sphericity is computed from the eigenvalues of the
quadratic momentum tensor

$\displaystyle M^{\alpha\beta}=\sum p_{i}^{\alpha}p_{i}^{\beta}\,,$ (2)

where $\alpha,\beta$ are the spatial momentum indices, and the sum runs over
the same particles as in Eq. 1. Sphericity is defined as
$\frac{3}{2}(\lambda_{2}+\lambda_{3})$ for eigenvalues $\lambda_{i}$ of the
$3\times 3$ matrix defined in Eq. 2 and
$\lambda_{3}\leq\lambda_{2}\leq\lambda_{1}$. Hadronization shifts event shapes
(see e.g., Ref. Abbate:2010xh ) and so these observables are sensitive to
hadronization modeling. Figure 5 shows that HADML agrees with Herwig within
10% across most of the spectra, which itself agrees with data at a similar
level. Individual particle spectra are shown in Fig. 6 for the transverse
momenta along the Thurst major and minor directions. The level of agreement is
similar to the event shapes where there is sufficient statistical power.

Figure 6:  Normalized, differential cross-sections of particle transverse
momenta along the Thrust major (left) and Thurst minor (right) axes for
Herwig, Herwig with our HADML, and for data from DELPHI at LEP. Error bars on
the predictions represent statistical uncertainties.

##  4 Summary and Outlook

In this paper, we have established a first step on the path towards a neural
network-based hadronization model. The cluster hadronization model from Herwig
has been emulated with a Generative Adversarial Network. This model is
designed to reproduce the two-body decay of clusters into pions. The GAN is
integrated into the full Herwig program by using all other hadronization
components from the Herwig default model. The kinematic properties of other
hadrons are emulated using the pion model and conservation of energy. We have
shown that the HADML is able to reproduce Herwig’s light cluster decays and
when integrated with the full Herwig simulation, is able to reproduce results
from $e^{+}e^{-}$ data as well.

The ultimate goal of this research direction is to train the ML model directly
on data to improve upon the existing hadronization models. A number of
technical and methodological steps are required to achieve this vision. First,
the deep generative model needs to be extended to directly accommodate
multiple hadron species and to model the relative probabilities of the various
final states. In this work, we have modeled different hadron species using
conservation of energy, but this means that the fragmentation is assumed
universal. Architectural modifications could allow for perturbations on
universality. Hyperparameter optimization, including the investigation of
alternative generative models, is an important component of future work. Once
the deep generative model has the capacity to reproduce all of the physics of
the Herwig cluster model, methodological innovation is required to explore how
to tune the model to data. Traditionally, $e^{+}e^{-}$ data are used for
tuning. Optimization with a large set of one-dimensional, binned measurements
will need to be explored. A non-trivial aspect of this optimization is that
while the hadronization model would be differentiable, the parton shower input
would not be. Building in a model of uncertainty would also be a central
aspect of model tuning. It may also be possible to tune with unbinned, and
higher-dimensional results from $ep$ and $pp$ data H1:2021wkz ;
Arratia:2021otl ; Vandegar:2020yvw ; Bellagente:2019uyp ; 1800956 ;
Andreassen:2019cjw .

While we have focused on hadronization in the context of collider physics, the
ideas and concepts described in this paper have broader implications. First of
all, hadronization is used across high energy particle and nuclear physics
(see e.g., Ref. Sjostrand:2021dal ) and perturbations on the collider model
may be required to accurately describe other systems. Second, there are other
physical systems where first-principles input is combined with
phenomenological models. For example, a complete description of observational
cosmology requires an $N$-body simulation of the dark matter to be combined
with a description of visible matter around dark matter halos (see e.g., Ref.
2016MNRASFeng ; Modi2021A&C ; Dai:2020ekz ; Bohm:2020ilt ; Modi:2021acq ).
While different applications call for domain-specific adaptations, some
components and core methodology is common. Further development in this
research area will enable important advances in simulation to improve
inference in high energy physics and beyond.

Note added: As this manuscript was being finalized, we became aware of the
recent work in Ref. Ilten:2022jfm , which has a similar goal. That study uses
a different Monte Carlo program (Pythia instead of Herwig) and uses a
different generative model (Variational Autoencoder instead of a GAN).
Reference Ilten:2022jfm  also focuses on the pion-only case.

## Acknowledgments

The work of AS was funded by grant no. 2019/34/E/ST2/00457 of the National
Science Centre, Poland and the Priority Research Area Digiworld under the
program Excellence Initiative – Research University at the Jagiellonian
University in Cracow. BN and XJ were supported by the Department of Energy,
Office of Science under contract number DE-AC02-05CH11231. AG is supported by
the U.S. Department of Energy (DOE), Office of Science under Grant No. DE-
SC0009920.

## References

  * (1) A. Buckley et al., _General-purpose event generators for LHC physics_ , _Phys. Rept._ 504 (2011) 145–233, [1101.2599].
  * (2) B. Andersson, G. Gustafson, G. Ingelman and T. Sjostrand, _Parton Fragmentation and String Dynamics_ , _Phys. Rept._ 97 (1983) 31–145.
  * (3) T. Sjostrand, _Jet Fragmentation of Nearby Partons_ , _Nucl. Phys. B_ 248 (1984) 469–502.
  * (4) T. Sjöstrand, S. Mrenna and P. Z. Skands, _A Brief Introduction to PYTHIA 8.1_ , _Comput. Phys. Commun._ 178 (2008) 852–867, [0710.3820].
  * (5) T. Sjöstrand, S. Mrenna and P. Z. Skands, _PYTHIA 6.4 Physics and Manual_ , _JHEP_ 05 (2006) 026, [hep-ph/0603175].
  * (6) D. Amati and G. Veneziano, _Preconfinement as a Property of Perturbative QCD_ , _Phys. Lett. B_ 83 (1979) 87–92.
  * (7) A. Bassetto, M. Ciafaloni and G. Marchesini, _Color Singlet Distributions and Mass Damping in Perturbative QCD_ , _Phys. Lett. B_ 83 (1979) 207–212.
  * (8) B. R. Webber, _A QCD Model for Jet Fragmentation Including Soft Gluon Interference_ , _Nucl. Phys. B_ 238 (1984) 492–528.
  * (9) G. Corcella, I. G. Knowles, G. Marchesini, S. Moretti, K. Odagiri, P. Richardson et al., _HERWIG 6: An Event generator for hadron emission reactions with interfering gluons (including supersymmetric processes)_ , _JHEP_ 01 (2001) 010, [hep-ph/0011363].
  * (10) M. Bahr et al., _Herwig++ Physics and Manual_ , _Eur. Phys. J._ C58 (2008) 639–707, [0803.0883].
  * (11) J. Bellm et al., _Herwig 7.0/Herwig++ 3.0 release note_ , _Eur. Phys. J._ C76 (2016) 196, [1512.01178].
  * (12) J. Bellm et al., _Herwig 7.2 release note_ , _Eur. Phys. J. C_ 80 (2020) 452, [1912.06509].
  * (13) T. Gleisberg, S. Hoeche, F. Krauss, M. Schonherr, S. Schumann, F. Siegert et al., _Event generation with SHERPA 1.1_ , _JHEP_ 02 (2009) 007, [0811.4622].
  * (14) Sherpa collaboration, E. Bothmann et al., _Event Generation with Sherpa 2.2_ , _SciPost Phys._ 7 (2019) 034, [1905.09127].
  * (15) ATLAS Collaboration, _Measurement of the Lund jet plane using charged particles in 13 TeV proton-proton collisions with the ATLAS detector_ , _Phys. Rev. Lett._ 124 (2020) 222002, [2004.03540].
  * (16) A. Buckley, H. Hoeth, H. Lacker, H. Schulz and J. E. von Seggern, _Systematic event generator tuning for the LHC_ , _Eur. Phys. J. C_ 65 (2010) 331–357, [0907.2973].
  * (17) P. Ilten, M. Williams and Y. Yang, _Event generator tuning using Bayesian optimization_ , _JINST_ 12 (2017) P04028, [1610.08328].
  * (18) A. Andreassen and B. Nachman, _Neural Networks for Full Phase-space Reweighting and Parameter Tuning_ , _Phys. Rev. D_ 101 (2020) 091901, [1907.08209].
  * (19) W. Wang, M. Krishnamoorthy, J. Muller, S. Mrenna, H. Schulz, X. Ju et al., _BROOD: Bilevel and Robust Optimization and Outlier Detection for Efficient Tuning of High-Energy Physics Event Generators_ , _SciPost Phys. Core_ 5 (3, 2021) 1, [2103.05751].
  * (20) I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair et al., _Generative adversarial nets_ , in _Proceedings of the 27th International Conference on Neural Information Processing Systems - Volume 2_ , NIPS’14, (Cambridge, MA, USA), pp. 2672–2680, MIT Press, 2014, http://dl.acm.org/citation.cfm?id=2969033.2969125.
  * (21) A. Creswell, T. White, V. Dumoulin, K. Arulkumaran, B. Sengupta and A. A. Bharath, _Generative adversarial networks: An overview_ , _IEEE Signal Processing Magazine_ 35 (Jan, 2018) 53.
  * (22) D. P. Kingma and M. Welling, _Auto-encoding variational bayes_ , 1312.6114.
  * (23) D. P. Kingma and M. Welling, _An Introduction to Variational Autoencoders_ , _Foundations and Trends in Machine Learning_ 12 (2019) 307.
  * (24) D. J. Rezende and S. Mohamed, _Variational inference with normalizing flows_ , _International Conference on Machine Learning_ 37 (2015) 1530.
  * (25) I. Kobyzev, S. Prince and M. Brubaker, _Normalizing Flows: An Introduction and Review of Current Methods_ , _IEEE Transactions on Pattern Analysis and Machine Intelligence_ (2020) 1.
  * (26) L. de Oliveira, M. Paganini and B. Nachman, _Learning Particle Physics by Example: Location-Aware Generative Adversarial Networks for Physics Synthesis_ , _Comput. Softw. Big Sci._ 1 (2017) 4, [1701.05927].
  * (27) M. Paganini, L. de Oliveira and B. Nachman, _CaloGAN : Simulating 3D high energy particle showers in multilayer electromagnetic calorimeters with generative adversarial networks_ , _Phys. Rev._ D97 (2018) 014021, [1712.10321].
  * (28) M. Paganini, L. de Oliveira and B. Nachman, _Accelerating Science with Generative Adversarial Networks: An Application to 3D Particle Showers in Multilayer Calorimeters_ , _Phys. Rev. Lett._ 120 (2018) 042003, [1705.02355].
  * (29) J. Bendavid, _Efficient Monte Carlo Integration Using Boosted Decision Trees and Generative Deep Neural Networks_ , 1707.00028.
  * (30) A. Butter, T. Plehn and R. Winterhalder, _How to GAN LHC Events_ , _SciPost Phys._ 7 (2019) 075, [1907.03764].
  * (31) Y. Alanazi et al., _AI-based Monte Carlo event generator for electron-proton scattering_ , 2008.03151.
  * (32) E. Bothmann, T. Janßen, M. Knobbe, T. Schmale and S. Schumann, _Exploring phase space with Neural Importance Sampling_ , 2001.05478.
  * (33) C. Gao, S. Höche, J. Isaacson, C. Krause and H. Schulz, _Event Generation with Normalizing Flows_ , _Phys. Rev. D_ 101 (2020) 076002, [2001.10028].
  * (34) C. Gao, J. Isaacson and C. Krause, _i-flow: High-Dimensional Integration and Sampling with Normalizing Flows_ , 2001.05486.
  * (35) A. Butter, T. Heimel, S. Hummerich, T. Krebs, T. Plehn, A. Rousselot et al., _Generative Networks for Precision Enthusiasts_ , 2110.13632.
  * (36) J. W. Monk, _Deep Learning as a Parton Shower_ , 1807.03685.
  * (37) S. Carrazza and F. A. Dreyer, _Lund jet images from generative and cycle-consistent adversarial networks_ , _Eur. Phys. J._ C79 (2019) 979, [1909.01359].
  * (38) R. Kansal, J. Duarte, B. Orzari, T. Tomei, M. Pierini, M. Touranakou et al., _Graph Generative Adversarial Networks for Sparse Data Generation in High Energy Physics_ , _34th Conference on Neural Information Processing Systems_ (11, 2020) , [2012.00173].
  * (39) Y. S. Lai, D. Neill, M. Płoskoń and F. Ringer, _Explainable machine learning of the underlying physics of high-energy particle collisions_ , 2012.06582.
  * (40) R. Kansal, J. Duarte, H. Su, B. Orzari, T. Tomei, M. Pierini et al., _Particle Cloud Generation with Message Passing Generative Adversarial Networks_ , 2106.11535.
  * (41) B. Orzari, T. Tomei, M. Pierini, M. Touranakou, J. Duarte, R. Kansal et al., _Sparse Data Generation for Particle-Based Simulation of Hadronic Jets in the LHC_ , in _38th International Conference on Machine Learning Conference_ , 9, 2021, 2109.15197.
  * (42) S. Tsan, R. Kansal, A. Aportela, D. Diaz, J. Duarte, S. Krishna et al., _Particle Graph Autoencoders and Differentiable, Learned Energy Mover’s Distance_ , in _35th Conference on Neural Information Processing Systems_ , 11, 2021, 2111.12849.
  * (43) M. Touranakou, N. Chernyavskaya, J. Duarte, D. Gunopulos, R. Kansal, B. Orzari et al., _Particle-based Fast Jet Simulation at the LHC with Variational Autoencoders_ , 2203.00520.
  * (44) S. Vallecorsa, F. Carminati and G. Khattak, _3D convolutional GAN for fast simulation_ , _Proceedings, 23rd International Conference on Computing in High Energy and Nuclear Physics (CHEP 2018): Sofia, Bulgaria, July 9-13, 2018_ 214 (2019) 02010.
  * (45) SHiP collaboration, C. Ahdida et al., _Fast simulation of muons produced at the SHiP experiment using Generative Adversarial Networks_ , 1909.04451.
  * (46) V. Chekalina, E. Orlova, F. Ratnikov, D. Ulyanov, A. Ustyuzhanin and E. Zakharov, _Generative Models for Fast Calorimeter Simulation.LHCb case_ , _CHEP 2018_ (2018) , [1812.01319].
  * (47) ATLAS collaboration, _Deep generative models for fast shower simulation in ATLAS_ , _ATL-SOFT-PUB-2018-001_ (Jul, 2018) .
  * (48) F. Carminati, A. Gheata, G. Khattak, P. Mendez Lorenzo, S. Sharan and S. Vallecorsa, _Three dimensional Generative Adversarial Networks for fast simulation_ , _Proceedings, 18th International Workshop on Advanced Computing and Analysis Techniques in Physics Research (ACAT 2017): Seattle, WA, USA, August 21-25, 2017_ 1085 (2018) 032016.
  * (49) S. Vallecorsa, _Generative models for fast simulation_ , _Proceedings, 18th International Workshop on Advanced Computing and Analysis Techniques in Physics Research (ACAT 2017): Seattle, WA, USA, August 21-25, 2017_ 1085 (2018) 022005.
  * (50) P. Musella and F. Pandolfi, _Fast and Accurate Simulation of Particle Detectors Using Generative Adversarial Networks_ , _Comput. Softw. Big Sci._ 2 (2018) 8, [1805.00850].
  * (51) M. Erdmann, L. Geiger, J. Glombitza and D. Schmidt, _Generating and refining particle detector simulations using the Wasserstein distance in adversarial networks_ , _Comput. Softw. Big Sci._ 2 (2018) 4, [1802.03325].
  * (52) K. Deja, T. Trzcinski and u. Graczykowski, _Generative models for fast cluster simulations in the TPC for the ALICE experiment_ , _Proceedings, 23rd International Conference on Computing in High Energy and Nuclear Physics (CHEP 2018): Sofia, Bulgaria, July 9-13, 2018_ 214 (2019) 06003.
  * (53) D. Derkach, N. Kazeev, F. Ratnikov, A. Ustyuzhanin and A. Volokhova, _RICH 2018_ , 1903.11788.
  * (54) M. Erdmann, J. Glombitza and T. Quast, _Precise simulation of electromagnetic calorimeter showers using a Wasserstein Generative Adversarial Network_ , _Comput. Softw. Big Sci._ 3 (2019) 4, [1807.01954].
  * (55) L. de Oliveira, M. Paganini and B. Nachman, _Tips and Tricks for Training GANs with Physics Constraints_ , 2017, https://dl4physicalsciences.github.io/files/nips_dlps_2017_26.pdf.
  * (56) L. de Oliveira, M. Paganini and B. Nachman, _Controlling Physical Attributes in GAN-Accelerated Simulation of Electromagnetic Calorimeters_ , _J. Phys. Conf. Ser._ 1085 (2018) 042017, [1711.08813].
  * (57) B. Hooberman, A. Farbin, G. Khattak, V. Pacela, M. Pierini, J.-R. Vlimant et al., _Calorimetry with Deep Learning: Particle Classification, Energy Regression, and Simulation for High-Energy Physics_ , 2017, https://dl4physicalsciences.github.io/files/nips_dlps_2017_15.pdf.
  * (58) D. Belayneh et al., _Calorimetry with Deep Learning: Particle Simulation and Reconstruction for Collider Physics_ , 1912.06794.
  * (59) E. Buhmann, S. Diefenbacher, E. Eren, F. Gaede, G. Kasieczka, A. Korol et al., _Getting High: High Fidelity Simulation of High Granularity Calorimeters with High Speed_ , 2005.05334.
  * (60) S. Diefenbacher, E. Eren, G. Kasieczka, A. Korol, B. Nachman and D. Shih, _DCTRGAN: Improving the Precision of Generative Models with Reweighting_ , _Journal of Instrumentation_ 15 (2020) P11004, [2009.03796].
  * (61) A. Maevskiy, F. Ratnikov, A. Zinchenko and V. Riabov, _Simulating the Time Projection Chamber responses at the MPD detector using Generative Adversarial Networks_ , 2012.04595.
  * (62) K. Deja, J. Dubiński, P. Nowak, S. Wenzel and T. Trzciński, _End-to-end sinkhorn autoencoder with noise generator_ , 2006.06704.
  * (63) K. Dohi, _Variational Autoencoders for Jet Simulation_ , 2009.04842.
  * (64) F. Rehm, S. Vallecorsa, V. Saletore, H. Pabst, A. Chaibi, V. Codreanu et al., _Reduced Precision Strategies for Deep Learning: A High Energy Physics Generative Adversarial Network Use Case_ , 2103.10142.
  * (65) F. Rehm, S. Vallecorsa, K. Borras and D. Krücker, _Validation of Deep Convolutional Generative Adversarial Networks for High Energy Physics Calorimeter Simulations_ , 3, 2021, 2103.13698.
  * (66) F. Rehm, S. Vallecorsa, K. Borras and D. Krücker, _Physics Validation of Novel Convolutional 2D Architectures for Speeding Up High Energy Physics Simulations_ , 5, 2021, 2105.08960.
  * (67) G. R. Khattak, S. Vallecorsa, F. Carminati and G. M. Khan, _Fast Simulation of a High Granularity Calorimeter by Generative Adversarial Networks_ , 2109.07388.
  * (68) L. Anderlini, _Machine Learning for the LHCb Simulation_ , 10, 2021, 2110.07925.
  * (69) C. Fanelli and J. Pomponi, _DeepRICH: Learning Deeply Cherenkov Detectors_ , _Mach. Learn. Sci. Tech._ 1 (11, 2019) 015010, [1911.11717].
  * (70) Y. Lu, J. Collado, D. Whiteson and P. Baldi, _SARM: Sparse Autoregressive Model for Scalable Generation of Sparse Images in Particle Physics_ , 2009.14017.
  * (71) E. Buhmann, S. Diefenbacher, E. Eren, F. Gaede, D. Hundhausen, G. Kasieczka et al., _Hadrons, Better, Faster, Stronger_ , 2112.09709.
  * (72) C. Krause and D. Shih, _CaloFlow: Fast and Accurate Generation of Calorimeter Showers with Normalizing Flows_ , 2106.05285.
  * (73) A. Hariri, D. Dyachkova and S. Gleyzer, _Graph Generative Models for Fast Detector Simulations in High Energy Physics_ , 2104.01725.
  * (74) E. Buhmann, S. Diefenbacher, E. Eren, F. Gaede, G. Kasieczka, A. Korol et al., _Decoding Photons: Physics in the Latent Space of a BIB-AE Generative Network_ , 2102.12491.
  * (75) ATLAS collaboration, G. Aad et al., _AtlFast3: the next generation of fast simulation in ATLAS_ , 2109.02551.
  * (76) C. Krause and D. Shih, _CaloFlow II: Even Faster and Still Accurate Generation of Calorimeter Showers with Normalizing Flows_ , 2110.11377.
  * (77) S. Bieringer, A. Butter, S. Diefenbacher, E. Eren, F. Gaede, D. Hundhausen et al., _Calomplification - The Power of Generative Calorimeter Models_ , 2202.07352.
  * (78) A. Butter et al., _Machine Learning and LHC Event Generation_ , in _2022 Snowmass Summer Study_ , 3, 2022, 2203.07460.
  * (79) A. Butter and T. Plehn, _Generative Networks for LHC events_ , 2008.08558.
  * (80) M. Feickert and B. Nachman, _A Living Review of Machine Learning for Particle Physics_ , 2102.02770.
  * (81) S. Forte, L. Garrido, J. I. Latorre and A. Piccione, _Neural network parametrization of deep inelastic structure functions_ , _JHEP_ 05 (2002) 062, [hep-ph/0204232].
  * (82) R. D. Ball et al., _The Path to Proton Structure at One-Percent Accuracy_ , 2109.02653.
  * (83) V. N. Gribov and L. N. Lipatov, _Deep inelastic e p scattering in perturbation theory_ , _Sov. J. Nucl. Phys._ 15 (1972) 438–450.
  * (84) Y. L. Dokshitzer, _Calculation of the Structure Functions for Deep Inelastic Scattering and e+ e- Annihilation by Perturbation Theory in Quantum Chromodynamics._ , _Sov. Phys. JETP_ 46 (1977) 641–653.
  * (85) G. Altarelli and G. Parisi, _Asymptotic Freedom in Parton Language_ , _Nucl. Phys. B_ 126 (1977) 298–318.
  * (86) J. Bai, F. Lu, K. Zhang et al., “Onnx: Open neural network exchange.”
  * (87) S. Gieseke, C. Rohr and A. Siodmok, _Colour reconnections in Herwig++_ , _Eur. Phys. J. C_ 72 (2012) 2225, [1206.0041].
  * (88) S. Gieseke, P. Kirchgaeßer, S. Plätzer and A. Siódmok, _Soft Gluon Evolution as Guiding Principle for Colour Reconnection_ , _Acta Phys. Polon. B_ 50 (2019) 1871–1879.
  * (89) J. Bellm, C. B. Duncan, S. Gieseke, M. Myska and A. Siódmok, _Spacetime colour reconnection in Herwig 7_ , _Eur. Phys. J. C_ 79 (2019) 1003, [1909.08850].
  * (90) M. Arratia et al., _Presenting Unbinned Differential Cross Section Results_ , 2109.13243.
  * (91) G. ’t Hooft, _A Planar Diagram Theory for Strong Interactions_ , _Nucl. Phys. B_ 72 (1974) 461.
  * (92) Particle Data Group collaboration, M. Tanabashi et al., _Review of Particle Physics_ , _Phys. Rev._ D98 (2018) 030001.
  * (93) F. Chollet, “Keras.” https://github.com/fchollet/keras, 2017.
  * (94) A. Y. H. Adnrew L Maas and A. Y. Ng, _Rectifier nonlinearities improve neural network acoustic models_ , in _in ICML Workshop on Deep Learning for Audio, Speech and Language Processing_ , 2013.
  * (95) D. Kingma and J. Ba, _Adam: A method for stochastic optimization_ , 1412.6980.
  * (96) L. N. Vaserstein, _Markov processes over denumerable products of spaces, describing large systems of automata_ , _Problemy Peredači Informacii_ 5 (1969) 64.
  * (97) L. V. Kantorovich, _Mathematical Methods of Organizing and Planning Production_ , _Management Science_ 6 (1939) 366.
  * (98) L. Lonnblad, _ThePEG, Pythia7, herwig++ and Ariadne_ , _Nucl. Instrum. Meth. A_ 559 (2006) 246–248.
  * (99) O. R. developers, “Onnx runtime.” https://onnxruntime.ai/, 2021.
  * (100) M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean et al., _Tensorflow: A system for large-scale machine learning._ , in _OSDI_ , vol. 16, pp. 265–283, 2016.
  * (101) DELPHI collaboration, P. Abreu et al., _Tuning and test of fragmentation models based on identified particles and precision event shape data_ , _Z. Phys. C_ 73 (1996) 11–60.
  * (102) A. Buckley, J. Butterworth, D. Grellscheid, H. Hoeth, L. Lonnblad, J. Monk et al., _Rivet user manual_ , _Comput. Phys. Commun._ 184 (2013) 2803–2819, [1003.0694].
  * (103) D. Reichelt, P. Richardson and A. Siodmok, _Improving the Simulation of Quark and Gluon Jets with Herwig 7_ , _Eur. Phys. J. C_ 77 (2017) 876, [1708.01491].
  * (104) E. Farhi, _A QCD Test for Jets_ , _Phys. Rev. Lett._ 39 (1977) 1587–1588.
  * (105) S. Brandt, C. Peyrou, R. Sosnowski and A. Wroblewski, _The Principal axis of jets. An Attempt to analyze high-energy collisions as two-body processes_ , _Phys. Lett._ 12 (1964) 57–61.
  * (106) R. Abbate, M. Fickinger, A. H. Hoang, V. Mateu and I. W. Stewart, _Thrust at $N^{3}LL$ with Power Corrections and a Precision Global Fit for $\alpha_{s}(mZ)$_, _Phys. Rev. D_ 83 (2011) 074021, [1006.3080].
  * (107) H1 collaboration, V. Andreev et al., _Measurement of lepton-jet correlation in deep-inelastic scattering with the H1 detector using machine learning for unfolding_ , 2108.12376.
  * (108) M. Vandegar, M. Kagan, A. Wehenkel and G. Louppe, _Neural Empirical Bayes: Source Distribution Estimation and its Applications to Simulation-Based Inference_ , in _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ (A. Banerjee and K. Fukumizu, eds.), vol. 130 of _Proceedings of Machine Learning Research_ , pp. 2107–2115, PMLR, 11, 2021, 2011.05836, https://proceedings.mlr.press/v130/vandegar21a.html.
  * (109) M. Bellagente, A. Butter, G. Kasieczka, T. Plehn and R. Winterhalder, _How to GAN away Detector Effects_ , 1912.00477.
  * (110) M. Bellagente, A. Butter, G. Kasieczka, T. Plehn, A. Rousselot and R. Winterhalder, _Invertible Networks or Partons to Detector and Back Again_ , 2006.06685.
  * (111) A. Andreassen, P. T. Komiske, E. M. Metodiev, B. Nachman and J. Thaler, _OmniFold: A Method to Simultaneously Unfold All Observables_ , _Phys. Rev. Lett._ 124 (2020) 182001, [1911.09107].
  * (112) T. Sjöstrand and M. Utheim, _Hadron interactions for arbitrary energies and species, with applications to cosmic rays_ , _Eur. Phys. J. C_ 82 (2022) 21, [2108.03481].
  * (113) Y. Feng, M.-Y. Chu, U. Seljak and P. McDonald, _FASTPM: a new scheme for fast simulations of dark matter and haloes_ , _Mon. Not. Roy. Astron. Soc._ 463 (Dec., 2016) 2273–2286, [1603.00476].
  * (114) C. Modi, F. Lanusse and U. Seljak, _FlowPM: Distributed TensorFlow implementation of the FastPM cosmological N-body solver_ , _Astronomy and Computing_ 37 (Oct., 2021) 100505, [2010.11847].
  * (115) B. Dai and U. Seljak, _Learning effective physical laws for generating cosmological hydrodynamics with Lagrangian Deep Learning_ , _PNAS_ 118 (10, 2020) e2020324118, [2010.02926].
  * (116) V. Böhm, Y. Feng, M. E. Lee and B. Dai, _MADLens, a python package for fast and differentiable non-Gaussian lensing simulations_ , _Astron. Comput._ 36 (2021) 100490, [2012.07266].
  * (117) C. Modi, F. Lanusse, U. Seljak, D. N. Spergel and L. Perreault-Levasseur, _CosmicRIM : Reconstructing Early Universe by Combining Differentiable Simulations with Recurrent Inference Machines_ , 2104.12864.
  * (118) P. Ilten, T. Menzo, A. Youssef and J. Zupan, _Modeling hadronization using machine learning_ , 2203.04983.