# Neural Network Decoders for Large-Distance 2D Toric Codes

Xiaotong Ni  QuTech, Delft University of Technology, P.O.Box 5046, 2600 GA
Delft, The Netherlands<EMAIL_ADDRESS>
(2020 April 8)

We still do not have perfect decoders for topological codes that can satisfy
all needs of different experimental setups. Recently, a few neural network
based decoders have been studied, with the motivation that they can adapt to a
wide range of noise models, and can easily run on dedicated chips without a
full-fledged computer. The later feature might lead to fast speed and the
ability to operate at low temperatures. However, a question which has not been
addressed in previous works is whether neural network decoders can handle 2D
topological codes with large distances. In this work, we provide a positive
answer for the toric code [1]. The structure of our neural network decoder is
inspired by the renormalization group decoder [2, 3]. With a fairly strict
policy on training time, when the bit-flip error rate is lower than $9\%$ and
syndrome extraction is perfect, the neural network decoder performs better
when code distance increases. With a less strict policy, we find it is not
hard for the neural decoder to achieve a performance close to the minimum-
weight perfect matching algorithm. The numerical simulation is done up to code
distance $d=64$. Last but not least, we describe and analyze a few failed
approaches. They guide us to the final design of our neural decoder, but also
serve as a caution when we gauge the versatility of stock deep neural
networks. The source code of our neural decoder can be found at [4].

##  1 Introduction

Before we can make the components of quantum computers as reliable as those of
classical computers, we will need quantum error correction so that we can
scale up the computation. The surface code and other topological codes are
popular choices for several qubit architectures because of their high
thresholds and low requirement on connectivity between qubits. However,
several good performing decoders have trouble to do real-time decoding for
qubits with fast error-correction cycles, such as superconducting qubits.
Moreover, as we are getting closer to the point where small size surface code
can be implemented in the lab, it is desirable that the decoders can adapt to
the noise models from the experimental setups. These considerations motivate
the study of decoders based on neural networks, which we will refer to as
neural decoders, for surface code and other topological codes [5, 6, 7, 8, 9,
10, 11, 12]. One question that has not been addressed so far is whether neural
networks can also be used for decoding 2D topological codes on a large lattice
with good performance. In this work, we will focus on answering this question
for the toric code. While it is the simplest topological code, it shares many
common features with others, which makes it a good test platform.

To design a neural decoder for large toric codes, a natural first step is to
use convolutional neural networks (CNNs) [13, 14], as the toric code and CNNs
are both translational-invariant on a 2D-lattice. Compared to normal neural
networks, the number of parameters in CNNs only scale with the depth of
networks. This gives an intuition that the training of the CNNs should remain
feasible for the lattice size of concern in the near future. We want the
decoder to be able to adapt to experimental noise, which we should assume to
be constantly changing, and thus the data for calibration is limited. The
structure of CNNs allows us to have a great control of how many parameters to
be re-trained during calibration so that we can avoid over-fitting (see
Appendix E for an example).

Interestingly, the renormalization group (RG) decoder [2, 3] for toric code
already has a structure very similar to the CNNs used in image classification.
Both of them try to keep the information needed for the output intact while
reducing the size of the lattice, by alternating between local computation and
coarse-grain steps. This similarity means that we should aim to achieve better
or similar performance with the neural decoder compared to the RG one. And in
case of bad performance, we can “teach” the neural decoder to use a similar
strategy as the RG decoder. This is indeed how we get good performance in the
end. Conceptually, this is similar to imitation learning (see [15] for an
overview). Even though we initialize the neural decoder by mimicking the RG
one, it can have the following advantages:

  * •

It can achieve a better performance than the RG decoder, as the latter one
contains some heuristic steps. On the other hand, the neural decoder can be
optimized to be a local minimum with respect to the parameters of the neural
network (strictly speaking, at least the gradient is very small). The idea of
improving belief propagation with neural networks is also used for decoding
classical linear codes [16].

  * •

It offers an additional way to adapt to experimental noise models, which is
simply training on experimental data.

It is tricky to evaluate the performance of neural decoders. As it stands, we
need to train the neural nets for different lattice sizes separately, and the
training process is not deterministic. Thus, we cannot define a threshold for
the decoder. This is fine if our main goal is to have an adaptable decoder for
near-future quantum devices. However, in order to know how optimal the neural
decoders are, we still make a “threshold plot” under a well-studied noise
model in Figure 4. Roughly speaking, the threshold benchmark is a good
indicator of how good the decoder can process syndrome information. At the
same time, we compare our neural decoder to the minimum-weight perfect
matching algorithm in Figure 5, and show in Appendix E that our neural decoder
can improve itself when trained on different error model. We hope these pieces
of information together can give us a first impression of neural decoders on
toric code.

The focus of this paper is not on how to obtain an optimal neural decoder.
Indeed, a lot of hyper-parameter optimizations can be done to further improve
the performance or reduce the amount of data needed for training. Instead, we
describe the key ideas that allow us to reliably obtain decent neural decoders
for the toric code. The knowledge we gained can help us design neural decoders
for other large codes.

##  2 Introduction to Toric code and the Renormalization Group Decoder

###  2.1 Toric Code

Figure 1: An illustration of the lattice and a $B_{p}$ stabilizer check. The
$2\times 2$ unit cells are marked by the red color. We also give an example of
a coarse-grained edge, which locates at the top-right and contains two edges
of the original lattice.

First, we give a brief introduction of toric code. Consider a $L\times L$
square lattice with periodic boundary condition, where a qubit lives on each
edge. The stabilizer group of toric code is generated by two types of operator
$A_{s}$ and $B_{p}$

$A_{s}=\bigotimes_{q\in n(s)}X_{q},\quad B_{p}=\bigotimes_{q\in n(p)}Z_{q},$
(1)

where $s$ and $p$ is any site and plaquette respectively, and $n(\cdot)$
consists of the 4 qubits neighboring $s$ or $p$. The logical-$Z$ operators
have the form

$\bar{Z}_{i}=\bigotimes_{q\in l_{i}}Z_{q},$ (2)

where $l_{1,2}$ are two shortest inequivalent non-contractible loops. The
toric code has a distance $d=L$.

In this paper, we will focus on the bit-flip noise model, i.e. only $X$ errors
can happen. We will also assume perfect measurements. Under this restriction,
the quantum states will stay in the $+1$ eigenspace of $A_{s}$. Therefore, we
only need to consider the expectation values of $B_{p}$ and $\bar{Z}_{i}$. For
simplicity, let us suppose in the beginning $\langle\bar{Z}_{i}\rangle=+1$.
And then a set of $X$ errors happened, which leads to the syndrome
$s=\\{\langle B_{p}\rangle\\}$. The goal of a decoder is to apply $X$ to the
qubits, such that $\langle B_{p}\rangle$ and $\langle\bar{Z}_{i}\rangle$
return to $+1$. Without going to detail, we claim it is enough to know the
parity of the number of $X$ errors that happened on the loops $l_{1,2}$. These
two parities will be the final training target for our neural decoder. We will
refer to the two parities as logical corrections.

###  2.2 Renormalization Group Decoder

Let us first set up some notation. We will use $e$ to denote an edge of the
original lattice or a coarse-grained edge. When we say $e$ is a coarse-grained
edge, we mean $e$ is an edge of a unit cell which consists of two (or more)
edges of the original lattice. We use $x(e)=1$ to denote an $X$ error happened
on edge $e$, and otherwise $x(e)=0$. When $e$ is a coarse-grained edge
consists of edges $\\{e_{i}\\}$, we set $x(e)$ to be

$x(e)=\sum_{i}x(e_{i})\mod 2.$ (3)

Lastly, the conditional marginal probability distribution
$p_{e}\left(x\left(e\right)|\text{syndrome}\right)$ of error on a coarse-
grained edge $e$ is denoted by $p_{e}$. Theoretically, $p_{e}$ can be computed
by enumerating all error configurations that have the given syndrome. However,
with renormalization group decoders and our neural decoders, we will only be
able to compute approximate distributions $p^{\prime}_{e}\approx p_{e}$.
Therefore, with a slight abuse of notation, we will use $p^{\prime}_{e}$ and
$p_{e}$ interchangeably. We will also use $p_{e}$ to denote the physical error
rate of an original edge $e$.

One renormalization stage consists of the following (see Figure 2 for an
example):

  1. 1.

Divide the lattice into $m\times m^{\prime}$ unit cells, where in this work
$m=m^{\prime}=2$.

  2. 2.

The ideal outputs of the renormalization step are $\\{p_{e}\\}$ for each
coarse-grained edge $e$ that is a border of a unit cell. However, we can only
compute $\\{p_{e}\\}$ approximately by using belief propagation, which is a
heuristic procedure for computing marginal probabilities (see Appendix A).
These approximate $\\{p^{\prime}_{e}\\}$ are treated as the error rate of the
coarse-grained edge $e$ for the next renormalization stage.

Figure 2: An illustration of the first renormalization stage. The numbers in
the figure are made up by the author and are likely not accurate. The left
side is the physical lattice, where the qubits have the same error rate. The
dots represent the locations of violated parity checks. The right side is the
coarse-grained lattice. The dots indicate that the corresponding coarse-
grained plaquettes contain an odd number of violated parity checks. The
numbers on coarse-grained edges are $\\{p^{\prime}_{e}(1)\\}$, which are
computed by belief propagation. With the help of $\\{p^{\prime}_{e}(1)\\}$, we
can more reliably figure out how the coarse-grained syndrome should be
matched.

At the end of the renormalization process, we obtain $p_{e}$ for $e$ being
either of the two non-contractible loops. For simplicity, we assume the two
non-contractible loops are $l_{1,2}$. Thus, we get an approximation of the
marginal probability distribution for logical correction.

##  3 Design and Training of the Neural nets

At a first glance, to build a neural decoder, we can simply train a
convolutional neural net with input-output pairs (syndrome, logical
correction). However, in practice, this does not allow us to get a good enough
performance. A detailed description of some simpler approaches and discussion
will be presented in Appendix C. Those failures eventually motivate us to
design and train the neural decoder in the following way.

###  3.1 Design of the network

The network follows the same structure as the renormalization decoder. Most of
the network is repetitively applying the renormalization block, which is
depicted in Figure 3. The belief propagation (BP) network, as its name
suggested, is intended to approximate the BP process (see Appendix A for an
introduction). More concretely, the first step of the training process is to
train the BP network with the data generated by a handcrafted BP algorithm.
This means initially the inputs to the BP network are syndromes and error
rates on each edge, and the outputs are supposed to approximate the error
rates on coarse-grained edges. However, later in the training process (i.e.
global training mentioned in subsection 3.2), the BP network can deviate from
this initial behavior. The post-processing has two steps. The first step is to
remove the superficial complexity from the coarse-grained lattice. Whenever
for a coarse-grained edge $e$ has $p_{e}(1)>0.5$, we apply an $X$ on $e$ and
switch $p_{e}(1)\leftrightarrow p_{e}(0)$. If $e$ is on either of the two non-
contractible loops $l_{1,2}$, then the desired logical correction will be
updated accordingly. Although this step only changes the representation of the
data, and in principle, neural nets can learn to do the same thing, it is a
quite costly step for neural nets as it can call the parity function multiple
times. The second step is coarse-graining. We need to reduce the lattice size
by half, and for convenience, this is done by the first layer of every belief
propagation network. We also compute the parity of four $\langle B_{p}\rangle$
in each unit cell and feed these parities to the next BP network as the
effective syndrome of the coarse-grained lattice.

Figure 3: Structure of the entire network, and the training order labeled by
the blue circles. After loading the pre-trained belief propagation network
into RN block, the first step is to train the dense layers, and the second
step is to train all the layers together. We will call the second step global
training.

In more detail, the input to the BP network can be packed in a tensor $I$ with
shape $(l,l,3)$, where $l$ is the initial lattice size or the output size of
the precedent BP network. For example, we can set $I(i,j,0)$ to be $\langle
B_{p}\rangle$ on plaquette $(i,j)$, and $I(i,j,1),I(i,j,2)$ to be the error
rates corresponding to the top and left qubits of the plaquette. Each BP
network consists of 13 convolution and 3 batch normalization layers. The
definition of convolution layers can be found in Appendix B, and batch
normalization is introduced in [17]. The first layer reduces the lattice size
$L$ by half. The reasoning is that the belief propagation is done based on
$2\times 2$ unit cells. The remaining layers keep the lattice size unchanged.
Among them, only four involve communication between unit cells, i.e. the
kernels of these four convolution layers have size $3\times 3$. They spread
evenly in the 13-layer network. Other layers only have kernels of size
$1\times 1$, which can then be viewed as computation inside unit cells. The
rationale behind this is that the messages likely need to be processed before
the next round of communication. The batch normalization layers also spread
evenly, with the hope that they can make the training more stable.

After the renormalization process reduces the lattice to a size of $2\times
2$, we apply 4 dense layers (a.k.a fully-connected layers). Note that the
dense layers conveniently break the translational symmetry imposed by the
convolution layers. In the end, we have a neural network with input shape
$(L,L,3)$ and output shape $(2)$111For efficient training, an additional
dimension called batch size will be added.. The input shape is $(L,L,3)$
because this is the input shape of BP networks. For $L=64$, the total number
of trainable layers in the network is around 60, which is very large compared
to early deep neural networks [14]. However, most of the computation cost and
the trainable parameters are concentrated in the 16 convolution layers with
kernel size $3\times 3$. Combining this and the careful training strategy we
describe below, we find that the training can be done very efficiently.

###  3.2 Training

In general, training neural networks becomes harder when the number of layers
increases. This is often attributed to the instability of gradient
backpropagation. Considering we have a very deep neural network, we should
find a way to train parts of the network first. The training is divided into
two stages. First, we train the belief propagation network to indeed do belief
propagation (BP). This corresponding to the blue circle with 0 in Figure 3. To
do this, we implement a BP algorithm and use it to generate training data for
the network. More concretely, we first assign a random error rate $e^{-k}$ to
each edge, where $k\in[0.7,7]$ from a uniform distribution. The choice of the
distribution is quite arbitrary. Then we sample error on each edge according
to its error rate and compute the syndrome. After that, we feed both the error
rates and syndrome into our handcrafted BP algorithm, which will output an
estimation of the error rates $p_{e}$ corresponding to the coarse-grained
edges. We can subsequently train the BP network with the same input-output
relation. An important detail is that we transform the error rates $p_{e}(1)$
in both input and output to $r_{e}=\log\left(p_{e}(1)/p_{e}(0)\right)$. The
reason behind this is described in Appendix C.

Next, we load the pre-trained belief propagation network into the decoder
network described in the previous subsection. To ensure $r_{e}$ stay bounded,
we perform a rescale $r_{e}\rightarrow 7r_{e}/\max_{e}|r_{e}|$ before feed it
into next RN block (the choice of 7 here is arbitrary). We can then train the
dense layers and afterward the whole network with input-output pairs
(syndrome, logical correction). These two trainings correspond to the blue
circle 1 and 2 in Figure 3, respectively. The training data is measurable in
experiments in these two training.

We train the decoders for different lattice sizes $L$ separately. Although
this makes the concept of threshold pointless, it is still useful to estimate
the “threshold” so that we can have a rough comparison of the neural decoder
with the existing ones. For this, we train the decoder for different $L$ with
the same amount of stochastic gradient steps, which also implies the optimizer
sees the same amount of training data for each $L$. In addition, the training
for each $L$ is done under 1 hour (on the year 2016 personal computer with 1
GPU). We consider this to be a fairly strict policy. The result is plotted in
Figure 4. We can also forgo this strict policy and spend more time in training
the neural decoder for $d=64$ toric code, which gives rise to Figure 5. The
training time is still under 2 hours. More details about the design and
training can be found in Appendix D and the source code [4], and more
discussion about the numerical results can be found in the following section.

##  4 Numerical results

Figure 4: Logical accuracy versus physical error rate. The neural decoders are
trained at physical error rate $9\%$. For the three solid lines, the decoder
has been trained globally, while the dashed lines it has not. The colors of
the dashed lines indicate the code distance they are evaluated on. The
vertical grid indicates the physical error rates for which we evaluate the
logical accuracy.  Figure 5: In this figure, we compare the performance of our
neural decoders to the MWPM algorithm. The solid line for $d=16$ and the
dashed line for $d=64$ decoder are using the strict training policy, while
more training has been done on the $d=64$ decoder corresponding to the solid
line. The “star” points are the performance of the minimum-weight perfect
matching algorithm. The colors of stars indicate the code distance they are
evaluated on. The vertical grid indicates the physical error rates for which
we evaluate the logical accuracy for the lines. We see the performance of
neural decoders can be almost as good as MWPM for a decent range of physical
error rates.

For the strict training policy, we plot the logical accuracies versus the
physical error rates in Figure 4. Logical accuracy is simply $(1-\text{logical
error rate})$ and is averaged over the two logical qubits. For the solid
lines, the decoders have been trained globally, i.e. have done both steps 1
and 2 in Figure 3. For the dashed lines, the decoders only did step 1, i.e.
only the dense layers are trained. The colors of the dashed lines indicate the
code distance they are evaluated on. The vertical grid indicates the physical
error rates for which we evaluate the logical accuracy, where for each point
we sample $10^{4}$ (syndrome, logical correction) pairs. We can see that the
solid lines cross around $p_{\text{physical}}=0.095$; therefore we might say
our neural decoder has an effective threshold around $9.5\%$. It can be seen
that the global training is crucial for getting a decent performance because
without it the effective threshold will be below $8\%$.

We can also spend more time to train the $d=64$ decoder, and then compare the
performance of the neural decoders to the minimum-weight perfect matching
algorithm (MWPM) in Figure 5. The “star” points are the logical accuracies of
MWPM, where each one is evaluated by $3000$ trials. The $d=16$ decoder
corresponding to the solid line and the $d=64$ decoder corresponding to the
dashed line are from the strict training policy. The $d=64$ decoder
corresponding to the solid line is obtained by doing more training while
having the same network architecture. We see that without the strict training
policy, the performance of the neural decoder is almost identical to MWPM for
a decent range of physical error rates. We can also compare to the
renormalization group (RG) decoder in [2], where the authors have shown a
threshold of $8.2\%$ when using $2\times 1$ unit cell, and claim a threshold
around $9.0\%$ if using $2\times 2$ unit cell. With the strict training
policy, our neural decoder is slightly better or at least comparable to the RG
decoder, while without the policy our neural decoder is clearly better for
$d\leq 64$.

##  5 Discussion

One obvious question is whether we can get a good neural decoder for surface
code or other topological codes on large lattices. In the case of surface
code, the major difference compared to the toric code is the existence of
boundaries. This means we have to inject some non-translational invariant
components into the network. For example, we can have a constant tensor $B$
with shape $(L,L,2)$ marks the boundary, e.g. $B(x,y,i)=1$ if $(x,y)$ is at
the smooth boundary and $i=0$, or if $(x,y)$ is at the rough boundary and
$i=1$; otherwise $B(x,y,i)=0$. We then stack $B$ with the old input tensor
before feed into the neural decoder. More generally, if a renormalization
group decoder exists for a topological code, we anticipate that a neural
decoder can be trained to have similar or better performance. For example,
neural decoders for surface code with measurement errors, for topological
codes with abelian anyons can be trained following the same procedure
described in this paper.

Another question we want to discuss is the viability of our neural decoder at
low physical error rates. On the one hand, we can train our neural decoders to
approximate the RG decoder, and therefore they can have similar performance at
low error rates. On the other hand, it will be much harder to improve neural
decoders just by training on experimental data, because it will take a long
time to encounter syndromes that are decoded incorrectly. Therefore, we should
expect neural decoders to gradually lose the ability to adapt to experimental
noise models as the physical error rates decrease.

We want to discuss a bit more about running neural networks on specialized
chips. It is straightforward to run our neural decoder on GPU or TPU [18] as
they are supported by Tensorflow [19], the neural network library used in this
work. There is software (e.g. OpenVINO) to compile common neural networks to
run on commercially available field-programmable gate arrays (FPGAs), but we
do not know how easy it is for our neural decoder222The only uncommon
component of our neural decoder is element-wise parity function.. Apart from
power efficiency, there is a study about operating FPGAs at 4K temperature
[20]. Overall, there is a possibility to run neural decoders at low
temperature. Note that for running on FPGAs or benchmarking the speed, it is
likely a good idea to first compress the neural networks, see [21].

##  6 Acknowledgement

The author wants to thank Ben Criger, Barbara Terhal, Thomas O’Brien for
useful discussion. The implementation of minimum-weight perfect matching
algorithm, including the one used in Appendix E, is provided by Christophe
Vuillot, which uses the backend from either Blossom V [22] or NetworkX [23].
The author acknowledge support through the ERC Consolidator Grant No. 682726.

## References

  * Kitaev [2003] A.Yu. Kitaev.  Fault-tolerant quantum computation by anyons.  _Annals of Physics_ , 303(1):2–30, jan 2003.  doi: 10.1016/s0003-4916(02)00018-0.
  * Duclos-Cianci and Poulin [2010] Guillaume Duclos-Cianci and David Poulin.  Fast decoders for topological quantum codes.  _Physical review letters_ , 104(5):050504, 2010\.  doi: 10.1103/PhysRevLett.104.050504.
  * Duclos-Cianci and Poulin [2014] Guillaume Duclos-Cianci and David Poulin.  Fault-tolerant renormalization group decoder for abelian topological codes.  _Quantum Information & Computation_, 14(9-10):721–740, 2014.
  * [4] https://github.com/XiaotongNi/toric-code-neural-decoder.
  * Baireuther et al. [2018] Paul Baireuther, Thomas E. O'Brien, Brian Tarasinski, and Carlo W. J. Beenakker.  Machine-learning-assisted correction of correlated qubit errors in a topological code.  _Quantum_ , 2:48, jan 2018.  doi: 10.22331/q-2018-01-29-48.
  * Varsamopoulos et al. [2017] Savvas Varsamopoulos, Ben Criger, and Koen Bertels.  Decoding small surface codes with feedforward neural networks.  _Quantum Science and Technology_ , 3(1):015004, nov 2017.  doi: 10.1088/2058-9565/aa955a.
  * Torlai and Melko [2017] Giacomo Torlai and Roger G Melko.  A neural decoder for topological codes.  _Physical Review Letters_ , 119(3):030501, 2017\.  doi: 10.1103/PhysRevLett.119.030501.
  * Breuckmann and Ni [2018] Nikolas P. Breuckmann and Xiaotong Ni.  Scalable neural network decoders for higher dimensional quantum codes.  _Quantum_ , 2:68, may 2018.  doi: 10.22331/q-2018-05-24-68.
  * Baireuther et al. [2019] Paul Baireuther, MD Caio, B Criger, Carlo WJ Beenakker, and Thomas E O’Brien.  Neural network decoder for topological color codes with circuit level noise.  _New Journal of Physics_ , 21(1):013003, 2019\.  doi: 10.1088/1367-2630/aaf29e.
  * Krastanov and Jiang [2017] Stefan Krastanov and Liang Jiang.  Deep neural network probabilistic decoder for stabilizer codes.  _Scientific Reports_ , 7(1), sep 2017.  doi: 10.1038/s41598-017-11266-1.
  * Maskara et al. [2019] Nishad Maskara, Aleksander Kubica, and Tomas Jochym-O’Connor.  Advantages of versatile neural-network decoding for topological codes.  _Physical Review A_ , 99(5):052351, 2019.  doi: 10.1103/PhysRevA.99.052351.
  * Chamberland and Ronagh [2018] Christopher Chamberland and Pooya Ronagh.  Deep neural decoders for near term fault-tolerant experiments.  _Quantum Science and Technology_ , 3(4):044002, jul 2018.  doi: 10.1088/2058-9565/aad1f7.
  * LeCun et al. [1989] Yann LeCun et al.  Generalization and network design strategies.  _Connectionism in perspective_ , pages 143–155, 1989.
  * Krizhevsky et al. [2012] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  In F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems 25_ , pages 1097–1105. Curran Associates, Inc., 2012.
  * Attia and Dayan [2018] Alexandre Attia and Sharone Dayan.  Global overview of imitation learning.  2018\.
  * Nachmani et al. [2016] Eliya Nachmani, Yair Be'ery, and David Burshtein.  Learning to decode linear codes using deep learning.  In _2016 54th Annual Allerton Conference on Communication, Control, and Computing (Allerton)_. IEEE, sep 2016.  doi: 10.1109/allerton.2016.7852251.
  * Ioffe and Szegedy [2015] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  2015\.
  * Jouppi et al. [2017] Norman P. Jouppi, Al Borchers, Rick Boyle, Pierre luc Cantin, Clifford Chao, Chris Clark, Jeremy Coriell, Mike Daley, Matt Dau, Jeffrey Dean, Ben Gelb, Cliff Young, Tara Vazir Ghaemmaghami, Rajendra Gottipati, William Gulland, Robert Hagmann, C. Richard Ho, Doug Hogberg, John Hu, Robert Hundt, Dan Hurt, Julian Ibarz, Nishant Patil, Aaron Jaffey, Alek Jaworski, Alexander Kaplan, Harshit Khaitan, Daniel Killebrew, Andy Koch, Naveen Kumar, Steve Lacy, James Laudon, James Law, David Patterson, Diemthu Le, Chris Leary, Zhuyuan Liu, Kyle Lucke, Alan Lundin, Gordon MacKean, Adriana Maggiore, Maire Mahony, Kieran Miller, Rahul Nagarajan, Gaurav Agrawal, Ravi Narayanaswami, Ray Ni, Kathy Nix, Thomas Norrie, Mark Omernick, Narayana Penukonda, Andy Phelps, Jonathan Ross, Matt Ross, Amir Salek, Raminder Bajwa, Emad Samadiani, Chris Severn, Gregory Sizikov, Matthew Snelham, Jed Souter, Dan Steinberg, Andy Swing, Mercedes Tan, Gregory Thorson, Bo Tian, Sarah Bates, Horia Toma, Erick Tuttle, Vijay Vasudevan, Richard Walter, Walter Wang, Eric Wilcox, Doe Hyun Yoon, Suresh Bhatia, and Nan Boden.  In-datacenter performance analysis of a tensor processing unit.  In _Proceedings of the 44th Annual International Symposium on Computer Architecture - ISCA '17_. ACM Press, 2017.  doi: 10.1145/3079856.3080246.
  * Abadi et al. [2015] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  URL https://www.tensorflow.org/.  Software available from tensorflow.org.
  * Lamb et al. [2016] I. D. Conway Lamb, J. I. Colless, J. M. Hornibrook, S. J. Pauka, S. J. Waddy, M. K. Frechtling, and D. J. Reilly.  An FPGA-based instrumentation platform for use at deep cryogenic temperatures.  _Review of Scientific Instruments_ , 87(1):014701, jan 2016.  doi: 10.1063/1.4939094.
  * Cheng et al. [2018] Yu Cheng, Duo Wang, Pan Zhou, and Tao Zhang.  Model compression and acceleration for deep neural networks: The principles, progress, and challenges.  _IEEE Signal Processing Magazine_ , 35(1):126–136, jan 2018.  doi: 10.1109/msp.2017.2765695.
  * Kolmogorov [2009] Vladimir Kolmogorov.  Blossom v: a new implementation of a minimum cost perfect matching algorithm.  _Mathematical Programming Computation_ , 1(1):43–67, apr 2009.  doi: 10.1007/s12532-009-0002-8.
  * Hagberg et al. [2004–] Aric Hagberg, Dan Schult, Pieter Swart, et al.  Networkx, 2004–.  URL https://networkx.github.io/.
  * Goh [2017] Gabriel Goh.  Why momentum really works.  _Distill_ , 2017.  doi: 10.23915/distill.00006.
  * Bengio and LeCun [2007] Yoshua Bengio and Yann LeCun.  Scaling learning algorithms towards ai.  _Large-scale kernel machines_ , 34(5):1–41, 2007\.
  * Erhan et al. [2010] Dumitru Erhan, Yoshua Bengio, Aaron Courville, Pierre-Antoine Manzagol, Pascal Vincent, and Samy Bengio.  Why does unsupervised pre-training help deep learning?  _Journal of Machine Learning Research_ , 11(Feb):625–660, 2010.
  * Maas et al. [2013] Andrew L Maas, Awni Y Hannun, and Andrew Y Ng.  Rectifier nonlinearities improve neural network acoustic models.  In _Proceedings of the 30th International Conference on Machine Learning_ , volume 28 of _JMLR Workshop and Conference Proceedings_ , Atlanta, Georgia, USA, 2013.
  * Kingma and Ba [2015] Diederik Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _3rd International Conference for Learning Representations, San Diego_ , 2015.  URL https://arxiv.org/abs/1412.6980.

##  Appendix A Implementation of Belief Propagation Algorithm

Belief propagation is a heuristic procedure for computing marginal
probabilities of graphical models. We choose to use a slightly different
belief propagation implementation compared to [2], as ours seems to be more
natural for the bit-flip noise model. We divide the lattice into $2\times 2$
unit cells. Let $G$ be a bipartite graph, where one part corresponds to unit
cells, and another part to coarse-grained edges. Two vertices in $G$ is
connected when the coarse-grained edge is adjacent to the unit cell. This
later decides how the messages flow in the graph. However, we assign two
variables $\\{x(e_{i}),x(e_{j})\\}$ to each vertex corresponding a coarse-
grained edge which $e_{i}$ and $e_{j}$ form. In this section, the symbol $e$
or $e_{i}$ will denote original edges of the lattice (i.e. not coarse-
grained). We define $E$ to be the set of all edges $e$, $E_{cg}\subset E$ to
be the set of $e$ which are components of coarse-grained edges (i.e. red edges
in Figure 1), and $\bar{E}_{cg}=E\setminus E_{cg}$. Given a syndrome $S$, the
unnormalized probability of an error configuration
$\vec{x}\equiv\\{x(e)\\}_{e\in E}$ can be written as

$p(\vec{x})=g(S,\vec{x})\prod_{e\in E}p_{e}(x(e)),$ (4)

where $g(S,\vec{x})=1$ if $\vec{x}$ has syndrome $S$, and otherwise
$g(S,\vec{x})=0$. It is obvious $g(S,\vec{x})$ can be factorized to local
terms. Thus, the marginal distribution for $e\in E_{cg}$ can then be
factorized according to $G$ as

$\sum_{\\{x(e),e\in\bar{E}_{cg}\\}}p(\vec{x})=\prod_{c}f_{c}(\\{x(e)\\}_{c}),$
(5)

where the product is taken over all unit cells $c$, and $\\{x(e)\\}_{c}$ are
the set of $x(e)$ such that $e\in E_{cg}$ is adjacent to $c$. We can then
apply the standard belief propagation to the graph $G$. Without further
explanation, we choose to use the following rule. A unit cell $c_{k}$ sends to
each of its adjacent cell $c_{n}$ a message containing 4 real numbers

$\\{m_{c_{k},c_{n}}(x(e_{i}),x(e_{j}))\\}\quad\text{for}\
x(e_{i}),x(e_{j})=0,1,$ (6)

where $e_{i}$ and $e_{j}$ form the coarse-grained edge between $c_{k}$ and
$c_{n}$. When we already fix an error configuration $\vec{x}$, we may use the
simplified notations

$p_{e}\equiv p_{e}(x(e)),\qquad m_{c_{k},c_{n}}\equiv
m_{c_{k},c_{n}}(x(e_{i}),x(e_{j})).$ (7)

To compute an out-going message from cell $c$, we take messages from the other
three directions, and consider them to be the probability of error
configuration on respective edges. We then sum over all error configurations
in the cell, which give the correct syndrome of the 4 plaquette stabilizer
checks.

Figure 6: An illustration of message passing for a unit cell.

More concretely, we define $\vec{x}_{c}$ to be $\vec{x}$ restricted to edges
in $c$ (i.e. all edges in Figure 6), and $g^{\prime}(S,\vec{x}_{c})$ checks
whether $\vec{x}_{c}$ is compatible with $S$ similar to $g$. We assume we want
to send the messages from $c$ to $c_{n}$, while the incoming messages are from
cells $I=\\{c_{k}\\}$. Then we have

$m_{c,c_{n}}=\sum_{\vec{x}_{c}}g^{\prime}(S,\vec{x}_{c})\prod_{c_{k}\in
I}m_{c_{k},c}\prod p_{e_{i}}.$ (8)

For the last term $\prod p_{e_{i}}$, the product is taken over the blue edges
in Figure 6, assuming $c_{n}$ is on the right of $c$.

In the end of the message passing, we can compute the marginal probability by

$P(x(e_{i}),x(e_{j}))=m_{c_{k},c_{n}}m_{c_{n},c_{k}}/\left(p_{e_{i}}p_{e_{j}}\right),$
(9)

where $e_{i}$ and $e_{j}$ are the edges between $c_{n}$ and $c_{k}$. From the
joint distribution $P(x(e_{i}),x(e_{j}))$ we can compute the distribution
$P\left(x(e_{i})+x(e_{j})\mod 2\right)$. It is not hard to see that the above
message passing rules will lead to the correct marginal probability when the
underlying graph is a tree (note this is not the case for the square lattices
we are considering). To generate training data for neural networks, we do 7
rounds of message passing defined above.

The key differences between our implementation and the one in [3] is

  * •

Ours utilizes all four stabilizer checks in each unit cell while in [3] only
three are used.

  * •

Each message contains 4 real numbers in our implementation while only 1 in
[3].

##  Appendix B Introduction to Neural Networks

A neural network, at the highest abstraction, can just be viewed as a black-
box function $f_{\text{nn}}(x,\vec{w})$ with many parameters $\vec{w}$ to be
tuned. We want $f$ to describe the input-output relation presented in a
dataset $D=\\{(x_{i},y_{i})\\}$. To do this333In this work, we choose to not
study the effect of overfitting, as we have the ability to generate infinite
data., we choose a (smooth) loss function $\mathcal{L}$, and then we do the
minimization

$\min_{\vec{w}}\sum_{i}\mathcal{L}(f(x_{i},\vec{w}),y_{i}).$ (10)

One important requirement is that $f$ is (almost-everywhere) differentiable
with respect to $\vec{w}$. This allows us to train the network with gradient
descent, for which a good introduction can be found in [24]. In general, we
can expect gradient descent will take us to a local minimum or some region
with very small gradients. This is the advantage of “end-to-end” training
compared to human-written heuristic algorithms, as the latter are unlikely to
be a local optimum (assuming we can add real number parameters to those
heuristic algorithms). A common loss function for classification problems is
cross-entropy. Assume we have a dataset $D=\\{(x_{i},y_{i})\\}$ where
$y_{i}\in\\{0,1\\}$, and the neural network output $y_{i}^{\prime}$ which
tries to approximate $\text{Prob}(y_{i}=1)$, the cross-entropy loss function
is then calculated as following:

$-\sum_{i}\left(y_{i}\log
y_{i}^{\prime}+\left(1-y_{i}\right)\log\left(1-y_{i}^{\prime}\right)\right).$
(11)

Note that when we use the notation $\text{Prob}(y_{i}=1)$, we implicitly
assume $D$ is obtained by sampling from an underlying probability
distribution.

More concretely, most neural networks consist of many layers. In this paper,
the two relevant types of layers are the dense and convolution layer. Dense
layers (a.k.a fully-connected layers) have the form $g(A\vec{x}+\vec{b})$,
where $g$ is some non-linear function applied entrywise, and the matrix $A$,
vector $\vec{b}$ are the trainable parameters. Assuming $A$ has a shape of
$n\times m$, we will say the output dimension of the dense layer is $n$. One
convolution layer, as the name suggests, contains a collection of discrete
convolutions. For this paper, the input to the layer resides on a
$2$-dimensional lattice of size $l^{2}$ with periodic boundary condition. On
each lattice site, there is a $d$-dimensional input vector
$x_{\vec{u}}\in\mathbb{R}^{d}$, where the subscript
$\vec{u}\in\mathbb{Z}_{l}^{2}$ (we use $\mathbb{Z}_{l}$ to denote integer in
range $[0,l-1]$). We define the kernel to be a tensor $K_{\vec{u},i}$, where
$\vec{u}\in\mathbb{Z}_{n}^{2}$ and $i\in\mathbb{Z}_{d}$. With a slight abuse
of notation, we will say such a kernel has size $n^{2}$. The convolution is
then

$y_{\vec{v}}=\sum_{\vec{u}\in\mathbb{Z}_{n}^{2}}\sum_{i\in\mathbb{Z}_{d}}x_{\vec{v}-\vec{u},i}K_{\vec{u},i},$
(12)

where $x_{\vec{v}-\vec{u},i}$ is the $i$th element of $x_{\vec{v}-\vec{u}}$,
and $\vec{v}-\vec{u}$ is calculated module $l$ because of the periodic
boundary condition. We will have a collection of kernels
$\\{K_{\vec{u},i,j}\\}_{j}$ for one convolution layer, which means we also
have a collection of outputs $\\{y_{\vec{v},j}\\}_{j}$. The cardinality of
$\\{K_{\vec{u},i,j}\\}_{j}$ is conventionally called the number of filters.
After Equation 12, we can also apply a non-linear function $g$ entrywise if
needed.

Before concluding this section, let us make one clarification. In this paper,
sometimes we only train part of the network, e.g. for blue circle 1 in Figure
3 we only train the dense layers. Assuming $\vec{w}_{1}$ are the parameters in
the part of the network we want to train and $\vec{w}_{2}$ are the rest, then
we are doing the optimization

$\min_{\vec{w}_{1}}\sum_{i}\mathcal{L}(f(x_{i},\vec{w}_{1},\vec{w}_{2}),y_{i})$
(13)

by using some gradient descent optimizer.

##  Appendix C Comparison to Simpler Approachers

In this section, we will show the performance of the neural net decoders when
trained with simpler approaches (more precisely, approaches with less human
involvement and prior knowledge of toric code decoding), and provide some
reasoning if possible. The neural nets will be the same as the ones we used in
the main text, except that they do not contain the “removing complexity” step
in the post-processing. We will see in general the performance gets much
worse, especially when the lattice size grows large. However, this does not
mean these simpler approaches will always fail. It just implies that a large
amount of training time / human involvement is needed, which could make them
impossible in practice.

The simplest approach is to train the whole network with input-output pairs
(syndrome, logical correction). During limited attempts, this approach does
not produce decoders much better than random guess for large toric codes. A
hand-waving explanation is the following. It is fair to assume a lot of parity
functions need to be evaluated during the decoding process. It is known that
the parity is not an easy function for neural nets to compute [25], and one
good way to approximate it is to increase the depth of the network. So let us
assume each renormalization stage needs 5 layers. This means to decode $L=32$
toric code, the network will have 25 layers, which exceeds the range where
neural nets can be reliably trained.

The problem of too many layers can be alleviated if we can pre-train the
earlier layers of the network. A similar strategy was used in training neural
nets for computer vision problems [26]. For the bit-flip noise model, we can
pre-train the earlier layers to mimic the renormalization group decoder.
Recall in section 2, we mentioned that for the renormalization group decoder,
the output corresponding to a coarse-grained edge $e$ is $p_{e}$. Since we
generate syndromes by first sampling $x(e)$ for all (not coarse-grained) edges
$e$, we also have the ability to generate pairs (syndrome, $\\{x(e)\\}$ for
coarse-grained edges $e$) for training. Although $\\{x(e)\\}$ in the above
pairs are binary numbers, with the cross-entropy as the cost function, in
theory the output will converge to $p_{e}(1)$. The pre-training is done one
renormalization block at a time. More concretely, with the network we are
using in the main text, we will train the output of the 12th layer with the
training target of first renormalization block, and the output of the 24th
layer with the second block, etc. We can try this method on $L=32$ toric code
and bit-flip error rate $0.08$. It does not work well, as we can see in Table
1.

Block | RN1 | RN2 | RN3 | RN4 | Dense
---|---|---|---|---|---
Cross-entropy | 0.16 | 0.22 | 0.28 | 0.4 | 0.5
Accuracy | 0.93 | 0.90 | 0.86 | 0.79 | 0.58
Table 1: Cross-entropy and accuracy after each renormalization block when we
do the training block by block. It is done with a bit-flip error rate $0.08$.
Although the accuracy after the first RN block is comparable to the input
error rate, it gradually decreases until barely above $0.5$

These numbers are not very accurate and coming from a single training
instance. However, the author has observed the same trend many times that the
loss and accuracy slowly degrade in the process of renormalization, even
though the error rate is way below the theoretical threshold. To diagnose the
reason, we first notice that the first RN block actually performs reasonably
well. This suggests that the optimizer is capable of training each RN block
alone. Assume this is indeed the case, the degrading of performance is likely
caused by the following two reasons.

First, later in the renormalization process, if we look at the coarse-grained
syndrome or $p_{e}$ alone, they behave more and more like white noise. While
it is possible for the neural nets to do the same post-processing described in
subsection 3.1, a few layers of the network will be occupied by this.
Therefore, the natural solution is to implement the post-processing ourselves.
By doing this, we suspect it is possible to reach a threshold of $8\%$, but
apparently $8\%$ is still not good enough.

The second reason is related to the convergence of $p_{e}$. When below the
threshold, we will encounter very often that $p_{e}$ is very close to $0$ or
$1$. For example, if $x(e)=1$ vs $x(e)=0$ corresponding to a weight $4$ vs
weight $1$ local configuration, then we will have $p_{e}\approx p_{0}^{3}$,
where $p_{0}$ is the initial error rate. This will become more prominent later
in the renormalization process, as $p_{e}$ from the previous renormalization
stage become the error rate in the next stage. It is important to know how
close to $0$ (or $1$) $p_{e}$ is on the logarithmic scale. Otherwise, in the
later stage, the information will not be accurate enough to deduce
configuration close to the minimum weight of errors. This poses the following
requirements:

  * •

When we pass $p_{e}$ to some intermediate layer of a neural net, it should be
able to distinguish between small $p_{e}$. However, recall that each layer
does the computation $f(Ax+b)$, where $f$ has a bounded derivative. Thus, to
distinguish a set of small $\\{p_{e}\\}$, we need $\|A\|\sim
1/\min\\{p_{e}\\}$. This will either not be achieved by training or cause
instability of the network. Another issue related to the minuscule nature of
$p_{e}$ is the cross-entropy loss function does not provide enough motivation
for $p_{e}$ to converge to the target value $q$ in log-scale. More accurately,
the derivative of the cross-entropy scales like $O(|p_{e}-q|)$ when
$p_{e}\approx q$, which will be too small before the convergence in log-scale.
A natural solution is we replace the appearance of $p_{e}$ with $\log
p_{e}-\log(1-p_{e})$.

  * •

Even with a good representation of $p_{e}$, we shall still be very cautious
about the training, as we are trying to estimate very small $p_{e}$ from
sampling. In the end, we decide to implement a belief propagation routine and
use the input-output pair from the routine to train the network. The advantage
is that belief propagation directly outputs the probability, which should be
reasonably accurate in the logarithmic scale. Therefore, we can get a much
stable training process.

##  Appendix D Technical Details

The objective of this section is to describe some technical details for people
who do not plan to read the source code.

For the majority of network, we use leaky ReLUs [27] as the activation
function, which has the form

$y=x\quad\text{if}\ x>0;\quad y=0.2x\quad\text{if}\ x\leq 0$ (14)

Apart from the last layer of each renormalization block and the last dense
layer, the number of filters in each convolution layer is 200, and the output
dimension of each dense layer is 50.

For training a belief propagation network, we generate a dataset of size
$80000$. The dataset consists of the input and output of the belief
propagation algorithm described in Appendix A when applied to $d=16$ toric
code. The optimizer we use is the ADAM [28]. The learning rate parameter of
the optimizer is set to $7\times 10^{-4}$ for training belief propagation
network, $10^{-3}$ for training dense layers, $7\times 10^{-5}$ for global
training of $L=16,32$ lattice, and $7\times 10^{-6}$ for $L=64$. The batch
size for training is 50. The training of the dense layers uses around 1000
batches. For the strong policy, the global training uses 3000 batches
regardless of the code distance. To see the potential of our neural decoder at
$d=64$, we also did a longer training and compared it to MWPM in Figure 5. In
total, it is trained using 18000 batches. The first 12000 batches are trained
on physical error rate $p=0.09$, and the last 6000 batches are on $p=0.095$.
The reason of switching to a higher error rate for late-stage training is that
the accuracy at $p=0.09$ is too close to 1 for effective training.

##  Appendix E Spatially Varying Error Rates

A natural use of the neural decoder is to train it with experimental data.
Naturally, the noise models in the experiments will not be translational
invariant. There are two simple ways to reconcile this fact with the
translational-invariance of convolutional neural nets:

  1. 1.

Allow the first few layers of the network to be non-translational invariant.

  2. 2.

Introduce site-dependent trainable variables to the networks.

In this section, we will consider the error model that has varying bit-flip
error rates across the lattice. For this, we can use the second approach,
where the site-dependent variables can in principle represents the varying
error rates. In fact, recall that the neural decoder has an input shape
$(L,L,3)$, which contains $2L^{2}$ numbers that are originally error rates
feeding into belief propagation. Thus we can simply feed the site-dependent
variables into the neural decoder and then train them.

However, there is still a complication regarding the starting point of this
training. A natural choice is to start with the trained neural decoder for
uniform error rate and only train the site-dependent variables. With this
route, there is a risk that if previously we trained the neural decoder under
uniform error rate for long enough, the neural decoder could learn to ignore
the error rate inputs as they are constant. In this case, only training the
site-dependent variables can fail.

Another route we can take is that we also reinitialize the first
renormalization block. More accurately, we do the following:

  1. 1.

We start with the trained neural decoder for uniform error rate.

  2. 2.

We reverse the belief propagation network in the first renormalization block
to the pre-trained weights. In other words, it now again approximates belief
propagation.

  3. 3.

We then train the site-dependent variables and the first renormalization block
together.

We test this scheme with a distance 16 toric code. For each qubit, there is a
$50\%$ chance it has bit-flip rate $0.16$, and a $50\%$ chance of rate $0$. In
other words, only around half of the qubits are noisy. With this new error
model, we train the site-dependent variables and the first renormalization
block for 4500 batches, where each batch contains 50 training data. The
learning rate of the ADAM optimizer is set to $2\times 10^{-4}$. By performing
this training, we increase the logical accuracy from $0.967$ to $0.993$, where
$0.967$ corresponds to the decoder trained on uniform bit-flip error model.
The accuracy is each evaluated by using $10^{5}$ (syndrome, logical
correction) pair. To provide some comparison, we run the minimum-weight
perfect matching (MWPM) algorithm for the same error model. Without providing
error rate information, the MWPM algorithm assigns an equal weight to each
qubit, and have a logical accuracy of $0.975$. If we provide the perfect
information about the error rates, we can assign weight 1 to noisy qubits, and
weight 100 to the noiseless ones. For each pair of violated parity checks, we
can select the path with the minimum total weight between them, and use this
weight for the MWPM algorithm. With this choice, the logical accuracy rises to
$1$, e.g. $100\%$ success. These accuracies are each evaluated by $10^{4}$
(syndrome, logical correction) pair.

Based on the thoughts above, it is likely better to not start with the neural
decoder trained on uniform error rate. Instead, we can train a neural decoder
with training data that has varying error rates, but otherwise the same
procedure as depicted in Figure 3. This way, the first renormalization block
will not learn to ignore the error rate inputs.


# Characterisation of Laser Wakefield Acceleration Efficiency with Octave
Spanning Near-IR Spectrum Measurements

M.J.V. Streeter<EMAIL_ADDRESS>The Cockcroft Institute, Keckwick
Lane, Daresbury, WA4 4AD, United Kingdom  Physics Department, Lancaster
University, Lancaster LA1 4YB, United Kingdom  The John Adams Institute for
Accelerator Science, Imperial College London, London, SW7 2AZ, UK  Y. Ma
Center for Ultrafast Optical Science, University of Michigan, Ann Arbor, MI
48109-2099, USA  The Cockcroft Institute, Keckwick Lane, Daresbury, WA4 4AD,
United Kingdom  Physics Department, Lancaster University, Lancaster LA1 4YB,
United Kingdom  B. Kettle  The John Adams Institute for Accelerator Science,
Imperial College London, London, SW7 2AZ, UK  S.J.D. Dann  The Cockcroft
Institute, Keckwick Lane, Daresbury, WA4 4AD, United Kingdom  Physics
Department, Lancaster University, Lancaster LA1 4YB, United Kingdom  E.
Gerstmayr  The John Adams Institute for Accelerator Science, Imperial College
London, London, SW7 2AZ, UK  F. Albert  Lawrence Livermore National Laboratory
(LLNL), P.O. Box 808, Livermore, California 94550, USA  N. Bourgeois  Central
Laser Facility, STFC Rutherford Appleton Laboratory, Didcot OX11 0QX, UK  S.
Cipiccia  Diamond Light Source, Harwell Science and Innovation Campus, Fermi
Avenue, Didcot OX11 0DE, UK  J.M. Cole  The John Adams Institute for
Accelerator Science, Imperial College London, London, SW7 2AZ, UK  I.G.
González  Department of Physics, Lund University, P.O. Box 118, S-22100, Lund,
Sweden  A. Higginbotham  York Plasma Institute, Department of Physics,
University of York, York YO10 5DD, UK  A.E. Hussein  Center for Ultrafast
Optical Science, University of Michigan, Ann Arbor, MI 48109-2099, USA
Department of Electrical and Computer Engineering, University of Alberta, 9211
116 Street NW Edmonton, Alberta, T6G 1H9, Canada  D.A. Jaroszynski  SUPA,
Department of Physics, University of Strathclyde, Glasgow G4 0NG, UK  The
Cockcroft Institute, Keckwick Lane, Daresbury, WA4 4AD, United Kingdom  K.
Falk  ELI-Beamlines, Institute of Physics, Academy of Sciences of the Czech
Republic, 18221 Prague, Czech Republic  K. Krushelnick  Center for Ultrafast
Optical Science, University of Michigan, Ann Arbor, MI 48109-2099, USA  N.
Lemos  Lawrence Livermore National Laboratory (LLNL), P.O. Box 808, Livermore,
California 94550, USA  N.C. Lopes  The John Adams Institute for Accelerator
Science, Imperial College London, London, SW7 2AZ, UK  GoLP/Instituto de
Plasmas e Fusão Nuclear, Instituto Superior Técnico, U.L., Lisboa 1049-001,
Portugal  C. Lumsdon  York Plasma Institute, Department of Physics, University
of York, York YO10 5DD, UK  O. Lundh  Department of Physics, Lund University,
P.O. Box 118, S-22100, Lund, Sweden  S.P.D. Mangles  The John Adams Institute
for Accelerator Science, Imperial College London, London, SW7 2AZ, UK  Z.
Najmudin  The John Adams Institute for Accelerator Science, Imperial College
London, London, SW7 2AZ, UK  P.P. Rajeev  Central Laser Facility, STFC
Rutherford Appleton Laboratory, Didcot OX11 0QX, UK  M. Shahzad  SUPA,
Department of Physics, University of Strathclyde, Glasgow G4 0NG, UK  The
Cockcroft Institute, Keckwick Lane, Daresbury, WA4 4AD, United Kingdom  M.
Smid  ELI-Beamlines, Institute of Physics, Academy of Sciences of the Czech
Republic, 18221 Prague, Czech Republic  R. Spesyvtsev  SUPA, Department of
Physics, University of Strathclyde, Glasgow G4 0NG, UK  The Cockcroft
Institute, Keckwick Lane, Daresbury, WA4 4AD, United Kingdom  D.R. Symes
Central Laser Facility, STFC Rutherford Appleton Laboratory, Didcot OX11 0QX,
UK  G. Vieux  SUPA, Department of Physics, University of Strathclyde, Glasgow
G4 0NG, UK  The Cockcroft Institute, Keckwick Lane, Daresbury, WA4 4AD, United
Kingdom  A.G.R. Thomas<EMAIL_ADDRESS>Center for Ultrafast Optical Science,
University of Michigan, Ann Arbor, MI 48109-2099, USA  The Cockcroft
Institute, Keckwick Lane, Daresbury, WA4 4AD, United Kingdom  Physics
Department, Lancaster University, Lancaster LA1 4YB, United Kingdom

###### Abstract

We report on high efficiency energy transfer in a GeV-class laser wakefield
accelerator. Both the transfer of energy from the laser to the plasma
wakefield, and from the plasma to the accelerated electron beam were diagnosed
experimentally by simultaneous measurement of the deceleration of laser
photons and the accelerated electrons as a function of acceleration length.
The extraction efficiency, which we define as the ratio of the energy gained
by the electron beam to the energy lost by the self-guided laser mode, was
maximised at $27\pm 2$% by tuning of the plasma density, plasma length and
incident laser pulse compression. At higher densities, the laser was observed
to fully redshift over an entire octave, from 800 nm to 1600 nm.

††preprint: APS/123-QED

Intense laser pulses can drive compact plasma-based electron accelerators
using a process known as Laser WakeField Acceleration (LWFA). As the laser
pulse propagates through a plasma, it drives electron oscillations that
produce large electrostatic fields, typically of order 100 GV/m. LWFA has been
successfully used to accelerate electrons to $>1$ GeV energy levels over
interaction distances on the order of a centimeter Leemans  _et al._ (2006);
Kneip  _et al._ (2009); Clayton  _et al._ (2010); Wang  _et al._ (2013);
Leemans  _et al._ (2014); Gonsalves  _et al._ (2019). A crucial consideration
for LWFAs is the efficiency of energy transfer from the laser to the
accelerated particle bunch. In radio-frequency (RF) linear accelerators,
efficient operation is achieved by storing the drive energy in a high quality-
factor cavity, which is then extracted by multiple electron beams in a bunch
train. In high-amplitude plasma-accelerators non-linearities eventually damp
out the plasma oscillations and so high efficiency energy transfer must be
achieved within a relatively small number of plasma oscillation periods.

In beam-driven plasma wakefield acceleration (PWFA), the extraction efficiency
is simply calculated as the ratio of the energy gained by the witness beam, to
the energy lost by the driver. Using this measure, an efficiency of $>30$% has
been observed experimentally Litos  _et al._ (2014). In LWFA, energy transfer
to the plasma wakefield occurs through redshifting of the driving laser pulse,
and so can be determined from spectral measurement of the post-interaction
laser pulse Esarey  _et al._ (2009); Shiraishi  _et al._ (2013). Combined with
measurement of the accelerated electron beam spectrum, it is possible to
simultaneously diagnose the efficiency with which the laser excites the plasma
wakefield, and the efficiency with which the electron beam extracts that
energy. Higher-order laser modes that are not guided in a central filament
will not drive strong plasma waves and therefore do not transfer significant
energy to the wake Vieira  _et al._ (2012). Consequently, the ‘extraction
efficiency’ for LWFA only includes energy transfer from the guided, and
therefore redshifted, laser mode. Measuring the electron beam properties and
the input laser energy alone does not provide enough information to infer this
quantity.

Regardless of the nature of the driver, 100% extraction efficiency would
require that the wake of the witness beam perfectly cancels the plasma wake
generated by the driver. With a suitably chosen trapezoidal electron beam
current profile Katsouleas  _et al._ (1987), the accelerating field over the
electron bunch can be kept constant at $E_{z}(\xi_{S})$ where $\xi_{S}$ is the
location of the head of the electron bunch in the co-moving frame $\xi=z-ct$.
Doing so allows for simultaneous high plasma wake extraction efficiency and
low energy spread for the accelerated beam in linear Katsouleas  _et al._
(1987) or non-linear blowout Gordienko and Pukhov (2005); Lu  _et al._ (2007);
Tzoufras  _et al._ (2008) regimes. If dephasing occurs, then this ideal beam-
loading condition can not be maintained, leading to increased energy spread
and lower overall efficiency. Any modification to the wakefield amplitude,
i.e. as the laser evolves, will affect both the ideal beam-loading condition
and the dephasing rate.

In this letter, we present experimental measurements of the transfer of laser
energy into a plasma wakefield and the efficiency with which that energy was
extracted by an electron beam as it was accelerated to $>1$ GeV. This required
measurement of the full laser pulse spectrum at the exit of the plasma, which
in this regime extended up to 1600 nm. Studying the energy transfer between
the laser, plasma and electron beam as a function of the plasma length was
used to reveal the dynamics responsible for this optimum.

For the driving laser of a LWFA, the energy loss per unit length is given by
$-\mathrm{d}W_{L}/\mathrm{d}z=-(W_{L0}/\omega_{0})\mathrm{d}\langle\omega\rangle/\mathrm{d}z$,
where $W_{L0}$ and $\omega_{0}$ are the initial laser pulse energy and
frequency, respectively. The electron-beam energy-gain per-unit-length,
$N_{B}m_{e}c^{2}d\langle\gamma\rangle/dz$, where $N_{B}$ is the number of
accelerated electrons and $\langle\gamma\rangle=\int
S(\gamma)\mathrm{d}\gamma/N_{B}$ is the average energy of the beam, can be
inferred from measurements of the electron spectrum $S(\gamma)$ as a function
of plasma length. Therefore, the instantaneous energy extraction efficiency
for a LWFA can be written as,

$\tilde{\eta_{b}}=-\frac{N_{B}m_{e}c^{2}\omega_{0}}{W_{L0}}\left[\frac{\mbox{d}\langle\gamma\rangle}{\mbox{d}z}\right]\left[\frac{\mbox{d}\langle\omega\rangle}{\mbox{d}z}\right]^{-1}\;.$
(1)

Due to contributions by dephasing, drive laser evolution, beam injection and
beam loss (changing $N_{B}$), the extraction efficiency is not a constant but
changes along the accelerator length. For the results of this paper, we
measure the accelerator averaged extraction efficiency $\eta_{b}$, i.e. the
ratio of the total energy gained by the electron beam to the energy lost by
the laser pulse.

An experiment was performed (see appendix for the experimental schematic) with
the Gemini laser at the Central Laser Facility. Each pulse contained $6.3\pm
0.6$ J in a pulse length of $\tau_{\rm FWHM}=45\pm 4$ fs, giving a peak power
of $P_{0}=134\pm 22$ TW. The linearly polarized pulse was focused to a spot
width of $40(\pm 2)$ $\mu$m $\times 50(\pm 2)$ $\mu$m (FWHM) with an $f/40$
parabolic mirror. A deformable mirror was used to optimize the wavefront,
giving a peak intensity in vacuum of $I_{0}=4.6(\pm 0.8)\times 10^{18}$ W cm-2
and a peak normalised vector potential $a_{0}=1.46\pm 0.12$ at focus.

The laser pulse was focused into a 3D printed two-stage gas cell Hussein  _et
al._ (2019), filled with a 2%/98% nitrogen/helium mix for the first ‘injector’
stage and pure helium in the second ‘accelerator’ stage. The cell walls had 1
mm wide vertical slits to allow for the gas cell to be translated vertically.
This enabled the accelerator length to be adjusted continuously as the exit
wall was angled at 45 ∘ to the vertical plane. The injector stage has an
internal length of 3 mm and the accelerator length was variable over 8-21 mm,
giving a total gas cell length 14-26 mm (including the cell boundaries). The
electron density in the gas cells was varied in the range $n_{e}=0-2.6\times
10^{18}$ cm-3, which was diagnosed by observing the spectrum of Raman side-
scattering from plasma waves generated by low intensity ($a_{0}<1$), long
duration ($\tau_{\mathrm{FWHM}}\approx 200$ fs) laser pulses Matsuoka  _et
al._ (2010).

After interaction with the plasma, the transmitted laser pulse was reflected
from two glass plates into a fiber coupler, sampling a 1 cm diameter region. A
fiber splitter directed the signal onto two spectrometers, one measuring
350–840 nm (Andor Shamrock) and one measuring 900–1700 nm (Ocean Optics
NIRQuest 512). The relative spectral sensitivities of the laser spectrometers
was calibrated using a pre-calibrated white light source. The electron beam
spectrum was measured using a magnetic dipole with integrated field strength
$\int B\mathrm{d}z=0.45$ Tm, which dispersed electrons in the energy range
450–3000 MeV onto a Lanex scintillator.

In order to determine the optimal conditions for electron generation, the gas
cell was positioned at its longest length and the plasma density $n_{e}$,
focus position $z_{f}$ and second order spectral phase $\beta^{(2)}$ of the
laser were independently scanned. The results of a gas cell density scan, for
an accelerator length of 21 mm (total plasma length 27 mm) and the optimal
values of $z_{f}$ and $\beta^{(2)}$, are plotted in figure 1a-b.

Figure 1:  a) Laser and b) electron spectrum at the exit of the 27 mm gas cell
as functions of electron density in the acceleration stage. c) Laser and d)
electron beam spectra at the exit of a gas cell as functions of accelerator
stage length at $n_{e}=1.25\pm 0.06\times 10^{18}$ cm-3. The red-bordered
region in the laser spectra plots indicates the gap between the measurement
ranges of the two spectrometers, which was filled by a Gaussian process
regression interpolation model. Each column of each image shows the average
from 3–5 shots at the same conditions. The width of each column is
proportional to the bin-width between measurements.

The laser spectrum, shown in figure 1a, was increasingly redshifted and
broadened at higher plasma density. For the highest plasma density
$n_{e}=2.6\times 10^{18}$ cm-3, the laser spectrum extended to the limit of
the spectrometer with a peak occurring at 1600 nm. Comparatively little laser
energy was blueshifted for the full density range demonstrating that
ionization blueshift or photon acceleration at the rear of the plasma wave
were not significant Schreiber  _et al._ (2010).

As shown in figure 1b, the highest electron energy occurred for
$n_{e}=1.5\times 10^{18}$ cm-3, where a peak in the spectrum was observed at
$1.4$ GeV. At higher densities, the maximum electron energies decreased, while
the total measured charge remained approximately constant at $eN_{B}\approx
250$ pC. The increasing laser redshift indicates that a strong plasma wave
continued to be driven at these high densities, but the injected electrons
experienced less acceleration.

Figure 1c-d show the result of scanning the accelerator length for a fixed
density of of $n_{e}=1.25\times 10^{18}$ cm-3. The laser (figure 1c redshifted
at an approximately linear rate as the acceleration length was increased. The
electron spectra (figure 1d, shows two distinctive electron bunches were
accelerated, with the higher energy component reaching 1.2 GeV. For plasma
length $z<18$ mm (acceleration cell length $z<13$ mm) the charge in the higher
energy component was much reduced, indicating that injection was sensitive to
small changes in the plasma profile as the cell was translated vertically.
Also for the longest plasma lengths ($z>25$ mm), the laser redshift was
reduced, indicating that less energy was coupled into the plasma. This was
likely due to obstructions to the laser path at the top of the gas cell slits.

Both the laser and electron beam energies were calculated by integrating the
measurements in figure 1 over the spectral axes. For the laser spectrum, there
was a small gap between the ranges of the two spectrometers. In order to
determine the interacting laser energy, Gaussian process regression (GPR) was
used to fit the observed signal and interpolate over this region. The relative
error of the interacting laser energy measurement introduced by this procedure
was calculated from the GPR model uncertainties as less than $3$%. Also, as
the laser spectrum was measured on a small sample region of the transmitted
beam ($\approx 1$ mrad), we have assumed that the spectrum is independent of
angle. Analysis of transmitted laser spectra from the PIC simulations
(described later in this paper) indicate that this assumption leads to the
laser depletion being overestimated by $\approx 4$%. Finally, it was assumed
that the total photon number of the driving laser pulse
($N_{ph}=W_{L0}/\hbar\omega_{0}$) was conserved throughout the interaction
Esarey  _et al._ (2009), such that the energy loss can be calculated from the
change in average laser frequency as in equation 1.

The laser pulse energy loss and electron energy gain are plotted as functions
of plasma density in figure 2a. The laser lost more energy for increasing
plasma density until reaching a plateau for $n_{e}>1.6\times 10^{18}$ cm-3.
The electron beam total energy reached a maximum of $0.23\pm 0.04$ J at
$n_{e}=1.6\times 10^{18}$ cm-3 and was lower for both higher and lower plasma
densities. The extraction efficiency $\eta_{b}=-W_{B}/\Delta W_{L}$, plotted
in figure 2b, reached a maximum of $19\pm 2$% at $n_{e}=1.05\times 10^{18}$
cm-3. Increasing the plasma density beyond this point resulted in a lower
extraction efficiency, even though the electron beam energy increased.

Figure 2: a) and c) Electron beam energy (red) and energy loss of the laser
(black) and b) and d) extraction efficiency as functions of plasma density and
length.

The laser energy loss and electron beam energy for the length scan are shown
in figure 2c. The laser energy loss was approximately linearly proportional to
the plasma length up to $z=25$ mm. The electron beam total energy increased
suddenly with the appearance of the higher energy feature in the electron
spectrum (visible in figure 1d), at $z=18$ mm. Over the range $(20\leq z\leq
26)$ mm, the electron charge was approximately constant at $eN_{B}=280\pm 90$
pC, while the electron beam total energy increased with $z$. The extraction
efficiency, shown in figure 2d), also increased with length for $z<20$ mm as
more charge was injected, before stabilizing at an average of $23\pm 4$% for
$(20\leq z\leq 26)$ mm with a maximum of $27\pm 2$%.

In order to explore the dynamics of the experiment, we performed quasi-3D PIC
simulations using FBPIC Lehe  _et al._ (2016), (for details see appendix). The
laser pulse was initialised using the experimentally measured temporal
profile, and a Gaussian approximation to the measured focal spot energy
distribution. Choosing the vacuum $a_{0}=1.46$ resulted in a larger redshift
than was measured experimentally, indicating that the combination of laser
pulse and target imperfections resulted in a lower proportion of the laser
energy being guided than for a pure Gaussian mode, in line with previous
observations Mangles  _et al._ (2012); Vieira  _et al._ (2012). Using a
Gaussian pulse with $a_{0}=1.3$ provided a good match between simulated and
measured laser spectra. This pulse contained 3.4 J of energy, indicating that
$54$ % of the laser energy contributed to the LWFA in the experiment.

The laser and electron spectra, as functions of propagation distance within
the PIC simulation, are shown in figure 3a and b. The laser spectrum was seen
to redshift and broaden by the same amount as the experimental measurements of
figure 1c, with a redshifted peak at $\omega=1.9\times 10^{15}$ rad fs-1. The
laser reached a peak normalised vector potential of $a_{0}=2.6$ at $z=4.6$ mm,
due to the effects of self-focusing and self-compression. For the rest of the
accelerator, $a_{0}<2.6$ and no self-injection was observed. The accelerated
electron beam resulted from trapping of inner shell electrons from nitrogen
(ionisation injection Rowlands-Rees  _et al._ (2008); Pak  _et al._ (2010);
McGuffey  _et al._ (2010); Chen  _et al._ (2012)) during the first 5 mm. A
total beam charge of 144 pC was observed in the accelerated beam at the end of
the simulation, with a peak in energy of 1.2 GeV, matching the experimental
observations.

Figure 3: a) Laser and b) electron spectra as functions of propagation
distance taken from a PIC simulation with a plateau density of
$n_{e}=1.25\times 10^{18}$ cm-3. c) The extraction efficiency, calculated as
the ratio of the electron beam energy to the laser energy loss, for
simulations with $n_{e}=(1.25,1.5)\times 10^{18}$ cm-3.

The extraction efficiency is shown in Figure 3c for $n_{e}=1.25\times 10^{18}$
cm-3 and for a slightly higher density $n_{e}=1.5\times 10^{18}$ cm-3. For the
higher density case, the laser intensity reached a larger value during the
initial self-focusing, resulting in the trapping of $\sim 50$% more charge.
Also for the higher density, the laser pulse maintained $a_{0}>2.5$ over the
full acceleration distance, thereby driving a higher amplitude wake. As a
result, beam-loading was less severe and so the extraction efficiency was
initially higher than for the lower density case. However, for the higher
density case the efficiency dropped significantly throughout the accelerator,
as the electron beam dephased and started to decelerate for $z>17$ mm.

Simulations using the same plasma profile with plateau densities of
$n_{e}=(1.1,1.2,1.3,1.4)\times 10^{18}$ cm-3, show that both the maximum final
electron energy and extraction efficiencies were optimised for
$n_{e}=1.25\times 10^{18}$ cm-3. Operating the LWFA at a higher density but
over a shorter length allows for higher efficiency, but results in a lower
maximum electron energy.

Figure 4 shows the electron beam position relative to the axial longitudinal
electric field in the lab frame with a Galilean coordinate transform to the
linear group velocity of the laser $\xi=z-v_{g}t$. The fields are shown both
with (figure 4a) and without (figure 4b) the contribution of the trapped
electron bunch. After the initial self-focusing phase, the laser pulse $a_{0}$
steadily dropped from $a_{0}=2$ to $a_{0}=1.2$ at $z=22$ mm, causing a gradual
reduction in the wavelength of the wakefield $\lambda_{p}$. In addition, as
the shape of the wakefield smoothly changed from the saw-tooth profile of the
highly non-linear regime to a quasi-linear sinusoidal profile, this also had
the effect of moving the position of peak accelerating field further forward.
As shown in figure 4b, these effects combined to balance the subluminal group
velocity of the driving laser pulse and so the trapped highly relativistic
electrons closely tracked the position of highest field strength. Figure 4a,
shows the effect of beam-loading, which reduced the accelerating field
experienced by the rear of the bunch to $\approx 15$% of the value at the
head.

Figure 4: a) Loaded and b) unloaded longitudinal electric field and injected
electron bunch position (green) over the propagation axis $z$ of PIC
simulations with $n_{e}=1.25\times 10^{18}$ cm-3. The unloaded fields were
extracted from a simulation without the ionization injection species. The
maximum accelerating field position in the first plasma period (black line) is
overlaid.

During injection, the normalised vector potential reached a peak value of
$a_{0}=2.6$, significantly above the threshold for ionisation injection.
Electrons injected at this point were trapped on _deep_ orbits Chen  _et al._
(2012) in which electrons obtained the wakefield phase velocity $v=v_{\phi}$
significantly before the back of the wake, and so never experienced the
maximum accelerating wakefield. From this point on, a non-evolving laser
driver would have caused the electron beam to dephase at $L_{\phi}=14.5$ mm
with a maximum energy of 800 MeV. However, the driver evolution acted to
mitigate dephasing, resulting in the higher observed electron beam energy of
1.2 GeV.

Methods for accelerating electrons to energies beyond the dephasing limit have
been explored, including non-uniform plasma profiles Pukhov and Kostyukov
(2008); Guillaume  _et al._ (2015); Ma  _et al._ (2018); Sadler  _et al._
(2020) or by using alternative laser focusing geometries with spatio-temporal
couplings Debus  _et al._ (2019); Palastro  _et al._ (2020); Caizergues  _et
al._ (2020). Phase-locked LWFA dynamics in a constant density plasma have been
observed in PIC simulations previously Li  _et al._ (2014) although in that
case it was attributed to pulse depletion. As depletion increases the
wakefield amplitude due to laser redshifting ($a_{0}\propto\omega^{-1/2}$),
this effect alone would actually increase the wakefield wavelength of and
cause the electron beam to dephase more rapidly. Here we show that through
careful management of the laser evolution it is possible to mitigate
dephasing, so that pulse depletion determines the electron energy limit for
the accelerator.

In conclusion, we have measured the extraction efficiency for an LWFA, which
reached a maximum of $27\pm 2$%, close to that previously observed in
electron-beam-driven PWFA Litos  _et al._ (2014). The measurements indicated
that only approximately 20% of the laser pulse energy was transferred to the
plasma wakefield, with approximately half of the laser energy wasted due to a
non-ideal focal spot. The overall efficiency could therefore be increased by
improving the spatial distribution of the laser pulse, as indicated by PIC
simulations with a Gaussian distribution as well as previous studies on the
effects of non-Gaussian focal spots Vieira  _et al._ (2012). Laser energy that
remains after the interaction could possibly be recovered to further improve
the total efficiency of a LWFA facility. For the highest plasma densities and
longest interaction length, the laser pulse spectrum was observed to span a
complete octave from 800–1600 nm. Further harnessing of these effects may open
up a route to relativistic intensity single cycled mid-IR laser pulses
Streeter  _et al._ (2018); Nie  _et al._ (2018, 2020).

###### Acknowledgements.

We acknowledge support from the UK STFC core grants ST/P002056/1 (Cockcroft
Institute), ST/P000835/1 (John Adams Institute), the National Science
Foundation (Grant No. 1804463) and the Air Force Office of Scientific Research
(Grant No. FA9550-16-1-0121). F.A. acknowledges funding from the DOE Early
Career research program (Fusion Energy Sciences SCW1575-1).

## References

  * Leemans  _et al._ (2006) W. P. Leemans, B. Nagler, A. J. Gonsalves, C. Tóth, K. Nakamura, C. G. R. Geddes, E. Esarey, C. B. Schroeder, and S. M. Hooker, Nature Physics 2, 696 (2006).
  * Kneip  _et al._ (2009) S. Kneip, S. R. Nagel, S. F. Martins, S. P. D. Mangles, C. Bellei, O. Chekhlov, R. J. Clarke, N. Delerue, E. J. Divall, G. Doucas, K. Ertel, F. Fiuza, R. Fonseca, P. Foster, S. J. Hawkes, C. J. Hooker, K. Krushelnick, W. B. Mori, C. A. J. Palmer, K. Ta Phuoc, P. P. Rajeev, J. Schreiber, M. J. V. Streeter, D. Urner, J. Vieira, L. O. Silva, and Z. Najmudin, Physical Review Letters 103, 035002 (2009).
  * Clayton  _et al._ (2010) C. E. Clayton, J. E. Ralph, F. Albert, R. A. Fonseca, S. H. Glenzer, C. Joshi, W. Lu, K. A. Marsh, S. F. Martins, W. B. Mori, A. Pak, F. S. Tsung, B. B. Pollock, J. S. Ross, L. O. Silva, and D. H. Froula, Physical Review Letters 105, 105003 (2010).
  * Wang  _et al._ (2013) X. Wang, R. Zgadzaj, N. Fazel, Z. Li, S. A. Yi, X. Zhang, W. Henderson, Y.-Y. Chang, R. Korzekwa, H.-E. Tsai, C.-H. Pai, H. Quevedo, G. Dyer, E. Gaul, M. Martinez, A. C. Bernstein, T. Borger, M. Spinks, M. Donovan, V. Khudik, G. Shvets, T. Ditmire, and M. C. Downer, Nature Communications 4, 1988 (2013).
  * Leemans  _et al._ (2014) W. P. Leemans, A. J. Gonsalves, H.-S. Mao, K. Nakamura, C. Benedetti, C. B. Schroeder, C. Tóth, J. Daniels, D. E. Mittelberger, S. S. Bulanov, J.-L. Vay, C. G. R. Geddes, and E. Esarey, Physical Review Letters 113, 245002 (2014).
  * Gonsalves  _et al._ (2019) A. J. Gonsalves, K. Nakamura, J. Daniels, C. Benedetti, C. Pieronek, T. C. H. de Raadt, S. Steinke, J. H. Bin, S. S. Bulanov, J. van Tilborg, C. G. R. Geddes, C. B. Schroeder, C. Tóth, E. Esarey, K. Swanson, L. Fan-Chiang, G. Bagdasarov, N. Bobrova, V. Gasilov, G. Korn, P. Sasorov, and W. P. Leemans, Physical Review Letters 122, 084801 (2019).
  * Litos  _et al._ (2014) M. Litos, E. Adli, W. An, C. I. Clarke, C. E. Clayton, S. Corde, J. P. Delahaye, R. J. England, A. S. Fisher, J. Frederico, S. Gessner, S. Z. Green, M. J. Hogan, C. Joshi, W. Lu, K. A. Marsh, W. B. Mori, P. Muggli, N. Vafaei-Najafabadi, D. Walz, G. White, Z. Wu, V. Yakimenko, and G. Yocky, Nature 515, 92 (2014).
  * Esarey  _et al._ (2009) E. Esarey, C. B. Schroeder, and W. P. Leemans, Reviews of Modern Physics 81, 1229 (2009).
  * Shiraishi  _et al._ (2013) S. Shiraishi, C. Benedetti, A. J. Gonsalves, K. Nakamura, B. H. Shaw, T. Sokollik, J. van Tilborg, C. G. R. Geddes, C. B. Schroeder, C. Tóth, E. Esarey, and W. P. Leemans, Physics of Plasmas 20, 063103 (2013).
  * Vieira  _et al._ (2012) J. Vieira, S. F. Martins, F. Fiúza, C. K. Huang, W. B. Mori, S. P. Mangles, S. Kneip, S. Nagel, Z. Najmudin, and L. O. Silva, Plasma Physics and Controlled Fusion 54, 055010 (2012).
  * Katsouleas  _et al._ (1987) T. Katsouleas, S. Wilks, P. Chen, J. M. Dawson, and J. Su, Particle Accelerators 22, 81 (1987).
  * Gordienko and Pukhov (2005) S. Gordienko and A. Pukhov, Physics of Plasmas 12, 043109 (2005).
  * Lu  _et al._ (2007) W. Lu, M. Tzoufras, C. Joshi, F. Tsung, W. Mori, J. Vieira, R. Fonseca, and L. Silva, Physical Review Special Topics - Accelerators and Beams 10, 061301 (2007).
  * Tzoufras  _et al._ (2008) M. Tzoufras, W. Lu, F. S. Tsung, C. Huang, W. B. Mori, T. Katsouleas, J. Vieira, R. A. Fonseca, and L. O. Silva, Physical Review Letters 101, 145002 (2008).
  * Hussein  _et al._ (2019) A. E. Hussein, N. Senabulya, Y. Ma, M. J. V. Streeter, B. Kettle, S. J. D. Dann, F. Albert, N. Bourgeois, S. Cipiccia, J. M. Cole, O. Finlay, E. Gerstmayr, I. G. González, A. Higginbotham, D. A. Jaroszynski, K. Falk, K. Krushelnick, N. Lemos, N. C. Lopes, C. Lumsdon, O. Lundh, S. P. D. Mangles, Z. Najmudin, P. P. Rajeev, C. M. Schlepütz, M. Shahzad, M. Smid, R. Spesyvtsev, D. R. Symes, G. Vieux, L. Willingale, J. C. Wood, A. J. Shahani, and A. G. R. Thomas, Scientific Reports 9, 3249 (2019).
  * Matsuoka  _et al._ (2010) T. Matsuoka, C. McGuffey, P. G. Cummings, Y. Horovitz, F. Dollar, V. Chvykov, G. Kalintchenko, P. Rousseau, V. Yanovsky, S. S. Bulanov, A. G. R. Thomas, A. Maksimchuk, and K. Krushelnick, Physical Review Letters 105, 034801 (2010).
  * Schreiber  _et al._ (2010) J. Schreiber, C. Bellei, S. P. D. Mangles, C. Kamperidis, S. Kneip, S. R. Nagel, C. A. J. Palmer, P. P. Rajeev, M. J. V. Streeter, and Z. Najmudin, Physical Review Letters 105, 235003 (2010).
  * Lehe  _et al._ (2016) R. Lehe, M. Kirchen, I. A. Andriyash, B. B. Godfrey, and J.-L. Vay, Computer Physics Communications 203, 66 (2016).
  * Mangles  _et al._ (2012) S. P. D. Mangles, G. Genoud, M. S. Bloom, M. Burza, Z. Najmudin, A. Persson, K. Svensson, A. G. R. Thomas, and C.-G. Wahlström, Physical Review Special Topics - Accelerators and Beams 15, 011302 (2012).
  * Rowlands-Rees  _et al._ (2008) T. P. Rowlands-Rees, C. Kamperidis, S. Kneip, A. J. Gonsalves, S. P. D. Mangles, J. G. Gallacher, E. Brunetti, T. Ibbotson, C. D. Murphy, P. S. Foster, M. J. V. Streeter, F. Budde, P. A. Norreys, D. A. Jaroszynski, K. Krushelnick, Z. Najmudin, and S. M. Hooker, Physical Review Letters 100, 105005 (2008).
  * Pak  _et al._ (2010) A. Pak, K. A. Marsh, S. F. Martins, W. Lu, W. B. Mori, and C. Joshi, Physical Review Letters 104, 025003 (2010).
  * McGuffey  _et al._ (2010) C. McGuffey, A. G. R. Thomas, W. Schumaker, T. Matsuoka, V. Chvykov, F. J. Dollar, G. Kalintchenko, V. Yanovsky, A. Maksimchuk, K. Krushelnick, V. Y. Bychenkov, I. V. Glazyrin, and a. V. Karpeev, Physical Review Letters 104, 025004 (2010).
  * Chen  _et al._ (2012) M. Chen, E. Esarey, C. B. Schroeder, C. G. R. Geddes, and W. P. Leemans, Physics of Plasmas 19, 033101 (2012).
  * Pukhov and Kostyukov (2008) A. Pukhov and I. Kostyukov, Physical Review E 77, 025401(R) (2008).
  * Guillaume  _et al._ (2015) E. Guillaume, A. Döpp, C. Thaury, K. Ta Phuoc, A. Lifschitz, G. Grittani, J. P. Goddet, A. Tafzi, S. W. Chou, L. Veisz, and V. Malka, Physical Review Letters 115, 155002 (2015).
  * Ma  _et al._ (2018) Y. Ma, D. Seipt, S. J. Dann, M. J. Streeter, C. A. Palmer, L. Willingale, and A. G. Thomas, Physics of Plasmas 25, 113105 (2018).
  * Sadler  _et al._ (2020) J. D. Sadler, C. Arran, H. Li, and K. A. Flippo, Physical Review Accelerators and Beams 23, 021303 (2020).
  * Debus  _et al._ (2019) A. Debus, R. Pausch, A. Huebl, K. Steiniger, R. Widera, T. E. Cowan, U. Schramm, and M. Bussmann, Physical Review X 9, 031044 (2019).
  * Palastro  _et al._ (2020) J. P. Palastro, J. L. Shaw, P. Franke, D. Ramsey, T. T. Simpson, and D. H. Froula, Physical Review Letters 124, 134802 (2020).
  * Caizergues  _et al._ (2020) C. Caizergues, S. Smartsev, V. Malka, and C. Thaury, Nature Photonics 14, 475 (2020).
  * Li  _et al._ (2014) W. Li, J. Liu, W. Wang, Z. Zhang, Q. Chen, Y. Tian, R. Qi, C. Yu, C. Wang, T. Tajima, R. Li, and Z. Xu, Applied Physics Letters 104, 093510 (2014).
  * Streeter  _et al._ (2018) M. J. V. Streeter, S. Kneip, M. S. Bloom, R. A. Bendoyro, O. Chekhlov, A. E. Dangor, A. Döpp, C. J. Hooker, J. Holloway, J. Jiang, N. C. Lopes, H. Nakamura, P. A. Norreys, C. A. J. Palmer, P. P. Rajeev, J. Schreiber, D. R. Symes, M. Wing, S. P. D. Mangles, and Z. Najmudin, Physical Review Letters 120, 254801 (2018).
  * Nie  _et al._ (2018) Z. Nie, C.-H. Pai, J. Hua, C. Zhang, Y. Wu, Y. Wan, F. Li, J. Zhang, Z. Cheng, Q. Su, S. Liu, Y. Ma, X. Ning, Y. He, W. Lu, H.-H. Chu, J. Wang, W. B. Mori, and C. Joshi, Nature Photonics 12, 489 (2018).
  * Nie  _et al._ (2020) Z. Nie, C.-H. Pai, J. Zhang, X. Ning, J. Hua, Y. He, Y. Wu, Q. Su, S. Liu, Y. Ma, Z. Cheng, W. Lu, H.-H. Chu, J. Wang, C. Zhang, W. B. Mori, and C. Joshi, Nature Communications 11, 2787 (2020).

##  I Appendix

###  I.1 Experimental Setup

Figure 5: Experimental schematic showing the gas cell geometry, the electron
spectrometer and the NIR (900–1700 nm) and VIS (350–840 nm) spectrometers
which measured the transmitted laser pulse.

###  I.2 Simulations

Simulations were performed using the quasi-3D Particle-In-Cell (PIC) code
FBPIC (https://fbpic.github.io/index.html) using cylindrical symmetry with
azimuthal mode decomposition. The simulation domain was $100\times 80$ $\mu$m
in the propagation ($z$) and radial $r$ axes, which were divided into
$2000\times 100$ cells with two azimuthal modes. The plasma electrons were
represented by $2\times 2\times 8$ macro-particles per cell in the $z$, $r$
and $\theta$ directions. The plasma was initialized with a density of
$1.25\times 10^{18}$ mm-3, with $500$ $\mu$m linear ramps either side of a 26
mm plateau. The temporal profile of the laser pulse was taken from FROG
measurements of the input laser pulses during the experiment. Measurements of
the laser focal spot were used to fit a 2D Gaussian transverse profile as
$a(x,y,t)=\hat{a}(t)\exp\left[-(x/\sigma_{x})^{2}-(y/\sigma_{y})^{2}\right]$
with $\sigma_{x}=40$ $\mu$m and $\sigma_{y}=32$ $\mu$m.