Hunting for QCD Instantons at the LHC
in Events with Large Rapidity Gaps

V. A. Khozea,b, V. V. Khozeb, D. L. Milneb and M. G. Ryskina,b

a Petersburg Nuclear Physics Institute, NRC Kurchatov Institute, Gatchina, St.
Petersburg, 188300, Russia

b IPPP, Department of Physics, Durham University, Durham, DH1 3LE, UK

###### Abstract

We outline a strategy of how to search for QCD instantons of invariant mass 20
– 60 GeV in diffractive events in low luminosity runs at the LHC. We show that
by imposing appropriate selection criteria on the final states, one can select
the kinematic regime where the instanton signal exceeds the background by a
factor of at least 8. In spite of the relatively strong cuts that we impose on
the total transverse energy and the number of charged tracks,
$\sum_{i}E_{T,i}>15$ GeV, $N_{\rm ch}>20$ measured within the $0<\eta<2$
interval and excluding events with high $p_{T}$ particles, the expected cross-
section is sufficiently large to study the instanton production in the events
with Large Rapidity Gaps at low luminosities, thus avoiding problems with
pile-up. The paper also includes an updated computation of instanton cross-
sections and other parameters relevant for the ongoing studies.

E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

Instantons are non-perturbative classical solutions of Euclidean equations of
motion in non-abelian gauge theories [1]. In the semi-classical limit,
instantons provide dominant contributions to the path integral and describe
quantum tunnelling between different vacuum sectors of the theory [2, 3, 4].
Over the years a number of phenomenological models and approaches were
developed where QCD instantons were either directly responsible for
generating, or at least contributed to many key aspects of non-perturbative
low-energy dynamics of strong interactions [5, 6, 7, 8, 9, 10]. These include
the role of instantons in the breaking of the $U(1)_{A}$ symmetry and the
spontaneous breakdown of the chiral symmetry, the formation of quark and gluon
condensates, topological susceptibility of the vacuum, and the non-
perturbative generation of the axion potential. Instanton-based
phenomenological models provide at least a qualitative description of the QCD
vacuum.111Their potential short-comings originate from attempting to describe
the strongly-coupled QCD vacuum dynamics using a semi-classical (weak-
coupling) approximation where the path integral is expanded around classical
instanton background fields. These predictions are in good agreement with
lattice calculations [11], see also [12, 13].

Our motivation here is different. We would like to search for QCD instanton
events directly in scattering experiments. This setup provides a more
controlled environment than QCD vacuum modelling, as the scattering event
kinematics will enable us to remain in a relatively weakly-coupled regime
where semiclassical instanton methods are reliable. In scattering processes
the instanton can be thought of as new non-perturbative multi-particle
vertices in a Feynman diagram, and for its computation at the microscopic
level we will use the formalism developed in the recent papers [14, 15]. QCD
instanton-generated processes are predicted to be produced with a large
scattering cross-section at small centre-of-mass partonic energies, although,
as shown in [15], discovering them at hadron colliders using conventional
high-$p_{T}$ trigger requirements, is a challenging task that requires
alternative dedicated search strategies. Developing such a strategy based on
the event selection with large rapidity gaps at low luminosities is the
purpose of this paper.

In the past, in the context of deep inelastic scattering, QCD instanton
contributions [16, 17, 18] were searched for (but not observed) [19, 20] at
the HERA collider.

In QCD, the instanton configuration consists of the gauge field,

$\displaystyle A_{\mu}^{a\,{\rm inst}}(x)$ $\displaystyle=$
$\displaystyle\frac{2\rho^{2}}{g}\frac{\bar{\eta}^{a}_{\mu\nu}(x-x_{0})_{\nu}}{(x-x_{0})^{2}((x-x_{0})^{2}+\rho^{2})}\,,$
(1)

along with the fermion components for light ($m_{f}<1/\rho$) fermions,

$\bar{q}_{Lf}=\psi^{(0)}(x)\,,\quad q_{Rf}=\psi^{(0)}(x)\,.$ (2)

The gauge field $A_{\mu}^{a\,{\rm inst}}$ is the Belavin-Polyakov-Schwartz-
Tyupkin (BPST) instanton solution [1] of the self-duality equations in the
singular gauge. Here $\rho$ is the instanton size and $x_{0}$ is the instanton
position. Constant group-theoretic coefficients $\bar{\eta}^{a}_{\mu\nu}$ are
the ’t Hooft eta symbols defined in [2]. The fermionic components $\psi^{(0)}$
are the corresponding normalised solutions of the Dirac equation
$\gamma^{\mu}D_{\mu}[A_{\mu}^{a\,{\rm inst}}]\psi^{(0)}\,=0$. These are the
fermion zero modes of the instanton. The instanton configuration is a local
minimum of the Euclidean action, and the action on the instanton is given by
$S_{I}=\frac{8\pi^{2}}{g^{2}}=\frac{2\pi}{\alpha_{s}}$.

The instanton configuration (1) has topological charge equal to one and thus,
due to the chiral anomaly, the instanton processes violate chirality. If the
instanton is produced by a two-gluon initial state, the final state of this
instanton-mediated process will have $N_{f}$ pairs of quarks and anti-quarks
with the same chirality,

$g+g\to n_{g}\times g+\sum^{N_{f}}_{f=1}(q_{Rf}+\bar{q}_{Lf})\ ,$ (3)

where $N_{f}$ is the number of light flavours relative to the inverse
instanton size, $m_{f}<1/\rho$. The instanton contribution to the amplitude
for this process comes from expanding the corresponding path integral in the
instanton field background. At leading order in the instanton perturbation
theory, the amplitude takes the form of an integral over the instanton
collective coordinates (see e.g. Refs. [14, 15] for more detail),

${\cal A}^{\rm\,L.O.}_{\,2\to\,n_{g}+2N_{f}}=\int
d^{4}x_{0}\int_{0}^{\infty}d\rho\,D(\rho)\,e^{-S_{I}}\,\prod_{i=1}^{n_{g}+2}A_{{\rm
LSZ}}^{\rm inst}(p_{i};\rho)\,\prod_{j=1}^{2N_{f}}\psi^{(0)}_{{\rm
LSZ}}(p_{j};\rho).$ (4)

The factors of $A_{{\rm LSZ}}^{\rm inst}(p_{i};\rho)$ and $\psi^{(0)}_{{\rm
LSZ}}(p_{j};\rho)$ are the standard insertions of the LSZ-reduced instanton
fields in the momentum representation and $D(\rho)$ is given by the known
expression for the instanton density [2].

From the point of view of Feynman graphs the leading order instanton amplitude
(4) reveals itself as a family of multi-particle vertices (with different
numbers of emitted gluons) integrated over the instanton position and size. It
describes the emission of a large number of gluons, $n_{g}\propto
E^{2}/\alpha_{s}$, together with a fixed number of quarks and anti-quarks, one
pair for each light flavour in accordance with (3). The semi-classical
suppression factor, $\exp(-S_{I})=\exp(-2\pi/\alpha_{s})$, will be partially
compensated by the growth with jet energy, $E$, of the high multiplicity
cross-section for the process (3). The fully factorised structure of the field
insertions on the right hand side of (4) implies that at leading order in
instanton perturbation theory there are no correlations between the momenta of
the external legs in the instanton amplitude. The momenta of individual
particles in the final state are mutually independent, apart from overall
momentum conservation.

Thus to discover the QCD instanton we have to observe in the final state a
multi-particle cluster or a fireball which contains in general a large number
of isotropically distributed gluon (mini)jets accompanied by $N_{f}$ pairs of
light quark jets generated by a subprocess such as (3).

It is quite challenging, however, to identify the instanton on top of the
underlying event. Recall that the instanton is not a particle and there will
be no peak in the invariant mass, $M_{\rm inst}$, distribution 222What we mean
by the instanton mass is the partonic energy $\sqrt{\hat{s}}$ of the initial
2-gluon state in the process (3). As we integrate over the Bjorken x variables
when computing hadronic cross-section, we sum over a broad range of instanton
masses. . The mean value of $M_{\rm inst}$ can at least in principle be
“measured” or reconstructed as the mass of the minjet system created by the
instanton fireball in each given event. Talking about the instanton we
actually mean a family of objects of different sizes, $\rho$ and different
orientation in the colour and the Lorentz spaces. The mean value of $M_{\rm
inst}$ depends on $\rho$, increasing when $\rho$ decreases. Since
experimentally it is impossible to measure the instanton size, $\rho$, below
we use the mass $M_{\rm inst}$ to characterise the properties of the instanton
production event.

The paper is organised as follows; in Section 2 we detail the main sources of
background and how their cross-section behaves with energy. Section 3 details
the physics of events with large rapidity gaps while Section 4 explains the
rationale behind the selection criteria employed in this paper. Section 5
updates the calculation of cross-sections presented in [15] to account for
small virtualities present in incoming gluons and provides details of the
Monte Carlo generation of our background. In Section 6 we give our results and
Section 7 provides a discussion theoretical uncertainties and additional
instanton sub-processes. In section 8 we present our conclusions.

##  2  Background

Let us consider the main sources of background which can mimic the instanton
signal.

The cross-section of instanton production falls steeply with $M_{\rm inst}$
mainly due to the factor $\exp(-S_{I})=\exp(-2\pi/\alpha_{s}(\rho))$ in the
amplitude. At smaller values of $\rho$ the instanton action
$S_{I}=2\pi/\alpha_{s}(\rho)$ increases since the QCD coupling $\alpha_{s}$
decreases. Calculations presented in Section 5 (cf. Table 1) show that the
elementary cross-section of the parton level subprocess (3) falls
approximately as

$\hat{\sigma}_{\rm inst}\propto M^{-6}_{\rm inst}\,,$ (5)

over a broad low-to-intermediate energy range (it becomes less steep,
$\hat{\sigma}_{\rm inst}\propto M^{-4}_{\rm inst}$, at lower energies 20 – 30
GeV).

On the other hand the dimensionless rate, $M^{2}\,\sigma_{\rm pQCD}$, of a
similar purely perturbative QCD subprocess, $gg\to N{\rm minijets}$, decreases
only logarithmically. For the perturbatively formed ‘hedgehog’ configuration
of $N$ final state jets we would expect

$\sigma_{\rm pQCD}(gg\to N\,{\rm
jets})\,\sim\,\frac{16\pi}{M^{2}}\left(\frac{N_{c}}{\pi}\,\alpha_{s}(M)\right)^{N}\,,$
(6)

where $M$ denotes the invariant energy of the perturbatively formed cluster of
minijets. Thus, at sufficiently large values of $M_{\rm inst}$ the instanton
signal (5) will become negligible relative to the purely perturbative QCD
production (6).

In the regime of interest of moderately small $M_{\rm inst}$ we face however
another problem that needs to be addressed. The instanton event can be
mimicked by the multiple parton interactions (MPI) illustrated in Fig. 1.
Indeed, the double (triple, …, $n$) parton scattering produces a few pairs of
jets (dijets) which would look like a fireball, thus obscuring the final state
signature of the genuine instanton signal.

Figure 1: Schematic diagram for the triple parton scattering process.

The cross-section of such a process can be evaluated as

$\frac{d\sigma}{dE^{2}_{1}...dE^{2}_{n}}\,\sim\,\left(\frac{d\sigma}{\sigma_{eff}dE^{2}_{1}}\
...\ \frac{d\sigma}{\sigma_{eff}dE^{2}_{n}}\right)\sigma_{eff}\ ,$ (7)

where $E_{i}$ denotes the transverse energy of a jet in the $i$ dijet system,
$\sigma_{eff}\sim 10$ mb (see Fig. 4 of [21] and references therein) and the
cross section $d\sigma/dE^{2}_{i}\sim\pi\alpha_{s}^{2}/E^{4}_{i}$. That is, at
low $M_{\rm inst}$ when the transverse energies of produced jets become small,
the MPI processes will dominate. It was shown in [22] that the events caused
by such multiple parton interactions manifest high sphericity S, similar to
the instanton signal. Our strategy to suppress them is outlined in the
following section and involves a final state selection with Large Rapidity
Gaps.

Finally, we would like to emphasise one additional point. Since our knowledge
on the dynamics of confinement and hadronization is quite limited, general
purpose Monte Carlo (MC) event generators introduce certain phenomenological
parameters that are tuned in order to better describe the experimental data
(for a review see e.g. Ref. [23]). But if a low-mass instanton contributes to
the scattering, its presence will affect the data, and after such tuning we
lose the possibility to distinguish the low $M_{\rm inst}$ instanton from
other hadronization effects.

##  3 Events with Large Rapidity Gaps

As discussed above, we expect a large ‘underlying event’ background both in
the region of large and small instanton masses. However the low $M_{\rm inst}$
background caused by multi-parton interactions can be effectively suppressed
by selecting events with Large Rapidity Gaps (LRG). Indeed, each additional
’parton-parton$\to$ dijet’ scattering is accompanied by the colour flow
created by the parton cascade needed to form the incoming partons (see Fig.
1). This colour flow produces secondaries which fill the LRG. The LRG survival
probability, $S^{2}$, (i.e. the probability not to destroy the LRG) is rather
small, $S^{2}\leq 0.1$, see e.g. [24]. Thus the probability to observe $n$
additional branches of parton-parton interactions in LRG events is suppressed
by the factor $(S^{2})^{n}$.

Moreover, recall that the events with an LRG mainly occur at large values of
$b_{t}$, (the separation between the two incoming protons in the transverse
plane), where the proton optical density (opacity) is small [24, 25]. On the
other hand, the MPI (Double/Triple/…/n Parton Scattering) processes proceed
dominantly at small values of $b_{t}$, where the proton opacity is much
larger. Thus, we expect that the suppression factor, $S^{2}$, for each
additional scattering should be smaller than the average value of $S^{2}$
discussed in [24]. Therefore in this paper when evaluating backgrounds for
relatively low mass instantons in the events with an LRG, to a first
approximation we neglect the MPI/nPS contributions.

##  4 Search strategy

We propose searching for the instanton as a multi-particle cluster/fireball
with a mass of about 20 – 60 GeV in the events with an LRG. The presence of an
LRG can be detected either by detecting the leading forward proton with beam
momentum fraction, $x_{L}=1-\xi$, very close to 1 ($\xi=x_{Pom}\leq 0.01$), or
by observing no hadron activity in the forward calorimeters.

Since the relatively heavy instanton produces a rather large number of jets we
are looking for high multiplicity events which:

  * •

do not contain high-$E_{T}$ jets, and

  * •

still have a large density of the transverse energy,
$\sum_{i}dE_{Ti}/d\eta\sim M_{\rm inst}/3$ (the sum is over all secondary
particles in the given $\eta$ interval).

Indeed, from Table 1 we see that the instanton of mass 30 GeV produces about
17 jets (9 gluons plus 4 light $\bar{q}q$ pairs). The energy of each jet
$E_{Ti}\sim 1/\rho\sim 2$ GeV. After hadronization in such an event we expect
about 40-60 particles. The large multiplicity can be used as the main (or
additional) trigger to select the events of interest.

Note that the elementary (parton level) instanton cross-section is predicted
to be rather large ($\sim 1\mu b$) while the background due to hedgehog-like
pQCD multijet production (6) is much smaller – already for the case of 6 final
jets the expected elementary (hedgehog QCD) cross-section
$\sim\mbox{few}\cdot$nb which can be safely neglected.

On the other hand, in the events with an LRG, the probability to observe such
large multiplicity in a limited rapidity interval falls steeply with
$N_{track}$ (see e.g. Fig. 13 of [25]). Moreover, since each jet from the
instanton cluster contains a leading hadron with $p_{T}\geq 1$ GeV, we can
select the events with a large multiplicity (say,
$N_{ch}(p_{T}>0.5\mbox{GeV})>20$) of high $p_{T}$ particles, in this way
strongly suppressing the soft QCD background.

Figure 2: Instanton production in a diffractive process with an LRG. The
pomeron exchange is shown by the thick doubled line. The red bar shows the
range of $\eta$ considered in this paper. $Y$ indicates the incoming proton
position in rapidity. As shown in the diagram secondaries will be produced
outside this range but will not be used when calculating $E_{T}$ or $N_{ch}$.

Since the LHC detectors never cover the whole ($4\pi$) rapidity interval,
there is no chance to adequately measure the value of $M_{\rm inst}$. To
select the events with appropriate $M_{\rm inst}$ we introduce the cut on the
total transverse energy measured within the given rapidity interval
$\sum_{i}E_{Ti}>M_{0}$.

Here we consider the instanton production in the pomeron-proton collision (in
terms of Regge theory the pomeron exchange is responsible for the presence of
the LRG) selecting events with a large multiplicity and relatively large
transverse energy 333The idea to search for the instanton in events with very
large multiplicity but not too large transverse energy and the first
evaluation of the instanton cross-section at collider energies was discussed
long ago in [26].. We expect that these events will be more or less
spherically symmetric, that is, in such events there should be a large
probability to observe the sphericity $S$ close to 1. Since we can not observe
the particles in the whole $4\pi$ sphere we consider the “transverse
sphericity” (or cylindricity) defined as
$S_{T}=2\lambda_{2}/(\lambda_{1}+\lambda_{2})$ where $\lambda_{i}$ are the
eigenvalues of the matrix

$S^{\alpha\beta}=\frac{\sum_{i}p^{\alpha}_{i}p^{\beta}_{i}}{\sum_{i}|\vec{p}^{2}_{i}|}\
,$ (8)

and $\lambda_{2}<\lambda_{1}$. Here $p_{i}^{\alpha}$ is the two-dimensional
transverse component of the momentum of the $i$-th particle and we sum over
all particles observed in the event within a given rapidity interval.

Finally recall that the decay of the instanton produces one additional pair of
each flavour of light ($m_{f}<1/\rho$) quark. So in the case of the signal we
expect to observe a larger number of strange and charm particles than in
background events.

##  5 Calculation

###  5.1 Parton-level instanton cross-sections

The instanton cross-sections of the elementary gluon-initiated process (3) are
calculated following the approach developed in the recent papers [14, 15]. In
this approach the total parton-level instanton cross-section for the process
$gg\to X$ is related by the optical theorem to the imaginary part of the
forward elastic scattering amplitude computed in the background of the
instanton–anti-instanton ($I\bar{I}$) configuration,

$\displaystyle\hat{\sigma}_{\rm inst}$ $\displaystyle=$
$\displaystyle\frac{1}{E^{2}}\,{\rm Im}\,{\cal
A}^{{I\bar{I}}}_{4}(p_{1},p_{2},-p_{1},-p_{2})\,.$ (9)

Here $E=\sqrt{(p_{1}+p_{2})^{2}}$ denotes the invariant mass $\sqrt{\hat{s}}$
of the 2-gluon initial state. The resulting instanton cross-section takes the
form of the finite-dimensional integral [15] over the instanton and anti-
instanton collective coordinates (scale-sizes, $\rho$ and $\bar{\rho}$, the
separation between the instanton and anti-instanton centres, $R$, and the
relative $I\bar{I}$ orientations, $\Omega$):

$\displaystyle\hat{\sigma}_{\rm tot}^{\rm inst}\,\simeq\,\frac{1}{E^{2}}\,{\rm
Im}\,\frac{\kappa^{2}\pi^{4}}{36\cdot
4}\int\frac{d\rho}{\rho^{5}}\int\frac{d\bar{\rho}}{\bar{\rho}^{5}}\int
d^{4}R\int
d\Omega\left(\frac{2\pi}{\alpha_{s}(\mu_{r})}\right)^{14}(\rho^{2}E)^{2}(\bar{\rho}^{2}E)^{2}\,{\cal
K}_{\rm ferm}(z)$
$\displaystyle(\rho\bar{\rho}\mu_{r}^{2})^{b_{0}}\,\exp\left(R_{0}E\,-\,\frac{4\pi}{\alpha_{s}(\mu_{r})}\,{\cal
S}(z)\,-\,\frac{\alpha_{s}(\mu_{r})}{16\pi}(\rho^{2}+\bar{\rho}^{2})\,E^{2}\,\log\frac{E^{2}}{\mu_{r}^{2}}\,-\,Q(\rho+\bar{\rho})\right).$
(10)

In the above integral $\kappa^{2}$ denotes the known normalization constant
and ${\cal K}_{\rm ferm}(z)$ is the factor arising from the overlap of fermion
zero modes of $I$ and $\bar{I}$ and the variable $z$ is a certain conformally-
invariant ratio of the collective coordinates $R$, $\rho$ and $\bar{\rho}$
[27, 28],

$z\,=\,\frac{R^{2}+\rho^{2}+\bar{\rho}^{2}+\sqrt{(R^{2}+\rho^{2}+\bar{\rho}^{2})^{2}-4\rho^{2}\bar{\rho}^{2}}}{2\rho\bar{\rho}}\,.$
(11)

The first factor in the exponent in (10) describes the energy input from the
initial state, the second factor is the action of the instanton-anti-instanton
configuration [28], $S_{I\bar{I}}(z)=\frac{4\pi}{\alpha_{s}}\,{\cal S}(z)$
where,

$\displaystyle{\cal
S}(z)\,=\,3\frac{6z^{2}-14}{(z-1/z)^{2}}\,-\,17\,-\,3\log(z)\left(\frac{(z-5/z)(z+1/z)^{2}}{(z-1/z)^{3}}-1\right)\,.$
(12)

The third term in the exponent,

$\exp\left(-\,\frac{\alpha_{s}(\mu_{r})}{16\pi}\,(\rho^{2}+\bar{\rho}^{2})\,E^{2}\,\log\frac{E^{2}}{\mu_{r}^{2}}\right)\,,$
(13)

is the quantum effect that takes into account resummed perturbative exchanges
between the hard initial state gluons and the instantons [29]. Inclusion of
these quantum effects is required in order to resolve the well-known infra-red
problem in the $\rho,\bar{\rho}\to\infty$ limit, as QCD instantons with
$\rho\gg\left(\frac{16\pi}{\alpha_{s}}\frac{1}{E^{2}\log E^{2}}\right)^{1/2}$
are automatically cut-off by these quantum corrections.

We note that the initial-state partons (in our case gluons) are not strictly
on mass-shell, but carry small order-GeV virtualities, $Q_{1}$ and $Q_{2}$.
For the gluon emitted from the pomeron we choose $Q_{1}=2$ GeV, and the second
gluon that originated from the proton has the virtuality $Q_{2}=1$ GeV.
444From the pomeron side the expected virtuality $Q^{2}_{1}\simeq q^{2}_{t}$
is close to the mean transverse momentum squared, $q^{2}_{t}$, of the gluon
inside the pomeron. Its value can be estimated as $q^{2}_{t}\sim
1/\alpha^{\prime}_{P}\sim 4$ GeV2 where we take the pomeron trajectory slope
$\alpha^{\prime}_{p}=0.25$ GeV-2 [30]. The gluon virtuality in the proton PDF,
$Q_{2}\sim 1$ GeV was chosen for the following reasons. As a rule the global
parton analysis based on DGLAP evolution starts the evolution at some
$Q=Q_{0}\geq 1$ GeV. We do not know the parton distributions at a smaller $Q$
(it is not even clear whether we can use the parton language at such small
$Q$). Note however that the low-$x$ gluon density increases as the scale $Q$
increases. So we consider the value of $Q_{2}=1$ GeV as a conservative
estimate. These virtualities introduce the form-factor555The form-factor is a
direct consequence of Fourier transforming the instanton field to momentum
space to obtain $A^{\rm inst}_{LSZ}(p_{i})$, where the momentum $p_{i}$ has
the virtuality, $Q_{i}^{2}$ [28], and was explored and used extensively in the
context of deep inelastic scattering in [16, 18]. $e^{-Q\rho}$ in the
instanton vertex, with $Q=Q_{1}+Q_{2}=3$ GeV in our case. This is the origin
of the final term,

$\exp\left(-Q(\rho+\bar{\rho})\right)\,,$ (14)

in the exponent of the cross-section in (10).

To further simplify the integrand we select a natural value for the
renormalization scale dictated by the inverse instanton size. This
prescription removes the large
$(\rho\mu_{r})^{b_{0}}(\bar{\rho}\mu_{r})^{b_{0}}$ factor in front of the
exponent in (10). Hence we choose,

$\mu_{r}\,=\,1/\langle\rho\rangle\,=\,1/\sqrt{\rho\bar{\rho}}\,.$ (15)

For the reference point of $\alpha_{s}$ we choose its value at the $\tau$
mass, $\alpha_{s}(m_{\tau})=0.32$ following [31, 32] that fits the
experimental data [33, 34]. The choice of $m_{\tau}$ as the reference scale
for $\alpha_{s}$ is conveniently close to the interesting (for us) regime of
relatively light instantons with sizes $1/\langle\rho\rangle\sim$ few GeV. For
the running coupling at energy scales above $m_{\tau}$, we use the 1-loop
expression, while for the lower scales we freeze the coupling at a fixed
critical value $\alpha_{s}=0.35$,

$\frac{4\pi}{\alpha_{s}(\langle\rho\rangle)}\,\simeq\,\begin{cases}\quad\frac{4\pi}{0.32}\,-\,2b_{0}\log\left(\langle\rho\rangle\,m_{\tau}\right)&:\,\,{\rm
for}\,\,\langle\rho\rangle^{-1}\geq 1.45\,{\rm GeV}\\\
\quad\frac{4\pi}{0.35}&:\,\,{\rm for}\,\,\langle\rho\rangle^{-1}<1.45\,{\rm
GeV}\,.\end{cases}$ (16)

At these energy-scales we are in the regime of

$N_{f}=4\,,$ (17)

active quarks, and this is the $N_{f}$ value we use in $b_{0}=11-2N_{f}/3$ and
in the instanton density expressions for $\kappa^{2}$ and ${\cal K}_{\rm
ferm}$ in (10).

The rationale for freezing the coupling at 0.35 is that the perturbative
evolution of strong coupling at the order $\sim 1$ GeV scale is known to be in
conflict with observations. The strong coupling values extracted from the fits
to charmonium spectrum give $\alpha_{s}\simeq 0.35$ at 1 to 1.2 GeV scale
[35], which is close to the measured value $\alpha_{s}(m_{\tau})$ at the
considerably higher scale $m_{\tau}=1.777$ GeV, and is about 50% below the
prediction of the 2-loop perturbative running. These considerations justify
freezing the strong coupling in the infrared at a critical value
$\alpha_{s}=0.35$, the approach we will follow in accordance with (16).

E [GeV] |  $\hat{\sigma}\left(gg\rightarrow I\right)$ [pb] |  $\langle 1/\rho\rangle$ [GeV] | $\alpha_{s}\left(1/\rho\right)$ | $\langle n_{g}\rangle$
---|---|---|---|---
20 | $2.01\times 10^{6}$ | 1.69 | 0.327 | 7.81
25 | $9.49\times 10^{5}$ | 1.98 | 0.306 | 8.58
30 | $4.64\times 10^{5}$ | 2.27 | 0.290 | 9.07
35 | $2.32\times 10^{5}$ | 2.52 | 0.279 | 9.61
40 | $1.25\times 10^{5}$ | 2.84 | 0.267 | 9.67
50 | $3.89\times 10^{4}$ | 3.38 | 0.251 | 10.56
60 | $1.38\times 10^{4}$ | 3.87 | 0.241 | 10.89
70 | $5.45\times 10^{3}$ | 4.33 | 0.232 | 11.38
80 | $2.36\times 10^{3}$ | 4.85 | 0.224 | 11.67
90 | $1.08\times 10^{3}$ | 5.24 | 0.219 | 12.31
100 | $5.44\times 10^{2}$ | 5.82 | 0.213 | 12.10
110 | $2.79\times 10^{2}$ | 6.21 | 0.209 | 12.62
120 | $1.53\times 10^{2}$ | 6.71 | 0.205 | 12.77
130 | $8.56\times 10^{1}$ | 7.13 | 0.201 | 13.04
140 | $4.99\times 10^{1}$ | 7.57 | 0.198 | 13.25
150 | $3.01\times 10^{1}$ | 8.00 | 0.195 | 13.45

Table 1: The instanton production cross-section for the sub-process (3) for a
range of partonic cms energies $E$ (i.e. instanton masses) and the mean number
of emitted gluons in addition to the mean value of $1/\rho$ and
$\alpha_{s}\left(1/\rho\right)$. For energies larger than 150 GeV, the
suppression coming from the 3 GeV virtuality is fairly negligible and one can
use the numbers given in [15]. Note that the mean value of $R/\rho$ was
$\sim$1.55 at all energies.

The authors of Ref. [15], evaluated the integral in (10) numerically using the
python package SciPy [36] to derive the parton-level instanton cross-sections
along with the mean number of gluons in the final state as functions of
partonic energy $E$. In this work we follow the computational approach of
[15], except that we now also account for the effect of $Q\simeq 3$ GeV gluon
virtualities (14) from the outset, and use the prescription (16) for the
running coupling.

The data in Table 1 shows that the instanton production cross-section
$\hat{\sigma}_{\rm inst}$ is sharply suppressed at high partonic energies
(high instanton masses) and becomes large at low values of $E$. This is also
consistent with the earlier calculations in [14] that used the same theory
set-up for the parton-level instanton cross-section (10), but evaluated the
integral in the steepest descent approximation. The sharp suppression of the
instanton cross-section at large partonic energies is the consequence of the
inclusion of the quantum exchanges between the hard initial gluons that
resulted in the exponential form-factor (13). It cuts off the (anti)-instanton
sizes more and more efficiently when $E$ increases, ultimately leading to the
integral being dominated by smaller and smaller instantons, thus exponentially
suppressing the overall cross-section, in agreement with the results shown in
Table. 1 and Fig. 6 in Ref. [14]. In the low-energy limit, contributions of
large instantons with $\rho\gg 1/(3{\rm GeV})$ are cut-off by gluon
virtualities in (14).

For more detail on the theory formalism leading to (10) and the resulting
evaluation of $\hat{\sigma}_{\rm inst}(E)$ and $\langle n_{g}\rangle(E)$ we
refer the reader to Ref. [15].

###  5.2 Hadronic cross-sections for signal and background

The cross-section of instanton production in proton-proton collisions reads,

$\sigma_{pp\rightarrow I}\,=\,\int_{\hat{s}_{min}}^{s_{pp}}dx_{1}dx_{2}\quad
f\left(x_{1},Q_{1}^{2}\right)f\left(x_{2},Q_{2}^{2}\right)\hat{\sigma}_{\rm
inst}\left(E^{2}=x_{1}x_{2}s_{pp}\right)\,,$ (18)

where ${s}_{pp}$ is the centre-of-mass energy of the hadron collider and
$\hat{\sigma}$ is the partonic instanton cross-section. Selecting the events
with an LRG or the leading proton we have to replace one parton distribution
(PDF), say $f(x_{1},Q_{1}^{2})$, which describes the probability to find an
appropriate parton (gluon) in the incoming proton, by the so-called
diffractive PDF, $f^{D}(x_{1},Q_{1}^{2})$, which describes the corresponding
distribution in the pomeron.

The lower limit $\hat{s}_{min}$ in the integral (18) is introduced for
technical reasons. Strictly speaking, one is supposed to integrate from
$\hat{s}_{min}=0$ to ${s}_{pp}$. However instantons of masses much smaller
than 20 GeV are not under theoretical control, the semiclassical dilute
instanton approximation becomes invalid for large instantons with $\rho\sim$
GeV, where the theory can no longer be considered weakly-
coupled.666Furthermore, when $\rho\gtrsim 0.3\,{\rm fm}\sim 1.5\,{\rm
GeV}^{-1}$ one can no longer distinguish between the instantons created in
collisions and the instanton-like configurations populating the QCD ground
state [11]. The role of $\hat{s}_{min}$ is to cut-off such large ground-state
field configurations. It is impossible to fix $\hat{s}_{min}$ experimentally,
instead we have to choose and impose experimental cuts on final states that
will suppress the contribution of the low mass instantons. In practice we will
set $\sqrt{\hat{s}_{min}}=20\,\text{GeV}$ in the integral,777Taking
$\hat{s}_{min}=\left(20\text{GeV}\right)^{2}$ gives a hadronic instanton
cross-section of 1.72 $\mu$b (after multiplying by $S^{2}$ and taking the
factorisation scale as $1/\langle\rho\rangle$). and then impose the cuts.

Instanton samples are generated using the RAMBO algorithm [37] and then
showered and hadronised using PYTHIA 8 [38]. The scale of the process was
taken to be $1/\langle\rho\rangle$ for the purposes of showering in PYTHIA. As
already mentioned, in the integral (18) we set the minimum instanton mass
generated to be $\hat{s}_{min}=20$ GeV. However we also note that the
inclusion of these low mass instantons is expected to give a negligible
contribution with our choice of selection criteria.

Background events were simulated by generating inelastic pion-proton
collisions (see Fig. 2) in PYTHIA for the proton energy $E_{p}=7$ TeV and the
pomeron energy $E_{Pom}=7\cdot x_{Pom}$ TeV (with a pion being used to
simulate a pomeron) and $x_{Pom}$ fixed at 0.003. For the proton PDF we use
the NNPDF3.1luxQED NNLO [39, 40, 41, 42] set while for the parton distribution
in the pomeron (mimicked in PYTHIA by the pion) the H1 fit B [43] is used. The
$Q^{2}$ evolution for the pomeron PDF is performed by QCDNUM [44].

When computing the signal, the instanton cross-section calculated in (18) was
multipled by the gap survival probabilty, $S^{2}$. For simplicity we took
$S^{2}=0.1$ in agreement with the experimental results [45, 46] for
diffractive high $E_{T}$ dijet production. The background cross-section was
normalised to the experimental yield of single dissociation
$d\sigma^{SD}/d\ln(1/x_{Pom})\simeq 0.15$ mb [47]. In generating both the
instanton and the background samples, multi-parton interactions were disabled
in PYTHIA.

Since in the case of soft QCD background only a very small fraction of events
satisfy the $\sum_{i}E_{Ti}>M_{0}$ GeV cut888Actually to get a large
$\sum_{i}E_{Ti}$ in a perturbative QCD event we have to produce a high $E_{T}$
jet or to consider the event with a large number of multiple parton
interactions. The last possibility is strongly suppressed by the LRG survival
factor. (which we had used further) we generate the ”Hard QCD” events with a
minimum transverse momentum of 10 GeV using PYTHIA. Using lower $p_{T}$
samples from PYTHIA, we concluded that the contribution of lower $p_{T}$
events is negligible.

##  6 Results

Figure 3: Multiplicity distribution of charged hadrons produced in the events
with the instanton (green) in comparison with the expected background (red).
The number of events is normalised to the integrated luminosity
$L=1\,\mbox{pb}^{-1}$ and $\Delta\ln(x_{Pom})=1$ interval. We required events
to have $\sum_{i}E_{T,i}>15$ GeV and $N_{ch}>20$, summing only over charged
particles in the region $0<\eta<2$ with $p_{T}>0.5$ GeV, with an additional
constraint that there is no charged particle in this region with $p_{T}>2$ GeV
(left figure), or no charged particle in the region $-2<\eta<2$ with $p_{T}>$
2.5 GeV (right figure).

As an example we consider instanton production at 14 TeV in the events with
the leading proton momentum fraction $x_{L}=1-x_{Pom}=0.997$. That is we
calculate the cross-section $d\sigma/d\ln(1/x_{Pom})$ at $x_{Pom}=0.003$. This
may be considered as the total observed cross-section when/if one integrates
over the interval of $x_{Pom}$ from 0.0018 to 0.005 (i.e. $\Delta\ln(x_{pom})$
interval equal to 1).999Note that in this region the dependence of single
dissociation cross-section $d\sigma/d\ln(1/x_{Pom})$ is practically flat as is
shown in [25, 47, 48].

Our results are presented in Figs. 3 and 4. To further suppress the background
we impose the following set of experimental selection criteria that look
realistic for the present ATLAS and CMS detectors. We exclude the very low
$p_{T}$ particles and consider only the secondary hadrons with transverse
momentum $p_{T}>0.5$ GeV produced within the $0<\eta<2$ rapidity interval, as
shown in Fig. 2. We select the events for which the total transverse energy of
the charged secondaries (with $p_{T}>0.5$ GeV and $0<\eta<2$)101010We shifted
the $\eta$ interval in the direction opposite to the LRG (or the leading
proton which has the negative rapidity) in order not to have gluons with too
large $x$ (where the gluon density rapidly decreases with increasing $x$) from
the pomeron PDF. Recall that for $x_{Pom}=0.003$ the energy of the ‘incoming’
pomeron is 21 GeV only. is $\sum_{i}E_{T,i}>15$ GeV and the number of
corresponding charged tracks $N_{ch}>20$. We also demand that there are no
charged particles in this region with $p_{T}>2$ GeV [left Figs. 3 and 4].
Plots on the right hand side in Figs. 3-4 implement the constraint that no
charged particles are present in the region $-2<\eta<2$ with $p_{T}>$ 2.5 GeV.

As is seen in Fig. 4 after the proposed cuts the instanton signal exceeds the
background by a factor of over (left) $8$ (right) $4$ and for $S_{T}>0.5$ the
signal to background ratio S/B is even higher. Notably, in the $N_{ch}>20$
region the expected cross-section is still rather large ($\sim$ 1 nb). This
allows one to observe the instanton production in low luminosity runs without
pile-up problems.

Figure 4: Distribution over the transverse sphericity $S_{T}$, Eq. (8), of the
charged hadrons produced in the events with the instanton (green) in
comparison with the expected background (red). The selection criteria used are
the same as those described in Fig. 3

.

Figure 5: Distribution over the number of neutral kaons produced in the events
with the instanton (green) in comparison with the expected background (red).
The selection criteria used are the same as those described in Fig. 3.

The sphericity distribution is shown in Fig. 4. It is clearly seen that while
the background sphericity falls for $S_{T}>0.5$ the instanton signal is mainly
concentrated at $S_{T}\sim 0.7-0.8$. For $S_{T}>0.8$ the expected instanton
contribution is a factor of 10 greater than the background (for both sets of
selection criteria).

Finally, Fig. 5 shows the distribution over the number of neutral kaons.
111111This can be observed via the $(K^{0})_{s}\to\pi^{+}\pi^{-}$ decay. The
selection criteria used are the same as those described in Fig. 3. As expected
the kaon multiplicity is larger in the events containing the instanton and for
$N_{K^{0}}>5$ the S/B ratio is about 5 (for both sets of selection criteria).

##  7 Discussion

The possibility of direct experimental observation of QCD instantons at hadron
colliders has recently attracted a fair amount of attention [14, 15, 49]. The
authors of [15] evaluated the sensitivity of the LHC and Tevatron high- and
low-luminosity runs, and [49] obtained a first limit based on existing Minimum
Bias data.

Extraction of the instanton signal on top of the underlying event at the LHC,
not relying on events with large rapidity gaps, is a challenging task. One
problem, not considered in [49], is how to select in the experiment the events
corresponding to the chosen $M_{\rm inst}$ interval. It is impossible to
distinguish the jet emitted by the instanton from the jet created during the
parton showering or due to the presence of a perturbative QCD multi-jet
‘hedgehog’ contribution (the latter is another problem that has not been
addressed until now).

The idea to observe instantons in the pomeron-induced process, namely in the
central exclusive production,

$pp~{}\to~{}p~{}+~{}X~{}+~{}p\ $ (19)

was first put forward in [50] (see also [51]). This was prompted by the
reported UA8 observation [52] of the large enhancement in the pomeron-pomeron
cross-section at invariant masses $2<M<5$ GeV with roughly isotropic
distribution of secondaries. Unfortunately, we (present authors) do not have
theoretical control over such low-mass instantons, or for that matter other
non-perturbative contributions in the infrared. In these settings, there is no
external mass parameter and it is not clear how to distinguish instanton
production from the background or, for example, from the glueball production.
On the other hand searching for a larger-mass instanton in pomeron-pomeron
collisions has no advantage in comparison with instanton production in events
with one LRG, but in the case of pomeron-pomeron collisions the expected
cross-section is much smaller.

###  7.1 Other instanton production sub-processes

Up to now we were only concentrating on the 2-gluon initiated process (3), but
one should also consider other quark, anti-quark and gluon initiated 2-parton
processes, such as,

$\displaystyle g+u_{L}$ $\displaystyle\to$ $\displaystyle n_{g}\times
g\,+\,u_{R}+\sum_{f=1}^{N_{f}-1}(q_{Rf}+\bar{q}_{Lf})\,,$ (20) $\displaystyle
u_{L}+\bar{u}_{R}$ $\displaystyle\to$ $\displaystyle n_{g}\times
g+\sum_{f=1}^{N_{f}-1}(q_{Rf}+\bar{q}_{Lf})\,,$ (21) $\displaystyle
u_{L}+d_{L}$ $\displaystyle\to$ $\displaystyle n_{g}\times
g\,+\,u_{R}+d_{R}+\sum_{f=1}^{N_{f}-2}(q_{Rf}+\bar{q}_{Lf})\,,$ (22)

and also include the contributions from multi-parton initial states.

We first consider the 2-parton initiated processes. Note that there is no
interference between all these sub-processes since the fermion content of the
external states in (3),(20)-(22) is different. Thus, the total contribution
comes from adding the individual cross-sections from these processes, in other
words, summing over the integrals (18) with the appropriate choice of pdfs for
each of the partons in the initial state. Since each of these contributions is
manifestly positive, our already computed contribution from the gluon fusion
process (3) cannot be reduced. It is also clear that one cannot expect any
sufficient enhancement to it from accounting for the additional channels
(20)-(22). This is because for the parton-level cross-section most of the
expression in (10), in particular its exponent, remains the same;121212Up to
now we have only included the initial-state corrections arising from the high-
energy limit of the gluon propagator in the instanton background, but as was
shown in [53], the same result (13) is expected to hold also for fermions in
the initial state. only the four field insertions for the external states
vary. On the other hand, as we have seen, the dominant contributions to
hadronic cross-sections come from the low-$x$ region, and as a result, our
gluon fusion instanton sub-process is expected to be the dominant one.

Next, we would like to discuss the effects of multi-parton initial states in
the instanton production. To do this we would like to distinguish between the
two cases: one, where multi-partons entering the instanton have originated
from just the two parent partons, and the second case, where the multi-parton
initial state has a genuine multi-parton origin, with all 3 or 4 partons being
emitted from the proton and/or the pomeron.

The first case, where the multi-partons in the initial state can be traced to
just two parent partons is already accounted for by the resummed quantum
corrections from the initial states. Such effects are precisely what is
included into the initial-initial state and the initial-final state
propagators in the instanton background computed in [29]. Our approach takes
into account the initial-initial (i.e. hard-hard) quantum effects by
explicitly including the Mueller form-factor term (13). The initial-final
state (or hard-soft interactions) were not directly included in our
approach,131313This is done to avoid double-counting, since the radiative
exchanges involving final-state partons are already accounted for by the
optical theorem. however, the form-factor arising from all such quantum
corrections is of the same generic exponential form [29],

$\exp\Bigl{(}-\,\frac{\alpha_{s}(\mu_{r})}{16\pi}\,\rho^{2}E^{2}\,\sum_{i,j}a_{ij}\log
a_{ij}\Bigr{)},\quad{\rm where}\quad a_{ij}:=p_{i}\cdot p_{j}/E^{2}\,,$ (23)

as the initial-initial state form-factor in (13). Hence we do not expect a
significant difference in the final result from the inclusion of these
effects.

The second case corresponds to genuine multi-parton initial states, and it is
accounted for by adding the contribution to the total cross-section from the
integrals of the type (18), now with multiple pdfs, one for each parton in the
initial state. Like in the two-particle case, discussed earlier, this is a
manifestly positive contribution that cannot destructively interfere with our
two-gluon contribution. In summary, we expect that a more careful inclusion of
all instanton production sub-processes could only enhance the gluon fusion
contribution (and only moderately).

###  7.2 Theory uncertainties

Note that all the present calculations should not be considered as the precise
theoretical predictions. There are a number of uncertainties.

First of all, one could expect sizeable corrections to the assumed
factorization of instanton and anti-instanton collective coordinate
integration measure in the forward elastic scattering amplitude (10), which
strictly holds only at infinitely large inter-instanton separations
$R^{2}/(\rho\bar{\rho})\gg 1$. At finite separations, the integration measure
was shown in [27] to satisfy,

$d\mu_{I\bar{I}}\,=\,d\mu_{I}\,d\mu_{\bar{I}}\left(1+{\cal
O}\left(\frac{1}{z^{2}}\right)\right)\,,$ (24)

where $z$ is given by (11). In our calculation the mean value of $R/\rho$ was
$\simeq$1.55 at all energies, which corresponds to $z\simeq 4.16$, or
$1/z^{2}\simeq 0.058$ which is a reassuringly small value to be optimistic
about the validity of the approximate factorization of instanton integration
measures (and thus the factorization of the determinants of quadratic
fluctuation operators in the instanton–anti-instanton background).

Hence, in terms of the conformal $z^{2}$ variable the individual instanton and
the anti-instanton are reasonably well-separated. For this reason, we are
similarly confident about the validity of the instanton–anti-instanton action
expression (12) that was computed for a particular exactly-solvable model [27,
28] of the instanton-anti-instanton valley configuration. What is more
difficult to estimate, however, is the role played by the unknown higher-order
perturbative corrections to the re-summed quantum exchanges from the initial
states in (13).

Furthermore, there is a rather large theoretical uncertainty in evaluation of
the gap survival factor $S^{2}$, which also depends on the details of the
particular subprocess (see [24]). Therefore we prefer to use the experimental
value measured in diffractive dijet production[45, 46]. We evaluate the
corresponding uncertainty as $\pm 50$%.

We should also note a potentially strong dependence of the predicted cross-
section on the values of the renormalization, $\mu_{R}$, and factorization,
$\mu_{F}$, scales. In particular, it was found in [54] that replacing the
natural scale $\mu_{F}=1/\rho$ by the value of $\mu_{F}=M_{\rm inst}$ enhances
the instanton production cross-section for instantons of mass larger than
$M_{\rm inst}=50$ GeV computed in [14] by approximately a factor of five.

Besides this, as we have already noted, the present calculations are based on
the simplest $gg\to Instanton$ subprocess neglecting the possibility to
produce the instanton in collision with quarks or with a pair of gluons, like
$g+{gg}\to Instanton$ and so on. Based on our earlier discussion in Section
7.1, we consider it unlikely that these new subprocesses would change the
qualitative features of the expected signal, but they may enlarge the final
cross-section.

##  8 Conclusions

We propose to search for the QCD instantons at the LHC in events with Large
Rapidity Gaps and in this paper we carried out the first preliminary
investigation of this search strategy. These LRG events are selected either by
detecting the forward leading proton with $x_{L}$ very close to 1 or by
observing no hadron activity in the forward calorimeter. We discussed the main
sources of background for the low and the high mass instantons and showed that
background for the low mass instantons, which comes mainly from multiple
parton interactions, is effectively suppressed in the LRG events by the small
gap survival factor $S^{2}\leq 0.1$. This allowed us to find the kinematical
domain where the signal from $M_{\rm inst}\sim 20-60$ GeV instantons reliably
exceeds background.

Even with these rather strong cuts in place, the expected instanton cross-
section remains sufficiently large ($\sim$ 1 nb) to effectively produce and
probe QCD instantons at the LHC, at low luminosity runs, avoiding pile-up
problems.

## Acknowledgments

We thank Valery Schegelsky and Michael Spannowsky for useful discussions. VVK
is supported by the STFC under grant ST/P001246/1. DLM is supported by an STFC
studentship.

## References

  * [1] A. A. Belavin, A. M. Polyakov, A. S. Schwartz and Y. S. Tyupkin, Phys. Lett. 59B (1975) 85.
  * [2] G. ’t Hooft, Phys. Rev. D 14 (1976) 343, Erratum: [Phys. Rev. D 18 (1978) 2199].
  * [3] C. G. Callan, Jr., R. F. Dashen and D. J. Gross, Phys. Lett. 63B (1976) 334.
  * [4] R. Jackiw and C. Rebbi, Phys. Rev. Lett. 37 (1976) 172.
  * [5] G. ’t Hooft, Phys. Rept. 142 (1986) 357.
  * [6] C. G. Callan, Jr., R. F. Dashen and D. J. Gross, Phys. Rev. D 17 (1978) 2717.
  * [7] V. A. Novikov, M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B 191 (1981) 301.
  * [8] E. V. Shuryak, Nucl. Phys. B 203 (1982) 116.
  * [9] D. Diakonov and V. Y. Petrov, Phys. Lett. 147B (1984) 351; Nucl. Phys. B 272 (1986) 457.
  * [10] T. Schäfer and E. V. Shuryak, Rev. Mod. Phys. 70 (1998) 323 hep-ph/9610451.
  * [11] A. Hasenfratz and C. Nieter, Phys. Lett. B 439 (1998) 366 hep-lat/9806026.
  * [12] T. A. DeGrand, A. Hasenfratz and T. G. Kovacs, Nucl. Phys. B 505 (1997) 417 hep-lat/9705009.
  * [13] D. A. Smith et al. [UKQCD Collaboration], Phys. Rev. D 58 (1998) 014505 hep-lat/9801008.
  * [14] V. V. Khoze, F. Krauss and M. Schott, JHEP 2004 (2020) 201 arXiv:1911.09726 [hep-ph].
  * [15] V. V. Khoze, D. L. Milne and M. Spannowsky, Phys. Rev. D 103 (2021) no.1, 014017 arXiv:2010.02287 [hep-ph].
  * [16] I. I. Balitsky and V. M. Braun, Phys. Lett. B 314 (1993) 237 hep-ph/9305269.
  * [17] S. Moch, A. Ringwald and F. Schrempp, Nucl. Phys. B 507 (1997) 134 hep-ph/9609445.
  * [18] A. Ringwald and F. Schrempp, Phys. Lett. B 438 (1998) 217 hep-ph/9806528.
  * [19] C. Adloff et al. [H1 Collaboration], Eur. Phys. J. C 25 (2002) 495 hep-ex/0205078.
  * [20] S. Chekanov et al. [ZEUS Collaboration], Eur. Phys. J. C 34 (2004) 255 hep-ex/0312048.
  * [21] M. Aaboud et al. [ATLAS Collaboration], Phys. Lett. B 790 (2019) 595 arXiv:1811.11094 [hep-ex].
  * [22] M. Sas and J. Schoppink, arXiv:2101.12367 [hep-ph].
  * [23] A. Buckley, J. Butterworth, S. Gieseke, D. Grellscheid, S. Hoche, H. Hoeth, F. Krauss, L. Lonnblad, E. Nurse and P. Richardson, et al. Phys. Rept. 504 (2011), 145-233 arXiv:1101.2599 [hep-ph].
  * [24] V. A. Khoze, A. D. Martin and M. G. Ryskin, J. Phys. G 45 (2018) no.5, 053002 arXiv:1710.11505 [hep-ph].
  * [25] V. A. Khoze, A. D. Martin and M. G. Ryskin, Eur. Phys. J. C 81 (2021) no.2, 175 arXiv:2012.07967 [hep-ph].
  * [26] I. I. Balitsky and M. G. Ryskin, Phys. Atom. Nucl. 56, 1106 (1993) [Yad. Fiz. 56N8, 196 (1993)]; Phys. Lett. B 296 (1992) 185.
  * [27] A. V. Yung, Nucl. Phys. B 297 (1988) 47.
  * [28] V. V. Khoze and A. Ringwald, Phys. Lett. B 259 (1991) 106;
“Valley trajectories in gauge theories,” CERN-TH-6082-91 (1991).

  * [29] A. H. Mueller, Nucl. Phys. B 348 (1991) 310; Nucl. Phys. B 353 (1991) 44.
  * [30] A. Donnachie and P. V. Landshoff, Nucl. Phys. B 231 (1984) 189.
  * [31] M. Beneke and M. Jamin, JHEP 0809 (2008) 044 arXiv:0806.3156 [hep-ph].
  * [32] G. Abbas, B. Ananthanarayan, I. Caprini and J. Fischer, Phys. Rev. D 87 (2013) no.1, 014008 arXiv:1211.4316 [hep-ph].
  * [33] S. Schael et al. [ALEPH Collaboration], Phys. Rept. 421 (2005) 191 hep-ex/0506072.
  * [34] K. Ackerstaff et al. [OPAL Collaboration], Eur. Phys. J. C 7 (1999) 571 hep-ex/9808019.
  * [35] A. M. Badalian and V. L. Morgunov, Phys. Rev. D 60 (1999) 116008 hep-ph/9901430.
  * [36] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser and J. Bright, et al. Nature Meth. 17 (2020), 261 arXiv:1907.10121 [cs.MS].
  * [37] R. Kleiss, W. J. Stirling and S. D. Ellis, Comput. Phys. Commun. 40 (1986), 359
  * [38] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen and P. Z. Skands, Comput. Phys. Commun. 191 (2015), 159-177 arXiv:1410.3012 [hep-ph].
  * [39] A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr and G. Watt, Eur. Phys. J. C 75 (2015), 132 arXiv:1412.7420 [hep-ph].
  * [40] V. Bertone et al. [NNPDF], SciPost Phys. 5 (2018) no.1, 008 arXiv:1712.07053 [hep-ph].
  * [41] A. V. Manohar, P. Nason, G. P. Salam and G. Zanderighi, JHEP 12 (2017), 046 arXiv:1708.01256 [hep-ph].
  * [42] A. Manohar, P. Nason, G. P. Salam and G. Zanderighi, Phys. Rev. Lett. 117 (2016) no.24, 242002 arXiv:1607.04266 [hep-ph].
  * [43] A. Aktas et al. [H1 Collaboration], Eur. Phys. J. C 48 (2006) 715 hep-ex/0606004.
  * [44] M. Botje, Comput. Phys. Commun. 182 (2011) 490 arXiv:1005.1481 [hep-ph].
  * [45] A. M. Sirunyan et al. [CMS and TOTEM Collaborations], Eur. Phys. J. C 80 (2020) no.12, 1164 arXiv:2002.12146 [hep-ex].
  * [46] G. Aad et al. [ATLAS Collaboration], Phys. Lett. B 754 (2016) 214 arXiv:1511.00502 [hep-ex].
  * [47] G. Aad et al. [ATLAS Collaboration], JHEP 2002 (2020) 042 Erratum: [JHEP 2010 (2020) 182] arXiv:1911.00453 [hep-ex].
  * [48] V. Khachatryan et al. [CMS Collaboration], Phys. Rev. D 92 (2015) no.1, 012003 arXiv:1503.08689 [hep-ex].
  * [49] S. Amoroso, D. Kar and M. Schott, arXiv:2012.09120 [hep-ph].
  * [50] E. Shuryak and I. Zahed, Phys. Rev. D 68 (2003), 034001 arXiv:hep-ph/0302231 [hep-ph].
  * [51] E. Shuryak and I. Zahed, arXiv:2102.00256 [hep-ph].
  * [52] A. Brandt et al. [UA8], Eur. Phys. J. C 25 (2002), 361-377 arXiv:hep-ex/0205037 [hep-ex].
  * [53] V. V. Khoze and D. L. Milne, Int. J. Mod. Phys. A 36 (2021) no.05, 2150032 arXiv:2011.07167 [hep-ph].
  * [54] M. Mangano, arXiv:2101.02719 [hep-ph].


# Simplicial convolutional neural networks

###### Abstract

Graphs can model networked data by representing them as nodes and their
pairwise relationships as edges. Recently, signal processing and neural
networks have been extended to process and learn from data on graphs, with
achievements in tasks like graph signal reconstruction, graph or node
classifications, and link prediction. However, these methods are only suitable
for data defined on the nodes of a graph. In this paper, we propose a
simplicial convolutional neural network (SCNN) architecture to learn from data
defined on simplices, e.g., nodes, edges, triangles, etc. We study the SCNN
permutation and orientation equivariance, complexity, and spectral analysis.
Finally, we test the SCNN performance for imputing citations on a coauthorship
complex.

Index Terms—  Simplicial complex, Hodge Laplacian, simplicial filter,
simplicial neural network.

##  1 Introduction

Graphs are powerful models to represent irregular data by encoding their
pairwise relationships. To process such networked data, signal processing
concepts have been extended to the graph domain, defining, for instance, the
graph Fourier transform and graph filters [1]. Meanwhile, graph neural
networks (GNNs) have achieved a good performance in tasks like rating
prediction in recommender systems, graph or node classification and link
prediction [2, 3, 4, 5, 6].

However, graph signal processing and GNNs are designed for data defined on
nodes of a graph. In real-world problems, we might also have data defined on
edges, triangles, etc, of a network, such as communication or traffic flow in
a data or road network, paper citations of a coauthorship network and so on
[7, 8, 9]. To model data defined on such higher-order network structures, we
can use a simplicial complex, which is a collection of simplices, i.e., nodes,
edges, triangles, etc. Recently, signal processing and neural networks on
simplicial complexes have emerged. In [7], simplicial data as well as the
simplicial Fourier transform have been defined. The authors of [10, 11]
provided an overview of some simplicial signal processing techniques to
process flow-typed data. Our previous work [12] analyzed the definition of
frequencies for simplicial signals and proposed two types of simplicial
filters based on the Hodge Laplacian.

Meanwhile, researchers have also attempted to develop neural networks on
simplicial complexes. In [8], a basic simplicial neural network (SNN) was
proposed with a convolutional layer composed of a basic simplicial filter [12]
and a nonlinearity. Message passing neural networks (MPNNs) have been
generalized to simplicial complexes in [13] where the aggregation and updating
functions consider in addition to the edge data also data defined on adjacent
simplices, i.e., nodes and triangles. The neural network architectures in [14,
15] are instances of [13] by specifying the aggregation functions as
simplicial filters. Another attempt in [16] considered recurrent architectures
in MPNNs for flow interpolation and graph classification tasks.

Motivated by the principle of the convolution operator, in this paper we
propose simplicial convolutional neural networks (SCNNs). Differently from the
earlier approach in [8], we build an SCNN with simplicial filters of higher
flexibility in exploiting the lower- and upper-neighbors of a simplex. And
differently from the MPNN, the proposed SCNN considers multihop information
exchange within a layer and enjoys spectral interpretability via the
simplicial Fourier transform. Our specific contributions are:  i) we propose
an inductive SCNN based on the more advanced simplicial filter of [12] and
discuss its connections to the related work;  ii) we analyze its permutation
and orientation equivariances as well as characterize the SCNN in the spectral
domain;  iii) we test the performance on citation data imputation in a
coauthorship complex outperforming the state-of-the-art.

##  2 Simplicial signal processing

In this section, we recall some important simplicial signal processing
concepts, including simplicial complexes and signals, the Hodge decomposition,
and simplicial filters.

Simplicial complexes and signals. Given a finite set of vertices
$\mathcal{V}$, a $k$-simplex $\mathcal{S}^{k}$ is a subset of $\mathcal{V}$
with cardinality $k+1$. A _face_ of $\mathcal{S}^{k}$ is a subset with
cardinality $k$ and thus a $k$-simplex has $k+1$ faces. A _coface_ of
$\mathcal{S}^{k}$ is a $(k+1)$-simplex that includes $\mathcal{S}^{k}$ [7,
17]. A simplicial complex of order $K$, $\mathcal{X}^{K}$, is a collection of
$k$-simplices $\mathcal{S}^{k}$, $k=0,\dots,K$, with an inclusion property–for
any $\mathcal{S}^{k}\in\mathcal{X}^{K}$, then
$\mathcal{S}^{k-1}\in\mathcal{X}^{K}$ if
$\mathcal{S}^{k-1}\subset\mathcal{S}^{k}$. We denote the number of
$k$-simplices in $\mathcal{X}^{K}$ by $N_{k}$. If two simplices share a common
face, then they are lower neighbours; if they share a common coface, they are
upper neighbours [12]. A graph is a simplicial complex where nodes are
0-simplices, and edges are 1-simplices.

In a simplicial complex, we define a $k$-simplicial signal
$\mathbf{x}^{k}=[x^{k}_{1},\dots,x^{k}_{N_{k}}]$ as a mapping from the
$k$-simplices to the real space $\mathbb{R}^{N_{k}}$. For example,
$\mathbb{R}^{N_{0}}$ is the graph signal space in GSP, and
$\mathbb{R}^{N_{1}}$ is the space of edge flows. For an edge flow
$\mathbf{x}^{1}\in\mathbb{R}^{N_{1}}$, the sign of its entry denotes the
direction of the flow relative to a chosen reference orientation [11, 17].

Hodge Laplacian and decomposition. We represent the relations between
$(k-1)$\- and $k$-simplices with the incidence matrix $\mathbf{B}_{k}$,
$k=1,\dots,K$. The rows of $\mathbf{B}_{k}$ are indexed by $(k-1)$-simplices
and the columns by $k$-simplices. E.g., matrix $\mathbf{B}_{1}$ is the node-
to-edge incidence matrix, and $\mathbf{B}_{2}$ is the edge-to-triangle
incidence matrix [11, 12].

We can also use the Hodge Laplacian matrices,
$\mathbf{L}_{k}=\mathbf{B}_{k}^{\top}\mathbf{B}_{k}+\mathbf{B}_{k+1}\mathbf{B}_{k+1}^{\top}$,
where $\mathbf{L}_{k,\rm{l}}\triangleq\mathbf{B}_{k}^{\top}\mathbf{B}_{k}$ and
$\mathbf{L}_{k,\rm{u}}\triangleq\mathbf{B}_{k+1}\mathbf{B}_{k+1}^{\top}$ are
the lower and the upper Laplacians, which encode lower and upper
neighbourhoods, respectively. When $k=0$, the Hodge Laplacian is the graph
Laplacian $\mathbf{L}_{0}=\mathbf{B}_{1}\mathbf{B}_{1}^{\top}$. Hodge
Laplacians admit a _Hodge decomposition_ , leading to three orthogonal
subspaces which the simplicial signal space can be decomposed into, i.e.,
$\mathbb{R}^{N_{k}}=\text{im}(\mathbf{B}_{k}^{\top})\oplus\text{im}(\mathbf{B}_{k+1})\oplus\text{ker}(\mathbf{L}_{k})$,
where $\oplus$ is the direct sum of vector spaces and $\text{im}(\cdot)$ and
$\text{ker}(\cdot)$ are the _image_ and _kernel_ of a matrix. For $k=1$, these
subspaces carry the following interpretations [7, 11].

Gradient space. By applying matrix $\mathbf{B}_{1}$ to an edge flow
$\mathbf{x}^{1}$, we compute its net flow at each node,
$\mathbf{B}_{1}\mathbf{x}^{1}$. The incidence matrix $\mathbf{B}_{1}$ is
called a _divergence operator_. Its adjoint $\mathbf{B}_{1}^{\top}$
differentiates a node signal $\mathbf{x}^{0}$ along the edges to induce an
edge flow $\mathbf{B}_{1}^{\top}\mathbf{x}^{0}$, i.e., it is the _gradient
operator_. As a result, any flow within $\text{im}(\mathbf{B}_{1}^{\top})$ can
be written as the gradient of a node signal $\mathbf{x}^{0}$, i.e.,
$\mathbf{x}^{1}=\mathbf{B}_{1}^{\top}\mathbf{x}^{0}$. We call
$\mathbf{x}^{1}\in\text{im}(\mathbf{B}_{1}^{\top})$ a _gradient flow_ and the
space $\text{im}(\mathbf{B}_{1}^{\top})$ the gradient space.

Curl space. We can induce an edge flow from a triangle signal $\mathbf{x}^{2}$
as $\mathbf{x}^{1}=\mathbf{B}_{2}\mathbf{x}^{2}$. The induced flow
$\mathbf{x}^{1}\in\text{im}(\mathbf{B}_{2})$ is called a _curl flow_ and the
space $\text{im}(\mathbf{B}_{2})$ is the curl space. The adjoint
$\mathbf{B}_{2}^{\top}$ is the curl operator. We can use it to compute the net
edge flow of $\mathbf{x}^{1}$ circulating along the triangles as
$\mathbf{B}_{2}^{\top}\mathbf{x}^{1}$.

Harmonic space. The remaining space $\text{ker}(\mathbf{L}_{1})$ is the
harmonic space. Any edge flow $\mathbf{x}^{1}\in\text{ker}(\mathbf{L}_{1})$
has zero divergence and curl, i.e., it is _divergence-_ and _curl-free_.

Due to the boundary condition $\mathbf{B}_{1}\mathbf{B}_{2}=\mathbf{0}$, any
gradient flow $\mathbf{x}^{1}\in\text{im}(\mathbf{B}_{1}^{\top})$ is curl-
free. The space orthogonal to the gradient space, i.e.,
$\text{ker}(\mathbf{B}_{1})=\text{im}(\mathbf{B}_{2})\oplus\text{ker}(\mathbf{L}_{1})$,
is called the cycle space, which consists of both the curl space and harmonic
space. Any flow in this space is divergence-free.

Simplicial filters. To process simplicial signals $\mathbf{x}^{1}$, we use a
simplicial convolutional filter of the following form:

$\mathbf{H}=\epsilon\mathbf{I}+\sum_{l_{1}=1}^{L_{1}}\alpha_{l_{1}}(\mathbf{B}_{1}^{\top}\mathbf{B}_{1})^{l_{1}}+\sum_{l_{2}=1}^{L_{2}}\beta_{l_{2}}(\mathbf{B}_{2}\mathbf{B}_{2}^{\top})^{l_{2}}$
(1)

where $\epsilon$,
$\boldsymbol{\alpha}=[\alpha_{1},\dots,\alpha_{L_{1}}]^{\top}$ and
$\boldsymbol{\beta}=[\beta_{1},\dots,\beta_{L_{2}}]^{\top}$ are the filter
coefficients [12]. For ease of exposition, we only discuss the simplicial
filter form (1) for the edge signal space $\mathbb{R}^{N_{1}}$, but similar
discussions apply to general simplicial filters with $\mathbf{B}_{k}$ and
$\mathbf{B}_{k+1}$.

Applying $\mathbf{H}$ to an input edge flow $\mathbf{x}^{1}$ consists of the
_simplicial shifting_ operations, $\mathbf{L}_{1,\rm{l}}\mathbf{x}^{1}$ and
$\mathbf{L}_{1,\rm{u}}\mathbf{x}^{1}$. The $i$th entry of
$\mathbf{L}_{1,\rm{l}}\mathbf{x}^{1}$ is
$[\mathbf{L}_{1,\rm{l}}\mathbf{x}^{1}]_{i}=\sum_{j\in\\{\mathcal{N}_{{\rm
l},i}\cup\;i\\}}[\mathbf{L}_{1,{\rm l}}]_{ij}[\mathbf{x}^{1}]_{j}$, which is a
local operation within the lower neighborhood of the $i$th edge, likewise for
$\mathbf{L}_{1,\rm{u}}\mathbf{x}^{1}$. Moreover, powers
$\mathbf{L}_{1,\rm{l}}^{k}\mathbf{x}^{1}=\mathbf{L}_{1,\rm{l}}(\mathbf{L}_{1,\rm{l}}^{k-1}\mathbf{x}^{1})$
can be recursively obtained by applying the local operation $k$ times [12].
This leads to a distributed implementation of simplicial filtering with a
complexity of order $\mathcal{O}(N_{1}D)$ for each shifting with $D$ being the
maximal number of neighbours.

As we can observe from (1), in the simplicial domain, different sets of
coefficients on $\mathbf{L}_{1,\rm{l}}$ and $\mathbf{L}_{1,\rm{u}}$ enable an
independent and flexible filtering within the lower and upper simplicial
neighbourhoods. When $L_{1}=L_{2}$ and
$\boldsymbol{\alpha}=\boldsymbol{\beta}$, filter $\mathbf{H}$ [cf. (1)]
reduces to the basic form $\mathbf{H}=\sum_{l=0}^{L}h_{l}\mathbf{L}_{1}^{l}$
at the cost of losing expressive power and flexibility [12].

##  3 Simplicial convolutional neural networks

Upon having a simplicial convolutional filter (1), we can build an SCNN by
composing filter banks with elementwise nonlinearities. This SCNN applies to
any $k$-simplicial signal but we again illustrate it for edge signals
$\mathbf{x}^{1}$ for the ease of intuition, and omit the superscript to avoid
overcrowded notation.

Consider a $P$-layer SCNN. In the first layer $p=1$, we apply $F$ filters
$\mathbf{H}_{1}^{f}$ [cf. (1)] and an elementwise nonlinearity $\sigma(\cdot)$
to the input $\mathbf{x}_{0}$ to get a collection of $F$ features
$\mathbf{x}_{1}^{f}$ as

$\mathbf{x}_{1}^{f}=\sigma[\mathbf{z}_{1}^{f}]=\sigma[\mathbf{H}_{1}^{f}\mathbf{x}_{0}],\quad
f=1,\dots,F$ (2)

which constitute the output feature matrix
$\mathbf{X}_{1}=[\mathbf{x}_{1}^{1},\dots,\mathbf{x}_{1}^{F}]$. In subsequent
intermediate layers $p=2,\dots,P-1$, we have
$\mathbf{X}_{p-1}=[\mathbf{x}_{p-1}^{1},\dots,\mathbf{x}_{p-1}^{F}]\in\mathbb{R}^{N_{1}\times
F}$ as input. Each input signal $\mathbf{x}_{p-1}^{g},g=1,\dots,F$ is passed
through a bank of filters $\mathbf{H}_{p}^{fg}$ to obtain $F$ intermediate
outputs
$\mathbf{z}_{p}^{fg}=\mathbf{H}_{p}^{fg}\mathbf{x}_{p-1}^{g},f=1,\dots,F$. To
avoid exponential filter growth, the intermediate outputs of the different
input signals $\mathbf{x}_{p-1}^{g}$ are summed, thus, the $p$th layer
generates $F$ features $\mathbf{x}_{p}^{f}$ as follows

$\mathbf{x}_{p}^{f}=\sigma\bigg{[}\sum_{g=1}^{F}\mathbf{z}_{p}^{fg}\bigg{]}=\sigma\bigg{[}\sum_{g=1}^{F}\mathbf{H}_{p}^{fg}\mathbf{x}_{p-1}^{g}\bigg{]}~{}f=1,\dots,F.$
(3)

The processing in (3) is repeated until the last layer $p=P$, where we
consider the output has a single feature, and hence each input is processed by
a single filter $\mathbf{H}^{g}$. Thus, the final output of the SCNN is given
by

$\mathbf{x}_{P}=\sigma\bigg{[}\sum_{g=1}^{F}\mathbf{z}_{P}^{g}\bigg{]}=\sigma\bigg{[}\sum_{g=1}^{F}\mathbf{H}_{P}^{g}\mathbf{x}_{P-1}^{g}\bigg{]}.$
(4)

Equations (2), (3) and (4) constitute the SCNN architecture based on the
simplicial filter form defined in (1). The lower and upper Hodge Laplacians
encode the lower and upper simplicial neighbourhoods, respectively. Simplicial
convolutions through filter (1) are performed independently within the lower
and upper neighbourhoods and controlled by different sets of coefficients. As
we shall show later on, this means that the gradient and curl components of
the input features are convolved independently, leading to more expressive
power. Next, we discuss the connections of the SCNN with current alternatives
and analyze its properties.

Links with related works. In [8], a similar convolutional layer was proposed
but based on filter $\mathbf{H}=\sum_{l=0}^{L}h_{l}\mathbf{L}_{1}^{l}$. This
SNN architecture is a particular case of the proposed SCNN with less
expressive power but also with less parameters. The message passing neural
network (MPNN) for simplicial complexes [13] aggregates and updates features
from direct simplicial neighbours and simplices of different orders, e.g.,
nodes and triangles. By considering an order-one simplicial convolution as the
message aggregation step, we then obtain the architectures in [14, Eq. 4] and
[15, Eq. 7]. When only edge features are available, such approaches are a
particular case of the SCNN with order $L_{1}=1$ and $L_{2}=1$. Recurrent
architectures are considered for flow interpolation and graph classification
in [16]. Compared to these works, the SCNN treats features from the lower and
upper neighbours differently and considers features from not only direct
neighbours but also for multihop neighbors.

Locality and complexity. The intermediate output $\mathbf{z}^{fg}_{p}$ at the
$p$th layer collects for each edge information from lower neighbours up to
$L_{1}$ hops away and upper neighbours $L_{2}$ hops away through filter
(1).This locality comes from the structure of the Hodge Laplacian, likewise
that of GNNs [2].

When only a single feature is available, we have $1+L_{1}+L_{2}$ parameters in
such layers. For layers with multiple input and output features, the number of
parameters grows $F^{2}$ times. The major complexity comes from the
convolutional filtering step, which as seen before it is a weighted linear
combination of different shifts of a simplicial signal; a local operation
within the simplicial neighbourhoods that can be computed recursively. Hence,
an SCNN layer performs the simplicial filtering for each edge with a cost of
order $\mathcal{O}((L_{1}+L_{2})D)$. Again, this complexity grows $F^{2}$
times when multiple features are used and $P$ times if $P$ layers are
considered.

Equivariance and invariance. In the following, we show that our SCNN is
equivalent with respect to a different simplex labeling and different
reference flow orientations. Consider the set of simplicial permutation
matrices

$\mathcal{P}=\\{\mathbf{P}_{k}\in\\{0,1\\}^{N_{k}\times
N_{k}}:\mathbf{P}_{k}\mathbf{1}=\mathbf{1},\mathbf{P}_{k}^{\top}\mathbf{1}=\mathbf{1},k\geq
0\\},$

where products $\mathbf{P}_{k}\mathbf{x}^{k}$ permute the $k$-simplicial
signal $\mathbf{x}^{k}$. Let
$\mathbf{P}=(\mathbf{P}_{0},\mathbf{P}_{1},\dots)$ denote a sequence of
permutations. Then, the following holds.

######  Proposition 1.

The SCNN is a permutation equivariant architecture. For an input edge flow
$\mathbf{x}$, the output of an edge space SCNN layer with a simplicial filter
$\mathbf{H}$, $\mathbf{y}=\sigma[\mathbf{H}\mathbf{x}]$, becomes
$\mathbf{y}^{\prime}=\mathbf{P}_{1}\mathbf{y}$ after a permutation sequence
$\mathbf{P}$.

###### Proof.

An SCNN layer with filter $\mathbf{H}$ gives the output
$\mathbf{z}=\mathbf{H}\mathbf{x}$. After a sequence of permutations
$\mathbf{P}$, the input edge flow $\mathbf{x}$ becomes
$\mathbf{P}_{1}\mathbf{x}$ and the boundary operators become
$\mathbf{P}_{0}\mathbf{B}_{1}\mathbf{P}_{1}^{\top}$ and
$\mathbf{P}_{1}\mathbf{B}_{2}\mathbf{P}_{2}^{\top}$. Thus, the Hodge
Laplacians become $\mathbf{P}_{1}\mathbf{L}_{1,\rm{l}}\mathbf{P}_{1}^{\top}$
and $\mathbf{P}_{1}\mathbf{L}_{1,\rm{u}}\mathbf{P}_{1}^{\top}$ due to
$\mathbf{P}_{k}^{\top}\mathbf{P}_{k}=\mathbf{I}$. Then we can express the
permuted intermediate output as

$\displaystyle\mathbf{z}^{\prime}$
$\displaystyle\\!=\\!\bigg{(}\\!\epsilon\mathbf{I}\\!+\\!\sum_{l_{1}=1}^{L_{1}}\\!\alpha_{l_{1}}(\mathbf{P}_{1}\mathbf{L}_{1,\rm{l}}\mathbf{P}_{1}^{\top})^{l_{1}}\\!+\\!\\!\\!\sum_{l_{2}=1}^{L_{2}}\\!\beta_{l_{2}}(\mathbf{P}_{1}\mathbf{L}_{1,\rm{u}}\mathbf{P}_{1}^{\top})^{l_{2}}\\!\bigg{)}\mathbf{P}_{1}\mathbf{x}$
$\displaystyle=\mathbf{P}_{1}\bigg{(}\epsilon\mathbf{I}+\sum_{l_{1}=1}^{L_{1}}\alpha_{l_{1}}\mathbf{L}_{1,\rm{l}}^{l_{1}}+\sum_{l_{2}=1}^{L_{2}}\beta_{l_{2}}\mathbf{L}_{1,\rm{u}}^{l_{2}}\bigg{)}\mathbf{x}=\mathbf{P}_{1}\mathbf{z}.$

Thus, the simplicial filter $\mathbf{H}$ is permutation equivariant.
Furthermore, since the nonlinearity $\sigma(\cdot)$ is elementwise, SCNNs are
permutation equivariant. ∎

In addition, in a simplicial complex, we have also set an arbitrary reference
orientation for a simplex. A new reference orientation can be modelled by
multiplying the rows and columns of the boundary matrices $\mathbf{B}_{k}$ and
$\mathbf{B}_{k+1}$ where that $k$-simplex appears and the corresponding
simplicial signal value by $-1$. Let then $\mathbf{D}_{k}$ be diagonal
matrices from the set

$\mathcal{D}=\\{\mathbf{D}_{k}=\text{diag}(\mathbf{d}_{k}):\mathbf{d}_{k}\in\\{\pm
1\\}^{N_{k}},k\geq 1,\mathbf{d}_{0}=\mathbf{1}\\},$

where $\mathbf{D}_{k}\mathbf{x}^{k}$ is the updated $k$-simplicial signal
$\mathbf{x}^{k}$. Let $\mathbf{D}=(\mathbf{D}_{0},\mathbf{D}_{1},\dots)$
denote a sequence of orientation changes. Then, the following holds.

######  Proposition 2.

The SCNN is orientation equivariant if the nonlinearity $\sigma(\cdot)$ is
odd. Without loss of generality, for an input flow $\mathbf{x}$, the output of
an edge space SCNN layer with a simplicial filter $\mathbf{H}$,
$\mathbf{y}=\sigma[\mathbf{H}\mathbf{x}]$ becomes
$\mathbf{y}^{\prime}=\mathbf{D}_{1}\mathbf{y}$ after a sequence of orientation
changes $\mathbf{D}$.

###### Proof.

After an orientation change, the edge flow $\mathbf{x}$ becomes
$\mathbf{D}_{1}\mathbf{x}$ and the boundary operators $\mathbf{B}_{1}$ and
$\mathbf{B}_{2}$ are updated as $\mathbf{D}_{0}\mathbf{B}_{1}\mathbf{D}_{1}$
and $\mathbf{D}_{1}\mathbf{B}_{2}\mathbf{D}_{2}$. Then, the Hodge Laplacians
become $\mathbf{D}_{1}\mathbf{L}_{1,\rm{l}}\mathbf{D}_{1}$ and
$\mathbf{D}_{1}\mathbf{L}_{1,\rm{u}}\mathbf{D}_{1}$. We can then express the
new filter output as
$\mathbf{z}^{\prime}=\mathbf{D}_{1}\mathbf{H}\mathbf{x}=\mathbf{D}_{1}\mathbf{z}$,
following similar steps as in (3). Thus, simplicial filter $\mathbf{H}$ is
orientation equivariant. When the nonlinearity $\sigma(\cdot)$ is an odd
function, we have
$\mathbf{y}^{\prime}=\sigma(\mathbf{z}^{\prime})=\mathbf{D}_{1}\sigma(\mathbf{z})=\mathbf{D}_{1}\mathbf{y}$
that completes the proof. ∎

Permutation and orientation equivariances preserve the output of an SCNN
regardless of the choices of the labeling and reference orientation of the
edges. In turn, they allow the SCNN to exploit the internal symmetries in the
complex.

Spectral analysis. For the spectral analysis, consider first that the
eigenvectors of the Hodge Laplacian $\mathbf{L}_{1}$ span the three spaces
given by the Hodge decomposition: (i) the gradient space
$\text{im}(\mathbf{B}_{1}^{\top})$ is spanned by a set of eigenvectors
$\mathbf{U}_{\rm{G}}$ of $\mathbf{L}_{1,\rm{l}}$ with positive eigenvalues;
(ii) the curl space $\text{im}(\mathbf{B}_{2})$ is spanned by a set of
eigenvectors $\mathbf{U}_{\rm{C}}$ of $\mathbf{L}_{1,\rm{u}}$ with positive
eigenvalues; and (iii) the harmonic space $\text{ker}(\mathbf{L}_{1})$ is
spanned by the eigenvectors $\mathbf{U}_{\rm{H}}$ of $\mathbf{L}_{1}$ with
zero eigenvalue. Moreover, we have
$\text{im}(\mathbf{L}_{1})=\text{im}(\mathbf{U}_{\rm{G}})\oplus\text{im}(\mathbf{U}_{\rm{C}})$,
i.e., gradient and curl spaces make up the image of $\mathbf{L}_{1}$ [12].
Then, we can eigendecompose $\mathbf{L}_{1}$ as
$\mathbf{L}_{1}=\mathbf{U}\mathbf{\Lambda}\mathbf{U}^{\top}$ where
$\mathbf{U}=[\mathbf{U}_{\rm{H}}\,\,\mathbf{U}_{\rm{G}}\,\,\mathbf{U}_{\rm{C}}]$
provides a simplicial Fourier basis, and
$\mathbf{\Lambda}=\text{diag}(\mathbf{\Lambda}_{\rm{H}},\mathbf{\Lambda}_{\rm{G}},\mathbf{\Lambda}_{\rm{C}})$
with $\mathbf{\Lambda}_{\rm{H}}=\text{diag}(\mathbf{0}_{N_{H}})$,
$\mathbf{\Lambda}_{\rm{G}}=\text{diag}(\lambda_{\rm{G},1},\dots,\lambda_{{\rm{G}},N_{\rm{G}}})$,
and
$\mathbf{\Lambda}_{\rm{C}}=\text{diag}(\lambda_{\rm{C},1},\dots,\lambda_{{\rm{C}},N_{\rm{C}}})$
collecting the harmonic, gradient, and curl frequencies, respectively; i.e.,
the simplicial frequencies [12].

For an edge flow $\mathbf{x}$, we can find its simplicial Fourier transform
(SFT) as $\tilde{\mathbf{x}}=\mathbf{U}^{\top}\mathbf{x}$. This further
defines three embeddings
$\tilde{\mathbf{x}}=[\tilde{\mathbf{x}}_{\rm{H}}^{\top}\,\,\tilde{\mathbf{x}}_{\rm{G}}^{\top}\,\,\tilde{\mathbf{x}}_{\rm{C}}^{\top}]$:
the harmonic embedding
$\tilde{\mathbf{x}}_{\rm{H}}=\mathbf{U}^{\top}_{\rm{H}}\mathbf{x}$, the
gradient embedding
$\tilde{\mathbf{x}}_{\rm{G}}=\mathbf{U}^{\top}_{\rm{G}}\mathbf{x}$, and the
curl embedding
$\tilde{\mathbf{x}}_{\rm{C}}=\mathbf{U}^{\top}_{\rm{C}}\mathbf{x}$, which
contain the weights of $\mathbf{x}$ at harmonic, gradient, and curl
frequencies, respectively.

(a) $10\%$ missing rate.

(b) $20\%$ missing rate.

Fig. 1: Training loss of SNN and SCNN. We see that our SCNN converges much
faster in the early training stage and results in a smaller training loss.

Via the eigendecomposition of $\mathbf{L}_{1}$, we can analyze the proposed
SCNNs in the spectral domain. First, the frequency response of a simplicial
convolutional filter $\mathbf{H}$ is given by

$\tilde{H}(\lambda_{i})=\begin{cases}\epsilon,&\text{ for }\lambda_{i}=0,\\\
\epsilon+\sum_{l_{1}=1}^{L_{1}}\alpha_{l_{1}}\lambda_{i}^{l_{1}},&\text{ for
}\lambda_{i}\in\mathcal{Q}_{\rm{G}},\\\
\epsilon+\sum_{l_{2}=1}^{L_{2}}\beta_{l_{2}}\lambda_{i}^{l_{2}},&\text{ for
}\lambda_{i}\in\mathcal{Q}_{\rm{C}},\\\ \end{cases}$ (5)

where $\mathcal{Q}_{\rm{G}}$ and $\mathcal{Q}_{\rm{C}}$ respectively collect
distinct gradient and curl frequencies. Thus, the $i$th entry of the SFT of
the intermediate output $\mathbf{z}$ of an SCNN layer can be expressed as
$\tilde{z}_{i}=\tilde{H}(\lambda_{i})\tilde{x}_{i}$. This spectral analysis
shows that the SCNN layers compute the high-level simplicial components as the
pointwise multiplication of the input embedding and the simplicial filter
frequency response, ultimately respecting the convolution theorem.
Furthermore, we have here different frequency responses for the gradient and
curl components, which corresponds to the different weights on lower and upper
Laplacians in filter $\mathbf{H}$. However, for a filter with
$\boldsymbol{\alpha}=\boldsymbol{\beta}$ as in [8], this would lead to
coupling between independent frequencies, and in turn to a limited learning
expressiveness.

##  4 Numerical results

We use the SCNN to impute missing citations in a coauthorship complex, in
which a paper with $k+1$ authors is represented by a $k$-simplex, and the
$k$-simplicial signal is the number of citations of the paper. We followed the
steps of [8], which lead to the citation dataset in Table 1. We compared the
SCNN with the SNN in [8] for the $k$-simplicial signals with $k=0,\dots,5$.
Missing data are generated randomly on the $k$-simplicial signals at 5 rates,
$10\%,20\%,\dots,50\%$. The input of the SCNNs is the $k$-simplicial signal
where missing citations are replaced by the median of known citations. As the
SNN in [8], our SCNN has 3 layers with 30 convolutional filters of total
length $5$ ($L_{1}=L_{2}=2$). We used LeakyReLU for $\sigma(\cdot)$ as in [8]
although not odd. The reference orientation didn’t seem to have much
influence. We used the $\ell_{1}$ norm to train the NNs over known citations
for 1000 iterations using the Adam optimizer with a learning rate $10^{-3}$.

We report the training loss of two instances in Fig. 1 and the mean
accuracy111A citation value is considered to be correct if the imputed value
is within $\pm 5\%$ of the true value. $\pm$ the standard deviation in Table 1
over 10 different experiments. The proposed SCNN approach achieves a smaller
training loss and a faster convergence than the SNN (Fig. 1) due to its better
expressive power. From Table 1, we observe that both NNs perform similarly for
dimensions 0 and 1. This is because for the former, the two NNs are the same,
and for the latter, the data dimension is rather small. However, the SCNN
gives consistently $1-2\%$ better accuracies for $k\geq 2$ with larger data
dimensions.

Table 1: Imputation accuracies for each dimension and missing rate by SNN
(first rows) and SCNN (second rows).

Order | 0 | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---|---
$N_{k}$ | 352 | 1474 | 3285 | 5019 | 5559 | 4547
$10\%$ | $0.91\pm 0.003$ | $0.91\pm 0.002$ | $0.90\pm 0.004$ | $0.91\pm 0.004$ | $0.90\pm 0.016$ | $0.90\pm 0.008$
$10\%$ | $0.91\pm 0.004$ | $0.91\pm 0.002$ | $0.91\pm 0.002$ | $0.92\pm 0.001$ | $0.92\pm 0.002$ | $0.92\pm 0.002$
$20\%$ | $0.81\pm 0.006$ | $0.82\pm 0.003$ | $0.82\pm 0.005$ | $0.83\pm 0.004$ | $0.82\pm 0.012$ | $0.83\pm 0.007$
$20\%$ | $0.81\pm 0.007$ | $0.82\pm 0.003$ | $0.83\pm 0.003$ | $0.83\pm 0.002$ | $0.84\pm 0.002$ | $0.84\pm 0.002$
$30\%$ | $0.72\pm 0.006$ | $0.73\pm 0.004$ | $0.73\pm 0.005$ | $0.75\pm 0.002$ | $0.75\pm 0.002$ | $0.75\pm 0.003$
$30\%$ | $0.72\pm 0.005$ | $0.73\pm 0.004$ | $0.74\pm 0.003$ | $0.75\pm 0.002$ | $0.76\pm 0.002$ | $0.77\pm 0.002$
$40\%$ | $0.63\pm 0.007$ | $0.64\pm 0.003$ | $0.65\pm 0.003$ | $0.66\pm 0.004$ | $0.67\pm 0.009$ | $0.67\pm 0.008$
$40\%$ | $0.63\pm 0.006$ | $0.64\pm 0.003$ | $0.65\pm 0.002$ | $0.66\pm 0.002$ | $0.67\pm 0.003$ | $0.69\pm 0.002$
$50\%$ | $0.54\pm 0.007$ | $0.55\pm 0.005$ | $0.56\pm 0.003$ | $0.57\pm 0.003$ | $0.59\pm 0.004$ | $0.60\pm 0.005$
$50\%$ | $0.54\pm 0.006$ | $0.55\pm 0.004$ | $0.56\pm 0.003$ | $0.58\pm 0.003$ | $0.59\pm 0.003$ | $0.61\pm 0.002$
|  |  |  |  |  |

##  5 Conclusion

This paper proposed a simplicial convolutional neural network architecture to
learn from data defined on higher-order structures of a network, i.e.,
simplices in a simplicial complex. We built an SCNN layer through a
composition of a simplicial filter and an elementwise nonlinearity. Due to the
use of an advanced simplicial filter, our SCNN is able to learn from
simplicial neighbours over multiple hops and process simplicial subcomponents
(e.g., gradient and curl) independently, compared with the current solutions.
The proposed SCNN applies to any $k$-simplicial signal case. We showed the
SCNN is equivariant to permutations in the topology and to orientations in the
flows, which allows it to exploit symmetries in a simplicial complex. In the
future, we plan to extend the SCNN to include also the data on adjacent
simplices.

## References

  * [1] A. Ortega, P. Frossard, J. Kovačević, J. M. Moura, and P. Vandergheynst, “Graph signal processing: Overview, challenges, and applications,” _Proceedings of the IEEE_ , vol. 106, no. 5, pp. 808–828, 2018.
  * [2] F. Gama, E. Isufi, G. Leus, and A. Ribeiro, “Graphs, convolutions, and neural networks: From graph filters to graph neural networks,” _IEEE Signal Processing Magazine_ , vol. 37, no. 6, pp. 128–138, 2020.
  * [3] E. Isufi, M. Pocchiari, and A. Hanjalic, “Accuracy-diversity trade-off in recommender systems via graph convolutions,” _Information Processing & Management_, vol. 58, no. 2, p. 102459, 2021.
  * [4] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” _arXiv preprint arXiv:1609.02907_ , 2016.
  * [5] E. Isufi, F. Gama, and A. Ribeiro, “Edgenets: Edge varying graph neural networks,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2021.
  * [6] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, and S. Y. Philip, “A comprehensive survey on graph neural networks,” _IEEE transactions on neural networks and learning systems_ , vol. 32, no. 1, pp. 4–24, 2020.
  * [7] S. Barbarossa and S. Sardellitti, “Topological signal processing over simplicial complexes,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 2992–3007, 2020.
  * [8] S. Ebli, M. Defferrard, and G. Spreemann, “Simplicial neural networks,” _arXiv preprint arXiv:2010.03633_ , 2020.
  * [9] G. Leus, M. Yang, M. Coutino, and E. Isufi, “Topological volterra filters,” in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 5385–5399.
  * [10] M. T. Schaub and S. Segarra, “Flow smoothing and denoising: Graph signal processing in the edge-space,” in _2018 IEEE Global Conference on Signal and Information Processing (GlobalSIP)_. IEEE, 2018, pp. 735–739.
  * [11] M. T. Schaub, Y. Zhu, J.-B. Seby, T. M. Roddenberry, and S. Segarra, “Signal processing on higher-order networks: Livin’on the edge… and beyond,” _arXiv preprint arXiv:2101.05510_ , 2021.
  * [12] M. Yang, E. Isufi, M. T. Schaub, and G. Leus, “Finite impulse response filters for simplicial complexes,” _arXiv preprint arXiv:2103.12587_ , 2021.
  * [13] C. Bodnar, F. Frasca, Y. G. Wang, N. Otter, G. Montúfar, P. Lio, and M. Bronstein, “Weisfeiler and lehman go topological: Message passing simplicial networks,” _arXiv preprint arXiv:2103.03212_ , 2021.
  * [14] E. Bunch, Q. You, G. Fung, and V. Singh, “Simplicial 2-complex convolutional neural nets,” _arXiv preprint arXiv:2012.06010_ , 2020.
  * [15] T. M. Roddenberry, N. Glaze, and S. Segarra, “Principled simplicial neural networks for trajectory prediction,” in _International Conference on Machine Learning_. PMLR, 2021, pp. 9020–9029.
  * [16] T. M. Roddenberry and S. Segarra, “Hodgenet: Graph neural networks for edge data,” in _2019 53rd Asilomar Conference on Signals, Systems, and Computers_. IEEE, 2019, pp. 220–224.
  * [17] L.-H. Lim, “Hodge laplacians on graphs,” _SIAM Review_ , vol. 62, no. 3, pp. 685–715, 2020.