# QCD in the cores of neutron stars††thanks: Presented at Quark Matter 2022

Oleg Komoltsev  Faculty of Science and Technology, University of Stavanger,
4036 Stavanger, Norway

###### Abstract

I discuss why state-of-the art perturbative QCD calculations of the equation
of state at large chemical potential that are reliable at asymptotically high
densities constrain the same equation of state at neutron-star densities. I
describe how these theoretical calculations affect the EOS at lower density. I
argue that the ab-initio calculations in QCD offer significant information
about the equation of state of the neutron-star matter, which is complementary
to the current astrophysical observations.

##  1 Introduction

The equation of state (EOS) of the dense matter at zero temperature is a
necessary input for the neutron-stars (NS) physics. Theoretical calculations
of the EOS can be done only at the two opposite (low- and high-density)
limits. At the low-density limit the matter can be described within the chiral
effective field theory (CET) [1, 2]. Those calculations are reliable up to
around nuclear saturation density $n_{s}=0.16/\textrm{fm}^{3}$. On the other
hand we can access the EOS using perturbative Quantum Chromodynamics (pQCD) at
the asymptotically high densities, above $\sim 40n_{s}$ [3, 4]. Central
densities of maximally massive neutron stars are around $4-8n_{s}$, which is
not reachable within CET or pQCD. Therefore, there are no tools in our
possession to compute EOS of the cores of NS from the first principles.

However, we can obtain an empirical access to the cores of NSs using recent
astrophysical observations. The most important probes of NS physics are the
discovery of massive NSs [5, 6, 7], mass - radius measurements [8, 9], and the
gravitational-wave and multi-messenger astronomy [10, 11]. Utilizing all
constraints coming from astrophysical observation as well as first principle
calculations narrows down dramatically the range of possible EOSs, which
allows us to use the densest objects in the Universe to test independently
various beyond standard model scenarios and/or general relativity.

Majority of the EOS studies extrapolate CET EOS up to NS densities 5-10$n_{s}$
and conditioning it with the observational inputs. The results differ from the
works that include high-density limit and interpolate between two orders of
magnitude. The qualitative difference is in the softening of the EOS happening
around $\epsilon\sim$750 MeV/fm-3, which can be interpreted as quark matter
cores inside the most massive NS [12].

In this work I answer why and how pQCD input offers significant information
about the EOS at NS densities. I find that pQCD input propagates non-trivial
constraints all the way down to 2.2$n_{s}$ just by using solely thermodynamic
stability, consistency and causality [13]. In addition the complementariness
of the pQCD input to the astrophysical observations was studied in [14]. I
show that pQCD is responsible for the softening of the EOS at the NS
densities. Therefore, it is essential to include pQCD input in any inference
study of the EOS.

##  2 Setup

All technical details as well as analytical formulas are presented in [13]. In
this section I describe the conditions I use, in particular stability,
consistency and causality, and the resulting propagation of the pQCD input
down to lower densities. Let us start with the baryon density $n$ as a
function of the chemical potential $\mu$ as shown in fig.1(a). The goal is to
find all possible lines that connect endpoint of CET results (dark blue line
in the bottom left corner) with the first point of pQCD calculations (purple
line in the upper right corner) using 3 conditions.

(a)

(b)

Figure 1: (a) - Baryon density as a function of chemical potential.
Simultaneous fulfilment of thermodynamic consistency, stability and causality
narrows down the allowed region (white area) of the EOS. (b) - Zero
temperature EOSs from CompOSE database plotted with the allowed area (gray
shape) arising from the new constraints in $\epsilon-p-n$ space. Consistent/in
tension/not consistent means that EOS is consistent with integral constraints
for all/some/none X values in a range [1,4].

The first condition is thermodynamic stability, which implies concavity of the
grand canonical potential $\partial^{2}_{\mu}\Omega(\mu)\leq 0$. At zero
temperature $\Omega(\mu)=-p(\mu)$, which implies that the number density is
monotonically increasing function of the chemical potential
$\partial_{\mu}n(\mu)\geq 0$.

The second condition is causality – the sound speed cannot exceed the speed of
light $c^{2}_{s}\leq 1$. This provides constraints on the first derivative of
the number density with respect to the chemical potential

$c^{-2}_{s}=\frac{\mu}{n}\frac{\partial n}{\partial\mu}\leq 1.$ (1)

For each point on the $\mu-n$ plane we can calculate the least allowed slope
coming from causality, which is represented by the arrows in the fig.1(a).
This cuts upper (lower) region of the plane, because any points from the area
above (below) orange line $c^{2}_{s}=1$ cannot be connected to pQCD (CET) in a
casual way.

The third condition is thermodynamic consistency. In addition to $n$ and $\mu$
we need to match pressure $p$ at the low- and high- density limits. The
pressure is giving by the integral of the number density

$\int^{\mu_{\rm QCD}}_{\mu_{\rm CET}}n(\mu)d\mu=p_{\rm QCD}-p_{\rm CET}=\Delta
p.$ (2)

This implies that the area under the curve for any EOS is fixed by our input
parameters. For each arbitrary point ${\mu_{0},n_{0}}$ we can construct the
EOS that maximize/minimize the area under the curve $\Delta
p_{max/min}(\mu_{0},n_{0})$ shown as a green/blue dashed line in the fig.1(a).
If $\Delta p_{max}(\mu_{0},n_{0})<\Delta p$ then any EOS that goes through the
point ${\mu_{0},n_{0}}$ does not have enough area under the curve. This
discards the region in the lower right corner in the fig.1(a) under the red
line called ”integral constraints”. If $\Delta p_{min}(\mu_{0},n_{0})>\Delta
p$ then any EOS that goes through the point ${\mu_{0},n_{0}}$ has too much
area under the curve. This cuts area in the upper left corner above the red
line. The integral constraints can be obtained without any assumptions of
interpolation function in a completely general and analytical way.

Figure 2: Constraints on the $\epsilon-p$ plane coming from low- and high-
density limits. Shapes outlined by solid black line are the allowed areas for
fixed number density without pQCD. Blue shapes are allowed regions after
imposing pQCD input.

We can map the allowed region from $\mu-n$ to $\epsilon-p$ plane. The results
of such mapping is shown in the fig.2. The green envelope corresponds to the
the white area in the fig.1(a) restricted by the causality and the integral
constraints. The shapes of allowed region with and without pQCD are shown for
the fixed number density $n$ = 2,3,5 and 10$n_{s}$. This explicitly shows how
pQCD input can propagate information down to lower density starting from
2.2$n_{s}$. And, strikingly, at 5$n_{s}$ it excludes 75% of otherwise allowed
area.

Using the new constraints we can check the consistency of publicly available
EOSs. Results for all zero temperature EoSs in $\beta$-equilibrium from the
public CompOSE database [15, 16] are shown in the fig.1(b). Almost all of the
EOSs start to be inconsistent with pQCD input at some density within the
provided range.

##  3 Bayesian inference of EOS

With the construction described above we can propagate information from ab-
initio QCD calculations down to NS densities, where we already have
constraints from astrophysical observations. To understand if the new
constraints from pQCD go beyond the constraints coming from the NS
measurements we construct a Bayesian-inference framework. This was done in
[4], where we generate a large ensemble of different EOSs using Gaussian-
process regression. We anchor the ensemble to CET calculations and extrapolate
it up to 10$n_{s}$, where we impose pQCD input as a blue shape from fig.2. We
condition ensemble sequentially with the astrophysical observations. With this
setup we can turn on and turn off pQCD input in order to study its effect on
our posterior after imposing astrophysical observation.

The results are present in fig.3. The reduction of the pressure (green arrow
on the right plot), which is caused by the QCD input, happens before the
density reaches its maximal central value. In another words, the prediction of
QCD input is the softening of the EOS that happens inside the most massive
neutron stars.

Figure 3: Left plot shows the sample of 10k EOSs. The coloring represents the
likelihood after imposing all observations as well as pQCD input. Right plot
shows 67 %-credible intervals conditioned with the different astrophysical
observations and high-density limit. The gray band shows 67 %-credible
interval for the maximal central energies density reached in NSs.

##  4 Conclusion

In this work, I show how QCD calculations at asymptotically high densities can
propagate information down to lower densities using solely thermodynamic
consistency, stability and causality. This information offers a significant
constraints to the EOS at NS density, which is complementary to the current
astrophysical observations. In addition, I show that the prediction of QCD
input is the softening of the EOS that happens in the most massive NSs. An
easy-to-use python script is provided to check consistency of the EOS with
pQCD input, available on Github [17].

In order to achieve accurate determination of the EOS it is crucial to utilize
all available controlled measurements and theoretical calculations. This
strategy either helps us to understand the matter of the densest objects in
the Universe or find a discrepancy between different inputs, which allows us
to use NS as a tool for fundamental discoveries.

## References

  * [1] I. Tews, T. Krüger, K. Hebeler, and A. Schwenk, “Neutron matter at next-to-next-to-next-to-leading order in chiral effective field theory,” _Phys. Rev. Lett._ , vol. 110, no. 3, p. 032504, 2013.
  * [2] C. Drischler, K. Hebeler, and A. Schwenk, “Chiral interactions up to next-to-next-to-next-to-leading order and nuclear saturation,” _Phys. Rev. Lett._ , vol. 122, no. 4, p. 042501, 2019.
  * [3] T. Gorda, A. Kurkela, P. Romatschke, M. Säppi, and A. Vuorinen, “Next-to-next-to-next-to-leading order pressure of cold quark matter: Leading logarithm,” _Phys. Rev. Lett._ , vol. 121, no. 20, p. 202701, 2018\.
  * [4] T. Gorda, A. Kurkela, R. Paatelainen, S. Säppi, and A. Vuorinen, “Soft Interactions in Cold Quark Matter,” _Phys. Rev. Lett._ , vol. 127, no. 16, p. 162003, 2021.
  * [5] P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E. Roberts, and J. H. T. Hessels, “A two-solar-mass neutron star measured using Shapiro delay,” _Nature_ , vol. 467, pp. 1081–1083, 2010.
  * [6] J. Antoniadis _et al._ , “A massive pulsar in a compact relativistic binary,” _Science_ , vol. 340, no. 6131, p. 1233232, 2013.
  * [7] E. Fonseca _et al._ , “Refined Mass and Geometric Measurements of the High-mass PSR J0740+6620,” _Astrophys. J. Lett._ , vol. 915, no. 1, p. L12, 2021.
  * [8] M. C. Miller _et al._ , “The radius of PSR J0740+6620 from NICER and XMM-Newton data,” _Astrophys. J. Lett._ , vol. 918, no. 2, p. L28, 2021\.
  * [9] T. E. Riley _et al._ , “A NICER view of the massive pulsar PSR J0740+6620 informed by radio timing and XMM-Newton spectroscopy,” _Astrophys. J. Lett._ , vol. 918, no. 2, p. L27, 2021.
  * [10] B. P. Abbott _et al._ , “GW170817: Observation of gravitational waves from a binary neutron star inspiral,” _Phys. Rev. Lett._ , vol. 119, no. 16, p. 161101, 2017.
  * [11] ——, “Multi-messenger observations of a binary neutron star merger,” _Astrophys. J. Lett._ , vol. 848, no. 2, p. L12, 2017.
  * [12] E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, “Evidence for quark-matter cores in massive neutron stars,” _Nat. Phys._ , vol. 16, no. 9, pp. 907–910, 2020.
  * [13] O. Komoltsev and A. Kurkela, “How Perturbative QCD Constrains the Equation of State at Neutron-Star Densities,” _Phys. Rev. Lett._ , vol. 128, no. 20, p. 202701, 2022.
  * [14] T. Gorda, O. Komoltsev, and A. Kurkela, “Ab-initio QCD calculations impact the inference of the neutron-star-matter equation of state,” 4 2022.
  * [15] S. Typel, M. Oertel, and T. Klähn, “CompOSE CompStar online supernova equations of state harmonising the concert of nuclear physics and astrophysics compose.obspm.fr,” _Phys. Part. Nucl._ , vol. 46, no. 4, pp. 633–664, 2015.
  * [16] M. Oertel, M. Hempel, T. Klähn, and S. Typel, “Equations of state for supernovae and compact stars,” _Rev. Mod. Phys._ , vol. 89, no. 1, p. 015007, 2017.
  * [17] “https://github.com/OKomoltsev/QCD-likelihood-function.”


where $\tilde{c}$ denotes the $|\mathcal{Z}|$ components of $c$ corresponding to the columns of $A$ put into in $\tilde{A}$, and $\tilde{c}_c$ are the remaining entries $c_g$ of $c$. Then $\alpha = {\tilde{A}}^{'-1}\tilde{c}$, which can be seen by left-multiplying the above equation by the $|\mathcal{Z}| \times |\mathcal{G}|$ matrix $[{\tilde{A}}^{'-1},\mathbf{0}^{|\mathcal{Z}| \times |\mathcal{G}|-|\mathcal{Z}|}]$. Intuitively, the system $A'\alpha = c$ is over-determined, so we only only need the components $\tilde{c}$ of $c$ to uniquely determine the vector $\alpha$.

Now consider the case in which $|\mathcal{G}| < |\mathcal{Z}|$, so that the system $A'\alpha = c$ is now undetermined. Suppose for now that the rank of $A$ is $|\mathcal{G}|$ so that it has full column rank. One solution $\alpha$ can then be obtained by  writing $A=\begin{bmatrix}
\tilde{A} \\ {\tilde{A}_c}
\end{bmatrix} $ where $\tilde{A}$ is an invertible $|\mathcal{G}| \times |\mathcal{G}|$ matrix representing $|\mathcal{G}|$ linearly independent rows of $A$. Now consider $\alpha = \begin{pmatrix} \tilde{A}^{-1}c\\ \mathbf{0}^{(|\mathcal{Z}|-|\mathcal{G}|) \times 1} \end{pmatrix}$ where note that $\tilde{A}^{-1}c$ is $|\mathcal{G}|-$component vector. This represents one solution to $A'\alpha = c$ because

$$A'\begin{pmatrix} \tilde{A}^{'-1}c\\ \mathbf{0}^{(|\mathcal{Z}|-|\mathcal{G}|) \times 1} \end{pmatrix} = [\tilde{A}, \tilde{A}_c] \begin{pmatrix} \tilde{A}^{-1}c\\ \mathbf{0}^{(|\mathcal{Z}|-|\mathcal{G}|) \times 1} \end{pmatrix} = c $$

We can combine the constructions in the two special cases considered above to relax any assumptions about the cardinality of $\mathcal{Z}$ and $\mathcal{G}$ or the rank of $A$. Let the rank of $A$ be $k \le \min\{|\mathcal{Z}|, |\mathcal{G}|\}$. Write $A = A_k[I_k,M]$ where $A_k$ is a $k \times |\mathcal{G}|$ matrix composed of $k$ linearly independent columns of $A$, and $M$ is $(|\mathcal{G}|-k) \times k$ matrix that expresses the remaining $(|\mathcal{G}|-k)$ columns of $A$ as linear combinations of the columns of $A$ represented in $A_k$. Write $c = \begin{pmatrix} \tilde{c}_k\\ \tilde{c}_c \end{pmatrix}$ where $\tilde{c}_k$ collects the corresponding $k$ components of $c$. Note that if $c' = \alpha'A$ has a solution, then $c' = \tilde{c}_k'[I_k,M]$, since $c' = (\alpha_k' A_k)[I,M]$ where the $k$ components of $c'$ corresponding to the columns in $A_k$ are $\alpha_k' A_k$, so $\tilde{c}_k'=\alpha_k' A_k$. Now split the rows of $A_k$ as
\tilde{A} \\ {\tilde{A}_c}
\end{bmatrix}$ where $\tilde{A}$ is a square invertible $k \times k$ matrix representing $k$ linearly independent rows of $A_k$ and $\tilde{A}_c$ is $(|\mathcal{Z}|-k) \times k$. Now $\alpha = \begin{pmatrix} \tilde{c}_k'\tilde{A}^{-1}\\ \mathbf{0}^{(|\mathcal{Z}|-k)\times 1} \end{pmatrix}$ represents a solution to $c'=\alpha'A$ because

$$[\tilde{c}_k' \tilde{A}^{-1},  \mathbf{0}^{1 \times (|\mathcal{Z}|-k)}] A = [\tilde{c}_k' \tilde{A}^{-1},  \mathbf{0}^{1 \times (|\mathcal{Z}|-k)}] \begin{bmatrix}
\tilde{A} \\ {\tilde{A}_c}
\end{bmatrix}[I_k,M] = \tilde{c}_k'[I_k,M] = c' $$

In all of the three cases considered above, we can write any non-zero elements $\alpha_z$ of a $\alpha$ yielding a binary combination as components $x_z$ of $x=M^{-1}b$, where $M$ is an invertible $n \times n$ binary matrix (i.e. having entries of $0$ or $1$), and $b$ an n-component binary vector. Equivalently, $x$ represents the unique solution to $Mx=b$. Cramer's rule for such a solution establishes that the $x_z$ can be written as

$$x_z = \frac{det(M_z)}{det(M)}$$

where $M_z$ is a matrix that replaces the column $z$ of the matrix $M$ with the vector $b$. Since both $M$ and $b$ are composed of binary entries, the matrix $M_z$ is always binary as well. The result now follows as stated in Proposition <ref> since $0$ is always a possible value of $det(M_z)$.

§.§ Proof of Theorem <ref>

§.§.§ Setup and notation

Let $\mathcal{Y} \subseteq \mathbbm{R}$ be the support of $Y$. For any $y \in \mathcal{Y}, z \in \mathcal{Z}$ and $t \in \mathcal{T}$, let $F_{(YD)|Z}(y,t|z) := E[\mathbbm{1}(Y_i \le y)\mathbbm{1}(T_i=t)|Z_i=z]$. The strategy will be based on the fact that knowing the distribution of $(Y_i,T_i,Z_i)$ is equivalent to knowing $F_{(YD)|Z}(y,t|z)$ for all $(y,t,z)$ along with the distribution of the instruments $\mathcal{P}_Z$. For any response type $g \in \mathcal{G}$, we use the notation $T_g(z)$ to denote the common counterfactual treatment functions among units $i$ of that response-type: $G_i=g$.

By the law of iterated expectations and Assumption 1:
\begin{align}F_{(YD)|Z}(y,t|z) &= E[\mathbbm{1}(Y_i \le y)\mathbbm{1}(T_i=t)|Z_i=z] = \sum_{g: T_g(z)=t} P(G_i=g)E[\mathbbm{1}(Y_i(t) \le y)|G_i=g] \nonumber\\
&= \sum_{g: T_g(z)=t} P(G_i=g) F_{Y(t)|G}(y|g): = \sum_{g \in \mathcal{G}} A^{[t]}_{zg} \cdot P(G_i=g) F_{Y(t)|G}(y|g) \label{Aeq}  \end{align}

where $A^{[t]}$ is the $|\mathcal{Z}| \times \mathcal{|\mathcal{G}|}$ matrix with typical entry $[A^{[t]}]_{zg} = \mathbbm{1}(T_g(z)=t)$.

Before proceeding, I state a Lemma that will be useful in what follows:

If $\mu_g^t$ is outcome-agnostic identified given instrument support $\mathcal{Z}$, it remains outcome-agnostic identified using data from $Z_i \in \mathcal{Z}_0$, where $\mathcal{Z}_0 \subseteq \mathcal{Z}$ is a subset of instrument values for which the rows of $A^{[t]}$ for $z \in \mathcal{Z}_0$ are linearly independent of one another.

Suppose that $A^{[t]}$ does not have full row rank. This implies that for some $\mathcal{Z}_0 \subset \mathcal{Z}$, the remaining rows of $A^{[t]}$ for $z \notin \mathcal{Z}_0$ can be written as a linear combination of the rows of $A^{[t]}$ for $z \in \mathcal{Z}_0$. Take one such $z^* \notin \mathcal{Z}_0$, and accordingly let

$$A^{[t]}_{z^*,g} = \sum_{z \in \mathcal{Z}_0} A^{[t]}_{z^*,g} \quad \textrm{ for all } g \in \mathcal{G}$$

Note then that Eq. (<ref>) implies that

\begin{align*}
F_{(YD)|Z}(y,t|z^*) &= P(Y_i \le y  \textrm{ and } T=t|Z_i=z^*)\\
&= P(Y_i(t) \le y \textrm{ and } A^{[t]}_{z^*,G_i}=1|Z_i=z^*)\\
&= P(Y_i(t) \le y \textrm{ and } A^{[t]}_{z^*,G_i}=1)\\
&= P(Y_i \le y \textrm{ and } \sum_{z \in \mathcal{Z}_0} A^{[t]}_{z^*,G_i}=1)
\end{align*}

where the RHS on the last line does not depend in any way on the distribution of observables for $i$ such that $Z_i=z^*$. Thus, $F_{(YD)|Z}(y,t|z^*)$ adds no information that is not contained in $F_{(YD)|Z}(y,t|z)$ for $z \in \mathcal{Z}_0$. If $\mu_g^t$ is outcome-agnostic identified, it must be using the distribution $\mathcal{P}_{YTZ|Z \in \mathcal{Z}_0}$ rather than the full unconditional distribution $\mathcal{P}_{YTZ}$.

Lemma <ref> implies that we can without loss of generality take $A^{[t}$ to have full row-rank. Lemma <ref> extends an observation of <cit.> that one can remove any rows of $A{[t]}$ that is an exact copy of another row (i.e. all selection groups behave the same for two instrument values), and there is hence a direct redundancy over instrument values.

§.§.§ Outcome-agnostic identification

Now define $\mathbf{F}_{(YD)|Z}(y)$ to be a $|\mathcal{T}|\cdot|\mathcal{Z}| \times 1$ vector of $F_{(YD)|Z}(y,t|z)$ over $z$ and $t$ and $\mathbf{G}^*(y)$ to be a $|\mathcal{T}|\cdot |\mathcal{G}|$-component vector of $P(G_i=g)\cdot F_{Y(t)|G}(y|g)$ over $g$ and $t$, for a fixed $y$. Now let $\mathbf{G}^*$ represent the whole vector-valued function $\mathbf{G}^*: \mathcal{Y} \rightarrow \mathbbm{R}^{|\mathcal{T}|\cdot|\mathcal{G}|}$, and define $\mathbf{F}_{(YD)|Z}$ similarly as the function $\mathcal{Y} \rightarrow \mathbbm{R}^{|\mathcal{T}|\cdot|\mathcal{Z}|}$ yielding the vector $F_{(YD)|Z}(y) $. In this notation $\mathcal{P}_Z$ and $\mathbf{F}_{(YD)|Z}$ encode the entire distribution of observables $(Y,T,Z)$ while $\mathcal{P}_Z$ and $\mathbf{G}^*_Y$ encode the entire distribution of model primitives $(Y(t),G,Z)$.

The relationship between the two can be characterized by writing Eq. (<ref>) as:

\begin{equation} \label{Aeq2}
\mathbf{F}_{(YD)|Z} = \mathcal{A} \circ \mathbf{G}^*
\end{equation}

where $\mathcal{A}$ is the linear map of functions  $\mathcal{Y} \rightarrow \mathbbm{R}^{|\mathcal{T}|\cdot |\mathcal{G}|}$ to functions $\mathcal{Y} \rightarrow\mathbbm{R}^{|\mathcal{T}|\cdot|\mathcal{Z}|}$ defined by:

$$ \left[\mathcal{A} \circ \bm{\mu}(y)\right]_{tz} = \sum_{g} A^{[t]}_{z,g} \cdot \bm{\mu}(y)_{t g} $$

holding for each $y$, for any vector-valued function $\bm{\mu}: \mathcal{Y} \rightarrow\mathbbm{R}^{|\mathcal{T}|\cdot|\mathcal{Z}|}$.

Let $\theta_c=\mathbbm{E}[Y_i(t)|c(G_i)=1]$ be the parameter of interest. Note that $\theta_c$ can also be written as a linear map applied to the function $\mathbf{G}^*$. In particular $\theta_c = \Theta_c \circ \mathbf{G}^*$, where for any function $\bm{\mu}: \mathcal{Y}$ to $\mathbbm{R}^{|\mathcal{T}|\cdot|\mathcal{G}|}$, $\Theta_c \circ \bm{\mu}$ is the scalar:

\begin{equation} \label{thetaeq} \sum_{g \in \mathcal{G}} \frac{c_g}{P(c(G_i)=1)} \cdot \int_\mathcal{Y} y\cdot  d\bm{\mu}(y)_{t,g}
\end{equation}

$$\mathcal{S}:= \{\bm{\mu}: \mathcal{A} \circ \bm{\mu} = \mathbf{F}_{(YD)|Z}\}$$
$$\mathcal{R}:= \{\bm{\mu}: [\bm{\mu}(y)]_{tg}/P(G_i=g) \textrm{ is a proper CDF  for each }t \in \mathcal{T} \textrm{ and }g \textrm{ for which } P(G_i=g)>0\}$$

Note that the sets $\mathcal{R}$ and $\mathcal{S}$ as well as the map $\Theta_c$ depend on the distribution $\mathcal{P}_{latent}$ (through $\mathbf{F}_{(YD)|Z}$, through the $P(G_i=g)$, and through $P(c(G_i)=1)$ respectively). Let us denote this dependence by $\mathcal{S}(\mathcal{P}_{obs})$, $\mathcal{R}(\mathcal{P}_{G})$ and $\Theta_c(\mathcal{P}_{G})$, though I will later leave this dependence implicit to ease notation.

Definition <ref> of outcome-agnostic identification from Section <ref>, translated into this notation, says that

$$\left\{\Theta_c(\mathcal{P}_G) \circ \bm{\mu}: \bm{\mu} \in \mathcal{R}(\mathcal{P}_G) \textrm{ and } \bm{\mu} \in \mathcal{S}(\mathscr{P}_{obs}) \textrm{ and } \mathcal{P}_G \in \mathscr{P}_G \right\} \textrm{ is a singleton for all } \mathcal{P}_{obs} \in \mathscr{P}_{obs}(\mathcal{G})$$

\begin{align*}
\mathscr{P}_{obs}(\mathcal{G}) &= \{\phi(\mathcal{P}_{latent} \times \mathcal{P}_Z): \mathcal{P}_{latent} \in \mathscr{P}(\mathcal{G}), \mathcal{P}_Z \in \mathscr{P}_Z\}\\
&= \{\Theta_c(\mathcal{P}_G) \circ \bm{\mu}: \bm{\mu} \in \mathcal{R}(\mathcal{P}_G), \mathcal{P}_{G} \in \mathscr{P}_G \textrm{ such that } supp\{\mathcal{P}_G\} \subseteq \mathcal{G}, \mathcal{P}_Z \in \mathscr{P}_Z\}
\end{align*}

Putting these together, we have that $\theta_c$ is outcome-agnostic identified if and only if

$$ \left\{\Theta_c \cdot \bm{\mu}\right\}_{\bm{\mu} \in (\mathcal{S} \cap \mathcal{R})} \textrm{ is a singleton for all } \mathcal{P}_{latent} \textrm{ with } supp\{\mathcal{P}_{G}(\mathcal{P}_{latent})\} \subseteq \mathcal{G}\vspace{.3cm}$$

§.§.§ A candidate for $\mathbf{G}^*$ that recovers observables

Consider the vector-valued function $\mathbf{G}$, where the $t,g$ component of $\mathbf{G}(y)$ is:

$$\left[\mathbf{G}(y)\right]_{t,g}:= \begin{cases}
P(G_i=g)\cdot F_{Y(t)|G}(y|g) & \textrm{ if } \max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=t)=0\\
\sum_z [(A^{[t]})^+]_{g,z} ]\cdot F_{(YD)|Z}(y,t|z) & \textrm{ if }  \max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=t)=1
\end{cases} $$

and $(A^{[t]})^+$ indicates the Moore-Penrose pseudoinverse of the matrix $A^{[t]}$. The reason for separating out the two cases is that if there exists a group $g$ that acts as a “never-taker” with respect to treatment $t$ such that $\max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=t)=0$, then this corresponds to a column of all zeros in $A^{[t]}$. A property of the Moore-Penrose inverse is that if column $g$ of $A^{[t]}$ is all zeros, then the corresponding row $g$ of $(A^{[t]})^+$ is also all zeros (see e.g. <cit.> which would leave $\left[\mathbf{G}(y)\right]_{t,g}=0$ for all $y$ if we simply defined $\left[\mathbf{G}(y)\right]_{t,g}=\sum_z [(A^{[t]})^+]_{g,z} F_{(YD)|Z}(y,t|z)$ for all $t,g$. This would make it impossible for $\mathbf{G}$ to represent a possible candidate for $\mathbf{G}^*$. The above construction avoids this problem by simply replacing such problematic combinations of $(g,t)$ by using the actual $[\mathbf{G}^*(y)]_{t,g}$. Note that if the first case holds for all $g \in \mathcal{G}$, then the matrix $A^{[t]}$ is simply the zero matrix, and outcome agnostic identification cannot hold, by Lemma <ref>. Thus, we can continue under the assumption that the second case holds for at least some $g \in \mathcal{G}$.

We first observe that $\mathbf{G}$ “recovers observables”, by which I mean that $\mathcal{A} \circ \bm{\mu} = \mathbf{F}_{(YD)|Z}$ and hence $\mathbf{G} \in \mathcal{S}$. To see this, note that:

\begin{align*}[\mathcal{A} \circ \mathbf{G}(y)]_{t,z} &= \sum_{g} A^{[t]}_{z,g}\left[\mathbf{G}(y)\right]_{t,g}\\
&=\cancel{\sum_{g: \max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=t)=0} A^{[t]}_{z,g}\cdot  P(G_i=g)\cdot F_{Y(t)|G}(y|g)}\\
&\hspace{1in} + \sum_{g: \max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=t)=1}\sum_{z'} A^{[t]}_{z,g}[(A^{[t]})^+]_{g,z'} F_{(YD)|Z}(y,t|z')\\
&=\sum_{g,z'} A^{[t]}_{z,g}[(A^{[t]})^+]_{g,z'} F_{(YD)|Z}(y,t|z')\\
&= \sum_{z'} [A^{[t]}(A^{[t]})^+]_{z,z'} F_{(YD)|Z}(y,t|z')= [F_{(YD)|Z}(y)]_{tz} \end{align*}

where the second and third equalities use that $A^{[t]}_{z,g} = 0$ for all $z$, if $g$ is such that $\max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=t)=0$. The final equality follows from $A^{[t]}(A^{[t]})^+ = I_{|\mathcal{Z}|}$, which in turn follows from $(A^{[t]})^+ = {A^{[t]}}'({A^{[t]}}{A^{[t]}}')^{-1}$ since we can by Lemma <ref> assume that $A^{[t]}$ has full row rank.

$\mathbf{G}$ may however not be in $\mathcal{R}$, as its definition above does not ensure that each $F_{Y(t)|G}(y|g)$ is necessarily weakly increasing in $y$ with a limit of unity as $y \uparrow \infty$. Note that $[\mathbf{G}]_{t,g}/P(G_i=g)$ does have the final two properties of a CDF: right-continuity and a left limit of zero. To see this, substitute Eq. (<ref>) into the definition of $\mathbf{G}$, to rewrite as:

\begin{equation} \label{eq:Fstarnew}
\left[\mathbf{G}(y)\right]_{t,g}:= \begin{cases}
P(G_i=g)\cdot F_{Y(t)|G}(y|g) & \textrm{ if } \max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=t)=0\\
\sum_{g'} [(A^{[t]})^+ A^{[t]}]_{g,g'}\cdot P(G_i=g')\cdot F_{Y(t)|G}(y|g') & \textrm{ if } \max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=t)=1
\end{cases}
\end{equation}

Right continuity of each element of $\mathbf{G}(y)$ in $y$ follows from right-continuity of the $F_{Y(t)|G}(y|g')$. Note that $\lim_{y \downarrow -\infty} \left[\mathbf{G}(y)\right]_{t,g} = 0$ follows from each of the CDFs $F_{(YD)|Z}$ approaching zero as $y \downarrow -\infty$, given that the components of $A^{[t]}$ and $P(G_i=g)$ are finite.

Let $\beta_{t,g} := \lim_{y \uparrow \infty} \left[\mathbf{G}(y)\right]_{t,g}$. For any $t,g$ such that $\max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=t)=0$, it follows from the definition of
$\mathbf{G}$ that $\beta_{t,g}=P(G_i=g)$, since each of the $F_{Y(t)|G}(y|g)$ are valid CDFs. For the other $t,g$, use (<ref>) to see that

\begin{align*}
\beta_{t,g}&=\lim_{y \uparrow \infty} \sum_{g'} [(A^{[t]})^+ A^{[t]}]_{g,g'} \cdot P(G_i=g')\cdot F_{Y(t)|G}(y|g') = \sum_{g'} [(A^{[t]})^+ A^{[t]}]_{g,g'} \cdot P(G_i=g')\\
&= [(A^{[t]})^+ A^{[t]}P]_{g}
\end{align*}

where $P$ is a vector of $P(G_i=g)$ for all $g \in \mathcal{G}$.

The following definition will be useful in the construction below. Pick any fixed $g^* \in \mathcal{G}$, and let us define a vector valued function $\mathbf{D}: \mathcal{Y}\rightarrow \mathbbm{R}^{|\mathcal{T}|\cdot |\mathcal{G}|}$ with components:

\begin{equation} \label{eq:D}
[\mathbf{D}(y)]_{t,g} := (P_g-\beta_{t,g})\cdot F_{Y(t)|G}(y|g^*)= [\left\{I-(A^{[t]})^+ A^{[t]}\right\}P]_{g}\cdot F_{Y(t)|G}(y|g^*)
\end{equation}

§.§.§ A class of alternative candidates that also recover observables

Now let us define for any $\lambda \in [0,1]$ the convex combination of $\mathbf{G}+\mathbf{D}$ and $\mathbf{G}^*$:

\begin{equation} \label{eq:Galpha}
\mathbf{G}^{\lambda} := \lambda\left(\mathbf{G} + \mathbf{D} \right) + (1-\lambda)\mathbf{G}^* = \mathbf{G}^*+\lambda\left\{\mathbf{G}-\mathbf{G}^*+\mathbf{D}\right\}
\end{equation}

Our first observation will be that $\mathcal{A} \circ \mathbf{G}^{\lambda} = \mathbf{F}_{(YD)|Z}$, and thus $\mathbf{G}^{\lambda} \in \mathcal{S}$. To see this, note that:

\begin{align*}[\mathcal{A} \circ \mathbf{G}^\lambda(y)]_{t,z} &= [\mathcal{A} \circ \mathbf{G}(y)]_{t,z} + \lambda \cdot [\mathcal{A} \circ \left\{\mathbf{G}-\mathbf{G}^*+\mathbf{D}\right\}(y)]_{t,g}\\
&=[F_{(YD)|Z}(y)]_{t,z} + \cancel{\lambda \cdot [\mathcal{A} \circ\mathbf{G}(y)]_{t,g}-\lambda \cdot [\mathcal{A} \circ\mathbf{G}^*(y)]_{t,g}} + \lambda \cdot [\mathcal{A} \circ\mathbf{D}(y)]_{t,g}\\
&=[F_{(YD)|Z}(y)]_{t,z}+\lambda \cdot \sum_{g,g'} A^{[t]}_{z,g} \cdot [(I-(A^{[t]})^+ A^{[t]})]_{g,g'}\cdot P(G_i=g')\cdot  F_{Y(t)|G}(y|g^*)\\
&=[F_{(YD)|Z}(y)]_{t,z}+\lambda \cdot \sum_{g'} [\cancel{A^{[t]}(I-(A^{[t]})^+ A^{[t]})}]_{z,g'}\cdot P(G_i=g')\cdot  F_{Y(t)|G}(y|g^*)\\
\end{align*}

since $\mathcal{A} \circ\mathbf{G}^*=\mathcal{A} \circ\mathbf{G}$ and $A^{[t]}(A^{[t]})^+ A^{[t]}=A^{[t]}$.

Now, we verify that for a small enough $\lambda$, $\mathbf{G}^{\lambda}$ yields $F_{Y(t)|G}(y|g)$ that satisfy the properties of a CDF and hence $\mathbf{G}^{\lambda} \in \mathcal{R}$.

First, note that $\left[\mathbf{G}^{\lambda}(y)\right]_{t,g}$ is right-continuous in $y$, since each of $[\mathbf{G}(y)]_{t,g}$, $[\mathbf{G}^*(y)]_{t,g}$, and $[\mathbf{D}(y)]_{t,g}$ are. We also have that $\lim_{y \downarrow -\infty} \left[\mathbf{G}^{\lambda}(y)\right]_{t,g}=0$, since
$$\lim_{y \downarrow -\infty} \left[\mathbf{G}(y)\right]_{t,g}=\lim_{y \downarrow -\infty} \left[\mathbf{G}^{*}(y)\right]_{t,g}=\lim_{y \downarrow -\infty} \left[\mathbf{D}(y)\right]_{t,g}=0$$

Note as well that

\begin{align*}
\lim_{y \uparrow \infty} \left[\mathbf{G}^{\lambda}(y)\right]_{t,g} &= \lim_{y \uparrow \infty} \left[\mathbf{G}^{*}(y)\right]_{t,g}+\lambda \cdot \lim_{y \uparrow \infty} \left[\left\{\mathbf{G}-\mathbf{G}^*+\mathbf{D}\right\}(y)\right]_{t,g} \\
&=P_g+\lambda \cdot \left\{\lim_{y \uparrow \infty} \left[\mathbf{G}(y)\right]_{t,g}-\lim_{y \uparrow \infty} \left[\mathbf{G}^{*}(y)\right]_{t,g}+\lim_{y \uparrow \infty} \left[\mathbf{D}(y)\right]_{t,g}\right\}\\
&=P_g+\lambda \cdot \left\{\beta_{t,g}-P_g+(P_g-\beta_{t,g})\cdot 1\right\} = P_g
\end{align*}

matching the correct normalization of the true $\lim_{y \uparrow \infty} \left[\mathbf{G}^{*}(y)\right]_{t,g} = P_g \cdot \lim_{y \uparrow \infty} F_{Y(t)|G=g}(y) = P_g$.

It only remains to be seen that for a small enough value of $\lambda$, $ \left[\mathbf{G}^{\lambda}(y)\right]_{t,g}$ is weakly increasing in $y$. This is where Assumption REG will prove useful.

Given Assumption REG, $\left[\mathbf{G}^{\lambda}(y)\right]_{t,g}$ is non-decreasing in $y$ for any $\lambda \in (0,\bar{\lambda}]$, where $\bar{\lambda} = \frac{\underbar{L}}{2|\mathcal{G}|\cdot \bar{L}} >0$.

Given Proposition <ref>, we have shown that for $\lambda \le \bar{\lambda}$, $\mathbf{G}^\lambda \in \mathcal{R}$ and hence $\mathbf{G}^\lambda \in (\mathcal{S} \cap \mathcal{R})$.

§.§.§ Outcome agnostic identification implies $c \in rs(A^{[t]})$

Thus, identification of $\theta_c$ requires that for any such $\lambda>0$ that is small enough that $\mathbf{G}^\lambda \in (\mathcal{S} \cap \mathcal{R}$), it holds that $\Theta_c \circ \mathbf{G}^\lambda = \Theta_c \circ \mathbf{G}^*$. This in turn requires, by Eq. (<ref>), that $\Theta_c \circ  \left\{\mathbf{G}-\mathbf{G}^*+\mathbf{D}\right\}=0$. Now:

\begin{align}
\Theta_c &\circ  \left\{\mathbf{G}-\mathbf{G}^*+\mathbf{D}\right\} =\frac{1}{P(c(G_i))} \sum_g c_g
\cdot \left\{\int_\mathcal{Y} y\cdot  d\mathbf{G}(y)_{t,g}-\int_\mathcal{Y} y\cdot  d\mathbf{G}^*(y)_{t,g}+\int_\mathcal{Y} y\cdot  d\mathbf{D}(y)_{t,g}\right\} \nonumber\\
&=\frac{1}{P(c(G_i))} \sum_g c_g \sum_{g'} [I-(A^{[t]})^+ A^{[t]}]_{g,g'} \cdot P(G_i=g')\cdot \mathbbm{E}[Y_i(t)|G_i=g'] \nonumber\\
&\hspace{.5in}+\frac{1}{P(c(G_i))} \sum_g c_g \sum_{g'} [(I-(A^{[t]})^+ A^{[t]})P]_{g,g'}\cdot P(G_i=g') \cdot \mathbbm{E}[Y_i(t)|G_i=g^*] \nonumber\\
&=\frac{1}{P(c(G_i))} \sum_{g'} [c'(I-(A^{[t]})^+ A^{[t]})]_{g'} \cdot P(G_i=g')\cdot \mathbbm{E}[Y_i(t)|G_i=g'] \nonumber\\
&\hspace{1in}+\frac{1}{P(c(G_i))} \sum_{g'} [c'(I-(A^{[t]})^+ A^{[t]})]_{g'}\cdot P(G_i=g')\cdot \mathbbm{E}[Y_i(t)|G_i=g^*] \nonumber\\
&=\frac{1}{P(c(G_i))} \sum_{g'} [c'(I-(A^{[t]})^+ A^{[t]})]_{g'} \cdot P(G_i=g')\cdot \left\{ \mathbbm{E}[Y_i(t)|G_i=g']-\mathbbm{E}[Y_i(t)|G_i=g^*]\right\} \label{eq:finalcondtion}
\end{align}

For this to be true for all joint distributions $\mathcal{P}_{latent}$ compatible with $\mathcal{G}$, we must have $[c'(I-(A^{[t]})^+ A^{[t]})]_{g'}=0$ for all $g' \in \mathcal{G}$, i.e. $c$ is in the rowspace of $A^{[t]}$. Otherwise, if $c'(I-(A^{[t]})^+ A^{[t]}) = \tilde{c}'$ for some non-zero vector $\tilde{c}$, we could always construct a $\mathcal{P}_{latent}$ yielding a non-zero vector of $P(G_i=g')\cdot \left\{ \mathbbm{E}[Y_i(t)|G_i=g']-\mathbbm{E}[Y_i(t)|G_i=g^*]\right\}$ across $g'$ that is parallel to $\tilde{c}$ in $\mathbbm{R}^{|\mathcal{G}|}$, such that $\sum_{g'} \tilde{c}'_{g'} \cdot P(G_i=g')\cdot \left\{ \mathbbm{E}[Y_i(t)|G_i=g']-\mathbbm{E}[Y_i(t)|G_i=g^*]\right\} > 0$.

Deriving Eq. (<ref>): Finally, let us see explicitly that if $c$ lies in the row-space of $A^{[t]}$, we can write $\theta_c$ in the form of Eq. (<ref>). Given the definition of $\mathbf{G}$, this allows us to express $\theta_c=\mathbbm{E}[Y_i(t)|c(G_i)=1]$ as:

\begin{align}
\theta_c &=\sum_{g} \frac{c_g}{\mathbbm{E}[c_{G_i}]} \int y d\left\{\sum_{z} [(A^{[t]})^+]_{g,z} F_{(YD_T)|Z}(y,t|z))\right\} \nonumber\\
&=\sum_{g} \frac{c_g}{\mathbbm{E}[c_{G_i}]}\sum_{z} [(A^{[t]})^+]_{g,z} P(T_i=t|Z_i=z)\int y dF_{Y|Z,T}(y|z,t) \nonumber\\
&=\sum_{z} \left(\sum_g [(A^{[t]})^+]_{g,z} \frac{c_g}{\mathbbm{E}[c_{G_i}]}\right) \cdot \mathbbm{E}[Y_iD^{[t]}_i|Z_i=z] \nonumber\\
&=\frac{\sum_{z} \left(\sum_g [(A^{[t]})^+]_{g,z} \cdot c_g \right) \cdot \mathbbm{E}[Y_iD^{[t]}_i|Z_i=z]}{\mathbbm{E}[c_{G_i}]} \label{eq:necc}
\end{align}

The numerator of (<ref>) takes the form of a linear combination in which $K = |\mathcal{Z}|$ with coefficients $\alpha_k = \sum_g [(A^{[t]})^+]_{g,z_k} \cdot c_g$, for an arbitrary ordering of the points $z_1 \dots z_K$ in $\mathcal{Z}$.

It only remains to be shown that $(t,\alpha)$ is a binary combination, and that furthermore $\mathbbm{E}[c_{G_i}] = \sum_{k} \alpha_{k}  \cdot \mathbbm{E}[D^{[t]}_i(z_k)] = \sum_{k} \alpha_{k}  \cdot \mathbbm{E}[D^{[t]}_i|Z_i=z_k]$. This by Proposition (<ref>) then establishes the final result that $\mathbbm{E}[Y_i(t)|c_{G_i}=1] = \mathbbm{E}[Y_i(t)|C_i^{[\alpha,t]}]$. For $(t,\alpha)$ to be a binary combination we must first verify that $\sum_{k} \alpha_{k} \cdot D^{[t]}_i(z_k) \in \{0,1\}$ for all $i$. Note that for any $z$, $D^{[t]}_i(z) = [A^{[t]}]_{z,G_i}$ by the definition of $A^{[t]}$. Thus:

\begin{align*}
\sum_{k} \alpha_{k}  \cdot D^{[t]}_i(z_k) &= \sum_{z} \left(\sum_g [(A^{[t]})^+]_{g,z} \cdot c_g \right) \cdot [A^{[t]}]_{z,G_i}\\
&= \sum_g c_g \cdot [(A^{[t]})^+(A^{[t]})]_{g,G_i} = [[(A^{[t]})^+(A^{[t]})]' c]_{G_i} = [[(A')(A')^+]c]_{G_i}
\end{align*}

where we let $A'$ denote the transpose of $A^{[t]}$ (dropping the $t$ index for clarity). Since $c$ lies in the row space of $A^{[t]}$ (equivalently $c$ lies in the column space of $A'$), we have that $[(A')(A')^+]c = c$. Thus have by the above that $\sum_{k} \alpha_{k}  \cdot D^{[t]}_i(z_k) = c_{G_i}=\mathbbm{1}[c(G_i)=1]$, establishing that $\sum_{k} \alpha_{k}  \cdot D^{[t]}_i(z_k) \in \{0,1\}$ for all $i$, and also that $P[c(G_i)=1] = \sum_{k} \alpha_{k}  \cdot \mathbbm{E}[D^{[t]}_i(z_k)] = \sum_{k} \alpha_{k}  \cdot \mathbbm{E}[D^{[t]}_i|Z_i=z_k]$.

§.§.§ Proof of Proposition <ref>

The key to monotonicity will be to choose $\lambda$ small enough that any decreases in the components of $\mathbf{G}^{\lambda}$ with $y$ is dominated by increases in the corresponding components of $\mathbf{G}^*$, so that each $\left[\mathbf{G}^{\lambda}\right]_{t,g}$ is monotonically increasing. Note that $\left[\mathbf{G}^\lambda(y)\right]_{t,g}$ is monotonically increasing in $y$ if for any $y'>y$:

$$\lim_{y' \downarrow y} \left[\mathbf{G}^\lambda(y')\right]_{t,g} - \left[\mathbf{G}^\lambda(y)\right]_{t,g}\ge 0$$

i.e. that

\begin{equation} \label{eq:ineq}
\lim_{y' \downarrow y} \left[\mathbf{G}^*(y')\right]_{t,g} - \left[\mathbf{G}^*\right]_{t,g}\ge \lambda\cdot \{\lim_{y' \downarrow y} \left[(\mathbf{G}^*-\mathbf{G})(y')-(\mathbf{G}^*-\mathbf{G})(y)\right]_{t,g} - (\lim_{y' \downarrow y} \left[\mathbf{D}(y')\right]_{t,g} - \left[\mathbf{D}\right]_{t,g})\}
\end{equation}

Let us turn first to $\left[(\mathbf{G}^*-\mathbf{G})(y)\right]_{t,g}$. Fix a $g$ and $t$, and any $y' > y$. Then, by (<ref>):

\begin{equation} \label{eq:Fstarnew2}
\left[\mathbf{G}(y')\right]_{t,g}-\left[\mathbf{G}^*(y)\right]_{t,g}= \begin{cases}
\sum_{g'} [(A^{[t]})^+ A^{[t]}]_{g,g'}\cdot  P(G_i=g')\cdot \left\{F_{Y(t)|G}(y'|g')-F_{Y(t)|G}(y|g')\right\}
\end{cases}
\end{equation}

where the first line indicates the case that $g$ is such that $\max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=t)=0$, and the second that $\max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=t)=1$. Thus $\left[(\mathbf{G}^*-\mathbf{G})(y')\right]_{t,g}-\left[(\mathbf{G}^*-\mathbf{G})(y)\right]_{t,g}$ is equal to $0$ if $\max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=t)=0$, and

$$\sum_{g'} [I-(A^{[t]})^+ A^{[t]}]_{g,g'} \cdot P(G_i=g')\cdot \left\{F_{Y(t)|G}(y'|g')-F_{Y(t)|G}(y|g')\right\}$$

if $\max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=t)=1$.

Thus we have by REG that if $\left[(\mathbf{G}^*-\mathbf{G})(y')\right]_{t,g}\ne \left[(\mathbf{G}^*- \mathbf{G})(y)\right]_{t,g}$:

\begin{align*}
&= \left|\sum_{g'} [I-(A^{[t]})^+ A^{[t]}]_{g,g'} \cdot P(G_i=g')\cdot \left\{F_{Y(t)|G}(y'|g')- F_{Y(t)|G}(y|g')\right\}\right|\\
&= \left\{F_{Y(t)|G}(y'|g^*)- F_{Y(t)|G}(y|g^*)\right\}\cdot \left|\sum_{g'} [I-(A^{[t]})^+ A^{[t]}]_{g,g'} \cdot P(G_i=g') \cdot \frac{F_{Y(t)|G}(y'|g')- F_{Y(t)|G}(y|g')}{F_{Y(t)|G}(y'|g^*)- F_{Y(t)|G}(y|g^*)}\right|\\
& \le \left\{F_{Y(t)|G}(y'|g^*)- F_{Y(t)|G}(y|g^*)\right\}\cdot  |\mathcal{G}|^{1/2} \cdot \sqrt{\sum_{g'} P(G_i=g')^2\cdot  \left(\frac{F_{Y(t)|G}(y'|g')- F_{Y(t)|G}(y|g')}{F_{Y(t)|G}(y'|g^*)- F_{Y(t)|G}(y|g^*)}\right)^2}\\
& \le \left\{F_{Y(t)|G}(y'|g^*)- F_{Y(t)|G}(y|g^*)\right\}\cdot |\mathcal{G}| \cdot \max_{g'} P(G_i=g') \cdot \max_{g'} \left|\frac{F_{Y(t)|G}(y'|g')- F_{Y(t)|G}(y|g')}{F_{Y(t)|G}(y'|g^*)- F_{Y(t)|G}(y|g^*)}\right|\\
& \le \left\{F_{Y(t)|G}(y'|g^*)- F_{Y(t)|G}(y|g^*)\right\}\cdot |\mathcal{G}| \cdot \max_{g'} \left|\frac{F_{Y(t)|G}(y'|g')- F_{Y(t)|G}(y|g')}{F_{Y(t)|G}(y'|g^*)- F_{Y(t)|G}(y|g^*)}\right|
\end{align*}

using that $[I-(A^{[t]})^+ A^{[t]}]$ is a projection (so that $|[I-(A^{[t]})^+ A^{[t]}]v| \le |v|$ for any vector $v \in \mathbbm{R}^{|\mathcal{G}|}$) and by the Cauchy-Schwarz inequality. Let $\delta_t^*(y):=\lim_{y' \downarrow y} \left\{F_{Y(t)|G}(y'|g^*)- F_{Y(t)|G}(y|g^*)\right\}$, which is guaranteed to exist by right continuity of $F_{Y(t)|G}(y|g^*)$. Then, by REG:

$$\lim_{y' \downarrow y} \left|\left[(\mathbf{G}^*-\mathbf{G})(y')\right]_{t,g}-\left[(\mathbf{G}^*-\mathbf{G})(y)\right]_{t,g}\right| \le \delta_t^*(y) \cdot |\mathcal{G}| \cdot \bar{L}$$

Now consider $\left[\mathbf{D}(y)\right]_{t,g}$. Fix a $g$ and $t$, and any $y' > y$. Similarly, we have that

\begin{align*}
\left|\left[\mathbf{D}(y')\right]_{t,g}-\left[\mathbf{D}(y)\right]_{t,g}\right|&= \left|\sum_{g'} [I-(A^{[t]})^+ A^{[t]}]_{g,g'} \cdot P(G_i=g')\cdot \{F_{Y(t)|G}(y'|g^*)-F_{Y(t)|G}(y|g^*)\}\right|\\
&\le \{F_{Y(t)|G}(y'|g^*)-F_{Y(t)|G}(y|g^*)\} \cdot |P|\\
&\le \{F_{Y(t)|G}(y'|g^*)-F_{Y(t)|G}(y|g^*)\} \cdot |\mathcal{G}|
\end{align*}

So, using Assumption REG:

$$ \lim_{y' \downarrow y} \left|\left[\mathbf{D}(y')\right]_{t,g}-\left[\mathbf{D}(y)\right]_{t,g}\right| \le \delta_t^*(y) \cdot |\mathcal{G}|\cdot \bar{L}$$

We can thus put an upper bound on the RHS of (<ref>)

$$\lambda\cdot \{\lim_{y' \downarrow y} \left[(\mathbf{G}^*-\mathbf{G})(y')-(\mathbf{G}^*-\mathbf{G})(y)\right]_{t,g} - (\lim_{y' \downarrow y} \left[\mathbf{D}(y')\right]_{t,g} - \left[\mathbf{D}\right]_{t,g})\} \le 2\lambda \cdot \delta_t^*(y) \cdot |\mathcal{G}|\cdot \bar{L}$$

Meanwhile, by REG:

\begin{align*}
\lim_{y' \downarrow y} \left\{\left[\mathbf{G}^*(y')\right]_{t,g} - \left[\mathbf{G}^*\right]_{t,g}\right\} & = \lim_{y' \downarrow y} \{F_{Y(t)|G}(y'|g^*)-F_{Y(t)|G}(y|g^*)\} \cdot \lim_{y' \downarrow y} \frac{F_{Y(t)|G}(y'|g')- F_{Y(t)|G}(y|g')}{F_{Y(t)|G}(y'|g^*)- F_{Y(t)|G}(y|g^*)}\\
\ge \delta_t^*(y) \cdot \underbar{L}
\end{align*}

Thus inequality (<ref>) holds provided that

$$\delta_t^*(y) \cdot \underbar{L} \ge  2\lambda \cdot \delta_t^*(y) \cdot |\mathcal{G}|\cdot \bar{L} \quad \iff \quad \lambda \le \frac{\underbar{L}}{2|\mathcal{G}|\cdot \bar{L}}$$

A visualization of the intuition behind this result is depicted in Figure <ref>.

Depiction of the monotonicity result. Blue sinusoidal function depicts an example of a $\left[(\mathbf{G}+\mathbf{D})^\lambda(y)\right]_{t,g}$ that is not weakly increasing. Orange curve depicts $\left[\mathbf{G}^*(y)\right]_{t,g}$ which is weakly increasing. Black curve depicts $\left[\mathbf{G}^\lambda(y)\right]_{t,g}$, which is a linear combination of the blue and orange functions with weights $\lambda=0.1$ and $1-\lambda = 0.9$, respectively. This value of $\lambda$ is small enough that the black curve is weakly increasing everywhere.

§ EXTENDED ANALYSIS OF IDENTIFICATION UNDER NSOG

It is known that unconditional means $\mathbbm{E}[Y_i(t)]$ of a given potential outcome $Y_i(t)$ can be point-identified, given an order condition on the instruments, under an assumption of “no-selection on gains” (NSOG) (see e.g. <cit.> for versions of this result).[<cit.> does not use the terminology of NSOG, but rather “constant average treatment effects”.] Note that identification of $\mathbbm{E}[Y_i(t)]$ and $\mathbbm{E}[Y_i(t')]$ implies identification of unconditional average treatment effects $\mathbbm{E}[Y_i(t')-Y_i(t)]$ as well.

NSOG says that treatment effects are mean independent of actual treatment, given any realization of the instruments:
[NOSG (no selection on gains)]
For any $t,t',t_1,t_2 \in \mathcal{T}$ and $z \in \mathcal{Z}$:

$$\mathbbm{E}[Y_i(t')-Y_i(t)|T_i=t_1,Z_i=z] = \mathbbm{E}[Y_i(t')-Y_i(t)|T_i=t_2,Z_i=z]$$

NSOG implies that if we consider any fixed treatment value $0 \in \mathcal{T}$, then $\mathbbm{E}[Y_i(t')-Y_i(0)|T_i=t,Z_i=z] = \mathbbm{E}[Y_i(t')-Y_i(0)|Z_i=z]$ for any $t,z$, which coupled with independence (<ref>) in turn implies that $\mathbbm{E}[Y_i(t')-Y_i(0)|T_i=t,Z_i=z] = \mathbbm{E}[Y_i(t')-Y_i(0)] := \Delta_{t'}$, where note that $\Delta_{t'}$ does not depend on $z$ or $t$. This normalization against an arbitrary treatment $0 \in \mathcal{T}$ allows us to carry around one less index in our expressions.

§.§ Identification under NSOG

This subsection shows that $\mathbbm{E}[Y_i(t)]$ can be identified for each $t \in \mathcal{T}$ under NSOG given rich enough support of the instruments. The proof essentially follows that of <cit.>, which adapts an argument from <cit.> to cases in which the treatments $\mathcal{T}$ are not necessarily ordered.

NSOG implies that with probability one:

$$\mathbbm{E}[Y_i(t)-Y_i(0)|T_i,Z_i=z] = \Delta_t$$

for any $t$ and thus letting $t=T_i$:

$$\mathbbm{E}[Y_i(T_i)-Y_i(0)|T_i,Z_i=z] = \Delta_{T_i} = \sum_{t \in \mathcal{T}} \mathbbm{1}(T_i=t) \cdot \Delta_t$$

again with probability one. Averaging over the conditional distribution of $T_i$ given $Z_i=z$, we have

$$\mathbbm{E}[Y_i-Y_i(0)|Z_i=z] = \sum_{t \in \mathcal{T}} P(T_i=t|Z_i=z) \cdot \Delta_t$$

To now see that $\mathbbm{E}[Y_i(t)]$ can be identified under NSOG given rich enough instrument support, let us assume that $|\mathcal{Z}| \ge |\mathcal{T}|$ and suppose that there exists a set of $|\mathcal{T}|$ instrument values $\tilde{\mathcal{Z}} \subseteq \mathcal{Z}$ such that the $|\mathcal{T}| \times |\mathcal{T}|$ matrix $\Sigma$ with entries $\Sigma_{zt} = P(Z_i=z, T_i=t)$ over all $z \in \tilde{\mathcal{Z}}$ is invertible, with $P(Z_i = z) > 0$ for each $z \in \tilde{\mathcal{Z}}$. The above equation can be re-written

\} \cdot \mathbbm{1}(Z_i=z)] = \sum_{t \in \mathcal{T}} \Sigma_{zt} \cdot \Delta_t$$

for each $z \in \tilde{\mathcal{Z}}$, or equivalently

\begin{align*}
\mathbbm{E}[Y_i \cdot \mathbbm{1}(Z_i=z)] &= P(Z_i=z)\cdot \mathbbm{E}[Y_i(0)] +  \sum_{t \in \mathcal{T}} \Sigma_{zt} \cdot \Delta_t\\
&= P(Z_i=z)\cdot \mathbbm{E}[Y_i(0)] +  \sum_{t \in \mathcal{T}, t \ne 0} \Sigma_{zt} \cdot \Delta_t\\
&= P(Z_i=z,T_i=0)\cdot \mathbbm{E}[Y_i(0)] + \sum_{t \in \mathcal{T}, t \ne 0} \Sigma_{zt} \cdot \mathbbm{E}[Y_i(0)] + \sum_{t \in \mathcal{T}, t \ne 0} \Sigma_{zt} \cdot \Delta_t\\
&= P(Z_i=z,T_i=0)\cdot \mathbbm{E}[Y_i(0)] + \sum_{t \in \mathcal{T}, t \ne 0} \Sigma_{zt} \cdot \mathbbm{E}[Y_i(t)]\\
& = \sum_{t \in \mathcal{T}} \Sigma_{zt} \cdot \mathbbm{E}[Y_i(t)]
\end{align*}

using that $\Delta_0 = 0$ in the second equality. This yields a system of $|\mathcal{T}|$ equations in the $|\mathcal{T}|$ unknowns $\mathbbm{E}[Y_i(t)]$ with identified coefficients $\Sigma_{zt}$. Given that $\Sigma^{-1}$ is invertible, we have then that

$$\mathbbm{E}[Y_i(t)] = \sum_{z \in \tilde{\mathcal{Z}}} \Sigma^{-1}_{tz} \cdot \mathbbm{E}[Y_i \cdot \mathbbm{1}(Z_i=z)]$$

§.§ How Theorem <ref> does not cover NSOG

Since the result of the last section makes no assumption about which response groups can show up in the population, it is compatible with any selection model $\mathcal{G} \subseteq \{0,1\}^{\mathcal{T}^\mathcal{Z}}$, including for example the full powerset $\{0,1\}^{\mathcal{T}^\mathcal{Z}}$ of possible selection groups $\mathcal{T}^\mathcal{Z}$.

Whatever $\mathcal{G}$ is, unconditional means like $\mathbbm{E}[Y_i(t)]$ correspond to the choice $c = (1, \dots, 1)'$ in $\mathbbm{R}^{|\mathcal{G}|}$. As long as $\mathcal{G}$ allows never-takers with respect to treatment $t$, this choice of $c$ will not lie in the rowspace of $A^{[t]}$ and hence $[c'(I-(A^{[t]})^+ A^{[t]})]_{g'}$ in the final step of the proof of Theorem <ref>, Eq. <ref>,  will be non-zero for at least some $g' \in \mathcal{G}$. The unrestricted selection model $\mathcal{G} = \{0,1\}^{\mathcal{T}^\mathcal{Z}}$ has the property that there exist never-takers with respect to any given $t \in \mathcal{T}$. Recall that Eq. (<ref>) says that if $\mu_c^t$ is outcome-agnostic identified, then:

\begin{equation} \label{eq:finalcondition2}
\sum_{g'} [c'(I-(A^{[t]})^+ A^{[t]})]_{g'} \cdot P(G_i=g')\cdot \left\{ \mathbbm{E}[Y_i(t)|G_i=g']-\mathbbm{E}[Y_i(t)|G_i=g^*]\right\} = 0
\end{equation}

where $g^*$ is an arbitrary fixed selection group in $\mathcal{G}$. Therefore, point identification of $\mathbbm{E}[Y_i(t)]$ under NSOG must be possible, even though the coefficients on $P(G_i=g')\cdot \left\{ \mathbbm{E}[Y_i(t)|G_i=g']-\mathbbm{E}[Y_i(t)|G_i=g^*]\right\}$ in the above are not all equal to zero.

Note that one way to satisfy (<ref>) without $[c'(I-(A^{[t]})^+ A^{[t]})]_{g'}=0$ for all $g'$ is to assume that $\mathbbm{E}[Y_i(t)|G_i=g']=\mathbbm{E}[Y_i(t)|G_i=g^*]$ for all $g' \in \mathcal{G}$. However, this is a much stronger assumption than NSOG, and in fact entirely rules out endogeneity of treatment, in the sense that a simple difference in means $\mathbbm{E}[Y_i|T_i=t'] - \mathbbm{E}[Y_i|T_i=t]$ would then be equal to the average treatment effect $\mathbbm{E}[Y_i(t')-Y_i(t)]$ between $t'$ and $t$.

In the case of NSOG, identification of the $\mathbbm{E}[Y_i(t)]$ remains consistent with the proof of Theorem <ref> in a less direct way. Recall the set $\mathcal{R}$ from the proof of Theorem <ref>, which consists of all vector valued functions $\bm{\mu}: \mathcal{Y} \rightarrow \mathbbm{R}^{|\mathcal{T}| \times |\mathcal{Z}|}$ that yield proper CDF functions $[\bm{\mu}(y)]_{tg}/P(G_i=g)$ for response groups $g$ that occur in the population. However, some of the $\bm{\mu} \in \mathcal{R}$ that may nevertheless violate NSOG, though the produce valid CDFs. Intuitively, the set $\mathcal{R}$ allowed in this proof may be “too big” if the assumption of NSOG is additionally maintained.

Let us see this possibility by considering the construction $\mathbf{G}^\lambda$ used in the proof of Theorem <ref> as a candidate value for $\mathbf{G}^*$ the true potential outcome CDFs, scaled by the selection group probabilities, i.e. $[\mathbf{G}^*(y)]_{t,g}:=P(G_i=g)\cdot F_{Y(t)|G}(y|g)$.

Consider any $g \in \mathcal{G}$ such that $\max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=t)=1$ and $\max_{z \in \mathcal{Z}} \mathbbm{1}(T_g(z)=0)=1$, where $0 \in \mathcal{T}$ is the arbitrary reference treatment. Then if $\mathbf{G}^*=\mathbf{G}$ (corresponding to the limit $\lambda \rightarrow 0$) we have by Eq. (<ref>) that:

$$P(G_i=g)\cdot \mathbbm{E}[Y_i(t)|G_i=g] = \sum_{g'} [(A^{[t]})^+A^{[t]}]_{g,g'}\cdot P(G_i=g')\cdot F_{Y(t)|G}(y|g')$$

and similarly for treatment $0$, so that

\begin{align*}
\mathbbm{E}[Y_i(t)|G_i=g, Z_i=z] = \mathbbm{E}[Y_i(t)|G_i=g]= \sum_{g'} \frac{P(G_i=g')}{P(G_i=g)}\cdot(A^{[t]})^+A^{[t]}]_{g,g'}\cdot \mathbbm{E}[Y(t)|G_i=g']
\end{align*}

using independence (<ref>), and similarly for treatment $0$. Therefore, for any $t_1 \in \mathcal{T}$:

\begin{align}
\mathbbm{E}[&Y_i(t)-Y_i(0)|T_i=t_1, Z_i=z] = \mathbbm{E}[Y_i(t)-Y_i(0)|A^{[t]}_{z,G_i} = 1, Z_i=z]  \nonumber \\
&= \sum_g P(G_i=g|A^{[t_1]}_{z,G_i} = 1)	\cdot \mathbbm{E}[Y_i(t)-Y_i(0)|G_i=g, Z_i=z] \nonumber \\
&= \frac{1}{P(T_i=t_1|Z_i=z)}\sum_{g:A^{[t_1]}_{zg}=1} P(G_i=g)\cdot \mathbbm{E}[Y_i(t)-Y_i(0)|G_i=g, Z_i=z] \nonumber \\
&= \frac{1}{P(T_i=t_1|Z_i=z)}\sum_{g} A^{[t_1]}_{zg}\cdot P(G_i=g)\cdot \{\mathbbm{E}[Y_i(t)|G_i=g, Z_i=z]-\mathbbm{E}[Y_i(0)|G_i=g, Z_i=z]\} \nonumber \\
&= \frac{1}{P(T_i=t_1|Z_i=z)} \sum_g A^{[t_1]}_{zg} \cdot \sum_{g'} P(G_i=g')\cdot \left\{[(A^{[t]})^+A^{[t]}]_{g,g'}\cdot \mathbbm{E}[Y(t)|G_i=g']\right. \nonumber \\
&\hspace{3.11in} \left.-[(A^{[0]})^+A^{[0]}]_{g,g'}\cdot \mathbbm{E}[Y(0)|G_i=g']\right\} \label{eq:nsogcondition}
\end{align}

where the third line uses that

\begin{align*}
P(A^{[t_1]}_{z,G_i} = 1) &= \sum_{g} P(G_i=g)\cdot A^{[t_1]}_{z,g} = \sum_{g} P(G_i=g)\cdot \mathbbm{1}(T_g(z)=t_1)\\
&= \sum_{g} P(G_i=g)\cdot P(T_i(z)=t_1|G_i=g,Z_i=z) = P(T_i=t_1|Z_i=z)
\end{align*}

by independence Eq. (<ref>). Defining $y^{[t]}$ to be a vector of $P(G_i=g) \cdot \mathbbm{E}[Y_i(t)|G_i=g]$ across all $g \in \mathcal{G}$, we can write the above as

$$\mathbbm{E}[Y_i(t)-Y_i(0)|T_i=t_1, Z_i=z] = \frac{[A^{[t_1]}\{(A^{[t]})^+A^{[t]} y^{[t]}-(A^{[0]})^+A^{[0]} y^{[0]}\}]_z}{[A^{[t_1]}P]_z}$$

under $\mathbf{G}^*=\mathbf{G}$, where $P$ is a vector of $P(G_i=g)$ for $g \in \mathcal{G}$ introduced previously in the proof of Theorem <ref>.

For NSOG to be satisfied, the RHS of Eq. (<ref>) must be equal to $\Delta_1$, regardless of the value of $t_1$ or $z$. This requires

$$ \frac{[A^{[t_1]}\{(A^{[t]})^+A^{[t]} y^{[t]}-(A^{[0]})^+A^{[0]} y^{[0]}\}]_z}{[A^{[t_1]}P]_z} = \Delta_t$$

for all $t_1$ and $z$. It is unclear if this equation can ever be satisfied unless $y^{[t]}_g=0$ for all $t$ and $g$ (in which case $\Delta_t=0$ and the equation is satisfied). Even if the strong condition that $\mathbbm{E}[Y_i(t)|G_i=g]$ be the same for all $g$ were to hold, then $y^{[t]} = \mathbbm{E}[Y_i(t)]\cdot P_g$ for all $g$ and
the equation above would reduce to

$$ \frac{[A^{[t_1]}(A^{[t]})^+A^{[t]}P]_z}{[A^{[t_1]}P]_z} \cdot \mathbbm{E}[Y_i(t)]-\frac{[A^{[t_1]}(A^{[0]})^+A^{[0]}P]_z}{[A^{[t_1]}P]_z} \cdot \mathbbm{E}[Y_i(0)] = \mathbbm{E}[Y_i(t)]-\mathbbm{E}[Y_i(0)]$$

This would be true if $A^{[t_1]} = A^{[t]} = A^{[0]}$, but that can never occur for $t \ne 0$ unless treatments $t$ or $0$ are never chosen by an selection group (since $[A^{[t]}+A^{[0]}]_{zg} \le 1$; a given individual cannot take two treatments under the same instrument value).


11institutetext: John Adams Institute for Accelerator Science, Blackett
Laboratory, Imperial College London, UK

# An Overview of Recent Progress in Laser Wakefield Acceleration Experiments

S.P.D. Mangles

###### Abstract

The goal of this paper is to examine experimental progress in laser wakefield
acceleration over the past decade (2004–2014), and to use trends in the data
to understand some of the important physical processes. By examining a set of
over 50 experiments, various trends concerning the relationship between plasma
density, accelerator length, laser power and the final electron beam energy
are revealed. The data suggest that current experiments are limited by
dephasing and that current experiments typically require some pulse evolution
to reach the trapping threshold.

Keywords
Laser wakefield accelerators; plasma accelerators; laser-plasma acceleration.

##  0.1 Introduction

This paper is a summary of a lecture given at the CERN Accelerator School 2015
on plasma-based wakefield accelerators. Its purpose is to provide an overview
of recent experimental progress in laser wakefield acceleration, concentrating
on the energy frontier and general trends that can be observed in the data
produced by the many groups around the world contributing to this growing
field. There are now more than 20 active laboratories performing laser
wakefield acceleration experiments. We will not detail the key results from
each of these laboratories in this paper, although we will examine data from
various published experiments and use the trends in the data to try to gain
some understanding of the underlying physical processes.

Laser wakefield accelerators were proposed 36 years ago in the seminal 1979
paper by Tajima and Dawson [1] and experiments in laser wakefield acceleration
were undertaken as soon as laser pulses of sufficiently short duration and
high power became available, thanks to the development of the laser technique
called ‘chirped pulse amplification’ [2]. Some of the early work in what we
call laser wakefield acceleration (where the pulse duration,
$\tau_{\mathrm{L}}$ is comparable to the plasma period
$2\pi/\omega_{\mathrm{p}}$) occurred in the 1990s (Refs.[3, 4]) using
picosecond glass lasers. The field passed a major milestone in the early 2000s
as high-power (${\sim}10\UW[T]$) femtosecond laser pulses, using titanium
sapphire laser systems, became available. One major result from that era
occurred when three groups from the UK, USA and France all demonstrated that
laser wakefield accelerators could produce electron beams with well-defined
energies [5, 6, 7]. The electron beams produced by these ${\sim}10\UW[T]$
lasers had energies of ${\sim}100\UMeV$ and were produced in plasmas of a
plasma density ${\sim}10^{19}\Ucm^{-3}$ that were only ${\sim}1\Umm$ long.

One of the key challenges that has driven progress in the field of laser
wakefield acceleration is the maximum achievable beam energy, and this paper
will concentrate on this challenge. There have, of course, been many other
significant areas of experimental progress including: improving beam
stability, especially by controlling injection (Refs.[8, 9]); diagnosis of
wakefield accelerators (Refs.[10, 11]) and the resulting improvements in our
understanding of the underlying processes; and the application of laser
wakefield accelerators for a range of applications, perhaps most notably their
use as novel sources of X-radiation (Refs.[12, 13, 14]). Details of the
progress in these areas are outside the scope of this paper.

Figure 1: Reported electron beam energies from laser wakefield experiments at
various laboratories over the last decade; data from Refs. [5, 6, 7, 15, 16,
17, 18, 19, 20, 21, 22, 8, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 9, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,
55, 56, 57, 58, 59, 60].

This paper will examine how the maximum achievable beam energy has progressed
over the last decade, and will use a set of published results [5, 6, 7, 15,
16, 17, 18, 19, 20, 21, 22, 8, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 9, 45, 46, 47, 48, 49, 50, 51, 52, 53,
54, 55, 56, 57, 58, 59, 60] to try to understand some of the physics behind
this progress. The progress in the maximum beam energy in laser wakefield
accelerator experiments has been rapid, as shown in energy-year, from a
maximum beam energy of $0.2\UGeV$ reported in 2002 [61] to the current record
of $4\UGeV$ from the group at the Lawrence Berkeley National Laboratory [58]
achieved in 2014 – an increase by a factor of 20 in just over a decade. It
should be noted that this is by no means an exhaustive list of all published
experiments in laser wakefield accelerators (there are just 52 publications in
this dataset, whereas a literature search for papers on ‘laser wakefield’ will
find over 1000 papers).

##  0.2 Overall trends in laser wakefield acceleration experiments

The rapid progress shown in energy-year is impressive. But how has it been
achieved? Over the same period of time, short-pulse (${\approx}30\Ufs$) laser
systems have become more powerful. Figure 2 shows that there is a clear trend:
higher-power lasers are capable of producing higher-energy electron beams.
However, these gains were not achieved by simply increasing the laser power,
the researchers behind these experiments have often found the optimum
conditions for their experiments. Key parameters involved in this optimization
include the operating plasma density, the length of the accelerator and the
laser intensity. This section will examine the data from various experiments
and compare them with predicted trends, to see whether they can confirm those
predictions and the underlying physical processes.

Figure 2: Variation of reported electron beam energy with laser power from
various experiments; data from Refs. [5, 6, 7, 15, 16, 17, 18, 19, 20, 21, 22,
8, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41,
42, 43, 44, 9, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59,
60].

###  0.2.1 Accelerator length and operating plasma density

Let us first consider the density of the plasma accelerator. The energy gained
by a particle of charge $q$ in an accelerating structure is simply
proportional to the product of the electric field and the length of the
accelerator, $d$, $W\simeq qEd$, where $E$ is the average accelerating
electric field experienced by the particle. One of the key physical
limitations in a laser wakefield accelerator is dephasing. The electrons
trapped in a wake are highly relativistic ($\gamma\gg 1$), so they travel at a
speed approaching that of light in vacuum ($v_{\mathrm{e}}\to c$), but the
phase speed of the wake is determined by the speed of the laser pulse that
drives the plasma wave. A simple expression for the speed of a laser pulse in
a plasma can be found by using the standard expression for the group velocity
of an electromagnetic wave in a plasma,

$\frac{v_{\mathrm{g}}}{c}=\sqrt{1-\frac{n_{\mathrm{e}}}{n_{\mathrm{c}}}}\approx
1-\frac{1}{2}\frac{n_{\mathrm{e}}}{n_{\mathrm{c}}}\ ,$ (1)

where $n_{\mathrm{e}}$ is the electron density of the plasma, $n_{\mathrm{c}}$
is the critical density for propagation of the electromagnetic (when the
plasma frequency $\omega_{\mathrm{p}}$ equals the frequency of the
electromagnetic wave, $\omega_{0}$) and it is assumed that $n_{\mathrm{e}}\ll
n_{\mathrm{c}}$. The wake’s phase speed is therefore slightly, but
significantly, less than $c$; crucially, the lower the plasma density, the
faster the phase velocity.

Because of this difference between the electron velocity and the wake phase
velocity, electrons in a laser-driven wake will outrun the wake111Note that
dephasing does not occur in wakefield accelerators driven by highly
relativistic charged particle beams as both the accelerated and driver beams
are highly relativistic.. If the electron is injected at the start of the
accelerating phase of the plasma wave and then outruns the wave by half a
plasma wavelength ($\lambda_{\mathrm{p}}/2=\pi c/\omega_{\mathrm{p}}$), it can
no longer gain energy from the plasma wave. If the electron has an initial
velocity $v_{\mathrm{e}}=\beta_{\mathrm{e}}c$ and the plasma wave has a phase
velocity $v_{\phi}=\beta_{\phi}c$, then the time it takes for this to occur is
$t_{\mathrm{d}}=\lambda_{\mathrm{p}}/(2c(\beta_{\mathrm{e}}-\beta_{\phi}))$.
The dephasing length is then the distance that the electron travels in this
time. Since $\beta_{\mathrm{e}}\to 1$ and $\beta_{\phi}\simeq
1-\frac{1}{2}({n_{\mathrm{e}}}/{n_{\mathrm{c}}})$, this reduces to

$L_{\rm
dephasing}\simeq\frac{n_{\mathrm{c}}}{n_{\mathrm{e}}}\lambda_{\mathrm{p}}\propto
n_{\mathrm{e}}^{-\frac{3}{2}}\ .$ (2)

It is interesting to see how the lengths of the accelerators in the set of
experiments vary, and how this compares with what we might expect if dephasing
is important. These data are shown in density-length. The top panel shows how
the reported electron beam energy varies with the length of the wakefield
accelerator. There is a clear correlation – the higher electron energies are
achieved with longer accelerators, as we might expect. The bottom panel of
density-length shows how the length of the accelerator and the plasma density
at which it was operating are related. The line on this curve is the simple
expression for the dephasing length (dephasing_length).

Figure 3: Top: Variation of reported electron beam energy with accelerator
length. Bottom: Relationship between operating plasma density and accelerator
length. The line shows the expression for the dephasing length
(dephasing_length). Data from Refs. [5, 6, 7, 15, 16, 17, 18, 19, 20, 21, 22,
8, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41,
42, 43, 44, 9, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59,
60].

The maximum electric field that a plasma wave can support increases with
plasma density, since it scales as

$E_{\mathrm{max}}\simeq
m_{\mathrm{c}}\omega_{\mathrm{p}}/e\propto\sqrt{n_{\mathrm{e}}}\ .$ (3)

The maximum energy that can be gained by an electron in a plasma wave as a
function of plasma density is therefore expected to be

$W(n_{\mathrm{e}})\simeq
E_{\max}L_{\mathrm{dephasing}}\propto\frac{1}{n_{\mathrm{e}}}\ .$ (4)

Figure 4 shows how the plasma density and beam energy vary in the set of
experiments. The line on energy-density is simply
$W/(m_{\mathrm{e}}c^{2})=\kappa\,n_{\mathrm{c}}/n_{\mathrm{e}}$ and shows good
agreement with the entire dataset for $\kappa=1$. The scaling laws in
Ref.[62], by Wei Lu et al., for the blow-out or ‘bubble’ regime of wakefield
accelerators, suggest that the scaling law should be

$W(n_{\mathrm{e}},a_{0})\simeq\frac{2}{3}a_{0}\frac{n_{\mathrm{c}}}{n_{\mathrm{e}}}m_{\mathrm{e}}c^{2}\propto\frac{a_{0}}{n_{\mathrm{e}}}\
,$ (5)

where $a_{0}=eA_{0}/(m_{\mathrm{e}}c)=eE_{0}/(m_{\mathrm{e}}\omega_{0}c)$ is
the normalized peak vector potential (or strength parameter) of the laser
pulse. This scaling predicts that the beam energy should not only be
proportional to $1/n_{\mathrm{e}}$ but also proportional to the laser
strength, $a_{0}$. The experiments shown correspond to a wide range of initial
laser intensities (corresponding to $a_{0}=0.5$–$4.0$), yet they do not appear
to show a dependence on $a_{0}$.

Figure 4: Variation of reported electron beam energy with the density in the
accelerator. The line shows the relation
$W/(m_{\mathrm{e}}c^{2})=\kappa\,n_{\mathrm{c}}/n_{\mathrm{e}}$ with
$\kappa=1$. Data from Refs. [5, 6, 7, 15, 16, 17, 18, 19, 20, 21, 22, 8, 23,
24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42,
43, 44, 9, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60].

One possible reason for this apparent discrepancy is that the initial value
for $a_{0}$ is not the value of $a_{0}$ that determines the wake amplitude. It
is well known that laser pulses can undergo significant evolution once they
enter the plasma. The processes of self-focusing, self-compression and photon
acceleration[63] can all act to change $a_{0}$ as the pulse propagates.
Together, these processes can be termed the ‘self-evolution’ of the laser
pulse. One interpretation of the experimental data is, therefore, that the
process of self-evolution has occurred until $a_{0}\approx 3$ for all of the
data shown. Why would this be the case? One reasonable hypothesis is that each
point in the dataset corresponds to the maximum energy achieved during a
particular experiment and that this will occur at (or at least close to) the
lowest density at which a particular experiment can trap and accelerate
electrons. Many experiments operate by fixing the laser power and plasma
length while varying the plasma density, for reasons of experimental
simplicity. When an experiment is conducted in this manner, there will be a
minimum density at which electron beams are trapped and accelerated (the
trapping threshold). Because self-evolution happens less quickly and severely
at lower densities, the maximum $a_{0}$ that is reached inside the accelerator
will decrease with decreasing plasma density. Therefore, the maximum achieved
electron beam energy will correspond to the minimum laser strength required to
produce trapping. The fact that the experimental dataset matches the non-
linear wakefield scaling but only if $a_{0}\approx 3$ suggests that the
minimum $a_{0}$ required for trapping is $a_{0}\approx 3$.

###  0.2.2 Laser spot size and matched guiding

The experimental data clearly show that higher-power lasers are required to
achieve higher electron beam energies. But what are the physical processes
behind this trend? It was argued that the experimental trends are consistent
with there being a minimum value of $a_{0}\approx 3$, which is needed for
trapping, and that this value is reached because of the way the pulse evolves
as it propagates.

Consider the relationship between the intensity, $I$ ($\propto a_{0}^{2}$),
and power, $P$, of a laser pulse. Since $I=P/A$, where $A$ ($\propto w^{2}$,
the laser spot size) is the focal spot area, we have

$P\propto a_{0}^{2}w^{2}\,.$ (6)

The fact that higher-power lasers are needed to reach the $a_{0}\approx 3$
threshold at lower densities therefore implies that the spot size that these
laser pulses produce after evolution is larger. Pulse evolution is a result of
the feedback between the refractive index gradient associated with the plasma
wave and the laser pulse, it is mediated by the plasma itself. Lower-density
plasmas, therefore, have a lesser effect on the laser pulse – resulting in
slower evolution – and crucially this affects the properties that the pulse
obtains as a result of self-evolution.

One important concept that arises from this is that of the _matched_ spot size
– one where the self-focusing caused by the plasma balances the natural
diffraction of the laser pulse and stable propagation occurs. In the blow-out
regime (where the laser pulse expels practically all the electrons from inside
the bubble), the transverse density profile of the bubble is approximately
zero and flat inside the bubble, with very steep walls at the edges. The
refractive index of an underdense plasma is $\eta\approx
1-n_{\mathrm{e}}/(2n_{\mathrm{c}})$, so the refractive index profile of this
bubble is similar to that of a single-mode optical fibre, we have a ‘core’ of
a certain diameter (the bubble diameter) surrounded by ‘cladding’ (the bubble
sheath) of a lower refractive index. The guided mode in such an optical fibre
has a transverse size that is approximately equal to the size of the core, so
we expect that in a laser wakefield accelerator we will get stable propagation
(no spot size oscillations) when the transverse size of the laser pulse is
approximately equal to that of the bubble.

Of course, one main difference between the bubble and an optical fibre is that
the size of the bubble is determined by the properties of the laser pulse
itself. If the laser spot is too small, then it will drive a bubble that is
larger than the laser spot. This over-sized bubble will support a larger mode
and the laser pulse will expand, which in turn reduces the bubble size. The
radius of the bubble can be found by balancing the ponderomotive force of the
laser pulse with the force due to the electric field of the bubble. If the
bubble is approximately the same size as the laser pulse then it will be
matched. The size of the matched spot can be found by balancing the forces on
an electron at the edge of the bubble (the laser’s ponderomotive force and the
force due to the electric field inside the bubble). An expression for the
matched spot size, $w_{\mathrm{m}}$ is

$w_{\mathrm{m}}\simeq 2\frac{c}{\omega_{\mathrm{p}}}\sqrt{a_{0}}\ ,$ (7)

where the numerical factor of two was found through particle-in-cell
simulations[62]. This expression is particularly useful for finding the
correct initial parameters of the accelerator, for a given laser system, one
should first determine the spot size at which the threshold $a_{0}\simeq 3$ is
reached; this expression can then be used to determine the correct operating
plasma density. However, when self-focusing plays an important role, we need
an expression for the matched spot size that depends not on the initial laser
intensity but rather the value of $a_{0}$ that it reaches after self-focusing.
We can use the fact that the ratio of the laser power, $P_{\mathrm{L}}$, to
the critical power for self-focusing, $P_{\mathrm{c}}$, can be written as

$\frac{P_{\mathrm{L}}}{P_{\mathrm{c}}}=\frac{1}{32}\frac{\omega_{p}^{2}}{c^{2}}\,{a_{0}^{2}w^{2}}\
,$ (8)

and the fact that $a_{0}^{2}w^{2}$ is constant during focusing (assuming that
self-focusing happens more rapidly than any pulse compression) to eliminate
$a_{0}$ from matched1. This results in the following expression for the
matched spot size:

$w_{\mathrm{m}}\simeq
2\sqrt{2}\frac{c}{\omega_{\mathrm{p}}}\left(\frac{P_{\mathrm{L}}}{P_{\mathrm{c}}}\right)^{\frac{1}{6}}\
.$ (9)

Note that $2w_{\mathrm{m}}\approx\lambda_{\mathrm{p}}$, as long as
$P_{\mathrm{L}}$ is not many times greater than $P_{\mathrm{c}}$.

It is interesting to examine how the spot size used in experiments compares
with this matched spot size. Figure 5 shows the variation in the initial
(vacuum) laser focal spot size with the operating plasma density in the
experiments in the dataset. There is a clear overall trend towards larger
initial spots at lower densities (and therefore higher electron beam
energies), and the spot sizes used are reasonably close to
$\lambda_{\mathrm{p}}$ (shown as a solid blue line in spot_size). However,
given that the evidence suggests that most of these experiments reached
$a_{0}\approx 3$ as a result of pulse evolution, it is also interesting to
compare the initial spot size with the expected matched spot size for
$a_{0}\approx 3$ (this is shown as a dashed red line in spot_size). Most of
the experiments are clearly operating at an initial spot size significantly
larger than this matched spot size, and none of the selected experiments
operate with a spot size below this. This suggests either that most
experiments are operating with too large an initial spot size, wasting
accelerator length and laser energy while the laser pulse self-focuses, or
that there is some experimental advantage in starting at a spot size larger
than the matched spot size and letting pulse evolution happen.

Figure 5: Variation of initial laser focal spot size with operating plasma
density in various laser wakefield acceleration experiments. Solid blue line,
plasma wavelength, $\lambda_{\mathrm{p}}$; red dashed line, matched spot size
assuming $a_{0}\approx 3$. Data from Refs. [5, 6, 7, 15, 16, 17, 18, 19, 20,
21, 22, 8, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39,
40, 41, 42, 43, 44, 9, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58,
59, 60].

###  0.2.3 To guide or not to guide?

It is well known that a tightly focused laser pulse will quickly diffract in a
vacuum. Effectively, the pulse can only remain intense over a distance of
about a Rayleigh length, $z_{\mathrm{R}}=\pi w_{0}^{2}/\lambda$. The
accelerator lengths used in laser wakefield experiments are typically much
longer than this and some sort of ‘guiding’ is therefore needed to keep the
laser intensity sufficiently high to drive a wake throughout the structure.
There are two principal techniques to achieve this guiding, which both rely on
creating a waveguide structure to counteract diffraction. This requires that
the transverse plasma density profile has a minimum on-axis. Such a density
minimum is naturally created by the laser pulse itself – that is, the bubble
itself acts as a waveguide. Alternatively, a preformed density minimum can be
formed, for example using a capillary discharge [24]; using a preformed
waveguide brings significant complexity to an experiment and restricts
diagnostic access (for wake imaging diagnostics based on ultrashort probes
[11]). However, the bubble cannot self-guide the very front slice of the laser
pulse, since the density minimum is not formed immediately. Because of this,
external channels are expected to be more efficient, they allow the laser to
propagate over greater distances at high $a_{0}$. With these points in mind,
it is interesting to see whether the experimental evidence supports the use of
external waveguides. In guiding_type, the data are sorted into self-guided and
externally guided experiments. The top panel of guiding_type appears to show
that there is no real advantage to using external guiding structures; the
electron beam energy, as a function of plasma density, follows the same trend
for both subsets of the data, they are both limited by dephasing. However, the
bottom panel of guiding_type reveals the distinct advantage that experiments
in externally guided structures have over self-guided ones. The highest
electron energy achieved at a given laser power is almost always from an
externally guided experiment; the self-guided experiments at the same power
tend to produce lower energy electron beams.

Figure 6: Variation of reported electron beam energy from various experiments.
Top: Variation of beam energy with plasma density. Bottom: Variation of beam
energy with laser power. Filled circles, experiments in preformed guiding
structures; open squares, experiments without preformed guiding structures.
Data from Refs. [5, 6, 7, 15, 16, 17, 18, 19, 20, 21, 22, 8, 23, 24, 25, 26,
27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 9, 45,
46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60].

Are these data sufficient to suggest that external guiding structures improve
the efficiency of laser wakefield accelerators, owing to reduced energy losses
as the pulse propagates? Or is something more subtle occurring? In one
experiment with an external waveguide using a ${\simeq}20\UW[T]$ laser [30],
differences between the plasma density measured during a high-intensity laser
shot and plasma density measurements made offline suggested that, under the
conditions in which electron beams were observed, the laser pulse caused
additional ionization of the plasma. As the plasma used in that experiment was
formed from hydrogen, which is readily ionized, this suggested that high-$Z$
impurities from the walls of the capillary were being ionized by the main
pulse and injected into the plasma wave. Since that work, a number of groups
have proceeded to exploit this ionization injection mechanism to reduce the
threshold for injection and increase the beam charge (Refs. [36, 39]). There
are also various other injection techniques, including injection due to
propagation in a density gradient [35, 9] and colliding pulse injection [8].
In injection_type, the variation of electron beam energy with laser power is
plotted again, this time with the data divided into subsets based on injection
type. As it is not known whether the majority of capillary-discharge-based
experiments rely on self-injection or whether ionization injection plays a
role, as in Ref. [30], these experiments have been placed in their own subset.
Some of the ionization injection experiments also produce higher electron beam
energies for a given laser power than self-injection experiments, and lie on
the upper curve of the entire dataset, just as the capillary discharge dataset
does. This might suggest that the injection mechanism is more important than
the guiding mechanism in determining the maximum energy that can be achieved
from a given laser power. This makes sense, as the value that $a_{0}$ reaches
after self-evolution decreases with decreasing laser power, but alternative
injection mechanisms should lower the value that $a_{0}$ needs to reach in
order for injection to occur. However, injection mechanisms other than
ionization injection seem to perform similarly to self-injection in terms of
the electron energy that can be obtained for a given laser power. At present,
the evidence is still inconclusive; this is clearly a matter that requires
further study.

Figure 7: Variation of reported electron beam energy from various experiments
as a function of laser power and for different injection mechanisms. Filled
black circles, capillary discharge experiments; black squares, ionization
injection; green diamonds, density down-ramp injection; cyan triangles,
colliding pulse injection; open black squares, self-injection. Data from Refs.
[5, 6, 7, 15, 16, 17, 18, 19, 20, 21, 22, 8, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 9, 45, 46, 47, 48, 49,
50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60].

##  0.3 Future directions

Having examined the trends in experimental data from laser wakefield
experiments over the last decade, we find that there is one overriding
message. To keep pushing the electron energy achievable from a single stage of
a laser wakefield to ever higher values clearly requires operation at lower
densities and over longer distances, and such experiments will need more
powerful lasers. However, it should be noted that laser power can be increased
in two ways, either by increasing the laser energy or by decreasing the pulse
duration. Most of the data presented were obtained using laser pulses with a
duration of ${\sim}30\Ufs$, so the trends in laser power are dominated by the
pulse energy. Simulations demonstrate that the pulse length in a laser
wakefield accelerator should not be too short. Energy is lost predominantly
from the front of the laser pulse (this is the part of the laser that is doing
work on the plasma), so the shorter the laser pulse, the more quickly it runs
out of energy. As new laser facilities come online around the world, it is
worth remembering that the shortest pulse is not necessarily the best for
driving laser wakefield accelerators. Pulse durations where
$c\tau_{\mathrm{L}}$ is comparable to the plasma wavelength,
$\lambda_{\mathrm{p}}$, drive efficient wakes and allow the pump depletion
length to match the dephasing length [62]. Pulses significantly shorter than
$\lambda_{\mathrm{p}}/2$ expend their energy before the maximum energy can be
reached.

The experimental trends show that most experiments are attaining electron beam
energies that are limited by dephasing. If the community is to be able to keep
up the impressive pace it has sustained over the last decade then methods to
overcome this limitation will become increasingly important, especially as
further increases in laser power become ever more expensive. As a result,
techniques to overcome dephasing, such as quasi-phase matching [64], staging
[65] and density tapering [66] will all become important areas of research.

The final point that we would like to make is that the laser systems currently
used to drive laser wakefield accelerators are woefully inefficient. A
titanium sapphire laser has a ‘wall plug’ efficiency of ${\sim}0.1\%$. They
also currently run at very low repetition rates compared with a conventional
accelerator. As laser wakefield accelerators push the energy frontier and
strive to become workhorses for applications, it will become increasingly
important that both repetition rate and efficiency are properly considered and
significantly improved. Much more efficient laser architectures are available,
including thin-disc [67] and fibre [68] lasers, that can also readily operate
at much higher repetition rates. However, these systems do not yet have the
capability to produce laser pulses with the high energy and short pulse
duration needed to drive current laser wakefield acceleration experiments.
Innovative solutions, including the use of coherent [69] and incoherent [70]
combinations of many low-energy laser pulses to drive wakefield experiments,
or the use of trains of many low-energy pulses to resonantly drive a high
amplitude wakefield [71], may well form the most promising routes to high
repetition rate, high-efficiency laser wakefield accelerators suitable for
particle physics experiments or light-source based applications.

## References

  * [1] T. Tajima and J.M. Dawson,  Phys. Rev. Lett. 43(4) (1979) 267. http://dx.doi.org/10.1103/PhysRevLett.43.267
  * [2] D. Strickland and G. Mourou,  Opt. Commun. 55(6) (1985), 447–449. http://dx.doi.org/10.1016/0030-4018(85)90151-8
  * [3] K. Nakajima et al.,  Phys. Scr. T52 (1994) 61. http://dx.doi.org/10.1088/0031-8949/1994/T52/009
  * [4] F. Amiranoff et al.,  Phys. Rev. Lett. 81(5) (1998) 995. http://dx.doi.org/10.1103/PhysRevLett.81.995
  * [5] S.P.D. Mangles et al.,  Nature 431 (2004) 535–538. http://dx.doi.org/10.1038/nature02939
  * [6] C.G.R. Geddes et al.,  Nature 431 (2004) 538–541. http://dx.doi.org/10.1038/nature02900
  * [7] J. Faure et al.,  Nature 431 (2004) 541–544. http://dx.doi.org/10.1038/nature02963
  * [8] J. Faure et al.,  Nature 444 (2006) 737–739. http://dx.doi.org/10.1038/nature05393
  * [9] A.J. Gonsalves et al.,  Nat. Phys. 7 (2011) 862–866. http://dx.doi.org/10.1038/nphys2071
  * [10] N.H. Matlis et al.,  Nat. Phys. 2 (2006) 749–753. http://dx.doi.org/10.1038/nphys442
  * [11] A. Sävert et al.,  Phys. Rev. Lett. 115(5) (2015) 055002. http://dx.doi.org/10.1103/PhysRevLett.115.055002
  * [12] A. Rousse et al.,  Phys. Rev. Lett. 93(13) (2004) 135005. http://dx.doi.org/10.1103/PhysRevLett.93.135005
  * [13] S. Kneip et al.,  Nat. Phys. 6 (2010) 980–983. http://dx.doi.org/10.1038/nphys1789
  * [14] M. Fuchs et al.,  Nat. Phys. 5 (2009) 826–829. http://dx.doi.org/10.1038/nphys1404
  * [15] E. Miura et al.,  Appl. Phys. Lett. 86(25) (2005) 251501. http://dx.doi.org/10.1063/1.1949289
  * [16] H. Kotaki et al.,  Laser Phys. 16(7) (2006) 1107–1110. http://dx.doi.org/10.1134/S1054660X06070140
  * [17] M. Mori et al.,  Phys. Lett. A 356(2) (2006) 146–151. http://dx.doi.org/10.1016/j.physleta.2006.06.001
  * [18] S. Masuda et al.,  J. Phys. IV France 133 (2006) 1127–1129. http://dx.doi.org/10.1051/jp4:2006133229
  * [19] C.-T. Hsieh et al.,  Phys. Rev. Lett. 96(9) (2006) 095001. http://dx.doi.org/10.1103/PhysRevLett.96.095001
  * [20] B. Hidding et al.,  Phys. Rev. Lett. 96(10) (2006) 105004. http://dx.doi.org/10.1103/PhysRevLett.96.105004
  * [21] T. Hosokai et al.,  Phys. Rev. E 73(3) (2006) 036407. http://dx.doi.org/10.1103/PhysRevE.73.036407
  * [22] S.P.D. Mangles et al.,  Phys. Rev. Lett. 96(21) (2006) 215001. http://dx.doi.org/10.1103/PhysRevLett.96.215001
  * [23] S.A. Reed et al.,  Appl. Phys. Lett. 89(23) (2006) 231107. http://dx.doi.org/10.1063/1.2400400
  * [24] W.P. Leemans et al.,  Nat. Phys. 2 (2006) 696–699. http://dx.doi.org/10.1038/nphys418
  * [25] S. Masuda et al.,  Phys. Plasmas 14(2) (2007) 023103. http://dx.doi.org/10.1063/1.2434248
  * [26] T. Ohkubo et al.,  Phys. Rev. ST Accel. Beams 10(3) (2007) 031301. http://dx.doi.org/10.1103/PhysRevSTAB.10.031301
  * [27] S. Karsch et al.,  New J. Phys. 9 (2007) 415. http://dx.doi.org/10.1088/1367-2630/9/11/415
  * [28] S.P.D. Mangles et al.,  Phys. Plasmas 14(5) (2007) 056702. http://dx.doi.org/10.1063/1.2436481
  * [29] A. Gamucci et al.,  IEEE Trans. Plasma Sci. 36(4) (2008) 1699–1706. http://dx.doi.org/10.1109/TPS.2008.2000898
  * [30] T.P. Rowlands-Rees et al.,  Phys. Rev. Lett. 100(10) (2008) 105005. http://dx.doi.org/10.1103/PhysRevLett.100.105005
  * [31] N.A.M. Hafz et al.,  Nat. Photonics 2 (2008) 571. http://dx.doi.org/10.1038/nphoton.2008.155
  * [32] K. Schmid et al.,  Phys. Rev. Lett. 102(12) (2009) 124801. http://dx.doi.org/10.1103/PhysRevLett.102.124801
  * [33] S. Kneip et al.,  Phys. Rev. Lett. 103(3) (2009) 035002. http://dx.doi.org/10.1103/PhysRevLett.103.035002
  * [34] D.H. Froula et al.,  Phys. Rev. Lett. 103(21) (2009) 215006. http://dx.doi.org/10.1103/PhysRevLett.103.215006
  * [35] K. Schmid et al.,  Phys. Rev. ST Accel. Beams 13(9) (2010) 091301. http://dx.doi.org/10.1103/PhysRevSTAB.13.091301
  * [36] A.E. Pak et al.,  Phys. Rev. Lett. 104(2) (2010) 025003. http://dx.doi.org/10.1103/PhysRevLett.104.025003
  * [37] T.P.A. Ibbotson et al.,  New J. Phys. 12 (2010) 45008. http://dx.doi.org/10.1088/1367-2630/12/4/045008
  * [38] C.E. Clayton et al.,  Phys. Rev. Lett. 105(10) (2010) 105003. http://dx.doi.org/10.1103/PhysRevLett.105.105003
  * [39] C. McGuffey et al.,  Phys. Rev. Lett. 104(2) (2010) 025004. http://dx.doi.org/10.1103/PhysRevLett.104.025004
  * [40] H. Lu et al.,  Appl. Phys. Lett. 99(9) (2011) 091502. http://dx.doi.org/10.1063/1.3626042
  * [41] S. Fourmaux et al.,  New J. Phys. 13 (2011) 033017. http://dx.doi.org/10.1088/1367-2630/13/3/033017
  * [42] J.S. Liu et al.,  Phys. Rev. Lett. 107(3) (2011) 035001. http://dx.doi.org/10.1103/PhysRevLett.107.035001
  * [43] B.B. Pollock et al.,  Phys. Rev. Lett. 107(4) (2011) 045001. http://dx.doi.org/10.1103/PhysRevLett.107.045001
  * [44] O. Lundh et al.,  Nat. Phys. 7 (2011) 219–222. http://dx.doi.org/10.1038/nphys1872
  * [45] P. Brijesh et al.,  Phys. Plasmas 19(6) (2012) 063104. http://dx.doi.org/10.1063/1.4725421
  * [46] M.Z. Mo et al.,  Appl. Phys. Lett. 100(7) (2012) 074101. http://dx.doi.org/10.1063/1.3685464
  * [47] S. Kneip et al.,  Phys. Rev. ST Accel. Beams 15(2) (2012) 021302. http://dx.doi.org/10.1103/PhysRevSTAB.15.021302
  * [48] R. Weingartner et al.,  Phys. Rev. ST Accel. Beams 15(11) (2012) 111302. http://dx.doi.org/10.1103/PhysRevSTAB.15.111302
  * [49] F. Albert et al.,  Phys. Rev. Lett. 111(23) (2013) 235004. http://dx.doi.org/10.1103/PhysRevLett.111.235004
  * [50] P.A. Walker et al.,  New J. Phys. 15 (2013) 045024. http://dx.doi.org/10.1088/1367-2630/15/4/045024
  * [51] M.Z. Mo et al.,  Appl. Phys. Lett. 102(13) (2013) 134102. http://dx.doi.org/10.1063/1.4799280
  * [52] S. Corde et al.,  Nat. Commun. 4 (1501) (2013) 1309.6364v1. http://dx.doi.org/10.1038/ncomms2528
  * [53] S. Chen et al.,  Phys. Rev. Lett. 110(15) (2013) 155003. http://dx.doi.org/10.1103/PhysRevLett.110.155003
  * [54] H.T. Kim et al.,  Phys. Rev. Lett. 111(16) (2013) 165002. http://dx.doi.org/10.1103/PhysRevLett.111.165002
  * [55] X. Wang et al.,  Nat. Commun. 4(1988) (2013). http://dx.doi.org/10.1038/ncomms2988
  * [56] G. Sarri et al.,  Phys. Rev. Lett. 113 (2014) 224801. http://dx.doi.org/10.1103/PhysRevLett.113.224801
  * [57] N.D. Powers et al.,  Nat. Photonics 8 (2014) 28–31. http://dx.doi.org/10.1038/nphoton.2013.314
  * [58] W.P. Leemans et al.,  Phys. Rev. Lett. 113(24) (2014) 245002. http://dx.doi.org/10.1103/PhysRevLett.113.245002
  * [59] K. Khrennikov et al.,  Phys. Rev. Lett. 114(19) (2015) 195003. http://dx.doi.org/10.1103/PhysRevLett.114.195003
  * [60] M. Schnell et al.,  J. Plasma Phys. 81(04) (2015) 475810401. http://dx.doi.org/10.1017/S0022377815000379
  * [61] V. Malka et al.,  Science 298 (5598) (2002) 1596–1600. http://dx.doi.org/10.1126/science.1076782
  * [62] W. Lu et al.,  Phys. Rev. ST Accel. Beams 10(6) (2007) 061301. http://dx.doi.org/10.1103/PhysRevSTAB.10.061301
  * [63] W.B. Mori,  IEEE J. Quant. Electron. 33(11) (1997) 1942–1953. http://dx.doi.org/10.1109/3.641309
  * [64] S.J. Yoon et al.,  Phys. Rev. ST Accel. Beams 15(8) (2012) 081305. http://dx.doi.org/10.1103/PhysRevSTAB.15.081305
  * [65] D. Kaganovich et al.,  Phys. Plasmas 12(10) (2005) 100702. http://dx.doi.org/10.1063/1.2102727
  * [66] P. Sprangle et al.,  Phys. Rev. E 63(5) (2001) 056405. http://dx.doi.org/10.1103/PhysRevE.63.056405
  * [67] A. Giesen and J. Speiser,  IEEE J. Sel. Top. Quant. Electron. 13(3) (2007) 598. http://dx.doi.org/10.1109/JSTQE.2007.897180
  * [68] C. Jauregui, J. Limpert and A. Tünnermann,  Nat. Photonics 7 (2013) 861–867. http://dx.doi.org/10.1038/nphoton.2013.273
  * [69] G. Mourou et al.,  Nat. Photonics 7 (2013) 258–261. http://dx.doi.org/10.1038/nphoton.2013.75
  * [70] C. Benedetti, et al.,  Phys. Plasmas 21(5) (2014) 056706\. http://dx.doi.org/10.1063/1.4878620
  * [71] S.M. Hooker et al.,  J. Phys. B: At. Mol. Opt. Phys 47(23) (2014) 234003. http://dx.doi.org/10.1088/0953-4075/47/23/234003