# Online Obstacle evasion with Space-Filling Curves

Ashay Wakode1, Arpita Sinha2 1 Ashay Wakode is with Mechanical Engineering
Department, Indian Institute Of Technology Bombay, Mumbai, India
<EMAIL_ADDRESS>2 Arpita Sinha is with System and Controls Department;
Indian Institute Of Technology Bombay, Mumbai, India
<EMAIL_ADDRESS>
###### Abstract

The paper presents a strategy for robotic exploration problems using Space-
Filling curves (SFC). The region of interest is first tessellated, and the
tiles/cells are connected using some SFC. A robot follows the SFC to explore
the entire area. However, there could be obstacles that block the systematic
movement of the robot. We overcome this problem by providing an evading
technique that avoids the blocked tiles while ensuring all the free ones are
visited at least once. The proposed strategy is online, implying that prior
knowledge of the obstacles is not mandatory. It works for all SFCs, but for
the sake of demonstration, we use Hilbert curve. We present the completeness
of the algorithm and discuss its desirable properties with examples. We also
address the non-uniform coverage problem using our strategy.

###### Index Terms:

Robotic Exploration, Space-Filling curve, Online Obstacle evasion, Non-uniform
coverage.

##  I Introduction

In 1878, George Cantor demonstrated that an interval $I=[0,1]$ can be mapped
bijectively onto $[0,1]\times[0,1]$. Later, G. Peano discovered one such
mapping that is also continuous and surjective; the image of such mapping when
parameterized in the interval $I$ to higher dimensions ($\mathbb{R}^{n}$) is
known as Space-Filling Curve (SFC). More SFCs were later discovered by E.
Moore, H. Lebesgue, W. Sierpinski, and G. Polya [1, 2]. See Fig 1. SFCs have
some interesting properties -

  * •

Self-Similar \- Each curve is made out of similar sub-curves

  * •

Surjective Map \- SFC passes through every point of $\mathbb{R}^{n}$

  * •

Locality Preserving \- Two points close by in $I$ map to close by points in
$\mathbb{R}^{n}$

Due to the above properties, SFCs have been used in many applications - data
collection from sensor network [3, 4]; ordering meshes of complex geometries
[2] and many more. An approximate solution to Travelling Salesman Problem
(TSP) can be found using Hilbert’s Space Filling curve [5]. Space Filling Tree
analogous to SFCs having tree-like structure have been proposed for sampling-
based path planning [6], as opposed to traditional methods like Rapid-
exploring Random Trees (RRTs) [7].

Figure 1: Hilbert curve and Sierpinski curve

One of the major applications of SFCs is robotic exploration. In robotic
exploration problem, a single or group of robotic agents are deployed to
search, survey or gather information about a specific region while avoiding
obstacles. Robotics exploration is one of the sub-problems of the larger
Coverage Planning Problem (CPP), wherein the agent is bestowed with the task
of visiting all points in a given region or volume while avoiding obstacles
[8, 9, 10]. Numerous approaches for CPP already exists - Graph based, Grid
based, Neural-Network based, Cellular decomposition based [8]. Each of these
approaches can be used for robotic exploration problem.

SFCs have been used for robotics exploration and the reasons making it
eligible are:

  * •

Complete and Robust : Coverage using SFCs is time complete and robust to
failure [11]

  * •

Extensible :

    * –

Dimension: The strategy developed for 2D can extended to 3D since similar
grammar exists for the construction of SFC in both dimensions [2]

    * –

Number of search agents : Coverage using SFCs have been shown to be complete
and robust when multiple agents are used [11]

    * –

Irregular Area : The generalized version of SFCs aka Generalized SFC (GSFC)
can span irregular quadrilateral and triangles as opposed to SFCs which map
regular shapes like square or isosceles right triangle [2]

  * •

Non-Uniform coverage : SFCs can be easily used in non-uniform coverage
scenarios requiring some parts to be searched more rigorously than others. On
top of that Hilbert curve has been shown to be more efficient in time/energy
than the lawnmower’s path [12]

A 2D area with multiple obstacles of arbitrary shape and size is explored
using a SFC. The area is tessellated (referred to as waypoints) and
tiles/cells are connected in the order determined by SFC. The tessellation
with obstacle is considered unreachable. A online strategy is suggested for a
search agent to explore the area using SFC and evade the obstacle on the go.
It suggests modified waypoint when unreachable waypoints are encountered. It
guarantees that all the waypoints reachable from the initial waypoint are
explored. The strategy works for all SFCs, but for the sake of demonstration
Hilbert’s curve is used in this paper.

The rest of this paper is organized into six sections: Section II elaborates
on the related work and the need for a new strategy. Section III discusses
preliminaries on SFCs. Section IV details the strategy, investigates its
completeness and discusses its useful properties. Section V compiles the
simulation result. Section VI concludes the paper and talks about the
limitations and possible future work.

##  II Related Work

There is sizable literature dealing with using SFC in the robotic exploration
problem. [11] suggests the use of a swarm of mobile robots for exploration,
each mobile robot following a SFC. This approach is efficient in terms of
energy, robust to failures and assures coverage in finite time, but considers
an obstacle-free environment. [13] proposes a new method for non-uniform
coverage planning for UAVs. [12] takes the work further and uses the Hilbert
curve for non-uniform coverage which is optimal as opposed to existing
methods. [14] introduced a UAV system to conduct aerial localization that uses
Moore’s SFC. However, [13, 12, 14] does not consider obstacles.

[15] introduced path planning approach for SFCs with obstacles and proved the
optimality for specific obstacle configurations, due to the existence of a
Hamiltonian path for such obstacle configurations. [16] introduced an
algorithm to construct SFCs for sensor networks with holes. The algorithm can
be used for motion planning with obstacles while using SFC. However, the
solutions proposed in [15, 16] require the knowledge of obstacles before
starting the exploration.

[17] formulated an online obstacle evasion strategy for Hilbert curve with
only one waypoint blocked by an obstacle. [18] carries forward the work and
suggests a strategy capable of evading two neighboring waypoints on a Hilbert
curve. [19] talks about online obstacle evasion for aerial vehicle covering a
region using the Sierpinski curve, but the obstacles need to block disjoint
cells.

Early work on the topic did not consider obstacles in the environment. Later
works considered obstacles but required knowledge about the obstacles or were
restricted to a particular class of obstacles and SFC used. Therefore an
online strategy capable of evading arbitrary obstacles while using any SFC is
the required next step.

##  III Preliminaries on SFCs

SFCs are defined as,
Definition : Given a mapping $f:I\rightarrow Q\subset\mathbb{R}^{n}$, with
$f_{*}\textit{(}I)$ as image. $f_{*}\textit{(}I)$ is called a SFC, if
$f_{*}\textit{(}I)$ has Jordan content (area, volume, ..) larger than 0 [2].

In this paper, approximations of SFCs are used. Approximate SFCs are
constructed by dividing $\mathbb{R}^{n}$ and connecting the centers of the
cells by a continuous piecewise straight line. The way the divisions are
connected depends on the grammar of the SFC being used.

Figure 2: Hilbert curve as mapped from $I$ Figure 3: Hilbert curve iterations
1 to 5; Colored circles represent specific translation + rotation rules for
creating $n^{th}$ iteration from $n-1^{th}$ iteration

Here in Fig 2, it is seen that a square is divided N times and the centers are
connected to generate Hilbert curve. In this case, hilbert curve is mapped
onto a square in $\mathbb{R}^{2}$. The degree to which the square is divided
is quantified by the term - ”iteration”. As iteration goes to infinity,
approximate SFC visits every point in the square; this curve is termed as SFC.
Thus explains the surjective mapping property of SFC. A Rigorous mathematical
explanation can be found in [1, 2]. The SFCs in Fig 1, 2 are approximate but
are referred to simply as SFC in literature.

Moreover, approximate SFCs are sufficient for the exploration given that the
robotic agent has a search radius enough to cover the entire division (area or
volume) while at its center.

Furthermore, it can be seen in Fig 3 that iteration $N$ is constructed by
rotation and translation of iteration $N-1$. Hence explaining the self-similar
property. The grammatical way of constructing SFCs uses this fact.

SFCs are Hölder continuous mappings. Given,

$x,y\in I$

and $f$ is SFC with $f:I\rightarrow Q\subset\mathbb{R}^{n}$, if,

$\lVert f(x)-f(y)\rVert_{2}\leq C|x-y|^{r}$

$f$ is said to be Hölder continuous with exponent $r$, if a constant $C$
exists for all $x,y$. The RHS can be interpreted as the distance between $x$
and $y$, while LHS is the distance between points in SFC. Here, the distance
between points on SFC is bounded by the distance between the points they are
mapped from, explaining the Locality Preservation.

Lastly, any problem in Q $\subset\mathbb{R}^{n}$ mapped by SFC can be
transformed into a problem in $I$. And often solving the problem in $I$ is
straightforward than solving it in the original space $Q$. This technique is
referred to as General Space-Filling Heuristics (GSFH).

##  IV Problem Formulation and Proposed Solution

A robotic agent with sensing radius $s$, is given a task to explore a 2D
region using SFC and avoiding the obstacles when detected. The region to be
explored can have any shape and size. But, SFCs only span regular shapes like
squares (Hilbert curve) and isosceles right triangle (Sierpinski curve). Here,
it is assumed that the shape to be explored is regular, since an irregular
shape can be decomposed into regular shapes and explored separately. Extension
of SFCs to General SFCs [2], which map arbitrary quadrilaterals and triangles,
also exists and can be used. For this paper, we look into SFCs only.

The strategy suggested can be used for any SFC. For demonstration the results
are worked out for Hilbert curve. A square with area $A$, has obstacles at
unknown locations. Without loss of generality, let the robotic agent start at
left-hand bottom corner of the square. It follows Hilbert curve to explore the
region while avoiding the static obstacles and get to the right-hand bottom
corner. The obstacle can have arbitrary shape and size. The iteration ($k\in
N$) of Hilbert curve used is such that the agent can scan an entire cell while
at its center. Mathematically,

$k\geq\lceil log_{2}(A/s\sqrt{2}-1)\rceil$ (1)

$Note$ : This relation has been found out by comparing sensing radius and
diagonal of cell for Hilbert curve. Similar relation for other SFCs can be
found out trivially.

Here, the robotic agent detects an obstacle just before going to the cell with
an obstacle. It is also assumed that the agent has the capability to memorize
the locations of detected obstacles.

###  IV-A Proposed Solution

Given an SFC mapping an area, the center of each cell is identified as
waypoint and are numbered from $0$ to $N-1$ based on SFC. The goal of the
robotic agent is to start at $0$ and visit all the unblocked reachable
waypoints, This means going to a waypoint with the maximum number possible.
Certain region may not be blocked by the obstacle but can still be unreachable
when it is blocked by obstacles on all sides, like in the case of annular
shaped obstacle.

Here, a graph theoritic solution is presented. Nomenclature used :

  * •

$G=$ Graph with waypoints of SFC as vertices; Adjacent vertices are connected
by the edge; No waypoint is assumed to be blocked by an obstacle in $G$. $G$
is a dual graph of SFC decomposition, similar to one used in [15].

  * •

$O=$ Set of vertices where obstacles have been detected

  * •

$V=$ Set of visited vertices

  * •

$A(V,G)=$ Set of vertices in $G$ adjacent to vertices in $V$ but not present
in $V$

The strategy is presented as a pseudocode in Algorithm 1. The agent uses the
strategy at each waypoint to decide the next waypoint. While at waypoint $c$,
the agent knows $G$ given the SFC, $V$ from the visited nodes and $O$ from the
detected obstacles till that time, and needs to decide about next waypoint to
visit. If all the waypoints adjacent to $V$ in $G$ are blocked, no waypoint
remains to be visited and the search can be terminated (line $4$). Suppose
there exists a set of waypoints that can be visited. In that case, the
strategy suggests minimum numbered waypoint $p$, and the shortest path $R$ to
$p$ is found using Dijkstra’s algorithm (any shortest path finding algorithm
can be used). $R$ refers to an array of waypoints starting with $c$ and
terminating at $p$. The agent checks if $p$ is blocked while at the pen-
ultimate waypoint (line $7$). If $p$ is found to be blocked, it is added to
$O$ and the strategy starts again at step $4$ (line $8$ to $10$). Otherwise,
$p$ is reached and added to $V$ (line $12$ and $13$). The algorithm is used
again to decide about the next waypoint. The strategy is explained using an
example in the next section.

Algorithm 1 Online Strategy for Obstacle evasion

1:Input : $G$, $V$, $O$

2:Output : Next waypoint

3:Initialization : Current waypoint ($c$)

4:if $A(V,G)-O\neq\emptyset$ then

5: $p=min(A(V,G)-O)$ $\triangleright$ min numbered element

6: $R=$ shortest route to $p$ starting from $c$

7: while Agent at $R[-2]$ do $\triangleright$ second last element of $R$

8: if $p$ is blocked $\lor$ $p\in O$ then

9: $O\leftarrow O\cup\\{p\\}$

10: Go back to step $4$

11: else

12: Next waypoint $=p$

13: $V\leftarrow V\cup\\{p\\}$

14: end if

15: end while

16:else

17: All the reachable waypoints visited

18:end if

######  Lemma 1.

An agent starting at waypoint $H$ and following the proposed strategy will
visit all the waypoints connected to $H$.

###### Proof.

Let us prove the lemma by contradiction. Assume here a waypoint $J$ (not
blocked by an obstacle) connected to $H$, but the agent terminated the
exploration without visiting $J$. This can happen only if,

$A(V,G)-O=\emptyset$

But,

$A(V,G)-O=\\{J\\}$

Since, $J$ is connected to $H$ and not visited. Therefore, a contradiction. ∎

Here, the agent will explore the area spanned by waypoints given the sensing
radius enough to cover the entire cell (Eq. 1). But, this does not guarantee
complete search.

Figure 4: Region blocked by tight space reachable through the use of higher
iteration (4) of Hilbert curve

In a situation with tight spaces, certain reachable regions may get blocked
due to the iteration of SFC used. In such cases, the iteration of SFC can be
increased. This would yield a better approximation of the obstacle and hence
better coverage. This is illustrated in Fig 4.
$Remark$ : There can be multiple strategies to select the next waypoint among
the connected waypoints. The proposed strategy ensures that SFC is followed in
the regions without obstacles.

###  IV-B Discussion

  * •

The strategy works well with both normal and sparse obstacles, which is not
the case with Lawnmower’s path which performs poorly with sparse obstacles.

  * •

The path modification occurs locally. Given a obstacle, the strategy suggests
path changes closer to it. Meaning SFC retains its original form in the rest
of the region.

  * •

Higher iterations of SFCs are preferred for sparse obstacles while lower
iterations are preferred for obstacles with larger size. Higher iterations
offer a more agile path and less occlusion than lower iterations when dealing
with sparse obstacles. At the same time, they offer a longer path and may not
be desired for larger obstacles.

  * •

The path generated using the strategy can be adapted for different parts of
the region to be explored, given the knowledge of the nature of obstacle
(Sparse or Normal). Further, the agent may want to explore certain parts with
more rigor. In such cases, higher iteration may give the agent more time and
focus the search on smaller cells/tiles [13, 12]. This can be easily achieved
with proposed strategy. An Example of such a situation is introduced in the
next section

Obstacles comparable to the sensing radius of the agent are considered normal.
While the obstacles much smaller than the sensing radius are considered
sparse. Here it is assumed that the agent can modify the sensing radius, as is
the case in [13, 12] with regards to aerial vehicles. Otherwise, the agent can
discard the data collected outside the required sensing radius.

##  V Results

The presented strategy was implemented in Python (version $3.8$). iGraph [21]
was used for graph operations (version $22.0.4$). Hilbert curve library [20]
was used for plotting the Hilbert curve (version $2.0.5$). The code for
implementation and generated figures are available at -
https://github.com/wakodeashay/SKC. Let’s look at examples,

  * •

Normal Obstacles : The strategy was implemented for a situation with normal
obstacles. Figure 5 shows two normal sized obstacles blocking the SFC. The
modified path is also shown.

Here, while at $0$,

$A(V,G)-O=\\{1,2,3\\}$

the agent searches for obstacle and goes to $1$. Next,

$V=\\{0,1\\},A(V,G)-O=\\{2,3,13\\}$

agent searches for obstacle and goes to $2$. When no obstacle is found, the
agent goes to the next waypoint in the SFC due to the minimum rule. This
happens till $21$, when obstacle is detected at $22$. At this point,

$V=\\{0,1,...,21\\}$ $A(V,G)-O=\\{23,29,30,31,32,53,54,57,58\\}$

So, the agent tries to reach $23$ but encounters obstacle while at $20$. Next,
it goes to $29$ from $20$ using shortest path. Now,

$V=\\{0,1,...,21,29\\}$ $A(V,G)-O=\\{24,28,30,31,32,53,54,57,58\\}$

The agent decides to go $24$ but detects an obstacle while at $29$. After
which the agent goes to $28$,

$V=\\{0,1,...,21,29,28\\}$ $A(V,G)-O=\\{27,30,31,32,53,54,57,58\\}$

Hence, the agent goes to $27$ and $26$ is added to $A(V,G)-O$, which is
reached finally. Now,

$V=\\{0,1,...,21,29,28,27,26\\}$ $A(V,G)-O=\\{25,30,31,32,53,54,57,58\\}$

The agent decides to go $25$ and detects obstacle while at $26$. Finally, the
agent goes to $30$, $31$ and so on following the SFC until next obstacle is
detected.

Figure 5: Hilbert curve with normal sized obstacles; Modified path, Brown
represent waypoints with detected obstacle

  * •

Sparse Obstacles : Figure 6 shows three scenarios with obstacle placed
sparsely in the area to be explored.

Figure 6: Obstacle evasion scenarios with sparse obstacles; 9.7% (100), 19.5%
(200) and 29.2% (300) of total waypoints (1024) blocked respectively

  * •

Non-Uniform Coverage : In certain cases user may want to search some parts
more rigorously. This can be done using higher iteration SFC, since the agent
spends more time with higher iteration and the search is focussed on smaller
radii. Nonetheless, if the user knows about the nature of obstacle (sparse or
normal), the iteration number can be increased or reduced on the go. This is
illustrated in 7.

$Remark$ : In certain situations, agent is unable to get to the last point of
the SFC. This is seen in the right-hand lower quadrant of fig 7. The agent
ends his journey at point $A$. Therefore, it does not know if it can directly
move to the first waypoint of the next SFC (could be blocked). In this case,
agent can be moved to a waypoint ($B$) on the shared edge closest to $A$.
Eventually, moving to the closest waypoint on the next SFC and start using the
proposed strategy. The agent will move to the first waypoint if it is
connected to $C$, in which case it, can start the exploration afresh if
required by the user or else continue using the strategy.

Figure 7: Inward red arrow shows starting position, Outward red arrow shows
terminal position

##  VI Conclusion

This paper discusses a strategy for online obstacle evasion on Space-Filling
curve. Preliminaries on Space-Filling curves were introduced, after which the
problem was formally put down. The strategy is presented and the completeness
is also touched upon. Later, it was validated on examples. The presented
strategy is a high-level waypoint planner. The path needs to be smoothened for
agents with non-holonomic constraints. Finally, The strategy can be optimized
to reduce the total number of visits to waypoints. Extension to 3D space,
environment with dynamic obstacles and Multi-agent scenarios could be
interesting future directions.

## References

  * [1] Sagan, H., 2012. Space-filling curves. Springer Science & Business Media.
  * [2] Bader, M., 2012. Space-filling curves: an introduction with applications in scientific computing (Vol. 9). Springer Science & Business Media.
  * [3] Yan, Y. and Mostofi, Y., 2014, December. An efficient clustering and path planning strategy for data collection in sensor networks based on space-filling curves. In 53rd IEEE Conference on Decision and Control (pp. 6895-6901). IEEE.
  * [4] Yan, Y. and Mostofi, Y., 2016. Efficient clustering and path planning strategies for robotic data collection using space-filling curves. IEEE Transactions on Control of Network Systems, 4(4), pp.838-849.
  * [5] Bartholdi III, J.J. and Platzman, L.K., 1982. An O (N log N) planar travelling salesman heuristic based on spacefilling curves. Operations Research Letters, 1(4), pp.121-125.
  * [6] Kuffner, J.J. and LaValle, S.M., 2011, September. Space-filling trees: A new perspective on incremental search for motion planning. In 2011 IEEE/RSJ International Conference on Intelligent Robots and Systems (pp. 2199-2206). IEEE.
  * [7] LaValle, S.M. and Kuffner Jr, J.J., 2001. Randomized kinodynamic planning. The international journal of robotics research, 20(5), pp.378-400.
  * [8] Galceran, E. and Carreras, M., 2013. A survey on coverage path planning for robotics. Robotics and Autonomous systems, 61(12), pp.1258-1276.
  * [9] Choset, H., 2001. Coverage for robotics–a survey of recent results. Annals of mathematics and artificial intelligence, 31(1), pp.113-126.
  * [10] Cabreira, T.M., Brisolara, L.B. and Paulo R, F.J., 2019. Survey on coverage path planning with unmanned aerial vehicles. Drones, 3(1), p.4.
  * [11] Spires, S.V. and Goldsmith, S.Y., 1998, July. Exhaustive geographic search with mobile robots along space-filling curves. In International Workshop on Collective Robotics (pp. 1-12). Springer, Berlin, Heidelberg.
  * [12] Sadat, S.A., Wawerla, J. and Vaughan, R., 2015, May. Fractal trajectories for online non-uniform aerial coverage. In 2015 IEEE international conference on robotics and automation (ICRA) (pp. 2971-2976). IEEE.
  * [13] Sadat, S.A., Wawerla, J. and Vaughan, R.T., 2014, September. Recursive non-uniform coverage of unknown terrains for uavs. In 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems (pp. 1742-1747). IEEE.
  * [14] Z. Liu, Y. Chen, B. Liu, C. Cao and X. Fu, ”HAWK: An Unmanned Mini-Helicopter-Based Aerial Wireless Kit for Localization,” in IEEE Transactions on Mobile Computing, vol. 13, no. 2, pp. 287-298, Feb. 2014, doi: 10.1109/TMC.2012.238.
  * [15] Tiwari, A., Chandra, H., Yadegar, J. and Wang, J., 2007. Constructing optimal cyclic tours for planar exploration and obstacle avoidance: A graph theory approach. In Advances in Cooperative Control and Optimization (pp. 145-165). Springer, Berlin, Heidelberg.
  * [16] Ban, X., Goswami, M., Zeng, W., Gu, X. and Gao, J., 2013, April. Topology dependent space filling curves for sensor networks and applications. In 2013 Proceedings IEEE INFOCOM (pp. 2166-2174). IEEE.
  * [17] S. H. Nair, A. Sinha and L. Vachhani, ”Hilbert’s space-filling curve for regions with holes,” 2017 IEEE 56th Annual Conference on Decision and Control (CDC), 2017, pp. 313-319, doi: 10.1109/CDC.2017.8263684.
  * [18] Joshi, A.A., Bhatt, M.C. and Sinha, A., 2019, December. Modification of Hilbert’s Space-Filling Curve to Avoid Obstacles: A Robotic Path-Planning Strategy. In 2019 Sixth Indian Control Conference (ICC) (pp. 338-343). IEEE.
  * [19] Wakode, A. and Sinha, A. (2022) ‘Online Evasive Strategy for Aerial Survey using Sierpinski curve’, IFAC-PapersOnLine, 55(22), pp. 129–134. doi: 10.1016/j.ifacol.2023.03.022.
  * [20] Altay, G. (2021). HilbertCurve: Python implementation of the Hilbert curve.
  * [21] Csardi, G., & Nepusz, T. (2006). The igraph software package for complex network research. InterJournal, Complex Systems, 1695.


In conclusion, we have shown that the dissipative toric code with single-shot recovery jumps will host a qubit steady state structure in the presence of dephasing perturbations if the strength of the recovery process scales with the linear system size. Moreover, our analysis suggests that any local perturbation will not destroy the qubit steady state for such a system. As we discuss in the main text, the perturbative analysis also misses a non-perturbative contribution for the system's lifetime, which will result in a quantitatively different scaling dependence on $\Delta$.

[1]
B. M. Terhal, “Quantum error correction for quantum memories,” Rev. Mod.
Phys., vol. 87, pp. 307–346, 2015.

[2]
D. Gottesman, Stabilizer codes and quantum error correction.
PhD thesis, California Institute of Technology, 1997.

[3]
J. P. Paz and W. H. Zurek, “Continuous error correction,” Proc. R.
Soc. A: Math. Phys. Eng. Sci., vol. 454, no. 1969, pp. 355–364, 1998.

[4]
C. Ahn, A. C. Doherty, and A. J. Landahl, “Continuous quantum error
correction via quantum feedback control,” Phys. Rev. A, vol. 65,
no. 4, p. 042301, 2002.

[5]
M. Sarovar and G. J. Milburn, “Continuous quantum error correction by
cooling,” Phys. Rev. A, vol. 72, no. 1, p. 012306, 2005.

[6]
H. Mabuchi, “Continuous quantum error correction as classical hybrid
control,” New J. Phys., vol. 11, no. 10, p. 105044, 2009.

[7]
O. Oreshkov, “Continuous-time quantum error correction,” arXiv:1311.2485, 2013.

[8]
V. V. Albert and P. Faist, eds., The Error Correction Zoo.

[9]
Z. Leghtas, S. Touzard, I. M. Pop, A. Kou, B. Vlastakis, A. Petrenko, K. M.
Sliwa, A. Narla, S. Shankar, M. J. Hatridge, M. Reagor, L. Frunzio, R. J.
Schoelkopf, M. Mirrahimi, and M. H. Devoret, “Confining the state of light
to a quantum manifold by engineered two-photon loss,” Science,
vol. 347, no. 6224, pp. 853–857, 2015.

[10]
S. Touzard, A. Grimm, Z. Leghtas, S. O. Mundhada, P. Reinhold, C. Axline,
M. Reagor, K. Chou, J. Blumoff, K. M. Sliwa, S. Shankar, L. Frunzio, R. J.
Schoelkopf, M. Mirrahimi, and M. H. Devoret, “Coherent oscillations inside a
quantum manifold stabilized by dissipation,” Phys. Rev. X, vol. 8,
p. 021005, 2018.

[11]
J. M. Gertler, B. Baker, J. Li, S. Shirol, J. Koch, and C. Wang, “Protecting a
bosonic qubit with autonomous quantum error correction,” Nature,
vol. 590, no. 7845, pp. 243–248, 2021.

[12]
R. Lescanne, M. Villiers, T. Peronnin, A. Sarlette, M. Delbecq, B. Huard,
T. Kontos, M. Mirrahimi, and Z. Leghtas, “Exponential suppression of
bit-flips in a qubit encoded in an oscillator,” Nat. Phys., vol. 16,
no. 5, pp. 509–513, 2020.

[13]
B. de Neeve, T.-L. Nguyen, T. Behrle, and J. P. Home, “Error correction of a
logical grid state qubit by dissipative pumping,” Nat. Phys., vol. 18,
no. 3, pp. 296–300, 2022.

[14]
P. Campagne-Ibarcq, A. Eickbusch, S. Touzard, E. Zalys-Geller, N. E. Frattini,
V. V. Sivak, P. Reinhold, S. Puri, S. Shankar, R. J. Schoelkopf, L. Frunzio,
M. Mirrahimi, and M. H. Devoret, “Quantum error correction of a qubit
encoded in grid states of an oscillator,” Nature, vol. 584, no. 7821,
pp. 368–372, 2020.

[15]
V. V. Sivak, A. Eickbusch, B. Royer, S. Singh, I. Tsioutsios, S. Ganjam,
A. Miano, B. L. Brock, A. Z. Ding, L. Frunzio, S. M. Girvin, R. J.
Schoelkopf, and M. H. Devoret, “Real-time quantum error correction beyond
break-even,” Nature, vol. 616, no. 7955, pp. 50–55, 2023.

[16]
B. J. Brown, D. Loss, J. K. Pachos, C. N. Self, and J. R. Wootton, “Quantum
memories at finite temperature,” Rev. Mod. Phys., vol. 88, p. 045005,

[17]
R. Peierls, “On ising's model of ferromagnetism,” Math. Proc. Cambridge
Philos. Soc., vol. 32, no. 3, p. 477–481, 1936.

[18]
R. B. Griffiths, “Peierls proof of spontaneous magnetization in a
two-dimensional ising ferromagnet,” Phys. Rev., vol. 136,
pp. A437–A439, 1964.

[19]
E. Dennis, A. Kitaev, A. Landahl, and J. Preskill, “Topological quantum
memory,” J. Math. Phys., vol. 43, pp. 4452–4505, 2002.

[20]
R. Alicki, M. Horodecki, P. Horodecki, and R. Horodecki, “On thermal stability
of topological qubit in kitaev's 4d model,” Open Syst. Inf. Dyn.,
vol. 17, no. 01, pp. 1–20, 2010.

[21]
F. Pastawski, L. Clemente, and J. I. Cirac, “Quantum memories based on
engineered dissipation,” Phys. Rev. A, vol. 83, p. 012304, 2011.

[22]
Y.-J. Liu and S. Lieu, “Dissipative phase transitions and passive error
correction,” arXiv:2307.09512, 2023.

[23]
S. Bravyi and B. Terhal, “A no-go theorem for a two-dimensional
self-correcting quantum memory based on stabilizer codes,” New
J. Phys., vol. 11, no. 4, p. 043029, 2009.

[24]
O. Landon-Cardinal and D. Poulin, “Local topological order inhibits thermal
stability in 2d,” Phys. Rev. Lett., vol. 110, p. 090502, 2013.

[25]
J. Haah and J. Preskill, “Logical-operator tradeoff for local quantum codes,”
Phys. Rev. A, vol. 86, p. 032308, 2012.

[26]
M. B. Hastings, “Topological order at nonzero temperature,” Phys. Rev.
Lett., vol. 107, p. 210501, 2011.

[27]
B. Yoshida, “Feasibility of self-correcting quantum memory and thermal
stability of topological order,” Ann. Phys., vol. 326, no. 10,
pp. 2566–2633, 2011.

[28]
J. Haah, “Commuting pauli hamiltonians as maps between free modules,” Comm. Math. Phys., vol. 324, no. 2, pp. 351–399, 2013.

[29]
F. Pastawski and B. Yoshida, “Fault-tolerant logical gates in quantum
error-correcting codes,” Phys. Rev. A, vol. 91, p. 012305, 2015.

[30]
S. Hallgren, D. Nagaj, and S. Narayanaswami, “The local hamiltonian problem on
a line with eight states is qma-complete,” Quantum Info. Comput.,
vol. 13, no. 9–10, p. 721–750, 2013.

[31]
D. Aharonov, D. Gottesman, S. Irani, and J. Kempe, “The power of quantum
systems on a line,” Comm. Math. Phys., vol. 287, pp. 41–65, 4 2009.

[32]
N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, “Computational
complexity of projected entangled pair states,” Phys. Rev. Lett.,
vol. 98, p. 140506, 2007.

[33]
S. Piddock and A. Montanaro, “The complexity of antiferromagnetic interactions
and 2d lattices,” Quantum Info. Comput., vol. 17, no. 7–8,
p. 636–672, 2017.

[34]
J. Haah, “Local stabilizer codes in three dimensions without string logical
operators,” Phys. Rev. A, vol. 83, no. 4, p. 042330, 2011.

[35]
A. Hamma, C. Castelnovo, and C. Chamon, “Toric-boson model: Toward a
topological quantum memory at finite temperature,” Phys. Rev. B,
vol. 79, p. 245122, 2009.

[36]
J. R. Wootton, “Topological phases and self-correcting memories in interacting
anyon systems,” Phys. Rev. A, vol. 88, p. 062312, 2013.

[37]
S. Chesi, B. Röthlisberger, and D. Loss, “Self-correcting quantum memory in
a thermal environment,” Phys. Rev. A, vol. 82, p. 022305, 2010.

[38]
E. Kapit, J. T. Chalker, and S. H. Simon, “Passive correction of quantum
logical errors in a driven, dissipative system: A blueprint for an analog
quantum code fabric,” Phys. Rev. A, vol. 91, p. 062324, 2015.

[39]
S. Lieu, Y.-J. Liu, and A. V. Gorshkov, “Candidate for a passively-protected
quantum memory in two dimensions,” arXiv:2205.09767, 2022.

[40]
D. Gottesman, “Fault-Tolerant Quantum Computation with Constant Overhead,”
arXiv:1310.2984, 2013.

[41]
T.-C. Lin and M.-H. Hsieh, “Good quantum ldpc codes with linear time decoder
from lossless expanders,” arXiv:2203.03581, 2022.

[42]
A. Leverrier and G. Zémor, “Quantum tanner codes,” in 2022 IEEE 63rd
Annual Symposium on Foundations of Computer Science (FOCS), pp. 872–883,
IEEE, 2022.

[43]
P. Panteleev and G. Kalachev, “Asymptotically good quantum and locally
testable classical ldpc codes,” in Proceedings of the 54th Annual ACM
SIGACT Symposium on Theory of Computing, pp. 375–388, 2022.

[44]
R. Alicki, M. Fannes, and M. Horodecki, “On thermalization in kitaev's 2d
model,” J. Phys. A: Math. Theor., vol. 42, no. 6, p. 065303, 2009.

[45]
A. Lucia, D. Pérez-García, and A. Pérez-Hernández, “Thermalization in
kitaev's quantum double models via tensor network techniques,” arXiv:2107.01628, 2023.

[46]
I. Bardet, A. Capel, L. Gao, A. Lucia, D. Pérez-García, and C. Rouzé,
“Rapid thermalization of spin chain commuting hamiltonians,” Phys.
Rev. Lett., vol. 130, p. 060401, 2023.

[47]
H.-P. Breuer, F. Petruccione, et al., The theory of open quantum
Oxford University Press on Demand, 2002.

[48]
V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, “Completely positive
dynamical semigroups of n-level systems,” J. Math. Phys., vol. 17,
no. 5, pp. 821–825, 1976.

[49]
G. Lindblad, “On the generators of quantum dynamical semigroups,” Comm. Math. Phys., vol. 48, no. 2, pp. 119–130, 1976.

[50]
E. Knill and R. Laflamme, “Theory of quantum error-correcting codes,” Phys. Rev. A, vol. 55, pp. 900–911, 1997.

[51]
A. L. Grimsmo, J. Combes, and B. Q. Baragiola, “Quantum computing with
rotation-symmetric bosonic codes,” Phys. Rev. X, vol. 10, p. 011058,

[52]
V. V. Albert, B. Bradlyn, M. Fraas, and L. Jiang, “Geometry and response of
lindbladians,” Phys. Rev. X, vol. 6, p. 041031, 2016.

[53]
V. V. Albert, “Lindbladians with multiple steady states: theory and
applications,” arXiv:1802.00010, 2018.

[54]
“NIST Digital Library of Mathematical Functions.”
<https://dlmf.nist.gov/>, Release 1.1.9 of 2023-03-15.
F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider,
R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and
M. A. McClain, eds.

[55]
A. Kossakowski, “On quantum statistical mechanics of non-hamiltonian
systems,” Rep. Math. Phys., vol. 3, no. 4, pp. 247–274, 1972.

[56]
D. Gottesman, “Fault-tolerant quantum computation with higher-dimensional
systems,” Chaos Solit. Fractals, vol. 10, no. 10, pp. 1749–1758,

[57]
J. Bierbrauer and Y. Edel, “Quantum twisted codes,” J. Comb. Des.,
vol. 8, no. 3, pp. 174–188, 2000.

[58]
A. Ketkar, A. Klappenecker, S. Kumar, and P. K. Sarvepalli, “Nonbinary
stabilizer codes over finite fields,” IEEE transactions on information
theory, vol. 52, no. 11, pp. 4892–4914, 2006.

[59]
T. Cormen, C. Leiserson, R. Rivest, and C. Stein, Introduction To
Mit Electrical Engineering and Computer Science, MIT Press, 2001.

[60]
A. Kitaev, “Anyons in an exactly solved model and beyond,” Ann. Phys.,
vol. 321, no. 1, pp. 2–111, 2006.

[61]
T. M. Stace, S. D. Barrett, and A. C. Doherty, “Thresholds for topological
codes in the presence of loss,” Phys. Rev. Lett., vol. 102, p. 200501,

[62]
D. K. Tuckett, A. S. Darmawan, C. T. Chubb, S. Bravyi, S. D. Bartlett, and
S. T. Flammia, “Tailoring surface codes for highly biased noise,” Phys. Rev. X, vol. 9, p. 041031, 2019.

[63]
A. R. Calderbank, E. M. Rains, P. W. Shor, and N. J. A. Sloane, “Quantum error
correction and orthogonal geometry,” Phys. Rev. Lett., vol. 78,
pp. 405–408, 1997.

[64]
C. Wang, J. Harrington, and J. Preskill, “Confinement-higgs transition in a
disordered gauge theory and the accuracy threshold for quantum memory,” Ann. Phys., vol. 303, no. 1, pp. 31–58, 2003.

[65]
J. Dengis, R. König, and F. Pastawski, “An optimal dissipative encoder for
the toric code,” New J. Phys., vol. 16, no. 1, p. 013023, 2014.

[66]
F. Ticozzi, G. Baggio, and L. Viola, “Quantum Information Encoding from
Stabilizing Dynamics,” in 2019 IEEE 58th Conference on Decision and
Control (CDC), pp. 413–418, IEEE, 2019.

[67]
M. H. Michael, M. Silveri, R. T. Brierley, V. V. Albert, J. Salmilehto,
L. Jiang, and S. M. Girvin, “New class of quantum error-correcting codes for
a bosonic mode,” Phys. Rev. X, vol. 6, p. 031006, 2016.

[68]
A. Kubica and J. Preskill, “Cellular-automaton decoders with provable
thresholds for topological codes,” Phys. Rev. Lett., vol. 123,
p. 020501, 2019.

[69]
M. A. Nielsen and I. L. Chuang, Quantum computation and quantum
Cambridge University Press, 10th anniversary ed ed.

[70]
B. Buča and T. Prosen, “A note on symmetry reductions of the lindblad
equation: transport in constrained open spin chains,” New J. Phys.,
vol. 14, no. 7, p. 073007, 2012.


[1]C. D. Arrowsmith

1]Department of Physics, University of Oxford, Parks Road, Oxford, OX1 3PU, UK

2]European Organization for Nuclear Research (CERN), CH-1211 Geneva 23,
Switzerland

3]GSI Helmholtzzentrum für Schwerionenforschung GmbH, Planckstraße 1 64291
Darmstadt, Germany

4]GoLP/Instituto de Plasmas e Fusão Nuclear, Instituto Superior Técnico,
Universidade de Lisboa, 1049-001, Lisboa, Portugal

5]Lawrence Livermore National Laboratory, 7000 East Ave, Livermore,
California, 94550, USA

6]Rutherford Appleton Laboratory, Chilton, Didcot, OX11 0QX, UK

7]University of Rochester Laboratory for Laser Energetics, Rochester,
NY,14623, USA

8]Science Institute, University of Iceland, Dunhaga 3, IS-107, Reykjavik,
Iceland

9]School of Electrical Engineering and Computer Science, KTH Royal Institute
of Technology, SE-100 44, Stockholm, Sweden

10]AWE, Aldermaston, Reading, Berkshire, RG7 4PR, UK

11]Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, D-69117, Heidelberg,
Germany

12]School of Applied Mathematics and Physical Sciences, National Technical
University of Athens, Athens 157 72, Greece

13]Department of Physics, University of Strathclyde, Glasgow, G4 0NG, UK

# Laboratory realization of relativistic pair-plasma beams

<EMAIL_ADDRESS>P. Simon  P. Bilbao  A. F. A. Bott  S.
Burger  H. Chen  F. D. Cruz  T. Davenne  I. Efthymiopoulos  D. H. Froula  A.
M. Goillot  J. T. Gudmundsson  D. Haberberger  J. Halliday  T. Hodge  B. T.
Huffman  S. Iaquinta  F. Miniati  B. Reville  S. Sarkar  A. A. Schekochihin
L. O. Silva  R. Simpson  V. Stergiou  R. M. G. M. Trines  T. Vieu  N.
Charitonidis  R. Bingham  G. Gregori  [  [  [  [  [  [  [  [  [  [  [  [  [

(Dated: )

Relativistic electron-positron (e±) plasmas are ubiquitous in extreme
astrophysical environments such as black holes and neutron star
magnetospheres, where accretion-powered jets and pulsar winds are expected to
be enriched with such pair plasmas. Their behaviour is quite different from
typical electron-ion plasmas due to the matter-antimatter symmetry of the
charged components and their role in the dynamics of such compact objects is
believed to be fundamental. So far, our experimental inability to produce
large yields of positrons in quasi-neutral beams has restricted the
understanding of electron-positron pair plasmas to simple numerical and
analytical studies which are rather limited. We present first experimental
results confirming the generation of high-density, quasi-neutral, relativistic
electron-positron pair beams using the 440 GeV/c beam at CERN’s Super Proton
Synchrotron (SPS) accelerator. The produced pair beams have a volume that
fills multiple Debye spheres and are thus able to sustain collective plasma
oscillations. Our work opens up the possibility of directly probing the
microphysics of pair plasmas beyond quasi-linear evolution into regimes that
are challenging to simulate or measure via astronomical observations.

## Main

Relativistic electron-positron (e±) pair plasmas are expected to be produced
around black holes [1] and neutron stars [2]. In these environments, pair
creation can occur due to intense, high-energy $\gamma$-ray fluxes (by the
Breit-Wheeler process [3]) or when the electromagnetic fields are comparable
to the Schwinger field: the critical field strength for vacuum breakdown
($E_{\rm c}=$1.3\text{\times}{10}^{18}\text{\,}\mathrm{V}$$/m, $B_{\rm
c}=$4.4\text{\times}{10}^{9}\text{\,}\mathrm{T}$$) [4, 5]. Because of the
symmetry of the charged components, electron-positron pair plasmas should
exhibit collective behaviour that is significantly different from typical
electron-ion plasmas [6]. Linear and non-linear wave processes can be affected
in both fluid and kinetic regimes because of the suppression of some wave
modes. This is important in a variety of astrophysical settings, with recent
attention focusing on fast radio burst generation and the stability of
astrophysical pair beam jets [7, 8, 9, 10, 11]. However, producing
sufficiently large yields and densities of e± pairs in the laboratory in order
to directly probe the relevant plasma microphysics has been challenging.
Presently, high flux laboratory sources of positrons include: (i) nuclear
reactors [12], (ii) electron accelerators [13, 14], and (iii) high-power
lasers [15, 16, 17, 18, 19, 20, 21]. All these approaches involve pair
production processes when sufficiently energetic $\gamma$-rays
($E_{\gamma}\geq 2m_{\rm e}c^{2}=$1.022\text{\,}\mathrm{MeV}$$) interact with
charged nuclei (so-called Trident and Bethe-Heitler processes [22]), with the
highest cross-section in high-Z materials. In the coming decade, it is
proposed to use magnetic chicanes at FACET-II (SLAC) to combine the
accelerator’s e+ and e- beams into a quasi-neutral jet [23]. The next
generation of ultra-intense lasers may also be able to produce pairs by
achieving the Schwinger limit for vacuum breakdown [24, 25, 26]. Meanwhile,
precision magnetic confinement techniques have been developed to trap low-
temperature e± pair plasmas [27, 28, 29], and relativistic laser-produced
plasmas [30, 31, 32]. However, despite significant efforts, none of these
approaches have so far been able to produce the yields and densities of pairs
needed to sustain collective modes in the plasma.

Figure 1: Experimental setup. Protons with 440 GeV/c momentum are extracted
from the SPS ring with maximum fluence exceeding $3\times 10^{11}$ protons in
a single bunch of duration $250\text{\,}\mathrm{ps}$ (1-$\sigma$), and
transverse size $\sigma_{r}=$1\text{\,}\mathrm{mm}$$. The transverse beam
profile of the secondary beam is imaged using a $70$ mm $\times\,50$ mm
$\times\,0.25$ mm chromium-doped (Chromox) luminescence screen positioned 10
cm downstream of the target, and a blocker foil
($50\text{\,}\mathrm{\SIUnitSymbolMicro m}$ aluminium) is used to minimize
stray optical light. The Chromox plate is oriented at 45∘ to the beam path and
viewed by a digital camera which has an exposure time of 24 ms to capture the
entire scintillation of the screen. The $3.8\text{\,}\mathrm{m}$ standoff
distance of the digital camera leads to image resolution of
$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$, however the actual resolution is
$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ due to the translucence of the
Chromox. At a distance $2\text{\,}\mathrm{m}$ downstream of the target,
electrons and positrons are separated from the secondary beam and spectrally
resolved using a magnetic spectrometer comprised of an electromagnet and a
pair of luminescence screens ($200$ mm $\times\,50$ mm $\times\,1$ mm) centred
at a distance $240\text{\,}\mathrm{mm}$ off-axis. 20-cm thick bricks of
concrete (not shown in the diagram) are placed at the entrance of the
electromagnet, leaving a $40\text{\,}\mathrm{mm}$ wide aperture. Concrete is
also placed to block the target from the direct view of the cameras to
minimize speckle background arising on the camera images from the impact of
high energy hadrons scattered around the experimental area.

Here we present a novel approach for producing quasi-neutral e± jets in which
a high-intensity, ultra-relativistic proton beam is converted into pairs via
hadronic and electromagnetic cascades with 2-3 orders of magnitude higher
yield than previously reported neutral beams [17, 19]. We performed our
experiment at the HiRadMat (High-Radiation to Materials) facility [33] in the
accelerator complex at CERN, Geneva. The experimental setup is shown in Figure
1. We have performed detailed Monte-Carlo simulations using the standard
computer code FLUKA [34, 35, 36] to characterize the e± pair production, as
well as the other secondary beam components (hadrons and $\gamma$-rays). The
predicted number of pairs produced with kinetic energy greater than
$1\text{\,}\mathrm{MeV}$ is
$N_{\pm}=\frac{1}{2}(N_{\rm{e}^{+}}+N_{\rm{e}^{-}})=1.5\times 10^{13}$, with
peak pair density
$n_{\pm}=$1.6\text{\times}{10}^{12}\text{\,}{\mathrm{cm}}^{-3}$$, and the
ratio of positrons to electrons $N_{\rm{e}^{+}}/N_{\rm{e}^{-}}=0.82$.
Downstream of the target, the positron ratio is even higher
($N_{\rm{e}^{+}}/N_{\rm{e}^{-}}\gtrsim 0.9$), as low-energy electrons present
in the beam at the target rear surface preferentially escape the beam due to
their higher divergence.

The large numbers of electron-positron pairs are generated using a single LHC-
type bunch of $3\times 10^{11}$ protons with momentum $440~{}\mbox{GeV/c}$ and
duration (1-$\sigma$) of $\tau=$250\text{\,}\mathrm{ps}$$. The protons are
extracted to the facility from the Super Proton Synchrotron, irradiating a
solid-target composed of a low-Z material (graphite) and a high-Z converter
(tantalum). The dominant process for producing electron-positron pairs is
hadronization of quarks and gluons inside the graphite section of the target.
This produces a shower of pions, kaons and other hadrons on scales comparable
to the nuclear interaction length in graphite [37]. A copious number of ultra-
relativistic neutral pions are produced which almost instantaneously undergo
electromagnetic decay to produce a highly-collimated flux of GeV-scale
$\gamma$-rays. Electromagnetic cascades are then generated with the
$\gamma$-rays producing pairs in the high-Z tantalum converter which is much
longer than the conversion length. Further pairs are created via subsequent
bremsstrahlung of electrons and positrons (Bethe-Heitler process). Secondary
$\gamma$-rays which do not convert into pairs can escape the target, along
with a much smaller (by orders of magnitude) number of protons and other
hadronic species. The effect of hadronic beam components (such as $\pi^{\pm}$
pairs) on e± pair plasma dynamics must be considered, but the effects are
expected to be negligible due to the much lower mobility and density of these
species. The choice of target material length constitutes a compromise between
the number of pairs produced and the emittance of the e± beam, with the
thickness of graphite and tantalum chosen in the current setup to maximize the
pair density ($n_{\pm}$) maintained over a $1\text{\,}\mathrm{m}$ length
downstream of the target. A choice of thicker target materials can produce an
even greater pair density and number of electron-positron pairs at the
immediate rear of the target [37].

Measurements have been performed to validate the predicted beam
characteristics by recording the transverse beam profile and particle energy
spectra (see Figure 1 for details). This is achieved using chromium-doped
alumina-ceramic (Chromox) luminescence screens [38, 39, 40] placed in the beam
path downstream of the target. When ionizing particles or radiation are
incident on the screen, it emits red visible light with a few ms decay time.
For our conditions, the intensity of light emitted from the screen is directly
proportional to the energy deposited by ionizing particles (see Methods).
Although it is not possible to distinguish which secondary particle species is
the cause of observed luminescence light, above relativistic energies (Lorentz
factor, $\Gamma\gtrsim 2$) charged particles deposit energy as ‘minimum-
ionizing’, a characteristic which results in relativistic charged particles
depositing an almost identical amount of energy into the screen (see
Supplementary Information). Given that the vast majority of secondary
particles incident on the screen are relativistic, the observed brightness is
thus proportional to the number of incident particles. The expectation from
FLUKA simulations is that the production of $10^{13}$ pairs will lead to about
100 times larger luminescence compared with a primary beam containing $\sim
10^{11}$ protons. Since the target is mounted onto a motorized stage, it can
be entirely removed from the proton beam path. By independently measuring the
incident proton beam fluence for each shot using upstream current monitors, an
absolute calibration of the charge incident on the screen can be made.

A common source of background in experiments producing energetic
electromagnetic cascades is the large number of scattered sub-MeV e- and
$\gamma$-rays which can flood the detectors. Given that our experiment is
carried out in an air environment, a large fraction of low-energy particles
and radiation are absorbed before they can reach the screens; for instance, e±
with energy $\lesssim$100\text{\,}\mathrm{keV}$$ are mostly absorbed by a few
centimetres of air. The air environment can provide additional sources of
background in the form of stray light arising from Cherenkov emission,
fluorescence of air molecules, and optical transition radiation (OTR) which is
produced by particles passing from different dielectric materials into the air
(such as the target-air interface). While the contribution of all of these
sources is small compared to the light collected due to luminescence, an
aluminium blocker foil is placed in the beam path before the luminescence
screen to reduce on-axis Cherenkov and OTR from illuminating the Chromox
plate.

Figure 2: Transverse beam profile imaged using a luminescence screen. (a)
Direct comparison of FLUKA Monte-Carlo simulations with raw image data
obtained when the target is irradiated and the secondary beam is produced
(‘Target in’), versus when the target is removed and only the primary proton
beam irradiates the screen (‘No target’). An absolute fluence calibration is
obtained using the known density profile of the primary proton beam. (b)
Integrated image intensity (total intensity) from 70 shots is converted to an
absolute particle number, showing the case where the target is irradiated (red
circles, 47 shots), and when it is removed (black diamonds, 23 shots). The
error bars reflect the standard errors of the fitted parameters for each shot.
FLUKA Monte-Carlo simulations of the predicted light yield are shown for both
cases (black-dashed and red-dashed lines), showing strong agreement with the
experimental data. The blue dot-dashed line shows the contribution from e± in
the FLUKA simulation, highlighting that this is the dominant contribution to
the enhanced signal.

The experimental results of the post-target in-beam luminescence screen are
summarized in Figure 2, comparing directly with FLUKA Monte-Carlo simulations.
Figure 2a shows the raw image data of the transverse beam profile when the
target is irradiated compared with when the target is removed from the proton
path. The image intensity is converted to an absolute particle fluence by
normalizing to the known beam density profile of the primary proton beam
(intensity $3\times 10^{11}$ protons, Gaussian width
$\sigma_{r}=$1\text{\,}\mathrm{mm}$$, temporal duration
$\tau=$250\text{\,}\mathrm{ps}$$). When the target is irradiated and the
secondary beam is produced, a 5 times increase in peak brightness is observed.
In addition, the transverse size of the beam broadens to a Lorentzian profile
with half-width $\Sigma_{r}=2.3$ mm due to its increased divergence from
Coulomb scattering of pairs with atomic nuclei in the target material.

The integrated image intensities (total intensities) obtained from 70 shots
with the target irradiated (red circles) and the target removed (black
diamonds) are shown in Figure 2b, again converted to an absolute particle
number. The total intensity scales with the number of protons in the primary
beam. Given that $90$% of the measured signal corresponds to e- and e+, the
number of e± captured on the screen is measured to be
$N_{\pm,\rm{exp}}=(1.02\pm 0.05)\times 10^{13}$. This measurement agrees with
the FLUKA expectation for a screen placed the same distance downstream of the
target, $N_{\pm,\rm{FLUKA}}=1.04\times 10^{13}$. We estimate the peak pair
density at the screen position to be
$n_{\pm,\rm{exp}}=$9\text{\times}{10}^{11}\text{\,}{\mathrm{cm}}^{-3}$$,
taking the longitudinal beam profile as identical to the primary beam (beam
elongation due to straggling in the target is calculated to be $\lesssim
2\%$).

Figure 3: Magnetic particle spectrometer. The energy spectra of e- (blue) and
e+ (red) are obtained from images of luminescence screens (dimensions $200$ mm
$\times\,50$ mm $\times\,1$ mm) centred $240\text{\,}\mathrm{mm}$ from the
beam axis on either side (see Figure 1). Electrons and positrons are deflected
onto the screens by the vertically-oriented dipole magnetic field of the
electromagnet, whilst hadrons with a higher momentum and uncharged
$\gamma$-rays are mostly absorbed by the beam dump behind the electromagnet.
The spectrum in the energy range $30\leq E~{}\mbox{[MeV]}\leq 220$ is
constructed by piecing together images from multiple shots using different
magnetic field strengths ($B=0.1-$$0.34\text{\,}\mathrm{T}$). The shaded
regions correspond to the error associated with the absolute calibration.
FLUKA simulations (histograms) are able to accurately predict the
experimentally obtained spectra.

An electromagnet is used to measure the particle energy spectra of the
electrons and positrons separately from the other secondaries. The current in
the magnet coils can be varied for different shots to sample different
portions of the e- and e+ energy spectra between $30-220$ MeV with an energy
resolution $10-20$%. In addition, the electromagnet can be switched off,
making it possible to characterize the $\gamma$-ray and hadron background.
FLUKA simulations predict that without the magnetic field, low energy
($\lesssim 10$ MeV) scattered e± can irradiate the screens, but when the
electromagnet is activated, these lower-momentum pairs are swept away and the
measured signals more closely resemble the true spectra. This is confirmed by
the overlapping of the segments of the spectrum measured in different energy
ranges when the full spectrum is pieced together (shown in Figure 3). At
positions on the luminescence screens sufficiently far from the beam axis
($\gtrsim 200$ mm) the number of particles irradiating the screens is enhanced
compared with when the electromagnet is deactivated (raw data is shown in the
Supplementary Information). This observation can only be explained by the
presence of large numbers of e± with momentum $\leq 220$ MeV/c deflected onto
the screen. Higher momentum hadrons or $\gamma$-rays are not deflected onto
the screens and are instead absorbed by the beam dump downstream of the
electromagnet. An absolute calibration of number of electrons and positrons
measured with the magnetic spectrometer is obtained by direct comparison with
the amount of light collected from the post-target luminescence screens,
taking the different optical setups into account (the shaded regions in Figure
3 account for the errors associated with this absolute calibration). The
results show that FLUKA simulations are able to accurately predict both the
absolute number and the spectral shape of the e± pairs that can reach the
screens. In this energy range, the number of pairs is measured to be
$N_{\pm,\rm{exp}}=(2.46\pm 0.62)\times 10^{11}$ and positron fraction is
$(N_{\rm{e}^{+}}/N_{\rm{e}^{-}})_{\rm{exp}}=0.92\pm 0.05$, compared with the
FLUKA simulation: $N_{\pm,\rm{FLUKA}}=2.45\times 10^{11}$ and
$(N_{\rm{e}^{+}}/N_{\rm{e}^{-}})_{\rm{FLUKA}}=0.89$. The results of both of
the luminescence screen diagnostics are summarized in Table 1.

Given the experimental validation of the pair beam characteristics predicted
by FLUKA simulations, we use the simulations to obtain the beam properties at
the rear surface of the target. At this position the pair density is at its
highest ($n_{\pm}=$1.64\text{\times}{10}^{12}\text{\,}{\mathrm{cm}}^{-3}$$).
We assess whether this density is sufficiently high for the pair beam to be
considered a plasma. A fundamental characteristic of a plasma is its ability
to screen electric potentials that are applied to it. This is the case if the
physical size of the plasma exceeds the characteristic Debye screening length,
$r_{\rm{scr}}$, or equivalently, if the number of particles per Debye
screening volume is much greater than one, i.e. $N_{\pm}/N_{\rm{D}}\gg 1$,
where $N_{\rm{D}}=n_{\pm}r_{\rm scr}^{3}$.

The pair beams produced are highly relativistic ($\langle\Gamma\rangle\gg 1$),
possessing a relativistic thermal spread ($\Gamma_{T}=k_{\rm B}T_{\pm}/m_{\rm
e}c^{2}\gg 1$). Therefore the appropriate screening length is that of a
relativistically hot pair plasma, defined by $r_{\rm scr}^{-2}=8\pi
n_{\pm}e^{2}/k_{\rm{B}}T_{\pm}$. This is derived assuming the electrons and
positrons have a relativistic Maxwellian (Jüttner-Synge) thermal distribution
with pair temperature $T_{\pm}$ (see Methods). In the derivation, the plasma
has no bulk flow so we evaluate $N_{\pm}/N_{\rm{D}}$ in an inertial frame
where the beam has zero net momentum (the ‘zero-momentum’ frame). The pair
temperature is obtained in the zero-momentum frame by performing a Lorentz
transformation ($\Gamma_{0}\sim 4.5$) on the FLUKA-simulated momentum
distributions, leading to an approximately-isotropic momentum distribution
which is fit by $T_{\pm}=4~{}\mbox{MeV}$. The pair density in the zero-
momentum frame is given by $n_{\pm}/\Gamma_{0}$, where $n_{\pm}$ is the pair
density in the laboratory frame. Using these parameters, the number of
particles in a Debye screening volume is $N_{\pm}/N_{\rm{D}}\sim 8$, with
scaling

$\frac{N_{\pm}}{N_{\rm D}}=8~{}\bigg{(}\frac{N_{\pm}}{1.5\times
10^{13}}\bigg{)}\bigg{(}\frac{n_{\pm}}{$1.6\text{\times}{10}^{12}\text{\,}{\mathrm{cm}}^{-3}$}\bigg{)}^{1/2}\bigg{(}\frac{\Gamma_{0}}{4.5}\bigg{)}^{-1/2}\bigg{(}\frac{T_{\pm}}{$4\text{\,}\mathrm{MeV}$}\bigg{)}^{-3/2}.$
(1)

The pair beams can be used to study plasma processes relevant to the
interpenetration of pair jets within an ambient electron-ion plasma. In this
case, the condition to observe collective plasma phenomena is less strict. The
beam size must exceed the characteristic growth length of plasma
instabilities, which is typically the plasma skin depth of the ambient plasma,
$\lambda_{\rm s}^{2}=m_{\rm{e}}c^{2}/4\pi n_{\rm{p}}e^{2}$. Comparing the pair
number with $N_{\rm s}=n_{\pm}\lambda_{\rm s}^{3}$, assuming equal density for
the pair beam and the ambient plasma, we obtain $N_{\pm}/N_{\rm s}\sim 132$.

In Figure 4 we plot the number of pairs and beam density achieved in our
experiment alongside results of high densities of pairs reported from previous
experiments. High pair densities and yields can be achieved in laser
experiments, but at the expense of beam dimensions or charge neutrality.

By these definitions, the electron-positron beam we produced at HiRadMat can
safely be assumed to behave as a relativistic pair plasma, making it possible
to study collective plasma phenomena and providing a laboratory benchmark of
several high-energy astrophysics processes.

Position | Pair yield, $N_{\pm}$ | Positron fraction, ($N_{\rm{e}^{+}}/N_{\rm{e}^{-}}$)
---|---|---
Simulation | Experiment | Simulation | Experiment
Target rear surface ($E\geq 1$ MeV) | $1.53\times 10^{13}$ | - | $0.82$ | -
Post-target screen ($E\geq 1$ MeV) | $1.04\times 10^{13}$ | $(1.02\pm 0.05)\times 10^{13}$ | $0.91$ | -
Spectrometer ($30\leq E\mbox{ [MeV]}\leq 220$) | $2.45\times 10^{11}$ | $(2.46\pm 0.62)\times 10^{11}$ | $0.89$ | $0.92\pm 0.05$

Table 1: Summary of the measured and simulated yields of electron-positron
pairs. The experimentally measured and FLUKA-simulated electron-positron pair
yield ($N_{\pm}$) and positron fraction ($N_{\rm{e}^{+}}/N_{\rm{e}^{-}}$) are
summarized for the following positions: (i) at the rear surface of the target,
(ii) on the post-target luminescence screen, and (iii) on the particle
spectrometer luminescence screens. Figure 4: Comparison of laboratory-produced
high-density pair beams. The peak number and density of pairs reported in this
study (red square) compared with previous experiments performed at high-power
laser facilities (black squares): Orion/OMEGA-EP [15], Texas-PW [16], Astra-
Gemini [17], OMEGA-EP(’21) [19], OMEGA-EP(’14) [30]. The data are labelled by
the facility where the experiment was performed, and the fill fraction of each
marker corresponds to the fraction of positrons to electrons in the
experiment, $N_{\rm{e}^{+}}/N_{\rm{e}^{-}}$ ($\blacksquare=100\%$,
$\square=0\%$, also see the key). The blue-shaded region corresponds to when
the number of particles per Debye sphere is fewer than one,
$N_{\pm}/N_{\rm{D}}<1$ (assuming the screening length used in this work).

## Methods

### Screening lengths of relativistic plasmas

For a relativistic plasma the screening length is defined through the static
limit of the dielectric tensor $\varepsilon_{ij}$ [41, 42, 43]:

$\lim_{k=0}\,\lim_{\omega/k=0}k^{2}(\varepsilon^{l}-1)=r_{\rm scr}^{-2},$ (2)

where the longitudinal component of the dielectric tensor is

$\varepsilon^{l}=1+\sum_{\alpha}\frac{4\pi
q_{\alpha}^{2}}{k^{2}}\lim_{\eta\rightarrow+0}\int\mathbf{k}\cdot\frac{\partial
f_{\alpha}(\mathbf{p})}{\partial\mathbf{p}}\frac{d^{3}p}{\omega+i\eta-\mathbf{k}\cdot\mathbf{v}},$
(3)

and $\alpha$ labels each species component of the plasma. Assuming a pair
plasma comprised of electrons and positrons which have a relativistic
Maxwellian (Jüttner-Synge) distribution:

$f(\mathbf{p})=\frac{n_{e}}{8\pi\Gamma_{T}^{3}}\exp\left(-\frac{\Gamma}{\Gamma_{T}}\right),$
(4)

with ultra-relativistic temperature
$\Gamma_{T}=k_{\rm{B}}T_{\pm}/m_{\rm{e}}c^{2}\gg 1$ (noting that this tends to
a Maxwellian in the limit $\Gamma\ll 1$, $\Gamma_{T}\ll 1$). The screening
length can be evaluated in the relativistic limit ($v=c$) as

$r_{\rm scr}^{-2}=\frac{8\pi n_{\pm}e^{2}}{k_{\rm{B}}T_{\pm}},$ (5)

where $n_{\pm}$ and $T_{\pm}$ refer to the density and temperature of pairs.
Identical treatment in the classical limit for an electron-ion plasma leads to
the well-known Debye screening, $r_{\rm scr,D}^{-2}=k_{\rm{D}}^{2}=4\pi
n_{\rm{e}}e^{2}/k_{\rm{B}}T_{\rm{e}}$. Since this is derived for a plasma
which has net-zero momentum flow, we obtain the plasma density and temperature
by fitting Jüttner-Synge distributions in a Lorentz-boosted inertial frame
(see the Supplementary Information for details of the fitting).

### Electron-positron pair production target

The target is designed such that a quasi-neutral e± beam is produced with pair
beam density maximized along a 1-m length downstream of the target. The target
can be irradiated by many (potentially hundreds or thousands) of single
bunches containing $O(10^{11})$ primary protons, without compromising its
structural integrity. In addition, the target must cool sufficiently fast to
return to room temperature in-between shots, at a maximum repetition rate of 1
bunch per minute. The target design has been optimized using simulations with
two codes: (i) FLUKA [34, 35] (along with the associated interface “Flair”
[36]), a particle transport Monte-Carlo scattering code capable of accurately
describing the hadronic and electromagnetic cascades when the target is
irradiated with 440 GeV/c protons, and (ii) Ansys${}^{\mbox{\textregistered}}$
Mechanical [44], a finite-element code used to estimate the target cooling and
the stress/strain induced by the energy deposition of the beam. A FLUKA-
simulated transient thermal field is used as the initial condition.

FLUKA uses a robustly bench-marked physics model. To achieve a good
statistical representation in the Monte Carlo method, $10^{5}$ iterations were
performed. The low-energy cutoff for particle transport in the simulation is
$10$ keV for $\rm{e}^{-}/\rm{e}^{+}/\gamma$ and $100$ keV for hadrons.

The chosen target design consisted of a 360 mm long cylinder of isostatic
graphite (SGL Carbon R6650,
$1.84\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$) and 10-mm thick disk of
tantalum, both having a $20\text{\,}\mathrm{mm}$ diameter. The graphite and
tantalum are housed inside a 400-mm length, 50-mm diameter cylinder of high-
strength T9 aluminium alloy that acts as both a confinement vessel and a heat
sink. The tantalum is press-fit to ensure maximal thermal contact. 2-mm thick
expanded graphite pieces (SGL Carbon Sigraflex,
$1\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$) separate the target
components to allow thermal expansion and reduce contact stresses during
irradiation, while 2-mm thickness Sigradur G glassy carbon beam windows are
clamped onto either end of the target by aluminium flanges with Viton O-rings
to hermetically seal the target materials. Using this design, the Ansys
simulations have shown that the thermal loading per $3\times 10^{11}$ protons
is highest inside the tantalum, reaching peak instantaneous temperatures of
$300\text{\,}\mathrm{\SIUnitSymbolCelsius}$. Radiative and convective cooling
via the outer surface of the target housing leads to cooling of the target to
room temperature within a few seconds following the beam impact, while the
beam-induced maximum strain of the tantalum remains in all cases well below
its plastic deformation limit.

### Chromium-doped luminescence screens

Chromium-doped alumina-ceramic luminescence screens (Chromox, Al2O3: 99.5%,
Cr2O3: 0.5% [38]) have been used to measure the particle beam intensity and
transverse profile during the experiment. In the Chromox screens, principal
luminescence is due to de-excitations of the lowest-excited state of Cr3+ when
energy is deposited in the screen by ionizing particles and radiation. Light
is emitted isotropically, strongest at wavelengths $\lambda_{1}=691$ nm, and
$\lambda_{2}=694$ nm with decay times $3-6$ ms.

The transverse beam profile of the secondary beam is imaged using a $70$ mm
$\times\,50$ mm $\times\,0.25$ mm screen positioned 10 cm downstream of the
target, and a blocker foil ($50\text{\,}\mathrm{\SIUnitSymbolMicro m}$
aluminium) is used to minimize stray optical light. The screen is oriented at
45∘ to the beam path and viewed directly by a digital camera (Basler
acA1920-40gm GigE camera with Sony IMX249 CMOS sensor and Canon EF 75-300 mm
f/4-5.6 III lens) at a standoff distance 3.8 m with an exposure time 24 ms. An
almost identical optical setup is used to image the screens in the magnetic
spectrometer, except larger screens are used ($200$ mm $\times\,50$ mm
$\times\,1$ mm, centred at a distance $240$ mm off-axis), and viewed through a
single mirror reflection at a standoff distance 6.2 m.

Given that relativistic particles in the energy range of interest exhibit
minimum-ionizing behaviour, the energy deposition of a particle passing
through the Chromox screen is expected to be approximately insensitive to
energy and constant between singly-charged particle species (a result
confirmed by FLUKA simulations, see the Supplementary Information). The
translucence of the Chromox screens to the luminescence light (attenuation
length, $\mu=0.8$ mm-1) limits the spatial resolution to
$\gtrsim$$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$, as the luminescence
light is not significantly attenuated as it is transmitted from a region where
energy is deposited deeper into the screen. The translucence of the screens
simplifies the analysis; as a first approximation we don’t consider the
different longitudinal energy deposition profiles, which are anyway shown in
simulations to be approximately uniform through the screen thickness for the
relativistic particles observed.

### Magnetic electron-positron spectrometer

Before the experiment, we have characterized the spatial magnetic field
profile of the electromagnet with currents supplied to the coils in the range
between $0-400$ A. This electromagnet, designated “MNPA” in the CERN internal
naming system, has a yoke length of $250$ mm (total length $544$ mm), an
aperture of width $260$ mm and height $202$ mm, and a maximum field at 400A of
0.34 T. The exact magnet geometry has been modelled in detail using the
finite-element code Opera 3D [45]. The field map has been calculated with
sufficient resolution (10 mm) to capture the magnetic field gradients inside
the magnet gap, given that magnetic field transitions from 10%-90% peak field
strength over a distance $200$ mm. The field map has been cross-checked
against direct measurements of the magnet field using a Hall probe and the
difference between the model and measured magnetic fields inside the magnet
gap are $\lesssim 2$%. Finally, an energy calibration for the spectrometer
(correlating screen position with e± energy) is obtained from particle ray-
tracing calculations using the magnetic field maps. The e± energy ranges
sampled by different magnet settings are $30-55$ MeV, $55-80$ MeV, $75-110$
MeV, and $90-220$ MeV.

An absolute calibration of the electron and positron numbers is made by using
the brightness (pixel counts per unit area) of the in-beam luminescence
screen. Specifically, we account for the difference between the amount of
light collected in the optical setups used for the in-beam screen luminescence
and for the spectrometer screens, considering (i) the different standoff
distances, which leads to different solid angle subtended ($d\Omega_{\rm
beam}/d\Omega_{\rm spect}=2.7\pm 0.4$); (ii) the different camera gain
settings used ($G_{\rm spect}=38\pm 2$, $G_{\rm beam}=1$); and (iii) the
different thickness of screen used, where thicker screens lead to larger
energy deposition per particle ($u_{\rm dep,spect}/u_{\rm dep,beam}=3.5\pm
0.5$). A small $1/\cos\theta$ geometric correction is applied to the energy
spectra to account for the additional path length of Chromox encountered by
obliquely-incident deflected particles, where $\theta=10^{\circ}$ at the
screen edge closest to the beam axis, and $\theta=25^{\circ}$ at the furthest
edge.

The digital cameras viewing the in-beam screen at a standoff distance of
$3.8\text{\,}\mathrm{m}$ can resolve features as small as
$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$ in size, whilst the cameras viewing
the spectrometer screens at a standoff distance $6.2\text{\,}\mathrm{m}$ can
resolve features $120\text{\,}\mathrm{\SIUnitSymbolMicro m}$ in size. However,
the resolution of the energy spectrum projected onto the spectrometer screens
is limited by the 20-cm thickness, 20-mm wide concrete aperture at the
entrance of the electromagnet.

In our experimental setup, the target has been placed on a vertically movable
high-precision stage, allowing us to acquire data with the target in-beam, as
well as in a ‘target-out’ position. In the latter case, the primary proton
beam continues at its full intensity through the luminescence screens and the
electromagnet towards the beam dump. Using the ‘target-out’ configuration, we
took measurements without the current supplied to the electromagnet to
characterize the hadron and lepton background produced as particles are back-
scattered by the proton beam impact on the beam dump. The remnant field of the
aforementioned electromagnet was measured extensively before the experiment
using a Hall-probe and was found reproducibly to be negligible (on the order
of the noise of the instrument, i.e. $B\lesssim$0.3\text{\,}\mathrm{mT}$$).

## References

  *   * [1] Bambi, C. (ed.) _Astrophysics of Black Holes: From Fundamental Aspects to Latest Developments_ (Springer, Berlin/Heidelberg, 2016).
  * [2] Arons, J.  Some problems of pulsar physics or I’m madly in love with electricity.  _Space Science Reviews_ 24, 437–510 (1979).
  * [3] Breit, G. & Wheeler, J. A.  Collision of Two Light Quanta.  _Physical Review_ 46, 1087 (1934).
  * [4] Schwinger, J.  On Gauge Invariance and Vacuum Polarization.  _Physical Review_ 82, 664 (1951).
  * [5] Erber, T.  High-Energy Electromagnetic Conversion Processes in Intense Magnetic Fields.  _Reviews of Modern Physics_ 38, 626 (1966).
  * [6] Tsytovich, V. & Wharton, C. B.  Laboratory electron-positron plasma – a new research object.  _Comments on Plasma Physics and Controlled Fusion_ 4, 91–100 (1978).
  * [7] Arons, J.  Pair creation above pulsar polar caps: Geometrical structure and energetics of slot gaps.  _The Astrophysical Journal_ 266, 215–241 (1983).
  * [8] Begelman, M. C., Blandford, R. D. & Rees, M. J.  Theory of extragalactic radio sources.  _Reviews of Modern Physics_ 56, 255 (1984).
  * [9] Blandford, R. D. & Levinson, A.  Pair cascades in extragalactic jets. I: Gamma rays.  _The Astrophysical Journal_ 441, 79–95 (1995).
  * [10] Turolla, R., Zane, S. & Watts, A. L.  Magnetars: the physics behind observations. A review.  _Reports on Progress in Physics_ 78, 116901 (2015).
  * [11] Lyubarsky, Y.  Emission Mechanisms of Fast Radio Bursts.  _Universe_ 7, 56 (2021).
  * [12] Hugenschmidt, C., Piochacz, C., Reiner, M. & Schreckenbach, K.  The NEPOMUC upgrade and advanced positron beam experiments.  _New Journal of Physics_ 14, 055027 (2012).
  * [13] Bernardini, C.  AdA: The first electron-positron collider.  _Physics in Perspective_ 6, 156–183 (2004).
  * [14] Blumer, P. _et al._ Positron accumulation in the GBAR experiment.  _Nuclear Instruments and Methods in Physics Research Section A_ 1040, 167263 (2022).
  * [15] Chen, H. _et al._ Scaling the Yield of Laser-Driven Electron-Positron Jets to Laboratory Astrophysical Applications.  _Physical Review Letters_ 114, 215001 (2015).
  * [16] Liang, E. _et al._ High e+/e– Ratio Dense Pair Creation with 1021 W cm-2 Laser Irradiating Solid Targets.  _Scientific Reports_ 5, 13968 (2015).
  * [17] Sarri, G. _et al._ Generation of neutral and high-density electron–positron pair plasmas in the laboratory.  _Nature Communications_ 6, 6747 (2015).
  * [18] Xu, T. _et al._ Ultrashort megaelectronvolt positron beam generation based on laser-accelerated electrons.  _Physics of Plasmas_ 23, 033109 (2016).
  * [19] Peebles, J. L. _et al._ Magnetically collimated relativistic charge-neutral electron–positron beams from high-power lasers.  _Physics of Plasmas_ 28, 074501 (2021).
  * [20] Jiang, S. _et al._ Enhancing positron production using front surface target structures.  _Applied Physics Letters_ 118, 094101 (2021).
  * [21] Chen, H. & Fiuza, F.  Perspectives on relativistic electron–positron pair plasma experiments of astrophysical relevance using high-power lasers.  _Physics of Plasmas_ 30, 020601 (2023).
  * [22] Bethe, H. & Heitler, W.  On the Stopping of Fast Particles and on the Creation of Positive Electrons.  _Proceedings of the Royal Society of London. Series A_ 146, 83–112 (1934).
  * [23] Yakimenko, V. _et al._ FACET-II facility for advanced accelerator experimental tests.  _Physical Review Accelerators and Beams_ 22, 101301 (2019).
  * [24] Bell, A. R. & Kirk, J. G.  Possibility of Prolific Pair Production with High-Power Lasers.  _Physical Review Letters_ 101, 200403 (2008).
  * [25] Ridgers, C. P. _et al._ Dense Electron-Positron Plasmas and Ultraintense $\gamma$ rays from Laser-Irradiated Solids.  _Physical Review Letters_ 108, 165006 (2012).
  * [26] Zhang, P., Bulanov, S. S., Seipt, D., Arefiev, A. V. & Thomas, A. G. R.  Relativistic plasma physics in supercritical fields.  _Physics of Plasmas_ 27, 050601 (2020).
  * [27] Greaves, R. G. & Surko, C. M.  An Electron-Positron Beam-Plasma Experiment.  _Physical Review Letters_ 75, 3846 (1995).
  * [28] Danielson, J. R., Dubin, D. H. E., Greaves, R. G. & Surko, C. M.  Plasma and trap-based techniques for science with positrons.  _Reviews of Modern Physics_ 87, 247 (2015).
  * [29] Stenson, E. V. _et al._ Lossless Positron Injection into a Magnetic Dipole Trap.  _Physical Review Letters_ 121, 235005 (2018).
  * [30] Chen, H. _et al._ Magnetic collimation of relativistic positrons and electrons from high intensity laser–matter interactions.  _Physics of Plasmas_ 21, 040703 (2014).
  * [31] von der Linden, J. _et al._ Confinement of relativistic electrons in a magnetic mirror en-route to a magnetized relativistic pair plasma.  _Physics of Plasmas_ 28, 092508 (2021).
  * [32] Stoneking, M. R. _et al._ A new frontier in laboratory physics: magnetized electron–positron plasmas.  _Journal of Plasma Physics_ 86, 155860601 (2020).
  * [33] Efthymiopoulos, I. _et al._ HiRadMat: A New Irradiation Facility for Material Testing at CERN, Technical Report CERN-ATS-2011-232, CERN (2011).
  * [34] Ahdida, C. _et al._ New Capabilities of the FLUKA Multi-Purpose Code.  _Frontiers in Physics_ 9, 788253 (2022).
  * [35] Battistoni, G. _et al._ Overview of the FLUKA code.  _Annals of Nuclear Energy_ 82, 10–18 (2015).
  * [36] Vlachoudis, V. _et al._ FLAIR: A Powerful but User Friendly Graphical Interface for FLUKA.  _Proc. Int. Conf. on Mathematics, Computational Methods & Reactor Physics (M&C 2009), Saratoga Springs, New York_ (2009).
  * [37] Arrowsmith, C. D. _et al._ Generating ultradense pair beams using 400 GeV/c protons.  _Physical Review Research_ 3, 023103 (2021).
  * [38] McCarthy, K. J. _et al._ Characterization of the response of chromium-doped alumina screens in the vacuum ultraviolet using synchrotron radiation.  _Journal of Applied Physics_ 92, 6541–6545 (2002).
  * [39] Burger, S., Biskup, B., Mazzoni, S., Turner, M. _et al._ Scintillation and OTR Screen characterization with a 440 GeV/c Proton Beam in air at the CERN HiRadMat Facility.  _Proceedings of IBIC 2016_ (2016).
  * [40] Gorgisyan, I. _et al._ Commissioning of beam instrumentation at the CERN AWAKE facility after integration of the electron beam line.  _Journal of Physics: Conference Series_ 1067, 072015 (2018).
  * [41] Lifshitz, E. M. & Pitaevskii, L. P.  _Physical Kinetics (Course of Theoretical Physics: Vol. 10)_ (Pergamon, Oxford, 1981).
  * [42] Silin, V. P.  On the electromagnetic properties of a relativistic plasma.  _Sov. Phys. JETP_ 11, 1136–1140 (1960).
  * [43] Braaten, E. & Segel, D.  Neutrino energy loss from the plasma process at all temperatures and densities.  _Physical Review D_ 48, 1478 (1993).
  * [44] _Ansys® Academic Research Mechanical, Release 18.1, Help System, Coupled Field Analysis Guide, ANSYS, Inc._
  * [45] _Opera 3D, Dassault Systèmes®_.

Acknowledgments We thank Prof. C. Joshi (UCLA), Dr. F. Albert (LLNL) and Dr.
C. Densham (STFC Rutherford Appleton Laboratory) for useful discussions, and
Dr. T. Ma (LLNL) for supporting this experiment. This project has received
funding from the European Union’s Horizon Europe Research and Innovation
programme under Grant Agreement No 101057511 (EURO-LABS). The work of G.G. was
partially supported by UKRI under grants no. ST/W000903/1 and EP/Y035038/1,
while A.F.A.B. was also supported by UKRI (grant number MR/W006723/1). The
work of D.H.F. and D.H. was supported by the U.S. Department of Energy under
Award Number DE-NA0004144. We also acknowledge funding from AWE plc., and the
Central Laser Facility (STFC). FLUKA simulations were performed using the STFC
Scientific Computing Department’s SCARF cluster. UK Ministry of Defence ©
Crown Owned Copyright 2023/AWE.

Author contributions This project was conceived by G.G. and R.B. The
experiment was designed by C.D.A., G.G., P.S., N.C. and R.B., and carried out
by C.D.A., P.S., N.C., G.G., T.H., R.S., J.H., P.B., S.B., F.D.C., A.M.G.,
D.H., S.I., V.S., T.V. and B.T.H. The data analysis was carried out by C.D.A.
The manuscript was written by C.D.A., with input from G.G., N.C. and R.B.
Numerical simulations were performed by C.D.A. and P.S. Further experimental
and theoretical support was provided by I.E., D.H.F., A.F.A.B., A.A.S.,
J.T.G., B.T.H., F.M., S.S., B.R., H.C., L.O.S., T.D. and R.M.G.M.T.

## Supplementary information

### Detailed secondary beam characteristics from FLUKA simulations

The yield, mean energy and beam radius of all secondary species produced by
the proton beam irradiation of the target are given in Table 2, obtained at
the rear surface of the target. The particle spectra of electrons, positrons
and protons are given as a function of relativistic Lorentz factor in Figure
5.

Species | Specified energy range | Yield (per primary) | Mean energy, $\langle E\rangle$ (GeV) | Beam radius, $\sigma_{r}$ (mm)
---|---|---|---|---
e- | $\geq$1\text{\,}\mathrm{MeV}$$ | 56 | 0.20 | 2.0
e+ | $\geq$1\text{\,}\mathrm{MeV}$$ | 46 | 0.27 | 1.9
p+ |  $\geq$100\text{\,}\mathrm{keV}$$, $<$430\text{\,}\mathrm{GeV}$$ | 1.4 | 31 | 1.2
p+ | $\geq$430\text{\,}\mathrm{GeV}$$ | 0.42 | 439 | 1.2
$\pi^{\pm}$ | $\geq$100\text{\,}\mathrm{keV}$$ | 6.2 | 6.6 | 1.7
k± | $\geq$100\text{\,}\mathrm{keV}$$ | 0.5 | 11 | 1.4
µ- | $\geq$100\text{\,}\mathrm{keV}$$ |  $<$0.05 | - | -
$\gamma$ | $\geq$10\text{\,}\mathrm{keV}$$ | 600 | 0.07 | 2.0
n | $\geq$100\text{\,}\mathrm{keV}$$ | 12.6 | 3.1 | 2.2

Table 2: Detailed secondary beam characteristics from FLUKA simulations. The
yield (per primary), mean energy ($\langle E\rangle=\langle\Gamma\rangle
mc^{2}$), and beam radius (half-width-half-maximum) are provided for the most
abundantly produced secondary species: electrons (e-), positrons (e+), protons
(p+), pions ($\pi^{\pm}$), kaons (k±), muons (µ-), $\gamma$-rays and neutrons
(n). The specified energy range defines the range of energies that the table
parameters refer to. In the case of electrons and positrons, for ease of
comparison with other studies, the low-energy cutoff used is
$1\text{\,}\mathrm{MeV}$, though there are even more pairs (up to 10% more) if
the energy cutoff is reduced to $10\text{\,}\mathrm{keV}$. For the remaining
species, the low energy cutoff is defined by the simulation. The yield per
primary can be multiplied by the primary proton beam intensity to obtain the
actual yield. For instance, with a primary proton beam intensity of $3\times
10^{11}$ protons, a total of $1.7\times 10^{13}$ electrons and $1.4\times
10^{13}$ positrons are produced. For the muons, the simulated number is too
small for reasonable statistics of the mean energy and beam radius.  Figure 5:
FLUKA simulations of particle spectra. Particle spectra of electrons (e-,
blue), positrons (e+, red) and protons (p+, green) at the exit of the target
are obtained from FLUKA simulations and plotted as a function of relativistic
Lorentz factor, $\Gamma=E/mc^{2}$. The electron and positron spectra are
matched except at energies $\lesssim$10\text{\,}\mathrm{MeV}$$, where
electrons can more easily escape the target than positrons. The e-e+ spectra
are characterized by extended tail distribution functions with energies up to
several tens of GeV. The population of primary protons which experience only
elastic scattering are identified by the peak in the spectra corresponding to
proton momentum $p=$440\text{\,}\mathrm{GeV}$$/c.

### Chromox luminescence as a function of particle type and energy

To further validate the assumption that particles deposit an approximately
equal amount of energy into the screens, FLUKA simulations are performed,
irradiating Chromox with different particle types and energies. The results
are shown in Figure 6. In the range $2\lesssim\Gamma\lesssim 500$, the
deposited energy per particle is approximately the same across energy scales
and between charged particle species. The mean energy deposited for
$\gamma$-rays in the range MeV-TeV is plotted (black-dotted), showing that the
energy deposition per $\gamma$-ray is approximately $300\times$ smaller.
Therefore, we expect the contribution of $\gamma$-rays to observed screen
luminescence to be negligible, despite the expected yield being $6\times$
greater than e+e-.

Figure 6: FLUKA simulations of energy deposited into Chromox screens. The
energy deposited per particle into a $250\text{\,}\mathrm{\SIUnitSymbolMicro
m}$-thickness Chromox screen at $45\text{\,}\mathrm{\SIUnitSymbolDegree}$
incidence is obtained from FLUKA simulations. Deposited energy is plotted as a
function of relativistic Lorentz factor for: electrons (e-, blue-solid),
positrons (e+, red-dashed), protons (p+, green-dot-dashed), charged pions
($\mathrm{\pi}^{\pm}$, orange-dashed) and charged kaons (k±, purple-dot-
dashed). The mean energy deposited per $\gamma$-ray in the range MeV-TeV is
also plotted (black-dotted), showing that the energy deposition of
$\gamma$-rays is approximately $300\times$ smaller.

### Magnetic particle spectrometer raw data

Raw image data of the Chromox screens in the magnetic particle spectrometer is
shown in Figure 7 (post-background subtraction). The increase in signal when
the electromagnet is turned on is evidence of positrons in the beam with
energy $E\leq$220\text{\,}\mathrm{MeV}$$. The horizontal lineouts from this
data in the central 25 mm are used to obtain the electron and positron spectra
shown in Figure 3.

Figure 7: Magnetic particle spectrometer raw data. Raw image data of the
chromium-doped luminescence screens used to measure the energy spectra of
electrons and positrons in the magnetic particle spectrometer. The screen
corresponding to positrons is shown when the target is irradiated and the
electromagnet is turned off (upper), and when the electromagnet is turned on
with peak magnetic field $B=$0.34\text{\,}\mathrm{T}$$ (lower). When the
magnet is turned on, positrons are deflected away from the beam axis onto the
screen. The secondary x-axis shows the corresponding positron energy according
to the energy calibration when the magnet is operating at the given setting.

### Fitting a relativistic Maxwellian to the pair beam momentum distribution

Given that the pair plasma beams produced in this work are characterized by
non-Maxwellian distribution functions, possessing relativistic thermal spreads
and bulk flow ($k_{\rm B}T_{\pm}\gg m_{\rm e}c^{2}$), comparisons with
lengthscales traditionally used to define collective behaviour in classical
plasmas are done with care. The pair temperature ($T_{\pm}$) is identified as
a key parameter which must be included in the consideration of the limiting
scales of collective plasma processes.

We derive the screening length for a relativistically hot plasma assuming a
relativistic Maxwellian (Jüttner-Synge) distribution (Eq. 4). This screening
lengthscale is analogous to the familiar Debye screening length for non-
relativistic Maxwellian plasmas. Since the screening length is derived
assuming zero bulk motion, to obtain the appropriate pair temperature we
perform a Lorentz transformation on the FLUKA-simulated momentum spectra and
fit the data to a Jüttner-Synge distribution. The number of pairs per Debye
sphere is calculated in the zero-momentum frame, using the beam density and
dimensions consistent with the Lorentz transformation from the laboratory
frame. The fitting of the Lorentz-transformed momentum distributions is shown
in Figure 8.

The number of pairs per Debye sphere is found to scale with:

$\frac{N_{\pm}}{N_{\rm D}}\propto
N_{\pm}~{}{n_{\pm}}^{1/2}~{}\Gamma_{0}^{-1/2}~{}\Gamma_{T}^{-3/2},$ (6)

where $\Gamma_{0}$ is the Lorentz factor required for a transformation from
the laboratory frame to the zero-momentum frame, and $\Gamma_{T}$ is the
Lorentz factor associated with the pair temperature fitted in the zero-
momentum frame ($k_{\rm{B}}T_{\pm}=\Gamma_{T}m_{\rm e}c^{2}$).

Alternative approaches of comparing the pair beam with limiting scales of
plasma collective processes do not consider the thermal spread of the beam in
the screening length. Instead, a comparison is made with the plasma skin
depth, using a density corresponding to a Lorentz boost to a frame co-moving
with the mean pair energy ($\langle E\rangle$ = $\langle\Gamma\rangle m_{\rm
e}c^{2}$) [21]. There are two issues with this approach. The first is that the
pair temperature is not considered, and so the treatment is as if the momentum
distribution is mono-energetic at the mean particle energy. Secondly, particle
spectra with extended tails can have mean Lorentz factor far-in-excess of the
median Lorentz factor, therefore misrepresenting the bulk. This alternative
approach leads to a different scaling for the number of particles per
screening volume: $(N_{\pm}/N_{\rm D})\propto
N_{\pm}~{}{n_{\pm}}^{1/2}~{}\langle\Gamma\rangle^{-3/2}$. However, the
approach presented in this work leads to approximately the same result if
$\Gamma_{0},\Gamma_{T}\sim\langle\Gamma\rangle^{3/4}$. Indeed this the case
given $\Gamma_{0}=4.5$ and $\Gamma_{T}=8$, and provided
$\langle\Gamma\rangle\sim 10-20$, which is generally true for the laser
experiments compared in Figure 4.

Figure 8: Fitting the momenta distributions to relativistic Maxwellian
distributions. (a) The longitudinal momentum distribution of positrons,
$p_{\rm{e}^{+}}^{\ell}$, corresponding to the laboratory frame
($\Gamma_{0}=1$) is obtained from FLUKA simulations (red). (b) Once a Lorentz
transformation ($\Gamma_{0}=4.5$) is applied to the particle vectors, the
longitudinal momentum spectrum is modified such that it becomes approximately
symmetric about zero-momentum. A relativistic Maxwellian distribution is
fitted with temperature parameter $\Gamma_{T}=8$ (black-dashed). (c) The
transverse momentum distribution of positrons, $p_{\rm{e}^{+}}^{t}$, is
plotted with the same fit.