# Zero-Shot Personalized Speech Enhancement
through Speaker-Informed Model Selection

###### Abstract

This paper presents a novel zero-shot learning approach towards personalized
speech enhancement through the use of a sparsely active ensemble model.
Optimizing speech denoising systems towards a particular test-time speaker can
improve performance and reduce run-time complexity. However, test-time model
adaptation may be challenging if collecting data from the test-time speaker is
not possible. To this end, we propose using an ensemble model wherein each
specialist module denoises noisy utterances from a distinct partition of
training set speakers. The gating module inexpensively estimates test-time
speaker characteristics in the form of an embedding vector and selects the
most appropriate specialist module for denoising the test signal. Grouping the
training set speakers into non-overlapping semantically similar groups is non-
trivial and ill-defined. To do this, we first train a Siamese network using
noisy speech pairs to maximize or minimize the similarity of its output
vectors depending on whether the utterances derive from the same speaker or
not. Next, we perform k-means clustering on the latent space formed by the
averaged embedding vectors per training set speaker. In this way, we designate
speaker groups and train specialist modules optimized around partitions of the
complete training set. Our experiments show that ensemble models made up of
low-capacity specialists can outperform high-capacity generalist models with
greater efficiency and improved adaptation towards unseen test-time speakers.

Index Terms—  Speech enhancement, deep learning, adaptive mixture of local
experts, model compression by selection

##  1 Introduction

Speech enhancement (SE) is a long-standing research area within signal
processing [1, 2, 3] which has experienced significant progress in the past
decade due to the pervasiveness of machine learning models and deep neural
networks (DNNs) [4, 5, 6, 7]. In this paper, we address the task of removing
non-stationary background noise from a monophonic recording of a single
speaker. Most studies in this area put forward models which perform
generalizable speech enhancement—that is, the denoising models do not make
assumptions about the test-time speaker or noisy environment. To address many
kinds of speaker types and noise types, most DNN-based general-purpose models
operate using millions of learnable parameters. However, given the ubiquity of
resource-constrained devices (i.e., smartphones or smart speakers), we focus
our interest on developing SE models which minimize run-time computational
complexity without sacrificing denoising quality.

As it is known that models trained specifically for a disjoint sub-problem
(i.e., specialists) outperform the ones that aim at the universal speech
enhancement task (i.e., generalists) [8], classifying which sub-problem the
test signal belongs to can help improve the performance of the SE problem. A
popular neural network design philosophy that handles those decomposed sub-
problem individually is the mixture of local experts (MLoE) [9]. An ensemble
model consists of independent expert modules, each of which addresses a subset
of all the training cases as a specialist. Hence, an auxiliary classifier
module must be trained such that it estimates the ensemble weights for all the
local experts based on their relevance towards a particular input example. The
weights are then used to emphasize more relevant specialists’ results.

Recent studies applied the MLoE paradigm directly towards speech enhancement.
In [10] local experts are associated with different hyperparameter choices,
such as different contextual window lengths, and then they are stacked up to
form multicontext stacking. MLoE is also extended to recurrent neural network
models to learn from the temporal structure of the speech enhancement problem
in [11]. While in these methods a local expert is not dedicated to solve a
specific sub-problem, in [12] MLoE-based SE system shows substantial
improvement by pre-defining two separate partitioning schemes: based on the
quality of input signal, i.e., in terms of signal-to-noise ratio (SNR), and
the gender of the speakers. Moreover, by introducing “sparseness” to the
ensemble weights, it performs test-time inference on only one most suitable
specialist. Compared to an ordinary generalist that needs to be very large to
achieve a certain level of speech denoising, each specialist can achieve the
same quality enhancement results using much smaller architecture. Hence, the
sparse ensemble of specialists is claimed as a model compression method in
[12]. Zezario et. al also proposed a similar partitioning strategy based on
the speech quality [13], but by employing a speech quality estimator in place
of the traditional gating module. In [13], it is also proposed to learn the
SNR-based partitions rather than pre-defining them.

In this paper, we propose to use MLoE as a means for personalizing the SE
model. In doing so, instead of partitioning the dataset into pre-defined
subsets, we take a more adaptive approach and learn the optimal speaker
grouping strategy from the data. With the speaker groups in place, the gating
module must estimate characteristics of the test-time speaker (i.e., by
computing an embedding), identify the most similar speaker group defined
within the training set, then forward the input signal to the appropriate
specialist network. This schema requires no training data from the test-time
speakers but more optimally denoises the test-time noisy utterances by
estimating the most suitable specialist. The idea of “zero-shot” speech
enhancement through model selection has seen preliminary assessment, however,
these works were limited to doing model selection based on the speaker-
agnostic aspects of the signals, such as quality of the test-time signal, i.e.
its signal-to-noise (SNR) level, [14, 13, 12], types of the noise sources
[15], or the gender [12].

Another important aspect of our work revolves around an open-ended question:
how do we cluster English speakers into appropriate groups? This question is
related to the task of learning speaker-characteristic embeddings, which has
been extensively researched in the last decade, largely for improving speaker
verification (SV) systems. Well-established embeddings include i-vectors,
which are calculated using a Gaussian mixture model [16], or x-vectors
computed from a time-delay neural network [17]. Autoencoders’ bottleneck
layers are also a popular choice for learning fixed-length embeddings, which
have been used to compress Mel-frequency cepstral coefficients (MFCCs) into
discriminative embeddings either through contrastive loss on pairs of inputs
[18] or by estimating subsequent frames for a single input signal [19]. While
these formulations are valid for learning speaker-identifying features, we
propose to customize this embedding learning models further, so they
effectively function as the gating module that is specialized for our SE
problem rather than verification. To this end, we will first employ a Siamese
network [20, 21] to learn the discriminative speaker embeddings, and then turn
it into a classifier as our gating module. Fine-tuning follows to further
evolve the learned embedding space into something more suitable for the SE
setup.

##  2 Proposed Methods

###  2.1 Ensemble Models

Given a large dataset containing many different speakers’ various utterances
$\mathbb{S}$, we postulate that there exists an optimal clustering within the
data based around the speaker’s identifying characteristics. Denoting $K$ to
be the number of clusters, one can create $K$ separate SE models trained only
to denoise utterances from each disjoint group of similar speakers. As
previously shown [12, 13], a sparsely active ensemble model is capable of
performing zero-shot adaptation because the gating module classifies the test-
time noisy input signals into one-of-$K$ groups.

Figure 1: The proposed sparse ensemble of specialists model.

We illustrate the architecture of a sparsely active ensemble model in Figure
1. An ensemble model is composed of one gating module and $K$ specialist
modules. The gating module processes a noisy speech input frame $\bm{x}$,
estimating a speaker-embedding first and then classifying it as belonging to
one-of-$K$ groups. The cluster probabilities vector $\bm{p}$ is used in two
ways—during training, all of the specialist modules outputs their own ideal
ratio mask (IRM) [22] estimates,
$\bm{m}^{(1)},\bm{m}^{(2)},\ldots,\bm{m}^{(K)}$, which are then combined in a
weighted sum using $\bm{p}$, i.e.,
$\hat{\bm{m}}=\sum_{k=1}^{K}p_{k}\bm{m}^{(k)}$. But during testing, only the
output from the $k^{*}$-th specialist, corresponding to the largest
probability, i.e., $k^{*}=\operatornamewithlimits{\arg\max}_{k}p_{k}$, is
chosen. This $\operatorname{argmax}$ operation selects a single specialist to
use during evaluation, making the ensemble sparsely active.

In the context of personalized speech enhancement, increasing hyper-parameter
$K$ can theoretically increase the level of specialization of each specialist
as well as the ensemble network’s capacity for personalization. However, there
is a trade-off with having too many models; a large $K$ can make the gating
module’s classification task too challenging, and may lead to specialist
modules becoming overfit on too small subset. In this paper, we investigate
three choices of $K$: $\mathrm{2}$, $\mathrm{5}$, and $\mathrm{1}\mathrm{0}$.
Determining the optimal number of clusters is an extended research topic
within unsupervised learning.

###  2.2 Discriminative Speaker-Specific Embeddings

The clustering of speakers is a significant matter when we build a successful
sparse ensemble model for SE. Although in theory all the specialists and the
gating module can be trained from scratch, training many modules
simultaneously is prone to result in suboptimal performance. Hence, we first
pre-train all the modules individually and then fine-tune them. The pre-
training step, therefore, requires sub-grouping of speakers.

To this end, we train a neural encoder that learns an an embedding function
$f$ which can characterize a noisy speech utterance with a low-rank embedding
vector. In order to train $f$, we formulate a speaker verification (SV)
upstream task. First, we sample utterances from a large training dataset
containing many speakers, $\bm{s}\in\mathbb{S}$, and noise signals from a
similarly large dataset of diverse noises, $\bm{n}\in\mathbb{N}$. Input
mixtures $\bm{x}$ are made by artificially mixing clean speech utterances
$\bm{s}$ with training noise signals $\bm{n}$; the amplitude of $\bm{n}$ is
scaled to simulate various signal-to-noise ratios (SNRs).

We can then generate pairs of noisy speech utterances, $\bm{x}_{i}$ and
$\bm{x}_{j}$. Once $f$ predicts the embeddings, i.e.,
$\bm{z}_{i}=f(\bm{x}_{i})$ and $\bm{z}_{j}=f(\bm{x}_{j})$, their inner product
serves as a measure of similarity. A sigmoid function follows to interpret it
as a probability $\hat{{y}}$. Our target is a binary value ${y}$, either
$\mathrm{1}$ or $\mathrm{0}$ depending on whether the utterances derive from
the same speaker or not. The embedding function $f$ is trained to minimize the
binary cross entropy loss between $\hat{y}$ and $y$.

This contrastive learning approach for deriving discriminative embeddings
follows the steps popularized by Siamese networks [20, 21] in that the same
embedding function $f$ is used for both input signals $\bm{x}_{i}$ and
$\bm{x}_{j}$. The rationale behind this embedding model is that the
discriminative nature of these embeddings can help the clustering process
prepare a semantically more meaningful partitioning of speakers.

(a) Speaker means obtained from SV with $K=2$ clustering.

(b) Speaker means obtained from SV with $K=10$ clustering.

(c) Speaker means derived uniquely by a fine-tuned $K=2$ ensemble.

(d) Speaker means derived uniquely by a fine-tuned $K=10$ ensemble.

Figure 2: Subplots comparing various choices of $K$ for using k-means
clustering on the speaker embeddings. The speaker verification (SV) pre-
training task creates a latent space of speaker embeddings $\mathbb{Z}$, from
which we can partition various groups, i.e. 2 in (a) and 10 in (b). After
fine-tuning an ensemble model, the gating network’s embedding function $f$
adjusts its parameters towards the speech enhancement (SE) objective. The
latent space is modified uniquely based on the ensemble’s configuration. In
(a) and (b), the class labels derive from k-means clustering, but in (c) and
(d) the class labels are estimated by the gating network’s classifier function
$g$.

###  2.3 Offline Speaker Clustering

Likewise, the gating module’s classification task and pre-training of
individual specialists rely on a reasonable clustering of speakers.
Determining how the $K$ groups are formulated, and which of the training set
speakers belong to each group, requires an offline clustering step. First, we
transform every utterance from the training corpus into the learned latent
space, i.e., $\bm{z}\leftarrow f(s)$. Embedding vectors from the same speaker
are averaged element-wise, which serves as the speaker-characteristic mean
vector. Finally, we apply k-means clustering to these mean vectors, which
learns $K$ speaker groups.

Figure 2 illustrates the clustering results with different choices of $K$.
Each of the $\mathrm{2}\mathrm{1}\mathrm{1}$ points represents one of the
Librispeech training set speakers, with marker style indicating speaker
gender. For plotting, the $\mathrm{3}\mathrm{2}$-dimensional embeddings
$\bm{z}$ are reduced to $\mathrm{2}$ dimensions using t-SNE (with
$\text{perplexity}=40$) [23]. These subplots show that the SV model succeeds
in learning a speaker embedding which can be clustered into loosely meaningful
groups, e.g., when $K=2$ the clusters implicitly form along the speaker gender
division. These speaker groups are used to pre-train our gating modules and
speaker-specific local experts.

###  2.4 Gating Module Pre-Training

The gating module must be able to classify the embedding vectors as belonging
to one of the $K$ speaker clusters. This neural network is a dense layer
followed by the softmax activation, which we denote by a parametric function
$\bm{p}=g(\bm{z};\mathcal{W}_{g})$, where $\mathcal{W}_{g}$ is its parameters.
The classifier function $g$ takes embeddings of noisy utterances $\bm{z}$ as
inputs, and outputs a vector of cluster probabilities $\hat{\bm{p}}$. As each
utterance belongs to a single cluster and the speaker IDs of the training set
speakers are known, we can encode the k-means clustering labels into one-hot
vector targets $\bm{p}$. These vectors are $K$-dimensional.

Note the discrepancy between the clustering done on embeddings of the clean
speech utterances and the actual use-case of the model that takes noisy
utterances. While the clustering results on clean data might be more reliable,
eventually it is always possible that a noisy test utterance can be
misclassified into a wrong speaker group, and then consequently assigned to a
suboptimal specialist. Moreover, since the embeddings are optimized for the SV
tasks, clustering on this representation may not be optimal for our SE
problem. We revisit this issue in Sec. 2.6 and propose a fine-tuning solution.

###  2.5 Specialist Pre-Training

The $K$ specialist modules are trained towards speech denoising in the
following manner. The large dataset of training noises $\mathbb{N}$ is
retained, but the large speech corpus $\mathbb{S}$ is partitioned into $K$
groups, based on the clustering results in Sec. 2.3. We denote the $k$-th
partition as $\mathbb{S}^{(k)}$. The $k$-th specialist module can be described
as a mapping function $h$ that updates its parameters $\mathcal{W}_{h}$ with
each iteration such that the distance $\mathcal{E}$ between the denoised
estimate signal $\hat{\bm{s}}$ and the target clean speech signal $\bm{s}$ is
minimized. We use the negative scale-invariant signal-to-distortion ratio (SI-
SDR) [24] as the loss function.

###  2.6 Ensemble Fine-Tuning

The ensemble model can now be used naïvely by assembling the pre-trained
specialist modules and a pre-trained gating module. However, we note that the
gating module may not classify all input signals with perfect accuracy.
Therefore, the ensemble model can benefit from fine-tuning (FT) in order to
adjust denoising performance for misclassified inputs. This potential co-
adaptation between gating and specialist modules can be found by adjusting the
parameters of all the underlying functions (i.e., embedding function $f$,
classifier function $g$, and denoising functions $h$ within each specialist).
In the fine-tuning phase, the ensemble model estimates the final ratio mask
$\hat{\bm{m}}$ by performing a normalized sum over the individual masks
$\bm{m}^{(k)}$ using the softmax vector, $\hat{\bm{p}}$, i.e.,
$\hat{\bm{m}}=\sum_{k=1}^{K}\hat{p}_{k}\bm{m}^{(k)}$. This ensures that the
ratio mask calculation is differentiable and can be seen as a “soft” gating
mechanism.

During testing, the weighted sum is replaced by a hard-decision, i.e.
$\hat{\bm{m}}=\bm{m}^{(k^{\ast})}$ where
$k^{\ast}=\operatorname*{argmax}_{k}p_{k}$. This switch in gating mechanism
between training- and evaluation-time is the essence of the ensemble scheme’s
efficiency: only one out of all the specialists is active during inference,
making the total used network parameters a fraction of the total learned. In
order to reduce the discrepancy between the hard and soft gating mechanisms, a
scaling parameter $\lambda=10$ is applied to the logit values to make the
softmax function harder.

Figure 2c and 2d show the speaker embedding vectors after fine-tuning. Note
that the comparison between the clustering on the SV embedding vectors and on
their fine-tuned version is not to argue that fine-tuning can improve the
clustering results. Instead, it is possible that fine-tuning with the speech
enhancement objective could in fact deteriorate the discriminative qualities
of the learned embedding vectors.

##  3 Experiment Setup

Mixtures are generated by combining randomly offset $5\text{\,}\sec$ segments
of utterances and noises. With every mixture, the noise signal is randomly
scaled such that the mixture SNR lies uniformly between
$-5\text{\,}\mathrm{t}$o $10\text{\,}\mathrm{dB}$. Utterances derive from the
LibriSpeech corpus [25] train-clean-100 folder, with
$\mathrm{2}\mathrm{1}\mathrm{1}$ speakers designated in the training set,
$\mathrm{2}\mathrm{0}$ in the validation set, and $\mathrm{2}\mathrm{0}$ in
the test set. Noises are selected from the MUSAN corpus [26], with
$\mathrm{6}\mathrm{2}\mathrm{8}$ noises from the free-sound folder used during
training and validation, and $\mathrm{5}\mathrm{4}$ noises from the sound-
bible folder used during test. Both LibriSpeech and MUSAN corpora are
resampled to $8\text{\,}\mathrm{kHz}$. When training the speaker verification
model, batches are made up of pairs of mixtures, with an equal chance of being
from the same speaker or not.

All mixture signals are processed in the time-frequency domain through STFT
using a frame size of $\mathrm{1}\mathrm{0}\mathrm{2}\mathrm{4}$ samples with
$75\text{\,}\%$ overlap. Throughout our experiment, every model performs
speech denoising by taking a series of magnitude spectra as input and
estimating IRM vectors $\bm{m}$. Masking is done element-wise onto the
complex-valued spectrum which possesses the noisy phase of the mixture signal.

Both the gating and specialist modules are composed of gated recurrent units
(GRU) cells [27]. The embedding function $f$ is built with $\mathrm{2}$ hidden
layers and $\mathrm{3}\mathrm{2}$ hidden units, with the output from last
frame becoming a fixed-length utterance-characteristic embedding $\bm{z}$. The
denoising functions $h$ are also built with 2 hidden layers but with a varied
number of hidden units. The baseline general-purpose SE model is constructed
in exactly the same manner as a specialist network, but is trained on the
entire speech corpus $\mathbb{S}$ instead of a personalized subset
$\mathbb{S}^{(k)}$. Throughout the experiment, we opt for a batch size of
$\mathrm{1}\mathrm{2}\mathrm{8}$, training all models using the Adam optimizer
[28] with learning rates of $10^{-3}$ for training and $10^{-4}$ for fine-
tuning.

##  4 Results

Figure 3: Comparison of speech enhancement performance between a baseline
general-purpose model against different configurations of speaker-informed
sparse ensemble models.

Figure 3 summarizes the findings of our experiments. The x-axis shows the
varying hidden sizes for the GRU layers. In order to calculate the total
number of parameters, the number in parenthesis must be multiplied by $K$,
i.e., for the ensemble model with $\mathrm{5}$ specialists of hidden size
$\mathrm{2}\mathrm{5}\mathrm{6}$, the total number of parameters equals
$5.6\text{\,}\mathrm{M}$. Because our ensemble models are sparsely active—that
is, one specialist is active at a time—the number of parameters effective at
run-time is only $1.1\text{\,}\mathrm{M}$, the amount listed on the x-axis.
Longitudinally, the baseline models share the same number of hidden units with
the specialist module, meaning the total number of model parameters for the
baseline is always $K$ times smaller than the ensemble model in comparison.
However their effective number of parameters is nearly equivalent. We note
that ensemble models are not fine-tuned for hidden sizes
$\geq$\mathrm{5}\mathrm{1}\mathrm{2}$$ due to GPU memory constraints. Larger
baseline models are trained and evaluated for comparison with the smaller
ensemble models.

Firstly, we see that across all configurations, our ensemble models
consistently yields a higher denoising performance when compared to a baseline
generalist model whose size is similar to one of the specialists. The naïve
ensemble models already show significant improvement (ranging from
$0.62\text{\,}\mathrm{t}$o $1.65\text{\,}\mathrm{dB}$), but different choices
of $K$ do not make a big difference. We also observe that fine-tuning the
ensemble models lift the performance even further (from
$1.24\text{\,}\mathrm{t}$o as much as $2.04\text{\,}\mathrm{dB}$. Furthermore,
fine-tuning introduces a larger gap in improvement when $K$ is larger;
intuitively, the more challenging classification task stands to benefit most
from fine-tuning.

The proposed method also performs model compression without sacrificing the
denoising performance. Overall, the smaller model architecture receives more
performance improvement, such as the $2.0\text{\,}\mathrm{dB}$ improvement in
the case of $\mathrm{6}\mathrm{4}$ hidden units. The model compression
benefits are made clear by comparing data points laterally. For example, as
circled in Figure 3, a generalist model requires at least
$\mathrm{5}\mathrm{1}\mathrm{2}$ hidden units in order to match the
performance of a fine-tuned ensemble model with $\mathrm{1}\mathrm{0}$
specialists each made up of GRUs with only $\mathrm{6}\mathrm{4}$ hidden
units. Including the cost of the gating module and all the other specialists
that are not chosen, this is still a $\mathrm{4}\mathrm{8}$% reduction in
terms of spatial complexity. Moreover, if we only count the gating module and
one chosen specialist, it is a $\mathrm{9}\mathrm{4}$% reduction in effective
parameters, greatly reducing arithmetic complexity during test-time.

Lastly, as hypothesized, we see that increasing the number of clusters results
can result in a more personalized speech enhancement so long as the ensemble
model is fine-tuned. The average SI-SDR improvement achieved with the ensemble
models increases along with $K$ from $\mathrm{2}$ to $\mathrm{5}$ to
$\mathrm{1}\mathrm{0}$ through fine-tuning.

##  5 Conclusion

In this paper, we investigated model adaptation through selection (the
“mixture of local experts” paradigm) as a means for zero-shot personalized
speech enhancement. Our method is zero-shot as the ensemble models are never
exposed to the test-time speaker or test-time environment; instead, the gating
module analyzes the test signal to determine the most similar training cases
and selects the most appropriate specialist for denoising. We obtain a
speaker-informed gating module by pretraining it with a constrastive speaker
verification task. The training cases are transformed to a learned latent
space where they are clustered using k-means. By identifying more clusters and
training more low-cost specialists, our ensemble models are able to adapt
better to unseen speakers. Our findings reinforce the idea that sparse
ensemble models can outperform general-purpose speech denoising models of a
similar architecture, additionally reducing run-time computational complexity.

## References

  * [1] S. F. Boll, “Suppression of acoustic noise in speech using spectral subtraction,” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. 27, pp. 113–120, 1979.
  * [2] Y. Ephraim and D. Malah, “Speech enhancement using a minimum-mean square error short-time spectral amplitude estimator,” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. 32, no. 6, pp. 1109–1121, 1984\.
  * [3] S. Gannot, D. Burshtein, and E. Weinstein, “Iterative and sequential Kalman filter-based speech enhancement algorithms,” _IEEE Transactions on Speech and Audio Processing_ , vol. 6, pp. 373–385, 1998.
  * [4] Y. Xu, J. Du, L.-R. Dai, and C.-H. Lee, “An experimental study on speech enhancement based on deep neural networks,” _IEEE Signal Processing Letters_ , vol. 21, no. 1, pp. 65–68, 2014.
  * [5] P. Huang, M. Kim, M. Hasegawa-Johnson, and P. Smaragdis, “Joint optimization of masks and deep recurrent neural networks for monaural source separation,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 23, no. 12, pp. 2136–2147, Dec 2015.
  * [6] F. Weninger, H. Erdogan, S. Watanabe, E. Vincent, J. Le Roux, J. R. Hershey, and B. Schuller, “Speech Enhancement with LSTM Recurrent Neural Networks and its Application to Noise-Robust ASR,” in _Proceedings of the International Conference on Latent Variable Analysis and Signal Separation (LVA/ICA)_ , Aug. 2015.
  * [7] D. L. Wang and J. Chen, “Supervised speech separation based on deep learning: An overview,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 26, no. 10, pp. 1702–1726, 2018.
  * [8] M. Kolbæk, Z. H. Tan, and J. Jensen, “Speech intelligibility potential of general and specialized deep neural network based speech enhancement systems,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 25, no. 1, pp. 153–167, Jan 2017.
  * [9] R. A. Jacobs, M. I. Jordan, S. J. Nowlan, and G. E. Hinton, “Adaptive mixtures of local experts,” _Neural Computation_ , vol. 3, no. 1, pp. 79–87, Mar. 1991.
  * [10] X.-L. Zhang and D. Wang, “A deep ensemble learning method for monaural speech separation,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 24, no. 5, pp. 967–977, 2016.
  * [11] S. Chazan, J. Goldberger, and S. Gannot, “Deep recurrent mixture of experts for speech enhancement,” in _Proceedings of the IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA)_. IEEE, 2017, pp. 359–363.
  * [12] A. Sivaraman and M. Kim, “Sparse Mixture of Local Experts for Efficient Speech Enhancement,” in _Proceedings of the Annual Conference of the International Speech Communication Association (Interspeech)_ , 2020, pp. 4526–4530.
  * [13] R. E. Zezario, C.-S. Fuh, H.-M. Wang, and Y. Tsao, “Speech enhancement with zero-shot model selection,” _arXiv preprint arXiv:2012.09359_ , 2020.
  * [14] R. E. Zezario, S.-W. Fu, X. Lu, H.-M. Wang, and Y. Tsao, “Specialized speech enhancement model selection based on learned non-intrusive quality assessment metric,” in _Proceedings of the Annual Conference of the International Speech Communication Association (Interspeech)_ , 2019.
  * [15] M. Kim, “Collaborative deep learning for speech enhancement: A run-time model selection method using autoencoders,” in _Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP)_ , 2017.
  * [16] N. Dehak, P. J. Kenny, R. Dehak, P. Dumouchel, and P. Ouellet, “Front-End Factor Analysis for Speaker Verification,” _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 19, no. 4, pp. 788–798, 2010.
  * [17] D. Snyder, D. Garcia-Romero, D. Povey, and S. Khudanpur, “Deep Neural Network Embeddings for Text-Independent Speaker Verification,” in _Proceedings of the Annual Conference of the International Speech Communication Association (Interspeech)_ , 2017, pp. 999–1003.
  * [18] K. Chen and A. Salman, “Learning speaker-specific characteristics with a deep neural architecture,” _IEEE Transactions on Neural Networks_ , vol. 22, no. 11, pp. 1744–1756, 2011.
  * [19] A. Jati and P. G. Georgiou, “Speaker2Vec: Unsupervised Learning and Adaptation of a Speaker Manifold Using Deep Neural Networks with an Evaluation on Speaker Segmentation,” in _Proceedings of the Annual Conference of the International Speech Communication Association (Interspeech)_ , 2017, pp. 3567–3571.
  * [20] J. Bromley, I. Guyon, Y. LeCun, E. Säckinger, and R. Shah, “Signature verification using a “siamese” time delay neural network,” in _Advances in Neural Information Processing Systems (NIPS)_ , 1994, pp. 737–744.
  * [21] D. Chicco, _Siamese Neural Networks: An Overview_. New York, NY: Springer US, 2021, pp. 73–94.
  * [22] A. Narayanan and D. L. Wang, “Ideal ratio mask estimation using deep neural networks for robust speech recognition,” in _Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP)_ , May 2013, pp. 7092–7096.
  * [23] L. V. der Maaten and G. Hinton, “Visualizing Data using t-SNE,” _Journal of Machine Learning Research_ , vol. 9, no. 11, 2008.
  * [24] J. L. Roux, S. Wisdom, H. Erdogan, and J. R. Hershey, “SDR - half-baked or well done?” _arXiv preprint arXiv:1811.02508_ , 2018.
  * [25] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur, “Librispeech: An ASR corpus based on public domain audio books,” in _Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP)_. IEEE, 2015, pp. 5206–5210.
  * [26] D. Snyder, G. Chen, and D. Povey, “MUSAN: A Music, Speech, and Noise Corpus,” _arXiv preprint arXiv:1510.08484_ , 2015.
  * [27] K. Cho, B. Van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio, “Learning phrase representations using RNN encoder-decoder for statistical machine translation,” _arXiv preprint arXiv:1406.1078_ , 2014.
  * [28] D. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _Proceedings of the International Conference on Learning Representations (ICLR)_ , 2015.


aainstitutetext: Department of Astronomy and Theoretical Physics, Lund
University, S-223 62 Lund, Swedenbbinstitutetext: Institute for Particle
Physics Phenomenology, Department of Physics, Durham University, Durham DH1
3LE, U.K.ccinstitutetext: High Energy Physics Group, Blackett Laboratory,
Imperial College, Prince Consort Road, London, SW7 2AZ, United Kingdom

# Collider Events on a Quantum Computer

Gösta Gustafson a Stefan Prestel b Michael Spannowsky c Simon Williams

###### Abstract

High-quality simulated data is crucial for particle physics discoveries.
Therefore, parton shower algorithms are a major building block of the data
synthesis in event generator programs. However, the core algorithms used to
generate parton showers have barely changed since the 1980s. With quantum
computers’ rapid and continuous development, dedicated algorithms are required
to exploit the potential that quantum computers provide to address problems in
high-energy physics. This paper presents a novel approach to synthesising
parton showers using the Discrete QCD method. The algorithm benefits from an
elegant quantum walk implementation which can be embedded into the classical
toolchain. We use the ibm_algiers device to sample parton shower
configurations and generate data that we compare against measurements taken at
the ALEPH, DELPHI and OPAL experiments. This is the first time a Noisy
Intermediate-Scale Quantum (NISQ) device has been used to simulate realistic
high-energy particle collision events.

IPPP/22/43 LU-TP-22-49

##  1 Introduction

Particle physics studies the smallest building blocks of matter and their
interactions by analysing the outcome of high-energy particle scattering
events produced in large-scale particle colliders. In such collision events,
the colliding particles annihilate, and the resulting energy often creates
highly energetic particles. These short-lived states quickly decay into high
occupation-number (i.e. large particle multiplicity), low average-energy
states through a process called “showering”. This phenomenon is particularly
important for colour-charged partons, such that developing “parton shower”
calculations Fox:1979ag ; Sjostrand:1985xi ; Gustafson:1986db  have become a
crucial aspect of research in Quantum Chromodynamics (QCD) Lonnblad:1992tz ;
Gieseke:2003rz ; Sjostrand:2004ef ; Schumann:2007mg ; Platzer:2009jq ;
Nagy:2014mqa ; Hoche:2015sya ; Fischer:2016vfv ; Hamilton:2020rcu .

All processes in particle physics are inherently quantum mechanical, making
particle physics an excellent driver for quantum algorithm development, as
pointed out early on by Feynman Feynman:1981tf . One major obstacle for such
applications is tackling environmental noise within the device, which can
cause decoherence effects, drastically reducing the fidelity of quantum
computers. Rapid progress in manufacturing and running Noisy Intermediate-
Scale quantum (NISQ) devices has led to renewed interest in particle physics
applications Jordan:2014tma ; Garcia-Alvarez:2014uda ; Arrighi:2018PRA ;
Marque-Martin:2018PRA ; Alexandru:2019nsa ; Jay:2019PRA ; Wei:2019rqy ;
Lamm:2019bik ; MottQuantum ; Bauer:2019qxa ; Alexandru:2019ozf ;
Blance:2020ktp ; Lamm:2019uyc ; Abel:2020ebj ; Abel:2020qzm ; Bepari:2020xqi ;
Williams:2021lvr ; DiMolfetta:2020QIP ; Araz:2022haf ; Matchev:2020wwx ;
DeJong:2020riy ; Ngairangbam:2021yma ; Bauer:2022hpo ; Agliardi:2022ghn ,
allowing to test the limits of current quantum hardware. However, these
applications are, so far, limited to proof-of-principle studies, often
focusing on implementations of low-dimensional field theories or simplified
models of quantum circuits.

The QCD parton showering process is crucial in modelling and analysing data
from particle colliders. This success of parton shower programs is based on
the ability to produce “synthetic” scattering event data, containing
physically meaningful particle final states that obey momentum conservation
Buckley:2011ms . Such simulated data can be compared to actual scattering
event data recorded by particle detectors, hinting at necessary improvements
in the theory prediction – or discoveries. Due to its relevance, the showering
process has also generated interest in the quantum computing community
Bauer:2019qxa ; Bepari:2020xqi ; Williams:2021lvr ; Deliyannis:2022uyh . These
attempts considered the showering of individual, isolated particles. The
generation of physically meaningful events was deemed impossible, in
particular, since it is assumed that each particle produced by the shower
samples had continuous momentum-space quantum numbers. However, the usefulness
of a semi-classical picture of factorising the model of scattering events into
the evolution of individual partons (inspired by the early work of Feynman and
Field Field:1977fa  and by factorisation proofs for specific, very simple one-
dimensional measurements Collins:1984kg ; Collins:1989gx ) has to be discarded
for the multi-dimensional modelling of scattering events that make current
parton showers successful.

This note proposes the first quantum algorithm capable of simulating realistic
scatting events that can be compared to data recorded at collider experiments,
achieved by reframing the classical parton shower algorithm as a quantum walk
PhysRevA.48.1687 ; QWProc ; Kempe . The reason to attempt implementing a
quantum algorithm for parton showers is twofold. First, on an abstract level,
quantum algorithms might eventually help overcome classical bottlenecks. These
might be related to raw computing power or overly complicated algorithms:
conventional parton showering programs implement very complicated code
structures when embedding quantum effects due to corrections to the
“factorisation” picture above, both for high-energy configurations and low-
energy “soft” gluons. Thus, the development of a quantum shower algorithm
allows to increase also the classical “toolbox” (which effectively consists of
a single algorithm, the “veto algorithm”, see e.g. Sjostrand:2006za ;
Bierlich:2022pfr  for details) and might provide more natural realisations of
quantum coherence effects. More pragmatically, we also aim to assess the
relevance of quantum noise on a real-world example – the description of data
collected by collider experiments.

In pursuing this goal, we are led to several improvements over previous work.
The most important aspect is the realistic handling of parton decay
kinematics, defined for any momentum-space configuration. This naturally
includes momentum-conservation constraints. Other important realisations are
that the energy-scale dependence of the strong coupling $\alpha_{s}$ – which
was omitted previously – and using soft-gluon (quantum) interference as a
guiding principle for quantum parton shower developments allows clarifying the
treatment of momentum-space sampling sufficiently to permit the generation of
event data.

##  2 The Parton Shower Algorithm

Parton shower processes evolve high-energy few-particle states to low-energy
multi-particle states by successively decaying particles into lower-energy
decay products. This process is dominated by collinear decay modes, in which
the decay products are almost parallel, or soft modes, which feature at least
one vanishing-energy decay product. Within the former limit, successive decay
steps factorise into independent quasi-classical steps, whereas interference
contributions only allow for partial factorisation in the latter.
Interestingly, the leading contributions to the decay rate in the collinear
limit are included in the “soft” decay rate. This makes the soft limit an
excellent starting point for developing parton showers Azimov:1985zta ;
Gustafson:1986db . In this limit, the decay of a high-energy state proceeds by
decaying colour-anticolour dipoles into lower-energy dipoles. Interpreting the
showering process as dipole decays allows to describe the most straightforward
interference patterns but also addresses momentum conservation within
sequential decay processes. Individual particles that fulfil relativistic on-
shell relations $p^{\mu}p_{\mu}=m^{2}$ cannot decay into on-shell decay
products that have non-vanishing angles or energies. Colour dipoles contain
two on-shell partons, which may readily decay into three on-shell particles
(forming two decay-product dipoles) without violating momentum conservation.
Thus, the dipole decay pattern is well-defined throughout momentum space, i.e.
the dipole picture allows for the generation of physical multi-particle data.
The ability of _all_ state-of-the-art conventional showers to synthesise data
events rests on dipole-shower-inspired strategies to generate the decay
product momenta Lonnblad:1992tz ; Gieseke:2003rz ; Sjostrand:2004ef ;
Schumann:2007mg ; Platzer:2009jq ; Nagy:2014mqa ; Hoche:2015sya ;
Fischer:2016vfv ; Hamilton:2020rcu .

Colour-charged states predominantly decay through (soft) gluon emission, the
inclusive decay probability of which is given by the eikonal interference
pattern

$\displaystyle
d\mathcal{P}\left(q(p_{\mathrm{I}})\bar{q}(p_{\mathrm{K}})\rightarrow
q(p_{i})g(p_{j})\bar{q}(p_{k})\right)\simeq\frac{ds_{ij}}{s_{\mathrm{IK}}}\frac{ds_{jk}}{s_{\mathrm{IK}}}C\frac{\alpha_{s}}{2\pi}\frac{2s_{ik}}{s_{ij}s_{jk}}\frac{d\phi}{2\pi}\simeq\frac{ds_{ij}}{s_{\mathrm{IK}}}\frac{ds_{jk}}{s_{\mathrm{IK}}}C\frac{\alpha_{s}}{2\pi}\frac{2s_{\mathrm{IK}}}{s_{ij}s_{jk}}\,,$
(1)

where $C$ is a colour-charge factor, the two-particle invariants $s_{ab}$ are
defined by $s_{ab}=2p_{a}\dot{p}_{b}$, and where the momentum conservation
condition $p_{\mathrm{I}}+p_{\mathrm{K}}=p_{i}+p_{j}+p_{k}$ implies the
relation $s_{\mathrm{IK}}=s_{ik}+s_{ij}+s_{jk}$ between pre- and post-
branching momenta. This decay process can equally be interpreted as the decay
of a dipole of invariant mass $s_{\mathrm{IK}}$ into two dipoles of invariant
masses $s_{ij}$ and $s_{jk}$, which are linked by the soft gluon. Dipole
showers employ this inclusive probability to calculate the contribution of
exclusive, fixed particle-multiplicity states $\Phi$ to an observable
$O(\Phi)$ according to the master equation

$\displaystyle\mathcal{F}_{n}(\Phi_{n},t_{n},t_{c};O)=\Delta(t_{n},t_{c})O(\Phi_{n})+\int\limits^{t_{n}}_{t_{c}}dtd\xi\frac{d\phi}{2\pi}C\frac{\alpha_{s}}{2\pi}\frac{2s_{ik}(t,\xi)}{s_{ij}(t,\xi)s_{jk}(t,\xi)}\Delta(t_{n},t)\mathcal{F}_{n}(\Phi_{n+1},t,t_{c};O)\,$
(2)

where the no-branching probability $\Delta$ is given by

$\displaystyle\Delta(t_{n},t)=\exp\left(-\int\limits^{t_{n}}_{t}dtd\xi\frac{d\phi}{2\pi}C\frac{\alpha_{s}}{2\pi}\frac{2s_{ik}(t,\xi)}{s_{ij}(t,\xi)s_{jk}(t,\xi)}\right)\,$
(3)

and where a particular choice of the functional form of $t$ and $\xi$ is
called the phase space parametrisation or momentum mapping. Finally, $t$
defines the “evolution variable”. Branchings are successively ordered in this
variable, with large-$t$ branchings occurring earlier in the evolution of the
state. Large evolution variables indicate transitions at high energy, while
low evolution variables indicate low-energy transitions. The first term in Eq.
2 models the probability of the $n$-body state $\Phi_{n}$ not decaying above
$t_{c}\simeq\Lambda_{\mathrm{QCD}}^{2}\simeq(300\mathrm{GeV})^{2}$, while the
second term describes one or more decays. This introduction highlights that
the development of a parton shower relies on $a)$ the choice of inclusive
decay probabilities $d\mathcal{P}$, $b)$ the choice of an evolution variable
$t$, and $c)$ the choice of a momentum mapping $s_{ij},s_{jk}\leftrightarrow
t,\xi$ which determines the relations between pre-and post-decay momenta.

It is worth noting that all conventional state-of-the-art parton showers use
slight variations of a single algorithm – the “veto algorithm” – to solve Eq.
2 numerically. This algorithm treats the variables $t$ and $\xi$ as continuous
degrees of freedom. It is thus unsuitable for (current) quantum devices. The
following section will develop other algorithmic solutions of Eq. 2, guided by
keeping in mind the feasibility of NISQ devices.

###  2.1 Reinterpreting classical parton shower algorithms as random walks

This section extends the classical shower algorithm toolbox by performing
several abstractions of the features of dipole showers. We are led to conclude
that the showering process can be described by creating and sampling from a
fixed set of primitive fractal structures, followed by a translation of the
chosen primitive structure into scattering event momenta. The first step has
an elegant implementation on intermediate-scale quantum devices.

Figure 1:  The phase space of effective gluon emission is discrete, since  A]1
gluons within a rapidity region $\delta y_{g}$ act coherently due to running-
coupling effects. The $\kappa$ (or equivalently the $k_{\perp}^{2}$) dimension
is also quantised, since  A]2 additional phase space folds opening due to
gluon emission are quantised in units of $\delta y_{g}$. See main text for
more details.

The first abstraction to consider is removing the independent treatment of
decay probability and momentum-space integration by absorbing the non-uniform
probability density in Eq. 1 into the integration measure. This can be
obtained by choosing a phase-space parametrisation in terms of the gluon’s
transverse momentum,

$\displaystyle
k_{\perp}^{2}=\frac{s_{ij}s_{jk}}{s_{\mathrm{IK}}}\quad\textnormal{and
rapidity}\quad y=\frac{1}{2}\ln\left(\frac{s_{ij}}{s_{jk}}\right)\,,$ (4)

which leads to

$\displaystyle
d\mathcal{P}\left(q(p_{\mathrm{I}})\bar{q}(p_{\mathrm{K}})\rightarrow
q(p_{i})g(p_{j})\bar{q}(p_{k})\right)\simeq=\frac{C\alpha_{s}}{\pi}d\kappa
dy\,\quad\textnormal{with}\quad\kappa=\ln\left(k_{\perp}^{2}/\Lambda^{2})\right)\,,$
(5)

where $\Lambda^{2}$ is an arbitrary mass scale. Within this phase space
parametrisation, allowed dipole decays are constrained to a triangular region
of height $L=\ln(s_{\mathrm{IK}}/\Lambda^{2})$ in the $(y,\kappa)$-plane, as
illustrated by the left-hand panel of Fig. 1. Due to the colour charge of an
emitted gluon, the rapidity span for subsequent dipole decays (at lower
$\kappa$) is increased with respect to the originally allowed range. This
feature of QCD can be interpreted as “folding out” smaller triangular regions
of additional phase space from the original triangle. Since these regions may
again contribute to allowed decays, this process may repeat, leading to a
fractal structure of triangles-attached-to-triangles. An example structure is
shown in the upper left panel of Fig. 2. This fractal picture is the starting
point of successful conventional dipole showers such as Ariadne
Lonnblad:1992tz ; Gustafson:1987rq .

These abstractions already simplify the treatment of parton showering since
complicated “splitting functions” have been subsumed into the phase space
parametrisation. Moreover, for fixed coupling $\alpha_{s}$, the inclusive
decay rate is uniformly distributed in the $(y,\kappa)$-plane, allowing for
straightforward sampling algorithms. Nevertheless, this continuous-variable
dipole decay picture is not yet suited for current universal quantum devices.

Examining the effect of a transverse-momentum-dependent running coupling – as
supported by higher-order QCD calculations – provides a path to an even
simpler picture. We may write

$\displaystyle\alpha_{s}(k_{\perp}^{2})=\frac{12\pi}{33-2n_{f}}\frac{1}{\ln\,(k_{\perp}^{2}/\Lambda_{\mathrm{QCD}}^{2})}\,,$
(6)

where the $n_{f}$-dependent term arises through $g\rightarrow q\bar{q}$
splitting. Neglecting the latter and combining the expression with Eq. 5 leads
to

$\displaystyle
d\mathcal{P}\left(q(p_{\mathrm{I}})\bar{q}(p_{\mathrm{K}})\rightarrow
q(p_{i})g(p_{j})\bar{q}(p_{k})\right)\simeq=\frac{d\kappa}{\kappa}\frac{dy}{\delta
y_{g}}\quad\textnormal{with}\quad\delta y_{g}=\frac{11}{6}\,,$ (7)

and where we have used that in the leading-colour limit $C\rightarrow
C_{\mathrm{A}}/2=3/2$ for any dipole decay. As argued in Andersson:1995jv ,
interpreting the running-coupling renormalisation group equation as gain-loss
equation means that gluons within a rapidity range $\delta y_{g}$ act
coherently as one “effective gluon”, as illustrated by  A]1 in Fig. 1. Thus,
the rapidity range of each triangular phase space region is quantised into
multiples of $\delta y_{g}$. Since the baseline of an additional triangle
extends to positive $y$ by $\ell/2$, the height $\ell$ at which to emit
effective gluons, i.e. their $\kappa$ value, is quantised into multiples of
$2\delta y_{g}$. Thus, we may model the parton shower by generating effective
gluons at the centre of discrete tiles covering the phase-space triangle.
These realisations form the basis of the “Discrete QCD algorithm” of
Andersson:1995jv .

Each rapidity slice can be treated independently of any other slice. Inserting
Eq. 7 into the exclusive decay probability (the second term in the shower
master equation, Eq. 2) shows that the exclusive rate for finding an effective
gluon in a fixed $y$-bin takes the straightforward form

$\displaystyle
dtd\xi\frac{d\phi}{2\pi}C\frac{\alpha_{s}}{2\pi}\frac{2s_{ik}(t,\xi)}{s_{ij}(t,\xi)s_{jk}(t,\xi)}\Delta(t_{n},t)=\frac{d\kappa}{\kappa}\exp\left(-\int\limits^{\kappa_{max}}_{\kappa}\frac{d\bar{\kappa}}{\bar{\kappa}}\right)=\frac{d\kappa}{\kappa_{max}}$
(8)

Thus, for a fixed $y$-bin, the $\kappa$-value of an effective gluon in Fig. 1
is, to first approximation, simply given by $(\textnormal{number of
tiles})^{-1}$. This is significantly simpler than in conventional showers,
which rely on sampling the no-emission probability $\Delta$ with the “veto
algorithm”.

Figure 2:  The Discrete QCD parton shower algorithm can be re-interpreted as a
one-dimensional random walk, since  A]3 the baseline of the folded structure
carries all necessary information. The “grove-like” baseline structure can
A]A be generated by a heavily constrained two-dimensional random walk. Due to
the low fractal dimension of the grove structure, a one-dimensional random-
walk algorithm  A]B is equally viable. For  A]B, the notation $n/n_{max}$
indicates that option $n$ was picked out of $n_{max}$ choices. The two-
particle invariants ($\ln(s_{ij})$) can be read off by following the path from
particle $i$ to particle $j$, and skipping segments whose colour was created
and reabsorbed along the way (e.g. skipping the pink segments when calculating
the invariant mass of the green and blue gluon tips).

Once an effective gluon has been selected, new triangles are folded out of the
parent region. Effective gluon positions in this fold are again quantised into
tiles of dimensions $\delta y_{g}\times 2\delta y_{g}$, as illustrated by the
upper left part of Fig. 2. However, a simpler interpretation is possible since
the height and the $y-$range of each fold are redundant: All information
necessary to calculate momentum invariants $s_{ij}$ can be read off the
baseline of the folded triangle structure, shown in the lower left part of
Fig. 2. We will call a specific baseline structure a “grove”. The shortest
distance (along the baseline) between two “tips” $i$ and $j$ can be shown to
equal $\ln(s_{ij}/\Lambda^{2})$. Together with the knowledge of the overall
centre-of-mass energy and uniformly sampled azimuthal decay angles $\phi$,
this information is sufficient to construct post-decay kinematics.

The Discrete QCD algorithm allows a simple method to produce groves with
correct rates. However, there are many ways to create the grove structures
apart from the Discrete QCD algorithm. One example is shown in the lower right
part of Fig. 2. Since the grove structure is contained in a triangular region
smaller than the original background triangle, and points on the baseline are
separated from their nearest neighbours by $\delta y_{g}$, a two-dimensional
fixed step-length random walk may be used to produce the grove structures.
However, such a two-dimensional random walk is heavily constrained since it
cannot extend above arbitrarily high or small $\kappa$ values, and the
$y$-value of each step cannot decrease. These constraints suggest that the
most straightforward algorithm to create the grove structure with correct
probabilities is a one-dimensional random walk with variable step length. The
step length would be determined by an auxiliary $y$-bin dependent choice,
which considers that small $|y|$ values allow for fewer $\kappa$ values. An
example of a one-dimensional random walk is given in the upper right part of
Fig. 2.

The random-walk algorithm can be summarised by

Start at the leftmost point on the baseline. Mark the position as quark.

Move to the next neighbour along the baseline with unit probability.

while  _walker has not reached the rightmost point on the primary baseline_ do

Calculate the number of tiles $n_{max}$ in this rapidity slice. Choose a new
baseline

height $h=n\delta y_{g}$, with $n\in\\{0,1,...,n_{max}-1\\}$ and probabilities
$P_{n}=1/n_{max}$.

if  _$h=0$_ then

move to the next neighbour along on current baseline with unit probability.

else

“Step into the fold”, i.e. take one step of length $\delta y_{g}$ along a new
baseline111Indicated by changing the “line colour” in upper right part of Fig.
2 . .

if  _the position after the step $=h$_ then

Mark the position as effective gluon222The position $h$ marks the “tips” of
the gluon fold in the upper left part of Fig. 2 . Positions $<h$ correspond to
one side of the gluon fold, positions $>h$ to the other side of the fold.

Add one step along the current baseline with unit probability.

else if  _the position after the step $=2h$_ then

“Step out of the fold”333Indicated by reverting to the previous colour in
upper right part of Fig. 2 . and move to next neighbour on the previous

baseline with unit probability.

Mark the rightmost point as antiquark.

This very compact algorithm is inspired by considering the constraints of a
quantum realisation on current NISQ devices. The moderate quantum volume of
these devices favours a less resource-intensive one-dimensional approach. This
compact classical parton shower algorithm is an essential finding of this
article and shows potential synergies between classical and quantum algorithm
development.

Note that all possible groves can be enumerated for fixed centre-of-mass
energy since the $(y,\kappa)$ plane is discretised in terms of $\delta y_{g}$.
This allows pre-computing all grove structures and their rates (e.g. with the
random walk algorithm above) before simulating event data from the information
stored in the grove. Event generation then proceeds by picking a grove
structure according to its probability and setting up the particle momentum
vectors defined by the grove. Such an algorithm is straightforward and
efficient. Moreover, during the event generation step, it is unnecessary to
know “how” the groves and their rates have been produced. In fact, we can use
an elegant quantum algorithm to achieve grove generation.

###  2.2 Generating scattering events from the grove structure

Event data can be generated once a grove structure has been selected. Within
the classical algorithm we use as a baseline for the quantum algorithm below,
we create post-decay momenta iteratively from the grove structure, creating
the highest-$\kappa$ effective gluons first. For each effective gluon $j$ that
has been emitted from a dipole $\mathrm{IK}$, we read off the values of
$s_{ij}$, $s_{jk}$ and $s_{\mathrm{IK}}$ from the grove, generate a uniformly
distributed azimuthal decay angle $\phi$, and then employ the momentum mapping
of Gehrmann-DeRidder:2011gkt  to produce post-branching momenta444Other
momentum mappings would, of course, be possible. Our choice reflects that our
dipole-shower algorithm does not distinguish individual partons, leading us to
believe that a momentum mapping that does not differentiate individual partons
is suitable.

There is some degree of freedom when reading off $s_{ij}$ and $s_{jk}$ from
the grove structure since the latter only defines which $(y,\kappa)$-tile (in
Fig. 1) is picked. In fact, as argued below Eq. 7, each point within the
$(y,\kappa)$-tile is equally likely. This freedom can be used to mitigate the
most glaring effects of phase-space discretisation. We follow the strategy of
Andersson:1995jv  and distribute the $y$ and $\kappa$-values of the effective
gluon in the highest and second-highest $\kappa$ tiles uniformly within the
tile. This ensures a good model for the highest-$\kappa$ gluons.

The simple grove selection probabilities rely on the fact that the
$(y,\kappa)$ plane is uniformly covered, by virtue of Eq. 5, making the
selection probabilities proportional to the area of the $(y,\kappa)$-tiles. If
a tile protrudes beyond the allowed phase space boundaries, its area is
reduced Andersson:1995jv . We incorporate this correction by weighting each
event with a factor

$w=\prod_{i_{s}}\frac{(\textnormal{area of tile picked in slice
$i_{s}$})\,/\,(\textnormal{area of slice $i_{s}$})}{(\textnormal{number of
tiles in slice $i_{s}$})^{-1}}~{}.$ (9)

##  3 The Quantum Parton Shower Algorithm

This section outlines the generation of the grove structures for the parton
shower described in Sec. 2.1 using a quantum device. The quantum algorithm has
been constructed using the quantum walk framework PhysRevA.48.1687 ; QWProc ;
Kempe . The quantum analogue of the classical random walk, the discrete-time
quantum walk defines the movement of a walker in position space,
$\mathcal{H}_{P}$, dictated by the result of a coin operation,
$\mathcal{H}_{C}$, such that the total Hilbert space of the quantum walk is

$\mathcal{H}=\mathcal{H}_{P}\otimes\mathcal{H}_{C}.$ (10)

A single step in the quantum walk is constructed from two components: the coin
flip operation, $C$, which encodes the coin state in $\mathcal{H}_{C}$ onto a
coin register; and the shift operation, $S$, which controls from the coin
register and moves the walker in $\mathcal{H}_{P}$ accordingly. The evolution
of the walk is entirely unitary and can thus be expressed as

$U=S\cdot(C\otimes\mathbb{I}),$ (11)

where $U$ is a unitary matrix that describes a single step of the walk in
$\mathcal{H}$, such that an $N$-step quantum walk is defined by $U^{N}$. The
quantum walk has drastically different characteristics from a classical random
walk due to possible interference effects introduced in the coin operation. A
thorough review of the quantum walk framework is available in Reference Kempe
and references therein.

Quantum walks are compact algorithms which have the potential to reduce the
Quantum Volume required on the device Williams:2021lvr . Consequently, the
quantum walk architecture is well suited to designing algorithms for NISQ
devices. Such quantum devices have a small number of controllable qubits and
require shallow circuit depths to limit decoherence effects on the device.
Appendix A discusses how algorithms can be streamlined to get practical
results from NISQ devices, highlighting some difficulties associated with such
devices.

Of particular interest to a quantum algorithm for parton shower simulation is
the quantum walk with memory McGettrick2010 ; Shakeel2014 . Quantum walks with
memory can simulate arbitrary dynamics by modifying the movement of the walker
based on the memory of the previous position Camilleri2014  and coin
operations PhysRevA.67.052317 ; PhysRevA.87.052302 . The amount of memory that
the quantum walk has affects the diffusive characteristics of the quantum
walk: from an ideal quantum distribution in the limit of minimal memory to an
ideal classical distribution in the limit of full memory PhysRevA.67.052317 .
Importantly, no decoherence effects are present within the algorithm, and the
evolution is entirely unitary, even within the limit of full memory. Moreover,
by reversing time, one can regain the initial state of the quantum walk, which
is impossible from a genuinely classical algorithm PhysRevA.87.052302 .
Therefore, quantum walks with memory provide the opportunity to simulate
arbitrary dynamics, including classical statistics, on a quantum device. For
this reason, quantum walks with memory are of great interest in computing
PhysRevA.93.042323 ; Roget2020  and high energy physics Williams:2021lvr .

In contexts outside of high energy physics, quantum walks have been shown to
exhibit speed-up over their classical counterparts for certain cases.
Simulated annealing PhysRevLett.101.130504 ; PhysRevA.78.042336  and search
algorithms PhysRevA.67.052307 ; Szegedy  are prominent examples, and it is
hoped that quantum walks will provide a speed-up for generic Markov Chain
Monte Carlo (MCMC) algorithms Montanaro . Improved speeds have been shown for
dedicated MCMC algorithms under specific conditions Montanaro ;
Lemieux2020efficientquantum . The critical characteristic that leads to the
speed of an MCMC algorithm is the Markov Chain mixing time, the time it takes
for the MCMC algorithm to reach the equilibrium distribution Levin2008 .
Quantum algorithms can result in a quadratic speed-up over classical analogues
in reaching equilibrium PhysRevA.78.042336 ; PhysRevA.76.042306 , though
proving that quantum walks lead to speed gains for MCMC algorithms in general,
is still an active field of research PhysRevA.76.042306 ;
Orsucci2018fasterquantummixing ; PhysRevA.104.032215 . Parton shower
algorithms are related to, but due to the state memory implied by shower
ordering conditions distinct from, MCMC algorithms. This makes developing a
quantum walk parton shower also interesting more abstractly: as an MCMC-
related class of algorithms that could benefit from quadratic speed-up. This
article does not consider an assessment of potential speed-ups, as we focus on
synthesising realistic data as a use case.

###  3.1 Implementation on a quantum device

Figure 3:  Schematic of the quantum Discrete QCD parton shower algorithm
circuit. The algorithm is a quantum walk with memory, constructed from maximum
five operations per step: the coin operation $C$, baseline shift $B$, the
$\lambda$ shift $S$, the gluon shift $S^{\prime}$, and the memory operation
$M$.

With their ability to simulate arbitrary dynamics and interference effects
whilst also having the potential to speed up MCMC algorithms, quantum walks
are an exciting test bed for developing quantum algorithms for the simulation
of parton showers. Here, we present a new approach to simulating a parton
shower on a quantum device, obtaining results comparable to data from collider
experiments. The quantum algorithm simulates the classical parton shower
outlined in Sec. 2.1. The algorithm is constructed from a quantum walk with
one dimension of memory (i.e. the memory register is constructed from one
memory qubit), where the Hilbert space of the walker is formed from the grove-
baseline ($\lambda$) space, $\mathcal{H}_{\lambda}$, defined in Fig. 2, and
the effective gluon space, $\mathcal{H}_{g}$, augmented by the coin space,
$\mathcal{H}_{C}$. A schematic of the quantum circuit is shown in Fig. 3. The
algorithm is constructed from a maximum of five operations per step: the coin
flip operation $C$, the baseline operation $B$, the $\lambda$-shift operation
$S$, the gluon-shift operation $S^{\prime}$, and the memory operation $M$.

The movement of the walker in the $\lambda$-space is schematically shown in
the bottom right of Fig. 2. Following the algorithm’s circuit shown in Fig. 3,
the walker starts in the left-most position of the initial fold. The coin-flip
operation $C$ constructs an equal state on the coin register, such that the
outcome of the coin gives an equal probability of selecting a tile in the
first slice, see App. A. As the walker will always move to the right every
step of the walk, the baseline operation $B$ is applied, irrespective of the
coin outcome, to increment the walker’s position in $\mathcal{H}_{\lambda}$.
The shift operations, $S$ and $S^{\prime}$, then control from the coin
register to move the walker correctly in $\mathcal{H}_{\lambda}$ and
$\mathcal{H}_{g}$ respectively, depending on the outcome of the coin. For
example, suppose the selected tile is at height $h=n\delta y_{n}$. In that
case, the walker is moved $2n$ steps in $\mathcal{H}_{\lambda}$, corresponding
to the walker moving along the baseline of the secondary fold created at $h$
(moving vertically down and then up in Fig. 2). The walker is then shifted in
$\mathcal{H}_{g}$ to simulate the emission of an effective gluon correctly.
The final operation of the step is the memory operation $M$. This allows for a
conditional coin to be applied to the register depending on the outcome of the
previous coin operation, equivalent to recycling the coin. The memory register
is constructed from a single qubit; thus, the walker only has a memory of the
previous coin operation. Furthermore, the memory operation is only needed once
for the simulation of events with $E_{\mathrm{CM}}=91.2~{}\mathrm{GeV}$, such
as events from LEP experiments, thus retaining quantum diffusive
characteristics. The algorithm is then repeated for any additional folds
created, omitting the baseline operation, which is accounted for in the parent
fold.

At the end of the algorithm, the gluon and $\lambda$ registers can be measured
to obtain a grove primitive. This can then be passed to the event generation
algorithm described in Sec. 2.2 to produce pseudodata from the grove structure
information returned from the quantum device. It should be noted that, for
exceptional cases, such as the example discussed here, the algorithm can be
streamlined by removing the $\mathcal{H}_{\lambda}$ space. Tailoring the
algorithm to simulate data from LEP experiments reduces the algorithm’s
extendability but decreases the circuit depth. Therefore, the algorithm can be
run on NISQ devices with higher fidelity than the full algorithm. A detailed
explanation of this modification is given in App. A. The ibm_algiers device
still returns noisy results despite using the streamlined algorithm.

###  3.2 Grove generation on a quantum device

(a) Probabilities of grove structures.

(b) Labelling of primitive grove structures.

Figure 4:  Generation of primitive grove structures for
$E_{\mathrm{CM}}=91.2~{}\mathrm{GeV}$ from the 32-qubit ibmq_qasm_simulator
and the 27-qubit ibm_algiers device, compared to the analytically calculated
classical distribution.

The parton shower algorithm presented in Sec. 3.1 is the first of its kind to
be comparable to real, archival collider data. To demonstrate this, the
algorithm has been constructed to simulate the evolution of a $q\overline{q}$
final state produced at $E_{\mathrm{CM}}=91.2~{}\mathrm{GeV}$, typical of
events at the LEP collider. For this centre-of-mass energy, there are a total
of 24 primitive grove structures which can be generated, with tertiary folds
being the maximum fold depth. The circuit consists of four registers, as shown
in Fig. 3, and requires 15 qubits. The circuit is compact, using 116 gate
operations (102 multi-qubit and 14 single qubit operations), and thus is
suited to NISQ devices.

Figure 4 shows the probability distribution for grove generation from the
32-qubit ibmq_qasm_simulator and the 27-qubit ibm_algiers quantum computer,
compared to the analytically calculated classical distribution. The simulator
has been run for 100,000 shots and shows good agreement with the classical
distribution. The ibmq_qasm_simulator simulates a fully fault tolerant device
and thus does not simulate noise effects or readout errors in NISQ devices.
Consequently, we get an exact agreement with the classical algorithm by
feeding the grove structures created by the simulator into the event
generation step outlined in Sec. 4. This comparison is shown in Fig. 5.

The ibm_algiers quantum computer has been run for 20,000 shots, using the
tailored circuit outlined in Sec. 3.1 and App. A. This reduces the number of
qubits in the circuit to 10 qubits and reduces the required number of gate
operations from 116 gate operations to 21 gate operations (12 multi-qubit and
nine single qubit gate operations). This dramatic decrease in circuit depth is
crucial for obtaining high-fidelity results from NISQ devices. Figure 4 shows
the uncorrected performance of the ibm_algiers device compared to the
ibmq_qasm_simulator for generating the 24 grove structures for
$E_{\mathrm{CM}}=91.2~{}\mathrm{GeV}$. The ibm_algiers device returns noisy
results, which, for some grove structures, greatly differ from the simulator.
This algorithm’s main source of error is a large amount of swap operations
needed to correctly implement the gluon shift, resulting in many multi-qubit
gate operations. Multi-qubit gates, such as the cnot gate, can introduce
environmental noise to the device, resulting in decoherence. Consequently, the
many swap operations needed in the computation directly affect the fidelity of
the results. Furthermore, the results tend to prefer states with fewer
effective gluons. This is likely due to the relaxation time of the qubits in
the gluon register, causing the qubits to return to the ground state before
measurement. A detailed discussion of the noise effects and how these affect
the results from the ibm_algiers device is given in App. B.

It is possible to employ quantum error mitigation schemes to suppress readout
errors in the results. However, it is shown in Sec. 4 that mitigation schemes
are not needed for event generation, and the algorithm is remarkably robust to
noise effects.

Figure 5:  Representative comparison the (classical) Discrete QCD model and
Quantum Parton Shower results from the 32-qubit ibmq_qasm_simulator, for the
transverse momentum and rapidity of all final-state particles in the synthetic
data. As expected from Figure 4, the two algorithms match exactly.

##  4 Event Generation using the Quantum algorithm

In this section, we compare the result of using the ibm_algiers device to
create data and compare the measurements of selected kinematic distributions
performed on these simulated data with archival collider data. We aim to
highlight different aspects of the simulation and thus compare both event
shape and jet observables. The value of the _thrust_ observable $T$
ALEPH:2003obs  classifies the event into “pencil-like” ($T\sim 1$) and
“planar-circular” events ($T\sim 0.5$). Thus, values of $T\rightarrow 1$
provide an excellent probe of the modelling of soft- and collinear parton
emissions in the parton shower, and the transition between partons and hadrons
for $T\sim 1$, while moderate $T$-values probe high-energy gluon emissions.
The _energy-energy correlation_ (EEC) is another event shape observable that
probes QCD evolution DELPHI:1996sen . Suppose the high-energy primary partons
mostly retain their direction (as would be the case for the asymptotically
free limit of QCD). In that case, one expects a particle at angle $\xi$ to be
accompanied by another particle at angle $\xi+\pi$. As the parton shower
evolution leads to collinear radiation, one expects several collimated
particles around a test particle at an angle $\xi$. At intermediate EEC values
from $-\frac{1}{2}\lesssim\cos\xi\lesssim\frac{1}{2}$, we expect soft-gluon
radiation as well as strong-coupling (non-perturbative) QCD effects to
contribute to form a plateau. The _Durham 2-jet rate_ $y_{23}$ JADE:1999zar
defines a “closeness” measure of the second and third-hardest jets in the
collision event666At LEP, each event is expected to have at least two jets,
producing two high-energy primary partons in the original
$e^{+}e^{-}$-annihilation.. Well-separated jets, meaning jets with high energy
and large angular separation, lead to high $y_{23}$ values. At such values,
higher-order real-emission corrections dominate the spectrum. At the
intermediate 2-jet rate, the jets start to coalesce. This region is
susceptible to parton shower evolution. Finally, at very small $y_{23}$
values, the jets are no longer separated, and non-perturbative (hadronisation)
effects dominate the distribution. Finally, the _jet mass difference_ $M_{D}$
ALEPH:2003obs  measures the (normalised) invariant mass of particles in the
left and right detector hemispheres. Small values of $M_{D}\sim 0$ correspond
to equally populated hemispheres, as would be the case for perfectly balanced
jets, while values of $M_{D}\lesssim\frac{1}{3}$ probe the modelling of the
hard- and multiple soft-gluon emission. Large $M_{D}$-values probe multiple
hard gluon emissions.

We embed the Quantum Parton Shower into a classical toolchain to facilitate
comparisons to data from the Large Electron-Positron collider (LEP). The
momentum distribution of quarks and antiquarks in the annihilation
$e^{+}e^{-}\rightarrow Z/\gamma^{*}\rightarrow q\bar{q}$ at
$E_{\mathrm{CM}}=91.2\,\mathrm{GeV}$ is generated classically. The parton
shower evolution of the $q\bar{q}$ final state is handled by the NISQ device,
which generates grove structures. We run the shower evolution algorithm on the
ibm_algiers device, which operates a Falcon r5.11 architecture. The grove
structures are then used to generate the desired kinematics, and the result is
stored in Les Houches Event (LHE) File format Alwall:2006yp ; Andersen:2014efa
. These LHE files are passed to the Pythia event generator Bierlich:2022pfr
for parton-to-hadron conversion by means of the Lund string model
Andersson:1983ia . Finally, we use Rivet analysis framework Buckley:2010ar ;
Bierlich:2019rhm  to compare our synthetic data to data taken with LEP
detectors.

We do not embark on a dedicated program to “tune” the parameters of the
simulation Buckley:2011ms  to improve the data description. The parton shower
mass scale $\Lambda$ is chosen as $\Lambda=0.4$ GeV since such a value allows
for a reasonably rich structure of possible groves at LEP collision energies.
As mentioned in Sec. 2.2, for the highest and second-highest $\kappa$ values
we distribute the y and $\kappa$ values for the effective gluons uniformly
within the tile. The parameters of the Lund symmetric fragmentation function
Andersson:1983ia  are set to $a=2$ and $b=3/5$, and width of the non-
perturbative string-$p_{\perp}$ is set to $\sigma=0.1$ GeV777A low value of
$\sigma$ is motivated by the fact that the mass scale $\Lambda=0.4$ GeV also
sets the shower cut-off. Thus, the shower covers lower $p_{\perp}$-values than
conventional methods and thus requires a smaller non-perturbative component to
the $p_{\perp}$ distribution..

Sample comparisons to experimental data are shown in Fig. 6. We find very
satisfactory agreement of the classical Dqcd shower results, in particular
given the simplicity of the shower model. As seen in the rightmost bins of the
jet mass difference $M_{D}$ and the Durham 2-jet rate $y_{23}$, the algorithm
tends to overpopulate states with multiple high-energy gluons. This is
expected since the soft-gluon emission rate overshoots the true matrix
elements in the hard region. Similarly, we expect that improving the
description outside the soft limit will improve the description at
intermediate thrust values $T\sim 0.75-0.9$. Both high $T$-values and moderate
EEC values indicate that the emission of soft gluons is well-described.

Comparing the results of running the quantum algorithm on the ibm_algiers
device to experimental data is again very favourable. Interestingly, the
device’s tendency to produce “too few” excited (multi-gluon) states slightly
improves the overall data description. Error correction would, of course,
mitigate this effect. Nevertheless, the sound data description of the
uncorrected quantum results may hint at improvements in the classical
algorithm.

Figure 6:  Sample comparisons of the Discrete QCD model and the Quantum Parton
Shower to data taken at the LEP collider ALEPH:2003obs ; DELPHI:1996sen ;
JADE:1999zar . The Quantum Parton Shower results are not corrected for errors
in the qubit evaluations.

##  5 Summary

In this article, we have synthesised realistic particle collision events by
sampling parton shower configurations on a quantum device. As an application,
we have compared the simulated data to data recorded at experiments at the
Large Electron-Positron (LEP) collider, finding favourable agreement.

To obtain these results, we have re-interpreted conventional parton showers –
which rely on an intricate form of rejection sampling – as one-dimensional
random walks in the $\lambda$-measure Andersson:1988ee , i.e. along the
baseline of the multigluon phase-space structure. The novel algorithm is a
reformulation of the Discrete QCD algorithm Andersson:1995jv . It allows
selecting the result of the parton shower evolution before multiparticle
momenta are generated while ensuring on-shell conditions and momentum
conservation. Developing this unique strategy was necessary to set a suitable
starting point for quantum algorithm development. Furthermore, the generation
and sampling of parton shower results can be performed on a quantum device,
whose outputs are then used to reconstruct kinematics and ultimately perform
comparisons to data.

This is the first time that the result of a quantum algorithm has been
compared to “real-life” particle physics data. The quantum algorithm is
constructed using the quantum walk framework and, consequently, is a compact
algorithm with a short circuit depth, an important aspect to consider to
obtain practical results from NISQ devices. We synthesise data using the
ibm_algiers device, which uses an IBM Falcon r5.11 processor. Figure 4 shows
that the results from the ibm_algiers device are subject to a non-negligible
amount of noise. The main noise source has been identified as cnot errors due
to gate decompositions and swap operations implemented in the transpilation
process888Transpilation is the process of rewriting a given input circuit to
match the topology of a specific quantum device, and/or to optimise the
circuit for execution on present day noisy quantum systems.. In addition to
cnot errors, it can be seen from the results in Fig. 4 and 6 that the
algorithm prefers final states will fewer soft emissions. From a close study
of the device output, this error has been attributed to the T1 (relaxation
time) of the qubits in the gluon register, causing decoherence during the
algorithm’s run time. A detailed discussion of the errors on the device has
been provided in App. B. Cumulatively, these errors result in a device output
far from ideal. However, since physical observables are supported by an
admixture of several parton shower structures, our results are remarkably
robust to noise, as demonstrated by the data comparison in Fig. 6. Overall,
our current quantum algorithm improves on previous quantum shower attempts in
several physics aspects (the inclusion of soft-gluon coherence, running-
coupling effects and realistic scattering data generation) as well as on
computational aspects: the required quantum volume to simulate data from LEP
experiments is less than that needed for all previously known quantum shower
algorithms Bauer:2019qxa ; Bepari:2020xqi ; Williams:2021lvr . Furthermore,
the necessity for fully fault tolerant devices has been drastically reduced.

A quantum event generator would also include quantum algorithms to sample the
highest-momentum transfer interactions and quantum algorithms for non-
perturbative hadronisation models. Progress on the former has been reported in
Bepari:2020xqi , while Klco:2018kyo  reported on modelling string breaking in
a simplified model on a NISQ device. Combined with these advances, our results
indicate that a programme of developing quantum event generators on near or
intermediate time scales is within reach. This programme should also include
several improvements to our quantum shower algorithm. An extension of the
classical algorithm to include $g\rightarrow q\bar{q}$ splittings, hard-
collinear corrections, and subleading colour terms may be inspired by previous
work Gustafson:1992uh , though detailed studies are necessary. A programme of
dedicated event generator parameter “tuning”, as is typical for conventional
event generators, should also be initiated for the quantum shower algorithm.
To obtain competetive predictions, appreciable work will be necessary to
incorporate precision QCD calculations into the Discrete QCD framework through
matching Frixione:2002ik ; Nason:2004rx ; Frixione:2007vw  or merging
Catani:2001cc ; Lonnblad:2001iq ; Mangano:2001xp ; Mrenna:2003if ;
Alwall:2007fs  methods.

Notwithstanding these challenges, it is exciting to see that the rapid
advances in algorithm design on NISQ devices allow us to think about such
technicalities already today.

## Acknowledgments

We thank Sarah Malik for valuable discussions. We acknowledge the use of IBM
Quantum services for this work and to advanced services provided by the IBM
Quantum Researchers Program. The views expressed are those of the authors, and
do not reflect the official policy or position of IBM or the IBM Quantum team.
We would like to specifically thank the IBM Support Centre for their help with
running on the ibm_cloud.

##  Appendix A Tailoring of the Quantum Implementation

Figure 7:  Schematic of the reduced quantum Discrete QCD parton shower
algorithm circuit. The algorithm is a quantum walk with memory, constructed
from maximum three operations per step: the coin operation $C$, the gluon
shift $S^{\prime}$, and the memory operation$M$.

Noisy Intermediate-Scale Quantum (NISQ) devices have a limited number of
qubits and, perhaps more importantly, require shallow circuit depths to return
high fidelity results. Furthermore, current quantum error mitigation schemes
are not to a high enough standard to suppress all errors from the device, such
as read-out and cnot errors. Therefore, to run practical algorithms on current
NISQ devices, it is essential to write efficient quantum algorithms which can
minimise the required quantum volume on the device, thus reducing noise
throughout the calculation.

The general and extendable quantum algorithm for the simulation of parton
showers outlined in this paper is the first quantum algorithm with the ability
to simulate a full parton shower with kinematics. The algorithm is implemented
using the quantum walk framework, allowing for compact and efficient circuit
implementation, as shown in Fig. 3. For a centre of mass energy typical of
collisions at the LEP collider, the algorithm can be run on the 27-qubit
ibm_algiers device and is constructed from 15 qubits and 116 gate operations
(102 multi-qubit and 14 single qubit gate operations). However, ibm_algiers is
a NISQ device, and even the compact circuit depth of the parton shower
algorithm introduces significant noise effects and produces low fidelity
results. For the specific example discussed in this paper, it is possible to
streamline the algorithm at the expense of the extendability of the algorithm
but at the benefit of drastically reducing the quantum volume required.

Removing the $\lambda$-walk in $\mathcal{H}_{\lambda}$ from the quantum walk,
restricting the walker to only $\mathcal{H}_{g}$ augmented by the coin space,
$\mathcal{H}_{C}$, allows for the circuit to be reduced in size to 10 qubits
and 21 gate operations (12 multi-qubit and nine single qubit gate operations).
Crucially, the number of multi-qubit gate operations is decreased by a factor
of approximately five, resulting in a decrease in the cnot errors within the
device. A schematic of this reduced circuit is shown in Fig. 7.

This algorithm has been run for 20,000 shots on the ibm_algiers device through
the ibm_cloud service. The graph for the grove generation is shown in Fig. 4,
and highlights the errors on the device compared to the output from the
ibmq_qasm_simulator. The error on the device remains large, despite the use of
the streamlined algorithm. Section 4 discusses the event generation using the
distribution of primitive grove structures produced by the quantum device.
Remarkably, the event generation algorithm is robust to noise on the quantum
device.

The streamlined algorithm uses a quantum walk with memory framework, discussed
in Sec. 3. For a centre of mass energy $E_{\mathrm{CM}}=91.2\,\mathrm{GeV}$,
the memory operation is only required once in the algorithm. This is applied
after the first step to implement the correct coin operation for the second
step. Consequently, the walker remains within the limit of quantum diffusive
characteristics.

The coin operation constructs an equal superposition of $n$ states on the coin
register of $m$ qubits, such that the outcome of the coin gives an equal
probability of obtaining each state,

$C|0\rangle^{\otimes
m}=\frac{1}{\sqrt{n}}\Big{(}|0\rangle+|1\rangle+\ldots+|n-1\rangle\Big{)}.$
(12)

This reflects the probability of selecting a tile in a slice of $n$ tiles, see
Fig. 1. As an example, consider the centre of mass energy
$E_{\mathrm{CM}}=91.2\,\mathrm{GeV}$. The phase space available is constructed
from slices of two and three tiles, requiring coin register states,

$\displaystyle\frac{1}{\sqrt{2}}\Big{(}|0\rangle+|1\rangle\Big{)},$ and
$\displaystyle\frac{1}{\sqrt{3}}\Big{(}|00\rangle+|10\rangle+|01\rangle\Big{)},$
(13)

respectively. These states can be constructed simply using the circuit
diagrams shown in Fig. 8.

Figure 8:  Coin operations for slices with $n=2$ and $n=3$ tiles.

##  Appendix B Identification of Quantum Errors

The use of Noisy Intermediate-Scale Quantum (NISQ) devices is impractical
without understanding the possible sources of quantum errors present within
the device RevModPhys.82.1155 . Unlike classical devices, quantum computers
are limited in qubit number and thus cannot rely on redundancy for error
correction FrontEng . For that reason, quantum error correction is a highly
active research field lidar2013quantum ; Devitt_2013 . However, a general
scheme for error correction, resulting in a fully fault tolerant device, does
not currently exist. Quantum devices are subject to environmental decoherence
effects, such as unwanted electric fields or hardware defects FrontEng .
Therefore, to obtain high fidelity results from a NISQ device, careful
considerations of circuit architecture have to be made, most importantly the
depth of the circuit and the number of multi-qubit operations. For this
reason, the quantum shower algorithm outlined in Sec. 3.1 has been streamlined
for the special case of $E_{\mathrm{CM}}=91.2\,\mathrm{GeV}$. The tailoring of
this special case is outlined in App. A.

The quantum parton shower algorithm has been run on the ibm_algiers devices
for 20,000 shots. A comparison to analytically calculated expected
probabilities and simulated quantum results from the ibmq_qasm_simulator is
given in Fig. 4. It is clear from the disparity between the ibm_algiers device
and the imbq_qasm_simulator that the real quantum device returns noisy
results, with the tendency to prefer grove structures with less effective
gluons produced, a visible effect in Fig. 4 and 6. The high number of multi-
qubit gate operations is the main source of error for the quantum shower
algorithm on the ibm_algiers device. Currently, on the ibm_algiers device,
multi-controlled gates such as the ccnot gate are decomposed into cnot gates.
Therefore, even if only a small amount of ccnot gates are used, including the
required swap operations, the number of cnot gates implemented can be an order
of magnitude higher. Consequently, this will greatly affect the results’
fidelity as cnot gates can introduce environmental noise leading to
decoherence.

Another primary source of error in the results from Fig. 4 comes from the
qubits’ T1 time (sometimes called the relaxation time). The T1 time of a qubit
is the time it takes the qubit to decay from the excited state, in this case,
$|1\rangle$, to the ground state, $|0\rangle$ on IBM devices QiskitT1 . The T1
time, therefore, dictates the time in which practical and meaningful
algorithms can be run. For the quantum shower algorithm, the gluon shift
records whether a gluon has been created by applying a shift in
$\mathcal{H}_{g}$ to the walker, as shown in Fig. 3. As more gluons are
created, more qubits are flipped from the $|0\rangle$ to the $|1\rangle$
state. If the first slice produces a gluon, then some qubits in the gluon
register will be flipped to the $|1\rangle$ state and can remain in this state
throughout the algorithm. These qubits could likely decay to the ground state
$|0\rangle$ whilst the algorithm is still running on the device and therefore
would lead to incorrect grove structures measured at the end of the algorithm.
This would explain the excess of events without soft emission. A thorough
overview of hardware noise and errors is given in QiskitT1  and references
therein.

Figure 9: Sample comparisons of the Discrete QCD model for varying values for
the mass scale $\Lambda$. This leads to non-negligible uncertainties for
observables that are dominated by perturbative QCD, as expected from a
leading-logarithmic model.

To mitigate these errors efficiently, an in-depth study of the transpilation
process needs to be performed using the IBM Q network IBMQ . This allows the
user to control the level of optimisation of the circuit to achieve a
compromise between the circuit depth and noise abatement. However,
sophisticated error mitigation is not required for the example presented in
this paper. Section 4 outlines the event generation from primitive grove
structures constructed using the quantum device. Figure 6 shows that the
algorithm is impressively robust to noise from the device, demonstrating the
usability and practicality of NISQ devices for exciting problems in high-
energy physics.

##  Appendix C Supplementary data comparisons

Figure 10: Comparative study for the affect of different tunes on the data
agreement. Observables dominated by non-perturbative dynamics show mild-
dependance on the mass scale $\Lambda$, but are sensitive to changes in the
tune.

The classical Discrete QCD algorithm introduced in Sec. 2.1 and 2.2 is one of
the main developments of this article. It constitutes an compact reformulation
of the orignal approach Andersson:1995jv . As such, additional data
comparisons can be helpful to assess the predictions provided by the method.

The Discrete QCD method provides an elegant model for perturbative aspects of
the showering process. Non-perturbative aspects are included by employing the
Lund string model implementation of Pythia Bierlich:2022pfr . The only
parameter of the perturbative model is the mass scale $\Lambda$. Figure 9
shows that variations of the mass scale do indeed yield non-negligible
uncertainties, as should be expected from a leading-logarithmic model. Changes
of the non-perturbative model (by varying the Pythia tune) lead to smaller
variations. Examples of observables that are dominated by non-perturbative
dynamics are shown in Fig. 10. These exhibit only very mild dependence on the
value of $\Lambda$, but are very sensitive to changes in the tune. Both the
default Pythia “Monash” tune Skands:2014pea  as well as the Dire parton shower
tune Isaacson:2018zdi  yield a better agreement with data than the nominal
parameters employed for the Dqcd shower. This confirms a dedicated tuning of
the Dqcd shower could improve the overall description of data. We have
refrained from this, to not give an artificially positive impression of the
Dqcd shower results.

## References

  * (1) G. C. Fox and S. Wolfram,  Nucl. Phys. B 168, 285 (1980).
  * (2) T. Sjostrand,  Phys. Lett. B 157, 321 (1985).
  * (3) G. Gustafson,  Phys. Lett. B 175, 453 (1986).
  * (4) L. Lonnblad,  Comput. Phys. Commun. 71, 15 (1992).
  * (5) S. Gieseke, P. Stephens, and B. Webber,  JHEP 12, 045 (2003), hep-ph/0310083.
  * (6) T. Sjostrand and P. Z. Skands,  Eur. Phys. J. C 39, 129 (2005), hep-ph/0408302.
  * (7) S. Schumann and F. Krauss,  JHEP 03, 038 (2008), 0709.1027.
  * (8) S. Platzer and S. Gieseke,  JHEP 01, 024 (2011), 0909.5593.
  * (9) Z. Nagy and D. E. Soper,  JHEP 06, 097 (2014), 1401.6364.
  * (10) S. Höche and S. Prestel,  Eur. Phys. J. C 75, 461 (2015), 1506.05057.
  * (11) N. Fischer, S. Prestel, M. Ritzmann, and P. Skands,  Eur. Phys. J. C 76, 589 (2016), 1605.06142.
  * (12) K. Hamilton, R. Medves, G. P. Salam, L. Scyboz, and G. Soyez,  JHEP 03, 041 (2021), 2011.10054.
  * (13) R. P. Feynman,  Int. J. Theor. Phys. 21, 467 (1982).
  * (14) S. P. Jordan, K. S. M. Lee, and J. Preskill,  (2014), 1404.7115.
  * (15) L. García-Álvarez et al.,  Phys. Rev. Lett. 114, 070502 (2015), 1404.2868.
  * (16) P. Arrighi, G. D. Molfetta, I. Márquez-Mártin, and A. Pérez,  Physical Review A 97, 062111 (2018), 1803.01015.
  * (17) I. Márquez-Mártin, P. Arnault, G. D. Molfetta, and A. Pérez,  Physical Review A 98, 032333 (2018), 1808.04488.
  * (18) NuQS, A. Alexandru et al.,  Phys. Rev. D 100, 114501 (2019), 1906.11213.
  * (19) G. Jay, F. Debbasch, and J. B. Wang,  Physical Review A 99, 032113 (2019), 1803.01304.
  * (20) A. Y. Wei, P. Naik, A. W. Harrow, and J. Thaler,  Phys. Rev. D 101, 094015 (2020), 1908.08949.
  * (21) NuQS, H. Lamm, S. Lawrence, and Y. Yamauchi,  Phys. Rev. D 100, 034518 (2019), 1903.08807.
  * (22) A. Mott, J. Job, J.-R. Vlimant, D. Lidar, and M. Spiropulu,  Nature 550, 375 (2017).
  * (23) C. W. Bauer, W. A. De Jong, B. Nachman, and D. Provasoli,  (2019), 1904.03196.
  * (24) NuQS, A. Alexandru, P. F. Bedaque, H. Lamm, and S. Lawrence,  Phys. Rev. Lett. 123, 090501 (2019), 1903.06577.
  * (25) A. Blance and M. Spannowsky,  JHEP 21, 170 (2020), 2103.03897.
  * (26) NuQS, H. Lamm, S. Lawrence, and Y. Yamauchi,  Phys. Rev. Res. 2, 013272 (2020), 1908.10439.
  * (27) S. Abel, N. Chancellor, and M. Spannowsky,  (2020), 2003.07374.
  * (28) S. Abel and M. Spannowsky,  PRX Quantum 2, 010349 (2021), 2006.06003.
  * (29) K. Bepari, S. Malik, M. Spannowsky, and S. Williams,  Phys. Rev. D 103, 076020 (2021), 2010.00046.
  * (30) S. Williams, S. Malik, M. Spannowsky, and K. Bepari,  (2021), 2109.13975.
  * (31) G. D. Molfetta and P. Arrighi,  Quantum Inf Process 19, 47 (2020), 1906.04483.
  * (32) J. Y. Araz and M. Spannowsky,  (2022), 2202.10471.
  * (33) K. T. Matchev, P. Shyamsundar, and J. Smolinsky,  (2020), 2003.02181.
  * (34) W. A. De Jong et al.,  Phys. Rev. D 104, 051501 (2021), 2010.03571.
  * (35) V. S. Ngairangbam, M. Spannowsky, and M. Takeuchi,  Phys. Rev. D 105, 095004 (2022), 2112.04958.
  * (36) C. W. Bauer et al.,  (2022), 2204.03381.
  * (37) G. Agliardi, M. Grossi, M. Pellen, and E. Prati,  Phys. Lett. B 832, 137228 (2022), 2201.01547.
  * (38) A. Buckley et al.,  Phys. Rept. 504, 145 (2011), 1101.2599.
  * (39) P. Deliyannis et al.,  (2022), 2203.10018.
  * (40) R. D. Field and R. P. Feynman,  Nucl. Phys. B 136, 1 (1978).
  * (41) J. C. Collins, D. E. Soper, and G. F. Sterman,  Nucl. Phys. B 250, 199 (1985).
  * (42) J. C. Collins, D. E. Soper, and G. F. Sterman,  Adv. Ser. Direct. High Energy Phys. 5, 1 (1989), hep-ph/0409313.
  * (43) Y. Aharonov, L. Davidovich, and N. Zagury,  Phys. Rev. A 48, 1687 (1993).
  * (44) D. Aharonov, A. Ambainis, J. Kempe, and U. Vazirani,  STOC ’01: Proceedings of the thirty-third annual ACM symposium on Theory of computing, pg. 50 (2001).
  * (45) J. Kempe,  Contemporary Physics 44, 307 (2003)
  * (46) T. Sjostrand, S. Mrenna, and P. Z. Skands,  JHEP 05, 026 (2006), hep-ph/0603175.
  * (47) C. Bierlich et al.,  (2022), 2203.11601.
  * (48) Y. I. Azimov, Y. L. Dokshitzer, V. A. Khoze, and S. I. Troian,  Phys. Lett. B 165, 147 (1985).
  * (49) G. Gustafson and U. Pettersson,  Nucl. Phys. B 306, 746 (1988).
  * (50) B. Andersson, G. Gustafson, and J. Samuelsson,  Nucl. Phys. B 463, 217 (1996).
  * (51) A. Gehrmann-De Ridder, M. Ritzmann, and P. Z. Skands,  Phys. Rev. D 85, 014013 (2012), 1108.6172.
  * (52) M. McGettrick,  Quant. Inf. and Comput. , 509 (2010).
  * (53) A. Shakeel, D. Meyer, and P. Love,  J. Math. Phys. 55, 122204 (2014).
  * (54) E. Camilleri, P. Rohde, and J. Twamley,  Sci Rep 4, 4791 (2014).
  * (55) T. A. Brun, H. A. Carteret, and A. Ambainis,  Phys. Rev. A 67, 052317 (2003).
  * (56) P. P. Rohde, G. K. Brennen, and A. Gilchrist,  Phys. Rev. A 87, 052302 (2013).
  * (57) D. Li, M. Mc Gettrick, F. Gao, J. Xu, and Q.-Y. Wen,  Phys. Rev. A 93, 042323 (2016).
  * (58) M. Roget, B. Herzog, and G. Di Molfetta,  Scientific Reports 10, 2045 (2020).
  * (59) R. D. Somma, S. Boixo, H. Barnum, and E. Knill,  Phys. Rev. Lett. 101, 130504 (2008).
  * (60) P. Wocjan and A. Abeyesinghe,  Phys. Rev. A 78, 042336 (2008).
  * (61) N. Shenvi, J. Kempe, and K. B. Whaley,  Phys. Rev. A 67, 052307 (2003).
  * (62) M. Szegedy,  Quantum speed-up of markov chain based algorithms,  in 45th Annual IEEE Symposium on Foundations of Computer Science, pp. 32–41, 2004.
  * (63) A. Montanaro,  Proc. R. Soc. 471 (2015).
  * (64) J. Lemieux, B. Heim, D. Poulin, K. Svore, and M. Troyer,  Quantum 4, 287 (2020).
  * (65) D. Levin, Y. Peres, and E. Wilmer,  Markov Chains and Mixing Times (American Mathematical Soc., 2008).
  * (66) P. C. Richter,  Phys. Rev. A 76, 042306 (2007).
  * (67) D. Orsucci, H. J. Briegel, and V. Dunjko,  Quantum 2, 105 (2018).
  * (68) Y. Atia and S. Chakraborty,  Phys. Rev. A 104, 032215 (2021).
  * (69) ALEPH, A. Heister et al.,  Eur. Phys. J. C 35, 457 (2004).
  * (70) DELPHI, P. Abreu et al.,  Z. Phys. C 73, 11 (1996).
  * (71) JADE, OPAL, P. Pfeifenschneider et al.,  Eur. Phys. J. C 17, 19 (2000), hep-ex/0001055.
  * (72) J. Alwall et al.,  Comput. Phys. Commun. 176, 300 (2007), hep-ph/0609017.
  * (73) J. R. Andersen et al.,  (2014), 1405.1067.
  * (74) B. Andersson, G. Gustafson, G. Ingelman, and T. Sjostrand,  Phys. Rept. 97, 31 (1983).
  * (75) A. Buckley et al.,  Comput. Phys. Commun. 184, 2803 (2013), 1003.0694.
  * (76) C. Bierlich et al.,  SciPost Phys. 8, 026 (2020), 1912.05451.
  * (77) B. Andersson, P. Dahlkvist, and G. Gustafson,  Phys. Lett. B 214, 604 (1988).
  * (78) N. Klco et al.,  Phys. Rev. A 98, 032331 (2018), 1803.03326.
  * (79) G. Gustafson,  Nucl. Phys. B 392, 251 (1993).
  * (80) S. Frixione and B. R. Webber,  JHEP 06, 029 (2002), hep-ph/0204244.
  * (81) P. Nason,  JHEP 11, 040 (2004), hep-ph/0409146.
  * (82) S. Frixione, P. Nason, and C. Oleari,  JHEP 11, 070 (2007), 0709.2092.
  * (83) S. Catani, F. Krauss, R. Kuhn, and B. R. Webber,  JHEP 11, 063 (2001), hep-ph/0109231.
  * (84) L. Lönnblad,  JHEP 05, 046 (2002), hep-ph/0112284.
  * (85) M. L. Mangano, M. Moretti, and R. Pittau,  Nucl. Phys. B632, 343 (2002), hep-ph/0108069.
  * (86) S. Mrenna and P. Richardson,  JHEP 05, 040 (2004), hep-ph/0312274.
  * (87) J. Alwall et al.,  Eur. Phys. J. C 53, 473 (2008), 0706.2569.
  * (88) A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf,  Rev. Mod. Phys. 82, 1155 (2010).
  * (89) N. A. of Engineering 2019,  Frontiers of Engineering: Reports on Leading-Edge Engineering from the 2018 Symposium (The National Academies Press, 2018).
  * (90) D. Lidar and T. Brun,  Quantum Error Correction (Cambridge University Press, 2013).
  * (91) S. J. Devitt, W. J. Munro, and K. Nemoto,  Reports on Progress in Physics 76, 076001 (2013).
  * (92) Calibrating qubits using qiskit pulse, https://learn.qiskit.org/course/quantum-hardware-pulses/calibrating-qubits-using-qiskit-pulse
  * (93) IBM Quantum, 2021.
  * (94) P. Skands, S. Carrazza and J. Rojo,  Eur. Phys. J. C 74 (2014) no.8, 3024
  * (95) J. Isaacson and S. Prestel,  Phys. Rev. D 99 (2019) no.1, 014021