# A free discontinuity model for smectic thin films

John M.Ball111J. M. Ball, School of Mathematical & Computer Sciences, Heriot-
Watt University, Edinburgh, EH14 4AS, United Kingdom. E-mail:
<EMAIL_ADDRESS>Giacomo Canevari222G. Canevari, Dipartimento di
Informatica, Università di Verona, Strada le Grazie 15, 37134 Verona, Italy.
E-mail<EMAIL_ADDRESS>Bianca Stroffolini333B. Stroffolini,
Dipartimento di Matematica e Applicazioni, Università di Napoli Federico II,
Via Cintia, 80126 Napoli, Italy. E-mail<EMAIL_ADDRESS>
###### Abstract

We attempt to describe surface defects in smectic A thin films by formulating
a free discontinuity problem — that is, a variational problem in which the
order parameter is allowed to have jump discontinuities on some (unknown) set.
The free energy functional contains an interfacial energy which penalizes
dislocations of the smectic layers at the jump. We discuss mathematical issues
related to the existence of minimizers and provide examples of minimizers in
some simplified settings.

##  1 Introduction

In this paper we report on an attempt to formulate a mathematical model,
capable of rigorous mathematical analysis, to describe aspects of the
interesting experiments on smectic A thin films of the group of E. Lacaze [17,
18, 24, 23, 10]. In these experiments thin films of 8CB liquid crystal are
deposited on various substrates, for example ${\rm MoS}_{2}$, mica and rubbed
PVA. Depending on the thickness of the films and on the nature of the
substrate different configurations of the smectic layers are observed [25].
The main features are illustrated in Figs 2, 2. Viewed from above (Fig. 2)
families of parallel ‘oily streaks’ are observed, each of which consists of a
flattened hemicylinder of smectic layers (Fig. 2). The representation in Fig.
2 of the configuration of layers in an oily streak on rubbed PVA substrate is
not a direct observation, but is deduced from X-ray diffraction and
ellipsometry; indeed recently some slight modifications to the likely
configuration have been suggested by Niyonzima et al [19].

Figure 1: AFM image of families of oily streaks of 8CB on ${\rm MoS}_{2}$
substrate [17].

Figure 2: Experimental representation of cross-sectional layer configuration
in an oily streak of 8CB on PVA substrate [10].

The observed layer configurations are a consequence of the antagonistic
boundary conditions, homeotropic on the upper surface in contact with air, and
unidirectional planar anchoring on the bottom surface parallel to the
substrate. Because for smectic A the director $\mathbf{n}$ is perpendicular to
the layers, this means that the layers prefer energetically to be tangent to
the upper free surface and perpendicular to the substrate.

A feature of the observed layer configurations is the presence of surface
(wall) defects across which the layer normal, and thus $\mathbf{n}$, jumps.
This suggests that a good mathematical framework to handle this problem is the
free discontinuity setting originated by De Giorgi & Ambrosio [12] for image
segmentation and subsequently used by Francfort & Marigo for Griffith fracture
models [13]. The relevance of such a setting for liquid crystal problems was
proposed in [6]. The free discontinuity formulation uses a free energy in
which there is a competition between bulk (volumetric) energy and interfacial
energy corresponding to unknown surfaces of discontinuity. In the case of
fracture mechanics the bulk energy is elastic energy and the interfacial
energy is located on the unknown crack surfaces. For smectic A thin films we
will take the bulk energy to be the elastic (Oseen-Frank) energy and the
discontinuity surfaces will correspond to walls across which $\mathbf{n}$
jumps.

##  2 The model

The thin film experiments described above would best be treated as a free
boundary problem for a given volume of liquid crystal deposited on the
substrate, with the aim of predicting the possible 3D configurations of oily
streaks as well as their internal structure. A key prediction of such a model
would be the width of the oily streaks. However we are not currently able to
give conditions under which the minimum of the corresponding energy is
attained and if so by what configurations of layers.

Instead we will work on a fixed domain $\Omega\subset\mathbb{R}^{2}$, to be
thought of as the cross-section of an oily streak, and seek to predict the
corresponding configurations of layers given suitable boundary conditions.

### The function space.

In order to specify in precise mathematical terms any energy minimization
problem, it is necessary to say in which function space the minimum is sought.
The function space describes the allowed singularities of the unknown function
or map and is part of the model. Making the function space larger, so that
worse singularities are allowed, may lead to different minimizers (this is
known as the Lavrentiev phenomenon, see [5, Section 3]). The main function
space used for free discontinuity problems, and developed by De Giorgi and
Ambrosio [12], is the space SBV of special functions of bounded variation. In
fact, in this work we will consider a slight variant of the space
$\operatorname{SBV}$, i.e. the space $\operatorname{SBV}^{2}$. Let $\Omega$ be
an open, bounded planar region. In technical terms, a map
$u\colon\Omega\to\mathbb{R}^{m}$ belongs to
$\operatorname{SBV}^{2}(\Omega,\,\mathbb{R}^{m})$ if its distributional
derivative $\mathrm{D}\mathbf{u}$ is a finite measure with no Cantor part, and
whose absolutely continuous part $\nabla\mathbf{u}$ is square integrable (see
e.g [12, 4] for more details). For the purposes of this paper, the key points
of the definition and theory are the following:

  1. 1.

for any ${\bf u}\in\operatorname{SBV}^{2}(\Omega,\,\mathbb{R}^{m})$ there is a
one-dimensional jump set $\mathrm{S}_{\bf u}$ consisting of the points
$\mathbf{x}\in\Omega$ at which ${\bf u}$ has a jump discontinuity,

  2. 2.

for (almost) any point $\mathbf{x}\in\mathrm{S}_{\bf u}$ there is a well-
defined unit normal
${\boldsymbol{\nu}}(\mathbf{x})={\boldsymbol{\nu}}_{\mathbf{u}}(\mathbf{x})$
to $\mathrm{S}_{\bf u}$ and well-defined limits ${\bf u}^{+}(\mathbf{x}),{\bf
u}^{-}(\mathbf{x})$ from either side of $\mathrm{S}_{\bf u}$. There may be an
exceptional set of points $\mathbf{x}\in\mathrm{S}_{\mathbf{u}}$ at which
${\boldsymbol{\nu}}(\mathbf{x})$, $\mathbf{u}^{+}(\mathbf{x})$ or
$\mathbf{u}^{-}(\mathbf{x})$ are not defined, but this must be a set of zero
length;

  3. 3.

the gradient $\nabla{\bf u}$ is defined in $\Omega\setminus\mathrm{S}_{\bf u}$
and $\displaystyle\int_{\Omega}|\nabla{\bf
u}(\mathbf{x})|^{2}\,\mathrm{d}\mathbf{x}<+\infty$.

The space $\operatorname{SBV}^{2}(\Omega,\,\mathbb{R}^{2\times 2})$, whose
elements are $2\times 2$ matrix fields on $\Omega$ in the class
$\operatorname{SBV}^{2}$, is defined analogously by identifying
$\mathbb{R}^{2\times 2}$ with $\mathbb{R}^{4}$.

### The elastic energy.

We use as a basic variable the director $\mathbf{n}$, $|\mathbf{n}|=1$, which
for smectic A is parallel to the layer normal $\mathbf{m}$, so that
$\mathbf{n}\otimes\mathbf{n}=\mathbf{m}\otimes\mathbf{m}$. The Oseen-Frank
energy density

(2.1) $\displaystyle
W(\mathbf{n},\nabla\mathbf{n})=\frac{1}{2}\left(K_{1}({\rm
div}\,\mathbf{n})^{2}+K_{2}(\mathbf{n}\cdot\operatorname{curl}\mathbf{n})^{2}+K_{3}|\mathbf{n}\wedge\operatorname{curl}\mathbf{n}|^{2}\right.$
$\displaystyle\left.+(K_{2}+K_{4})(\operatorname{tr}(\nabla\mathbf{n})^{2}-({\rm
div}\,\mathbf{n})^{2})\right),$

with Frank constants $K_{i}$, is invariant to changing the sign of
$\mathbf{n}$, and thus can be expressed in terms of
$\mathbf{M}:=\mathbf{n}\otimes\mathbf{n}$ and $\nabla\mathbf{M}$. For example,
we have that

(2.2) $\displaystyle|\mathbf{n}\wedge\operatorname{curl}\mathbf{n}|^{2}$
$\displaystyle=$ $\displaystyle|(\mathbf{n}\cdot\nabla)\mathbf{n}|^{2}$
$\displaystyle=$ $\displaystyle\frac{1}{2}M_{rs}M_{ij,r}M_{ij,s},$

where $M_{ij}=n_{i}n_{j}$ and we have used
$2n_{i}n_{i,j}=(|\mathbf{n}|^{2})_{,j}=0$. We make the assumption of constant
layer thickness (or locally parallel layers), which for sufficiently smooth
$\mathbf{m}$ is equivalent to the condition $\operatorname{curl}\mathbf{m}=0$,
and thus (assuming $\mathbf{n}$ sufficiently smooth so that we can take
$\mathbf{n}=\mathbf{m}$) to the condition $\operatorname{curl}\mathbf{n}=0$.

We assume two-dimensional symmetry for $\mathbf{n}$, so that for
$\mathbf{x}=(x_{1},x_{2})\in\Omega$

(2.3) $\mathbf{n}(\mathbf{x})=({\bf n}(x_{1},x_{2}),0)$

with ${\bf n}(x_{1},x_{2})\in\SS^{1}$, where $\SS^{1}$ is the unit circle. For
such two-dimensional director fields, an explicit computation shows that the
term $\operatorname{tr}(\nabla{\bf n})^{2}-(\operatorname{div}{\bf n})^{2}$ is
identically equal to zero so that, taking into account also the constraint
$\operatorname{curl}\mathbf{n}=0$, the Oseen-Frank energy (2.1) reduces to

(2.4) $W({\bf n},\nabla{\bf n})=\frac{1}{2}K_{1}|\nabla{\bf n}|^{2},$

where we assume that $K_{1}>0$. (This is the original model for smectics
proposed by Oseen [20].) It can be viewed as a special case
$\mathbf{m}=\mathbf{n}$ of the model of [16]. There are various other models
for smectics allowing variable layer thickness ([11, 9, 15, 14]), that
typically introduce the molecular number density $\rho=\rho(\mathbf{x})$ or
fluctuations about it as a new macroscopic variable, with the smectic layers
being seen as density waves. These models can describe dislocations in smectic
layers, for example, although it is unclear how to understand the macroscopic
variable $\rho$ varying over a molecular length-scale.

To allow for nonorientable configurations, we reformulate the problem in terms
of the two-dimensional $\mathbf{Q}$-tensor

(2.5) $\mathbf{Q}:=\frac{1}{\sqrt{2}}\left({\bf n}\otimes{\bf n}-\frac{\bf
I}{2}\right)$

where $\mathbf{I}$ is the $2\times 2$ identity matrix. The constant
$1/\sqrt{2}$ is a normalization factor, which plays no essential role in our
analysis. At each point $\mathbf{x}\in\Omega$, $\mathbf{Q}(\mathbf{x})$ is a
symmetric trace-free $2\times 2$ matrix which belongs to the set

(2.6) $\mathscr{N}:=\left\\{\frac{1}{\sqrt{2}}\left({\bf n}\otimes{\bf
n}-\frac{\bf I}{2}\right)\colon{\bf n}\in\SS^{1}\right\\}.$

If $\mathbf{n}\colon\Omega\to\SS^{1}$ is a smooth vector field and
$\mathbf{Q}$ is defined as in (2.5), an explicit computation shows that
$\left|\nabla{\bf
n}\right|^{2}=\left|\nabla\mathbf{Q}\right|^{2}:=Q_{ij,k}\,Q_{ij,k}$,
expressing the elastic energy (2.4) in terms of $\mathbf{Q}$.

We also need to express the constraint $\operatorname{curl}\mathbf{n}=0$ in
terms of $\mathbf{Q}$. We define

${\bf
A}(\mathbf{Q})(\nabla\mathbf{Q},\nabla\mathbf{Q}):=\left(\sqrt{2}Q_{hk}+\frac{1}{2}\delta_{hk}\right)Q_{ij,h}\,Q_{ij,k}$

This is a quadratic form in $\nabla\mathbf{Q}$, reminiscent of the cubic term
in the Landau-de Gennes elastic energy (see e.g. [5, Section 4]). Thanks to
(2.2), (2.5), we find that

$\begin{split}{\bf
A}(\mathbf{Q})(\nabla\mathbf{Q},\nabla\mathbf{Q})=|\mathbf{n}\wedge\operatorname{curl}\mathbf{n}|^{2},\end{split}$

and since
$\left|\operatorname{curl}\mathbf{n}\right|^{2}=\left|\mathbf{n}\wedge\operatorname{curl}\mathbf{n}\right|^{2}+(\mathbf{n}\cdot\operatorname{curl}\mathbf{n})^{2}$,
and keeping in mind that a two-dimensional vector field is orthogonal to its
curl, we obtain
$\mathbf{A}(\mathbf{Q})(\nabla\mathbf{Q},\,\nabla\mathbf{Q})=\left|\operatorname{curl}\mathbf{n}\right|^{2}$.
Therefore, we impose the constraint

(2.7) $\mathbf{A}(\mathbf{Q})(\nabla\mathbf{Q},\,\nabla\mathbf{Q})=0,$

which expresses constant layer thickness in terms of $\mathbf{Q}$.

As we are considering a free-discontinuity problem, the map $\mathbf{Q}$ is
allowed to have jump discontinuities on a one-dimensional set
$\mathrm{S}_{\mathbf{Q}}$, and the condition (2.7) loses its meaning at points
of $\mathrm{S}_{\mathbf{Q}}$. Therefore, the constraint (2.7) is only enforced
in the complement $\Omega\setminus\mathrm{S}_{\mathbf{Q}}$. In technical
terms, (2.7) only involves the absolutely continuous part $\nabla\mathbf{Q}$
of the distributional derivative of $\mathbf{Q}$.

### The jump energy.

Let

(2.8)
$\operatorname{SBV}^{2}(\Omega,\,\mathscr{N}):=\left\\{\mathbf{Q}\in\operatorname{SBV}^{2}(\Omega,\,\mathbb{R}^{2\times
2}),\ \mathbf{Q}(x)\in\mathscr{N}\ \textrm{ for }x\in\Omega\right\\}$

be the set of $\mathscr{N}$-valued $\mathbf{Q}$-tensors in the class
$\operatorname{SBV}^{2}$. For
$\mathbf{Q}\in\operatorname{SBV}^{2}(\Omega,\,\mathscr{N})$, we consider the
jump energy

(2.9)
$\int_{\mathrm{S}_{\mathbf{Q}}}\varphi(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}})\,\mathrm{d}\mathscr{H}^{1}$

where the integral is taken over the (one-dimensional) jump set
$\mathrm{S}_{\mathbf{Q}}$, with respect to the length measure
$\mathrm{d}\mathscr{H}^{1}$, and
$\varphi\colon\mathscr{N}\times\mathscr{N}\times\SS^{1}\to\mathbb{R}$ is a
continuous function. We want to choose the jump energy density $\varphi$ so
that (2.9) is a good model for the energy of a domain wall. Natural conditions
on $\varphi$ are:

  1. (C1)

$\varphi(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}})=0$ if
$\mathbf{Q}^{+}=\mathbf{Q}^{-}$;

  2. (C2)

invariance with respect to the orientation of the jump set, that is

$\varphi(\mathbf{Q}^{-},\,\mathbf{Q}^{+},\,-{\boldsymbol{\nu}})=\varphi(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}})$

for any $\mathbf{Q}^{+}\in\mathscr{N}$, $\mathbf{Q}^{-}\in\mathscr{N}$,
${\boldsymbol{\nu}}\in\SS^{1}$;

  3. (C3)

frame-indifference, that is

$\varphi(\mathbf{R}\mathbf{Q}^{+}\mathbf{R}^{\mathsf{T}},\,\mathbf{R}\mathbf{Q}^{-}\mathbf{R}^{\mathsf{T}},\mathbf{R}{\boldsymbol{\nu}})=\varphi(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}})$

for all $\mathbf{Q}^{+}\in\mathscr{N}$, $\mathbf{Q}^{-}\in\mathscr{N}$,
${\boldsymbol{\nu}}\in\SS^{1}$ and all rotations
$\mathbf{R}\in\mathrm{SO}(2)$.

Another condition we impose is that $\varphi$ should penalize dislocations of
the smectic layers. Generically, it may not be possible to define the smectic
layers consistently across the jump set, because there may be dislocations. A
geometric condition for the layers to match at the jump set is that, at each
jump point, the normal to the jump set ${\boldsymbol{\nu}}$ bisects any of the
angles between the smectic layers. For smectic A liquid crystals, the layers
are orthogonal to the the molecular directors $\mathbf{n}^{+}$,
$\mathbf{n}^{-}$ on either side of the jump, so bisecting any of the angles
between the layers is equivalent to bisecting any of the angles between
$\mathbf{n}^{+}$ and $\mathbf{n}^{-}$. This can be written as

(2.10)
$(\mathbf{n}^{+}\cdot{\boldsymbol{\nu}})^{2}=(\mathbf{n}^{-}\cdot{\boldsymbol{\nu}})^{2}$

or equivalently, in terms of the $\mathbf{Q}$-tensor, as

(2.11)
$\mathbf{Q}^{+}{\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}=\mathbf{Q}^{-}{\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}.$

In order to penalize dislocations, we therefore impose the following condition
on $\varphi$:

  1. (C4)

for given $\mathbf{Q}^{+},\mathbf{Q}^{-}\in\mathscr{N}$, the function
${\boldsymbol{\nu}}\in\SS^{1}\mapsto\varphi(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}})$
is minimised for a ${\boldsymbol{\nu}}$ satisfying (2.11).

Given distinct $\mathbf{Q}^{+}$ and $\mathbf{Q}^{-}$, there are four such
${\boldsymbol{\nu}}$ that bisect the angles between the corresponding
directors $\mathbf{n}^{+}$, $\mathbf{n}^{-}$; we label them as
${\boldsymbol{\nu}}_{1}$, ${\boldsymbol{\nu}}_{2}$,
${\boldsymbol{\nu}}_{3}=-{\boldsymbol{\nu}}_{1}$ and
${\boldsymbol{\nu}}_{4}=-{\boldsymbol{\nu}}_{2}$, as in Fig. 3. By taking
$\mathbf{R}=-\mathbf{I}$ in condition (C3), we obtain that
$\varphi(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}}_{1})=\varphi(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}}_{3})$
and
$\varphi(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}}_{2})=\varphi(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}}_{4})$.
However, we must account for the possibility that, in general,
$\varphi(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}}_{1})\neq\varphi(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}}_{3})$.
Indeed, configurations whose jump sets are oriented according to the unit
normal ${\boldsymbol{\nu}}_{1}$ or ${\boldsymbol{\nu}}_{3}$ have (generically)
different geometric properties, because
$\mathbf{Q}^{+}{\boldsymbol{\nu}}_{1}\cdot{\boldsymbol{\nu}}_{1}\neq\mathbf{Q}^{+}{\boldsymbol{\nu}}_{3}\cdot{\boldsymbol{\nu}}_{3}$
unless the directors $\mathbf{n}^{+}$, $\mathbf{n}^{-}$ corresponding to
$\mathbf{Q}^{+}$, $\mathbf{Q}^{-}$ are orthogonal to each other. The condition
(C4) is compatible with an energy density that, as a function of
${\boldsymbol{\nu}}$, has two global minima at ${\boldsymbol{\nu}}_{1}$,
${\boldsymbol{\nu}}_{3}$ (or, respectively, at ${\boldsymbol{\nu}}_{2}$,
${\boldsymbol{\nu}}_{4}$) and, say, two local minima at
${\boldsymbol{\nu}}_{2}$, ${\boldsymbol{\nu}}_{4}$ (respectively, at
${\boldsymbol{\nu}}_{1}$, ${\boldsymbol{\nu}}_{3}$), at a possibly higher
energy value.

Figure 3: Left: given $\mathbf{Q}^{+}$ and $\mathbf{Q}^{-}$, there are four
unit vectors ${\boldsymbol{\nu}}_{1}$, ${\boldsymbol{\nu}}_{2}$,
${\boldsymbol{\nu}}_{3}$, ${\boldsymbol{\nu}}_{4}$ that satisfy (2.10). Centre
and right: two piecewise constant configurations that satisfy (2.10). The
black lines represent the smectic layers, which are orthogonal to the
molecular director, and the thick red line represents the jump.

A singular jump energy satisfying all the conditions (C1)–(C4) is given by

(2.12)
$\zeta_{\alpha}(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}}):=\begin{cases}\left|\mathbf{Q}^{+}-\mathbf{Q}^{-}\right|^{\alpha}&\textrm{if
}\mathbf{Q}^{+}{\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}=\mathbf{Q}^{-}{\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}\\\
+\infty&\textrm{otherwise,}\end{cases}$

where $\alpha$ is a parameter such that $0<\alpha<1$. Choosing $\alpha>0$
guarantees that the condition (C1) is satisfied. On the other hand, taking
$\alpha<1$ is required for reasons of mathematical consistency; had we taken
$\alpha\geq 1$, then there would be no guarantee that a minimizer for the
energy functional exists in the space
$\operatorname{SBV}^{2}(\Omega,\,\mathscr{N})$. (We do not discuss this issue
in detail and refer to e.g. [4, Sections 4.1–2].) However, even if the
parameter $\alpha$ is chosen carefully, the functional associated with (2.12)
suffers from a mathematical pathology, which is illustrated in Fig. 4: there
exist sequences $\mathbf{Q}_{j}\in\operatorname{SBV}(\Omega,\,\mathscr{N})$
that converge to a limit map
$\mathbf{Q}\in\operatorname{SBV}(\Omega,\,\mathscr{N})$ in a suitable sense,
yet the energy of $\mathbf{Q}$ is strictly larger than the limit, as
$j\to+\infty$, of the energy of $\mathbf{Q}_{j}$. This behaviour, known as
lack of lower semicontinuity, is particularly evident in Fig. 4, because the
energy of the limit configuration is infinite. However, even if we considered
a modified energy density that takes only finite values, this pathological
behaviour could persist. The main issue is that the energy density
$\zeta_{\alpha}$ is not _BV-elliptic_. BV-ellipticity, which was introduced by
Ambrosio and Braides [3], is a necessary condition for the lower
semicontinuity of the energy functional, and is an important assumption to
ensure the existence of minimizers of free-discontinuity problems [3, 1].

Figure 4: A mathematical pathology associated with the energy density
$\zeta_{\alpha}$, given by (2.12). The black lines represent the smectic
layers, while the thick red line represents the jump set. On the left, a
piecewise constant configuration in a rectangle, with a zig-zag interface and
zig-zag angles of $45$ degrees; on the right, a piecewise constant
configuration with a horizontal jump set. If we choose the energy density as
in (2.12), the energy of the configuration on the left is $\sqrt{2}b\mu$,
where $b$ is the length of the bottom base of the rectangle, irrespective of
the spacing of the zig-zag. As the latter tends to zero, the configuration on
the left converges (in a suitable sense) to that on the right, which costs an
infinite energy.

We briefly explain the notion of BV-ellipticity. Let
$\mathbf{Q}^{+}\in\mathscr{N}$, $\mathbf{Q}^{-}\in\mathscr{N}$,
${\boldsymbol{\nu}}\in\SS^{1}$ be given. Let $C$ be a (closed) unit square,
centred at the origin, whose sides are parallel to ${\boldsymbol{\nu}}$ and to
${\boldsymbol{\nu}}^{\perp}:=(-\nu_{2},\,\nu_{1})$. Let $C^{\prime}$ be a
slightly larger such square, such that $C$ is contained in the interior of
$C^{\prime}$. We consider a piecewise constant map
$\mathbf{Q}^{*}\in\operatorname{SBV}(C^{\prime},\,\mathscr{N})$, defined as

(2.13) $\mathbf{Q}^{*}(\mathbf{x}):=\begin{cases}\mathbf{Q}^{+}&\textrm{if
}\mathbf{x}\cdot{\boldsymbol{\nu}}>0\\\ \mathbf{Q}^{-}&\textrm{if
}\mathbf{x}\cdot{\boldsymbol{\nu}}<0\end{cases}$

The map $\mathbf{Q}^{*}$ jumps along a straight line orthogonal to
${\boldsymbol{\nu}}$. The notion of BV-ellipticity is defined in terms of a
suitable minimization problem. Let
$\mathscr{C}[\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}}]$ be the
class of all maps
$\mathbf{P}\in\operatorname{SBV}^{2}(C^{\prime},\,\mathscr{N})$ that satisfy
the following properties:

  1. (i)

$\mathbf{P}=\mathbf{Q}^{*}$ in $C^{\prime}\setminus C$;

  2. (ii)

$\mathbf{P}$ has finite range (i.e., $\mathbf{P}$ is a piecewise constant
configuration that only takes finitely many values).

The condition (i) plays the rôle of a boundary condition. As for condition
(ii), restricting our attention to piecewise constant configurations with
finite range allows us to neglect all elastic contributions for the moment and
focus on the jump energy.

######  Definition 2.1 (BV-ellipticity).

We say that a function
$\varphi\colon\mathscr{N}\times\mathscr{N}\times\SS^{1}\to\mathbb{R}$ is BV-
elliptic if and only if, for any
$(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}})\in\mathscr{N}\times\mathscr{N}\times\SS^{1}$,
we have

$\inf_{\mathbf{P}\in\mathscr{C}[\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}}]}\int_{\mathrm{S}_{\mathbf{P}}\cap
C}\varphi(\mathbf{P}^{+},\,\mathbf{P}^{-},\,{\boldsymbol{\nu}}_{\mathbf{P}})\,\mathrm{d}\mathscr{H}^{1}=\varphi(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}})$

Equivalently, $\varphi$ is BV-elliptic if and only if, for any
$(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}})\in\mathscr{N}\times\mathscr{N}\times\SS^{1}$,
the map $\mathbf{Q}^{*}$ defined in (2.13) minimizes the jump energy
functional associated with $\varphi$ among all competitors in
$\mathscr{C}[\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}}]$. Roughly
speaking, a jump energy density $\varphi$ is BV-elliptic if the corresponding
jump energy functional favours the jump set to be (locally) a straight line.
This indeed rules out pathological behaviour, such as the one discussed in
Fig. 4 (see e.g. [3, 4]). In practice, deciding whether a given function is
BV-elliptic or not may not an easy task. Sufficient and necessary conditions
for BV-ellipticity have been proposed in the literature (see e.g. [3, 2, 4,
8]), but even these might not be immediately applicable to concrete examples.

The function $\zeta_{\alpha}$ is not BV-elliptic, because it fails to satisfy
some necessary conditions for BV-ellipticity (such as convexity in the
variable ${\boldsymbol{\nu}}$ — see e.g. [4, Theorems 5.11 and 5.14]).
However, a natural candidate for a BV-elliptic function that satisfies
(C1)–(C4) is the _BV-elliptic envelope_ of $\zeta_{\alpha}$, that is, the
largest BV-elliptic function $\varphi$ such that $\varphi\leq\zeta_{\alpha}$.
As it turns out, the BV-elliptic envelope of $\zeta_{\alpha}$ can be
calculated explicitly, and is given by

(2.14)
$\varphi_{\alpha}(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}}):=\left|\mathbf{Q}^{+}-\mathbf{Q}^{-}\right|^{\alpha}\left(1+\frac{\sqrt{2}\left|\mathbf{Q}^{+}{\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}-\mathbf{Q}^{-}{\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}\right|}{\left|\mathbf{Q}^{+}-\mathbf{Q}^{-}\right|}\right)^{1/2}$

if $\mathbf{Q}^{+}\neq\mathbf{Q}^{-}$, and
$\varphi_{\alpha}(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}})=0$ if
$\mathbf{Q}^{+}=\mathbf{Q}^{-}$. The function $\varphi_{\alpha}$ can be
expressed in terms of angular variables, which is sometimes convenient for
computations. If $\mathbf{n}^{+}=(\cos\beta_{+},\,\sin\beta_{-})$,
$\mathbf{n}^{-}=(\cos\beta_{-},\,\sin\beta_{-})$ are the directors
corresponding to $\mathbf{Q}^{+}$, $\mathbf{Q}^{-}$ respectively and
${\boldsymbol{\nu}}=(\cos\gamma,\,\sin\gamma)$, with $\beta_{+}$, $\beta_{-}$,
$\gamma$ real numbers, then

(2.15)
$\begin{split}\varphi_{\alpha}(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}})&=\left|\sin(\beta_{+}-\beta_{-})\right|^{\alpha}\left(1+\left|\sin\left(\beta_{+}+\beta_{-}-2\gamma\right)\right|\right)^{1/2}\\\
&=\left|\sin(\beta_{+}-\beta_{-})\right|^{\alpha}\left(\left|\cos\left(\frac{\beta_{+}+\beta_{-}}{2}-\gamma\right)\right|+\left|\sin\left(\frac{\beta_{+}+\beta_{-}}{2}-\gamma\right)\right|\right)\end{split}$

Equation (2.15) is obtained from (2.14), by substituting (2.5) for
$\mathbf{Q}^{+}$, $\mathbf{Q}^{-}$ and applying trigonometric identities.

The proof that $\varphi_{\alpha}$ is indeed BV-elliptic (and, in fact, is the
BV-elliptic envelope of $\zeta_{\alpha}$) is rather technical, and we will
present it in a forthcoming paper [7]444 In [7], we will show that
$\varphi_{\alpha}$ is not only BV-elliptic, but also jointly convex, in the
sense of [4, Definition 5.17]. Joint convexity is a sufficient condition for
BV-ellipticity, but it is not known whether it is a necessary condition as
well.. However, from (2.14) and (2.15), it is immediate to check that
$\varphi_{\alpha}$ satisfies (C1)–(C4). For given values
$\mathbf{Q}^{+}\neq\mathbf{Q}^{-}$, the function
${\boldsymbol{\nu}}\mapsto\varphi_{\alpha}(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}})$
has four global minima, at the same energy value, corresponding exactly to the
four unit vectors ${\boldsymbol{\nu}}$ that satisfy (2.11) (i.e.
$\gamma=(\beta_{+}+\beta_{-})/2+k\pi/2$, for integer $k$). We do not have an
example of an energy density that satisfies (C1)–(C4), is BV-elliptic, and has
non-equal local minima at the four vectors ${\boldsymbol{\nu}}$ that satisfy
(2.11).

##  3 Existence of minimizers and examples

When the jump energy density is chosen as in (2.14), we can prove existence of
minimizers for the energy functional, subject to appropriate boundary
conditions. For instance, strong anchoring at the boundary is often
represented mathematically by imposing Dirichlet boundary conditions, of the
form $\mathbf{Q}=\mathbf{Q}_{{\mathrm{bd}}}$ on $\partial\Omega$, where
$\mathbf{Q}_{{\mathrm{bd}}}\colon\partial\Omega\to\mathscr{N}$ is a boundary
datum. However, these Dirichlet boundary conditions are not suited to free-
discontinuity problems. Instead, one has to define the boundary conditions
carefully in order to cater for the possibility that $\mathbf{Q}$ jumps at the
boundary. This is done by ‘thickening the boundary’. In other words, we
consider an exterior neighbourhood $\Gamma$ of the boundary $\partial\Omega$
and impose the condition

(3.1) $\mathbf{Q}=\mathbf{Q}_{\mathrm{bd}}\qquad\textrm{on }\Gamma,$

where the datum $\mathbf{Q}_{{\mathrm{bd}}}$ is now defined in $\Gamma$. The
map $\mathbf{Q}_{{\mathrm{bd}}}$ needs to be compatible with our setting
(i.e.,
$\mathbf{Q}_{{\mathrm{bd}}}\in\operatorname{SBV}^{2}(\Gamma,\,\mathscr{N})$
and it must satisfy
$\mathbf{A}(\mathbf{Q}_{{\mathrm{bd}}})(\nabla\mathbf{Q}_{{\mathrm{bd}}},\,\mathbf{Q}_{{\mathrm{bd}}})=0$
in $\Gamma$).

Using the direct methods in the Calculus of Variations, we could prove the
following existence theorem:

######  Theorem 3.1.

For any $K_{1}>0$, $\mu>0$ and $0<\alpha<1$, the energy functional:

(3.2)
$I(\mathbf{Q})=\frac{K_{1}}{2}\int_{\Omega}|\nabla\mathbf{Q}|^{2}\mathrm{d}x+\mu\int_{\mathrm{S}_{\mathbf{Q}}}\varphi_{\alpha}(\mathbf{Q}^{+},\mathbf{Q}^{-},{\bf\nu})\,\mathrm{d}\mathscr{H}^{1}$

with $\varphi_{\alpha}$ as in (2.14), attains a minimum among
$\mathbf{Q}\in\operatorname{SBV}^{2}(\Omega,\mathscr{N})$ satisfying ${\bf
A}(\mathbf{Q})(\nabla\mathbf{Q},\nabla\mathbf{Q})=0$ and the boundary
conditions (3.1).

The proof of Theorem 3.1 will be given in [7]. To try to get some
understanding of the behaviour of minimizers, we consider a simplified (or
over-simplified) problem, for which we can find the minimizer explicitly. We
consider a rectangular domain, $\Omega=(-L,\,L)\times(0,\,H)$, with $L>0$ and
$H>L/2$. We focus attention on a restricted class of configurations, whose
jump set can be described in polar coordinates as the curve

$\mathcal{C}=\left\\{(\rho(\theta)\cos\theta,\,\rho(\theta)\sin\theta)\colon
0\leq\theta\leq\pi\right\\}$

where $\rho\colon[0,\,\pi]\to\mathbb{R}$ is a scalar function to be
determined, subject to the conditions $\rho(0)=\rho(\pi)=L$. Moreover, we
assume that $\mathbf{Q}$ is given by

(3.3)
$\mathbf{Q}(\mathbf{x}):=\frac{1}{\sqrt{2}}\left(\mathbf{n}(\mathbf{x})\otimes\mathbf{n}(\mathbf{x})-\frac{\mathbf{I}}{2}\right)$

where the director $\mathbf{n}$ is defined in polar coordinates as

$\mathbf{n}(r\cos\theta,\,r\sin\theta):=\begin{cases}(\cos\theta,\,\sin\theta)&\textrm{if
}r<\rho(\theta)\\\ (0,\,1)&\textrm{otherwise.}\end{cases}$

This configuration $\mathbf{Q}$ is uniquely determined by the function $\rho$,
$\mathbf{Q}=\mathbf{Q}[\rho]$; it has horizontal layers above the curve
$\mathcal{C}$ and concentric circular layers below $\mathcal{C}$. In
particular, $\mathbf{Q}$ satisfies planar anchoring conditions at the boundary
$x_{2}=0$, homeotropic anchoring conditions at $x_{2}=H$ and periodic boundary
conditions at $x_{1}=\pm L$. To simplify the problem further, we neglect the
elastic energy. (Heuristically, we expect this approximation to be relevant in
the limit as $\mu/K_{1}\to+\infty$.)

Figure 5: The minimizing configuration given by Proposition 3.2. The black
lines represent the smectic layers, while the thick red line is the jump set
$\mathcal{C}$.

######  Proposition 3.2.

For any $0<\alpha<1$, the unique minimizer of the functional

$\rho\mapsto\int_{\mathcal{C}}\varphi_{\alpha}(\mathbf{Q}^{+}[\rho],\,\mathbf{Q}^{-}[\rho],\,{\boldsymbol{\nu}}_{\mathbf{Q}[\rho]})\,\mathrm{d}\mathscr{H}^{1},$

among all $($Lipschitz continuous$)$ functions
$\rho\colon[0,\,\pi]\to\mathbb{R}$ such that $\rho(0)=\rho(\pi)=L$, is given
by

(3.4) $\rho(\theta):=\frac{L}{1+\sin\theta}$

for $\theta\in[0,\,\pi]$.

The proof of Proposition 3.2 will be given in [7]. When $\rho$ is given by
(3.4), the jump set $\mathcal{C}$ can be described in Cartesian coordinates as
the graph of

$x_{2}=\frac{L}{2}-\frac{x_{1}^{2}}{2L}\qquad\textrm{for }-L\leq x_{1}\leq L.$

In particular, $\mathcal{C}$ is a parabolic arc (see Fig. 5). It can be shown
that, at each point, the normal to the curve $\mathcal{C}$ bisects the angle
between the smectic layers.

Next we consider a domain which is a quarter disk of unit radius,
$\Omega:=\\{\mathbf{x}=(x_{1},\,x_{2})\in\mathbb{R}^{2}\colon
x_{2}^{2}+x_{2}^{2}<1,\ x_{1}>0,\ x_{2}>0\\}$. Again, we minimize within a
restricted class of configurations, whose jump set has the form

(3.5)
$\mathcal{C}=\left\\{(\rho(\theta)\cos\theta,\,\rho(\theta)\sin\theta)\colon
0\leq\theta\leq\frac{\pi}{2}\right\\}$

The function $\rho$ is unknown but, in contrast to the previous example, the
boundary values $\rho(0)$, $\rho(\pi/2)$ are unspecified. We consider maps
$\mathbf{Q}$ of the form (3.3), where the director $\mathbf{n}$ is given by

(3.6) $\mathbf{n}(r\cos\theta,\,r\sin\theta):=\begin{cases}(0,\,1)&\textrm{if
}r<\rho(\theta),\ 0<\theta<\pi/2\\\ (\cos\theta,\,\sin\theta)&\textrm{if
}r>\rho(\theta),\ 0<\theta<\pi/2.\end{cases}$

In particular, the layers are horizontal on the left side of $\mathcal{C}$ and
circular on the right side. Moreover, we impose planar anchoring conditions
near the boundary at $x_{1}=0$, defined in terms of the director as
$\mathbf{n}(x_{1},\,x_{2})=(1,\,0)$ for $x_{2}<0$. As a result, the jump set
of $\mathbf{Q}$ contains an additional component on the boundary of
$\partial\Omega$, i.e. the straight line segment of endpoints $(0,\,0)$ and
$(0,\,\rho(0))$, where the director jumps from the value $(1,\,0)$ to $(0,1)$.
In the energy functional, we account for this additional contribution as well.
Moreover, in this example we do _not_ neglect the elastic energy.

(a) $\mu=1$, $\alpha=0.5$ (b) $\mu=1$, $\alpha=0.2$ (c) $\mu=2$, $\alpha=0.5$
(d) $\mu=2$, $\alpha=0.2$

Figure 6: Numerics for the problem on a quarter circle. All the pictures have
$K_{1}=2$.

Fig. 6 shows a numerical approximation of the minimizer of (3.2) within this
restricted class, for $K_{1}=2$ and different values of $\alpha$, $\mu$. The
simulation is based on a MATLAB code and the details are provided in Appendix
6. The numerical method we use is _not_ guaranteed to converge to a global
minimum of the problem. However, we repeated the simulations using several
different initial guesses for $\rho$, including random ones, and obtained
qualitatively similar profiles for the jump set. The numerically found jump
set very nearly agrees with a union of two parabolic arcs that bisect the
smectic layers at each point, given explicitly as

(3.7)
$x_{1}=\min\left(\left(a^{2}+2ax_{2}\right)^{1/2},\,\left(b^{2}-2bx_{2}\right)^{1/2}\right)$

for some positive numbers $a$, $b$. In Fig. 6, we compare the numerical
solution with the parabolic arcs given by (3.7), where we have chosen the
parameters $a$, $b$ so as to match the boundary values of the numerically
found solution; the difference between the two curves is almost unnoticeable.
Numerical tests show that the jump set tends to shrink towards the centre of
the circle as $\mu$ increases or $\alpha$ decreses, although the dependence on
$\alpha$ seems to be weaker. This is consistent with the fact that the jump
energy density is monotonically increasing as a function of $\mu$ and
decreasing as a function of $\alpha$. Therefore, as $\mu$ increases or
$\alpha$ decreases, energy minimization requires the jump set to become
shorter, in order to compensate for the larger values of the jump energy
density.

##  4 Discussion and future directions

In this paper we have shown how a free-discontinuity model has the potential
to describe the configurations of smectic A layers observed in thin film
experiments. An advantage of free-discontinuity models in this context is that
they represent domain walls as sharp interfaces with precise locations.
However, our analysis is just a beginning and much remains to be done. In
particular, we are currently unable to formulate the thin film problem as a
well-posed free-boundary problem in either two or three dimensions, so that
the profile of the upper free surface and the width of the oily streaks can be
predicted. Analytic work on these problems would need to be supplemented by
numerical studies, and there is a need to develop appropriate numerical
methods for free-discontinuity (and free boundary) problems associated with
jump energy densities of the type (2.14). Such jump energies also need further
study, and it would be useful to broaden the class of BV-elliptic jump
energies satisfying (C1)–(C4), in particular to allow different minima at the
two angle bisectors.

Finally, it is important to understand the relation between models based on a
molecular density and our ‘sharp interface’ model. In this context there are
interesting numerical computations of Xia, Maclachlan, Atherton & Farrell [22]
using a modification of a model of Ball & Bedford [6], which is in turn based
on that of Pevnyi, Selinger & Sluckin [21].

##  5 Acknowledgement.

We would like to thank Emmanuelle Lacaze for many very helpful discussions. We
also thank Giacomo Albi and Marco Caliari from the University of Verona for
their help with the numerical simulations.

## References

  * [1] L. Ambrosio.  Existence theory for a new class of variational problems.  Arch. Rational Mech. Anal., 111(4):291–322, 1990.
  * [2] L. Ambrosio and A. Braides.  Functionals defined on partitions in sets of finite perimeter ii: semicontinuity, relaxation and homogenization.  J. Math. Pures Appl., 69:307–333, 1990.
  * [3] L. Ambrosio and A. Braides.  Functionals defined on partitions of sets of finite perimeter, I: integral representation and ${\Gamma}$-convergence.  J. Math. Pures. Appl., 69:285–305, 1990.
  * [4] L. Ambrosio, N. Fusco, and D. Pallara.  Functions of bounded variation and free discontinuity problems.  Oxford Mathematical Monographs. The Clarendon Press, Oxford University Press, New York, 2000.
  * [5] J. M. Ball.  Mathematics of liquid crystals.  Molecular Crystals and Liquid Crystals, 647(1):1–27, 2017.
  * [6] J. M. Ball and S. J. Bedford.  Discontinuous order parameters in liquid crystal theories.  Molecular Crystals and Liquid Crystals, 612(1):1–23, 2015.
  * [7] J. M. Ball, G. Canevari, and B. Stroffolini.  Analysis of a free discontinuity model for smectic thin films.  To appear.
  * [8] D. G. Caraballo.  BV-ellipticity and lower semicontinuity of surface energy of Caccioppoli partitions of $\mathbb{R}^{n}$.  Journal of Geometric Analysis, 23:202– 220, 2013.
  * [9] J.-H. Chen and T. C. Lubensky.  Landau-Ginzburg mean-field theory for the nematic to smectic-C and nematic to smectic-A phase transitions.  Physical Review A, 14(3):1202, 1976.
  * [10] D. Coursault, B. Zappone, A. Coati, A. Boulaoued, L. Pelliser, D. Limagne, N. Boudet, Bicher H. Ibrahim, A. De Martino, M. Alba, M. Goldmann, Y. Garreau, B. Gallasa, and E. Lacaze.  Self-organized arrays of dislocations in thin smectic liquid crystal films.  Soft Matter, 12(3):678–688, 2016.
  * [11] P. G. de Gennes.  An analogy between superconductors and smectics A.  Solid State Communications, 10(9):753–756, 1972.
  * [12] E. De Giorgi and L. Ambrosio.  Un nuovo funzionale nel calcolo delle variazioni.  Atti Accad. Naz. Lincei Rend. Cl. Sci. Fis. Mat. Nat. (8), 82(2):199–210 (1989), 1988.
  * [13] G. A. Francfort and J.-J. Marigo.  Revisiting brittle fracture as an energy minimization problem.  Journal of the Mechanics and Physics of Solids, 46(8):1319–1342, 1998.
  * [14] J. Han, Y. Luo, W. Wang, P. Zhang, and Z. Zhang.  From microscopic theory to macroscopic theory: a systematic study on modeling for liquid crystals.  Archive for Rational Mechanics and Analysis, 215(3):741–809, 2015\.
  * [15] M. Kléman and O. Parodi.  Covariant elasticity for smectics A.  Journal de Physique, 36(7-8):671–681, 1975.
  * [16] F. M. Leslie, I. W. Stewart, and M Nakagawa.  A continuum theory for smectic C liquid crystals.  Molecular Crystals and Liquid Crystals, 198(1):443–454, 1991.
  * [17] J.-P. Michel, E. Lacaze, M. Alba, M. De Boissieu, M. Gailhanou, and M. Goldmann.  Optical gratings formed in thin smectic films frustrated on a single crystalline substrate.  Physical Review E, 70(1):011709, 2004.
  * [18] J.-P. Michel, E. Lacaze, M. Goldmann, M. Gailhanou, M. de Boissieu, and M. Alba.  Structure of smectic defect cores: X-ray study of 8CB liquid crystal ultrathin films.  Phys. Rev. Lett., 96:027803, Jan 2006.
  * [19] J. D. Niyonzima et al.  Topological defects in smectic thin films.  Talk at ILCC 2022.
  * [20] C. W. Oseen.  The theory of liquid crystals.  Trans. Faraday Soc., 29:883–899, 1933.
  * [21] M. Y. Pevnyi, J. V. Selinger, and T. J. Sluckin.  Modeling smectic layers in confined geometries: Order parameter and defects.  Phys. Rev. E, 90:032507, Sep 2014.
  * [22] J. Xia, S. MacLachlan, T. J. Atherton, and P. E. Farrell.  Structural landscapes in geometrically frustrated smectics.  Physical Review letters, 126(17):177801, April 2021.
  * [23] B. Zappone, C. Meyer, L. Bruno, and E. Lacaze.  Periodic lattices of frustrated focal conic defect domains in smectic liquid crystal films.  Soft Matter, 8:4318–4326, 2012.
  * [24] B. Zappone and E. Lacaze.  Surface-frustrated periodic textures of smectic-A liquid crystals on crystalline surfaces.  Phys. Rev. E, 78:061704, Dec 2008.
  * [25] B. Zappone and E. Lacaze.  One-dimensional patterns and topological defects in smectic liquid crystal films.  Liquid Crystals Reviews, 0(0):1–18, 2022.

##  6 Appendix: numerical simulations

In this section, we present the details of the numerical simulations for the
minimization problem on a quarter circle, described in Section 3. The
$\mathbf{Q}$-tensor is defined in terms of a scalar function
$\rho\colon[0,\,\pi/2]\to\mathbb{R}$, as in (3.3), (3.6). However, we find it
convenient to introduce a new variable $u=u(\theta)$, defined by

(6.1) $\rho(\theta)=e^{-u(\theta)}\qquad\textrm{for any
}0<\theta<\frac{\pi}{2}.$

The original variable $\rho$ is subject to the geometric constraints
$0<\rho<1$, because the domain is a quarter circle of unit radius. However, we
minimize numerically the energy functional subject to _no_ constraints on $u$.
Equation (6.1) guarantees that $\rho>0$, but for some values of $\alpha$,
$\mu$ (not shown in Fig. 6), we did find numerical solutions that are not
admissible, because they do not satisfy $\rho<1$.

The energy can be written as a functional of $u$ by a direct computation,
using (3.3), (3.6) and (6.1). The energy consists of a sum of three terms:

(6.2)
$I(u)=I^{\mathrm{el}}(u)+I^{\mathrm{jump,int}}(u)+I^{\mathrm{jump,bd}}(u)$

The first term, $I^{\mathrm{el}}(u)$, is the elastic energy:

(6.3)
$I^{\mathrm{el}}(u):=\frac{K_{1}}{2}\int_{\Omega}\left|\nabla\mathbf{Q}(\mathbf{x})\right|^{2}\,\mathrm{d}\mathbf{x}=\frac{K_{1}}{2}\int_{0}^{\pi/2}\left(\int_{e^{-u(\theta)}}^{1}\frac{\mathrm{d}\rho}{\rho}\right)\mathrm{d}\theta=\frac{K_{1}}{2}\int_{0}^{\pi/2}u(\theta)\,\mathrm{d}\theta$

The second term, $I^{\mathrm{jump,int}}(u)$, is the contribution to the jump
energy from the curve $\mathcal{C}$, given in (3.5):

(6.4)
$\begin{split}I^{\mathrm{jump,int}}(u)&:=\mu\int_{\mathscr{C}}\varphi(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}})\,\mathrm{d}\mathscr{H}^{1}\\\
&=\mu\int_{0}^{\pi/2}\left|\cos\theta\right|^{\alpha}\,e^{-u(\theta)}\left({u^{\prime}(\theta)}^{2}+1+\left|f(\theta,\,u^{\prime}(\theta))\right|\right)^{1/2}\mathrm{d}\theta\end{split}$

with

(6.5)
$f(\theta,\,u^{\prime}(\theta)):=({u^{\prime}(\theta)}^{2}-1)\cos\theta+2u^{\prime}(\theta)\,\sin\theta$

The expression (6.4)–(6.5) is deduced from (2.15), by applying trigonometric
identities to simplify the form of the integrand. The integrand in (6.4) is
not differentiable at the points where $f(\theta,\,u^{\prime}(\theta))=0$. As
the numerical minimization is based on a quasi-Newton method, which works best
for differentiable functions, we regularize the integrand by introducing a
parameter $\varepsilon>0$:

(6.6)
$\begin{split}I^{\mathrm{jump,int}}(u)\approx\mu\int_{0}^{\pi/2}\left|\cos\theta\right|^{\alpha}\,e^{-u(\theta)}\left({u^{\prime}(\theta)}^{2}+1+\left(\varepsilon+f(\theta,\,u^{\prime}(\theta))^{2}\right)^{1/2}\right)^{1/2}\mathrm{d}\theta\end{split}$

The pictures in Fig. 6 are obtained for $\varepsilon=10^{-12}$. Finally,
$I^{\mathrm{jump,bd}}(u)$ is the contribution to the jump energy from jumps
that are located on the boundary of $\Omega$ — more precisely, on the line
segment $S$ of endpoints $(0,\,0)$ and $(0,\,\rho(0))$:

(6.7)
$\begin{split}I^{\mathrm{jump,bd}}(u)&:=\mu\int_{S}\varphi(\mathbf{Q}^{+},\,\mathbf{Q}^{-},\,{\boldsymbol{\nu}})\,\mathrm{d}\mathscr{H}^{1}=\sqrt{2}\mu\,e^{-u(0)}\end{split}$

This term can also be written in integral form, by considering a smooth
function $g\colon[0,\,\pi/2]\to\mathbb{R}$ such that $g(0)=1$, $g(\pi/2)=0$
and writing

(6.8)
$\begin{split}I^{\mathrm{jump,bd}}(u)&=-\sqrt{2}\mu\int_{0}^{\pi/2}\frac{\mathrm{d}}{\mathrm{d}\theta}\left(e^{-u(\theta)}\,g(\theta)\right)\mathrm{d}\theta=\sqrt{2}\mu\int_{0}^{\pi/2}e^{-u(\theta)}\left(u^{\prime}(\theta)\,g(\theta)-g^{\prime}(\theta)\right)\mathrm{d}\theta\end{split}$

We found that, when writing the boundary term in the form (6.7), the numerical
solution deviates from the parabolic arcs (3.7) near $\theta=0$. However, the
thickness of this ‘boundary layer’ is mesh-dependent, so this feature is
probably a numerical artifact. The pictures in Fig. 6 are obtained by
considering the integral form (6.8) of $I^{\mathrm{jump,bd}}$, with
$g(\theta):=1-2\theta/\pi$, and present no boundary layer. Other choices of
the function $g$, e.g. $g(\theta):=\cos\theta$, produce qualitatively similar
profiles for the jump set.

We discretize the functional (6.2) on a uniform mesh of $m$ points in
$[0,\,\pi/2]$. We approximate the derivative $u^{\prime}(\theta)$ by second-
order central finite differences and we approximate the integrals by the
trapezoid rule. The number $m$ of mesh points is increased gradually, from
$m=50$ to $m=100$. At each step of the iteration over $m$, we call the built-
in MATLAB function fminunc, which applies the Broyden-Fletcher-Goldfarb-Shanno
Quasi-Newton algorithm to minimize the functional (6.2) under no constraints
on $u$. The initial guess for the minimization process is defined by the
numerical minimizer found at the previous step.


Plasma Injection Schemes for Laser–Plasma Accelerators

J. Faure

Laboratoire d'Optique Appliquée, ENSTA-CNRS École Polytechnique, Palaiseau, France

Plasma injection schemes are crucial for producing high-quality electron beams in laser–plasma accelerators. This article introduces the general concepts of plasma injection. First, a Hamiltonian model for particle trapping and acceleration in plasma waves is introduced; ionization injection and colliding-pulse injection are described in the framework of this Hamiltonian model. We then proceed to consider injection in plasma density gradients.

Laser-plasma accelerators; laser wakefield; electron injection; ionization injection; colliding pulse injection; gradient injection.

§ INTRODUCTION

In accelerator physics, the starting point of a high-energy machine is always the injector. The injector can be extremely important; it usually delivers a particle beam at the$ \UMeV$  level, and its characteristics—such as emittance, bunch duration and energy spread—strongly impact the parameters of the final beam. By analogy, it seems only natural that laser–plasma accelerators should also include an injection stage; this would permit (i) decoupling of the injection and acceleration stages, (ii) more control over the beam parameters and the possibility to tune them independently, and (iii) better stability of the beam. However, injecting electrons into the plasma wakefield accelerating structure is not an easy task: the wavelength of the accelerating structure is on the order of $\lambda_{\rm p}\simeq10\:\mu$m; here, the plasma wavelength is $\lambda_{\rm p}=2\pi c/\omega_{\rm p}$ (where $c$ is the speed of light and $\omega_{\rm p}=(n_{\rm e} e^2/\epsilon_0m_{\rm e})^{1/2}$ is the plasma frequency), so that $\lambda_{\rm p}=10$–$30\:\mu$m for electron plasma densities in the $10^{18}$–$10^{19}\:$cm$^{-3} $ range. The production of monoenergetic electron bunches in such a micrometre structure requires that the injected beam have a duration shorter than $\lambda_{\rm p}/c=30$–$100\:$fs and that it be focused to a few microns. In addition to these stringent requirements, the injected particle beam should be synchronized at the femtosecond level with the laser pulse driving the wakefield, which also poses a considerable experimental challenge. The production of such short and well-synchronized bunches is at the edge of conventional radio-frequency (r.f.) accelerator technology, and experiments using external injectors have not yet succeeded in providing high-quality electron beams with narrow energy spreads [1, 2].

In many laser wakefield acceleration experiments, electrons are self-injected into the accelerating structure: when the plasma wave amplitude reaches very high levels (close to the wave-breaking threshold), background plasma electrons can be injected into the plasma waves [3, 4, 5, 6, 7]. This is somewhat analogous to dark current in a r.f. accelerator: at high field strengths, the cavities start to release electrons. This self-injection mechanism can lead to the production of high charge and, on occasion, narrow energy spreads in the range of a few percent to 10% can be achieved. However, it has proven difficult to obtain high-quality, i.e. narrow, energy spreads in a stable manner using this method. In addition, electron injection results from a succession of nonlinear effects, such as relativistic self-focusing [8, 9] of the laser pulse, spectral broadening and self-steepening [10, 11, 12, 13]. This is the reason that self-injection is rather difficult to control precisely and does not allow tuning of the injected beam parameters. In this context, many researchers in the field of laser–plasma accelerators have started to develop injection methods in which the injected electrons originate from the plasma itself (as opposed to external injection requiring an existing electron source).

While several methods have been proposed for injecting plasma electrons into an existing plasma wakefield, they can be summarized by the following general principles.

* Create electrons at the right phase in the wakefield. Even if an electron is created at rest, it can end up being trapped provided that it is born at the appropriate phase. This is the idea behind ionization injection.

* Give an initial kick to electrons so that their initial longitudinal velocity is high enough for trapping. This is analogous to surfers paddling in order to gain momentum and catch the wave. Several methods are based on this idea, including injection by colliding pulses.

* Slow down the wakefield to facilitate trapping. Decreasing the phase velocity of the wakefield locally can be achieved by tailoring the plasma density, for example, and it leads to a controlled injection mechanism. Injection in downward density ramps is based on this idea.

In this article we will review some of these injection mechanisms. In Section <ref>, a Hamiltonian formalism based on  [14, 15] is derived in order to find the trapping threshold for an electron interacting with an intense laser pulse and a plasma wave. In Section <ref>, the findings from the Hamiltonian model will be applied to the particular case of ionization injection where high-$Z$ atoms provide the source of injected electrons. In Section <ref> we show how electrons can be injected using colliding and counter-propagating laser pulses. Finally, in Section <ref>, we show how electrons can also be trapped in wakefields with decreasing phase velocities.

§ ONE-DIMENSIONAL HAMILTONIAN MODEL

§.§ Assumptions of the model

The goal of this section is to derive
analytically the trajectories of electrons in a laser field and a plasma wave [15]. In
particular, some electrons can be trapped and accelerated in the
plasma wakefield. The theoretical framework is the following.

* We start with a 1D model: we consider only the motion of electrons along the longitudinal coordinate $z$, and neglect the role of the radial electric fields. In
this case, the wakefield potential $\phi$ depends only on $z$ and $t$.

* For simplicity, we assume that the laser pulse does not
change during its propagation. Consequently,
the plasma wakefield is also stationary along the propagation.
This is important because it will allow us to use a
conservation-of-energy law.

* The plasma is modelled by an electron fluid. This fluid is
described by macroscopic quantities such as its density $n(\mathbf{r},t)$ and
velocity $\mathbf{v}(\mathbf{r},t)$. Let us note that in
such a model, kinetic effects (e.g. trapping, wave-breaking) are not
taken into account.

* Plasma ions are assumed to be immobile. This is
justified when the typical time for ion motion ($\omega_{\rm pi}^{-1}$) is large
compared to the driver pulse duration
(i.e. $\tau\ll\omega_{\rm pi}^{-1}$).

* The electron fluid is cold. In the case of a laser driver, this is justified
when the quiver velocity of electrons in the laser field is orders of magnitude greater
than the thermal velocity: $v_{\rm osc}\simeq~eE_{\rm laser}/(m_{\rm e}\omega_0)\gg
v_{\rm th}=(k_{\rm B} T_{\rm e}/m_{\rm e})^{1/2}$.

The laser pulse is represented by its normalized vector potential
$\mathbf{a}=e\mathbf{A}/m_{\rm e}c$, where $\mathbf{A}$ is the laser vector potential.
For a pulse propagating along the $z$-axis and polarized
along the $x$-axis,

\begin{equation}
\mathbf{a}=\hat{a}(z,t)\cos(k_0z-\omega_0t)\,\mathbf{e}_x,
\end{equation}

where $k_0$ and $\omega_0$ are the wavevector and frequency of
the laser electromagnetic field, and $\hat{a}$ is an envelope
function that represents the longitudinal shape of the pulse.
We assume a Gaussian shape for $\hat{a}$
so that $\hat{a}^2(\zeta)=a_0^2 \exp(-\zeta^2/L_0^2)$, where $\zeta=z-v_{\rm g}t$ with $v_{\rm g}$ being the laser group velocity, $L_0$ is the laser pulse length, and $a_0$
is the laser peak amplitude, $a_0=8.6\times10^{-10}\sqrt{I\,[\wcm]\,\lambda^2\,[\mathrm{\mu m}]}$, with $I$ and $\lambda$ being the laser intensity and the central wavelength, respectively. The plasma wakefield is represented by its normalized potential $\phi=e\Phi/m_{\rm e}c^2$, which is obtained by solving the equation

\begin{equation}\label{eqphi}
\frac{\partial^2\phi}{\partial
\zeta^2}=k_{\rm p}^2\gamma_{\rm p}^2\left[\beta_{\rm p}\Bigg(1-
\frac{1+a^2}{\gamma_{\rm p}^2(1+\phi)^2}\Bigg)^{-1/2} -1\right].
\end{equation}

Here, $\gamma_{\rm p}$ is the Lorentz factor corresponding to the plasma wave phase velocity $v_{\rm p}$: $\gamma_{\rm p}=(1-v_{\rm p}^2/c^2)^{-1/2}$. It is also assumed that the plasma wave velocity is given by the group velocity of the excitation laser: $v_{\rm p}\simeq v_{\rm g}$ and $\gamma_{\rm p}\simeq \gamma_{\rm g}$. Finally, the plasma wavevector is $k_{\rm p}=\omega_{\rm p}/v_{\rm p}$. eqphi is the nonlinear plasma wave equation, which is valid even for $a>1$, i.e. for relativistic laser intensities with $I>10^{18}\:\wcm$. Figure <ref> shows the result for a 20 fs laser pulse with $a_0=2$ and a plasma density of $n_{\rm e}=7\times~10^{18}\:\cc$. Note that the laser field was averaged over the fast time-scale at $\omega_0$ for clarity (keeping the fast frequency does not affect our conclusions). Notice the characteristic nonlinear shape of the plasma density perturbation, which in turn causes a nonlinear longitudinal electric field with a sawtooth shape.

The top panel plots the laser amplitude $a$ (red) and the wakefield potential $\phi$ (blue); the middle panel shows the corresponding normalized electric field $E_z/E_0$ (where $E_0=m_{\rm e}c\,\omega_{\rm p}/e$ is the wave-breaking field). The bottom panel shows the electron density perturbation $n(\zeta)/n_0-1$.

§.§ Hamiltonian for an electron interacting with a laser and a plasma wave

Let us consider an electron with Lorentz factor $\gamma$ in a 1D wakefield, represented by the normalized potential $\phi$. Its Hamiltonian reads $H =\gamma-\phi(z-v_{\rm g}t)=\sqrt{1+u_\bot^2+u_z^2}-\phi(z-v_{\rm g}t)$, where $u_{\bot}=p_{\bot}/(m_{\rm e}c)$ and $u_{z}=p_{z}/(m_{\rm e}c)$ are the transverse and longitudinal normalized linear momenta, respectively. The Hamiltonian depends on $z$ and $t$ in a particular manner: it depends only on the variable $\zeta=z-v_{\rm g}t$. We can use this
fact to change variables using a canonical transformation
$(z,u_z)\rightarrow (\zeta,u_z)$. We use a second-type generating
function $F_2(z,u_z)=u_z(z-v_{\rm g}t)$; the transformation satisfies $\zeta=\partial F_2/\partial
u_z=z-v_{\rm g}t$ and $u_z=\partial F_2/\partial z$, and the new Hamiltonian
is given by $H'=H+\frac{1}{c}\frac{\partial F_2}{\partial t}$. In this case, the new Hamiltonian (denoted simply by $H$ from now on) is

\begin{equation}\label{EqHbasic}
H=\sqrt{1+u_\bot^2+u_z^2}-\phi(\zeta)-\beta_{\rm p}u_z
\end{equation}

where $\beta_{\rm p}=v_{\rm p}/c$.

This Hamiltonian has several constants of motion. We now introduce the canonical momentum $\mathbf{P}=\mathbf{p}+q\mathbf{A}$. In normalized units, the canonical momentum is denoted by $\mathbf{U}$ and we see that $U_z=u_z$ and $U_\bot=u_\bot-a$ ($q=-e$ for electrons), because in 1D the laser field has only a transverse component. Writing the Hamiltonian with just the transverse canonical momentum gives

\[
H =\sqrt{1+(U_\bot+a)^2+u_z^2} -\phi(\zeta)-\beta_{\rm p} u_z.
\]

From Hamilton's equations, one finds that in 1D, the transverse canonical momentum is conserved:

\begin{equation*}
\dot{U}_\bot=-\frac{\partial H}{\partial r_\bot}=0 \,\implies \, u_\bot(\zeta)-a(\zeta)=\mathrm{const}.
\end{equation*}

In the case of an electron initially at rest in front of the laser pulse, $\zeta_{\rm i}=+\infty$, this gives the important result that $u_\bot(\zeta)=a(\zeta)$. The other constant of motion is simply the energy: this Hamiltonian does not depend on time but only on
$\zeta$; as a consequence, the energy of the system is conserved and
$H$ is conserved along an electron trajectory. Hence, for an electron with initial energy $H_0$, one can solve for its longitudinal momentum $u_z$ to get

\begin{equation}\label{EqTraj}
u_z=\beta_{\rm p}\gamma_{\rm p}^2(H_0+\phi)\pm\gamma_{\rm p}\sqrt{\gamma_{\rm p}^2(H_0+\phi)^2-\gamma_\bot^2}
\end{equation}

with $\gamma_\bot^2=1+u_\bot^2$. Once $a(\zeta)$ and $\phi(\zeta)$ are known, this equation gives the electron trajectory in $(\zeta,u_z)$ phase space. Figure <ref> shows various electron trajectories that were obtained for different initial conditions.
For electrons initially at rest in front of the laser pulse, i.e. $\zeta_{\rm i}=+\infty$ and $u_z(\zeta_{\rm i})=u_\bot(\zeta_{\rm i})=~0$, the Hamiltonian is $H_0=1$. The trajectory of such electrons is referred to as the fluid orbit and corresponds to the trajectory of plasma background electrons that contribute to the formation of the plasma wakefield (black line in lower panel of figsep). These electrons are not trapped and oscillate in the plasma wakefield with low energies.
The separatrix is the special trajectory which separates the trapped from the untrapped orbits. It can be found by considering an electron moving along $z$ with $v_z=v_{\rm p}$ (or $u_z(\zeta_{\min})=\beta_{\rm p}\gamma_{\rm p}$) and located initially at a minimum of the potential $\phi(\zeta_{\min})=\phi_{\min}<0$, i.e. at a node of the electric field. Conservation of canonical momentum reads $u_\bot(\zeta_{\min})=a(\zeta_{\min})$, so that the Hamiltonian on the separatrix is

\begin{equation}
H_{\rm sep}=\frac{\sqrt{1+a^2(\zeta_{\min})}}{\gamma_{\rm p}}-\phi_{\min}.
\end{equation}

This trajectory is represented by a red line in the lower panel of figsep. The dashed blue lines represent some typical trapped orbits for which electrons gain significant energy in the plasma wave. Trapped orbits are obtained when the Hamiltonian satisfies $H_{\rm i}\leqslant H_{\rm sep}$. To prove this, let us consider a trapped electron behind the laser pulse at phase $\zeta_{\rm i}$ and moving with  velocity $v_z(\zeta_{\rm i})=v_{\rm p}$. At this point in phase space, the plasma wave potential is such that $\phi(\zeta_{\rm i})=\phi_{\rm i}\geqslant\phi_{\min}$. In addition, following EqHbasic, conservation of the Hamiltonian can be written  as $H_{\rm i}=(1+\beta_{\rm p}^2\gamma_{\rm p}^2)^{1/2}-\phi_{\rm i}-\beta_{\rm p}^2\gamma_{\rm p}=1/\gamma_{\rm p}-\phi_{\rm i}$. Since
$\phi_{\rm i}\geqslant\phi_{\min}$, we get $H_{\rm i}\leqslant H_{\rm sep}$, which is a necessary and sufficient condition for trapping.

The top panel plots the laser amplitude $a^2$ (red), the wakefield potential (blue), and the longitudinal electric field (black dashed); the bottom panel shows the associated electron trajectories in phase space. The black line in the bottom panel represents the fluid orbit, the red line the separatrix, and the dashed blue lines various trapped orbits. The $y$-axis is in log-scale and, to avoid negative values, we have plotted $u_z+1-\min(u_{\rm fluid})$, where $u_{\rm fluid}(\zeta)$ is the fluid orbit. The horizontal black dashed line represents $u_z=0$, while the horizontal red dashed line represents $u_z=\beta_{\rm p}\gamma_{\rm p}$. The parameters are $a_0=2$, $n_{\rm e}/n_{\rm c}=0.44\%$, $\lambda=0.8\:\mathrm{\mu m}$ and $\tau=20\:\mathrm{fs}$.

§.§ Trapping thresholds

Now that the basic trajectories have been established, it is possible to determine the trapping threshold of plasma electrons in a plasma wave. To do so, we consider that an electron located in front of the laser pulse will be trapped in the wakefield if its longitudinal momentum satisfies $u_z(+\infty)>u_{z}^{\rm sep}(+\infty)$, which simply means that the initial momentum has to be larger than the initial momentum on the separatrix, $u_{z}^{\rm sep}(+\infty)$. According to EqTraj, this is simply $u_{z}^{\rm sep}(+\infty)=\beta_{\rm p}\gamma_{\rm p}^2H_{\rm sep}-\gamma_{\rm p}\sqrt{\gamma_{\rm p}^2H_{\rm sep}^2-1}$. Therefore, an electron with initial energy $E>E_{\rm trap}$ will be trapped and accelerated in the wakefield, where

\begin{equation}
E_{\rm trap}=m_{\rm e}c^2\Bigl(\sqrt{1+\{u_z^{\rm sep}(+\infty)\}^2}-1\Bigr).
\end{equation}

Figure <ref> shows the variation of the trapping threshold with the plasma wave amplitude $\phi_{\min}$, as well as its variation with the plasma wave Lorentz factor $\gamma_{\rm p}$. Clearly, trapping is easier for high-amplitude plasma waves and/or for small phase velocities (small values of $\gamma_{\rm p}$). Note that as $\phi_{\min}\rightarrow -1$, the trapping threshold $E_{\rm trap}$ tends toward zero. This is the onset of wave-breaking: the plasma wave amplitude becomes so high that all plasma electrons that were initially at rest get injected into the plasma wave. As $\phi_{\min}\rightarrow -1$, the longitudinal electric field reaches the cold wave-breaking limit, $E_{\rm WB}=\sqrt{2}(\gamma_{\rm p}-~1)^{1/2}E_0$, where $E_0=m_{\rm e}c\,\omega_{\rm p}/e$ was defined earlier. Thus, this model is able to give a 1D picture of self-injection, although one has to keep in mind that in experiments self-injection is a 3D process. More complicated models are necessary to fully capture the physics of self-injection in 3D [16].

(a) Trapping thresholds plotted as a function of wake amplitude (where $\phi_{\min}$ represents the minimum of the wake potential), for three different wake phase velocities. (b) Trapping thresholds plotted as a function of $\gamma_{\rm p}$, the Lorentz factor associated with the wakefield velocity, for three different wake amplitudes.

§ IONIZATION INJECTION

As we have just seen, unless the wakefield reaches wave-breaking amplitudes, it is not possible for plasma electrons to be injected unless a different method is used. Ionization injection is probably the easiest method for injecting electrons in a wakefield [17, 18]. The idea is to use a high-$Z$ gas so that the first levels of ionization occur at low intensity (typically below $10^{16}\:\wcm$). These electrons are born at rest in a region where the laser intensity is relatively low and provide the electrons which form the plasma wave and follow fluid orbits. On the other hand, ionization from inner shells occurs at higher intensities (typically for $I>10^{18}\:\wcm$), so these electrons are born at rest in regions of strong fields. It follows that they are born inside the plasma wave, at a totally different phase from fluid electrons. The whole challenge is to put these electrons on trapped orbits so that they will be injected and accelerated. Figure <ref> illustrates the different classes of trajectories of  electrons born at the front of the pulse and in the middle of the pulse. In this section, we follow the approach of  [19, 17].

Principle of ionization injection. (a) Ionization of nitrogen by an intense laser pulse at $10^{19}\:\wcm$; the first five electrons, from N$^+$ to N$^{5+}$, appear at the front of the laser pulse, while electrons from N$^{6+}$ appear in the middle of the laser pulse. (b) Schematic showing the different classes of trajectories for electrons born at the front of the pulse (fluid trajectories) or in the middle of the pulse (trapped trajectories).

We now determine the condition(s) under which these electrons will be trapped. We first make a few assumptions concerning ionization. In our model,  ionization is assumed to occur via barrier suppression ionization [20]. Considering an ionization level with energy $E_{\rm i}$, the threshold intensity required to ionize the electron is given by

\[
I_{\rm thresh}\,[\!\wcm]=4\times 10^9\,E_{\rm i}^4\,{\rm [eV]}/Z^{*2},
\]

where $Z^*$ is the charge of the resulting ion. Furthermore, we assume that electrons are born at rest (which is a good approximation when considering the energy that the electron will gain in the laser and plasma field). Technically, if the electron is born at phase $\zeta_{\rm ion}$, it should witness the corresponding laser amplitude $a(\zeta_{\rm ion})$. However, in linear polarization, ionization occurs mostly at the peak of the electric field, i.e. $a(\zeta_{\rm ion})\simeq 0$. Therefore, in the case of an ionized electron born at rest, conservation of canonical momentum reads $u_\bot(\zeta)= a(\zeta)-a(\zeta_{\rm ion})\simeq a(\zeta)$. The initial Hamiltonian for such electrons can be found from EqHbasic as

\[
H_{\rm ion}=1-\phi_{\rm ion},
\]

and the ionized electron trajectory can be computed by inserting $H_{\rm ion}$ into EqTraj. The conditions for trapping are simply that the intensity should be high enough  for ionizing a given electron level, i.e. $a(\zeta_{\rm ion})>a_{\rm thresh}$, and that the electron should be born on a trapped orbit, i.e. $H_{\rm ion} < H_{\rm sep}$. These conditions define the region of phase space where ionized electrons are trapped and further accelerated in the wakefield. Typically, these conditions are fulfilled only in the setting of a high-intensity laser pulse ($a>1$) and a large-amplitude plasma wave. Figure <ref> illustrates this discussion: the green area in the upper panel defines the region of phase space where electrons can be injected; the corresponding trajectories are plotted in green in the lower panel (note that the fast laser frequency was considered for this case, as electrons are born at zeros of the vector potential). As discussed previously, electrons are born at rest, oscillate in the laser field, and gain energy as they get on a trapped orbit.

Ionization injection—a phase space picture. The upper panel plots the laser amplitude $a^2$ (red) and the wakefield potential (blue); the dashed light-blue horizontal line indicates the ionization threshold $a_{\rm thresh}=1.37$, corresponding here to $I_{\rm thresh}=4\times 10^{18}\:\wcm$, and the region of phase space where the injection conditions are satisfied is shown in green. The lower panel shows the fluid orbit (black), the separatrix (red), and the trajectories followed by ionized electrons at the front and back of the injection region. The parameters are $a_0=2$, $n_{\rm e}/n_{\rm c}=~0.44\%$, $\lambda=0.8\:\mathrm{\mu m}$ and $\tau=20\:\mathrm{fs}$.

The injection volume in phase space depends on the wakefield amplitude and the ionization threshold, so in principle it can be controlled via the laser intensity. The injection volume is responsible for the injected charge as well as the energy spread. Obviously, large injection volumes lead to large injected charge and larger energy spreads. On the other hand, smaller energy spreads can be obtained by reducing the injection volume at the expense of charge. Experiments have demonstrated the concept of ionization injection using nitrogen and argon [17, 18]. In practice, it is quite difficult to control the injection volume experimentally: ionization injection occurs all along the propagation, as long as the laser intensity exceeds the threshold intensity; as a result, electrons are injected continuously along the propagation and the energy distributions are quite broad [17, 18]. Narrower energy spreads have been obtained by using a two-stage laser–plasma accelerator [21]: in the first stage, a short nitrogen gas jet is used to inject electrons while a second gas jet is used to boost the acceleration. In Ref. [22], by restricting the first stage to a $500\mic$ jet, the injection was kept relatively local, resulting in energy spreads at the 10% level.

To summarize our discussion on ionization injection, this method is easy to implement experimentally as it requires simply the use of a high-$Z$ gas. It is a straightforward way to increase the injected charge without reaching wave-breaking amplitudes. It usually results in increased injected beam loads, but controlling the beam quality can be rather difficult, as (i) the injection volume is directly related to the laser intensity, and (ii) injection tends to occur continuously along the propagation.

§ COLLIDING-PULSE INJECTION

We now present the colliding pulse injection scheme, where an auxiliary laser pulse is used to trigger a very localized injection ($\simeq 10\mic$), resulting in beams with small energy spread. The idea of using an additional laser pulse for injecting electrons was first proposed by Umstadter et al. [23]. It was further developed by Esarey et al. [24], who proposed a scheme based on counter-propagating laser pulses. In its simplest form, the scheme uses two counter-propagating ultra-short laser pulses with the same central wavelength and polarization. The first laser pulse, the `pump' pulse, creates a wakefield, whereas the second laser, the `injection' pulse, is used only for injecting electrons. The laser pulses collide in the plasma, and their interference creates a laser beat-wave pattern which can pre-accelerate plasma background electrons. If the laser intensities are high enough, this pre-acceleration permits the injection and trapping of electrons in the wakefield and their further acceleration to relativistic energies. The principle of the method is depicted in figcollid. Analytical work [25] and simulations [25, 26] have shown that this two-stage acceleration mechanism can lead to the production of high-quality electron bunches with narrow energy spread, small divergence and ultra-short duration, even when using relatively modest lasers (e.g. with normalized vector potential $a_0=1$ for the pump pulse and $a_1=0.3$ for the injection pulse).

Principle of colliding-pulse injection

A theoretical description of electron injection by colliding laser
pulses was first given by Esarey et al. in Ref. [24].
We will briefly review the principles of their
fluid model. First, the pump laser, with normalized vector potential
$a_0$, excites a plasma wave. Before the arrival of the laser pulse, the plasma
background electrons are initially at rest. In the absence of
wave-breaking and self-trapping, these electrons are not trapped in
the plasma wave. As we have seen earlier, they simply oscillate back and forth in the
plasma wave potential along the fluid orbit, and their momenta
oscillate between $u_{z,\min}^{\rm fluid}$ and $u_{z,\max}^{\rm fluid}$, where

\begin{equation}\label{eq-fluid}
u_{z,\max,\min}^{\rm fluid}=\beta_{\rm p}\gamma_{\rm p}^2(1+\phi_{\min,\max})-\gamma_{\rm p}\sqrt{\gamma_{\rm p}^2(1+\phi_{\min,\max})^2-1},
\end{equation}

with $\phi_{\min}$ and $\phi_{\max}$ representing the minimum and maximum of the
potential, respectively.
In order to be accelerated in the
wakefield, electrons need to move along trapped orbits, inside the separatrix.
When electrons move along the
separatrix, their momenta oscillate between $u_{z,\min}^{\rm sep}$ and
$u_{z,\max}^{\rm sep}$, where

\begin{equation}\label{eq-fluid-sep}
u_{z,\max,\min}^{\rm sep}=\beta_{\rm p}\gamma_{\rm p}(1+\gamma_{\rm p}\Delta\phi)\pm\gamma_{\rm p}\sqrt{(1+\gamma_{\rm p}\Delta\phi)^2-1}
\end{equation}

with $\Delta\phi=\phi_{\max}-\phi_{\min}$.

Therefore, one needs to find a way to push plasma background
electrons, which follow the fluid trajectories, into trapped
trajectories. The interference of two laser pulses generates a
beat-wave which is able to heat electrons and provide just such a
mechanism. The laser pulses are represented by
their normalized vector potentials $\mathbf{a}_{0,1}$, where
the subscripts 0 and 1 represent, respectively, the pump pulse and the
injection pulse. Assuming that the lasers are counter-propagating
along the $z$ direction, one can write the normalized vector
potential as

\begin{equation}\label{eq-defa}
\mathbf{a}_{0,1}=\frac{a_{0,1}(\zeta_{0,1})}{\sqrt{1+\sigma}}\bigl[\cos(k_{0,1}\zeta_{0,1})\mathbf{e}_x+\sigma\sin(k_{0,1}\zeta_{0,1})\mathbf{e}_y\bigr],
\end{equation}

where $\sigma=0$ for linear polarization and $\sigma=1$ for
circular polarization, $k_0=-k_1$ are the wavevectors, and
$\zeta_{0,1}=z\pm v_{\rm g}t$. Although in general the two waves can have different frequencies, all experiments to date have used identical frequencies, so here we assume that the two laser pulses have the same frequency $\omega_0$.

When the two pulses overlap, they interfere, and the resulting
squared electromagnetic field can be written as
$\mathbf{a}^2=(\mathbf{a}_0^2+\mathbf{a}_1^2+2\mathbf{a}_0\cdot \mathbf{a}_1)/(1+\sigma)$. The last term is
the beat-wave term, and it cancels out for crossed polarizations. In
the case where the two polarizations are parallel, it generates a
standing wave (i.e. the beat-wave has zero phase velocity) with a
spatial scale of $\lambda_0/2$. The ponderomotive force in
the beat-wave is very large (because $F_{\rm beat}\simeq 2 k_0 a_0 a_1$), and
it can pre-accelerate the plasma electrons. Neglecting the plasma
potential, electron trajectories in the beat-wave are governed by
the Hamiltonian

\begin{equation}\label{eq-Hb}
H_{\rm b}=\sqrt{1+u_z^2+\mathbf{a}^2}.
\end{equation}

In order to obtain an analytical estimate, we will assume circular
polarization ($\sigma=1$) for the pump and injection beams, so that
$\mathbf{a}^2=(a_0^2+a_1^2)/2+a_0a_1\cos(2k_0z)$. The separatrix in
this beat-wave pattern is then given by

\begin{equation}\label{eq-sepb}
u_{z}^{\rm beat}=\pm\sqrt{a_0a_1(1-\cos{2k_0z})};
\end{equation}

so on the beat-wave separatrix, electrons oscillate between
$u_{z,\min}^{\rm beat}$ and $u_{z,\max}^{\rm beat}$ (see figbeatsep), where

\begin{equation}\label{uzmin-b}
u_{z,\max,\min}^{\rm beat}=\pm\sqrt{2a_0a_1}.
\end{equation}

For instance, an electron trapped in the beat-wave created by two laser pulses with $a_0=2$ and $a_1=0.3$ can gain about $200$ keV, which is sufficient for getting trapped in the wakefield.

Trajectories of electrons following the beat-wave separatrix at the collision of the two circularly polarized laser pulses.

One can define an approximate threshold for injection into the wakefield by applying a phase space separatrix overlap condition [24, 25]. Specifically, island overlap requires the following conditions:

\begin{align}
u_{z,\max}^{\rm beat} & >   u_{z,\min}^{\rm sep},\nonumber\\[-7pt]
& \label{eqcond} \\[-7pt]
u_{z,\min}^{\rm beat} & <  u_{z,\min}^{\rm fluid}.\nonumber
\end{align}

These conditions are illustrated in figcollidthresh(a). Using criterion (<ref>), we find that the injection threshold can be reached with $a_0\simeq 1$ and $a_1\simeq 0.1$ (for this estimation, we have used pulse durations of $\tau=30$ fs to calculate the plasma wave amplitude $\phi$), for plasma density in the range of $n_{\rm e}=10^{18}$–$10^{19}\:\mathrm{cm^{-3}}$; see the results in figcollidthresh(b).

(a) Phase space representation of the injection conditions in (<ref>) in the case of colliding pulses: the blue lines represent the amplitude of the beat-wave separatrix, while the red lines are the wakefield separatrix and the fluid orbit; the area shaded in green represents in injection volume. (b) Plots of the injection thresholds for various laser and plasma parameters (a 30 fs and 800 nm laser pulse was assumed).

When the pulse polarizations are linear, the Hamiltonian in eq-Hb is time-dependent and no longer integrable. It has been shown that in this case, electron motion is stochastic [27, 28, 29] and the simple calculations above do not hold. In fact, stochastic heating turns out to be an efficient mechanism for heating and subsequently injecting electrons [28, 30]. Other effects, such as plasma wave inhibition [31], make the physics of colliding pulses more complicated than what this simple model tends to suggest. Full 3D particle-in-cell simulations are necessary to fully capture all the physical ingredients [30], but the concepts presented here still hold.

The colliding-pulse injection method requires some heavy experimental investment, as a second intense laser pulse is needed. The two laser pulses need to be precisely synchronized at the femtosecond level and spatially overlapped at the micron level. Despite its complexity, colliding-pulse injection has intrinsic advantages, some of which are summarized as follows.

* The injection is very localized as it occurs only during the collision of the two laser pulses, i.e. $10\mic$ for 30 fs pulses. This helps to reduce the energy spread and to inject very short electron bunches. The injected electron beams were measured to be only a few femtoseconds in duration [32].

* The physics of injection is largely determined by the heating at the collision, so injection is less sensitive to intensity fluctuations of the main laser pulse (though the intensity of the main pulse still plays a role, as it determines the wakefield amplitude). First experiments have shown that stable beams with narrow energy spreads can be achieved in this way [33, 34].

* The location of the collision and hence the injection can be easily controlled. Therefore it is possible to tune the acceleration length. Experiments have shown that by controlling this parameter, the beam energy can be controlled over a wide range, typically 10–200 MeV [33].

* The injection volume can be controlled by tuning the intensity of the injection laser pulse. For high $a_1$, the electron heating at the collision is large and results in large injected beam loads as well as larger energy spreads. The energy spread can be tuned and reduced at the 1% level by decreasing the injection laser amplitude, as demonstrated experimentally in Ref. [35].

Finally, based on theoretical studies, a three-pulse scheme has been proposed where the main beam is  used only to generate the wakefield [24]; two auxiliary beams are used for injection only. Their polarization is orthogonal to the main beam polarization so that the two injection pulses do not interfere with the main beam. The advantage of this method is that it provides an additional tuning knob: the injection phase can now be chosen by tuning the distance between the main pulse and the collision point. More recently, several variations of the colliding-pulse scheme have been proposed; the basic concepts remain similar but the analysis extends to the 3D case [36, 37].

§ INJECTION IN DENSITY GRADIENTS

We now focus on another scheme, in which the plasma needs to be engineered: by tailoring the plasma density, one can gain some control over the plasma wave phase velocity. As we have seen earlier, the lower the phase velocity, the lower the trapping threshold. Therefore, by setting up a local decrease in the phase velocity, one can trigger injection in a local manner as well [38, 39, 40]. This can be achieved by sending the laser pulse through a downward density ramp, which causes the wakefield to slow down. A full Hamiltonian description of this problem is beyond the scope of this article; in a density transition, the wakefield potential also depends on $z$, as $\phi(z,\zeta)$, and the constants of motion that we have derived in Section <ref> no longer hold in this case. Consequently, we will develop a simple fluid model that will provide some physical intuition of this process.

We start with the fluid equation describing the excitation of the wakefield by an intense laser pulse, restricting ourselves to the weakly relativistic case where $a^2\ll1$. We consider a gentle density gradient, $\frac{1}{n_0}\frac{\partial n_0}{\partial z}\ll k_{\rm p}$, or $L_{\rm g}k_{\rm p}\ll1$, where $L_{\rm g}$ is the gradient scale length. In this case, the plasma wave equation reads

\begin{equation}
\left(\frac{\partial^2}{\partial t^2}+\omega_{\rm p}^2(z)
\right)\phi=\omega_{\rm p}^2(z)\frac{\langle a^2\rangle}{2}.
\end{equation}

Note the spatial dependence of the plasma frequency in the gradient, $\omega_{\rm p}(z)\propto \sqrt{n_0(z)}$. We can now perform the usual change of variables to solve the problem in the moving frame: $\zeta=z-v_{\rm g}t$ and $\tau=t$. By applying the quasi-static approximation, one can then neglect the $\partial /\partial\tau$ terms. The previous equation, written in the co-moving variables, becomes

\begin{equation}
\left(\frac{\partial^2}{\partial \zeta^2}+k_{\rm p}^2(z)
\right)\phi=k_{\rm p}^2(z)\frac{\langle a^2\rangle}{2}
\end{equation}

where $k_{\rm p}(z)=\omega_{\rm p}(z)$. This equation can be integrated, and the solution behind the laser pulse has the form

\begin{equation}\label{eqphigrad}
\phi(\zeta,z)=\phi_0(z)\sin[k_{\rm p}(z)(z-v_{\rm g}t)],
\end{equation}

where the wakefield amplitude is $\phi_0(z)=-\frac{\sqrt{\pi}}{4}a_0^2(z)k_{\rm p}(z)L_0\exp\{-k_{\rm p}(z)^2L_0^2/4\}$ and its phase is $\varphi=k_{\rm p}(z)(z-v_{\rm g}t)$, so that one can compute the local oscillation frequency and wavevector:

\begin{align}
\omega & = -\partial \varphi/\partial t=k_{\rm p}(z)/v_{\rm g}=\omega_{\rm p}(z),\\
k & = \partial \varphi/\partial z=k_{\rm p}(z)+\partial k_{\rm p}/\partial z (z-v_{\rm g} t).
\end{align}

It follows that in a downward density gradient ($\partial k_{\rm p}/\partial z<0$ and $z-v_{\rm g}t <0$ behind the laser pulse), the wavevector increases with time (alternatively, the plasma wavelength decreases with time). In contrast, the plasma frequency does not depend on time: $\omega=\omega_{\rm p}(z)$. As a result of this time-varying wavevector, the phase velocity $v_{\rm p}(z,t)=\omega_{\rm p}(z)/k(z,t)$ is

\begin{equation}\label{eqvp}
v_{\rm p}(z,t) = \frac{v_{\rm g}}{1+(z-v_{\rm g}t)\dfrac{1}{k_{\rm p}}\dfrac{\rmd k_{\rm p}}{\rmd z}}.
\end{equation}

Consequently, as the wavevector increases with time, the phase velocity decreases. Injection occurs behind the laser pulse when the wakefield becomes slow enough to trap plasma background electrons.

These effects are illustrated in figgrad. In panel (a), we display the wakefield potential in the case where a laser pulse is tightly focused into a density gradient. The density gradient has a Gaussian shape and $k_{\rm p0}L_{\rm g}=50$, where $k_{\rm p0}=\omega_{\rm p0}/v_{\rm g0}$ is defined using the maximum density at the top of the downward ramp. One can see clearly that the plasma frequency decreases with $z$ as the density decreases, as expected. In addition, while the initial phase velocity is close to $c$, the bending of the wakefield in the $(t,z)$ plane indicates that the phase velocity decreases with time. In panel b), the phase velocity is plotted for various density scale lengths $L_{\rm g}$. It can be seen that the phase velocity decreases faster for sharper gradients, indicating that trapping is likely to occur more rapidly in this case.

(a) Wakefield potential map $\phi(z,t)$ in the case where a laser pulse is focused in a density gradient; eqphigrad was used to compute the potential; the dashed line represents the speed of light and one can see clearly the bending of the wakefield in the $(t,z)$ plane, indicating that the phase velocity is decreasing with time. (b) Illustration of the decrease in the phase velocity with time for two different gradient scale lengths $L_{\rm g}$ (computed using eqvp).

Numerous experiments have demonstrated that trapping in density gradients is an efficient method for injecting electrons into wakefields [41, 42, 43, 44]. This controlled injection scheme has resulted in more stable beams [41, 43] with energy spreads in the $10\%$ range. Some experiments have utilized the density down ramp at the exits of short gas jets, resulting in the production of low-energy beams of a few hundred keV [41, 44]. Other research groups have engineered a sharp density gradient by creating a shock in the gas flow [43, 45] or by using another laser pulse to create a density perturbation [42]. In this case, the injection location can be controlled and the resulting beam energy can be tuned [45, 46].

§ CONCLUSION

In conclusion, plasma injection schemes are at the forefront of current research in laser–plasma accelerators. Researchers continue to propose new injection schemes to test new ideas and to provide high-quality beams while maintaining a relatively simple experimental set-up [47, 48]. The idea of combining several methods might prove useful in the future; for example, performing colliding-pulse injection in a density gradient could increase the trapped charge while still preserving the beam quality[49]. Similarly, combining ionization injection and density gradient injection has the potential to yield interesting results.

§ ACKNOWLEDGEMENTS

J. Faure's current research on laser–plasma acceleration is supported by the European Research Council under contract no. 306708, ERC Starting Grant FEMTOELEC.

[1]
M. Everett, A. Lal, D. Gordon, C. Clayton, K. Marsh and C. Joshi,
Nature 368 (1994) 527. http://dx.doi.org/10.1038/368527a0

[2]
F. Amiranoff, S. Baton, D. Bernard, B. Cros, D. Descamps, F. Dorchies,
F. Jacquet, V. Malka, G. Matthieussent, J.R. Marquès et al.,
Phys. Rev. Lett. 81 (1998) 995.

[3]
A. Modena, A. Dangor, Z. Najmudin, C. Clayton, K. Marsh, C. Joshi, V. Malka,
C. Darrow, D. Neely and F. Walsh,
Nature 377 (1995) 606. http://dx.doi.org/10.1038/377606a0

[4]
D. Umstadter, S.-Y. Chen, A. Maksimchuk, G. Mourou and R. Wagner,
Science 273 (1996) 472. http://dx.doi.org/10.1126/science.273.5274.472

[5]
J. Faure, Y. Glinec, A. Pukhov, S. Kiselev, S. Gordienko, E. Lefebvre, J.-P.
Rousseau, F. Burgy and V. Malka,
Nature 431 (2004) 541. http://dx.doi.org/10.1038/nature02963

[6]
C.G.R. Geddes, C.Tóth, J. van Tilborg, E. Esarey, C.B. Schroeder,
D. Bruhwiler, C. Nieter, J. Cary and W.P. Leemans,
Nature 431 (2004) 538. http://dx.doi.org/10.1038/nature02900

[7]
S.P.D. Mangles, C.D. Murphy, Z. Najmudin, A.G.R. Thomas, J.L. Collier,
A.E. Dangor, E.J. Divall, P.S. Foster, J.G. Gallacher, C.J. Hooker
et al.,
Nature 431 (2004) 535.

[8]
G. Sun, E. Ott, Y.C. Lee and P. Guzdar,
Phys. Fluids 30 (1987) 526.

[9]
A.B. Borisov, A.V. Borovskiy, O.B. Shiryaev, V.V. Korobkin, A.M.
Prokhorov, J.C. Solem, T.S. Luk, K. Boyer and C.K. Rhodes,
Phys. Rev. A45 (1992) 5830.

[10]
E. Esarey, C.B. Schroeder, B.A. Shadwick, J.S. Wurtele and W.P. Leemans,
Phys. Rev. Lett. 84 (2000) 3081. http://dx.doi.org/10.1103/PhysRevLett.84.3081

[11]
D.F. Gordon and B. Hafizi,
Phys. Rev. Lett. 90 (2003)  215001.

[12]
J. Faure, Y. Glinec, J.J. Santos, F. Ewald, J.-P. Rousseau, S. Kiselev,
A. Pukhov, T. Hosokai and V. Malka,
Phys. Rev. Lett. 95 (2005) 205003.

[13]
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,
Phys. Rev. Lett. 105 (2010) 235003.

[14]
E. Esarey and M. Pilloff,
Phys. Plasmas 2 (1995) 1432. http://dx.doi.org/10.1063/1.871358

[15]
E. Esarey, P. Sprangle, J. Krall and A. Ting,
IEEE Trans. Plasma Sci. 24 (1996) 252.

[16]
I. Kostyukov, A. Pukhov and S. Kiselev,
Phys. Plasmas 11 (2004) 5256.

[17]
A. Pak, K.A. Marsh, S.F. Martins, W. Lu, W.B. Mori and C. Joshi,
Phys. Rev. Lett. 104 (2010) 025003. http://dx.doi.org/10.1103/PhysRevLett.104.025003

[18]
C. McGuffey, A.G.R. Thomas, W. Schumaker, T. Matsuoka, V. Chvykov, F.J.
Dollar, G. Kalintchenko, V. Yanovsky, A. Maksimchuk, K. Krushelnick
et al.,
Phys. Rev. Lett. 104 (2010) 025004. http://dx.doi.org/10.1103/PhysRevLett.104.025004

[19]
M. Chen, E. Esarey, C.B. Schroeder, C.G.R. Geddes and W.P. Leemans,
Phys. Plasmas 19 (2012) 033101. http://dx.doi.org/10.1063/1.3689922

[20]
M.V. Ammosov, N.B. Delone and V.P. Krainov,
Sov. Phys. JETP 64 (1986) 1191.

[21]
B.B. Pollock, C.E. Clayton, J.E. Ralph, F. Albert, A. Davidson, L. Divol,
C. Filip, S.H. Glenzer, K. Herpoldt, W. Lu et al.,
Phys. Rev. Lett. 107 (2011) 045001.

[22]
G. Golovin, S. Chen, N. Powers, C. Liu, S. Banerjee, J. Zhang, M. Zeng,
Z. Sheng and D. Umstadter,
Phys. Rev. ST Accel. Beams 18 (2015) 011301.

[23]
D. Umstadter, J.-K. Kim and E. Dodd,
Phys. Rev. Lett. 76 (1996) 2073.

[24]
E. Esarey, R.F. Hubbard, W.P. Leemans, A. Ting and P. Sprangle,
Phys. Rev. Lett. 79 (1997) 2682. http://dx.doi.org/10.1103/PhysRevLett.79.2682

[25]
G. Fubiani, E. Esarey, C. Schroeder and W. Leemans,
Phys. Rev. E70 (2004) 016402.

[26]
H. Kotaki, S. Masuda, M. Kando, J. Koga and K. Nakajima,
Phys. Plasmas 11 (2004) 3296.

[27]
J.T. Mendonça,
Phys. Rev. A28 (1983) 3592. http://dx.doi.org/10.1103/PhysRevA.28.3592

[28]
Z.-M. Sheng, K. Mima, Y. Sentoku, M.S. Jovanović, T. Taguchi, J. Zhang
and J.M. ter Vehn,
Phys. Rev. Lett. 88 (2002) 055004. http://dx.doi.org/10.1103/PhysRevLett.88.055004

[29]
A. Bourdier, D. Patin and E. Lefebvre,
Physica D206 (2005) 1.

[30]
X. Davoine, E. Lefebvre, J. Faure, C. Rechatin, A. Lifschitz and V. Malka,
Phys. Plasmas 15 (2008) 113102. http://dx.doi.org/10.1063/1.3008051

[31]
C. Rechatin, J. Faure, A. Lifschitz, V. Malka and E. Lefebvre,
Phys. Plasmas 14 (2007) 060702. http://dx.doi.org/10.1063/1.2741387

[32]
O. Lundh, J. Lim, C. Rechatin, L. Ammoura, A. Ben-Ismail, X. Davoine,
G. Gallot, J.-P. Goddet, E. Lefebvre, V. Malka et al.,
Nat. Phys. 7 (2011) 219. http://dx.doi.org/10.1038/nphys1872

[33]
J. Faure, C. Rechatin, A. Norlin, A. Lifschitz, Y. Glinec and V. Malka,
Nature 444 (2006) 737.

[34]
J. Faure, C. Rechatin, A. Norlin, F. Burgy, A. Tafzi, J.-P. Rousseau and
V. Malka,
Plasma Phys. Control. Fusion 49 (2007) B395. http://dx.doi.org/10.1088/0741-3335/49/12B/S36

[35]
C. Rechatin, J. Faure, A. Ben-Ismail, J. Lim, R. Fitour, A. Specka, H. Videau,
A. Tafzi, F. Burgy and V. Malka,
Phys. Rev. Lett. 102 (2009) 164801. http://dx.doi.org/10.1103/PhysRevLett.102.164801

[36]
X. Davoine, E. Lefebvre, C. Rechatin, J. Faure and V. Malka,
Phys. Rev. Lett. 102 (2009) 065001. http://dx.doi.org/10.1103/PhysRevLett.102.065001

[37]
R. Lehe, A.F. Lifschitz, X. Davoine, C. Thaury and V. Malka,
Phys. Rev. Lett. 111 (2013) 085005. http://dx.doi.org/10.1103/PhysRevLett.111.085005

[38]
S. Bulanov, N. Naumova, F. Pegoraro and J. Sakai,
Phys. Rev. E58 (1998) R5257.

[39]
H. Suk, N. Barov, J.B. Rosenzweig and E. Esarey,
Phys. Rev. Lett. 86 (2001) 1011.

[40]
A.V. Brantov, T.Z. Esirkepov, M. Kando, H. Kotaki, V.Y. Bychenkov and S.V.
Phys. Plasmas 15 (2008) 073111. http://dx.doi.org/10.1063/1.2956989

[41]
C.G.R. Geddes, K. Nakamura, G.R. Plateau, C. Tóth, E. Cormier-Michel,
E. Esarey, C.B. Schroeder, J.R. Cary and W.P. Leemans,
Phys. Rev. Lett. 100 (2008) 215004.

[42]
J. Faure, C. Rechatin, O. Lundh, L. Ammoura and V. Malka,
Phys. Plasmas 17 (2010) 083107.

[43]
K. Schmid, A. Buck, C.M.S. Sears, J.M. Mikhailova, R. Tautz, D. Herrmann,
M. Geissler, F. Krausz and L. Veisz,
Phys. Rev. ST Accel. Beams 13 (2010) 091301.

[44]
Z.-H. He, B. Hou, J.H. Easter, J. Faure, K. Krushelnick, J.A. Nees and
A.G.R. Thomas,
New J. Phys. 15 (2013) 053016. http://dx.doi.org/10.1088/1367-2630/15/5/053016

[45]
A. Buck, J. Wenz, J. Xu, K. Khrennikov, K. Schmid, M. Heigoldt, J.M.
Mikhailova, M. Geissler, B. Shen, F. Krausz et al.,
Phys. Rev. Lett. 110 (2013) 185006.

[46]
P. Brijesh, C. Thaury, K.T. Phuoc, S. Corde, G. Lambert, V. Malka, S.P.D.
Mangles, M. Bloom and S. Kneip,
Phys. Plasmas 19 (2012) 063104. http://dx.doi.org/10.1063/1.4725421

[47]
N. Bourgeois, J. Cowley and S.M. Hooker,
Phys. Rev. Lett. 111 (2013) 155004.

[48]
L.-L. Yu, E. Esarey, C.B. Schroeder, J.-L. Vay, C. Benedetti, C.G.R. Geddes,
M. Chen and W.P. Leemans,
Phys. Rev. Lett. 112 (2014) 125001. http://dx.doi.org/10.1103/PhysRevLett.112.125001

[49]
G. Fubiani, E. Esarey, C.B. Schroeder and W.P. Leemans,
Phys. Rev. E73 (2006) 026402.