$[0,L]\times\cdots\times[0,L]\setminus\bigcup_{i\in[d]}J_{i}$, the zero
Lebesgue measure sets will be mapped to zero Lebesgue measure sets by $f_{I}$.
Therefore,

$\mu_{L}\left(f_{I}\left(S\cap([0,L]^{d}\setminus\cup_{i\in[d]}J_{i})\right)\right)=0.$

Consequently,

$\displaystyle\mu_{L}\left(S^{\frac{1}{2}}_{I}\cap\left[-\sqrt{L},\sqrt{L}\right]\right)$
$\displaystyle=\mu_{L}\left(f_{I}\left(S\cap([0,L]^{d}\setminus\cup_{i\in[d]}J_{i})\right)\right)+\mu_{L}\left(f_{I}\left(S\cap(\cup_{i\in[d]}J_{i})\right)\right)$
$\displaystyle\leq 0+d\frac{\varepsilon}{d}=\varepsilon.$

Since $\varepsilon$ is arbitrary,
$\mu_{L}(S^{\frac{1}{2}}_{I}\cap\left[-\sqrt{L},\sqrt{L}\right]^{d})=0$,
completing the proof.

Using this observation, we can prove the following useful lemma.

###### Lemma 15

If there are $Q\subseteq\mathbb{R}^{d}$ and $c\in\mathbb{R}$, such that
$\mu_{L}(Q)>0$ (Lebesgue measure) and

$\|\boldsymbol{x}\|^{2}_{{\boldsymbol{M}}_{1}}=\|\boldsymbol{x}\|^{2}_{{\boldsymbol{M}}_{2}}+c\quad\quad\text{for
each $\boldsymbol{x}\in Q$},$

then ${\boldsymbol{M}}_{1}={\boldsymbol{M}}_{2}$ and $c=0$.

Proof  For each $\boldsymbol{z}\in Q$, we clearly have
$\boldsymbol{x}_{i}^{t}({\boldsymbol{M}}_{1}-{\boldsymbol{M}}_{2})\boldsymbol{x}_{i}=c$.
Since ${\boldsymbol{M}}={\boldsymbol{M}}_{1}-{\boldsymbol{M}}_{2}$ is a real
value symmetric matric,
${\boldsymbol{M}}={\boldsymbol{U}}^{t}{\boldsymbol{D}}{\boldsymbol{U}}$ where
${\boldsymbol{U}}$ is an orthonormal matrix and ${\boldsymbol{D}}$ is a
diagonal $d\times d$ matrix whose $(i,i)$ element is $a_{i}$. To prove the
assertion, it suffices to show that ${\boldsymbol{D}}=0$. For a contradiction,
suppose that ${\boldsymbol{D}}\neq 0$. Set
$Q^{\prime}=\\{{\boldsymbol{U}}\boldsymbol{x}:\boldsymbol{x}\in Q\\}$ and
$S=\left\\{\boldsymbol{y}\in\mathbb{R}_{\geq
0}^{{}^{d}}:\langle\boldsymbol{y},(a_{1},\ldots,a_{d})\rangle=c\right\\}$. For
each $\boldsymbol{z}=(z_{1},\ldots,z_{d})\in Q^{\prime}$,

$\displaystyle\sum_{i=1}^{d}z_{i}^{2}a_{i}$
$\displaystyle=\boldsymbol{z}^{t}{\boldsymbol{D}}\boldsymbol{z}$
$\displaystyle=\boldsymbol{x}^{t}{\boldsymbol{U}}^{t}{\boldsymbol{D}}{\boldsymbol{U}}\boldsymbol{x}$
$\displaystyle=\boldsymbol{x}^{t}{\boldsymbol{M}}\boldsymbol{x}=c,$

which concludes that

$Q^{\prime}\subseteq\left\\{\boldsymbol{z}\in\mathbb{R}^{d}:\langle\boldsymbol{z}^{2},(a_{1},\ldots,a_{d})\rangle=c\right\\}=S^{\frac{1}{2}}.$

Using Observation 2, since we know $\mu_{L}(S)=0$ we obtain
$\mu_{L}(S^{\frac{1}{2}})=0$ and thus $\mu_{L}(Q^{\prime})=0$. On the other
hand, $\mu_{L}(Q)=\mu_{L}(Q^{\prime})$ which implies $\mu_{L}(Q)=0$, a
contradiction.

##  Appendix C $\varepsilon$-cover ($\varepsilon$-net) for
$\mathcal{M}\times[0,B]$

As we are going to prove a uniform convergence theorem between empirical and
true losses, we need to define the $\varepsilon$-cover of a metric space. For
a metric space $(\mathcal{X},d)$, an $\varepsilon$-cover $\mathcal{E}$ is a
subset of $\mathcal{X}$ such that for each $x\in\mathcal{X}$, there is some
$y\in\mathcal{E}$ with $d(x,y)<\varepsilon$. In the following, we introduce an
$\varepsilon$-cover for $\mathcal{M}\times[0,B]$. However, we should first
define a metric over $\mathcal{X}=\mathcal{M}\times[0,B]$. For each
$({\boldsymbol{M}}_{1},\tau_{1}),({\boldsymbol{M}}_{2},\tau_{2})\in\mathcal{X}$,
we define

${\rm
d}\left(({\boldsymbol{M}}_{1},\tau_{1}),({\boldsymbol{M}}_{2},\tau_{2})\right)=\|{\boldsymbol{M}}_{1}-{\boldsymbol{M}}_{2}\|_{2}+|\tau_{1}-\tau_{2}|.$

###### Lemma 16

There exists an $\varepsilon$-cover $\mathcal{E}$ of $\mathcal{M}\times[0,B]$
under metric ${\rm d}$ of size

$\frac{B}{\varepsilon}\left(\frac{4\beta
d\sqrt{d}}{\varepsilon}\right)^{d^{2}}.$

Proof  The inequality

$\|{\boldsymbol{M}}\|_{F}\leq\sqrt{d}\|{\boldsymbol{M}}\|_{2}\leq\beta\sqrt{d}$

indicates that $\mathcal{M}$ is a subset of the $d^{2}$ dimensional Euclidean
ball of the radius $\beta\sqrt{d}$ centered at the origin, i.e.,

$\mathcal{M}=\left\\{{\boldsymbol{M}}_{d\times d}\colon{\boldsymbol{M}}\text{
is p.s.d. and
}\|{\boldsymbol{M}}\|_{2}\leq\beta\right\\}\subseteq\\{{\boldsymbol{M}}\in\mathbb{R}^{d^{2}}\colon\|{\boldsymbol{M}}\|_{F}\leq\beta\sqrt{d}\\}.$

It is known that a $k$-dimentional Euclidean ball of radius $r$ can be covered
by at most $\left(\frac{2r\sqrt{k}}{\varepsilon}\right)^{k}$ number of balls
of radius $\varepsilon$. So, $\mathcal{M}$ has an
$\frac{\varepsilon}{2}$-cover $\mathcal{E}_{2}$ of size at most
$\left(\frac{4\beta d\sqrt{d}}{\varepsilon}\right)^{d^{2}}$. As an
$\frac{\varepsilon}{2}$-cover $\mathcal{E}_{2}$ for $[0,B]$ (with respect to
$L_{1}$-norm), we can partition $[0,B]$ into $\frac{B}{\varepsilon}$ intervals
of length $\varepsilon$ and consider the end points of those intervals as the
$\frac{\varepsilon}{2}$-cover. Now the cartesian product of these two
$\frac{\varepsilon}{2}$-covers, $\mathcal{E}_{1}\times\mathcal{E}_{2}$, is an
$\varepsilon$-cover of size

$\frac{B}{\varepsilon}\left(\frac{4\beta
d\sqrt{d}}{\varepsilon}\right)^{d^{2}}.$

for $\mathcal{X}=\mathcal{M}\times[0,B]$ with respect to metric ${\rm d}$.

##  Appendix D Uniform convergence of $R_{N}$ to $R$

Although we have proved in Theorem 2 that the true loss
$R({\boldsymbol{M}},\tau)$ is uniquely minimized at
$({\boldsymbol{M}}^{*},\tau^{*})$ , in reality, we do not have the true loss.
Indeed, we only have access to the empirical loss
$R_{N}({\boldsymbol{M}},\tau)$. In this part, broadly speaking, we will show
that $R_{N}({\boldsymbol{M}},\tau)$ is uniformly close to
$R({\boldsymbol{M}},\tau)$ as $N$ gets large, and then, we conclude, instead
of minimizing $R_{N}({\boldsymbol{M}},\tau)$, we can minimize
$R_{N}({\boldsymbol{M}},\tau)$ to approximately find
$({\boldsymbol{M}}^{*},\tau^{*})$.

In the next lemma, we will see that if the two p.s.d. matrices are close via
spectrum norm, then the Mahalanobis norm defined based on these two matrices
are also close.

######  Observation 3 (Equation (5))

For given two p.s.d. ${\boldsymbol{M}}_{1}$ and ${\boldsymbol{M}}_{2}$,

$|\|\boldsymbol{x}\|^{2}_{{\boldsymbol{M}}_{1}}-\|\boldsymbol{x}\|^{2}_{{\boldsymbol{M}}_{2}}|\leq\|{\boldsymbol{M}}_{1}-{\boldsymbol{M}}_{2}\|_{2}\|\boldsymbol{x}\|^{2}.$

Proof  Using Cauchy–Schwarz inequality and the definition of the spectral
norm, we obtain

$\displaystyle|\|\boldsymbol{x}\|^{2}_{{\boldsymbol{M}}_{1}}-\|\boldsymbol{x}\|^{2}_{{\boldsymbol{M}}_{2}}|$
$\displaystyle=|\boldsymbol{x}^{\top}({\boldsymbol{M}}_{1}-{\boldsymbol{M}}_{2})\boldsymbol{x}|$
$\displaystyle=\langle\boldsymbol{x}^{\top},\boldsymbol{x}^{\top}({\boldsymbol{M}}_{1}-{\boldsymbol{M}}_{2})\rangle$
$\displaystyle\leq\|\boldsymbol{x}\|\|({\boldsymbol{M}}_{1}-{\boldsymbol{M}}_{2})\boldsymbol{x}\|$
$\displaystyle\leq\|({\boldsymbol{M}}_{1}-{\boldsymbol{M}}_{2})\|_{2}\|\boldsymbol{x}\|^{2},$

concluding the inequality.

As $-\log\Phi_{\rm Noise}(\cdot)$ is a decreasing function, using Equation 2,
we have

$0\leq-\log\Phi_{\rm
Noise}\left(\ell_{i}(\|\boldsymbol{z}_{i}\|^{2}_{{\boldsymbol{M}}}-\tau)\right)\leq-\log\Phi_{\rm
Noise}(-\beta A)=T,$

which indicates that the random variables $\boldsymbol{z}_{i}$’s are bounded.
Whenever we are dealing with a summation of bounded i.i.d. random variables,
one strong concentration inequality to use is Chernoff-Hoeffding bound. This
inequality states if $X_{1},...,X_{N}$ are $N$ independent random variables
such that ${\displaystyle X_{i}\in[a_{i},b_{i}]}$ almost surely for all $i$,
and $\displaystyle S_{N}=\frac{X_{1}+\cdots+X_{N}}{N}$, then

$\displaystyle\operatorname{P}\left(\left|S_{n}-\mathrm{E}\left[S_{n}\right]\right|\geq\alpha\right)\leq
2\exp\left(-{\frac{2N^{2}\alpha^{2}}{\sum_{i=1}^{N}(b_{i}-a_{i})^{2}}}\right).$

Since,

$\displaystyle R_{N}({\boldsymbol{M}},\tau)$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}-\log\Phi_{\rm
Noise}\left(\ell_{i}(\|\boldsymbol{z}_{i}\|^{2}_{{\boldsymbol{M}}}-\tau)\right)$

and $\mathbb{E}(R_{N}({\boldsymbol{M}},\tau))=R({\boldsymbol{M}},\tau)$, we
can use Chernoff-Hoeffding bound to control the probability that
$|R_{N}({\boldsymbol{M}},\tau)-R({\boldsymbol{M}},\tau)|$ is large.

###### Lemma 17

If $E=\\{({\boldsymbol{M}}_{i},\tau_{i})\colon i=1,\ldots,m=m(\alpha)\\}$ is
an $\alpha$-cover for $\mathcal{M}\times[0,B]$, then

$\displaystyle
P\big{(}|R_{N}({\boldsymbol{M}}_{i},\tau_{i})-R({\boldsymbol{M}}_{i},\tau_{i})|\geq\alpha\text{
for some $i\in[m]$}\big{)}$ $\displaystyle\leq
2me^{-\frac{2N\alpha^{2}}{T^{2}}}.$

Proof  Consider a fixed $i\in[m]$. For simplicity of notation, set
$Z=-\log\sigma\left(\ell(\|\boldsymbol{x}\|_{{\boldsymbol{M}}_{i}}-\tau_{i})\right)$
and
$Z_{j}=-\log\sigma\left(\ell_{j}(\|\boldsymbol{x}_{j}\|_{{\boldsymbol{M}}_{i}}-\tau_{i})\right)$.
As we explained above, $Z\in[0,T]$, see Table 1, and, using Chernoff-Hoeffding
Inequality, we obtain

$\displaystyle
P\big{(}|R_{N}({\boldsymbol{M}}_{i},\tau_{i})-R({\boldsymbol{M}}_{i},\tau_{i})|\geq\alpha\big{)}$
$\displaystyle=P\Big{(}\Big{|}\frac{1}{N}\sum_{j=1}^{N}Z_{j}-E(Z)\Big{|}\geq\alpha\Big{)}$
$\displaystyle\leq 2e^{-\frac{2N\alpha^{2}}{T^{2}}}.$

Now, using union bound, we have the desired inequality.
In the next theorem, we prove that, with high probability, the empirical loss
$R_{N}$ is everywhere close to the true loss $R$.

###### Theorem 18

For any $\varepsilon,\delta>0$, assume parameters $B$, $F$, and $\beta$ are
constants, define

$N_{d}(\varepsilon,\delta)=O\left(\frac{1}{\varepsilon^{2}}\Big{[}\log\frac{1}{\delta}+d^{2}\log\frac{d}{\varepsilon}\Big{]}\right).$

If $N>N_{d}(\varepsilon,\delta)$, then with probability at least $1-\delta$,

$\sup_{({\boldsymbol{M}},\tau)\in\mathcal{M}\times[0,B]}\left|R_{N}({\boldsymbol{M}},\tau)-R({\boldsymbol{M}},\tau)\right|<\varepsilon.$

Proof  To prove the assertion, we can equivalently prove

$P\left(\sup_{({\boldsymbol{M}},\tau)\in\mathcal{M}\times[0,B]}\left|R_{N}({\boldsymbol{M}},\tau)-R({\boldsymbol{M}},\tau)\right|\geq\varepsilon\right)<\delta.$

Set $\alpha=\frac{\varepsilon}{3\zeta(F+1)}$. Consider
$\mathcal{E}=\\{({\boldsymbol{M}}_{i},\tau_{i});i=1,\ldots,m=m(\alpha)\\}$ as
an $\alpha$-cover for $\mathcal{M}\times[0,B].$ For an arbitrary
$({\boldsymbol{M}},\tau)$, there is an index $i\in[m]$ such that

$d\left(({\boldsymbol{M}},\tau)-({\boldsymbol{M}}_{i},\tau_{i})\right)<\alpha.$

Consequently, using Lemma 3,

$|R({\boldsymbol{M}},\tau)-R({\boldsymbol{M}}_{i},\tau_{i})|<(F+1)\zeta\alpha=\frac{\varepsilon}{3}$

and

$\displaystyle|R_{N}({\boldsymbol{M}},\tau)-R_{N}({\boldsymbol{M}}_{i},\tau_{i})|<(F+1)\zeta\alpha=\frac{\varepsilon}{3}.$

So far, we have proved that for every $({\boldsymbol{M}},b)$, there exists an
index $i\in[m]$ such that

$|R({\boldsymbol{M}},\tau)-R({\boldsymbol{M}}_{i},\tau_{i})|<\frac{\varepsilon}{3}\quad\quad\text{and}\quad\quad|R_{N}({\boldsymbol{M}},\tau)-R_{N}({\boldsymbol{M}}_{i},\tau_{i})|<\frac{\varepsilon}{3}.$

Using triangle inequality, it concludes

$\displaystyle\left|R_{N}({\boldsymbol{M}},\tau)-R({\boldsymbol{M}},\tau)\right|\leq$
$\displaystyle\left|R_{N}({\boldsymbol{M}},\tau)-R_{N}({\boldsymbol{M}}_{i},\tau_{i})\right|+\left|R_{N}({\boldsymbol{M}}_{i},\tau_{i})-R({\boldsymbol{M}}_{i},\tau_{i})\right|$
$\displaystyle+\left|R({\boldsymbol{M}}_{i},\tau_{i})-R({\boldsymbol{M}},\tau)\right|$
$\displaystyle\leq$
$\displaystyle\frac{2\varepsilon}{3}+\left|R_{N}({\boldsymbol{M}}_{i},\tau_{i})-R({\boldsymbol{M}}_{i},\tau_{i})\right|.$

Via Lemma 16, there is an $\alpha$-cover of size

$\displaystyle m(\alpha)$ $\displaystyle=\frac{B}{\alpha}\left(\frac{4\beta
d\sqrt{d}}{\alpha}\right)^{d^{2}}$
$\displaystyle=\frac{3\zeta(F+1)B}{\varepsilon}\left(\frac{12\zeta(F+1)\beta
d\sqrt{d}}{\varepsilon}\right)^{d^{2}}$

for $\mathcal{X}=\mathcal{M}\times[0,B]$ with respect to metric ${\rm d}$. On
the other hand,

$\displaystyle\frac{T^{2}}{2\alpha^{2}}\log\frac{2m}{\delta}$
$\displaystyle=\frac{9\zeta(F+1)^{2}T^{2}}{2\varepsilon^{2}}\Big{[}\log\frac{6\zeta(1+F)B}{\varepsilon\delta}+d^{2}\log\frac{12\zeta(1+F)\beta}{\varepsilon}+\frac{3}{2}d^{2}\log
d\Big{]}$ (12)
$\displaystyle=O\left(\frac{1}{\varepsilon^{2}}\Big{[}\log\frac{1}{\delta}+d^{2}\log\frac{d}{\varepsilon}\Big{]}\right)$

As setting

$\displaystyle N$
$\displaystyle>\frac{T^{2}}{2\alpha^{2}}\log\frac{2m}{\delta}$
$\displaystyle=O\left(\frac{1}{\varepsilon^{2}}\Big{[}\log\frac{1}{\delta}+d^{2}\log\frac{d}{\varepsilon}\Big{]}\right)$

implies $2me^{-\frac{2N\alpha^{2}}{T^{2}}}<\delta$, using Lemma 17, we obtain,
with probability at least $1-\delta$,

$\left|R_{N}({\boldsymbol{M}}_{j},\tau_{j})-R({\boldsymbol{M}}_{j},\tau_{j})\right|<\alpha=\frac{\varepsilon}{3\zeta(F+1)}<\frac{\varepsilon}{3}\quad\quad\text{for
all $j\in[m]$}.$

Therefore, if
$N>O\left(\frac{1}{\varepsilon^{2}}\Big{[}\log\frac{1}{\delta}+d^{2}\log\frac{d}{\varepsilon}\Big{]}\right)$,
with probability at least $1-\delta$, for all $({\boldsymbol{M}},\tau)$ we
have

$\left|R_{N}({\boldsymbol{M}},\tau)-R({\boldsymbol{M}},\tau)\right|\leq\varepsilon,$

as desired.

##  Appendix E Simple Noise Properties in Table 1

In Subsections 3.1, 3.2, 3.3, 3.4, we derived some properties of $\Phi_{\rm
Noise}(\cdot)$ when noise is one of the simple noises listed in Table 1. This
section can be seen as a complementary section for those sections. For noises
listed in Table 1, one can verify that each of those noise distributions is
simple. Here, we verify some other information listed in that table. When
${\rm Noise}(\eta)$ is simple, $-\log\Phi_{\rm Noise}(\eta)$ is a decreasing
convex function which implies that ${{\rm d}\over{\rm
d}\eta}\left(-\log\Phi_{\rm Noise}(\eta)\right)=-\frac{{\rm
Noise}(\eta)}{\Phi_{\rm Noise}(\eta)}$ is a negative increasing function.
Therefore, $-\log\Phi_{\rm Noise}(\eta)$ is $\zeta$-Lipschitz over $[-\beta
F,\beta B]$ for

$\zeta=\frac{{\rm Noise}(-\beta F)}{\Phi_{\rm Noise}(-\beta F)}.$

Also, from Equation 6, we know

$T=-\log\Phi_{{\rm Normal}}(-\beta F).$

In what follows, we approximate $\zeta$ and $T$ for each noise choice.

  * •

Logistic noise. In this case, it is easy to see that

$\zeta=\frac{\sigma(-\beta F)\sigma(\beta F)}{\sigma(-\beta F)}=\sigma(\beta
F)<1$

and $T=-\log\sigma(-\beta F)=\log(1+e^{\beta F})\leq 1+\beta F$. Thus,
$\Phi_{\rm Logistic}(\eta)$ in $1$-log-Lipschitz and $T=O(\beta F).$

  * •

Normal noise. We start with an approximation of $\Phi_{\rm Normal}(\eta)$ (the
lower bound comes from a webpage Cook ). We set $\Phi^{c}_{\rm
Normal}(\eta)=1-\Phi_{\rm Normal}(\eta)=\Phi_{\rm Normal}(-\eta)$.

Set $g(t)=\Phi^{c}(t)-\frac{1}{\sqrt{2\pi}}\frac{t}{t^{2}+1}e^{-t^{2}/2}$. As
$g(0)>0$,
$g^{\prime}(t)=-\frac{2}{\sqrt{2\pi}}\frac{e^{-t^{2}/2}}{(t^{2}+1)^{2}}<0$ and
$\lim\limits_{t\rightarrow+\infty}g(t)=0$, we obtain
$\Phi^{c}(t)\geq\frac{1}{\sqrt{2\pi}}\frac{t}{t^{2}+1}e^{-t^{2}/2}$. Using the
above formula for $\zeta$, we have

$\displaystyle\zeta$ $\displaystyle={1\over\sqrt{2\pi}}\frac{e^{-\frac{(\beta
F)^{2}}{2}}}{\Phi(-\beta F)}$
$\displaystyle={1\over\sqrt{2\pi}}\frac{e^{-\frac{(\beta
F)^{2}}{2}}}{\Phi^{c}(\beta F)}$ $\displaystyle\leq\frac{(\beta
F)^{2}+1}{\beta F}=O(\beta F).$

Also,

$T=-\log\Phi_{{\rm Normal}}(-\beta F)=-\log\Phi^{c}_{{\rm Normal}}(\beta
F)=O((\beta F)^{2}).$

  * •

Laplace noise. In this setting,

$\zeta=\frac{{\rm Noise}(-\beta F)}{\Phi_{\rm Noise}(-\beta
F)}=\frac{\frac{1}{2}e^{-\beta F}}{\frac{1}{2}e^{-\beta F}}=1$

and $T=-\log\Phi_{{\rm Normal}}(-\beta F)=\beta F\log 2=O(\beta F)$.

  * •

Hyperbolic secant noise.

Remind that $\Phi^{c}(t)=\int_{t}^{\infty}\frac{1}{2}{\rm
sech}(\frac{\pi}{2}\eta){\rm d}\eta$. Set

$\displaystyle g(t)$ $\displaystyle=\Phi^{c}(t)-\frac{1}{\pi}{\rm
sech}(\frac{\pi}{2}t)$

Also, note that $g(0)>0$, $\lim\limits_{t\rightarrow+\infty}g(t)=0$, and

$\displaystyle g^{\prime}(t)$ $\displaystyle=-\frac{1}{2}{\rm
sech}(\frac{\pi}{2}t)+\frac{1}{2}\tanh(\frac{\pi}{2}t)\operatorname{sech}(\frac{\pi}{2}t)$
$\displaystyle=\frac{1}{2}{\rm
sech}(\frac{\pi}{2}t)\left(-1+\tanh(\frac{\pi}{2}t)\right)<0\quad\quad\forall
t\in\mathbb{R}.$

This implies that, for each $t\in\mathbb{R}$,

$\Phi^{c}(t)>\frac{1}{\pi}{\rm sech}(\frac{\pi}{2}t).$

Using this inequality, we obtain

$\displaystyle\zeta$ $\displaystyle=\frac{{\rm Noise}(-\beta F)}{\Phi_{\rm
Noise}(-\beta F)}$ $\displaystyle=\frac{\frac{1}{2}{\rm
sech}(-\frac{\pi}{2}\beta F)}{\Phi_{\rm Noise}^{c}(\beta F)}$
$\displaystyle\leq\frac{\pi}{2}.$

Furthermore,

$\displaystyle T$ $\displaystyle=-\log\Phi_{{\rm Normal}}(-\beta
F)\leq-\log\frac{1}{\pi}{\rm sech}(\frac{\pi}{2}\beta F)$
$\displaystyle=\log\pi+\log\cosh(\frac{\pi}{2}\beta F)=O(\beta F).$

##  Appendix F Connection Between $L_{1}(f)$ Norm and Spectral Norm

In Theorem 7 and Corollary 8, we study the connection between $L_{f}(f)$ norm
and the sample complexity of our problem. However, as the $L_{1}(f)$-metric is
dependent on the distribution $f(\boldsymbol{z})$, which is unavoidable, it is
not very intuitive. We indeed prefer some more informative norms such as
spectral norm. However, to this end, we must restrict the distribution
$f(\boldsymbol{z})$. We here assume that

  * •

We here assume that there is a positive constant $c$, such that
$f(\boldsymbol{z})\geq c$ for each $\boldsymbol{z}\in B^{d}(1);$ recall we
assume almost surely $\|\boldsymbol{z}\|^{2}\leq F=1$ in this section, and
$B^{d}(1)=\\{\boldsymbol{z}\in\mathbb{R}^{d}\mid\|\boldsymbol{z}\|\leq 1\\}$.

We next prove a statement similar to Corollary 8 but in terms of the ${\rm
d}$-metric instead of the $L_{1}(f)$-metric. The following definition is also
needed for the following two results. For $0\leq a<b\leq 1$, where $z_{1}$ is
the first coordinate of $\boldsymbol{z}$, define

${\rm Cone}(a,b)=\left\\{\boldsymbol{z}=(z_{1},\ldots,z_{d})\mid
3z_{1}^{2}\geq 2\|\boldsymbol{z}\|_{2}^{2}\ \&\ a\leq z_{1}\leq b\right\\}.$
(13)

###### Lemma 19

If $f(\boldsymbol{z})\geq c>0$ for each $\boldsymbol{z}\in B^{d}(1)$, then for
all $({\boldsymbol{M}},\tau)\in\mathcal{M}\times[0,B]$

$\|({\boldsymbol{M}},\tau)-({\boldsymbol{M}}^{*},\tau^{*})\|_{L_{1}(f)}\geq\frac{c\pi^{d/2}}{20\Gamma(d/2+1)}(\frac{1}{18})^{d}{\rm
d}\left(({\boldsymbol{M}},\tau),({\boldsymbol{M}}^{*},\tau^{*})\right).$

In particular, if $f(\boldsymbol{z})$ is uniform on unit disk, then for all
$({\boldsymbol{M}},\tau)\in\mathcal{M}\times[0,B]$

$\|({\boldsymbol{M}},\tau)-({\boldsymbol{M}}^{*},\tau^{*})\|_{L_{1}(f)}\geq\frac{1}{20}(\frac{1}{18})^{d}{\rm
d}\left(({\boldsymbol{M}},\tau),({\boldsymbol{M}}^{*},\tau^{*})\right).$

Proof  We remind that

$\displaystyle\|({\boldsymbol{M}},\tau)-({\boldsymbol{M}}^{*},\tau^{*})\|_{L_{1}(f)}$
$\displaystyle=\int
f(\boldsymbol{z})\left|(\|\boldsymbol{z}\|^{2}_{{\boldsymbol{M}}}-\tau)-(\|\boldsymbol{z}\|^{2}_{{\boldsymbol{M}}^{*}}-\tau^{*})\right|{\rm
d}\boldsymbol{z}$ $\displaystyle=\int
f(\boldsymbol{z})\Big{|}\boldsymbol{z}^{\top}(\underbrace{\hat{{\boldsymbol{M}}}-{\boldsymbol{M}}^{*}}_{=\bar{{\boldsymbol{M}}}})\boldsymbol{z}-(\underbrace{\tau-\tau^{*}}_{=\bar{\tau}})\Big{|}{\rm
d}\boldsymbol{z}$ $\displaystyle=\int
f(\boldsymbol{z})\left|\boldsymbol{z}^{\top}\bar{{\boldsymbol{M}}}\boldsymbol{z}-\bar{\tau}\right|{\rm
d}\boldsymbol{z}.$

Note that $\bar{{\boldsymbol{M}}}$ is a symmetric matrix. So, there are a real
value orthonormal matrix ${\boldsymbol{U}}$ and a real value diagonal matrix
${\boldsymbol{D}}$ such that
$\bar{{\boldsymbol{M}}}={\boldsymbol{U}}^{\top}{\boldsymbol{D}}{\boldsymbol{U}}$.
Let the vector ${\boldsymbol{\lambda}}$ denote the diagonal of
${\boldsymbol{D}}$. Without loss of generality, we assume that
$\lambda_{1}=\max\\{|\lambda_{1}|,\ldots,|\lambda_{d}|\\}.$ One can verify
that $\lambda_{1}=\|\bar{{\boldsymbol{M}}}\|_{2}.$ As ${\boldsymbol{U}}$ is
orthonormal, we have

$\displaystyle\int
f(\boldsymbol{z})\left|\boldsymbol{z}^{\top}\bar{{\boldsymbol{M}}}\boldsymbol{z}-\bar{\tau}\right|{\rm
d}\boldsymbol{z}$ $\displaystyle=\int
f(\boldsymbol{z})\left|({\boldsymbol{U}}\boldsymbol{z})^{\top}{\boldsymbol{D}}({\boldsymbol{U}}\boldsymbol{z})-\bar{\tau}\right|{\rm
d}\boldsymbol{z}$ $\displaystyle=\int
f({\boldsymbol{U}}^{\top}\boldsymbol{z})\left|\boldsymbol{z}^{\top}{\boldsymbol{D}}\boldsymbol{z}-\bar{\tau}\right|{\rm
d}\boldsymbol{z}$ $\displaystyle=\int
f({\boldsymbol{U}}^{\top}\boldsymbol{z})\left|\sum_{i=1}^{d}z_{i}^{2}\lambda_{i}-\bar{\tau}\right|{\rm
d}\boldsymbol{z}$

Note that
$\sum_{i=1}^{d}z_{i}^{2}\lambda_{i}\leq\lambda_{1}\|\boldsymbol{z}\|^{2}$ and,
for each $\boldsymbol{z}\in{\rm Cone}(0,1)$,

$\displaystyle\sum_{i=1}^{d}z_{i}^{2}\lambda_{i}$
$\displaystyle\geq\lambda_{1}z_{1}^{2}-\sum_{i=2}^{d}\lambda_{1}z_{i}^{2}$
$\displaystyle=\lambda_{1}\left(2z_{1}^{2}-\sum_{i=1}^{d}z_{i}^{2}\right)$
$\displaystyle=\lambda_{1}(2z_{1}^{2}-\|\boldsymbol{z}\|_{2}^{2})\geq\frac{\lambda_{1}}{2}z_{1}^{2}$

Set $q=\lambda_{1}+|\bar{\tau}|={\rm
d}\left((\bar{{\boldsymbol{M}}},\bar{\tau}),({\boldsymbol{M}}^{*},\tau^{*})\right)$.
We next consider two cases based on $|\bar{\tau}|$ and $q$.

  * •

$|\bar{\tau}|\leq 0.1q$. It implies $\lambda_{1}\geq 0.9q$ and thus, if
$\boldsymbol{z}\in{\rm Cone}(\frac{1}{\sqrt{3}},1)$, then

$\displaystyle\sum_{i=1}^{d}z_{i}^{2}\lambda_{i}-\bar{\tau}$
$\displaystyle\geq\frac{\lambda_{1}}{2}z_{1}^{2}-\bar{\tau}$
$\displaystyle\geq q\left(\frac{9}{20}z_{1}^{2}-\frac{1}{10}\right)$
$\displaystyle\geq\frac{1}{20}q.$

Therefore,

$\displaystyle\int
f({\boldsymbol{U}}^{\top}\boldsymbol{z})\Big{|}\sum_{i=1}^{d}z_{i}^{2}\lambda_{i}-\bar{\tau}\Big{|}{\rm
d}\boldsymbol{z}$ $\displaystyle\geq\frac{q}{20}\int_{\boldsymbol{z}\in{\rm
Cone}(\frac{1}{\sqrt{3}},1)}f({\boldsymbol{U}}^{\top}\boldsymbol{z}){\rm
d}\boldsymbol{z}$
$\displaystyle=\frac{q}{20}\mu_{f}({\boldsymbol{U}}^{\top}{\rm
Cone}(\frac{1}{\sqrt{3}},1))$ (14)
$\displaystyle\geq\frac{cq}{20}\int_{\boldsymbol{z}\in{\rm
Cone}(\frac{1}{\sqrt{3}},1)\cap B^{d}(1)}{\rm d}\boldsymbol{z}$
$\displaystyle=\frac{cq}{20}\times{\rm Volume}\left({\rm
Cone}(\frac{1}{\sqrt{3}},1)\cap B^{d}(1)\right)$
$\displaystyle\geq\frac{cq}{20}\times{\rm Volume}\left({\rm
Cone}(\frac{1}{\sqrt{3}},\sqrt{\frac{2}{3}})\right)$
$\displaystyle=\frac{cq}{20d}\left[\sqrt{\frac{2}{3}}V^{d-1}(\sqrt{\frac{1}{3}}))-\sqrt{\frac{1}{3}}V^{d-1}(\sqrt{\frac{1}{6}})\right]$
$\displaystyle=\frac{cq}{20}\frac{\pi^{\frac{d-1}{2}}}{d3^{\frac{d}{2}}\Gamma(\frac{d+1}{2})}\left[\sqrt{2}-\frac{1}{2^{\frac{d-1}{2}}}\right].$
(15)

  * •

$|\bar{\tau}|>0.1q$. It implies $\lambda_{1}<0.9q$ and thus

$\displaystyle\Big{|}\bar{\tau}-\sum_{i=1}^{d}z_{i}^{2}\lambda_{i}\Big{|}$
$\displaystyle\geq|\bar{\tau}|-\lambda_{1}\|\boldsymbol{z}\|_{2}^{2}$
$\displaystyle\geq
q\left(\frac{1}{10}-\frac{9}{10}\|\boldsymbol{z}\|_{2}^{2}\right)$
$\displaystyle\geq\frac{q}{20}\quad\text{ if
$\|\boldsymbol{z}\|_{2}^{2}\leq\frac{1}{18}$}$

Therefore,

$\displaystyle\int
f({\boldsymbol{U}}^{\top}\boldsymbol{z})\Big{|}\sum_{i=1}^{d}z_{i}^{2}\lambda_{i}-\bar{\tau}\Big{|}{\rm
d}\boldsymbol{z}$ $\displaystyle\geq\frac{q}{20}\int_{\boldsymbol{z}\in
B^{q}(\frac{1}{18})}f({\boldsymbol{U}}^{\top}\boldsymbol{z}){\rm
d}\boldsymbol{y}$ $\displaystyle=\frac{q}{20}\mu_{f}(B^{q}(\frac{1}{18}))$
(16) $\displaystyle\geq\frac{cq}{20}\times{\rm
Volume}\left(B^{q}(\frac{1}{18})\right)$
$\displaystyle=\frac{cq}{20}\frac{\pi^{d/2}}{\Gamma(d/2+1)}(\frac{1}{18})^{d}$
(17)

Now Lower bounds (15) and (17) implies the proof.
If $f(\boldsymbol{z})$ is a rotationally symmetric pdf, then
$\mu_{f}({\boldsymbol{U}}^{\top}B)=\mu_{f}(B)$ for any measurable set $B$.
Therefore, combining two Lower bounds (14) and (16), we obtain the following
lemma

###### Lemma 20

If $f(\boldsymbol{z})$ is a rotationally symmetric pdf, then for all
$({\boldsymbol{M}},\tau)\in\mathcal{M}\times[0,B]$

${\|({\boldsymbol{M}},\tau)-({\boldsymbol{M}}^{*},\tau^{*})\|_{L_{1}(f)}\over{\rm
d}\left(({\boldsymbol{M}},\tau),({\boldsymbol{M}}^{*},\tau^{*})\right)}\geq\frac{1}{20}\max\left(\mu_{f}({\rm
Cone}(\frac{1}{\sqrt{3}},1),\mu_{f}(B^{d}(\frac{1}{18}))\right).$

##  Appendix G Further Experimental Study

This section can be seen as a complementary section for Section 5.

##### Loss Function Behavior.

In Subsection 5.2, we experimentally studied the logistic model with different
noises. In Figures 2 and 3, we evaluate the model for different noise in terms
of eigenvalue recover and the accuracies. As a complementary information to
these figures, in Figure 9, as the iteration increases, we compare the values
of the loss function $R_{N}$ on $(\hat{{\boldsymbol{M}}},\hat{\tau})$,
compared with $({\boldsymbol{M}}^{*}/s,\tau^{*}/s)$ which we do not expect to
surpass. Observe that the loss on $({\boldsymbol{M}}^{*}/s,\tau^{*}/s)$ is a
constant red line at the bottom at around $0.23$. When considering Logistic
noise (blue) we reach this loss around 700 iterations, and nearly do when
considering Gaussian noise. For other types of noise, the method does worse;
Noisy labeling only achieves a loss value around $0.5$.

Figure 9: Loss history for Logistic model and different noises

##### Larger Dimension.

In Section 5, we dealt with small value for dimension $d$ ($d=10$ for
synthetic and $d=9,24$ for real data). In this section, following the same
approach as in Subsection 5.1, we generate synthetic data with $d=100$ and
${\rm rank}({\boldsymbol{M}}^{*})=30$. We also set the level of noise at
$20\%$. In Figure 10, we observe that how the sample complexity can be
affected by the dimension. When the sample size is less than $30K$, the model
overfits, which is completely natural as our model has $d^{2}+1$ parameters.
But for the larger values, we can see that the model starts to neutralize the
noise and the non-noisy accuracy approaches to $1$ (blue and magenta curves).
We have $97.59\%,97.56\%$ non-noisy train and test accuracy for $600K$ sample
size.

Figure 10: Sample complexity for $d=100$


# Scalable Non-Cartesian Magnetic Resonance Imaging with R2D2

Yiwei Chen1, Chao Tang1,2, Amir Aghabiglou1, Chung San Chu1, Yves Wiaux12 The
work was supported by EPSRC under grants EP/T028270/1 and ST/W000970/1.
Computing resources came from the Cirrus UK National Tier-2 HPC Service at
EPCC (http://www.cirrus.ac.uk) funded by the University of Edinburgh and EPSRC
(EP/P020267/1). 1Institute of Sensors, Signals and Systems, Heriot-Watt
University, Edinburgh EH14 4AS, United Kingdom
2EPCC, University of Edinburgh, Edinburgh EH8 9BT, United Kingdom  Email:
<EMAIL_ADDRESS>
###### Abstract

We propose a new approach for non-Cartesian magnetic resonance image
reconstruction. While unrolled architectures provide robustness via data-
consistency layers, embedding measurement operators in Deep Neural Network
(DNN) can become impractical at large scale. Alternative Plug-and-Play (PnP)
approaches, where the denoising DNNs are blind to the measurement setting, are
not affected by this limitation and have also proven effective, but their
highly iterative nature also affects scalability. To address this scalability
challenge, we leverage the “Residual-to-Residual DNN series for high-Dynamic
range imaging (R2D2)” approach recently introduced in astronomical imaging.
R2D2’s reconstruction is formed as a series of residual images, iteratively
estimated as outputs of DNNs taking the previous iteration’s image estimate
and associated data residual as inputs. The method can be interpreted as a
learned version of the Matching Pursuit algorithm. We demonstrate R2D2 in
simulation, considering radial k-space sampling acquisition sequences. Our
preliminary results suggest that R2D2 achieves: (i) suboptimal performance
compared to its unrolled incarnation R2D2-Net, which is however non-scalable
due to the necessary embedding of NUFFT-based data-consistency layers; (ii)
superior reconstruction quality to a scalable version of R2D2-Net embedding an
FFT-based approximation for data consistency; (iii) superior reconstruction
quality to PnP, while only requiring few iterations.

##  I Introduction

Magnetic Resonance Imaging (MRI) enables high-precision reconstruction of
structures and organs within the human body via the interaction of magnetic
fields and radio waves. In MRI, signals are measured in the spatial frequency
domain, namely k-space, using multiple receiver coils. The complex-valued MR
image is reconstructed based on the inverse Fourier transform. For this proof
of concept and without loss of generality, similarly to [1], we consider real-
valued positive images, which would correspond to the magnitude of complex-
valued MR images. We also restrict the study to a single-coil setting. We
consider non-Cartesian k-space sampling, more precisely radial sampling
sequences, which are commonly used and possess advantages such as reduced
motion artifacts in dynamic scenes [2]. The acquisition process can thus be
formulated as

$\bm{y}=\bm{\mathsf{\Phi}}\bm{\bar{x}}+\bm{n},$ (1)

where $\bm{y}\in\mathbb{C}^{M}$ denotes the k-space measurements,
$\bm{\mathsf{\Phi}}\colon\mathbb{R}^{N}\to\mathbb{C}^{M}$ is a non-uniform
Fourier sampling measurement operator, $\bm{\bar{x}}\in\mathbb{R}^{N}_{+}$ is
the Ground Truth (GT) image and $\bm{n}\in\mathbb{C}^{M}$ represents a
complex-valued Gaussian random noise following $\mathcal{N}(0,\tau^{2})$.
$\bm{\mathsf{\Phi}}$ can be implemented via the Non-Uniform Fast Fourier
Transform (NUFFT), _i.e._ as the product of $\bm{\mathsf{U,F,Z}}$, where
$\bm{\mathsf{Z}}$ is a zero-padding operator, $\bm{\mathsf{F}}$ denotes the
Fast Fourier Transform (FFT), and $\bm{\mathsf{U}}$ is an interpolation
operator.

To accelerate the acquisition process, reducing the number of k-space
measurements is common. DNN models like PnP algorithms [3] and unrolled
networks [4] have gained attention for their promising performance on
undersampled MRI data. PnP algorithms, bridging optimization theory and deep
learning, train denoising DNNs deployed within optimization algorithms to
replace handcrafted regularization operators. While effective and unaffected
by measurement settings, they face scalability challenges due to their highly
iterative nature. Unrolled DNNs ensure the consistency between the
reconstructions and measurements by mirroring optimization algorithm
iterations across layers, achieving high imaging precision in non-Cartesian
MRI. However, in large-dimensional settings, typically when a high number of
acquisition coils are involved or in 3D or 4D dynamic MRI [5], embedding NUFFT
operators into DNN architectures entails large computational cost and memory
requirements at both training and inference stages. These can rapidly become
prohibitive due to corresponding limitation of GPU hardware [4].

Recently, the R2D2 approach, a learned variant of the Matching Pursuit
approach [6], has been demonstrated to deliver a new regime of joint image
reconstruction quality and speed in radio astronomical imaging [7, 8, 9].
R2D2’s reconstruction is formed as a series of residual images, iteratively
estimated as outputs of DNNs taking the previous iteration’s image estimate
and associated back-projected data residual as inputs. In contrast with an
unrolled DNN approach, R2D2 thus utilizes a series of DNNs without embedding
the measurement operator into their individual architectures. Interestingly,
the number of required networks (iterations) was shown to be extremely small
compared to the typical number of iterations of a PnP algorithm.

We propose two unrolled variants of the R2D2 model, named R2D2-Net (FFT) and
R2D2-Net (NUFFT), embedding distinct data-consistency layers, respectively
based on the NUFFT and an FFT approximation. Preliminary results suggest that
R2D2 naturally exhibits suboptimal performance compared to the R2D2-Net
(NUFFT), which, however, is non-scalable due to the essential embedding of
NUFFT-based data-consistency layers. However, R2D2 largely outperforms the
scalable R2D2-Net (FFT) variant. With respect to the state of the art, R2D2
also largely outperforms the (non-scalable) unrolled network named NC-PDNet,
also embedding the NUFFT [4], and a (non-scalable) PnP benchmark [10].

##  II Methodology

###  II-A R2D2 approach

Given an image estimate $\bm{x}$, the back-projected data residual denoted by
$\bm{r}$, also named the data-consistency term [11], is given by

$\bm{r}=\bm{x}_{\text{d}}-\kappa\operatorname{Re}\\{\bm{\mathsf{\Phi}^{\dagger}}\bm{\mathsf{\Phi}}\bm{x}\\},$
(2)

where
$\bm{x}_{\text{d}}=\kappa\operatorname{Re}\\{\bm{\mathsf{\Phi}^{\dagger}}\bm{y}\\}$
is the data back-projected from the raw measurements. Here,
$\kappa=1/\max(\operatorname{Re}\\{\bm{\mathsf{\Phi}^{\dagger}}\bm{\mathsf{\Phi}}\bm{\delta}\\})$
is the normalization parameter with $\bm{\delta}$ being the Dirac image with
the central pixel value as 1 and others as 0.

The R2D2 approach aims to utilize the previous image estimate $\bm{x}^{i-1}$
and corresponding $\bm{r}^{i-1}$ to reduce the discrepancies between the image
estimate and GT image by a series of DNNs, denoted as
$\\{\bm{\mathsf{G}}_{\bm{\theta_{i}}}\\}_{1\leq i\leq I}$, with learnable
parameters denoted as $\\{\bm{\theta_{i}}\\}_{1\leq i\leq I}$. The iteration
structure reads

$\bm{x}^{i}=\bm{x}^{i-1}+\bm{\mathsf{G}}_{\bm{\theta_{i}}}(\bm{r}^{i-1},\bm{x}^{i-1}),$
(3)

with the initialized image estimate and back-projected data residual as
$\bm{x}^{0}=\bm{0}$ and $\bm{r}^{0}=\bm{x}_{\text{d}}$, respectively. DNNs are
trained sequentially. Specifically, for the $i$-th DNN, the loss function is
defined as:

$\underset{\bm{\theta}_{i}}{\min}\frac{1}{K}\sum_{k=1}^{K}\lVert\bm{\bar{x}}_{k}-[\bm{x}^{i-1}_{k}+\bm{\mathsf{G}}_{\bm{\theta_{i}}}(\bm{r}^{i-1}_{k},\bm{x}^{i-1}_{k})]_{+}\rVert_{1},$
(4)

where $\lVert\cdot\rVert_{1}$ is the L1 norm, $K$ is the number of training
samples, and $[\cdot]_{+}$ denotes the projection onto the positive orthant
$\mathbb{R}^{N}_{+}$.

At the inference stage, the last output of the DNN series denoted as
$\bm{x}^{I}$ gives the final reconstruction. It is noteworthy that the first
iteration of R2D2 is equivalent to a standard end-to-end learning approach.

###  II-B DNN structure

The standard U-Net [12] serves as the foundational structure for DNNs,
featuring both a contracting and expanding path. The contracting path
integrates multiple convolutional and pooling layers to gradually reduce the
spatial dimensions of the input image. Conversely, the expanding path
incorporates convolutional and upsampling layers to progressively upsample the
feature maps, resulting in an output with identical dimensions to the input
image. Skip connections are added to link layers between the contracting and
expanding paths.

###  II-C R2D2-Net

Formally, an unrolled variant of R2D2, named R2D2-Net [8], can be developed by
unrolling the R2D2 approach itself, with a predetermined number of internal
iterations. Two implementations, distinguished by data-consistency
computation, are considered. The first, denoted by R2D2-Net (NUFFT), uses the
exact measurement operator based on the NUFFT to calculate the data
consistency. To promote scalability, we employ an FFT approximation of the
exact measurement operator to estimate the data consistency. Firstly, the
point spread function image is calculated as
$\bm{\mathbf{h}}=\kappa\operatorname{Re}\\{\bm{\mathsf{\Phi}^{\dagger}}\bm{\mathsf{\Phi}}\bm{\delta}\\}\in\mathbb{R}^{N}$.
Secondly, we replace
$\kappa\operatorname{Re}\\{\bm{\mathsf{\Phi}^{\dagger}}\bm{\mathsf{\Phi}}\\}$
by
$\operatorname{Re}\\{\bm{\mathsf{F}^{\dagger}}(\bm{\mathsf{F}}\bm{\mathbf{h}})\bm{\mathsf{F}}\\}$,
where the FFT operator and its adjoint are denoted by $\bm{\mathsf{F}}$ and
$\bm{\mathsf{F}^{\dagger}}$, respectively. Note that
$\bm{\mathsf{F}}\bm{\mathbf{h}}$ is precomputed once and stored. Although this
mapping to the back-projected data space is approximate, it is fast and
memory-efficient. This substitution results in a scalable unrolled DNN denoted
by R2D2-Net (FFT).

###  II-D Normalization procedures

Normalization procedures are employed to avoid generalizability issues arising
from large variations in pixel value ranges and stabilize the training
process. An iteration-specific normalization factor for the $i$-th network is
denoted as $\alpha^{i-1}$ and given by the mean pixel value of the previous
estimate $\bm{x}^{i-1}$. Note that $\alpha^{0}$ is the mean pixel value of
$\bm{x}_{\text{d}}$. In training, $\bm{\bar{x}}$, $\bm{x}^{i-1}$ and
$\bm{r}^{i-1}$ are divided by $\alpha^{i-1}$. At the inference stage, the
normalization mapping:
$\bm{\mathsf{G}}_{\bm{\theta_{i}}}(\cdot)\mapsto\alpha^{i-1}\bm{\mathsf{G}}_{\bm{\theta_{i}}}(\cdot/\alpha^{i-1})$
is applied to each subnetwork to make the input normalized and the output
denormalized accordingly.

###  II-E Related works

Firstly, compared to PnP algorithms, DNNs in R2D2 aim to learn high-level
features rather than being trained solely as denoisers [3]. Secondly, in
contrast to unrolled DNNs [11, 4], R2D2 externalizes the data consistency
calculation from the network structure. This strategy relieves the significant
computational cost of NUFFT in training, leading to better scalability.
Formally, R2D2-Net shares the same network architecture as NC-PDNet [4], with
the key distinction being that the core subnetwork architecture of R2D2-Net is
a U-Net instead of a shallow convolution neural network. Furthermore, R2D2-Net
introduces a novel variant employing the FFT approximation to enhance
scalability.

##  III Data Simulation

###  III-A Denoised GT images

We adopt the magnitude MR images with a size of $320\times 320$ from the
FastMRI single-coil knee dataset [13]. As depicted in Fig. 1 (a), the raw GT
image exhibits noticeable noise. To improve the image quality and generate
k-space measurements with diverse noise levels, the raw images undergo
preprocessing by the denoising network SCUNet [14], followed by soft-
thresholding to remove residual backgrounds. The denoised image is then
normalized with intensities ranging from 0 to 1. We split denoised images into
training and validation datasets consisting of 25743 and 6605 images,
respectively.

###  III-B Radial sampling

Radial sampling starts at the centre of k-space and extends outward in an
angular fashion. The radial sampling trajectory in k-space is formulated as
$k_{x}(n)=r\cos\left(\alpha_{n}\right)$,
$k_{y}(n)=r\sin\left(\alpha_{n}\right)$, where $k_{x}(n)$ and $k_{y}(n)$ are
the Cartesian coordinates of the sampled point at the $n$-th spoke and $r$
represents the radial distance from the k-space origin with
$r\in\mathbb{Z},-R\leq r\leq R$. Generally, for a square image, $R$ is the
unilateral dimension of the image [4]. $\alpha_{n}$ is the angle of the $n$-th
spoke with respect to the $k_{x}$ axis. The golden angle,
$\alpha_{g}=111.25^{\circ}$ for the 2D case [2], is a common choice for the
specific angular increment to achieve an efficient and relatively uniform
angular distribution for the radial sampling trajectory, resulting in
$\alpha_{n}=n\alpha_{g}$.

Figure 1: Data simulation process. (a) Denoising the raw GT image; (b)
Generating the simulated measurements and back-projected image.

###  III-C Density compensation

During the measurement process of radial sampling, large intensities are
allocated to the densely sampled region at the centre of the k-space, which
leads to a back-projected image with abnormally large values. Density
Compensation (DC), as proposed in [15], calculates factors in k-space that
evenly weigh different sample locations by iteratively applying the
interpolation matrix and the adjoint interpolation matrix over $m$ iterations.
In our implementation, $m$ is set to be 10 following [4]. After obtaining the
DC weights denoted by $\bm{d}$, we multiply the k-space data by $\bm{d}$
before the back-projection.

###  III-D Acceleration factor

In data acquisition, oversampling can be applied along each spoke without
increasing the scan time. The number of spokes, denoted by
$N_{\text{s}}\in\mathbb{Z}_{+}$, determines the sparsity of the k-space
sampling, directly related to the speed of data acquisition, but with sparser
sampling implying a more challenging inverse problem for image reconstruction.
Following [4], we define an Acceleration Factor (AF) for 2D radial sampling
based on $N_{\text{s}}$ as ${\rm AF}=\sqrt{N}/N_{\text{s}}$. In previous
studies, AF was set to be a constant, _e.g._ 4 or 8. Aiming to develop a model
adaptable for a range of AFs, and to demonstrate R2D2’s performance at high
AFs, the number of spoke is left variable during training.

###  III-E Noise

In our dataset, the highest intensity of GT images is 1 due to the
normalization. The faintest image intensity, denoted by $\sigma$, is the
standard deviation of the background Gaussian noise in the image domain, which
is predicated on the assumption that all predictable intensities surpass the
noise level. The Dynamic Range (DR) is defined as the ratio between the
intensities of the maximum and faintest features, which in our dataset is
specifically the reciprocal of $\sigma$. We configured the DR to span from
$10$ to $10^{4}$, covering both low- and high-noise samples, facilitating
robust evaluation of imaging methods. Following [16], we model the
relationship between the additive noise in the k-space and its back-projection
in the image domain as $\tau=\sqrt{2{\rm L^{2}/L_{p}}}\sigma_{n}$, where $\rm
L$ and $\rm L_{\text{p}}$ are the spectral norm of the measurement operator
when the weighting is applied once and twice, respectively. Setting
$\sigma_{n}=\sigma$ ensures that the Gaussian noise in the synthetic
measurements preserves the DR of the GT images.

|  |  |  |  |  |
---|---|---|---|---|---|---
|  |  |  |  |  |
GT and Back-projected |  AIRI: (15.55, 18.96) |  U-Net: (14.93, 17.88) |  NC-PDNet: (14.40, 17.65) |  R2D2-Net (FFT): (17.46, 22.85) |  R2D2: (18.44, 23.15) |  R2D2-Net (NUFFT): (19.39, 25.10)
Figure 2: Reconstructed images in the first row and corresponding differences
compared to the GT image in the second row by different methods for one of the
validation samples. Here, DR = 167 and AF = 16. SNR and logSNR values (dB) are
displayed below.

###  III-F Details of the data generation

In training, we randomly selected $N_{\text{s}}$ from the integer set
$\\{8,9,\dots,79,80\\}$, corresponding to the AFs from 4 to 40, to generate
k-space trajectories. We simulated k-space measurements using these radial
trajectories and added corresponding noise as shown in Fig. 1 (b). For each GT
image, we generated one back-projected image with a randomly selected sampling
pattern to construct pairs of samples for supervised learning. In testing, we
randomly selected 20 GT images from the validation dataset and then applied
sampling patterns that vary the number of spokes from 10 to 80 with a step of
10, resulting in a total of 160 inverse problems.

##  IV Experimental Results

###  IV-A Comparison methods

We compare the proposed methods, including R2D2, R2D2-Net (NUFFT), and
R2D2-Net (PSF), to the state-of-the-art imaging methods, from the vanilla DNN
baseline U-Net [12], to the unrolled DNN baseline NC-PDNet [4], and to the
advanced PnP algorithm AIRI111AIRI is a PnP algorithm based on a Forward-
Backward optimisation structure. Its denoisers are trained to satisfy firm a
nonexpansiveness constraint necessary to PnP converge [17]. It leverages a
shelf of denoisers trained for a range noise levels, which are used adaptively
across the iteration process. The combination of its convergence guarantees
and its adaptive noise level functionality was shown to deliver precision and
robustness superior to more basic PnP approaches. For the purpose of this
comparison, the denoisers were trained on the same GT dataset as the DNNs of
the R2D2 family., also recently introduced in astronomical imaging [10, 18].

###  IV-B Implementation details

The DNN models were implemented using PyTorch [19]. We used the PyTorch
implementation of NUFFT from [20], which benefits from GPU acceleration. We
utilized the Adam optimizer [21] with a basic learning rate of 0.0001. In line
with the previous work [8], the number of output channels of the first
convolution layer for U-Net is set as 64. All DNN models, including the
denoisers in AIRI, were trained using a single NVIDIA Tesla V100-SXM2-16GB GPU
with a batch size of 4 on Cirrus 4, a high-performance computing system. They
were initialized randomly and trained until convergence. The same GPU was
utilized in reconstruction for DNN models, while AIRI’s iterative
reconstruction process employed the GPU combined with the dual Intel 18-core
Xeon E5-2695 processor on Cirrus 4.

(a) SNR

(b) logSNR

(c) SNR

(d) logSNR

Figure 3: The reconstruction performances of considered methods across (a),
(b) the number of iterations, and (c), (d) the number of spokes. The scalable
methods are denoted by continuous lines, while non-scalable methods are
denoted by dotted lines.

###  IV-C Evaluation metrics

We adopt linear and logarithmic signal-to-noise ratio metrics to evaluate
imaging quality. Firstly, Signal-to-Noise Ratio (SNR) is defined as
$\operatorname{SNR}(\bm{x},\bm{\bar{x}})=20\log_{10}\left(\|\bm{\bar{x}}\|_{2}/\|\bm{\bar{x}}-\bm{x}\|_{2}\right)$,
where $\|\cdot\|_{2}$ is the L2 norm. Secondly, based on the logarithmic
mapping of the images parameterized by $a>0$ is
$\operatorname{rlog}\colon\bm{x}\mapsto\operatorname{log}_{a}(a\bm{x}+1)$, the
logarithmic SNR (logSNR) to evaluate the ability of recovering faint signals
is defined as
$\operatorname{logSNR}(\bm{x},\bm{\bar{x}})=\operatorname{SNR}(\operatorname{rlog}(\bm{x}),\operatorname{rlog}(\bm{\bar{x}}))$,
where $a$ is set as the DR of the image.

###  IV-D Results

Fig. 3 (a) and (b) depict the results in terms of SNR and logSNR,
respectively, averaged from the 160 inverse problems. We categorize the
considered methods into two groups: scalable and non-scalable, based on their
capability to be applied to large-dimensional scenarios. As the number of
iterations increases, the performance exhibits a rising and converging trend
for R2D2, surpassing all other algorithms aside R2D2-Net (NUFFT). Fig. 3 (c)
and (d) show the results across the AFs, averaged from the 20 inverse problems
for each sampling pattern. NC-PDNet performs poorly in high-AF scenarios due
to the simple subnetwork structure. As expected, R2D2-Net (NUFFT) delivers
superior outcomes compared to R2D2 due to the joint training of its networks.
However, its non-scalability to large dimensions renders it akin to an oracle.

TABLE I: Additional acceleration ratios for R2D2. Algorithm$\backslash$SNR | 16 | 21 | 26 | Scalability
---|---|---|---|---
U-Net | 1.92 | 2.17 | - | Yes
AIRI | 1.69 | 1.21 | 1.35 | No
R2D2-Net (FFT) | 1.15 | 1.14 | 1.35 | Yes
NC-PDNet | 1.54 | 1.14 | 0.98 | No
R2D2-Net (NUFFT) | 0.77 | 0.76 | 0.79 | No

Based on Fig. 3 (c), we provide the additional acceleration ratio for R2D2 to
show its capability of reducing scan time compared to other methods in Tab. I.
The ratio is defined as $\text{AF}_{\text{R2D2}}/\text{AF}_{\text{Algo}}$,
where $\text{AF}_{\text{R2D2}}$ is the AF of R2D2 and
$\text{AF}_{\text{Algo}}$ is the AF of another algorithm at the same imaging
quality (SNR). In Fig. 2, R2D2 and R2D2-Net (NUFFT) provide the best
reconstruction results in terms of metrics and visual performance, while only
subtle differences between R2D2 and R2D2-Net (NUFFT) can be observed in the
zoomed areas.

TABLE II: The number of parameters, TTs, ITs, and the number of iterations for the different methods considered. Algorithm | Par. (M) | TT (h) | IT (s) | Iteration
---|---|---|---|---
AIRI | $0.6$ | $48$ | $1857\pm 387$ | $616\pm 138$
U-Net | $31$ | $52$ | $0.053\pm 0.007$ | $1$
R2D2-Net (FFT) | $248$ | $152$ | $0.129\pm 0.004$ | $1$
R2D2 | $248$ | $140$ | $0.237\pm 0.005$ | $8$
NC-PDNet | $1.6$ | $230$ | $0.263\pm 0.011$ | $1$
R2D2-Net (NUFFT) | $248$ | $315$ | $0.237\pm 0.005$ | $1$

Tab. II provides the number of trainable parameters, Training Times (TTs),
Inference Times (ITs) and number of iterations, averaged from the 160 inverse
problems. We classify the considered methods into three categories: the PnP
algorithm, scalable DNN models and non-scalable DNN models, delineated by a
horizontal line in the table. The TTs of scalable DNNs are notably shorter
compared to non-scalable ones. Specifically, when comparing the TT of R2D2-Net
(NUFFT) to R2D2-Net (FFT), we observe an increase of approximately 163 hours,
primarily due to the NUFFT-related computation during training. In terms of
ITs, R2D2 is four orders of magnitude faster than AIRI, attributed to its
smaller number of iterations.

##  V Conclusion

We have introduced the R2D2 deep learning-based image reconstruction paradigm
to MRI, leveraging a series of end-to-end DNNs in a “Matching Pursuit”
flavour. Each network in the series utilizes the back-projected data residual
and previous estimate to enhance the reconstruction. Two unrolled variants
were proposed: R2D2-Net (NUFFT), utilizing NUFFT-based data-consistency
layers, and R2D2-Net (FFT), employing the FFT approximation for data
consistency. For this proof of concept, and without loss of generality, we
have concentrated on the single-coil scenario to facilitate a comparison
between scalable and non-scalable methods. Preliminary simulations on
magnitude MR images with single-coil radial sampling demonstrate that R2D2
achieves: (i) suboptimal performance compared to its unrolled counterpart
R2D2-Net (NUFFT), which is non-scalable due to embedded NUFFT-based layers;
(ii) superior reconstruction quality to R2D2-Net (FFT); (iii) superior
reconstruction quality to state-of-the-art methods, including AIRI and NC-
PDNet.

Future research includes: (i) studying more evolved incarnations than those
considered, in particular investigating advanced DNN architectures in lieu of
R2D2’s U-Net core module; (ii) accounting for the complex-valued nature of
real-world MRI data; (iii) confirming R2D2’s practical performance in large-
dimensional scenarios, including multi-coil settings, and for 3D and 4D MRI.
It is worth mentioning that preliminary experiments suggest that embedding the
NUFFT into R2D2-Net already becomes impractical for training in a 2D multi-
coil setting with 32 coils at the same image sizes as those considered here.

## References

  * [1] Chang Min Hyun, Hwa Pyung Kim, Sung Min Lee, Sungchul Lee, and Jin Keun Seo,  “Deep learning for undersampled mri reconstruction,”  Physics in Medicine & Biology, vol. 63, no. 13, pp. 135007, 2018.
  * [2] Katherine L Wright, Jesse I Hamilton, Mark A Griswold, Vikas Gulani, and Nicole Seiberlich,  “Non-cartesian parallel imaging reconstruction,”  Journal of Magnetic Resonance Imaging, vol. 40, no. 5, pp. 1022–1040, 2014.
  * [3] Ulugbek S Kamilov, Charles A Bouman, Gregery T Buzzard, and Brendt Wohlberg,  “Plug-and-play methods for integrating physical and learned models in computational imaging: Theory, algorithms, and applications,”  IEEE Signal Processing Magazine, vol. 40, no. 1, pp. 85–97, 2023.
  * [4] Zaccharie Ramzi, GR Chaithya, Jean-Luc Starck, and Philippe Ciuciu,  “Nc-pdnet: A density-compensated unrolled network for 2d and 3d non-cartesian mri reconstruction,”  IEEE Transactions on Medical Imaging, vol. 41, no. 7, pp. 1625–1638, 2022.
  * [5] Zoran Stankovic, Bradley D Allen, Julio Garcia, Kelly B Jarvis, and Michael Markl,  “4d flow imaging with mri,”  Cardiovascular Diagnosis and Therapy, vol. 4, no. 2, pp. 173, 2014.
  * [6] Stéphane G Mallat and Zhifeng Zhang,  “Matching pursuits with time-frequency dictionaries,”  IEEE Transactions on Signal Processing, vol. 41, no. 12, pp. 3397–3415, 1993.
  * [7] Amir Aghabiglou, Matthieu Terris, Adrian Jackson, and Yves Wiaux,  “Deep network series for large-scale high-dynamic range imaging,”  in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2023, pp. 1–5.
  * [8] Amir Aghabiglou, Chung San Chu, Adrian Jackson, Arwa Dabbech, and Yves Wiaux,  “Ultra-fast high-dynamic range imaging of cygnus a with the r2d2 deep neural network series,”  arXiv preprint arXiv:2309.03291, 2023.
  * [9] Amir Aghabiglou, Chung San Chu, Arwa Dabbech, and Yves Wiaux,  “The r2d2 deep neural network series paradigm for fast precision imaging in radio astronomy,”  arXiv preprint arXiv:2403.05452, 2024.
  * [10] Matthieu Terris, Arwa Dabbech, Chao Tang, and Yves Wiaux,  “Image reconstruction algorithms in radio interferometry: From handcrafted to learned regularization denoisers,”  Monthly Notices of the Royal Astronomical Society, vol. 518, no. 1, pp. 604–622, 09 2022.
  * [11] Anuroop Sriram, Jure Zbontar, Tullie Murrell, Aaron Defazio, C Lawrence Zitnick, Nafissa Yakubova, Florian Knoll, and Patricia Johnson,  “End-to-end variational networks for accelerated mri reconstruction,”  in International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI). Springer, 2020, pp. Part II 23, 64–73.
  * [12] Olaf Ronneberger, Philipp Fischer, and Thomas Brox,  “U-net: Convolutional networks for biomedical image segmentation,”  in International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI). Springer, 2015, pp. Part III 18, 234–241.
  * [13] Jure Zbontar, Florian Knoll, Anuroop Sriram, Tullie Murrell, Zhengnan Huang, Matthew J Muckley, Aaron Defazio, Ruben Stern, Patricia Johnson, Mary Bruno, et al.,  “fastmri: An open dataset and benchmarks for accelerated mri,”  arXiv preprint arXiv:1811.08839, 2018.
  * [14] Kai Zhang, Yawei Li, Jingyun Liang, Jiezhang Cao, Yulun Zhang, Hao Tang, Deng-Ping Fan, Radu Timofte, and Luc Van Gool,  “Practical blind image denoising via swin-conv-unet and data synthesis,”  Machine Intelligence Research, vol. 20, no. 6, pp. 822–836, 2023.
  * [15] James G Pipe and Padmanabhan Menon,  “Sampling density compensation in mri: rationale and an iterative numerical solution,”  Magnetic Resonance in Medicine, vol. 41, no. 1, pp. 179–186, 1999.
  * [16] Amanda G Wilber, Arwa Dabbech, Adrian Jackson, and Yves Wiaux,  “Scalable precision wide-field imaging in radio interferometry: I. usara validated on askap data,”  Monthly Notices of the Royal Astronomical Society, vol. 522, no. 4, pp. 5558–5575, 2023.
  * [17] Jean-Christophe Pesquet, Audrey Repetti, Matthieu Terris, and Yves Wiaux,  “Learning maximally monotone operators for image recovery,”  SIAM Journal on Imaging Sciences, vol. 14, no. 3, pp. 1206–1237, 2021.
  * [18] Matthieu Terris, Chao Tang, Adrian Jackson, and Yves Wiaux,  “Plug-and-play imaging with model uncertainty quantification in radio astronomy,”  arXiv preprint arXiv:2312.07137, 2023.
  * [19] Adam Paszke, Sam Gross, Francisco Massa, Lerer, et al.,  “Pytorch: An imperative style, high-performance deep learning library,”  Advances in neural information processing systems (NIPS), pp. 8024–8035, 2019.
  * [20] M. J. Muckley, R. Stern, T. Murrell, and F. Knoll,  “Torchkbnufft: A high-level, hardware-agnostic non-uniform fast fourier transform,”  in ISMRM Workshop on Data Sampling & Image Reconstruction, 2020.
  * [21] Diederik P Kingma and Jimmy Ba,  “Adam: A method for stochastic optimization,”  arXiv preprint arXiv:1412.6980, 2014.


# Bright X-ray radiation from plasma bubbles in an evolving laser wakefield
accelerator

M.S. Bloom  The John Adams Institute for Accelerator Science, Imperial College
London, London, SW7 2AZ, UK  M.J.V. Streeter  The John Adams Institute for
Accelerator Science, Imperial College London, London, SW7 2AZ, UK  S. Kneip
The John Adams Institute for Accelerator Science, Imperial College London,
London, SW7 2AZ, UK  R.A. Bendoyro  GoLP/Instituto de Plasmas e Fusão Nuclear,
Instituto Superior Técnico, U.L., Lisboa 1049-001, Portugal  O. Cheklov
Central Laser Facility, STFC Rutherford Appleton Laboratory, Didcot OX11 0QX,
UK  J. M. Cole  The John Adams Institute for Accelerator Science, Imperial
College London, London, SW7 2AZ, UK  A. Döpp  The John Adams Institute for
Accelerator Science, Imperial College London, London, SW7 2AZ, UK  Fakultät
für Physik, Ludwig-Maximilians-Universität München, D-85748 Garching, Germany
Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Str. 1, D-85748
Garching, Germany  C.J. Hooker  Central Laser Facility, STFC Rutherford
Appleton Laboratory, Didcot OX11 0QX, UK  J. Holloway  Department of Physics
and Astronomy, University College London, London WC1E 6BT, UK  J. Jiang
GoLP/Instituto de Plasmas e Fusão Nuclear, Instituto Superior Técnico, U.L.,
Lisboa 1049-001, Portugal  N.C. Lopes  The John Adams Institute for
Accelerator Science, Imperial College London, London, SW7 2AZ, UK
GoLP/Instituto de Plasmas e Fusão Nuclear, Instituto Superior Técnico, U.L.,
Lisboa 1049-001, Portugal  H. Nakamura  The John Adams Institute for
Accelerator Science, Imperial College London, London, SW7 2AZ, UK  P.A.
Norreys  Central Laser Facility, STFC Rutherford Appleton Laboratory, Didcot
OX11 0QX, UK  P.P. Rajeev  Central Laser Facility, STFC Rutherford Appleton
Laboratory, Didcot OX11 0QX, UK  D.R. Symes  Central Laser Facility, STFC
Rutherford Appleton Laboratory, Didcot OX11 0QX, UK  J. Schreiber  The John
Adams Institute for Accelerator Science, Imperial College London, London, SW7
2AZ, UK  Fakultät für Physik, Ludwig-Maximilians-Universität München, D-85748
Garching, Germany  Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Str.
1, D-85748 Garching, Germany  J. C. Wood  The John Adams Institute for
Accelerator Science, Imperial College London, London, SW7 2AZ, UK  M. Wing
Department of Physics and Astronomy, University College London, London WC1E
6BT, UK  Z. Najmudin  The John Adams Institute for Accelerator Science,
Imperial College London, London, SW7 2AZ, UK  S.P.D. Mangles  The John Adams
Institute for Accelerator Science, Imperial College London, London, SW7 2AZ,
UK<EMAIL_ADDRESS>
###### Abstract

We show that the properties of the electron beam and bright X-rays produced by
a laser wakefield accelerator can be predicted if the distance over which the
laser self-focuses and compresses prior to self-injection is taken into
account. A model based on oscillations of the beam inside a plasma bubble
shows that performance is optimised when the plasma length is matched to the
laser depletion length. With a 200 TW laser pulse this results in an X-ray
beam with median photon energy of 20 keV, $>10^{9}$ photons per shot and a
peak brightness of $4\times
10^{23}\,\mathrm{photons~{}s^{-1}mrad^{-2}mm^{-2}(0.1\%BW)^{-1}}$.

††preprint: APS/123-QED

‘

Laser wakefield acceleratorsTajima and Dawson (1979) have gathered increasing
interest since it was first shown that they were capable of producing high
quality electron beamsMangles  _et al._ (2004); Geddes  _et al._ (2004); Faure
_et al._ (2004). Development has continued apace, and laser wakefield
accelerators can now produce ultra-short bunches of electrons, down to a few
femtosecondsLundh  _et al._ (2011), and reach multi-GeV beam energiesLeemans
_et al._ (2014); Wang  _et al._ (2013); Kim  _et al._ (2013). One of the
principal near-term applications of laser wakefield accelerators is the
production of bright, femtosecond duration pulses of broadband X-raysRousse
_et al._ (2004); Kneip  _et al._ (2010), that are suitable for a range of
applicationsAlbert  _et al._ (2014).

A laser wakefield accelerator is formed when an intense, short-duration laser
pulse is fired into a moderate density plasma. The ponderomotive force
associated with the laser pushes plasma electrons out of its way as it
propagates. Due to the much heavier–effectively stationary–ions the electrons
rush back towards their equilibrium position once the laser has passed,
setting up a charge oscillation with a relativistic phase velocity in the wake
of the laser pulse. When driven by a sufficiently intense laser pulse, almost
all of the electrons can be expelled from an approximately spherical cavity
behind the drive, known as the plasma bubblePukhov and Meyer-ter Vehn (2002).
The electric fields inside this bubble can be capable of accelerating
electrons to 1 GeV in just 1 cmKneip  _et al._ (2008).

Electrons can be self-injected into the bubble from the background plasma if
the wave has a sufficiently high amplitudeMangles  _et al._ (2012). The three-
dimensional structure of the bubble means that the injected electron beam
undergoes betatron oscillations with a short wavelength ($\sim$100 µm for
$\sim 1\,\mathrm{GeV}$) and so generates synchrotron X-rays in the multi-keV
spectral rangeRousse  _et al._ (2004).

In this letter we report on the experimental optimisation of the X-rays
generated by a laser wakefield accelerator driven by a 200 TW laser. We show
that the electron and X-ray properties are well described by an analytical
model that includes the fact that self-injection only occurs after the pulse
has self-focused and self-compressed to a sufficiently high intensity. We also
show that this source outperforms previous reports of X-ray emission from
wakefield acceleratorsKneip  _et al._ (2010); Wang  _et al._ (2013).

The experiment was performed using the Astra Gemini laser at the Rutherford
Appleton Laboratory, which delivered pulses of 12 J and duration 55 fs fwhm
onto a supersonic gas jet target. The laser was focused with an $f$/20 off-
axis parabolic mirror to a 22 µm fwhm spot containing $\approx 30\%$ of the
energy. The peak intensity at focus in vacuum was $I\approx 2.2\times
10^{19}\,\mathrm{Wcm^{-2}}$, corresponding to a normalised vector potential
$a_{\rm 0}=eA_{\rm 0}/(m_{\rm e}c^{2})\simeq 3$.

Two different helium gas jet targets with exit diameters of (10, 15) mm were
used, producing approximately uniform density plasmas of length (8.5, 13) mm
with electron densities up to $(8.0,\,4.0)\times 10^{18}$ cm-3 respectively.
The laser was focused onto the front edge of the gas flow. The generated
electron beam was analysed using a magnetic spectrometer consisting of a 30 cm
long 1 T permanent dipole magnet and two scintillating (Kodak Lanex regular)
screens.

The magnet also swept the electron beam away from an indirect detection X-ray
ccd camera (Princeton Instruments PIXIS) placed on the laser axis. This was
mounted outside the vacuum chamber behind a 180 µm thick Be window. An array
of 16 metallic filters mounted on a thin Al/mylar substrate was placed
directly in front of the camera’s CsI scintillator.

The X-ray spectrum was found by performing a least squares fit to the signal
detected behind each filter, taking into account the transmission through each
filter and the absolutely calibrated detector response under the assumptions
that the spectrum was spatially uniform over the detector and has a
synchrotron-like shape given by:
$d^{2}I/(dEd\Omega)_{\theta=0}\propto\xi^{2}\mathcal{K}^{2}_{2/3}(\xi/2)$,
where $\mathcal{K}_{2/3}(x)$ is a modified Bessel function of order $2/3$ and
$\xi=E/E_{\rm c}$. The shape of this spectrum is characterized by a single
parameter, the critical energy, $E_{\rm c}$ 111This definition of $E_{\rm c}$
is different to that used in refs Esarey  _et al._ (2002); Kneip  _et al._
(2008) but consistent with refs Rousse  _et al._ (2004, 2007); Jackson
(1999)). Gaps between the X-ray filters and repeated filters allowed gradients
in the X-ray beam profile to be taken into account.

Figure 1: Variation of electron and X-ray beams as a function of $n_{\rm e}$
for 10 mm nozzle. Top: Dispersion-corrected spectrally-dispersed images of the
electron beam for a selection of shots. Bottom: X-ray CCD camera images for
the same shots.

Fig. 1 shows the variation of the electron beam spectrum and corresponding
X-ray ccd images with plasma density for the 10 mm nozzle. Shots in this data
sub-set correspond to the brightest X-ray image recorded at each density. The
data shows that there is an optimum density for acceleration of $n_{\rm
e}\approx 3.8\times 10^{18}\,\mathrm{cm^{-3}}$, and that the X-ray signal is
correlated with the electron beam energy. Above the optimum density, the
electron beam begins to develop transverse structure and increased divergence.
This is consistent with the electron beam interacting with the plasma, driving
its own wake and becoming susceptible to propagation instabilitiesHuntington
_et al._ (2011). As the electron beam dephases, it can also interact with the
laser fieldMangles  _et al._ (2006) which can also increase the X-ray
fluxCipiccia  _et al._ (2011). However, under these conditions no major
enhancement in X-ray flux is evident due to these effects.

Figure 2: Variation of $W_{\rm max}$ with $n_{\rm e}$ for 10 (squares) and 15
mm (triangles) gas jets. Each point is the mean from $N$ = 1–5 shots, bars
represent the combined statistical and measurement errors. Curves show maximum
energy predicted using: eqn. 1 (dot-dash); eqn. 3 (dash); eqn. 4 for the 10 mm
(blue); and 15 mm (red) nozzles. The regions (I, II and III) described in the
text are indicated for each nozzle. Region II for each nozzle is also
indicated as a shaded region.

Fig. 2 shows the variation of the electron energy, $W_{\rm max}$ with $n_{\rm
e}$ for both the 10 and 15 mm gas nozzles. Both datasets show a similar trend;
above a threshold density, $W_{\rm max}$ rapidly increases until it reaches a
maximum after which it decreases approximately $\propto\frac{1}{n_{\rm e}}$.
Both the threshold and optimum densities occur at lower densities for the 15
mm nozzle. As the X-ray signal is optimised when the electron beam energy is
highest, optimisation of the X-ray generation relies on understanding how to
optimise the electron energy. The maximum energy predicted by Lu et al.,Lu
_et al._ (2007):

$\displaystyle W_{\rm max}/{\left(m_{\rm e}c^{2}\right)}$
$\displaystyle={\textstyle\frac{2}{3}}a_{\rm 0}({n_{\rm c}/{n_{\rm e}}})\,,$
(1)

is plotted in fig. 2 (dot-dash line) using the $a_{\rm 0}$ for vacuum. Here
$n_{\rm c}=\epsilon_{\rm 0}m_{\rm e}\omega^{2}/e^{2}$ is the critical plasma
density for radiation with angular frequency $\omega$. Eqn. 1 consistently
underestimates the electron energy at high densities and does not predict an
optimum density.

Eqn. 1 assumes a non-evolving laser pulse. However, self-focusing and pulse
compression occur as the pulse propagates in its self-generated wakefield,
causing $a_{\rm 0}$ to increase. The maximum value of $a_{\rm 0}$ can be
calculated by assuming that the pulse evolves to a matched spot size that
satisfies $a_{\rm 0}\approx 2(P_{\rm f}/P_{\rm c})^{1/3}$ Lu  _et al._ (2007),
where $P_{\rm c}=2m_{\rm e}c^{3}n_{\rm c}/(r_{\rm e}n_{\rm e})$ is the
critical power for relativistic self-focusing and $r_{\rm e}$ is the classical
electron radius. The final power, $P_{\rm f}$, depends on pulse compression
and the amount of laser energy that is trapped in the bubble. For a fraction
$\alpha$ of the total laser energy $E_{\rm L}$ compressed to a pulse duration
$\tau_{\rm f}$, the resulting expression is:

$\displaystyle a_{\rm max}$ $\displaystyle\simeq 2\left(\alpha E_{\rm
L}/{\tau_{\rm f}P_{\rm c}}\right)^{\frac{1}{3}}\,.$ (2)

Even though this expression does not include the effects of pump depletion and
photon deceleration, it has been shown to be sufficiently accurate to predict
the self-injection thresholdMangles  _et al._ (2012). The pulse compression
can be quantified by assuming that the front of the pulse travels at the
linear group velocity in the plasma while the rear of the pulse, which sits in
the significantly reduced plasma density inside the bubble, travels at $c$.
For an initial pulse duration $\tau_{\rm 0}$, the pulse duration after
propagation length $l$ is then: $\tau_{\rm f}(l)\approx\tau_{\rm
0}-\frac{n_{\rm e}l}{2cn_{\rm c}}$Schreiber  _et al._ (2010). The maximum
propagation length, $L_{\rm max}$, will be limited by pump depletion $L_{\rm
pump}\approx c\tau_{\rm 0}\,n_{\rm c}/n_{\rm e}$ Lu  _et al._ (2007) or by the
length of the target, $L_{\rm target}$, if $L_{\rm target}<L_{\rm pump}$. So
accounting for pulse evolution, the beam energy varies as:

$\displaystyle W_{\rm max}^{\prime}/(m_{\rm e}c^{2})$
$\displaystyle\approx{\textstyle\frac{4}{3}}\left(\alpha E_{\rm L}/\tau_{\rm
f}P_{\rm c}\right)^{\frac{1}{3}}(n_{\rm c}/n_{\rm e})$ (3)

Eqn. 3 is shown as the dashed line in fig. 2. This model overestimates the
observed energy gain, only approaching the data at high densities. It also
still fails to explain the initial increase in beam energy with increasing
density. These features can be explained by including a distance over which
the $a_{\rm 0}$ amplification occurs before self-injection. We call this the
pre-injection pulse evolution length, $L_{\rm PIPE}$. This pipe length will
decrease at higher densities as the pulse evolution rates increaseMori (1997).
The variation of electron energy with density can be split into three regions.
At low density (region I), the pipe length is longer than the gas jet and so
no electrons are injected.

As the density is increased, the evolution becomes fast enough that injection
occurs before the end of the gas jet, resulting in low energy electron beams.
This region II is marked by the shaded area in fig. 2. Increasing the density
further reduces the pipe length and brings the injection point earlier in the
gas jet. But in region II, the density is low enough that the laser has not
depleted by the end of the gas jet, so earlier injection leads to an increase
in acceleration length. This coupled with the increase in the accelerating
field ($\propto\sqrt{n_{\rm e}}$) results in a rapid increase of beam energy
with increasing density.

Once the density is high enough that the pump depletion length is less than
the gas jet length (region III), increasing the density actually decreases the
length over which the electron beam is accelerated. Therefore in region III,
despite the continuing increase in electric field strength, higher density
results in a decrease in the electron energy. At the highest densities, the
pipe length is very short and, because the pump depletion length is longer
than the dephasing length, the electron beam energy approaches the dephasing
limit (eqn. 3).

To model this behaviour, we assume that the pipe length $L_{\rm PIPE}=Sn_{\rm
c}/n_{\rm e}$, where $S$ is a constant of proportionality to be determined.
The fraction of the dephasing length, $L_{\phi}=\frac{4}{3}a_{\rm
0}^{1/2}(n_{\rm c}/n_{\rm e})c/\omega_{\rm p}$, that the acceleration occurs
over is therefore $\Delta_{\rm acc}=(L_{\rm max}-L_{\rm PIPE})/L_{\phi}$.

In the bubble regime, the electric field strength varies linearly inside the
bubble, so that the beam energy varies quadratically according to:

$\displaystyle W_{\rm max}^{\prime\prime}=W_{\rm
max}^{\prime}\left(2\Delta_{\rm acc}-\Delta_{\rm acc}^{2}\right).$ (4)

Eqn. 4 is plotted in fig. 2 (solid red and blue lines) with a single fitting
parameter, $S=$11.9 µm, chosen to best reproduce the experimental trend using
a least squares fit. A detailed study of the laser pulse properties after the
interaction show that the length over which power amplification due to pulse
compression occurs is closely matched to this pipe lengthStreeter  _et al._
(2017). The pipe model reproduces all of the main features of the
experimentally observed variation in electron energy including: the rapid rise
to an optimum density; the reduction in optimum density for the longer nozzle;
and the slower fall-off of maximum energy at higher densities. The model
performs equally well for both the 10 and 15 mm nozzles with the same value
for $S$, indicating that the precise shape of the gas density plays a minor
role in determining the evolution as compared to the initial laser parameters.
One of the main results of this model is that the density which maximises the
electron beam energy for a fixed length target occurs when the length of the
target equals the pump depletion length.

Fig. 3 plots how the X-ray critical energy, $E_{\rm c}$, and peak brightness,
$\mathcal{B}_{\rm 0}$, vary with $n_{\rm e}$ for the 10 mm gas nozzle. Both
the critical energy and brightness show similar behaviour; rapidly increasing
as the density is increased before turning over above $n_{\rm e}\approx
3.8\times 10^{18}\,\mathrm{cm^{-3}}$. The critical energy reaches $E_{\rm
c}\approx 30\,\mathrm{keV}$. The brightness, calculated assuming a constant
duration matching that of the laser pulse (55 fs) and an r.m.s. source radius
consistent with previous measurements (1 µm), reaches $4\times
10^{23}\,\mathrm{photons~{}s^{-1}mrad^{-2}mm^{-2}(0.1\%BW)^{-1}}$. This is
significantly higher than previous results at lower laser power, e.g. Kneip
_et al._ (2010), primarily due to the ability to accelerate electrons to
$\sim$ GeV energies, thus increasing their ability to radiate. Under optimum
conditions, a total photon yield of $N_{\rm X}\approx 3\times 10^{9}$ per shot
was measured.

Figure 3:  Observed X-ray: a) $E_{\rm c}$ and b) $\mathcal{B}_{\rm 0}$ as a
function of $n_{\rm e}$ for the 10 mm nozzle. Each point is the mean of $N=$
1–5 shots, bars represent the combined statistical and measurement errors.
Curves show predictions using eqn. 6.

Calculations of the radiation emitted by self-injected electrons oscillating
in the plasma bubble as they are accelerated and then decelerated from the
back of the bubble to the front, i.e. over a distance $L=2L_{\rm\phi}$, have
been presented in Thomas (2010). The number of betatron oscillations is found
to be, $N_{\beta}=\gamma_{\rm p}$, where $\gamma_{\rm p}$ is the Lorentz
factor associated with the bubble motion. The resulting X-ray spectrum is
synchrotron-like with an enhanced high-energy tail. To model the X-ray
emission for an arbitrary acceleration length less than $2L_{\phi}$, the
expressions from Thomas (2010) can readily be modified. The modified spectrum
for a beam of $N_{\rm b}\approx 3.1\times 10^{8}\lambda_{\rm 0}\sqrt{P_{\rm
f}}$ electronsKneip  _et al._ (2012) is:

$\displaystyle\frac{{\rm d}^{2}I}{{\rm d}\omega{\rm d}\Omega}$
$\displaystyle=\gamma_{\rm p}\frac{3e^{2}}{\pi^{3}\epsilon_{\rm 0}c}N_{\rm
b}\left(\frac{W_{\rm max}^{\prime}}{m_{\rm
e}c^{2}}\right)^{2}\mathcal{A}\left(\frac{\omega}{2\omega_{\rm
c0}},\Delta_{\rm acc}\right)\,,$ (5)

where $E_{\rm c0}=\hbar\omega_{\rm c0}=7.75\times 10^{-4}a_{\rm
0}^{5/2}(n_{\rm c}/n_{\rm e})^{9/8}$ is the critical energy of a synchrotron
spectrum corresponding to that emitted by an electron at the maximum energy
$W_{\rm max}^{\prime}$. The spectral shape function,
$\mathcal{A}(\xi,\Delta_{\rm acc})$ is:

$\displaystyle\mathcal{A}$ $\displaystyle=\xi^{2}\int_{-1}^{\Delta_{\rm
acc}-1}(1-x^{2})^{-\frac{3}{2}}\mathcal{K}_{2/3}^{2}\left(\xi(1-x^{2})^{-\frac{7}{4}}\right){\rm
d}x\,.$ (6)

The calculated variation in $E_{\rm c}$ and $\mathcal{B}_{\rm 0}$ have been
overlaid on the experimental data in fig. 3a & b, respectively, for a plasma
length of 8.5 mm with $L_{\rm PIPE}=Sn_{\rm c}/n_{\rm e}$ and $S=$11.9 µm. The
model, which has no fitting parameters other than than the pipe length
determined from the electron data, is consistent with the data. It exhibits
the same the sharp initial rise and slow fall off for both $E_{\rm c}$ and
$\mathcal{B}_{\rm 0}$. The model shows that the X-ray emission is also
optimised for $L_{\rm PIPE}+L_{\rm acc}\simeq L_{\rm target}$.

Fig. 4 compares our measured $E_{\rm c}$ and $\mathcal{B}_{\rm 0}$ with
previous experiments. Kneip et al. Kneip  _et al._ (2012) calculated the
following scaling laws in terms of the laser power, $P$ (in TW) for both
$E_{\rm c}$ (in keV) and $\mathcal{B}_{\rm 0}$ (in
$\mathrm{photons~{}s^{-1}mrad^{-2}mm^{-2}(0.1\%BW)^{-1}}$):

$\displaystyle E_{\rm c}$ $\displaystyle=3.4\times 10^{-2}\left(\beta\alpha
P\right)^{\frac{5}{6}}\left(n_{\rm c}/n_{\rm e}\right)^{\frac{7}{24}}\,,$ (7)
$\displaystyle\mathcal{B}_{\rm 0}$ $\displaystyle=2.7\times
10^{19}\left(\beta\alpha P\right)^{\frac{2}{3}}\left(n_{\rm c}/n_{\rm
e}\right)^{\frac{19}{12}}\,,$ (8)

where a power amplification factor, $\beta$ is included to account for pulse
compression. Compression over a depletion length at the rate in Schreiber  _et
al._ (2010) gives $\beta=2$.

Eqns. 7 and 8 treat $P$ and $n_{\rm e}$ as independent parameters. However,
increasing $P$ decreases the threshold density for self-injection, $n_{\rm
th}$, so experiments at higher power typically operate at lower densityMangles
(2016). Using Mangles  _et al._ (2012) to find $n_{\rm th}(P)$, we can
therefore eliminate $n_{\rm e}$ from the above expressions. Furthermore
experiments are typically optimised just above the threshold, so the curves
for $E_{c}$ and $\mathcal{B}_{\rm 0}$ (plotted in fig. 4) span the range
$[n_{\rm th},3n_{\rm th}]$.

The scaling laws and experimental data both show that higher power lasers
produce higher energy and brighter X-ray sources. However it is important to
note that experiments can differ significantly from the scaling law. For
example, higher photon energies can be produced if asymmetries are present in
the wake which increase the betatron oscillation amplitude and so increase the
source size and can decrease the brightnessMangles  _et al._ (2009).

The photon source in our experiments is much brighter, and at higher photon
energy than those at lower powerKneip  _et al._ (2010) and comparable to an
unoptimised experiment at five times higher laser powerWang  _et al._ (2013).
Our experimental data and the model developed here suggest that a significant
increase in the flux of X-ray radiation from a laser wakefield accelerator can
be achieved by careful optimisation of the length and plasma density for a
given laser system, and so greatly broaden the range of applications that are
accessible with theses sourcesAlbert  _et al._ (2014).

Figure 4: Scaling of X-ray radiation with laser power from various experiments
Mangles  _et al._ (2009); Rousse  _et al._ (2004); Schnell  _et al._ (2012);
Kneip  _et al._ (2010); Wenz  _et al._ (2015); Fourmaux  _et al._ (2011);
Albert  _et al._ (2013): a) $E_{\rm c}$; b) $\mathcal{B}_{\rm 0}$. The shaded
region corresponds to eqns. 7 and 8 for $n_{\rm th}\leq n_{\rm e}\leq 3n_{\rm
th}$.

The authors thank the staff at the Central Laser Facility, Rutherford Appleton
Laboratory for their assistance. This research was supported by STFC
(ST/J002062/1, ST/P000835/1), EPSRC (EP/I014462/1), and the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (Grant Agreement No. 682399). The authors confirm that the all data
used in this study are available without restriction. Data can be obtained by
contacting<EMAIL_ADDRESS>
## References

  * Tajima and Dawson (1979) T. Tajima and J. M. Dawson, Physical Review Letters 43, 267 (1979).
  * Mangles  _et al._ (2004) S. P. D. Mangles, C. D. Murphy, Z. Najmudin, A. G. R. Thomas, _et al._ , Nature 431, 535 (2004).
  * Geddes  _et al._ (2004) C. G. R. Geddes, C. Toth, J. van Tilborg, E. Esarey, C. B. Schroeder, D. Bruhwiler, C. Nieter, J. Cary, and W. P. Leemans, Nature 431, 538 (2004).
  * Faure  _et al._ (2004) J. Faure, Y. Glinec, A. Pukhov, S. Kiselev, S. Gordienko, E. Lefebvre, J. P. Rousseau, F. Burgy, and V. Malka, Nature 431, 541 (2004).
  * Lundh  _et al._ (2011) O. Lundh, J. Lim, C. Rechatin, L. Ammoura, A. Ben-Ismaïl, X. Davoine, G. Gallot, J.-P. Goddet, E. Lefebvre, V. Malka, and J. Faure, Nature Physics 7, 219 (2011).
  * Leemans  _et al._ (2014) W. P. Leemans, A. J. Gonsalves, H. S. Mao, K. Nakamura, C. Benedetti, C. B. Schroeder, C. Tóth, J. Daniels, D. E. Mittelberger, S. S. Bulanov, J. L. Vay, C. G. R. Geddes, and E. Esarey, Physical Review Letters 113, 245002 (2014).
  * Wang  _et al._ (2013) X. Wang, R. Zgadzaj, N. Fazel, Z. Li, _et al._ , Nature Communications 4, 1988 (2013).
  * Kim  _et al._ (2013) H. T. Kim, K. H. Pae, H. J. Cha, I. J. Kim, T. J. Yu, J. H. Sung, S. K. Lee, T. M. Jeong, and J. Lee, Physical Review Letters 11, 165002 (2013).
  * Rousse  _et al._ (2004) A. Rousse, K. T. Phuoc, R. C. Shah, A. Pukhov, E. Lefebvre, V. Malka, S. Kiselev, F. Burgy, J. P. Rousseau, D. P. Umstadter, and D. Hulin, Physical Review Letters 93, 135005 (2004).
  * Kneip  _et al._ (2010) S. Kneip, C. McGuffey, J. L. Martins, Martins, _et al._ , Nature Physics (2010), 10. 1038/nphys1789.
  * Albert  _et al._ (2014) F. Albert, A. G. R. Thomas, S. P. D. Mangles, S. Banerjee, S. Corde, A. Flacco, M. Litos, D. Neely, J. Vieira, Z. Najmudin, R. Bingham, C. Joshi, and T. Katsouleas, Plasma Physics and Controlled Fusion 56, 084015 (2014).
  * Pukhov and Meyer-ter Vehn (2002) A. Pukhov and J. Meyer-ter Vehn, Applied Physics B: Lasers and Optics 74, 355 (2002).
  * Kneip  _et al._ (2008) S. Kneip, S. R. Nagel, C. Bellei, N. Bourgeois, _et al._ , Physical Review Letters 100, 105006 (2008).
  * Mangles  _et al._ (2012) S. P. D. Mangles, G. Genoud, M. S. Bloom, M. Burza, Z. Najmudin, A. Persson, K. Svensson, A. G. R. Thomas, and C. G. Wahlström, Physical Review Special Topics - Accelerators and Beams 15, 011302 (2012).
  * Note (1) This definition of $E_{\rm c}$ is different to that used in refs Esarey  _et al._ (2002); Kneip  _et al._ (2008) but consistent with refs Rousse  _et al._ (2004, 2007); Jackson (1999)).
  * Huntington  _et al._ (2011) C. M. Huntington, A. G. R. Thomas, C. McGuffey, T. Matsuoka, V. Chvykov, G. Kalintchenko, S. Kneip, Z. Najmudin, C. Palmer, V. Yanovsky, A. Maksimchuk, R. P. Drake, T. Katsouleas, and K. Krushelnick, Physical Review Letters 106, 105001 (2011).
  * Mangles  _et al._ (2006) S. P. D. Mangles, A. G. R. Thomas, M. C. Kaluza, O. Lundh, F. Lindau, A. Persson, F. S. Tsung, Z. Najmudin, W. B. Mori, C. G. Wahlström, and K. Krushelnick, Physical Review Letters 96, 215001 (2006).
  * Cipiccia  _et al._ (2011) S. Cipiccia, M. R. Islam, B. Ersfeld, R. P. Shanks, _et al._ , Nature Physics 7, 867 (2011).
  * Lu  _et al._ (2007) W. Lu, M. Tzoufras, C. Joshi, F. S. Tsung, W. B. Mori, J. Vieira, R. Fonseca, and L. Silva, Physical Review Special Topics - Accelerators and Beams 10, 061301 (2007).
  * Schreiber  _et al._ (2010) J. Schreiber, C. Bellei, S. P. D. Mangles, C. Kamperidis, S. Kneip, S. R. Nagel, C. A. J. Palmer, P. P. Rajeev, M. J. V. Streeter, and Z. Najmudin, Physical Review Letters 105, 235003 (2010).
  * Mori (1997) W. Mori, IEEE Journal of Quantum Electronics 33, 1942 (1997).
  * Streeter  _et al._ (2017) M. J. V. Streeter, S. Kneip, M. S. Bloom, R. A. Bendoyro, _et al._ , ArXiv e-prints (2017), arXiv:1710.05337 [physics.plasm-ph] .
  * Thomas (2010) A. G. R. Thomas, Physics of Plasmas 17, 056708 (2010).
  * Kneip  _et al._ (2012) S. Kneip, Z. Najmudin, and A. G. R. Thomas, High Energy Density Physics 8, 133 (2012).
  * Mangles (2016) S. P. D. Mangles, CERN Yellow Reports 1, 289 (2016).
  * Mangles  _et al._ (2009) S. P. D. Mangles, G. Genoud, S. Kneip, M. Burza, K. Cassou, B. Cros, N. P. Dover, C. Kamperidis, Z. Najmudin, A. Persson, J. Schreiber, F. Wojda, and C. G. Wahlström, Applied Physics Letters 95, 181106 (2009).
  * Schnell  _et al._ (2012) M. Schnell, A. Sävert, B. Landgraf, M. Reuter, M. Nicolai, O. Jäckel, C. Peth, T. Thiele, O. Jansen, A. Pukhov, O. Willi, M. C. Kaluza, and C. Spielmann, Physical Review Letters 108, 075001 (2012).
  * Wenz  _et al._ (2015) J. Wenz, S. Schleede, K. Khrennikov, M. Bech, P. Thibault, M. Heigoldt, F. Pfeiffer, and S. Karsch, Nature Communications 6, 7568 (2015).
  * Fourmaux  _et al._ (2011) S. Fourmaux, S. Corde, K. Ta Phuoc, P. M. Leguay, S. Payeur, P. Lassonde, S. Gnedyuk, G. Lebrun, C. Fourment, V. Malka, S. Sebban, A. Rousse, and J. C. Kieffer, New Journal of Physics 13, 033017 (2011).
  * Albert  _et al._ (2013) F. Albert, B. B. Pollock, J. L. Shaw, K. A. Marsh, J. E. Ralph, Y. H. Chen, D. Alessi, A. Pak, C. E. Clayton, S. H. Glenzer, and C. Joshi, Physical Review Letters 111, 235004 (2013).
  * Esarey  _et al._ (2002) E. Esarey, B. A. Shadwick, P. Catravas, and W. P. Leemans, Physical Review E 65, 56505 (2002).
  * Rousse  _et al._ (2007) A. Rousse, K. T. Phuoc, R. C. Shah, R. Fitour, and F. Albert, European Physical Journal D 45, 391 (2007).
  * Jackson (1999) J. D. Jackson,  _Classical Electrodynamics_ , 3rd ed. (Wiley, New York, 1999).