# A Programmatic and Semantic Approach to Explaining and Debugging
Neural Network Based Object Detectors

Edward Kim
UC Berkeley
<EMAIL_ADDRESS>Gopinath, Divya
SGT Inc, NASA AMES
<EMAIL_ADDRESS>Corina Pasareanu
NASA AMES, Carnegie Mellon University
<EMAIL_ADDRESS>Sanjit A. Seshia
UC Berkeley
<EMAIL_ADDRESS>
###### Abstract

Even as deep neural networks have become very effective for tasks in vision
and perception, it remains difficult to explain and debug their behavior. In
this paper, we present a programmatic and semantic approach to explaining,
understanding, and debugging the correct and incorrect behaviors of a neural
network-based perception system. Our approach is semantic in that it employs a
high-level representation of the distribution of environment scenarios that
the detector is intended to work on. It is programmatic in that scenario
representation is a program in a domain-specific probabilistic programming
language which can be used to generate synthetic data to test a given
perception module. Our framework assesses the performance of a perception
module to identify correct and incorrect detections, extracts rules from those
results that semantically characterizes the correct and incorrect scenarios,
and then specializes the probabilistic program with those rules in order to
more precisely characterize the scenarios in which the perception module
operates correctly or not. We demonstrate our results using the Scenic
probabilistic programming language and a neural network-based object detector.
Our experiments show that it is possible to automatically generate compact
rules that significantly increase the correct detection rate (or conversely
the incorrect detection rate) of the network and can thus help with
understanding and debugging its behavior.

Figure 1: Overview of the workflow proposed in this paper. The green and red
bounding boxes are ground truth and prediction, respectively.

##  1 Introduction

Models produced by Machine Learning (ML) algorithms, especially deep neural
networks (DNNs), have proved very effective at performing various tasks in
computer vision and perception. Moreover, ML models are being deployed in
domains where trustworthiness is a big concern, such as automotive systems
[18], health care [2], and cyber-security [6]. Research in adversarial machine
learning [11], verification [8, seshia-arxiv16], and testing [27] has shown
that DNN-based vision/perception systems are not always robust and can be
fooled, sometimes leading to unsafe situations for the overall system (e.g.,
autonomous vehicle).

Given this lack of robustness and potential for unsafe behavior, it is crucial
that we develop methods to better understand, debug, and characterize
scenarios where DNN-based perception components fail and where they perform
correctly. The emerging literature on explaining and understanding ML models
provides one approach to address this concern. However, while there are
several techniques proposed to explain the behavior of ML-based perception
(e.g. [5, 16, 17, 19, 24]), almost all of them operate on the concrete input
feature space of the network. For example, attribution-based methods(e.g. [25,
30, 23]) indicate pixels in an input image that are associated with the output
of a DNN on that input. These methods, while very useful, do not directly
identify the higher-level “semantic” features of the scene that are associated
with that decision; they require a human to make that judgment. Additionally,
in many cases it is important to generate “population-level” explanations of
correct/incorrect behavior on such higher-level features. For example, it
would be useful to identify whether the perception module of an autonomous
vehicle generally misses cars of a certain model or color, or on a particular
region of a road, and leverage this knowledge to describe a high-level
success/failure scenario of a perception module without the bottleneck of
human intervention.
In this paper, we present a programmatic and semantic approach to explaining
and debugging DNN-based perception module, with a focus on object detection.
In this approach, we begin by formalizing the semantic feature space as a
distribution over a set of scenes, where a scene is a configuration of objects
in the three dimensional space and semantic features are features of the scene
that capture its semantics (e.g., the position and orientation of a car, its
model and color, the time of day, weather, etc.). We then represent the
semantic feature space using a program in a domain-specific programming
language – hence the term programmatic. Given such a representation and
generated data corresponding to correct and incorrect behaviors of an object
detector, we seek to compute specializations of the program corresponding to
those correct/incorrect behaviors. The specialized programs serve as
interpretable representations of environment scenes that result in those
correct/incorrect behaviors, enabling us to debug failure cases and to
understand where the object detector succeeds.

We implement our approach using the Scenic [1, 9] probabilistic programming
language. Probabilistic programming has already been demonstrated to be
applicable to various computer vision tasks (see, e.g., [14]). Scenic is a
domain-specific language used to model semantic feature spaces, i.e.,
distributions over scenes. It has a generative back-end that allows one to
automatically produce synthetic data when it is connected to a renderer or
simulator, such as the Grand Theft Auto V (GTA-V) video game. It is thus a
particularly good fit for our approach. Using Scenic, we implement the
workflow shown in Fig. 1. We begin with a Scenic program $P$ that captures a
distribution that we would like our DNN-based detector to work on. Generating
test data from $P$, we evaluate the performance of the detector, partitioning
the test set into correct and incorrect detections. For each partition, we use
a rule extraction algorithm to generate rules over the semantic features that
are highly correlated with successes/failures of the detector. Rule extraction
is performed using decision tree learning and anchors [22]. We further propose
a novel white-box approach that analyzes the neuron activation patterns of the
neural network to get insights into its inner workings. Using these activation
patterns, we show how to derive semantically understandable rules over the
high-level input features to characterize scenarios.

The generated rules are then used to refine $P$ yielding programs $P^{+}$ and
$P^{-}$ that characterize more precisely the correct and incorrect feature
spaces, respectively. Using this framework, we evaluate DNN-based object
detector for autonomous vehicles, using data generated using Scenic and GTA-V.
We demonstrate that our approach is very effective, producing rules and
refined programs that significantly increase the correct detection rate (from
$65.3\%$ to $89.4\%$) and incorrect detection rate (from $34.7\%$ to $87.2\%$)
of the network and can thus help with understanding, debugging and retraining
the network.

In summary, we make the following contributions:

  * $\bullet$

Formulation of a programming language-based semantic framework to characterize
success/failure scenarios for an ML-based perception module as programs that
help delineate its performance boundaries and generate new data in a
principled way;

  * $\bullet$

An approach based on anchors and decision tree learning for deriving rules for
refining scenario programs;

  * $\bullet$

A novel white-box technique that uses activation patterns of convolutional
neural networks to enhance scenario feature space refinement;

  * $\bullet$

A data generation platform enabling research into debugging and explaining
DNN-based perception, and

  * $\bullet$

Experimental results demonstrating that our framework is effective for a
complex convolutional neural network used in autonomous driving.

##  2 Background

Feature | Range
---|---
Weather |  Neutral, Clear, Extrasunny, Smog, Clouds,
Overcast, Rain, Thunder, Clearing, Xmas,
Foggy,Snowlight, Blizzard, Snow
Time | [00:00, 24:00)
Car Model |  Blista, Bus, Ninef, Asea, Baller, Bison, Buffalo, Bobcatxl, Dominator, Granger, Jackal, Oracle, Patriot, Pranger
Car Color | R = [0, 255], G = [0, 255], B =[0, 255]
Car Heading | [0, 360) deg
Car Position | Anywhere on a road on GTA-V’s map
Table 1: Environment features and their ranges in GTA-V

Scenic [9, 1] is a probabilistic programming language for scenario
specification and scene generation. The language can be used to describe
environments for autonomous systems, i.e. autonomous cars or robots. Such
environments are scenes, i.e. configurations of objects and agents. Scenic
allows assigning distributions to the features of the scenes, as well as hard
and soft mathematical constraints over the features in the scene. Generating
scenes from a Scenic program requires sampling from the distributions defined
in the program. Scenic comes with efficient sampling techniques that take
advantage of the structure of the Scenic program, to perform aggressive
pruning. The generated scenes are rendered into images with the help of a
simulator. In this paper (and similar to [9]) we use the Grand Theft Auto V
(GTA-V) game engine [10] to create realistic images with a case study that
uses SqueezeDet [29], a convolutional neural network for object detection in
autonomous cars. Note that the framework we put forth is not specific to
SqueezeDet, and can be used with other object detectors as well.

The semantic features that we use in our case study are described in Table 1.
These features are determined and limited by the environment parameters that
are controllable by users in the simulator. If distributions over these
environment features are not specified in a Scenic program, then, by default,
they are uniformly randomly selected from ranges shown in Table 1. Note that
for a different application domain, we would have a different set of features.

Scenic is designed to be easily understandable, with simple and intuitive
syntax. We illustrate it via an example, shown in Figure 2. The formal syntax
and semantics can be found in [9].

As shown in Figure 2, the program describes a rare situation where a car is
illegally intruding over a white striped traffic island to either cut in or
belatedly avoid entering elevated highway. In line 1, ”param time = (6*60,
18*60)” means that time of the day is uniformly randomly sampled from 6:00 to
18:00. In line 2, an ego car is placed at specific x @ y coordinate on GTA-V’s
map. In line 4, a spot on a traffic island that is within a visible region
from a camera mounted on ego car is selected. Of all visible region of the
traffic island, a spot is uniformly randomly sampled. In line 7 and 8,
otherCar is placed on the spot facing -90 to 90 degree off of where ego car is
facing, simulating cases when

Figure 2: Example Scenic program

a car may be protruding into a traffic flow. Lastly, scenic allows users to
define hard and soft constraints using ”require” statements. In this scenario,
all four require statements define hard constraints. In line 10, entire
surface of the otherCar must be within the view region of the ego car. So, a
scene where only front half of the otherCar is visible is not allowed. In line
11, the otherCar must be positioned in the right half of the ego car’s visible
region. In line 12 and 13, the distance of the otherCar from ego car should be
5 to 20 meters.

##  3 Related Work

Most techniques that aim to provide explainability and interpretability for
deep neural networks (DNNs) in the field of computer vision focus on
attributing the network’s decisions to portions of the input images( [16, 19,
24, 25, 30]). GradCAM [23] is a popular approach for interpreting CNN models
that visualizes how parts of the image affect the neural network’s output by
looking into class activation maps (CAM). Other techniques focus on
understanding the internal layers by visualizing their activation patterns [5,
17]. Our approach, on the other hand, aims to provide characterizations at a
higher level than raw image pixels, namely at the level of abstract features
defined in a Scenic program.

Rule extraction techniques either aim to represent the entire functionality of
the network as a set of rules making it too complex [31] or require the
presence of pre-mined set of rules [15] which would be difficult to obtain for
the object detection scenario. Anchors [22], which improves on LIME [21], is
closest to our work (and we discuss it in more detail later).

Recent work aims to explain the decisions of DNNs in terms of higher-level
concepts. The technique in [13] introduces the idea of concept activation
vectors, which provide an interpretation of a neural network’s internal state
in terms of human-friendly concepts. Feature Guided Exploration [28] aims to
analyze the robustness of networks used in computer vision applications by
applying perturbations over high-level input features extracted from raw
images. They use object detection techniques (such as SIFT – Scale Invariant
Feature Transform) to extract the features from an image. In contrast to these
techniques we directly leverage Scenic which defines the high-level features
in a way that is already understandable for humans. Existing approaches
typically use classification networks whose output directly corresponds to the
decision being made and rely on the derivative of the output with respect to
the input to calculate importance. In our application, there is no direct
correlation between the output of the object detector network and the validity
of the bounding boxes. Furthermore, unlike all previous work, we can use the
synthesized rules to automatically generate more input instances, by refining
the original Scenic program and then using it to generate data. These
instances can be used to test, debug and retrain the network.

##  4 Approach

Figure 3: Analysis Pipeline

The key idea of our approach is to leverage the high-level semantic features
formally encoded in a Scenic program to derive rules (sufficient conditions)
that explain the behavior of a detection module in terms of those features.
Our hypothesis is that since these features describe the important
characteristics that should be present in an image and furthermore they are
much fewer than the raw, low-level pixels, they should lead to small, compact
rules that have a clear meaning for the developer.

The problem that our technique aims to address can be formalized as follows.
Suppose a function g defines a mapping from a feature vector,
$[f_{1},f_{2},...,f_{n}]\in D_{1}\bigtimes D_{2}\bigtimes...\bigtimes D_{n}$,
to a matrix of pixels, $m\in M$, of an image, where each $D_{i}$ represents
the feature domain of feature $f_{i}$ and $M$ is a domain of $m$. Let function
$h$ denote the given perception module. Finally, let $e$ be an evaluation
function which compares the perception module’s prediction to the ground
truths, and outputs a boolean class (correct or incorrect) based on a certain
performance threshold. Given a Scenic program, according to its feature
dependencies and hard and soft constraints, the feature space, $D_{1}\bigtimes
D_{2}\bigtimes...\bigtimes D_{n}$, is defined. The problem is to find the
subset feature space, $d_{1}\bigtimes d_{2}\bigtimes...\bigtimes
d_{n}\subseteq D_{1}\bigtimes D_{2}\bigtimes...\bigtimes D_{n}$ such that when
we sample a certain number of features $[f_{1},f_{2},...,f_{n}]\in
d_{1}\bigtimes d_{2}\bigtimes...\bigtimes d_{n}$, the probability that
$e(h(g([f_{1},f_{2},...,f_{n}])))$ is equal to a target class (correct or
incorrect) is maximized.

A high-level overview of our analysis pipeline is illustrated in Figure 3. We
start with a Scenic program that encodes constraints (and distributions) over
high-level semantic features that are relevant for a particular application
domain, in our case object detection for autonomous driving. Intuitively, the
program (henceforth called scenario) encodes the environments that the user
wants to focus on in order to test the module. Based on this scenario, Scenic
generates a set of feature vectors by sampling from the specified
distributions. A simulator is then used to generate a set of realistic,
synthetic images (i.e. raw low-level pixel values) based on those features.

The images are fed to the object detector. Each image is assigned a binary
label, correct or incorrect, based on the performance of the object detector
on the image (see Section 4.1). The labels obtained for the images are mapped
back to the feature vectors that led to the generation of the respective
images. The result is a labeled data set that maps each high-level feature
vector to the the respective label.

We then use off-the-shelf methods to extract rules from this data set. The
rule extraction is described in more detail in Sec. 4.2. The result is a set
of rules encoding the conditions on high-level features that lead to likely
correct or incorrect detection. The obtained rules can be used to refine the
Scenic program, which in turn can be sampled to generate more images that can
be used to test, debug or retrain the detection module. This iterative process
can continue until one obtains refined rules, and Scenic programs, of desired
precision. In the following we provide more details about our approach.

###  4.1 Labelling

Obtaining the label (correct/incorrect) for an image is performed using the F1
score metric (harmonic mean of the precision and recall). This metric is
commonly used in statistical analysis of binary classification.

The F1 score is computed in the following way. For each image, the true
positive (TP) is the number of ground truth bounding boxes correctly predicted
by the detection module. Correctly predicted here means intersection-over-
union (IoU for object detection) is greater than 0.5. The false positive (FP)
is the number of predicted bounding boxes that falsely predicted ground
truths. This false prediction includes duplicate predictions on one ground
truth box. The false negative (FN) is the number of ground truth boxes that is
not detected correctly. We computed the F1 score for each image, and if it is
greater than a threshold, we assigned correct label; if not, incorrect. The
threshold used in our experiments was 0.8.

###  4.2 Rule Extraction

#### Methods:

We experimented with two methods, decision tree (DT) learning for
classification [20] and anchors [22], to extract rules capturing the subspace
of the feature space defined in the given Scenic program.

Decision tree learning is commonly used to extract rules explaining the global
behavior of a complex system while the anchors method is a state-of-the art
technique for extracting explanation rules that are locally faithful.

Decision trees encode decisions (and their consequences) in a tree-like
structure. They are highly interpretable, provided that the trees are short.
One can easily extract rules for explaining different classes, by simply
following the paths through the trees and conjuncting the decisions encoded in
the tree nodes. We used the rpart [26] package in R software, which implements
corresponding algorithm in [4], with default parameters.

The anchor method is a state-of-the art technique that aims to explain the
behavior of complex ML models with high-precision rules called anchors,
representing local, sufficient conditions for predictions. The system can
efficiently compute these explanations for any black-box model with high-
probability guarantees. We used the code from
https://github.com/marcotcr/anchor with the default parameters. Applying the
method to the object detector directly would result in anchors describing
conditions on low-level pixel values, which would be difficult to interpret
and use. Instead what we want is to extract anchors in terms of high-level
features. While one can use the simulator together with the object detector as
the black-box model, this would be very inefficient. Instead we built a
surrogate model mapping high-level Scenic features to output labels; we used a
random forest learning algorithm for this purpose as in the code. This
surrogate model was then passed to the anchor method to extract the rules.

#### Blackbox vs Whitebox Analysis:

So far we explained how we can obtain rules when treating the detection module
as a black box. We also investigated a white-box analysis, to determine
whether we can exploit the information about the internal workings of the
module to improve the rule inference. The white-box analysis is one of our
novel contributions in this paper. We leverage recent work [12] which aims to
infer likely properties of neural networks. The properties are in terms of
on/off activation patterns (at different internal layers) that lead to the
same predictions. These patterns are computed by applying decision-tree
learning over the activations observed during the execution of the network on
the training or testing set.

We analyzed the architecture of the SqueezeDet network and we determined that
there are three maxpool layers which provide a natural decomposition of the
network. Furthermore they have relatively low dimensionality making them a
good target for property inference.

We consider activation patterns over maxpool neurons based on whether the
neuron output is greater or equal to zero. A decision tree can then be learned
over these patterns to fit the prediction labels. For our experiments we
selected patterns from the maxpool layer 5, which turned out to be highly
correlated to images that lead to correct/incorrect predictions.

Then, we augmented the assigned correct and incorrect labels with
corresponding decision pattern in the following way. For example, using a
decision pattern for correct labels (i.e. the decision pattern that most
correlated to images with correct label), we created two sub-classes for
correct class. By feeding in only images with correct label to the perception
module, the images satisfying the decision pattern is re-labelled as ”correct-
decision-pattern,” otherwise, ”correct-unlabelled.” Likewise, the incorrect
class is augmented using a decision pattern that is most correlated to images
with incorrect label. It is our intuition that the decision pattern captures
more focused properties (or rules) among images belonging to a target class.
Hence, we hypothesize that this label augmentation would help anchor and
decision tree methods to better identify rules.

#### Rule Selection Criteria:

Once we extracted rules with either DT or anchors, we selected the best rule
using following criteria. To best achieve our objective, first, we chose the
rule with highest precision on a held-out testset of feature vectors. If there
are more than one rule with equal high precision, then we chose the rule with
the highest coverage (i.e. the number of feature vectors satisfying the
rules). Finally, if there is still more than one rule left, then we broke the
tie by choosing the most compact rule which has the least number of features.
The last two criteria are established to select the most general high-
precision rule.

##  5 Experiments

Scenario # | Rules
---|---
(Baseline$\rightarrow$Rule Precision) |
Scenario 1 |  x coordinate $\geq$ -198.1
($65.3\%\rightarrow 89.4\%$) |
|  hour $\geq$ 7.5 $\wedge$
|  weather = all except neutral $\wedge$
Scenario 2 |  car0 distance from ego $\geq$ 11.3m $\wedge$
($72.3\%\rightarrow 82.3\%$) |  car0 model = {Asea, Bison, Blista, Buffalo, Dominator, Jackal, Ninef, Oracle}
Scenario 3 |  car0 red color $\geq$ 74.5 $\wedge$
($61.7\%\rightarrow 79.4\%$) |  car0 heading $\geq$ 220.3 deg
|  car0 model = {Asea, Baller, Blista,
Scenario 4 |  Buffal, Dominator, Jackal, Ninef,
($89.6\%\rightarrow 96.2\%$) |  Oracle}

Table 2: Rules for correct behaviors of the detection module with the highest precision from Table 6 Scenario # | Rules
---|---
(Baseline$\rightarrow$Rule Precision) |
|  x coordinate $\leq$ -200.76 $\wedge$
Scenario 1 |  distance $\leq$ 8.84 $\wedge$
($34.7\%\rightarrow 87.2\%$) |  car model = PRANGER
|  hour $\geq$ 7.5 $\wedge$
Scenario 2 |  weather = all except Neutral $\wedge$
($27.7\%\rightarrow 44.9\%$) |  car0 distance from ego $<$ 11.3
|  weather = neutral $\wedge$
Scenario 3 |  agent0 heading = $\leq$ 218.08 deg $\wedge$
($38.3\%\rightarrow 83.4\%$) |  hour $\leq$ 8.00 $\wedge$
|  car2 red color $\leq$ 95.00
|  car0 model = PATRIOT $\wedge$
|  car1 model = NINEF $\wedge$
Scenario 4 |  car2 model = BALLER $\wedge$
($10.4\%\rightarrow 57.3\%$) |  92.25 $<$ car0 green color $\leq$ 158 $\wedge$
car0 blue color $\leq$ 84.25 $\wedge$
178.00 $<$ car2 red color $\leq$ 224

Table 3: Rules for incorrect behaviors of detection module with the highest
precision from Table 7

In this section we report on our experiments with the proposed approach on the
object detector. We investigate whether we can synthesize rules that are
effective in generating test inputs that increase the probability of
correct/incorrect detection, thus explaining the correct/incorrect behavior of
the analyzed module. We evaluate the proposed techniques along the following
dimensions: decision tree (DT) vs anchor, black-box (BB) vs white-box (WB).

###  5.1 Scenarios

We experimented with our approach on four different scenarios. Scenario 1
(Figure 2) describes the situation where a car is illegally intruding over a
white striped traffic island at the entrance of an elevated highway. Scenario
2 describes two-car scenario where one car occludes the ego car’s view of
another car at a T-junction intersection on an elevated road. describes scenes
where other cars are merging into ego car’s lane. The location in this
scenario is carefully chosen such that the sun rises in front of ego car,
causing a glare. describes a set of scenes when nearest car is abruptly
switching into ego car’s lane while another car on the opposite traffic
direction lane is slightly intruding over the middle yellow line into ego
car’s lane.

Scenario 2 describes two-car scenario where a car occludes the ego car’s view
of another car at a T-junction intersection . In the Scenic program, to cause
occlusion in scenes, we place an agent0 car within a visible region from ego
car. Then, we place agent1 car within a close vicinity, defined by small
horizontal (i.e. leftRight) and vertical (i.e. frontback) perturbations in the
program, to agent0 car. The metric of these perturbations are in meters. The
images from Scenario 2 are shown in Figure 8.

Scenario 3 describes scenes where three cars are merging into the ego car’s
lane. The location for Scenario 3 is carefully chosen such that the sun rises
in front of the ego car, causing a glare. The Scenic program describes three
cars merging in a platoon-like manner where one car is following another car
in front with some variations in distance between front and rear cars. The
metric for distance perturbation is in meters. The images from Scenario 3 are
shown in Figure 9.

Finally, Scenario 4 describes a set of scenes when the nearest car is abruptly
switching into ego car’s lane while another car on the opposite traffic
direction lane is slightly intruding over the middle yellow line into the ego
car’s lane. Failure to detect these two cars out of four cars may potentially
be safety-critical. The images from Scenario 4 are shown in Figure 10. The
locations of all four cars, in Scenario 4 Scenic program, are hard-coded with
respect to ego car’s location. The Scenic program would have become much more
interpretable had we described car locations with respect to lanes. The reason
we had to code in this undesirable manner is due to the simulator as
illustrated later in this section.

##  6 Success and Failure Scenario Descriptions

The refined Scenic programs characterizing success/failure scenarios are shown
in Figure 11, 12, 13, and 14. The red/green parts of programs represent the
rules automatically generated by our technique, which are cut and pasted to
original Scenic programs. These success/failure inducing rules are shown in
Table 2 and 3. As aforementioned, we generated new images using Scenic
programs that characterize failure scenarios. Examples of these images from
failure scenarios are shown in Figure 15, 16, 17, and 18.

###  6.1 Setup

The object detector was trained on a separate set of 10,000 GTA images with
one to four cars in various locations of the map producing different
background scenes. The GTA-V simulator provided images, ground truth boxes,
and values of the environment features.

For each scenario, we generated 950 new images as a train set and another 950
new images as a test set. We denote the labels corresponding to the maxpool
layer 5 decision pattern as p5c(correct) and p5_ic(incorrect) and the
remaining as correct_unlabelled and incorrect_unlabelled, respectively. We
augmented the feature vector with some extra features that are not part of the
feature values provided by the simulator but could help with extracting
meaningful rules. For example, in Scenario 1, the distance from ego to
otherCar is not part of the feature values provided by GTA-V. However, it can
be computed with Euclidean distance metric using (x,y) location coordinates of
ego and otherCar. Also, the difference in heading angle between ego and
otherCar is also added as extra feature to represent “badAngle” variable in
the program.

From the train set, we extracted rules to predict each label based on the
feature vectors. These rules were evaluated on the test set based on
precision, recall, and F1 score metrics. For DT learning we adjusted the label
weight to account for the uneven ratio among labels for both black-box and
white-box labels. For the Anchors method, we applied it on each instance of
the training set until we had covered a maximum of 50 instances for every
label (correct, incorrect for Black Box, and p5c, p5_ic, correct_unlabelled,
incorrect_unlabelled for White Box). The best anchor rule for every label is
selected based on the rule selection criteria mentioned in section 4.2.

Figure 4: The cumulative ratio of incorrectly detected images generated from refined Scenic programs (using incorrect rules) stabilizes over 500 samples. Each color has four graphs representing four different rule extraction methods Scenario # | 1 | 2 | 3 | 4
---|---|---|---|---
Correct DP | 0.626 | 0.651 | 0.514 | 0.824
Incorrect DP | 0.276 | 0.175 | 0.234 | 0.212

Table 4: Support for correct and incorrect decision patterns Scenario # | 1 | 2 | 3 | 4
---|---|---|---|---
BB Decision Tree | 0.723 | 0.342 | 0.631 | 0.622
WB Decision Tree | 0.727 | 0.696 | 0.601 | 0.778
BB Anchor | 0.361 | 0.457 | 0.302 | 0.438
WB Anchor | 0.520 | 0.188 | 0.149 | 0.438

Table 5: F1 score of correct rules on testset Scenario # | 1 | 2 | 3 | 4
---|---|---|---|---
Original Program | 0.653 | 0.723 | 0.617 | 0.896
BB Decision Tree | 0.843 | 0.778 | 0.787 | 0.950
WB Decision Tree | 0.826 | 0.823 | 0.788 | 0.962
BB Anchor | 0.727 | 0.811 | 0.652 | 0.928
WB Anchor | 0.894 | 0.817 | 0.794 | 0.928

Table 6: Precision of correct rules on testset

###  6.2 Results

Tables 2 and 3 show the best rules (wrt. precision) extracted with our
proposed framework, along with the baseline correct/incorrect detection rate
for each given scenario and the detection rate for the generated rules. The
results indicate that indeed our framework can generate rules that increase
significantly the correct and incorrect detection rate of the module.
Furthermore, the generated rules are compact and easily interpretable.

For example, the rule for correct behavior for Scenario 1 is ”x coordinate
$\geq-198.1$.” In GTA-V, at ego car’s specific location, the condition on x
coordinate was equivalent to the otherCar’s distance from ego being greater
than 11m. On the other hand, the rule for incorrect behavior for Scenario 1
requires the otherCar to be within 8.84m and its car model to be PRANGER.
These rules, counter-intuitively, indicate that the object detector fails when
the otherCar is close by, and performs well when located further away.

Scenario # | 1 | 2 | 3 | 4
---|---|---|---|---
Original Program | 0.347 | 0.277 | 0.383 | 0.104
BB Decision Tree | 0.703 | 0.418 | 0.506 | 0.375
WB Decision Tree | 0.73 | 0.449 | 0.494 | 0.099
BB Anchor | 0.872 | 0.357 | 0.834 | 0.573
WB Anchor | 0.674 | 0.422 | 0.365 | 0.176

Table 7: Incorrectly detected image ratio among five hundred new data generated from each refined Scenic program Scenario # | 1 | 2 | 3 | 4
---|---|---|---|---
Feature Space Coverage | 0.692 | 0.956 | 0.898 | 0.871

Table 8: The proportion of original feature space covered by the best
incorrect rule for each scenario from Table 7

#### Results for Correct Behavior:

Tables 5 and 6 summarize the results for the rules explaining correct
behavior. The results indicate that there are clear signals in the heavily
abstracted feature space and they can be used effectively for scenario
characterization via the generated high-precision rules.

The results also indicate that DT learning extracts rules with better F1
scores for all scenarios as compared to anchors. This could be attributed to
the difference in the nature of the techniques. The anchor approach aims to
construct rules that have high precision in the locality of a given instance.
Decision-trees on the other hand aim to construct global rules that
discriminate one label from another. Given that a large proportion of
instances were detected correctly by the analyzed module, the decision tree
was able to build rules with high precision and coverage for correct behavior.

The results also highlight the benefit of using white-box information to
extract rules for correct behavior.

Table 4 shows the support for the decision pattern is significant (greater
than 65% on average for all scenarios). The support is defined as a
correlation of the decision pattern to a specific label. Using this
information to augment the labels of the dataset helped to improve the
precision and F1 score of the rules (w.r.t. Scenic features) for both DT
learning and anchor method.

#### Results for Incorrect Behavior:

Tables 7 summarize the results for the rules explaining incorrect behavior.
Rule derivation for incorrect behavior is more challenging than for correct
behavior due to the low percentage of inputs that lead to the incorrect
detection for a well trained network.

In fact the F1 scores (computed on the test set) for rules predicting
incorrect behavior were too low due to very low (in some cases 0) recall
values.

To properly validate the efficacy of the generated rules, we refined the
Scenic programs by encoding the rules as constraints and we generated 500 new
images. We then evaluated our module’s performance on these new datasets.
Figure 4 justifies our choice of 500 as the number of new images that we
generate for evaluation.

All four methods contributed to more precisely identifying the subset features
spaces in which the module performs worse. Specifically, Table 7 illustrates
that the black-box anchor method enhanced the generation rate of incorrectly
detected images by 48% on average in Scenarios 1, 3, and 4 compared to the
baseline. This is a significant increase in the ratio of incorrectly labelled
images generated from the program, providing evidence that the refined
programs are more precisely characterizing the failure scenarios.

We also note that the anchor method outperforms DT learning. This is expected,
because the anchor method extracts rules that are highly precise within a
local feature space. The exception is Scenario 2. We conjecture that the
reason that the anchor method did not perform better than DT learning is due
to uncontrollable non-determinism in GTA-V, which generated pedestrians in
close vicinity to the camera of ego car even though its Scenic program did not
have any pedestrian. GTA-V non-deterministically instantiated these
pedestrians, and the perception module often incorrectly predicted the
pedestrians as cars. This is an issue with the GTA-V which originally was not
built for data generation purpose. GTA-V does not allow users to control or
eliminate these pedestrians and it does not provide features related to
pedestrians during data collection process. In future work, we plan to
incorporate simulators that allows a deterministic control (such as CARLA [7])
for further experimentation.

Unlike the results for correct behavior, the whitebox approach tends to
perform worse than blackbox when focusing on incorrect behavior. This outcome
can be attributed to very low support for decision patterns computed for
incorrect behavior, with maximum of 27.6% among the four scenarios as shown in
Table 4.

However, we do observe that the white-box approach for both DT learning and
anchors does, in general, enhance the ratio of incorrectly detected images as
shown in Table 7, compared to those of the original programs.

#### Limitations:

Our technique relies on abstracting an image with a high resolution (for
instance 1920 x 1200 in our example) to a vector of a small set of semantic
features. In our experiments we were able to derive compact rules with high
precision and coverage. However, we do note that in other application domains,
other than autonomous driving, the abstraction may lead to under-determined
representation, which may not yield any noticeable patterns. Therefore,
appropriate selection of a subset of essential features for a given
application domain (facilitated by an appropriate definition using Scenic), is
essential. We also note that all the Scenic programs we experimented with
contained only uniform distributions. Also, for each of the scenario programs
that we analyzed, we fixed the location and heading angle of the camera. In
these restricted settings, we were able to extract rules that distinguished
correctly detected scenes from the incorrect ones.

##  7 Conclusion and Future Work

We presented a semantic and programmatic framework for characterizing success
and failure scenarios of a given perception module in the form of programs.
The technique leverages the Scenic language to derive rules in terms of high-
level, meaningful features and generates new inputs that conform with these
rules. For future work, we plan on applying this approach to other domains, by
looking into more general input distributions and transformations.

##  8 Acknowledgment

We thank Daniel Fremont for his help on our use of Scenic, Jinkyu Kim and
Taesung Park for their thorough comments, and Xiangyu Yue for interfacing
GTA-V with Scenic. This project is supported by an NSF graduate fellowship
(Grant#: DGE1752814), NSF grants CNS-1545126 (VeHICaL), CNS-1739816, and
CCF-1837132, by the DARPA Assured Autonomy program, by Berkeley Deep Drive,
and by Toyota under the iCyPhy center. This work was partially done under the
NASA Ames Internship Program, 2019.

## References

  * [1] The SCENIC Probabilistic Programming Language.  https://github.com/BerkeleyLearnVerify/Scenic.
  * [2] Babak Alipanahi, Andrew Delong, Matthew T Weirauch, and Brendan J Frey.  Predicting the sequence specificities of DNA-and RNA-binding proteins by deep learning.  Nature biotechnology, 2015.
  * [3] Nese Alyuz and Tevfik Metin Sezgin.  Interpretable machine learning for generating semantically meaningful formative feedback.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR) Workshops, June 2019.
  * [4] L. Breiman, J. Friedman, R. Olshen, and C. Stone.  Classification and Regression Trees.  Wadsworth and Brooks, Monterey, CA, 1984.
  * [5] Shan Carter, Zan Armstrong, Ludwig Schubert, Ian Johnson, and Chris Olah.  Activation atlas.  Distill, 2019.  https://distill.pub/2019/activation-atlas.
  * [6] George E Dahl, Jack W Stokes, Li Deng, and Dong Yu.  Large-scale malware classification using random projections and neural networks.  In Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3422–3426. IEEE, 2013\.
  * [7] Alexey Dosovitskiy, German Ros, Felipe Codevilla, Antonio Lopez, and Vladlen Koltun.  CARLA: An open urban driving simulator.  In Proceedings of the 1st Annual Conference on Robot Learning, pages 1–16, 2017.
  * [8] Tommaso Dreossi, Alexandre Donze, and Sanjit A. Seshia.  Compositional falsification of cyber-physical systems with machine learning components.  In Proceedings of the NASA Formal Methods Conference (NFM), pages 357–372, May 2017.
  * [9] Daniel J. Fremont, Tommaso Dreossi, Shromona Ghosh, Xiangyu Yue, Alberto L. Sangiovanni-Vincentelli, and Sanjit A. Seshia.  Scenic: A language for scenario specification and scene generation.  In Proceedings of the 40th annual ACM SIGPLAN conference on Programming Language Design and Implementation (PLDI), June 2019.
  * [10] Rockstar Games.  Grand theft auto v.  Windows PC version, 2015.
  * [11] Ian Goodfellow, Patrick McDaniel, and Nicolas Papernot.  Making machine learning robust against adversarial inputs.  Communications of the ACM, 61(7):56–66, 2018.
  * [12] Divya Gopinath, Hayes Converse, Corina S. Pasareanu, and Ankur Taly.  Property inference for neural networks.  CoRR, abs/1904.13215, 2019.
  * [13] Been Kim, Martin Wattenberg, Justin Gilmer, Carrie J. Cai, James Wexler, Fernanda B. Viégas, and Rory Sayres.  Interpretability beyond feature attribution: Quantitative testing with concept activation vectors (tcav).  In Jennifer G. Dy and Andreas Krause, editors, ICML, volume 80 of Proceedings of Machine Learning Research, pages 2673–2682. PMLR, 2018\.
  * [14] Tejas D Kulkarni, Pushmeet Kohli, Joshua B Tenenbaum, and Vikash Mansinghka.  Picture: A probabilistic programming language for scene perception.  In Proceedings of the ieee conference on computer vision and pattern recognition, pages 4390–4399, 2015.
  * [15] Himabindu Lakkaraju, Stephen H. Bach, and Jure Leskovec.  Interpretable decision sets: A joint framework for description and prediction.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, San Francisco, CA, USA, August 13-17, 2016, pages 1675–1684, 2016.
  * [16] Scott M Lundberg and Su-In Lee.  A unified approach to interpreting model predictions.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30, pages 4765–4774. Curran Associates, Inc., 2017.
  * [17] Alexander Mordvintsev, Michael Tyka, and Christopher Olah.  DeepDream,https://github.com/google/deepdream.
  * [18] NVIDIA.  Nvidia tegra drive px: Self-driving car computer, 2015.
  * [19] Zhongang Qi, Saeed Khorram, and Fuxin Li.  Visualizing deep networks by optimizing with integrated gradients.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR) Workshops, June 2019.
  * [20] J. Ross Quinlan.  Induction of decision trees.  Machine learning, 1(1):81–106, 1986.
  * [21] Marco Túlio Ribeiro, Sameer Singh, and Carlos Guestrin.  ”why should I trust you?”: Explaining the predictions of any classifier.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, San Francisco, CA, USA, August 13-17, 2016, pages 1135–1144, 2016.
  * [22] Marco Túlio Ribeiro, Sameer Singh, and Carlos Guestrin.  Anchors: High-precision model-agnostic explanations.  In Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence, (AAAI-18), the 30th innovative Applications of Artificial Intelligence (IAAI-18), and the 8th AAAI Symposium on Educational Advances in Artificial Intelligence (EAAI-18), New Orleans, Louisiana, USA, February 2-7, 2018, pages 1527–1535, 2018.
  * [23] Ramprasaath R. Selvaraju, Michael Cogswell, Abhishek Das, Ramakrishna Vedantam, Devi Parikh, and Dhruv Batra.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  In IEEE International Conference on Computer Vision, ICCV 2017, Venice, Italy, October 22-29, 2017, pages 618–626, 2017.
  * [24] Daniel Smilkov, Nikhil Thorat, Been Kim, Fernanda B. Viégas, and Martin Wattenberg.  Smoothgrad: removing noise by adding noise.  CoRR, abs/1706.03825, 2017.
  * [25] Mukund Sundararajan, Ankur Taly, and Qiqi Yan.  Axiomatic attribution for deep networks.  CoRR, abs/1703.01365, 2017.
  * [26] Terry Therneau, Beth Atkinson, and Brian Ripley.  rpart: Recursive partitioning and regression trees.
  * [27] Yuchi Tian, Kexin Pei, Suman Jana, and Baishakhi Ray.  Deeptest: automated testing of deep-neural-network-driven autonomous cars.  In Proceedings of the 40th International Conference on Software Engineering, ICSE 2018, Gothenburg, Sweden, May 27 - June 03, 2018, pages 303–314, 2018.
  * [28] Matthew Wicker, Xiaowei Huang, and Marta Kwiatkowska.  Feature-guided black-box safety testing of deep neural networks.  In Tools and Algorithms for the Construction and Analysis of Systems - 24th International Conference, TACAS 2018, Held as Part of the European Joint Conferences on Theory and Practice of Software, ETAPS 2018, Thessaloniki, Greece, April 14-20, 2018, Proceedings, Part I, pages 408–426, 2018.
  * [29] Bichen Wu, Forrest N. Iandola, Peter H. Jin, and Kurt Keutzer.  Squeezedet: Unified, small, low power fully convolutional neural networks for real-time object detection for autonomous driving.  In Conference on Computer Vision and Pattern Recognition Workshops, CVPR Workshops, pages 446–454, 2017.
  * [30] Bolei Zhou, David Bau, Aude Oliva, and Antonio Torralba.  Interpreting deep visual representations via network dissection.  CoRR, abs/1711.05611, 2017.
  * [31] Jan Zilke, Eneldo Mencía, and Frederik Janssen.  Deepred – rule extraction from deep neural networks.  pages 457–473, 10 2016.

(a) Scenario 1 Scenic program

(b) Scenario 2 Scenic program

Figure 5: Scenario 1 and Scenario 2

(a) Scenario 3 Scenic program

(b) Scenario 4 Scenic program

Figure 6: Scenario 3 and Scenario 4  Figure 7: Images generated using Scenario
1 Figure 8: Images generated using Scenario 2 Figure 9: Images generated using
Scenario 3 Figure 10: Images generated using Scenario 4

(a) Refined Success Scenario 1 Scenic program

(b) Refined Failure Scenario 1 Scenic program

Figure 11: Refined Success and Failure Scenic programs for Scenario 1, with
red parts representing failure inducing rules and green parts representing the
success inducing rules as shown in Table 2 and 3

(a) Refined Success Scenario 2 Scenic program

(b) Refined Failure Scenario 2 Scenic program

Figure 12: Refined Success and Failure Scenic programs for Scenario 2, with
red parts representing failure inducing rules and green parts representing the
success inducing rules as shown in Table 2 and 3

(a) Refined Success Scenario 3 Scenic program

(b) Refined Failure Scenario 3 Scenic program

Figure 13: Refined Success and Failure Scenic programs for Scenario 3, with
red parts representing failure inducing rules and green parts representing the
success inducing rules as shown in Table 2 and 3

(a) Refined Success Scenario 4 Scenic program

(b) Refined Failure Scenario 4 Scenic program

Figure 14: Refined Success and Failure Scenic programs for Scenario 4, with
red parts representing failure inducing rules and green parts representing the
success inducing rules as shown in Table 2 and 3 Figure 15: Failure images
generated with Failure Scenic Scenario 1 as shown in Figure 11 (b), with
ground truth bounding boxes marked in green and prediction bounding boxes in
red Figure 16: Failure images generated with Failure Scenic Scenario 2 as
shown in Figure 12 (b), with ground truth bounding boxes marked in green and
prediction bounding boxes in red Figure 17: Failure images generated with
Failure Scenic Scenario 3 as shown in Figure 13 (b), with ground truth
bounding boxes marked in green and prediction bounding boxes in red Figure 18:
Failure images generated with Failure Scenic Scenario 4 as shown in Figure 14
(b), with ground truth bounding boxes marked in green and prediction bounding
boxes in red


[1,2]Marcel Lamač

[1]ELI Beamlines Facility, Extreme Light Infrastructure ERIC, Za Radnicí 835,
Dolní Břežany, 25241, Czechia

2]Faculty of Mathematics and Physics, Charles University, Ke Karlovu 3, Prague
2, 12116, Czechia

3]Faculty of Nuclear Science and Physical Engineering, Czech Technical
University in Prague, Břehová 7, Prague 1, 11519, Czechia

4]Kansai Photon Science Institute, National Institutes for Quantum and
Radiological Science and Technology, 8-1-7 Umemidai, Kizugawa, Kyoto,
619-0215, Japan

# Bright coherent attosecond X-ray pulses from beam-driven relativistic
mirrors (Supplementary Information)

<EMAIL_ADDRESS>Petr Valenta  Jaroslav Nejdl  Uddhab Chaulagain  Tae
Moon Jeong  Sergey V. Bulanov  *  [  [  [

###### Contents

  1. 1 Theory of weak reflection from relativistic mirrors
  2. 2 Beam-driven nonlinear plasma wave equation
  3. 3 Lorentz factor of plasma electrons and wave breaking threshold
  4. 4 Wavelength of interior plasma waves
  5. 5 Optimal driver length
  6. 6 Mirror propagation distance due to relativistic beam-plasma instability
  7. 7 Structure of beam-driven interior plasma wave breaking
  8. 8 Reflection coefficient of beam-driven relativistic mirror at the wave breaking threshold
  9. 9 Energy spectrum of reflected coherent radiation
  10. 10 Peak spectral brightness

##  1 Theory of weak reflection from relativistic mirrors

Combining Maxwell’s equations and using the Lorenz gauge condition
$\partial_{\mu}A^{\mu}=\nabla\cdot\textbf{A}+\partial_{t}\phi/c^{2}=0$, we
obtain the inhomogenous wave equation

$\left(\frac{\partial^{2}}{\partial
t^{2}}-c^{2}\nabla^{2}\right)\textbf{A}=\frac{\textbf{j}}{\epsilon_{0}},$ (1)

where A is the magnetic vector potential, $\epsilon_{0}$ is the vacuum
permittivity and j is the current density, acting as a source of radiation. We
consider only the one-dimensional geometry of a planar mirror, with a counter-
propagating incident laser. Therefore, we have for the transverse derivatives
$\partial_{y}\textbf{A}=\partial_{z}\textbf{A}=0$. To further simplify the
analysis, we consider only linear polarization in the $y$ direction and define
the dimensionless vector potential $\textbf{a}=e\textbf{A}/m_{e}c$, reducing
Eq. 1 into a scalar form. We assume that the radiative currents are due to a
propagating electron distribution oscillating in electromagnetic field,
therefore $j=-en_{e}v=-ecn_{e}a/\gamma$, where $\gamma$ is the Lorentz factor
of the plasma electrons, which gives us the one-dimensional Klein-Gordon
equation describing propagation of electromagnetic waves in plasma

$\left(\frac{\partial^{2}}{\partial t^{2}}-c^{2}\frac{\partial^{2}}{\partial
x^{2}}+\omega_{p}^{2}(x,t)\right)a(x,t)=0,$ (2)

where $\omega_{p}=\sqrt{e^{2}n_{e}(x,t)/(m_{e}\epsilon_{0}\gamma)}$ is the
relativistic plasma frequency and $n_{e}(x,t)$ is the electron density
distribution propagating with normalized velocity
$\beta=v/c=\sqrt{1-(1/\gamma)^{2}}$. Counter-propagating radiation scattered
from such a relativistic electron distribution is Doppler shifted in
frequency, which makes analysis of Eq. 2 a non-trivial task in the laboratory
frame. We therefore perform a Lorentz boost into the mirror reference frame
where both incident and reflected radiation have the same frequency
$\omega^{\prime}$. The $x$ coordinate then transforms as
$x^{\prime}=\gamma(x-vt)$. Since the operator
$\Box=(\partial_{t}^{2}-c^{2}\partial_{x}^{2})$ is Lorentz-invariant, the form
of the equation doesn’t change and it can be written in the mirror reference
frame as

$\left(\frac{\partial^{2}}{\partial t^{\prime
2}}-c^{2}\frac{\partial^{2}}{\partial x^{\prime
2}}+\omega_{p}^{\prime}\text{}^{2}(x^{\prime})\right)a(x^{\prime},t^{\prime})=0,$
(3)

where the apostrophe denotes quantities transformed to the mirror rest frame.
The normalized vector potential is transformed only in coordinates, because
the normalized amplitude ($a_{0}\propto E_{0}/\omega_{0}$) is Lorentz-
invariant.

Figure 1: Analytical solution of Eq. 9 for $k_{p}^{2}(x)=e^{-2x^{2}}$ showing
reflected radiation with light rays of constant phase propagating along the
light-cone variables $x-ct=const$.

For clarity of the following rest frame analysis, we drop apostrophes for all
the transformed variables and return them when necessary. The solution to Eq.
3 can be written generally as a superposition of the incident and reflected
waves

$a(x,t)=a_{0}e^{-i(\omega t+kx-\frac{\pi}{2})}+a_{r}(x,t),$ (4)

where $\omega$ is the radiation frequency and $k=\omega/c$ is the wave number,
both evaluated in the mirror rest frame. We include the phase factor $\pi/2$
so that the real part of the incident wave corresponds to a sine wave, such
that the real part of the incident wave is zero at spacetime origin $x=ct=0$.
Plugging (4) into (3) while assuming small reflection, $a_{r}(x,t)\ll
a_{0}\approx const.$, gives us the following equation for the reflected wave

$\left(\frac{\partial^{2}}{\partial t^{2}}-c^{2}\frac{\partial^{2}}{\partial
x^{2}}\right)a_{r}(x,t)=-i\omega_{p}^{2}(x)a_{0}e^{-i(\omega t+kx)}.$ (5)

We note that the approximation of weak reflection is equivalent to the
physically reasonable assumption that the electrons oscillate with the same
phase and normalized transverse momentum as the incident wave,
$j=-(ecn_{e}a_{0}/\gamma)e^{-i(\omega t+kx-\frac{\pi}{2})}$. The incident wave
$a_{0}(x,t)$ then represents a homogeneous solution to the differential
equation 1, while the particular solution $a_{r}(x,t)$ is obtained by solving
Eq. 1 with the approximated current.

The natural coordinate system for this problem are the light-cone coordinates
$\xi=x-ct,\eta=x+ct$, which correspond to a clockwise $\pi/4$ rotation of the
spacetime coordinate system. The derivatives transform as
$\partial_{t}+c\partial_{x}=2c\partial_{\eta}$ and
$\partial_{t}-c\partial_{x}=-2c\partial_{\xi}$, which using
$\partial_{t}^{2}-c^{2}\partial_{x}^{2}=(\partial_{t}+c\partial_{x})(\partial_{t}-c\partial_{x})$
gives us the reflected wave equation in terms of the light-cone coordinates as

$\frac{\partial^{2}a_{r}}{\partial\eta\partial\xi}(\xi,\eta)=\frac{i}{4}k_{p}^{2}\left(\frac{\eta+\xi}{2}\right)a_{0}e^{-ik\eta},$
(6)

where $k_{p}(x)=\omega_{p}(x)/c$ is the inhomogeneous relativistic plasma wave
number corresponding to the electron distribution. Eq. 6 can be solved by
direct integration and application of boundary conditions. We first apply the
boundary condition that the reflection occurs only for $t>0$, therefore at
$t=0$ we have $a_{r}(x,0)=0$. This also implies that the solution exists when
$x+ct\geq x-ct$ for all $x$. Integrating the equation over $\eta$, where
$\eta\geq\xi$, gives us the derivative of the reflected wave as

${\frac{\partial
a_{r}}{\partial\xi}(\xi,\eta)=\int_{\xi}^{\eta}\frac{i}{4}k_{p}^{2}\left(\frac{\zeta+\xi}{2}\right)a_{0}e^{-ik\zeta}d\zeta.}$
(7)

Note that, at $t=0$, Eq. 7 satisfies $\partial_{\xi}a_{r}(\eta=x,\xi=x)=0$,
i.e. no radiation is initially emitted along the light cone coordinate $\xi$.
To clearly show that Eq. 8 describes the reflected radiation, we substitute
$\zeta+\xi=2s$ in the integrand, giving us

$\frac{\partial
a_{r}}{\partial\xi}(\xi,\eta)=a_{0}\frac{ie^{ik\xi}}{2}\int_{\xi}^{\frac{\eta+\xi}{2}}k_{p}^{2}\left(s\right)e^{-i2ks}ds.$
(8)

If we assume the radiation condition $\partial_{\xi}a_{r}\approx ika_{r}$, we
obtain the spatiotemporal profile of the reflected radiation in spacetime
coordinates of the mirror rest frame as

${a_{r}(x,t)=a_{0}e^{-i(\omega
t-kx)}\frac{1}{2k}\int_{x-ct}^{x}k_{p}^{2}\left(s\right)e^{-i2ks}ds}.$ (9)

Eq. 9 satisfies the radiation boundary condition $\partial_{\xi}a_{r}=ika_{r}$
exactly for $t\rightarrow\infty$, since then $(\eta+\xi)/2=x$ and
$\xi=x-ct\rightarrow-\infty$, and only the plane wave phase keeps the $\xi$
dependence. This is equivalent to the assumption that, compared to the phase,
the amplitude of the wave described by the integral term varies slowly with
$\xi$. A more precise solution near could be obtained near the spacetime
origin by direct integration of Eq. 7, which leads to a more complicated form.
Since we are mostly interested in the far-field behavior, we only keep the
term given by Eq. 9 which well describes radiation behavior for large $t$ and
$x<ct$.

Finally, the spatiotemporal profile of the reflected wave is obtained by
taking the real part of Eq. 9, which is shown in Fig. 2 where a sine wave is
reflecting from a Gaussian electron profile. We see that the error of our far-
field approximation is small and leads to a small phase shift and reflected
amplitude reduction near the line $x=ct$, while the behavior at $x<ct$ shows
propagation of reflected radiation with amplitude corresponding to the exact
limit value. Precise integration without any approximations of Eq. 7 would
give us correct phase behavior near the spacetime origin. The absolute value
of Eq. 9, corresponding to the amplitude of the reflected wave, reduces in the
limit $x\rightarrow\infty$ to the Fourier transform of electron distribution
in the mirror rest frame, which was obtained in previous work on relativistic
mirrors [1], and the complex amplitude reflection coefficient $r=a_{r}/a_{0}$
in the far-field limit can be therefore expressed as

$r=\frac{1}{2k}\int_{-\infty}^{\infty}k_{p}^{2}\left(s\right)e^{-i2ks}ds,$
(10)

where $k=\sqrt{\frac{1+\beta}{1-\beta}}k_{0}$ is the wave number in the mirror
rest frame. The reflected amplitude in the laboratory frame is then simply
$E_{r}/(rE_{0})=(1+\beta)/(1-\beta)$, where $E_{0}$ is the electric field
amplitude of the incident laser in the laboratory frame.

In the following sections, we first analyze and identify properties of beam-
driven relativistic mirrors to consequently obtain the electron density
distribution required to calculate Eqs. 9-10.

##  2 Beam-driven nonlinear plasma wave equation

We proceed to analyze plasma waves driven by a charged particle beam in one-
dimensional approximation, which is valid for broad drivers satisfying
$k_{p}r_{b}\gg 1$, where $k_{p}=\omega_{p}/c$ is the linear plasma wave
number, $\omega_{p}=\sqrt{e^{2}n_{0}/\epsilon_{0}m_{e}}$ is the linear
electron plasma frequency, $n_{0}$ is the electron plasma density and $r_{b}$
is the radial size of the beam. The Poisson equation and the continuity
equation for plasma electrons are then given as

$\displaystyle\frac{\partial^{2}\phi}{\partial
x^{2}}=\frac{e}{\epsilon_{0}}(n_{e}-n_{0}-qn_{b}),\;$ (11)
$\displaystyle\frac{\partial j_{x}}{\partial x}-e\frac{\partial
n_{e}}{\partial t}=0,$ (12)

where $n_{b}$ is the number density of the driving beam and $q$ is the charge
number of the driving particles, e.g. it is positive for positrons and protons
($q=1$) or ions ($q\geq 1$) and it is negative for electrons ($q=-1$). We
assume that the driver density satisfies $n_{b}(x,t)=n_{b}(x-vt)$, which means
that it propagates in the $x$-direction and its shape does not evolve in time,
an assumption that is valid over a relatively long time for ultra-relativistic
beams with $\gamma\gg 1$, where $\gamma=1/\sqrt{1-(v/c)^{2}}$ is the
relativistic Lorentz factor of the driving beam and $c$ is the speed of light
in vacuum. We now assume that all the plasma quantities also depend only on
the Lagrangian coordinate of the driving beam $\xi=x-vt$. This gives us the
following transformations for the spatial and temporal derivatives,
respectively
$\partial_{x}f(\xi(x,t))=(\partial_{\xi}/\partial_{x})\partial_{\xi}f(\xi)=\partial_{\xi}f(\xi)$
and
$\partial_{t}f(\xi(x,t))=(\partial_{\xi}/\partial_{t})\partial_{\xi}f(\xi)=-v\partial_{\xi}f(\xi)$.
Using these transformations and introducing the following normalized units
$\omega_{p}\xi/c=\xi,\text{ }e\phi/m_{e}c^{2}\rightarrow\phi,\text{
}eE/m_{e}c\omega_{p}\rightarrow E,\text{ }n_{e}/n_{0}\rightarrow n_{e},\text{
}qn_{b}/n_{0}=\rho,\text{ }v/c=\beta,\text{ }j/en_{0}c\rightarrow j,\text{
}p/m_{e}c\rightarrow p,$ equations 11-12 become

$\displaystyle\frac{\partial^{2}\phi}{\partial\xi^{2}}=(n_{e}-1-\rho),\;$ (13)
$\displaystyle\frac{\partial}{\partial\xi}\left(n_{e}\beta-
n_{e}\beta_{e}\right)=0,$ (14)

where $\beta_{e}$ is the normalized velocity of the plasma electrons.
Integrating Eq. 14 with the initial conditions $n_{e}(\xi=0)=1\text{ and
}\beta(\xi=0)=0$, which are satisfied by plasma initially at rest, we obtain
$n_{e}=1/(1-\beta_{e}/\beta)$, which shows that extremely large density waves
are produced when $\beta_{e}\rightarrow\beta$, and Eq. 13 becomes

$\frac{\partial^{2}\phi}{\partial\xi^{2}}=\frac{\beta}{\beta-\beta_{e}}-1-\rho.$
(15)

It can be shown that the motion of an electron in a one-dimensional nonlinear
plasma wave is given by the normalized one-dimensional Hamiltonian,
$H(\xi,p_{x})=\gamma-\beta p_{x}-\phi(\xi)$. For plasma electrons which are
initially at rest, the following quantity is conserved, $\gamma-\beta
p_{x}=1+\phi$. Since only longitudinal fields are considered, the Lorentz
factor is simply $\gamma=\sqrt{1+p_{x}^{2}}$ and Eq. 15 can be cast as

$\frac{\partial^{2}\phi}{\partial\xi^{2}}=\gamma^{2}\left[\beta\left(1-\frac{1}{\gamma^{2}(1+\phi)^{2}}\right)^{-1/2}-1\right]-\rho,.$
(16)

Figure 2: Poisson equation (16) numerical solutions in the limit
$\gamma^{2}\gg\rho^{2}$. Wakefield-optimal solutions with bunch length given
by Eq. 25 for (a) $\rho=-0.25$ and (b) $\rho=1$. Single-cycle interior-wave-
optimal solution with bunch length with given by Eq. 24 for $\rho=1.5$.

To obtain analytical results, we approximate the driver as a flat-top beam,
$\rho(\xi)=\rho=const.$ for $-L_{b}<\xi<0$ and zero otherwise, where $L_{b}$
is the normalized driver length. We now proceed to analyze the plasma waves in
the region of the ultra-relativistic particle beam, where $-L_{b}<\xi<0$.
Multiplying Eq. 16 by $\phi^{\prime}$ and integrating with the initial
conditions $\phi(\xi=0)=\phi^{\prime}(\xi=0)=0$, we get the following equation
for the electric field of the plasma wave within the driving particle beam,

$\begin{split}\left(\frac{\partial\phi}{\partial\xi}\right)^{2}=&-2(\rho+\gamma^{2})\phi\\\
+&2\beta\gamma^{2}\left(\sqrt{(1+\phi)^{2}-\frac{1}{\gamma^{2}}}-\beta\right).\end{split}$
(17)

The existence of plasma oscillations requires the existence of electric field
stationary points, which is stated through Eq. 16 as $\phi^{\prime\prime}=0$.
From this condition we obtain that real solutions exist only if
$\rho>-1/(1+\beta)\approx-1/2$, which tells us that for dense, negatively-
charged drivers with $\rho\ll-1/2$, the background plasma cannot support
oscillations in the region of the driver due to the strong electrostatic field
of the charged beam. For $\rho>-1/(1+\beta)$, equations 16-17 allow us to find
the amplitude of the nonlinear plasma wave which exists in the interior region
of the driver, which we call interior plasma wave from now on. From the
electric field null-point condition $\phi^{\prime}=0$, we obtain the solutions
$\phi_{0}=0$ and
$\phi_{min}=-2\rho(1-1/\gamma^{2})/(2\rho+1+\rho^{2}/\gamma^{2})$, where the
zero and non-zero solutions correspond respectively to the interior electron
density wave trough and crest, respectively, in the case when $\rho>0$ and
vice versa for $\rho<0$.

##  3 Lorentz factor of plasma electrons and wave breaking threshold

Using the conserved Hamiltonian of plasma electrons initially at rest, the
electron Lorentz factor can be expressed explicitly in terms of the normalized
potential as
$\gamma_{e}=\gamma^{2}(1+\phi)\left[1-\beta\sqrt{1-1/(\gamma(1+\phi))^{2}}\right]$,
which gives us the Lorentz factor of the electrons in the density peak
(corresponding to $\phi_{min}$) of the nonlinear interior wave as

$\gamma_{e}=1+\frac{2\rho^{2}\beta^{2}}{2\rho+1+(\rho/\gamma)^{2}}.$ (18)

The wavebreaking limit of interior plasma waves is given by
$\gamma_{e}\rightarrow\gamma$. Solving Eq. 18 in this limit gives us the
following normalized driver density for which wavebreaking occurs immediately
upon beam entry into the plasma,

$\rho_{wb,\pm}=\begin{cases}\gamma,&\text{if $q>0$}.\\\
-\frac{\gamma}{2\gamma+1},&\text{if $q<0$}.\end{cases}$ (19)

The driver density corresponding to wavebreaking (19) has two solutions, which
correspond to positively and negatively charged drivers, respectively.

##  4 Wavelength of interior plasma waves

Inverting Eq. 17 and integrating from $\xi(\phi)<0$ to $\xi=0$, we get an
explicit dependence of the normalized Lagrangian coordinate $\xi$ on the
normalized potential and driver velocity as

$\begin{split}\xi(\phi,\gamma)&=\int_{\phi}^{0}\left[\vphantom{\int_{1}^{2}}-2(\rho+\gamma^{2})\phi\right.\\\
&\left.+2\beta\gamma^{2}\left(\sqrt{(1+\phi)^{2}-\frac{1}{\gamma^{2}}}-\beta\right)\vphantom{\int_{1}^{2}}\right]^{-1/2}d\phi.\end{split}$
(20)

Finally, to obtain the wavelength of the interior plasma waves, we integrate
Eq. 20 between the two null-points of the electric field given by $\phi_{0}=0$
and $\phi_{min}=-2\rho(1-1/\gamma^{2})/(2\rho+1+\rho^{2}/\gamma^{2})$, which
gives us the half-period of the interior wave. The wavelength of the nonlinear
interior wave normalized to the linear plasma wavelength is then
$\lambda_{i}/\lambda_{p}=-2\xi(\phi_{min})/2\pi$. Therefore

$\begin{split}\frac{\lambda_{i}}{\lambda_{p}}(\gamma,\rho)&=\frac{1}{\pi}\int_{\frac{-2\rho(1-1/\gamma^{2})}{(2\rho+1+\rho^{2}/\gamma^{2})}}^{0}\left[\vphantom{\int_{1}^{2}}-2(\rho+\gamma^{2})\phi^{\prime}\right.\\\
&\left.+2\beta\gamma^{2}\left(\sqrt{(1+\phi^{\prime})^{2}-\frac{1}{\gamma^{2}}}-\beta\right)\vphantom{\int_{1}^{2}}\right]^{-1/2}d\phi^{\prime},\end{split}$
(21)

which is valid when the condition for interior plasma wave existence is
satisfied, $\rho>-1/(1+\beta)$. Evaluating the integral 21 in the limit
$\gamma\gg 1$, we get the following expression for the nonlinear interior
plasma wave wavelength

$\frac{\lambda_{i}}{\lambda_{p}}=\frac{2}{\pi}\frac{E\left(\frac{2\rho}{2\rho+1}\right)}{\sqrt{2\rho+1}},$
(22)

where $E(k)$ is the complete elliptic integral of the second kind. For
negatively-charged drivers with $-1/2<\rho\leq 0$, the plasma wavelength grows
with increasing negative charge density and diverges at $\rho\rightarrow-1/2$.
For positively-charged drivers with $\rho\gg 1$, expansion of Eq. 22 yields
$\lambda_{i}/\lambda_{p}\approx(2/\rho)^{1/2}/\pi$, which reveals a unique
mechanism of relativistic nonlinear plasma waves with extremely short-
wavelengths, driven by positively-charged particle beams.

At the wavebreaking limit $\rho\rightarrow\gamma\gg 1$, the wavelength of the
nonlinear wave slightly decreases from Eq. 22, which can be estimated from Eq.
20 in this limit as
$\pi\lambda_{i}/\lambda_{p}\approx\left(2/\gamma\right)^{1/2}$. This stands in
stark contrast to plasma wakefields driven by either negatively-charged beams
or laser pulses, in which case the plasma wave elongates with increasing
driver velocity as
$\pi\lambda_{i}/\lambda_{p}\approx\left(2\gamma\right)^{1/2}$ [2, 3].

##  5 Optimal driver length

Now we proceed to analyze the plasma waves behind the driver, where
$\xi<-L_{b}$. The Poisson equation (16) with $\rho=0$ is satisfied with
initial conditions given by $\phi(-L_{b}),\phi^{\prime}(-L_{b})$, which
correspond to the solutions obtained in the driver region $-L_{b}\leq\xi\leq
0$ due to the continuity of the electric field. Integrating Eq. 16 for
$\rho=0$ with non-zero initial conditions, Eq. 17 takes the following form [4]

$\gamma_{e}+\frac{1}{2}\left(\frac{E_{wake}}{E_{0}}\right)^{2}=C_{0},$ (23)

where $C_{0}=\gamma_{e}(-L_{b})+(-\phi^{\prime}(-L_{b}))^{2}/2$ is a constant
given by initial conditions $\phi(-L_{b}),\phi^{\prime}(-L_{b})$ at the driver
rear $\xi=-L_{b}$.

When the driver length is exactly the length of a single cycle of the
nonlinear plasma wave, $L_{b}=\lambda_{i}$, the potential after a full-cycle
gives, as shown above, the initial conditions as
$\phi(-L_{b})=\phi^{\prime}(-L_{b})=0$, and we have $C_{0}=1$ and the maximum
amplitude of the wake is $E_{wake}/E_{0}=\sqrt{2(C_{0}-\gamma_{e})}=0$, since
$1\leq\gamma_{e}\leq C_{0}=1$. Therefore, a single-cycle plasma wave is driven
optimally and no plasma oscillations occur outside of the region occupied by
the ultra-relativistic particle beam. The same analysis holds for optimal
multi-cycle waves, in such a case the initial conditions for the Poisson
equation are the same and the driver length condition for optimal multi-cycle
nonlinear plasma wave generation becomes

$L_{b}=n\lambda_{i},$ (24)

where $n$ is a natural number and $\lambda$ is given by Eq. 22. Any deviation
from Eq. 24 gives $C_{0}>1$ and therefore non-zero plasma wakefield amplitude.
The wakefield amplitude is maximized when $\phi^{\prime}(-L_{b})=0$ and
$C_{0}=\gamma_{e}$, where $\gamma_{e}$ is given by Eq. 18. These boundary
conditions correspond to a driver whose length $L_{b}$ ends exactly at the
density peak of the interior plasma wave, as shown in Figs. 2a,b. The
strongest wakefield is therefore created when the driver length satisfies the
following resonance condition,

$L_{b}=\left(n-\frac{1}{2}\right)\lambda_{i},$ (25)

i.e. for the single-cycle wave $m=1$, the driver length fits the first half-
cycle of the interior plasma wave, $L_{b}=\lambda_{i}/2$, and the plasma
wakefield amplitude is maximized for given driver density and velocity, as
shown in Fig. 2a.

##  6 Mirror propagation distance due to relativistic beam-plasma instability

As a relativistic charged particle beam propagates through plasma, an
instability develops which modulates its initially uniform charge density.
This modulation eventually breaks the coherent nonlinear oscillations of
plasma electrons which sustain the relativistic mirror. To calculate the
distance in plasma a beam-driven relativistic mirror propagates before it
breaks, we need to find the instability growth rate.

Let us consider a relativistic charged particle beam with normalized charge
density $\rho_{b}=qn_{b}/n_{0}$, mass $M$ and velocity $v$ propagating along
the positive direction of the x-axis in plasma composed of electrons and ions
with particle number density $n_{0}$. However, since the beam is relativistic,
and its longitudinal size is of the order of the electron plasma wave
$L_{b}\propto\lambda_{p}$, the influence of slow ions is negligible compared
to the fast plasma electrons. The continuity equation and the Euler momentum
equation are then given for both relativistic particle beam and plasma
electron species as

$\displaystyle\frac{\partial n_{i}}{\partial
t}+\frac{\partial(n_{i}v_{i})}{\partial x}=0,$ (26) $\displaystyle
m\frac{\partial(\gamma_{i}v_{i})}{\partial
t}+mv_{i}\frac{\partial(\gamma_{i}v_{i})}{\partial x}=qE_{x},$ (27)

where $i=b,e$ is the species index and the plasma electrons are initially at
rest, $\gamma_{e}\approx 1$. Now, we consider the small-amplitude modes which
perturb the the relativistic beam, as $n_{i}\approx n_{0}+\delta n_{i0}$ and
$v_{i}\approx v_{i0}+\delta v_{i}$, where $\delta n_{i},\delta E_{x},\delta
v_{i}\propto e^{i(kx-\omega t)}$. Expanding the momentum in terms of the small
perturbation, we get
$p_{i}=\gamma_{i}v_{i}\approx\gamma_{i0}v_{i0}+\gamma_{i0}^{3}\delta v_{i}$.
Combining the linearized Eqs. 26,27 with the Poisson equation,
$\epsilon_{0}ik\delta E_{x}=\delta\rho_{b}+\delta\rho_{e}$, we get the
following dispersion relation for the relativistic two-stream instability

$\frac{\omega_{b}^{2}}{\gamma_{b}^{3}(\omega-
kv_{b})^{2}}+\frac{\omega_{p}^{2}}{\omega^{2}}=1,$ (28)

where $\omega_{p}^{2}=e^{2}n_{0}/(m_{e}\epsilon_{0})$ and
$\omega_{b}^{2}=(qe)^{2}n_{b}/(M\epsilon_{0})=(q\rho_{b}m/M)\omega_{p}^{2}$.
The $\gamma_{b}^{3}$ factor in the dispersion relation can be intuitively
understood as the well-known concept of the longitudinal mass of a particle
oscillating along the direction of its relativistic motion,
$M\rightarrow\gamma_{b}^{3}M$. Eq. 28 can be rewritten as

$\frac{1}{\omega-
kv_{b}}=\pm\frac{\gamma_{b}^{3/2}}{\omega_{b}}\sqrt{1-\left(\frac{\omega_{p}}{\omega}\right)^{2}}.$
(29)

Using $\omega=\omega_{r}+i\Gamma$, where $\omega_{r}=\mathcal{R}(\omega)$, and
expanding by $\omega_{r}-kv_{b}-i\Gamma$, we get

$\frac{\omega_{r}-kv_{b}-i\Gamma}{\omega_{r}^{2}-2\omega_{r}kv_{b}+(kv_{b})^{2}+\Gamma^{2}}=\pm\frac{\gamma_{b}^{3/2}}{\omega_{b}}\sqrt{1-\left(\frac{\omega_{p}}{\omega}\right)^{2}}.$
(30)

For even moderately relativistic beams satisfying $v_{b}\approx c$, unstable
copropagating electrostatic modes must be also relativistic, with constant
phase velocity given by $\omega_{r}\approx ck$. Subsequently, the growth rate
must be very small, $\omega=\omega_{r}+i\Gamma\approx\omega_{r}<\omega_{p}$,
which can be verified by solving Eq. 28 numerically, as shown in Fig. 3. Using
this in Eq. 30, we get the following estimate for the growth rate of the
relativistic non-resonant two-stream unstable modes

$\frac{\omega_{p}}{\Gamma}\approx\gamma_{b}^{3/2}\left(\frac{M}{qm_{e}\rho_{b}}\right)^{1/2}\frac{\omega_{p}}{kv_{b}}.$
(31)

Figure 3: Real (blue) and imaginary (red) parts of solution of Eq. 28,
calculcated for a positively-charged relativistic particle beam with
$qn_{b}/n_{0}=\rho_{wb,+}=\gamma$, $M=m_{e}$ and $\gamma=20$ (positron). The
real part corresponding to the growing unstable modes is given by linearly
growing branch, $\omega_{r}\approx kv_{b}\approx kc$.

Finally, the unstable spectrum observed in PIC simulations is broadband.
Typically, the modulation occurs only within a very small region surrounding
the copropagating electron density peak of the nonlinear wave (Fig. 2c), where
the electrons pile-up and modulate the driving particle beam. This makes
single mode interpretation complicated.

To proceed, we make a generalizing heuristic estimate. Requiring that the
mirror breaks immediately, $\tau_{wb}\rightarrow 0$, when $Q\geq Q_{wb}$,
where $Q$ the total charge of the relativistic beam and $Q_{wb}$ is the total
charge required to produce wave breaking, we get the following estimate that
correctly reproduces the timescale ($\tau_{wb}\propto 1/\Gamma$) of the
unstable broadband spectra observed in the PIC simulations, written as
$\omega_{p}/kv_{b}\rightarrow\sqrt{1-Q/Q_{wb}}$. This yields the following
expression for the beam-driven relativistic mirror propagation distance

$k_{p}L_{wb}\approx
N_{e}\gamma_{b}^{3/2}\left(\frac{M}{qm_{e}\rho_{b}}\right)^{1/2}\left(1-\frac{Q}{Q_{wb}}\right)^{1/2},$
(32)

where the factor $N_{e}$ corresponds to the number of e-folds observed in PIC
simulations before the mirror breaks. Typically, $N_{e}\approx 1$. For
positively-charged drivers with $\rho_{b}\lesssim\rho_{wb,+}\approx\gamma$, we
have $Q_{b}/Q_{wb}\approx\rho_{b}/\rho_{wb,+}$ and the mirror propagation
distance, for $N_{e}\approx 2\pi/5$, takes the form presented in the main text

$\frac{L_{wb,+}}{\lambda_{p}}\approx\frac{\gamma_{b}}{5}\left(\frac{M}{qm_{e}}\right)^{1/2}\left(1-\frac{\rho_{b}}{\rho_{wb}}\right)^{1/2}.$
(33)

For negatively-charged drivers with
$\rho_{b}\approx\rho_{wb,-}=-1/(2+(1/\gamma))\approx-1/2$ and
$Q_{b}/Q_{wb}\approx L_{b}/L_{wb}<1$, the mirror propagation distance can be
estimated as

$\frac{L_{wb,-}}{\lambda_{p}}\approx\frac{\gamma_{b}^{3/2}}{5}\left(\frac{2M}{|q|m_{e}}\right)^{1/2}\left(1-\frac{L_{b}}{L_{wb}}\right)^{1/2}.$
(34)

##  7 Structure of beam-driven interior plasma wave breaking

As we have shown in S1, evaluating reflection of radiation requires electron
plasma density distribution of the breaking wave. To gain a foothold, we
return to the Poisson equation (15) to analyze the structure of interior
plasma waves at the wavebreaking threshold, $\gamma_{e}\rightarrow\gamma$. We
use the integral of motion $\gamma_{e}-\beta p_{x}-1=\phi$ and the identity
$1/(1-\beta_{e}/\beta)-1=p_{x}/(\beta\gamma_{e}-p_{x})$ to rewrite Eq. 15 as

$(\gamma_{e}-\beta
p_{x})^{\prime\prime}=\frac{p_{x}}{\beta\gamma_{e}-p_{x}}-\rho.$ (35)

We now expand the momentum around the peak of the nonlinear plasma wave,
$p=p_{m}+\delta p+\mathcal{O}(\delta p^{2})$, where
$p_{m}=\sqrt{\gamma_{e}^{2}-1}$, where $\gamma_{e}$ is given by Eq. 18. At the
wavebreaking threshold we have $\gamma_{e}\rightarrow\gamma$, therefore we
proceed to use the following power series expansion for left-hand side of Eq.
35, $\gamma-\beta p_{x}=\sqrt{1+p_{x}^{2}}-\beta p_{x}\approx 1/\gamma+\delta
p^{2}/(2\gamma^{3})+\mathcal{O}(\delta p^{3})$. For the right-hand side of Eq.
35 we use $\gamma_{e}\rightarrow\gamma$ and therefore $\beta\gamma-
p_{x}=\beta\sqrt{1+p_{x}^{2}}-p_{x}\approx-\delta p/\gamma^{2}+\beta\delta
p^{2}/(2\gamma^{3})+\mathcal{O}(\delta p^{3})$. Keeping only the leading order
terms in Eq. 35, we get

$(\delta p^{2})^{\prime\prime}=-2\gamma^{6}\frac{\beta}{\delta p}-\gamma.$
(36)

Multiplying Eq. 36 by $(\delta p^{2})^{\prime}$ and integrating over $\xi$, we
get

$(\delta p^{\prime}\delta p)^{2}+2\gamma^{6}\beta\delta p+\gamma\delta
p^{2}/2=C,$ (37)

where $C$ is an integration constant. In the limit $\delta p\rightarrow 0$,
the solution of Eq. 37 is given by the value of the integration constant $C$.
If $\delta p^{\prime}\delta p\rightarrow 0$ when $\delta p\rightarrow 0$, then
we have $C=0$ and up to the leading order Eq. 37 takes the following form

$(\delta p^{\prime})^{2}=-2\gamma^{6}\frac{\beta}{\delta p}.$ (38)

In the limit $\xi\rightarrow 0$, where $\xi$ is the distance from the position
of the density peak, a solution of Eq. 37 is given as

$\delta p=-\left(\frac{3}{\sqrt{2}}\gamma^{3}\beta^{1/2}\xi\right)^{2/3}.$
(39)

The velocity can be expanded up to the first order as
$\beta_{e}\approx\beta+\frac{\delta p}{\gamma^{3}}+\mathcal{O}(\delta p^{2})$,
therefore
$\beta_{e}\approx\beta-\frac{1}{\gamma}\left(\frac{3}{\sqrt{2}}\beta^{1/2}\xi\right)^{2/3}$
and the normalized electron density becomes

$\frac{n_{e}}{n_{0}}=\frac{\beta}{\beta-\beta_{e}}\approx\gamma\left(\frac{\sqrt{2}}{3}\frac{\beta}{\xi}\right)^{2/3}.$
(40)

Eq. 40 shows that the density structure of a breaking interior plasma wave has
the form of a cusp caustic [5], which has an integrable singularity at
$\xi\rightarrow 0$, therefore the total number of electrons remains finite
even for the asymptotic leading-order term given by Eq. 40. It was previously
shown that such a cusp singularity is a general feature of electron plasma
wave breaking [6, 1, 7]. The structure of beam-driven wave breaking is not
different, since the effect of driver density represented by the third term of
Eq. 37 only contributes in the higher order corrections.

##  8 Reflection coefficient of beam-driven relativistic mirror at the wave
breaking threshold

In the laboratory frame, the electron distribution given by the asymptotic
cusp distribution (40) can be written as

$\frac{n_{e}}{n_{0}}=\gamma\left(\frac{\sqrt{2}}{3}\frac{\beta}{k_{p}(x-vt)}\right)^{2/3}.$
(41)

We now transforming into the mirror rest frame using
$x^{\prime}=\gamma(x-vt)$. The electron charge density is part of the four-
current which satisfies the continuity equation $\partial_{\mu}j^{\mu}=0$, and
it must therefore transform as $\rho^{\prime}=\gamma(\rho-j_{x}\beta/c)$ due
to charge conservation. Transforming the spacetime coordinates, the rest frame
electron density distribution then becomes
$n_{e}^{\prime}/n_{0}^{\prime}=\gamma(\sqrt{2}\beta\gamma/3)^{2/3}(x^{\prime})\text{}^{-2/3}$,
where $n_{0}^{\prime}=n_{0}/\gamma$ is the background plasma density in the
mirror rest frame. We therefore have the following inhomogenous relativistic
plasma wave number (using $\beta\approx 1$)

$k_{p}^{\prime
2}(x^{\prime})=\left(\frac{2}{9}\right)^{1/3}\frac{k_{p}^{4/3}\gamma^{2/3}}{x^{\prime
2/3}},$ (42)

where $k_{p}=\omega_{p0}/c$ is the linear plasma wave number calculated in the
laboratory frame. Using Eq. 42 to calculate the amplitude reflection
coefficient given by Eq. 10, we get

$r=\kappa\left(\frac{n_{0}}{\gamma n_{c}}\right)^{2/3},$ (43)

where $\kappa=\pi/(2^{4/3}3^{2/3}\Gamma\left(\frac{2}{3}\right))\approx 0.44$
and $\Gamma(z)$ is the Euler gamma function. Complex amplitude reflection
coefficient (43) was previously obtained in the context of laser-driven
relativistic mirrors [1, 6]. Calculating Eq. 9 in the limit
$t\rightarrow\infty$, we can further obtain the spatial profile of the
reflected amplitude around the mirror position, which is shown in Fig. 4.

Figure 4: Spatial profile of complex amplitude reflection coefficient near the
origin of the mirror rest frame ($x=0$). Here $\lambda$ represents the rest
frame radiation wavelength. The reflected wave amplitude is modulated at the
second harmonic of the incident radiation near the mirror origin. In the limit
$x\rightarrow\infty$, the amplitude stabilizes at the value given by Eq. 43.

In PIC simulations, we observe that Eq. 43 overestimates amplitude of the
reflected radiation. Indeed, Eq. 43 is derived using the singular cusp
distribution. However, in fully kinetic simulations electron density is finite
both before, during and after wavebreaking, which is the expected behavior due
to the discrete number of electrons (and macroparticles in the case of
particle-in-cell simulations), which compose nonlinear plasma waves.

To amend this, we introduce the finite distribution scale length $l$, also
known as sharpness parameter [1], and generalize the electron distribution in
the laboratory frame, given by Eq. 41, as

$\frac{n_{e}}{n_{0}}=\frac{\gamma
2^{1/3}}{3^{2/3}k_{p}^{2/3}}\left(\frac{\beta^{2}}{l^{2}+(x-vt)^{2}}\right)^{1/3},$
(44)

where the limit $l\rightarrow 0$ reproduces the dominant asymptotic cusp
structure obtained by solving Eq. 35. Transforming Eq. 44 into the rest frame,
we get the following rest frame inhomogeneous plasma wave number

$k_{p}^{\prime
2}(x^{\prime})=\left(\frac{2}{9}\right)^{1/3}\frac{k_{p}^{4/3}\gamma^{2/3}}{l^{\prime
2}+x^{\prime 2}},$ (45)

where $l^{\prime}=\gamma l$ is the finite distribution length scale in the
mirror rest frame. The rest frame reflection coefficient evaluated in terms of
laboratory frame parameters then becomes

$r\approx\frac{\pi^{1/2}}{2^{5/6}3^{2/3}\Gamma\left(\frac{1}{3}\right)}\frac{k_{p}^{4/3}l^{1/6}}{\gamma^{1/3}k_{0}^{7/6}}K_{1/6}\left(4\gamma^{2}k_{0}l\right),$
(46)

where $K_{n}(x)$ is the modified Bessel function of the second kind. By
measuring the finite distribution scale length $l$ from simulations, the
reflection coefficient is precisely given by Eq. 46. Fully-coherent emission
will occur when a distribution of radiating particles satisfies $k_{r}l\ll 1$,
where $k_{r}$ is the wave number of the reflected radiation in the laboratory
frame, $k_{r}\approx 4\gamma^{2}k_{0}$.

In the PIC simulations, we observe finite distribution length scales of the
order $k_{r}l_{sim}\approx 3.5\times 10^{-2}$, therefore we estimate
$4\gamma^{2}k_{0}l\approx 1/(9\pi)$ and the rest frame amplitude reflection
coefficient can be written in the following simple form

$r\approx\alpha\left(\frac{n_{0}}{n_{c}\gamma}\right)^{2/3},$ (47)

where
$\alpha=\frac{\pi^{1/3}}{2^{1/6}3^{2}\Gamma(1/3)}K_{1/6}\left(\frac{1}{3^{2}\pi}\right)\approx
0.3$ is the proportionality constant obtained by introducing finite
distribution length scale $l$. The scaling follows that of the cusp reflection
coefficient given by Eq. 43, but the coefficient is roughly $30\%$ smaller due
to the electron density finiteness. As shown with the results presented in the
main text, Eq. 47 perfectly describes the particle-in-cell simulation results
of relativistic mirrors close to the wave breaking threshold in the parameter
range $\gamma\in\\{5,10,15,20\\}$ and
$n_{0}/n_{c}\in\\{0.01,0.02,0.04,0.08\\}$.

##  9 Energy spectrum of reflected coherent radiation

If we assume that the radiation is reflected with constant amplitude
reflection coefficient, the electric field in the spectral domain of a signal
with duration $\tau_{r}$ can be expressed in the laboratory frame as

$\displaystyle\begin{split}E_{r}(x,\omega)=&E_{r}\int_{-\infty}^{\infty}e^{-i(\omega_{r}t-k_{r}x)}e^{i\omega
t}dt\\\
=&E_{r}e^{ik_{r}x}\int_{0}^{\tau_{r}}e^{i(\omega-\omega_{r})t}dt,\end{split}$
(48)

where $E_{r}=rE_{0}(1+\beta)/(1-\beta)$ is the reflected amplitude and $E_{0}$
is the incident electric field amplitude. Therefore, after directly evaluating
the last integral in Eq. 48, the normalized energy spectrum can be expressed
as

$\left|\frac{E(\omega)}{E_{p}}\right|^{2}=\text{sinc}^{2}\left(\frac{(\omega-\omega_{r})\tau_{r}}{2}\right),$
(49)

where $\text{sinc}(x)=\text{sin}(x)/x$ is the cardinal sine function,
$E_{p}^{2}=\left|E_{r}\right|^{2}\tau_{r}^{2}$ is the value of the spectrum
peak at $\omega=\omega_{r}$, which grows quadratically in time due to temporal
coherence, and $\tau_{r}$, $\omega_{r}$ are respectively the double Doppler
shifted reflected angular frequency and pulse duration.

##  10 Peak spectral brightness

The value of peak spectral brightness for a monochromatic light source
reflected from a homogeneous mirror in terms of the conventional unit of
photons per second per relative spectral bandwidth of $10^{-3}$ per transverse
radiation phase space can be written as

$\mathcal{B}=\frac{N_{ph,r}(N/1000)}{2\pi^{2}\theta^{2}\sigma^{2}\tau_{r}},$
(50)

where $N$ is the number of reflected cycles corresponding to the relative
spectral bandwidth as $\Delta\omega/\omega=1/N$, $N_{ph,r}$ is the total
number of reflected photons, $\theta$ is the half-angle divergence of
reflected radiation and $\rho$ is the waist of the incident laser. The
reflected divergence is $\theta\approx\theta_{0}/(2\gamma)^{2}$, where
$\theta_{0}$ is the half-angle divergence of the incident laser. If we
consider a diffraction limited incident laser and a homogeneous mirror, the
beam parameter product characterizing transverse phase space of incident
radiation becomes $\rho\theta_{0}=\lambda_{0}/\pi$, where $\lambda_{0}$ is the
incident laser wavelength. The peak spectral brightness is then

$\mathcal{B}=\frac{8\gamma^{4}N_{ph,r}(N/1000)}{\lambda_{0}^{2}\tau_{r}}.$
(51)

Finally, considering the reflected pulse duration satisfies
$\tau_{r}\approx\tau_{0}/(2\gamma)^{2}\approx(N\lambda_{0}/c)/(2\gamma)^{2}$,
we obtain the final formula for the peak spectral brightness of radiation
reflected from a homogeneous beam-driven relativistic mirror,

$\mathcal{B}=\frac{32c\gamma^{6}N_{ph,r}}{1000\lambda_{0}^{3}}.$ (52)

The ratio of reflected to incident number of photons can be calculated using
the rest frame reflection coefficient as
$N_{ph,r}/N_{ph,i}=\mathinner{\\!\left\lvert r\right\rvert}^{2}$. Finally, to
convert to the conventional synchrotron unit of $\text{photons/(s
mm}^{2}\text{ mrad}^{2}\text{ }0.1\%\text{BW)}$, we use
m2rad${}^{2}=10^{12}$mm2mrad2 to obtain

$\mathcal{B}\left[\frac{\text{photons}}{\text{s mm${}^{2}$ mrad${}^{2}$
0.1$\%$ BW}}\right]=\frac{32c\gamma^{6}N_{ph,i}\mathinner{\\!\left\lvert
r\right\rvert}^{2}}{10^{15}\lambda_{0}^{3}},$ (53)

where the reflection coefficient is given by Eq. 47.

## References

  *   * Bulanov et al. [2013] Bulanov, S.V., Esirkepov, T.Z., Kando, M., Pirozhkov, A.S., Rosanov, N.N.: Relativistic mirrors in plasmas. novel results and perspectives. Physics-Uspekhi 56(5), 429 (2013)
  * Esarey et al. [2009] Esarey, E., Schroeder, C.B., Leemans, W.P.: Physics of laser-driven plasma-based electron accelerators. Reviews of modern physics 81(3), 1229 (2009)
  * Rosenzweig [1987] Rosenzweig, J.: Nonlinear plasma dynanics in the plasma wakefield accelerator. IEEE transactions on plasma science 15(2), 186–191 (1987)
  * Akhiezer and Polovin [1956] Akhiezer, A.I., Polovin, R.: Theory of wave motion of an electron plasma. Soviet Phys. JETP 3 (1956)
  * [5] Arnold, V.I.: Catastrophe Theory. Berlin Springer-Verlag, 1986
  * Panchenko et al. [2008] Panchenko, A., Esirkepov, T.Z., Pirozhkov, A., Kando, M., Kamenets, F., Bulanov, S.: Interaction of electromagnetic waves with caustics in plasma flows. Physical Review E 78(5), 056402 (2008)
  * Bulanov et al. [1998] Bulanov, S., Naumova, N., Pegoraro, F., Sakai, J.: Particle injection into the wave acceleration phase due to nonlinear wake wave breaking. Physical Review E 58(5), 5257 (1998)