#  Equivariant, Safe and Sensitive
— Graph Networks for New Physics

Akanksha Bhardwaj  0000-0001-8695-972X b Christoph Englert
0000-0003-2201-0667 b Wrishik Naskar  0000-0002-4357-899 c Vishal S.
Ngairangbam  0000-0002-7143-715X c and Michael Spannowsky  0000-0002-8362-0576
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Physics, Oklahoma State
University, Stillwater, OK, 74078, USASchool of Physics and Astronomy,
University of Glasgow, Glasgow G12 8QQ, United KingdomInstitute for Particle
Physics Phenomenology, Department of Physics,
Durham University, Durham DH1 3LE, United Kingdom

###### Abstract

This study introduces a novel Graph Neural Network (GNN) architecture that
leverages infrared and collinear (IRC) safety and equivariance to enhance the
analysis of collider data for Beyond the Standard Model (BSM) discoveries. By
integrating equivariance in the rapidity-azimuth plane with IRC-safe
principles, our model significantly reduces computational overhead while
ensuring theoretical consistency in identifying BSM scenarios amidst Quantum
Chromodynamics backgrounds. The proposed GNN architecture demonstrates
superior performance in tagging semi-visible jets, highlighting its potential
as a robust tool for advancing BSM search strategies at high-energy colliders.

###### Keywords:

Jet substructure, Equivariance, Dark showers

††preprint:  IPPP/24/07

##  1 Introduction

The application of machine learning algorithms to jet classification
Cogan:2014oua ; deOliveira:2015xxd ; Kasieczka:2017nvn ; Louppe:2017ipp ;
Dreyer:2018nbf ; Heimel:2018mkt ; Farina:2018fyg ; Qu:2019gqs ;
Kasieczka:2019dbj ; Bernreuther:2020vhm ; pmlr-v162-qu22b ; Kasieczka:2020nyd
; Dreyer:2020brq ; Khosa:2021cyk ; Canelli:2021aps ; Dreyer:2021hhr ;
Cavallini:2021vot ; Butter:2022xyj ; Dillon:2022mkq ; Faucett:2022zie ;
Keicher:2023mer ; Finke:2023veq ; Lu:2023gjk ; Bardhan:2023mia ;
Metodiev:2023izu ; Favaro:2023xdl ; Romero:2023hrk ; Furuichi:2023vdx ;
Gaertner:2023ycs  provides an ideal environment to gauge the interplay between
performance and interpretability. On the one hand, we need highly performant
algorithms utilising the wealth of experimental data recorded at the Large
Hadron Collider (LHC), while on the other, we want to ascertain the reason
behind the algorithm’s efficiency. Architectures which are Lorentz equivariant
Bogatskiy:2020tje ; Gong:2022lye ; Li:2022xfc ; Hao:2022zns ;
Bogatskiy:2022czk , or are infra-red and collinear (IRC) safe Choi:2018dag ;
Komiske:2018cqr ; Konar:2021zdg ; Atkinson:2022uzb ; Athanasakos:2023fhq ;
Shen:2023ofd ; Konar:2023ptv  have been shown to enhance the physical biases
of the algorithms with generally minimal loss, or at times, comparatively
increased performance Gong:2022lye ; Bogatskiy:2022czk  with respect to
physics-opaque constructions. In this work, we combine IRC safety with
Euclidean equivariance in the rapidity-azimuth plane. With well-known examples
of jet-shape Almeida:2008yp ; Thaler:2010tr ; Thaler:2011gf ; Larkoski:2013eya
; Larkoski:2014pca ; Moult:2016cvt ; Dasgupta:2015lxh ; Caletti:2021oor
variables, constructed as functions of distance in the rapidity-azimuth plane,
this pushes the IRC-safe features closer to human-engineered QCD features
without losing out on performance, even with two orders of magnitude lower
number of model parameters.

Symmetry in a different sense, namely the degeneracy over nearly degenerate
final states, is fundamental to the theoretical robustness of information that
can be gained from phenomenological analyses. The Kinoshita-Lee-Nauenberg
theorem Kinoshita:1962ur ; Lee:1964is  guarantees cancellations of soft and
(final-state) collinear divergences when algorithms are infrared (IR) and
collinear safe; any departure from IRC-safe methods ultimately implies
uncorrectable interpretation shortfalls. This technically elevates IRC safety
to a comparable level as the space-time symmetries to perturbatively model and
interpret scattering processes.

The recent surge in the application of machine learning for enhancing the
sensitivity to new physics in available and projected particle physics data
sets has led to a departure from observables that transparently reflect
symmetry properties and IRC safety. Nonetheless, any such serious attempt
should critically feature IRC safety whilst exploiting symmetries of the data
sets in the corresponding architecture. This work aims to demonstrate that
this is not only feasible but directly applicable to the determination of BSM
parameter regions in case a discovery can be claimed via these techniques. We
construct an IRC-safe $\mathbf{E(2)}$-equivariant network and demonstrate the
enhanced optimisation capability of such an algorithm compared to approaches
that do not exploit symmetry. We then apply this approach to a classification
of jet substructure analyses (we consider Hidden Valley Models
Strassler:2006im ; Carloni:2010tw ; Carloni:2011kk ; Cohen:2015toa ;
Kar:2020bws ; Bernreuther:2020xus ; Barron:2021btf ; Albouy:2022cin  as a
representative example), where IRC safety plays a critical role in obtaining a
theoretically robust output score when the substructure becomes sparse or
soft.

This work is organised as follows: In Sec. 2, we discuss the inductive biases
relevant to phenomenology at the LHC, namely IRC safety and equivariance in
the rapidity-azimuthal plane. Then, in Sec. 3, we apply this network to a
relevant physics case where it discriminates between QCD jets and dark shower
jets (semi-visible jets). We discuss the network architecture and its
performance. Finally, we provide a summary and conclusion in Sec. 4.

##  2 Inductive biases and QCD

The success of deep learning algorithms relies on reducing the search space of
possible functions to a predetermined subset by imposing additional structures
on the architecture. These inductive biases favour the extraction of relevant
features which are not too specific (like high-level variables) but follow the
intuition of the application domain. Here, we discuss the two inductive biases
we impose on the Message Passing Neural Network utilised in our study: IRC
safety and equivariance in the rapidity-azimuth plane.

###  2.1 Infrared and Collinear Safety

From a QCD perspective, infra-red and collinear safe feature extraction
promises to provide better interpretability while not losing out on the high
classification capabilities. Here, we briefly describe the IRC-safe Energy-
weighted Message Passing Network (EMPN) for jet classification first proposed
in reference Konar:2021zdg . The message-passing operation has the structure

$\mathbf{H}^{(l+1)}_{i}=\sum_{i\in\mathcal{N}[i]}\;\omega^{(\mathcal{N}[i])}_{j}\;\;\hat{\Phi}^{(l+1)}(\mathbf{H}^{(l)}_{i},\mathbf{H}^{(l)}_{j})\quad$
(1)

where the scope-dependent factors are an analogue of
$z_{i}=p^{i}_{T}/\sum_{k}p^{k}_{T}$, defined as

$\omega_{j}^{(\mathcal{N}[i])}=\frac{p_{T}^{j}}{\sum_{k\in\mathcal{N}[i]}\;p_{T}^{k}}\quad,$

with $\mathcal{N}[i]$ denoting the IRC-safe local neighbourhood111It is
important to note that IRC safety always demands a closed neighbourhood, where
$\mathcal{N}[i]$ should always contain $i$. of particle $i$. The function
$\hat{\Phi}^{(l+1)}$, in general, can be any well-behaved neural network which
we take to be an edge convolution network DBLP:journals/corr/abs-1801-07829 .
For the initial input, we use the relative coordinates in the rapidity-azimuth
plane from the jet axis, $\mathbf{H}^{(0)}_{i}=(\Delta
y_{iJ},\Delta\phi_{iJ})$. Given the neighbourhood condition and
$\mathbf{H}^{(0)}_{q}=\mathbf{H}^{(0)}_{r}=\mathbf{H}^{(0)}_{s}$, Eq. (1)
“conserves collinearity” of the updated features
$\mathbf{H}^{(l)}_{q}=\mathbf{H}^{(l)}_{r}=\mathbf{H}^{(l)}_{s}$ for any $l$
when $z_{q}=z_{r}+z_{s}$. In the following discussions, we will discuss the
conservation of collinearity when these conditions are satisfied. We will not
explicitly state the additional assumptions on the energy weights and the
neighbourhood. The neighbourhood $\mathcal{N}[i]$ is constructed by imposing
the step function $\Theta(\Delta R_{ij}<R_{0})$ for each particle $i$, $R_{0}$
being a parameter that determines the density of the constructed graph. The
IRC-safe features are obtained after an energy-weighted graph readout as

$\mathbf{G}_{l}=\sum_{i}z_{i}\;\mathbf{H}_{i}^{(l)}\quad.$ (2)

The graph features for $l>0$ can be concatenated up to the total number of
such operations $L$, which is a hyperparameter. For classification, a
downstream network takes the graph representation as an input and gives a
classification score.

###  2.2 Equivariance in rapidity-azimuth plane

Equivariance is at the core of the powerful learning ability shown by modern
deep-learning algorithms such as Convolutional Neural Networks (CNNs) and
Graph Neural Networks (GNNs). They reduce the possible set of functions that a
neural network can approximate by assuming symmetries of the input data. For
instance, CNNs assume translational symmetry in $\mathbb{R}^{2}$ while GNNs
assume permutation symmetry of the nodes. Equivariance enforces group
algebraic structures in the hidden representations in a neural network. This
reduces the set of possible functions that a neural network can approximate to
those which follow the group equivariant property.

Mathematically, for a group $\mathcal{G}$ with an element denoted as $g$, let
$\mathbf{X}$ and $\mathbf{Y}$ be two sets with $\mathbf{x}\in\mathbf{X}$ and
$\mathbf{y}\in\mathbf{Y}$ that admit group actions222An action of a group
$\mathcal{G}$ is defined as
$\mathbf{T}_{\mathbf{X}}:\mathcal{G}\times\mathbf{X}\to\mathbf{X}$ such that
for any $\mathbf{x}\in\mathbf{X}$, we have
$\mathbf{T}_{\mathbf{X}}(e,\mathbf{x})=\mathbf{x}$ for the neutral element
$e\in\mathcal{G}$ and
$\mathbf{T}_{\mathbf{X}}(g_{1},\mathbf{T}_{\mathbf{X}}(g_{2},\mathbf{x}))=\mathbf{T}_{\mathbf{X}}(g_{1}g_{2},\mathbf{x})$,
for all $g_{1},g_{2}\in\mathcal{G}$. $\mathbf{T}_{\mathbf{X}}(g,\mathbf{x})$
and $\mathbf{T}_{\mathbf{Y}}(g,\mathbf{y})$, respectively. A function
$f:\mathbf{X}\to\mathbf{Y}$ is $\mathcal{G}$-equivariant if

$f(\mathbf{T}_{\mathbf{X}}(g,\mathbf{x}))=\mathbf{T}_{\mathbf{Y}}(g,f(\mathbf{x}))\quad.$
(3)

An equivariant neural network Lim2023  approximates such a function $f$ for a
given group, and $f$ is generally non-linear. From an approximation
perspective, the neural network being able to learn only those functions that
follow Eq. (3) rather than any general function necessarily reduces the
expressive power. However, the reduction generally results in higher
efficiency of the optimisation process in terms of network complexity, on top
of increasing the interpretability as the particular group $\mathcal{G}$ is
determined from domain knowledge. If the action $\mathbf{T}_{\mathbf{Y}}$ is
trivial, i.e., $\mathbf{T}_{\mathbf{Y}}(g,\mathbf{y})=\mathbf{y}\;\forall
g\in\mathcal{G}$, then the function $f$ is said to be $\mathcal{G}$-invariant.
GNNs are permutation equivariant in the node update stage, while for graph-
level purposes, the global readout renders the graph representation
permutation invariant.

For our purpose, the group which preserves the separation $\Delta R$ in the
$(y,\phi)$ plane is the Euclidean group $\mathbf{E}(2)$. To incorporate
$\mathbf{E}(2)$-equivariance in an EMPN, we modify the general
$\mathbf{E}(n)$-equivariant message-passing operation of reference
pmlr-v139-satorras21a . Given $\mathbf{E}(n)$ scalar node features
$\mathbf{h}^{(l)}_{i}$ and vectors $\mathbf{x}^{(l)}_{i}$, it is given as

$\displaystyle\mathbf{m}^{(l+1)}_{ij}$
$\displaystyle=\Phi^{(l+1)}_{e}(\mathbf{h}^{(l)}_{i},\mathbf{h}^{(l)}_{j},|\mathbf{x}^{(l)}_{i}-\mathbf{x}^{(l)}_{j}|)\quad,$
(4a) $\displaystyle\mathbf{x}_{i}^{(l+1)}$
$\displaystyle=\mathbf{x}^{(l)}_{i}+\sum_{j\neq
i}\,\;(\mathbf{x}^{(l)}_{i}-\mathbf{x}^{(l)}_{j})\;\Phi^{(l+1)}_{x}(\mathbf{m}_{ij})\quad,$
(4b) $\displaystyle\mathbf{m}^{(l+1)}_{i}$ $\displaystyle=\sum_{j\neq
i}\,\;\mathbf{m}^{(l+1)}_{ij}\quad,$ (4c)
$\displaystyle\mathbf{h}^{(l+1)}_{i}$
$\displaystyle=\Phi_{h}(\mathbf{h}^{(l)}_{i},\mathbf{m}^{(l+1)}_{i})\quad,$
(4d)

where ${\Phi}^{(l+1)}_{e}$, ${\Phi}^{(l+1)}_{x}$, and ${\Phi}^{(l+1)}_{h}$ are
multilayer perceptrons (MLPs), with $\Phi^{(l+1)}_{x}$ giving a one-
dimensional output with a sigmoid activation, thereby being interpreted as a
scalar weight in the aggregation. The restriction $j\neq i$ in the aggregation
steps can be relaxed as the equivariance is satisfied in its absence.
Moreover, as the neighbourhood construction modifies the message functions
(${\Phi}^{(l+1)}_{e}$ and ${\Phi}^{(l+1)}_{x}$) by a product with a step
function $\Theta(\Delta R_{ij}<R_{0})$ which is itself $\mathbf{E}(2)$
invariant, we can incorporate a local aggregation scheme without impeding
global equivariance. We note that modifying Eq. (1) to

$\mathbf{H}^{(l+1)}_{i}=\mathbf{F}^{(l+1)}\left(\mathbf{H}_{i}^{(l)},\sum_{i\in\mathcal{N}[i]}\;\omega^{(\mathcal{N}[i])}_{j}\;\;\hat{\Phi}^{(l+1)}(\mathbf{H}^{(l)}_{i},\mathbf{H}^{(l)}_{j})\right)\,,$
(5)

with a well-behaved (to be learned) function $\mathbf{F}^{(l+1)}$ still
conserves collinearity of the updated node features. Equipped with this
observation, it is straightforward to modify Eq. (4) to satisfy IRC safety by
incorporating the energy-weighted sum structure in the aggregation as

$\displaystyle\mathbf{m}^{(l+1)}_{ij}$
$\displaystyle=\Phi^{(l+1)}_{e}(\mathbf{h}^{(l)}_{i},\mathbf{h}^{(l)}_{j},|\mathbf{x}^{(l)}_{i}-\mathbf{x}^{(l)}_{j}|^{2})\quad,$
(6a) $\displaystyle\mathbf{x}_{i}^{(l+1)}$
$\displaystyle=\mathbf{x}^{(l)}_{i}+\sum_{j\in\mathcal{N}[i]}\omega^{(\mathcal{N}[i])}_{j}\;(\mathbf{x}^{(l)}_{i}-\mathbf{x}^{(l)}_{j})\;\Phi^{(l+1)}_{x}(\mathbf{m}_{ij})\quad,$
(6b) $\displaystyle\mathbf{m}^{(l+1)}_{i}$
$\displaystyle=\sum_{j\in\mathcal{N}[i]}\omega^{(\mathcal{N}[i])}_{j}\;\mathbf{m}^{(l+1)}_{ij}\quad,$
(6c) $\displaystyle\mathbf{h}^{(l+1)}_{i}$
$\displaystyle=\Phi^{(l+1)}_{h}(\mathbf{h}^{(l)}_{i},\mathbf{m}^{(l+1)}_{i})\quad.$
(6d)

Conservation of collinearity of the updated scalar features
$\mathbf{h}^{(l+1)}_{i}$, and vector features $\mathbf{x}^{(l+1)}_{i}$ follows
from comparing their corresponding aggregation equations to Eq. (5). We have
modified the input Euclidean norm to the squared value as we have self-loops,
and the square root is not differentiable at zero.

The general $\mathbf{E}(n)$ group consists of translations and rotations.
While there are some instances where the translational symmetry is utilised by
defining vector observables like jet pull Gallicchio:2010sw ; ATLAS:2018olo ;
Larkoski:2019fsm , most substructure observables are confined to a
rotationally symmetric definition by using only the scalar separation $\Delta
R$. Therefore, we also study an $\mathbf{O}(2)$-equivariant architecture by
forgoing translational equivariance. This can be done by modifying Eqs. (6a)
and (6b) to

$\displaystyle\mathbf{m}^{(l+1)}_{ij}$
$\displaystyle=\Phi^{(l+1)}_{e}(\mathbf{h}^{(l)}_{i},\mathbf{h}^{(l)}_{j},|\mathbf{x}^{(l)}_{i}-\mathbf{x}^{(l)}_{j}|^{2},\mathbf{x}^{(l)}_{i}.\,\mathbf{x}^{(l)}_{j})\quad,$
(7a) $\displaystyle\mathbf{x}_{i}^{(l+1)}$
$\displaystyle=\sum_{j\in\mathcal{N}[i]}\omega^{(\mathcal{N}[i])}_{j}\;\mathbf{x}^{(l)}_{j}\;\Phi^{(l+1)}_{x}(\mathbf{m}_{ij})\quad,$
(7b)

and keeping Eqs. (6c) and (6d) unchanged. We have added the inner product
$\mathbf{x}^{(l)}_{i}.\,\mathbf{x}^{(l)}_{j}$ as it is invariant to rotations,
which is not the case for translations.

Note that for the final IRC safe feature extraction, we need to perform an
energy-weighted summed global readout as given in Eq. (2) for both
$\mathbf{x}^{(l)}_{i}$ and $\mathbf{h}^{(l)}_{i}$. Moreover, the input scalars
$\mathbf{h}^{(0)}_{i}$ cannot contain $z_{i}$ from the requirement of IRC
safety, even though it is invariant under $\mathbf{E}(2)$. A schematic
representation of the equivariant message passing is shown on the left of
figure 1.

Figure 1: A schematic representation of $\mathbf{E}(2)$-EMPN and
$\mathbf{O}(2)$-EMPN network architecture.

##  3 A Physics application: Tagging semi-visible jets

The radiation pattern arising from perturbative QCD emissions is known to be
approximately uniform in the rapidity-azimuth plane. This guides the design of
most jet-shape observables Almeida:2008yp  to use the relative distance
between $\Delta R$, between pairs of particles Larkoski:2013eya ;
Moult:2016cvt , or from well-defined axes Thaler:2010tr ; Thaler:2011gf . The
implicit assumption underlying such definitions is the invariance of the
observables under the two-dimensional Euclidean group $\mathbf{E}(2)$ in the
rapidity-azimuth plane. It is, therefore, a natural extension of IRC-safe
feature extraction to incorporate $\mathbf{E}(2)$-equivariance to better
connect to known jet substructure observables in literature.

Although equivariance, in principle, would help segregate different QCD
radiation patterns with n-prong structures, a particularly challenging case
would be when the parton shower itself is altered with non-trivial physics.
For example, a QCD-like dark shower in Hidden Valley models can produce a
mixture of BSM and normal SM particles, generally modelled via interleaved
parton showers. In such a case, the differences that arise from the QCD
radiation are not in the presence of the hard prongs but rather at an
intricate level in the relative distribution of the particle’s position in the
$(y,\phi)$ plane, with possibly some missing patches left by undetectable
stable or quasi-stable particles within the jet. Generally, the energy values
would be similar to a one-prong QCD jet. At the same time, the relative
position of the particles in the $(y,\phi)$ plane would be the distinguishing
feature in such semi-visible jets. Therefore, $\mathbf{E}(2)$ equivariance
would be more prominent in identifying such semi-visible jets.

Hidden Valley Models offer an intriguing perspective in particle physics,
suggesting the existence of secluded, almost decoupled sectors that interact
weakly with the SM. Such a scenario can be brought about by extending the SM
gauge group with a new group $G_{v}$. The choices of the gauge group and new
fermions and their relation to the SM particles offer a rich phenomenology
with possibly distinct signatures like emerging jets or semi-visible jets not
covered by traditional searches at the LHC, like monojet searches CMS:2021far
; ATLAS:2022ckd . In general, the SM particles remain neutral under the new
group $G_{v}$, while there can be two classes of new particles: those charged
under the SM gauge group and the $G_{v}$ and those which solely have a $G_{v}$
charge. Interactions between the new sector and SM particles can occur through
intermediary mechanisms, such as a TeV scale $Z^{\prime}$ or through higher-
dimensional operators and loops involving heavy particles possessing non-
trivial $G_{\text{SM}}$ and $G_{v}$ charges. The case of a confined dark
sector decoupled from SM in the lower energies is of particular interest.

###  3.1 Simulation details

Strongly coupled hidden-sector particles produced at the LHC could produce
invisible dark mesons along with stable mesons. In such a case, novel LHC
signatures such as semi-visible jets could be produced, almost
indistinguishable from QCD backgrounds, therefore posing an exciting challenge
for jet classification. This study considers resonant and non-resonant
production of semi-visible jets. For the latter, we consider both Electroweak
(EW) and QCD productions separately. We simulate these dark shower jets events
at 14 TeV centre-of-mass energy proton-proton collisions using the Hidden
Valley (HV) module of Pythia8.310 Bierlich:2022pfr . They are generated as:

  * •

Resonant production: A new leptophobic vector boson $Z^{\prime}$ is produced
as a BSM mediator, decaying to a pair of dark quarks. Events are generated
setting with the process HiddenValley:ffbar2Zv = on. We set the mass of
$Z^{\prime}$ = 1 TeV with a width of 10 MeV. The $Z^{\prime}$ decays to a pair
of dark quarks ($q_{d}\,\bar{q}_{d})$, which has no SM charge but can decay to
another dark sector particle and an SM quark. The subsequent showering and
hadronisation produce dark mesons and SM particles.

  * •

Non-resonant EW production: We produce a pair of $U_{v}$ quarks using the
process HiddenValley:ffbar2UvUvbar = on. These $U_{v}$’s are charged under the
fundamental representation of the dark SU(N) and mirror the SM charges of the
$u$-quark. It, therefore, radiates into both SM and dark sector particles.

  * •

Non-resonant QCD production: The same $U_{v}$ pair of quarks are produced
using the process HiddenValley:gg2UvUvbar = on. Unlike the previous channel,
which occurs via an s-channel colour singlet propagator, this channel has an
SM colour connection between the initial and the final state lines, resembling
QCD dijet productions more closely.

The fraction of stable to total dark mesons produced in the process is
controlled by the probVector parameter in the HV Module. We have chosen
probVector=0.5 in our analysis. This choice is motivated by higher values
producing events with low missing transverse energy, which is covered by dijet
bump hunts. In comparison, lower values correspond to signatures where the
missing transverse energy’s $\phi$ direction is away from the jet axes and,
hence, are covered by more traditional dark matter searches.

All final state particles detectable at the LHC are clustered into microjets
of radius $R=0.1$ using the anti-$k_{t}$ algorithm Cacciari:2008gp , mimicking
the calorimeter resolution. Inclusive microjets with transverse momentum
$p_{T}>1$ GeV are then clustered into jets of radius $R=0.8$ with the
anti-$k_{t}$ algorithm. The jet clusterings are performed using the FastJet
(v3.4.1) Cacciari:2011ma  package. Events are required to have at least two
jets within $|\eta|<3$, the leading jet is required to have
$p_{T}>150~{}\text{GeV}$, and other sub-leading jets are required to have
$p_{T}>120~{}\text{GeV}$.

For the background, QCD dijet events are generated by setting HardQCD:all =
on. To have better efficiency in the event selection, we set
PhaseSpace:pTHatMin= 100, which puts a lower bound of 100 GeV on the
transverse momentum on each of the two final state legs in the $2\to 2$
processes for all classes. We generate 300k events passing the selection
criteria for each class, which are used to create three binary classification
datasets for each of the three signals vs the background. The leading jet’s
constituents in these events are used to construct jet graphs. We use
$R_{0}=0.5$ to construct the radius graph and extract
$\mathbf{h}^{(0)}_{i}=\Delta R_{iJ}$, and $\mathbf{x}^{(0)}_{i}=(\Delta
y_{iJ},\Delta\phi_{iJ})$ along with $z_{i}$ and
$\omega^{(\mathcal{N}[i])}_{j}$. For the ordinary EMPN architecture, we use
$\mathbf{H}^{(0)}_{i}=\mathbf{x}^{(0)}_{i}$. Since all three networks process
the node features to extract relevant edge features as inputs to the learnable
functions, no edge features are added. The combined datasets of the signal and
background (consisting of 600k samples) are segregated into 60% training, 20%
testing, and 20% validation datasets.

Figure 2: Representative Feynman diagrams for the dark shower event generation
processes.

###  3.2 Network Architecture and training

The network analysis is carried out using the PyTorch-Geometric fey2019fast
package. We consider an IRC-safe EdgeConv architecture similar to the one
studied in reference Konar:2021zdg  as an example of a non-equivariant
architecture. For ease of reading, we name the three architecture as
$\mathbf{E(2)}$-EMPN, $\mathbf{O(2)}$-EMPN, and EdgeConv-EMPN. As one of the
hallmarks of equivariant architectures is the efficiency of feature extraction
with a low number of tunable parameters, we consider a small variant
containing about 4k parameters and a large variant containing about 150k
parameters. The $\mathbf{E(2)}$-EMPN and $\mathbf{O(2)}$-EMPN have the same
structure of MLPs, and the difference arises from the edge input to
$\Phi^{(l)}_{e}$ and the aggregation of the vectors. Since the EdgeConv-EMPN
has only one MLP, while the equivariant ones have three, we reduce the
dimensionality of the hidden (scalar) node representations in the latter to
keep the number of parameters comparable. We apply the message-passing
operation twice for all architectures and extract the IRC-safe graph
representation after each operation to feed to a classifier MLP. A schematic
diagram for the equivariant ones is shown on the right side of figure 1.

For the small parameter size, the equivariant architectures have a
12-dimensional updated scalar node feature after each message operation, while
the EdgeConv-EMPN has 16-dimensional updated node features. For
$l\in\\{1,2\\}$, the MLPs $\Phi^{(l)}_{e}$, $\Phi^{(l)}_{x}$, and
$\Phi_{h}^{(l)}$ have two hidden layers with the same dimension as the updated
scalar feature dimension. Similarly, the edge MLPs have two hidden layers with
the same dimensions as the updated node features. An analogous architecture is
repeated for the large parameter case, with the equivariant ones updating
80-dimensional scalar node features and the EdgeConv-EMPN updating
128-dimensional ones. All hidden layers have ReLU activation. Except for
$\Phi_{x}^{(l)}$’s, which have a one-dimensional output with Sigmoid
activation, we do not put any output activation for the other MLPs. The graph
representation for the $\mathbf{E(2)}$-EMPN and $\mathbf{O(2)}$-EMPN cases is
obtained from the two scalar node representations and two vector node
representations.

In total, we have three binary classification datasets and six networks,
thereby giving us eighteen training cases. We train the network from random
initialisation five times for each of these cases with a batch size of 300 for
200 epochs. The networks are optimised with the Adam
DBLP:journals/corr/KingmaB14  optimiser with an initial learning rate of
$10^{-3}$. A decay-on-plateau criterion is utilised for the learning rate,
triggering a decay by a factor of 0.5 if the validation loss remains stagnant
for three consecutive epochs.

###  3.3 Results

Figure 3: The network performance is shown for two scenarios: one with a low
number of neural network parameters $(\approx 4k)$ and another with a high
number of parameters $(\approx 150k)$, for $\mathbf{E(2)}$-EMPN (solid),
$\mathbf{O(2)}$-EMPN (dashed) and EdgeConv-EMPN (dotted) for the three cases
of semi-visible jet tagger.

We extract the prediction on the test dataset from the epoch with the lowest
validation lost for each of the five training of the eighteen instances. To
obtain a summary representation of the receiver-operator characteristics (ROC)
curve, we extract the background acceptance $\epsilon_{B}$, at fixed signal
efficiency $\epsilon_{S}$ over the five training instances of training using a
custom implementation.333This was verified against the scikit-learn scikit-
learn  implementation with AUCs being equal up to four decimal points.
Extracting the mean value of the background acceptance $\bar{\epsilon}_{B}$,
over these five instances, we plot the ROC curve between
$1/\bar{\epsilon}_{B}$ and $\epsilon_{S}$ in figure 3. The average of the area
under the curve (AUC) and the standard deviation are listed in table 1. Before
comparing the relative differences for the considered architectures, we see
that the classification efficiency follows the usual physics intuition. The
resonant production of a heavy colour singlet $Z^{\prime}$ sets a hard scale,
which induces significant changes in the jets, while the non-resonant
productions do not have such a hard scale. In the colour-singlet s-channel
non-resonant production via the $\gamma$ and $Z$, there is no colour flow from
the initial quarks to the final state (even though the $U_{v}$ quark carries
an SM QCD charge). These jets are less QCD-like than the non-resonant
production via gluons, where colour flows between the initial and final states
and hence closely resembles QCD jet production.

NN parameters 4$k$ 150$k$ Process $g\,g\rightarrow U_{v}\,\bar{U}_{v}$
$q\,\bar{q}\rightarrow\gamma^{*}/Z^{*}\rightarrow U_{v}\,\bar{U}_{v}$
$q\,\bar{q}\rightarrow Z^{\prime}\rightarrow q_{d}\,\bar{q}_{d}$
$g\,g\rightarrow U_{v}\,\bar{U}_{v}$
$q\,\bar{q}\rightarrow\gamma^{*}/Z^{*}\rightarrow U_{v}\,\bar{U}_{v}$
$q\,\bar{q}\rightarrow Z^{\prime}\rightarrow q_{d}\,\bar{q}_{d}$
$\mathbf{E(2)}$-EMPN 0.807±0.051 0.896±0.003 0.913±0.002 0.855±0.002
0.906±0.001 0.923±0.001 $\mathbf{O(2)}$-EMPN 0.827±0.036 0.872±0.037
0.904±0.028 0.848±0.003 0.903±0.002 0.922±0.001 EdgeConv-EMPN 0.767±0.041
0.817±0.042 0.845±0.006 0.848±0.000 0.904±0.000 0.921±0.001

Table 1: The network performance is illustrated for two scenarios: one with a
low number of neural network parameters $(\approx 4k$) and another with a high
number of parameters $(\approx 150k$), for the $\mathbf{E(2)}$-EMPN,
$\mathbf{O(2)}$-EMPN, and EdgeConv-EMPN.

For the large network on the right side of figure 3, the three architectures
have virtually overlapping ROCs with a very minute decrease for the non-
equivariant EdgeConv-EMPN. The values of the AUCs in table 1 follow the same
trend. This suggests that with sufficiently large parametrisation, the
EdgeConv-EMPN can catch the distinguishing features as well as the
$\mathbf{E}(2)$ and $\mathbf{O}(2)$ equivariant ones, even though the
extracted features are not equivariant. In stark contrast, there is a very
noticeable dip in performance for the EdgeConv-EMPN of 4k parameters for all
three signal cases compared to the two equivariant architectures, with the
$\mathbf{E(2)}$-EMPN having the best classification power for the resonant and
electroweak non-resonant processes, while the $\mathbf{O(2)}$-EMPN bettering
the $\mathbf{E(2)}$-EMPN for the non-resonant production via gluons. Moreover,
the best equivariant AUCs in the low parameter case are not that far off from
their corresponding best ones in the large parameter cases, with the highest
difference for the hardest case of distinguishing gluon-induced production of
semi-visible jets. This clearly shows that the inductive bias of equivariance
in the rapidity-azimuth plane helps optimise feature extraction and reduces
the need for many parameters.

##  4 Summary and Conclusions

The future roadmap of particle physics crucially depends on maximising the
discovery potential of the LHC and its high luminosity phase. Recent
developments in optimising searches for new physics using machine learning
have demonstrated the enormous opportunity created by less traditional
approaches to data analysis. Interfacing such highly adapted strategies with
theoretical predictions to make robust statements about parameter
measurements, whether related to SM or BSM effects, highlights infrared and
collinear safety as a central theme of neural network architecture design,
e.g. through energy-weighted message passing. Adding such features to the
neural network design can practically increase computational demands, which
requires optimisation improvements. In this work, we have considered
$\mathbf{E(2)}$ equivariance as an avenue to remove computational overhead. In
particular, we have considered jet-substructure-based analysis to show that an
$\mathbf{E(2)}$ equivariant Graph Neural Network optimises the network’s
learning capability by directly exploiting the physical symmetry that the data
set possesses. As most jet substructure observables deal with rotationally
invariant quantities, we also considered $\mathbf{O(2)}$-equivariance which is
a subset of $\mathbf{E(2)}$.

Such advances directly impact physics applications. To highlight this, we have
considered a representative (but not exclusive) scenario of semi-visible jets.
In theories underpinning this BSM signature, extracting, e.g., dark shower
parameters critically depends on an IRC-safe implementation of the tagging
algorithm as the model navigates between visibly hard and invisibly soft final
states. This asks for the correct modelling of the QCD null hypothesis to draw
theoretically consistent conclusions. The architecture presented in this work
provides excellent signal vs. background discrimination with these desired
properties. This opens up the possibility of employing the architectures
proposed in this work as a new tool for BSM discovery in challenging QCD-
dominated final states.

### Acknowledgements

A.B. is supported by the U.S. Department of Energy under grant number DE-SC
0016013. C.E. is supported by the UK Science and Technology Facilities Council
(STFC) under grant ST/X000605/1 and the Leverhulme Trust under RPG-2021-031.
W.N. is funded by a University of Glasgow College of Science and Engineering
Scholarship. V. S. N. and M.S. are supported by the STFC under grant
ST/P001246/1. Parts of the computation utilised the Param Vikram-1000 High
Performance Computing Cluster at the Physical Research Laboratory.

## References

  * (1) J. Cogan, M. Kagan, E. Strauss and A. Schwarztman, _Jet-Images: Computer Vision Inspired Techniques for Jet Tagging_ , _JHEP_ 02 (2015) 118 [1407.5675].
  * (2) L. de Oliveira, M. Kagan, L. Mackey, B. Nachman and A. Schwartzman, _Jet-images — deep learning edition_ , _JHEP_ 07 (2016) 069 [1511.05190].
  * (3) G. Kasieczka, T. Plehn, M. Russell and T. Schell, _Deep-learning Top Taggers or The End of QCD?_ , _JHEP_ 05 (2017) 006 [1701.08784].
  * (4) G. Louppe, K. Cho, C. Becot and K. Cranmer, _QCD-Aware Recursive Neural Networks for Jet Physics_ , _JHEP_ 01 (2019) 057 [1702.00748].
  * (5) F. A. Dreyer, G. P. Salam and G. Soyez, _The Lund Jet Plane_ , _JHEP_ 12 (2018) 064 [1807.04758].
  * (6) T. Heimel, G. Kasieczka, T. Plehn and J. M. Thompson, _QCD or What?_ , _SciPost Phys._ 6 (2019) 030 [1808.08979].
  * (7) M. Farina, Y. Nakai and D. Shih, _Searching for New Physics with Deep Autoencoders_ , _Phys. Rev. D_ 101 (2020) 075021 [1808.08992].
  * (8) H. Qu and L. Gouskos, _ParticleNet: Jet Tagging via Particle Clouds_ , _Phys. Rev. D_ 101 (2020) 056019 [1902.08570].
  * (9) A. Butter et al., _The Machine Learning landscape of top taggers_ , _SciPost Phys._ 7 (2019) 014 [1902.09914].
  * (10) E. Bernreuther, T. Finke, F. Kahlhoefer, M. Krämer and A. Mück, _Casting a graph net to catch dark showers_ , _SciPost Phys._ 10 (2021) 046 [2006.08639].
  * (11) H. Qu, C. Li and S. Qian, _Particle transformer for jet tagging_ , in _Proceedings of the 39th International Conference on Machine Learning_ , K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu and S. Sabato, eds., vol. 162 of _Proceedings of Machine Learning Research_ , pp. 18281–18292, PMLR, 17–23 Jul, 2022, https://proceedings.mlr.press/v162/qu22b.html.
  * (12) G. Kasieczka, S. Marzani, G. Soyez and G. Stagnitto, _Towards Machine Learning Analytics for Jet Substructure_ , _JHEP_ 09 (2020) 195 [2007.04319].
  * (13) F. A. Dreyer and H. Qu, _Jet tagging in the Lund plane with graph networks_ , _JHEP_ 03 (2021) 052 [2012.08526].
  * (14) C. K. Khosa and S. Marzani, _Higgs boson tagging with the Lund jet plane_ , _Phys. Rev. D_ 104 (2021) 055043 [2105.03989].
  * (15) F. Canelli, A. de Cosa, L. L. Pottier, J. Niedziela, K. Pedro and M. Pierini, _Autoencoders for semivisible jet detection_ , _JHEP_ 02 (2022) 074 [2112.02864].
  * (16) F. A. Dreyer, G. Soyez and A. Takacs, _Quarks and gluons in the Lund plane_ , _JHEP_ 08 (2022) 177 [2112.09140].
  * (17) L. Cavallini, A. Coccaro, C. K. Khosa, G. Manco, S. Marzani, F. Parodi et al., _Tagging the Higgs boson decay to bottom quarks with colour-sensitive observables and the Lund jet plane_ , _Eur. Phys. J. C_ 82 (2022) 493 [2112.09650].
  * (18) A. Butter, B. M. Dillon, T. Plehn and L. Vogel, _Performance versus resilience in modern quark-gluon tagging_ , _SciPost Phys. Core_ 6 (2023) 085 [2212.10493].
  * (19) B. M. Dillon, L. Favaro, T. Plehn, P. Sorrenson and M. Krämer, _A normalized autoencoder for LHC triggers_ , _SciPost Phys. Core_ 6 (2023) 074 [2206.14225].
  * (20) T. Faucett, S.-C. Hsu and D. Whiteson, _Learning to identify semi-visible jets_ , _JHEP_ 12 (2022) 132 [2208.10062].
  * (21) P. Keicher, _Machine Learning in Top Physics in the ATLAS and CMS Collaborations_ , in _15th International Workshop on Top Quark Physics_ , 1, 2023, 2301.09534.
  * (22) T. Finke, M. Krämer, A. Mück and J. Tönshoff, _Learning the language of QCD jets with transformers_ , _JHEP_ 06 (2023) 184 [2303.07364].
  * (23) C.-T. Lu, H. Lv, W. Shen, L. Wu and J. Zhang, _Probing dark QCD sector through the Higgs portal with machine learning at the LHC_ , _JHEP_ 08 (2023) 187 [2304.03237].
  * (24) D. Bardhan, Y. Kats and N. Wunch, _Searching for dark jets with displaced vertices using weakly supervised machine learning_ , _Phys. Rev. D_ 108 (2023) 035036 [2305.04372].
  * (25) E. M. Metodiev, J. Thaler and R. Wynne, _Anomaly Detection in Collider Physics via Factorized Observables_ , 2312.00119.
  * (26) L. Favaro, M. Krämer, T. Modak, T. Plehn and J. Rüschkamp, _Semi-visible jets, energy-based models, and self-supervision_ , 2312.03067.
  * (27) A. Romero and D. Whiteson, _Jet Rotational Metrics_ , 2311.06686.
  * (28) A. Furuichi, S. H. Lim and M. M. Nojiri, _Jet Classification Using High-Level Features from Anatomy of Top Jets_ , 2312.11760.
  * (29) T. Gaertner and J. Reiten, _The simplicial substructure of jets_ , 2312.06596.
  * (30) A. Bogatskiy, B. Anderson, J. T. Offermann, M. Roussi, D. W. Miller and R. Kondor, _Lorentz Group Equivariant Neural Network for Particle Physics_ , 2006.04780.
  * (31) S. Gong, Q. Meng, J. Zhang, H. Qu, C. Li, S. Qian et al., _An efficient Lorentz equivariant graph neural network for jet tagging_ , _JHEP_ 07 (2022) 030 [2201.08187].
  * (32) C. Li, H. Qu, S. Qian, Q. Meng, S. Gong, J. Zhang et al., _Does Lorentz-symmetric design boost network performance in jet physics?_ , 2208.07814.
  * (33) Z. Hao, R. Kansal, J. Duarte and N. Chernyavskaya, _Lorentz group equivariant autoencoders_ , _Eur. Phys. J. C_ 83 (2023) 485 [2212.07347].
  * (34) A. Bogatskiy, T. Hoffman, D. W. Miller and J. T. Offermann, _PELICAN: Permutation Equivariant and Lorentz Invariant or Covariant Aggregator Network for Particle Physics_ , 2211.00454.
  * (35) S. Choi, S. J. Lee and M. Perelstein, _Infrared Safety of a Neural-Net Top Tagging Algorithm_ , _JHEP_ 02 (2019) 132 [1806.01263].
  * (36) P. T. Komiske, E. M. Metodiev and J. Thaler, _Energy Flow Networks: Deep Sets for Particle Jets_ , _JHEP_ 01 (2019) 121 [1810.05165].
  * (37) P. Konar, V. S. Ngairangbam and M. Spannowsky, _Energy-weighted message passing: an infra-red and collinear safe graph neural network algorithm_ , _JHEP_ 02 (2022) 060 [2109.14636].
  * (38) O. Atkinson, A. Bhardwaj, C. Englert, P. Konar, V. S. Ngairangbam and M. Spannowsky, _IRC-Safe Graph Autoencoder for Unsupervised Anomaly Detection_ , _Front. Artif. Intell._ 5 (2022) 943135 [2204.12231].
  * (39) D. Athanasakos, A. J. Larkoski, J. Mulligan, M. Płoskoń and F. Ringer, _Is infrared-collinear safe information all you need for jet classification?_ , 2305.08979.
  * (40) W. Shen, D. Wang and J. M. Yang, _Hierarchical high-point Energy Flow Network for jet tagging_ , _JHEP_ 09 (2023) 135 [2308.08300].
  * (41) P. Konar, V. S. Ngairangbam and M. Spannowsky, _Hypergraphs in LHC phenomenology — the next frontier of IRC-safe feature extraction_ , _JHEP_ 01 (2024) 113 [2309.17351].
  * (42) L. G. Almeida, S. J. Lee, G. Perez, G. F. Sterman, I. Sung and J. Virzi, _Substructure of high- $p_{T}$ Jets at the LHC_, _Phys. Rev. D_ 79 (2009) 074017 [0807.0234].
  * (43) J. Thaler and K. Van Tilburg, _Identifying Boosted Objects with N-subjettiness_ , _JHEP_ 03 (2011) 015 [1011.2268].
  * (44) J. Thaler and K. Van Tilburg, _Maximizing Boosted Top Identification by Minimizing N-subjettiness_ , _JHEP_ 02 (2012) 093 [1108.2701].
  * (45) A. J. Larkoski, G. P. Salam and J. Thaler, _Energy Correlation Functions for Jet Substructure_ , _JHEP_ 06 (2013) 108 [1305.0007].
  * (46) A. J. Larkoski, J. Thaler and W. J. Waalewijn, _Gaining (Mutual) Information about Quark/Gluon Discrimination_ , _JHEP_ 11 (2014) 129 [1408.3122].
  * (47) I. Moult, L. Necib and J. Thaler, _New Angles on Energy Correlation Functions_ , _JHEP_ 12 (2016) 153 [1609.07483].
  * (48) M. Dasgupta, L. Schunk and G. Soyez, _Jet shapes for boosted jet two-prong decays from first-principles_ , _JHEP_ 04 (2016) 166 [1512.00516].
  * (49) S. Caletti, O. Fedkevych, S. Marzani, D. Reichelt, S. Schumann, G. Soyez et al., _Jet angularities in Z+jet production at the LHC_ , _JHEP_ 07 (2021) 076 [2104.06920].
  * (50) T. Kinoshita, _Mass singularities of Feynman amplitudes_ , _J. Math. Phys._ 3 (1962) 650.
  * (51) T. D. Lee and M. Nauenberg, _Degenerate Systems and Mass Singularities_ , _Phys. Rev._ 133 (1964) B1549.
  * (52) M. J. Strassler and K. M. Zurek, _Echoes of a hidden valley at hadron colliders_ , _Phys. Lett. B_ 651 (2007) 374 [hep-ph/0604261].
  * (53) L. Carloni and T. Sjostrand, _Visible Effects of Invisible Hidden Valley Radiation_ , _JHEP_ 09 (2010) 105 [1006.2911].
  * (54) L. Carloni, J. Rathsman and T. Sjostrand, _Discerning Secluded Sector gauge structures_ , _JHEP_ 04 (2011) 091 [1102.3795].
  * (55) T. Cohen, M. Lisanti and H. K. Lou, _Semivisible Jets: Dark Matter Undercover at the LHC_ , _Phys. Rev. Lett._ 115 (2015) 171804 [1503.00009].
  * (56) D. Kar and S. Sinha, _Exploring jet substructure in semi-visible jets_ , _SciPost Phys._ 10 (2021) 084 [2007.11597].
  * (57) E. Bernreuther, J. C. Mejia, F. Kahlhoefer, M. Krämer and P. Tunney, _On the challenges of searching for GeV-scale long-lived particles at the LHC_ , _JHEP_ 04 (2021) 210 [2011.06604].
  * (58) J. Barron, D. Curtin, G. Kasieczka, T. Plehn and A. Spourdalakis, _Unsupervised hadronic SUEP at the LHC_ , _JHEP_ 12 (2021) 129 [2107.12379].
  * (59) G. Albouy et al., _Theory, phenomenology, and experimental avenues for dark showers: a Snowmass 2021 report_ , _Eur. Phys. J. C_ 82 (2022) 1132 [2203.09503].
  * (60) Y. Wang, Y. Sun, Z. Liu, S. E. Sarma, M. M. Bronstein and J. M. Solomon, _Dynamic graph CNN for learning on point clouds_ , _CoRR_ abs/1801.07829 (2018) [1801.07829].
  * (61) L.-H. Lim and B. J. Nelson, _What is… an equivariant neural network?_ , _Notices of the American Mathematical Society_ 70 (2023) 1.
  * (62) V. G. Satorras, E. Hoogeboom and M. Welling, _E(n) equivariant graph neural networks_ , in _Proceedings of the 38th International Conference on Machine Learning_ , M. Meila and T. Zhang, eds., vol. 139 of _Proceedings of Machine Learning Research_ , pp. 9323–9332, PMLR, 18–24 Jul, 2021, https://proceedings.mlr.press/v139/satorras21a.html.
  * (63) J. Gallicchio and M. D. Schwartz, _Seeing in Color: Jet Superstructure_ , _Phys. Rev. Lett._ 105 (2010) 022001 [1001.5027].
  * (64) ATLAS collaboration, _Measurement of colour flow using jet-pull observables in $t\bar{t}$ events with the ATLAS experiment at $\sqrt{s}=13\,\hbox{TeV}$_, _Eur. Phys. J. C_ 78 (2018) 847 [1805.02935].
  * (65) A. Larkoski, S. Marzani and C. Wu, _Safe Use of Jet Pull_ , _JHEP_ 01 (2020) 104 [1911.05090].
  * (66) CMS collaboration, _Search for new particles in events with energetic jets and large missing transverse momentum in proton-proton collisions at $\sqrt{s}$ = 13 TeV_, _JHEP_ 11 (2021) 153 [2107.13021].
  * (67) ATLAS collaboration, _Search for new phenomena in final states with photons, jets and missing transverse momentum in pp collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector_, _JHEP_ 07 (2023) 021 [2206.06012].
  * (68) C. Bierlich et al., _A comprehensive guide to the physics and usage of PYTHIA 8.3_ , _SciPost Phys. Codeb._ 2022 (2022) 8 [2203.11601].
  * (69) M. Cacciari, G. P. Salam and G. Soyez, _The anti- $k_{t}$ jet clustering algorithm_, _JHEP_ 04 (2008) 063 [0802.1189].
  * (70) M. Cacciari, G. P. Salam and G. Soyez, _FastJet User Manual_ , _Eur. Phys. J. C_ 72 (2012) 1896 [1111.6097].
  * (71) M. Fey and J. E. Lenssen, _Fast graph representation learning with pytorch geometric_ , 2019.
  * (72) D. P. Kingma and J. Ba, _Adam: A method for stochastic optimization_ , in _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , Y. Bengio and Y. LeCun, eds., 2015, http://arxiv.org/abs/1412.6980.
  * (73) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel et al., _Scikit-learn: Machine learning in Python_ , _Journal of Machine Learning Research_ 12 (2011) 2825.


###### Abstract

We investigate experimentally and analytically the coalescence of
reflectionless (RL) states in symmetric complex wave-scattering systems. We
observe RL-exceptional points (EPs), first, with a conventional Fabry-Perot
system for which the scattering strength within the system is tuned
symmetrically, and then with single- and multi-channel symmetric disordered
systems. We identify that an EP of the parity-time (PT)-symmetric RL-operator
is obtained when the spacing between central frequencies of two natural
resonances of the system is equal to the decay rate into incoming and outgoing
channels. Finally, we leverage the transfer functions associated with RL and
RL-EPs states to implement first- and second-order analog differentiation.

††preprint: APS/123-QED

Exceptional points (EP) are spectral singularities in non-Hermitian systems at
which two or more eigenvalues and eigenstates coalesce [1, 2, 3, 4, 5]. EPs
have mainly been explored for resonant states that are the poles of the
scattering matrix and scattering states. In systems with losses only, finding
an EP between two resonant modes requires that the mutual coupling between
resonances satisfy a critical relation with their loss factors [6]. EPs are
however more easily realizable in systems with PT-symmetry obtained by
balancing gain and loss in symmetric regions [7, 6]. Many unconventional
features of EPs have been demonstrated theoretically and experimentally such
as the design of unidirectional invisibility [8, 9], asymmetric mode switching
[10] or directional lasing [11], to cite a few (see Refs. [12, 6] for
reviews). These spectral singularities are also interesting for sensing
applications since the energy (or frequency) splitting between $n$ degenerate
eigenstates scales as the $n$th root of the perturbation [13, 14].

Recently, a new kind of EPs associated with reflectionless (RL) scattering
states rather than resonances has been investigated. RL states are eigenstates
of a non-Hermitian operator $H_{\mathrm{RL}}$ based on the wave equation with
incoming channels connected to a scattering system modeled as gain and
outgoing channels modeled as losses [15, 16, 17]. The eigenvalues
$\tilde{\omega}_{R}$ of $H_{\mathrm{RL}}$ are distinct from the resonance
spectra related to the poles of the scattering matrix $S(\omega)$. When an
eigenvalue is tuned to the real axis, the corresponding RL-state enable
reflectionless coupling of incoming channels. The excitation of multi-channel
RL-states in disordered matter requires non-trivial wavefront shaping [18].

A special case of RL-EPs found for purely incoming boundary conditions is the
coalescence on the real axis of two zeros of the complete scattering matrix
$S(\omega)$ known as a perfectly absorbing EP [19, 20, 21]. The phenomenon of
coherent perfect absorption (CPA) occurs when absorption within a scattering
medium balances the excitation rate of incoming channels [22, 23, 24, 25]. CPA
is the time-reverse of lasing at threshold [26, 27] and a generalization of
the critical coupling condition [28, 24]. CPA has been implemented in a wide
range of regular scattering systems [29], as well as in disordered matter [30,
31, 32, 33, 34, 35]. The latter has been achieved in particular through
purposeful perturbations of the scattering system, which can include
additional constraints, such as on frequency or CPA wavefront [30, 32, 33,
34]. When two CPA-states coalesce, the broadened lineshape of the absorption
spectrum indicates the existence of perfectly absorbing EPs [19, 20, 21].

In flux-conserving systems, zero reflection and therefore perfect transmission
of an incoming wavefront occurs when the decay rates into incoming and
outgoing channels are equal. The probability of finding RL scattering modes is
naturally enhanced in scattering systems with mirror-symmetry for which the
RL-operator is PT-symmetric and RL-eigenvalues feature an exciting property.
An RL-eigenvalue remains on the real axis under continuous perturbation of the
system until it coalesces with a second one at an EP before splitting into a
complex-conjugate pair [15, 19, 17]. The spectral broadening of the
transmission peak corresponding to RL-EPs has been observed in multimirror
cavities [36, 37, 38] and in numerical simulations of quantum dots [39] and
atomic wires [40], even though not interpreted as an EP. However, a clear
experimental demonstration of RL-EPs in disordered systems and an analysis of
RL-states in term of natural resonances has to date not been reported.

Figure 1: (a) Photography of the experimental setup which is a rectangular
waveguide supporting only the fundamental transverse-electric mode between 7
and 11 GHz. The penetration depth $z$ of two rods inserted symmetrically near
the middle of the waveguide can be finely controlled with a translation stage.
Inset: a channel (coax-to-waveguide transition) with an empty opening or with
an alumina slab placed in front of it to increase the reflectivity at the
interfaces of the waveguide. (b-c) Colorscale representation of the reflection
on a dB-scale, $10\mathrm{log}[R(\nu,z)]$, measured with a vector network
analyzer between 8 and 10 GHz for a symmetric (b) and a asymmetric (c)
perturbation. RL-EPs appear only in the symmetric case when two branches with
$R\sim 0$ shown in yellow coincide. (d-f) $R(\nu,z)$ in the underperturbed
$z<z_{\mathrm{EP}}$ (d), critically perturbed $z=z_{\mathrm{EP}}$ (e) and
overperturbed $z>z_{\mathrm{EP}}$ (f) regimes as a function of rescaled
frequency $\nu-\nu_{0}$ ($\nu_{0}=8.2$ GHz). The spectra are shown for zero
(green line), one (red line) and two (blue line) dielectric slabs at each
interface. The corresponding RL-eigenvalues found from Eq. (2) are shown on
the top in the complex plane.

In this article, we observe the existence of RL-EPs in symmetric regular and
disordered systems that are tuned toward this scattering anomaly by inserting
a symmetric defect. We first start with a multimirror Fabry-Perot cavity for
which the reflectivity at the center is progressively increased. Our harmonic
analysis reveals that an RL-EP is obtained when the spacing between the
central frequencies of two quasi-normal modes (QNMs) is equal to their
linewidth. Our experimental results in the microwave range are in excellent
agreement with a two-level model. We then demonstrate RL-EPs in disordered
single- and multi-channel waveguides. Finally, we show that designing systems
operating at an RL-EP is relevant to analog computations of derivatives.

We measure spectra of the reflection coefficient $r(\nu)$ and the transmission
coefficient $t(\nu)$ through a single-mode rectangular waveguide (length
$L=400$ mm, width $W=22.86$ mm and height $H=10.16$ mm) with two coax-to-
waveguide transitions attached to the openings (see Fig. 1(a)). A dielectric
alumina slab of reflectivity $R_{0}=0.25$ (see SM) is positioned between each
transition and the waveguide to increase the internal reflection at the
interfaces. A perturbation is introduced symmetrically with respect to the
center of the waveguide by inserting two aluminium rods (diameter $2$ mm)
through two holes drilled into the top plate and spaced by $12$ mm. The
penetration depth $z$ varies from 0 to $8$ mm in steps of $\Delta z=0.02$ mm
($\Delta z\sim\lambda/1666$). This system is therefore equivalent to a multi-
mirror Fabry-Perot interferometer with tunable reflectivity at the center.

A colormap of the reflection $R(\nu,z)=|r(\nu,z)|^{2}$ is shown in Fig. 1(b)
on a log-scale. In absence of the perturbation ($z=0$), the frequencies
corresponding to extremely small reflection (bright area) are regularly spaced
as expected for a Fabry-Perot interferometer. As the penetration depth $z$
increases, the RL-frequencies move closer until they coalesce for a critical
perturbation ($z=z_{\mathrm{EP}}$) to form RL-EPs. The two peaks on the
spectrum of $R(\nu)$ found in the underperturbed regime ($z<z_{\mathrm{EP}}$)
transform into a single broadband one at $z_{\mathrm{EP}}$ with a flattened
quartic line shape which is characteristic of EPs [19, 20, 41, 21]. Once two
RL-eigenvalues have collapsed, they leave the real axis as complex-conjugate
pairs. In the overperturbed regime $z>z_{\mathrm{EP}}$, the minimum of
$R(\nu)$ increases with $z$ as the imaginary part of RL-eigenvalues moves away
from the real axis. In contrast, for an asymmetrical perturbation of the
Fabry-Perot cavity (a single rod inserted at $x=L/4$), the zero-reflection
frequencies do not collapse but move independently in the complex plane (see
Fig. 1(c)).

The same procedure is then repeated for samples with zero and two alumina
slabs at each interface, yielding $R_{0}=0.002$ and $R_{0}=0.46$. For the
smallest $R_{0}$, the reflectivity is solely due to the small impedance
mismatch between the transitions and the waveguide. In each case, a symmetric
perturbation leads to the formation of RL-EPs. The linewidth of resonances
decreases with increasing $R_{0}$ and the spectral dips corresponding to RL-
states narrow as seen in Fig. 1(d-f).

Figure 2: (a) Modal analysis of the transmission coefficient
$t(\nu,z_{\text{EP}})$. Two modes with lines shown in orange and cyan mainly
contribute to $T=|t(\nu,z_{\text{EP}})|^{2}$ (blue line) which is perfectly
reconstructed from harmonic inversion (red dashed line). (b) the DOS present
two peaks corresponding to the maxima of the two resonances. (c,d) Real (c)
and Imaginary (d) parts of the reflectionless eigenvalues (red circles) and
eigenfrequencies (blue crosses) with respect to $z$. The black lines are the
theoretical values of $\tilde{\omega}_{R}$ found from Eq. (2) in which the
values of $\omega_{n}(z)$ and $\gamma(z)$ are injected.

We now analyze RL-eigenvalues found by tracking the local minima of
$R(\omega,z)$ in terms of natural resonances of the medium. The QNMs are
solutions of the wave equation with purely outgoing boundary conditions (i.e.
for both left and right channels). The eigenstates $\psi_{n}(x)$ are
associated with complex frequencies
$\tilde{\omega}_{n}=\omega_{n}-i\Gamma_{n}/2$ with central frequency
$\omega_{n}$ and linewidth $\Gamma_{n}$. We extract the set of complex
frequencies $\tilde{\omega}_{n}$ from a modal analysis of transmission spectra
between 7 and 11 GHz using the harmonic inversion method [42, 2] (see SM). In
absence of the perturbation, the central frequencies almost coincide with RL-
frequencies. However, for $z=z_{EP}$, two overlapping resonances mainly
contribute to the flattened transmission peak as seen in Fig. 2(a). Two peaks
are observed on the spectrum of the density of states (DOS)
$\rho=(1/2\pi)\Sigma_{n}(\Gamma_{n}/2)/[(\omega-\omega_{n})^{2}+(\Gamma_{n}/2)^{2}]$
in Fig. 2(b) with a spacing between central frequencies approximately given by
the linewidths. Even for the highest penetration depth of the rods, the
resonances do not coalesce.

We now derive analytical expressions of $\tilde{\omega}_{n}$ and RL-
eigenvalues $\tilde{\omega}_{\mathrm{R}}$ using a two-level model. The
effective Hamiltonian $H_{\mathrm{eff}}$ of the open system is expressed in
the basis of two successive modes of a closed one-dimensional waveguide,
$H_{\text{eff}}=H_{0}-iV_{0}V_{0}^{T}/2-iV_{1}V_{1}^{T}/2$, with $H_{0}$ being
a diagonal matrix with elements $\omega_{0}-\delta\omega_{0}/2$ and
$\omega_{0}+\delta\omega_{0}/2$. The left and right vectors $V_{0}$ and
$V_{1}$ account for the coupling of the closed system to the channels and are
real for systems with time-reversal symmetry. The scattering matrix
$S(\omega)$ writes
$S(\omega)=1-iV^{T}[\omega\mathbb{1}-H_{\text{eff}}]^{-1}V$, where
$V=[V_{0}\;V_{1}]$. For an empty one-dimensional cavity of length $L$, the
$n$th eigenfunction of the closed systems is
$\psi_{n}(x)=(2/L)\mathrm{sin}(n\pi x/L)$. We consider here an even resonance
($n\equiv 0\;(\bmod\;2)$) and an odd resonance ($n\equiv 1\;(\bmod\;2)$). The
coupling vectors are related to the derivative of $\psi_{n}(x)$ at $x=0$ and
$x=L$. $V_{0}$ is therefore symmetric, $V_{0}=\sqrt{\gamma/2}(1~{}1)$ and
$V_{1}$ is anti-symmetric $V_{1}=\sqrt{\gamma/2}(1~{}-1)$. The coupling rate
$\gamma/2$ of each channel to the cavity depends on the reflectivity at the
interfaces of the waveguide. The effective Hamiltonian is therefore a diagonal
matrix with eigenvalues
$\tilde{\omega}_{M\pm}=\omega_{0}\pm\delta\omega_{0}/2-i\gamma/2$ (see SM).

Inserting a rod in the middle of the waveguide leads to a local change
$\Delta\epsilon(x=L/2)$ of the permittivity. For a small penetration depth,
the frequency shift of the two resonances can be calculated using standard
perturbation theory. At first order, the shift is
$\Delta\omega_{n}\propto-\Delta\epsilon(x)|\psi_{n}(x)|^{2}$ [44, 45]. For an
even resonance, the field cancels at the center of the waveguide and
$\Delta\omega_{n}\sim 0$. For odd resonances, $|\psi_{n}(x)|^{2}\neq 0$ and
the central frequency shifts towards smaller values. As the penetration depth
increases, resonances come together in pairs with central frequencies and
spacings that are functions of $z$, $\omega_{0}(z)$ and $\delta\omega_{0}(z)$.

RL-modes are eigenvalues of the PT-symmetric operator
$H_{\mathrm{R}}=H_{0}+iV_{0}V_{0}^{T}/2-iV_{1}V_{1}^{T}/2$, where an effective
gain is associated to incoming channels at the left interface. The coupling of
the channels now results in anti-diagonal terms $i\gamma$:

$H_{\mathrm{R}}=\omega_{0}(z)\mathbb{1}+\frac{1}{2}\begin{pmatrix}-\delta\omega_{0}(z)&i\gamma\\\
i\gamma&\delta\omega_{0}(z)\end{pmatrix}$ (1)

The RL-eigenvalues $\tilde{\omega}_{R}$ are then found from a diagonalization
of $H_{\mathrm{RL}}$:

$\tilde{\omega}_{R\pm}=\omega_{0}(z)\pm\frac{1}{2}\sqrt{\delta\omega_{0}(z)^{2}-\gamma^{2}}.$
(2)

Uniform absorption within the waveguide is incorporated by adding an imaginary
part $-i\gamma_{a}$ to complex frequencies. We estimate from the modal
analysis that $\gamma_{a}=30$ MHz (see SM).

When the coupling of channels is small compared to the spacing between
resonances, $\gamma\ll\delta\omega_{0}$, RL-eigenvalues are real and coincide
with the central frequencies of resonances $\omega_{n}$. However, the spacing
between RL-eigenvalues decreases more rapidly than $\delta\omega_{0}(z)$ as
the perturbation strength is symmetrically tuned. A reflectionless EP is found
when losses through channels are equal to the spacing between the two
resonances, $\delta\omega_{0}(z_{\text{EP}})=\gamma$. For larger
perturbations, the two RL-eigenvalues become a complex-conjugate pair.
Assuming that $\delta\omega_{0}(z)$ scales linearly with the penetration depth
$z$ near the EP, $\delta\omega_{0}(z)=\gamma+\kappa/2(z_{\text{EP}}-z)$ (see
SM), gives a splitting between RL-eigenvalues for $z<z_{\text{EP}}$ near the
EP
$\tilde{\omega}_{R+}-\tilde{\omega}_{R-}\sim\sqrt{\gamma\kappa(z_{\text{EP}}-z)}$
which is characteristic of the square-root detuning behavior of EPs under
small perturbations. The same splitting is found on imaginary parts of
$\tilde{\omega}_{R}$ for $z>z_{EP}$.

The experimental results in Fig. 2(c) for $\text{Re}[\tilde{\omega}_{R}]$ are
in excellent agreement with the prediction of Eq. (2) in which we inject the
values of $\omega_{0}(z)$, $\gamma(z)$ and $\delta\omega_{0}(z)$ extracted
from the modal analysis. This agreement also highlights the effectiveness of
the two-level model. For $z<z_{\mathrm{EP}}$, the imaginary part of
$\tilde{\omega}_{R}$ is equal to the absorption decay rate
$\Gamma_{R}=-\gamma_{a}$. For $z>z_{\mathrm{EP}}$, we fit
$R(\omega=\omega_{0})$ using the following analytical expression for the
reflection coefficient (see SM for derivation):

$\begin{split}r&=\frac{[\omega-\omega_{0}(z)]^{2}+[\gamma^{2}-\delta\omega_{0}(z)^{2}]/4}{[\omega-\omega_{0}(z)]^{2}-[\delta\omega_{0}(z)^{2}+\gamma^{2}]/4+i\gamma[(\omega-\omega_{0}(z)]}\\\
&=\frac{(\omega-\tilde{\omega}_{R+})(\omega-\tilde{\omega}_{R-})}{(\omega-\tilde{\omega}_{M+})(\omega-\tilde{\omega}_{M-})}.\end{split}$
(3)

and obtain an excellent agreement for $\text{Im}[\tilde{\omega}_{R}]$ in Fig.
2(d).

A signature of an EP is the flattening of $R(\omega)$ close to
$\omega=\omega_{0}(z_{\text{EP}})$ [20, 21]. Eq. (S9) demonstrates that
$R(\omega)$ scales as
$R(\omega)\sim[\omega-\omega_{0}(z_{\text{EP}})]^{4}/\gamma^{4}$ at an EP and
therefore features a quartic line shape. This is confirmed in Fig. 1(e) in
reflection and Fig. 2(a) in transmission. We emphasize that $\gamma$ decreases
with increasing internal reflection so that the line shape is narrower for the
sample with two dielectric slabs.

Figure 3: (a) Colorscale representation of $R(\nu,z)$ for a symmetric single-
channel disordered medium. The inset shows a picture of metallic spheres
inserted within the waveguide. (b) Spectrum of reflection eigenvalues in a
multichannel disordered cavity in the underperturbed regime (left) and tuned
at an RL-EP (right). The arrows indicate the frequencies of the coalescing
zeros of $\tau_{4}(\nu)$. At the EP, the smallest eigenvalue displays the
expected broadband quartic line shape.

RL-EPs are not restricted to regular systems but exist in any complex
multiple-scattering system. We now add a symmetric disorder made of 28
metallic spheres of diameter 5 mm within the empty waveguide (see Fig. 3(a)).
The spacing between two spheres on the left side is drawn from a uniform
random distribution and we replicate this disorder with a mirror symmetry on
the right side. The spectrum $R(\nu)$ still presents clear dips corresponding
to RL-states but the spacing between two dips is now random. As the
perturbation is symmetrically inserted, pairs of these dips collapse to give
rise to RL-EPs. $z_{\text{EP}}$ is also a random variable which reflects the
random field distribution within the waveguide. In contrast to the Fabry-Perot
cavity case, however, not all RL-eigenvalues are on the real axis for $z=0$.
An interesting case is the variation of RL-states around 8.73 GHz. Two
eigenvalues first split at an EP for $z=1$ mm, move away along the real axis
and then come back together to form a second EP for $z=3.8$ mm.

We also realize an RL-EP in a symmetric multichannel systems. Two arrays of
$N=4$ transitions operating between 11 and 16 GHz are attached to an
effectively two-dimensional rectangular cavity with metallic boundary
conditions of length $L=0.5$ m, width $W=0.25$ m and height $h=8$ mm (see
Refs. [46, 47] and SM for details). Spectra of the $N\times N$ reflection
matrices $r(\nu)$ are measured both at the left and right interfaces of the
cavity. The cavity is made disordered with 6 aluminum rods symmetrically
placed on each side of the cavity. By tuning the penetration depth of two rods
in the same way as for the single channel case, we identify an RL-EP on the
last eigenvalue $\tau_{4}(\nu,z)$ of $r^{\dagger}r$ at 11.15 GHz. The two
zeros of $\tau_{4}(\nu,z)$ coalesce at the RL-EP with a characteristic quartic
line shape (see Fig. 3(b) and SM for the colorscale representation). We have
verified that the same RL-EP is found on the reflection matrix at the other
side of the cavity.

The incident wavefront $v_{R}$ corresponding with an RL-state is non-trivial
in complex scattering systems and can be given by the eigenvector of
$r^{\dagger}r$ with the smallest (near-zero) eigenvalue. When the spacing
$\delta$ between two zeros is large, the corresponding eigenvectors are
independent wavefronts. The degree of correlation between eigenvectors
$C=|v_{R+}v_{R-}^{\dagger}|$ is equal to 0.55 for $\delta=0.03$ GHz. As we
approach an RL-EP, the eigenvectors become strongly correlated with
$C\rightarrow 1$ as the two RL-states coalesce. A single eigenvalue
$\tau_{4}(\nu,z)$ with a flattened shape is close to zero and no decrease is
observed on $\tau_{3}(\nu,z)$. Note that the case of diabolical point with
non-degenerate states [4] would be dramatically different. The orthogonality
of the two states within the sample would lead to two reflection eigenvalues
$\tau_{3}$ and $\tau_{4}$ being simultaneously close to zero with independent
eigenvectors.

Figure 4: Analog differentiation of first (a) and second (b) order based on
the transfer function in reflection associated with a RL state and a RL-EP
state, respectively, of the 1D disordered waveguide. The figures display the
output signal envelopes calculated for the injection of a Gaussian pulse (45
MHz bandwidth). The carrier frequencies are 8.72 GHz (a) and 10.04 GHz (b)
with $z=5.9$ mm (a) and $z=4.49$ mm (b), respectively.

Finally, we note that the transfer functions $r(\nu)$ associated with an RL
state and an RL-EP state in reflection coincides with those of a first- and
second-order differentiator, respectively. Wave-based signal processing holds
the promise of being fast and energy-efficient. Implementations of higher-
order derivatives have been proposed mainly in carefully engineered static
optical fiber systems [48, 49, 50, 51, 52]. Recently, Ref. [41] proposed to
leverage the unprecedented flexibility of complex scattering systems tuned to
CPA in order to perform meta-programmable analog differentiation and
implemented higher-order differentiation by cascading two systems tuned to CPA
in order to emulate the transfer function of a CPA-EP. Here, by controlling
the depth of the perturber penetration, we can tune our system to either an RL
or RL-EP state. Thereby, we implement the second-order differentiator transfer
function without using any non-linear components. For the prototypical example
of a Gaussian input pulse, we calculate the envelopes of the output signal
based on the measured transfer function for these two states. The results are
displayed in Fig. 4 and in agreement with the analytically expected
derivative. We attribute the imperfections to the fact that absorption
prevents our RL-EP from lying directly on the real axis.

In conclusion, we have observed the coalescence of RL-states into EPs in
regular and disordered single- and multi-channel scattering systems. We have
shown that an RL-EP in a symmetric medium requires that the spacing between
two resonances is perfectly balanced by the coupling rate of symmetric and
anti-symmetric modes. An interesting extension of our work would consist in
observing EPs related to transmissionless states [53]. These states can be
found in two and three-dimensional scattering systems which theoretically
provide the same behavior (broadened lineshape of the transmission spectrum
and enhanced sensitivity of the zeros under perturbation of the system).

##  I Acknowledgments

This publication was supported by the European Union through the European
Regional Development Fund (ERDF), by the French region of Brittany and Rennes
Métropole through the CPER Project SOPHIE/STIC & Ondes. M. D. acknowledges the
Institut Universitaire de France. C. F. acknowledges funding from the French
“Ministère de la Défense, Direction Générale de l’Armement”.

##  II Data availability

The data that support the plots within this paper and other findings of this
study are available from the corresponding authors on reasonable request.

## References

  * Heiss [2000] W. Heiss, Phys. Rev. E 61, 929 (2000).
  * Berry [2004] M. V. Berry, Czechoslov. J. Phys. 54, 1039 (2004).
  * Rotter [2009] I. Rotter, J. Phys. A 42, 153001 (2009).
  * Moiseyev [2011] N. Moiseyev,  _Non-Hermitian quantum mechanics_ (Cambridge University Press, 2011).
  * Cao and Wiersig [2015] H. Cao and J. Wiersig, Rev. Mod. Phys. 87, 61 (2015).
  * Miri and Alù [2019] M.-A. Miri and A. Alu, Science 363, eaar7709 (2019).
  * Özdemir  _et al._ [2019] Å. K. Özdemir, S. Rotter, F. Nori, and L. Yang, Nat. Mater. 18, 783 (2019).
  * Lin  _et al._ [2011] Z. Lin, H. Ramezani, T. Eichelkraut, T. Kottos, H. Cao, and D. N. Christodoulides, Phys. Rev. Lett. 106, 213901 (2011).
  * Fleury  _et al._ [2015] R. Fleury, D. Sounas, and A. Alù, Nat. Comm. 6, 5905 (2015).
  * Doppler  _et al._ [2016] J. Doppler, A. A. Mailybaev, J. Böhm, U. Kuhl, A. Girschik, F. Libisch, T. J. Milburn, P. Rabl, N. Moiseyev, and S. Rotter, Nature 537, 76 (2016).
  * Peng  _et al._ [2016] B. Peng, Å. K. Özdemir, M. Liertzer, W. Chen, J. Kramer, H. Yilmaz, J. Wiersig, S. Rotter, and L. Yang, Proc. Natl. Acad. Sci. 113, 6845 (2016).
  * El-Ganainy  _et al._ [2018] R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Nat. Phys. 14, 11 (2018).
  * Chen  _et al._ [2017] W. Chen, Å. Kaya Özdemir, G. Zhao, J. Wiersig, and L. Yang, Nature 548, 192 (2017).
  * Hodaei  _et al._ [2017] H. Hodaei, A. U. Hassan, S. Wittek, H. Garcia-Gracia, R. El-Ganainy, D. N. Christodoulides, and M. Khajavikhan, Nature 548, 187 (2017).
  * Bonnet-Ben Dhia  _et al._ [2018] A.-S. Bonnet-Ben Dhia, L. Chesnel, and V. Pagneux, Proc. R. Soc. A 474, 20180050 (2018).
  * Sweeney  _et al._ [2020] W. R. Sweeney, C. W. Hsu, and A. D. Stone, Phys. Rev. A 102, 063511 (2020).
  * Stone  _et al._ [2020] A. D. Stone, W. R. Sweeney, C. W. Hsu, K. Wisal, and Z. Wang, Nanophotonics 1 (2020).
  * Pichler  _et al._ [2019] K. Pichler, M. Kuhmayer, J. Böhm, A. Brandstatter, P. Ambichl, U. Kuhl, and S. Rotter, Nature 567, 351 (2019).
  * Sweeney  _et al._ [2019] W. R. Sweeney, C. W. Hsu, S. Rotter, and A. D. Stone, Phys. Rev. Lett. 122, 093901 (2019).
  * Suwunnarat  _et al._ [2021] S. Suwunnarat, Y. Tang, M. Reisner, F. Mortessagne, U. Kuhl, and T. Kottos, arXiv:2103.03668 (2021).
  * Wang  _et al._ [2021] C. Wang, W. R. Sweeney, A. D. Stone, and L. Yang, Science 373, 1261 (2021).
  * Li  _et al._ [2017] H. Li, S. Suwunnarat, R. Fleischmann, H. Schanz, and T. Kottos, Phys. Rev. Lett. 118, 044101 (2017).
  * Fyodorov  _et al._ [2017] Y. V. Fyodorov, S. Suwunnarat, and T. Kottos, J. Phys. A Math 50, 30LT01 (2017).
  * Baranov  _et al._ [2017] D. G. Baranov, A. Krasnok, T. Shegai, A. Alù, and Y. Chong, Nat. Rev. Mater. 2, 17064 (2017).
  * Krasnok  _et al._ [2019] A. Krasnok, D. Baranov, H. Li, M.-A. Miri, F. Monticone, and A. Alù, Adv. Opt. Photonics 11, 892 (2019).
  * Chong  _et al._ [2010] Y. D. Chong, L. Ge, H. Cao, and A. D. Stone, Phys. Rev. Lett. 105, 053901 (2010).
  * Roger  _et al._ [2015] T. Roger, S. Vezzoli, E. Bolduc, J. Valente, J. J. F. Heitz, J. Jeffers, C. Soci, J. Leach, C. Couteau, N. I. Zheludev, and D. Faccio, Nat. Commun. 6, 7031 (2015).
  * Piper  _et al._ [2014] J. R. Piper, V. Liu, and S. Fan, Appl. Phys. Lett. 104, 251110 (2014).
  * Wan  _et al._ [2011] W. Wan, Y. Chong, L. Ge, H. Noh, A. D. Stone, and H. Cao, Science 331, 889 (2011).
  * F. Imani  _et al._ [2020] M. F. Imani, D. R. Smith, and P. del Hougne, Adv. Funct. Mater. 30, 2005310 (2020).
  * Chen  _et al._ [2020] L. Chen, T. Kottos, and S. M. Anlage, Nat. Commun. 11, 1 (2020).
  * Frazier  _et al._ [2020] B. W. Frazier, T. M. Antonsen Jr, S. M. Anlage, and E. Ott, Phys. Rev. Research 2, 043422 (2020).
  * del Hougne  _et al._ [2021] P. del Hougne, K. B. Yeo, P. Besnier, and M. Davy, Laser Photonics Rev. 15, 2000471 (2021).
  * Del Hougne  _et al._ [2021] P. Del Hougne, K. B. Yeo, P. Besnier, and M. Davy, Phys. Rev. Lett. 126, 193903 (2021).
  * Chen  _et al._ [2021] L. Chen, S. M. Anlage, and Y. V. Fyodorov, Phys. Rev. E 103, L050203 (2021).
  * Van de Stadt and Muller [1985] H. Van de Stadt and J. M. Muller, J. Opt. Soc. Am. A 2, 1363 (1985).
  * Stone  _et al._ [1990] J. Stone, L. Stulz, and A. Saleh, Electron. Lett. 26, 1073 (1990).
  * Stephen  _et al._ [2017] M. Stephen, M. Fahey, and I. Miller, Applied optics 56, 2636 (2017).
  * Joe  _et al._ [2000] Y. S. Joe, D. S. Ikeler, R. M. Cosby, A. M. Satanin, and C. Sub Kim, J. App. Phys. 88, 2704 (2000).
  * Lee and Kim [2001] H.-W. Lee and C. Kim, Phys. Rev. B 63, 075306 (2001).
  * Sol  _et al._ [2021] J. Sol, D. R. Smith, and P. del Hougne, arXiv:2108.06178 (2021).
  * Kuhl  _et al._ [2008] U. Kuhl, R. Höhmann, J. Main, and H. J. Stöckmann, Phys. Rev. Lett. 100, 254101 (2008).
  * Davy and Genack [2018] M. Davy and A. Z. Genack, Nat. Comm. 9, 4714 (2018).
  * [44] M. Cotrufo, E. Verhagen, and A. Fiore, in  _Quantum Sensing and Nano Electronics and Photonics XIV_ , Vol. 10111 (International Society for Optics and Photonics) p. 1011128.
  * Cogn’ee  _et al._ [2019] K. G. Cognée, W. Yan, F. La China, D. Balestri, F. Intonti, M. Gurioli, A. F. Koenderink, and P. Lalanne, Optica 6, 269 (2019).
  * Davy  _et al._ [2015] M. Davy, Z. Shi, J. Park, C. Tian, and A. Z. Genack, Nat. Commun. 6 (2015).
  * Davy  _et al._ [2021] M. Davy, M. Kuhmayer, S. Gigan, and S. Rotter, Commun. Phys 4, 85 (2021).
  * Ngo  _et al._ [2004] N. Ngo, S. F. Yu, S. Tjin, and C. Kam, Opt. Commun. 230, 115 (2004).
  * Hsue  _et al._ [2004] C.-W. Hsue, L.-C. Tsai, and K.-L. Chen, IEEE Trans. Microw. Theory Tech. 52, 1443 (2004).
  * Kulishov and Azaña [2005] M. Kulishov and J. Azaña, Opt. Lett. 30, 2700 (2005).
  * Park  _et al._ [2007] Y. Park, J. Azaña, and R. Slavík, Opt. Lett. 32, 710 (2007).
  * Li  _et al._ [2009] M. Li, D. Janner, J. Yao, and V. Pruneri, Opt. Express 17, 19798 (2009).
  * Kang and Genack [2021] Y. Kang and A. Z. Genack, Phys. Rev. B 103, L100201 (2021).

Supplementary Material for ’Exceptional points of PT-symmetric reflectionless
states in complex scattering systems”

##  I Experimental Setups

###  I.1 Multi-Mirror Fabry-Perot Interferometer

Figure S1: Colorscale representation of the reflection on a dB-scale,
$10\mathrm{log}[R(\nu,z)]$, measured with a Vector Network Analyzer between 8
and 10 GHz with zero (a) and two (b) alumina slabs at each interface.

In Fig. 1(b,c) of the main text, we have presented a colormap of the
reflection $R(\nu,z)=|r(\nu,z)|^{2}$ with respect to the frequency $\nu$ and
the depth $z$ of the two rods inserted within the waveguide when a single
dielectric alumina slab was placed between each coax-to-waveguide transition
and the waveguide. The aim of this slab is to enhance the reflectivity at the
interface and the Q-factor of the resonances. The collapse of RL states into
RL-EPs however does not depend on this reflectivity and is a more general
phenomenon. In Fig. S1, we show the same colormap for 0 and 2 dielectric slabs
at each interface.

Remarkably, the system still acts as a Fabry-Perot interferometer even without
the slabs with regularly spaced RL-frequencies in absence of the perturbation.
The reflectivity is very small but non-vanishing as the leads (coax-to-
waveguide transitions) are well but not perflectly matched to the waveguide.
The reflection average over the frequency range and the Q-factor of the
resonances are small leading to broad zeros in reflection. In contrast, for 2
dielectric slabs, the resonances with large Q-factors are peaked and the EPs
are narrower. In both cases, As $z$ increases, RL-EPs can be oberseved as $z$
increases.

Figure S2: Spectrum of the numerical (blue) and experimental (red) reflection
for an empty waveguide with one alumina slab at each interface.

We estimate the reflectively $R_{0}$ at the interfaces provided by 0, 1 or 2
dielectric slabs by fitting the spectrum of the reflection $R(\nu,z=0)$ for an
empty waveguide with its theoretical formula for a Fabry-Perot cavity

$R=1-\frac{(1-R_{0})^{2}}{1+R_{0}^{2}+2R_{0}\cos{\Delta\phi}},$ (S1)

with $\Delta\phi=2kL$ quantifying the phase shift of the wave propagating from
one Fabry-Perot mirror to the other. Here $k=\frac{2\pi}{\lambda}$ is the wave
number and $L$ the distance between the mirrors. The spectrum of the
reflection and its best fit for a single dielectric slab placed at each
interface is shown in Fig. S2. We estimate that $R_{0}=0.002$, $R_{0}=0.25$
and $R_{0}=0.46$ for 0, 1 and 2 slabs, respectively.

###  I.2 Disordered one-dimensional medium

As mentioned in the main text, RL-EPs can exist in any complex multiple-
scattering systems. We add inside the waveguide a random single-channel
disorder made of the same 28 metallic spheres as those presented in the main
text. The reflection $R(\nu,z)$ is shown in Fig. S3. Without perturbation
($z=0$), the spectrum presents randomly spaced dips. In contrast to the Fig.
3(a), one can see less RL-EPs with the random disorder than the symmetric one.
Indeed, with a symmetrical disorder, the RL-eigenvalues can leave the real
axis only by coalescing with another RL-eigenvalue whereas with a random
disorder, a RL-eigenvalue can randomly leave the real axis by itself.

Figure S3: Colorscale representation of $R(\nu,z)$ for a random disordered
medium without mirror symmetry.

##  II Theoretical analysis

In this section, we provide more details on the calculations using the
effective two-level Hamiltonian approach. The effective Hamiltonian writes

$H_{\text{eff}}=H_{0}-\frac{i}{2}V_{0}V_{0}^{T}-\frac{i}{2}V_{1}V_{1}^{T},$
(S2)

where the left and right vectors $V_{0}$ and $V_{1}$ account for the coupling
of the channels to the closed system. In the basis of two modes of the closed
waveguide, $H_{0}$ is a diagonal matrix with elements
$\omega_{0}-\delta\omega_{0}/2$ and $\omega_{0}+\delta\omega_{0}/2$. We have
identified that the coupling vectors are symmetric,
$V_{0}=\sqrt{\gamma/2}(1~{}1)$ and anti-symmetric
$V_{1}=\sqrt{\gamma/2}(1~{}-1)$, respectively (see main text). $\gamma$ is the
coupling rate of the channels to the cavity and therefore on the reflectivity
at the interfaces of the waveguide. This yields

$H_{\text{eff}}=\omega_{0}\mathbb{1}+\frac{1}{2}\begin{pmatrix}-\delta\omega_{0}-i\gamma&0\\\
0&\delta\omega_{0}-i\gamma\end{pmatrix}$ (S3)

Its eigenvalues are then

$\tilde{\omega}_{M\pm}=\omega_{0}\pm\delta\omega_{0}/2-i\gamma/2$ (S4)

The scattering matrix $S(\omega)$ can then be derived from the effective
Hamiltonian using $S(\omega)=1-iV^{T}[\omega\mathbb{1}-H_{\text{eff}}]^{-1}V$,
where $V=[V_{0}V_{1}]$. The reflection coefficient $r(\omega)$ at the left
interface is especially given by
$r=1-iV_{0}^{T}[\omega\mathbb{1}-H_{\text{eff}}]^{-1}V_{0}$. Using that

$[\omega\mathbb{1}-H_{\text{eff}}]^{-1}=\begin{pmatrix}\frac{1}{\omega-\omega_{0}+\delta\omega_{0}/2+i\gamma/2}&0\\\
0&\frac{1}{\omega-\omega_{0}-\delta\omega_{0}/2+i\gamma/2}\end{pmatrix}$ (S5)

we find

$r(\omega)=1-\frac{i\gamma/2}{\omega-\omega_{0}+\frac{1}{2}[\delta\omega_{0}+i\gamma]}-\frac{i\gamma/2}{\omega-\omega_{0}+\frac{1}{2}[-\delta\omega_{0}+i\gamma]}.$
(S6)

This equation can then be simplified to

$r(\omega)=\frac{[\omega-\omega_{0}]^{2}+[\gamma^{2}-\delta\omega_{0}^{2}]/4}{[\omega-\omega_{0}]^{2}-[\delta\omega_{0}^{2}+\gamma^{2}]/4+i\gamma[(\omega-\omega_{0}]}.$
(S7)

which is Eq. (3) of the main text.

Simple calculations show that the expression of $r(\omega)$ can also be
factorized in terms of eigenfrequencies $\tilde{\omega}_{M\pm}$ and RL-
eigenvalues
$\tilde{\omega}_{R\pm}=\omega_{0}(z)\pm\frac{1}{2}\sqrt{\delta\omega_{0}(z)^{2}-\gamma^{2}}$
as

$r(\omega)=\frac{(\omega-\tilde{\omega}_{R+})(\omega-\tilde{\omega}_{R-})}{(\omega-\tilde{\omega}_{M+})(\omega-\tilde{\omega}_{M-})}.$
(S8)

As expected, $r(\omega)$ vanishes when $\omega=\tilde{\omega}_{R\pm}$.

##  III Modal Analysis

###  III.1 Harmonic Inversion

To extract the resonances $\tilde{\omega}_{n}=\omega_{n}-i\Gamma_{n}/2$ of the
system, we perform a modal analysis between 7 and 11 GHz of the transmission
coefficient. We seek to decompose $t(\omega)$ as

$t(\omega)=\Sigma_{n}\frac{t_{n}}{\omega-\omega_{n}+i\Gamma_{n}/2}.$ (S9)

The coefficients $t_{n}$ are the modal transmission coefficients that are
associated with each resonance. We use the Harmonic Inversion (HI) method
applied to the inverse Fourier transform of $t(\omega)$ in the time domain [1,
2]. An excellent agreement with experimental results is found with 19 modes.

Figure S4: The transmission spectrum (blue line) and its reconstruction using
the modal analysis (red dashed-line). The corresponding modes are shown in
blue and orange. By compensating the absorption strength,
$\Gamma_{n}\rightarrow\Gamma_{n}-\Gamma_{a}$ with $\Gamma_{a}=30$ MHz, the
maximum of the reconstructed transmission (black line) reaches unity. The
dashed blue and orange lines are the contributions of the two modes when the
absorption is removed.

Uniform losses within the sample broaden the linewidth,
$\Gamma_{n}\rightarrow\Gamma_{n}+\Gamma_{a}$, and reduces transmission through
the sample. We now estimate the absorption rate $\Gamma_{a}$ at a RL-EP in the
case of a perturbed Fabry-Perot interferometer with one dielectic slab placed
at each interface. The maximum transmission found experimentally at the RL-EP
is 0.865. We compensate the linewidth in the reconstructed spectrum of
$T(\omega))$ and find that the transmission reaches unity for $\Gamma_{a}=30$
MHz, as shown in Fig.S4.

We observe that a slight asymmetry on the flattened line shape of the
transmission spectrum found experimentally at the EP. The strengths of the two
modes contributing to the transmission are indeed not equal. Their maxima are
$T_{n}=|t_{n}|^{2}=1.536$ and $T_{n+1}=1.398$. We attribute this difference to
the joint effects of a non-perfect symmetry of the system due to inevitable
fabrication errors and a non-uniform absorption strength. Note that the modal
strength exceeds unity as a consequence of the bi-orthogonality of quasi-
normal modes in open systems [3]. In symmetric systems, $T_{n}$ is equal to
the Petermann factor of the modes, $T_{n}=K_{n}$ which is a measure of the
degree of complexness of the eigenfunctions [4]. For isolated resonances,
$K_{n}~{}sim1$ but $K_{n}$ increases with the modal overlap between
resonances. Here, two resonances are overlapping yielding $T_{n}=K_{n}>1$.

###  III.2 Spacing between resonant frequencies

Figure S5: Spacing $\delta\omega_{0}(z)$ between central frequencies of the
two resonances close to 8.2 GHz for a multimirror cavity with one dielectric
slab on each interface. Near the EP at $z_{\mathrm{EP}}=6.2$ mm, the
experimental data (blue line) are in good agreement with a linear fit (red
line) as a function of $z$.

In Fig. S5, the spacing between two resonant frequencies $\delta\omega_{0}(z)$
is shown as a function of the penetration depth $z$ for a Fabry-Perot cavity
with a a single alumina slab placed at each interface. Near the EP at
$z_{\mathrm{EP}}=6.8$ mm, $\delta\omega_{0}(z)$ is seen to be well
approximated with a linear function. Using that
$\delta\omega_{0}(z_{\mathrm{EP}})=\gamma$, we can write
$\delta\omega_{0}(z)=\gamma+\kappa/2(z_{\mathrm{EP}}-z)$, where $\kappa$ is a
positive coefficient. The spacing between real RL-eigenvalues
$\delta\tilde{\omega}_{R}=\sqrt{\delta\omega_{0}(z)^{2}-\gamma^{2}}$ gives
$\delta\tilde{\omega}_{R}=\sqrt{\kappa\gamma(z_{\mathrm{EP}}-z)+\kappa^{2}/4(z_{\mathrm{EP}}-z)^{2}}$.
In the vicinity of the EP, the square root behavior
$\delta\tilde{\omega}_{R}\sim\sqrt{\kappa\gamma(z_{\mathrm{EP}}-z)}$ reflects
the square-root sensitivity of EPs to an external perturbation of the system.
The same is obtained for $z>z_{\mathrm{EP}}$ on the splitting of the imaginary
parts of the two complex-conjugate RL-eigenvalues.

##  IV Multichannel cavity

Figure S6: (a) Photography of the Multi-channel cavity. (b) Sketch of the
cavity for which the top has been removed to see the symmetric disorder inside
the cavity. Figure S7: Colorscale representation of the last reflection
eigenvalue $\tau_{4}(\nu,z)$ on a dB-scale for the multichannel case. The
minimal reflection is shown for the left (top) and right (bottom) reflection
matrices.

The reflection coefficient $r(\nu)$ and the transmission coefficient $t(\nu)$
are measured through a rectangular multi-channel cavity (length $L=500$ mm,
width $W=250$ mm and height $H=8$ mm) between two arrays of $N=4$ single-
channel. Note that $8$ transitions are observed at each side of the cavity on
Fig. S6. However, only 4 of them are used to measure the reflection matrix.
The others are not connected to the VNA. This open-circuit condition is
similar to a metallic boundary conditions for the cavity.

An array of $14\times 20$ holes are drilled symmetrically with respect to the
center on the top plate of the cavity with a spacing of $12$ mm between each
hole (see Fig.S6(a)). Ten aluminium rods are inserted symmetrically through
these holes thus creating a symmetric disorder inside the cavity (see Fig.
S6(b)). As with the single mode waveguide presented in the main text, a
perturbation is then introduced in the cavity by inserting symmetrically two
aluminium rods spaced by $12$ mm near to the center. The penetration depth of
these two rods is tuned to observe the RL-EPs.

Due to its large dimensions and inevitable fabrication errors, the cavity is
not perfectly symmetric. However, we verify that the RL-EP presented in the
main text is found on the both the left and right reflection matrices. In Fig.
S7, we present the colorscale representation of the smallest reflection
eigenvalues $\tau_{4}(\nu,z)$ for both matrices. The multichannel RL-EP can be
clearly observed on the two figures as expected.

## References

  * Mandelshtam and Taylor [1997] V. A. Mandelshtam and H. S. Taylor, J. Chem. Phys. 107, 6756 (1997).
  * Davy and Genack [2018] M. Davy and A. Z. Genack, Nat. Comm. 9, 4714 (2018).
  * Davy and Genack [2019] M. Davy and A. Z. Genack, Phys. Rev. Res. 1, 033026 (2019).
  * Poli  _et al._ [2010] C. Poli, O. Legrand, and F. Mortessagne, Phys. Rev. E 82, 055201(R) (2010).