A multidimensional theory for electron trapping by a plasma wake generated in the bubble regime

We present a theory for electron self-injection in nonlinear, multidimensional plasma waves excited by a short laser pulse in the bubble regime or by a short electron beam in the blowout regime. In these regimes, which are typical for electron acceleration in the last impressive experiments, the laser radiation pressure or the electron beam charge pushes out plasma electrons from some region, forming a plasma cavity or a bubble with a huge ion charge. The plasma electrons can be trapped in the bubble and accelerated by the plasma wakefields up to a very high energy. We derive the condition of the electron trapping in the bubble. The developed theory predicts the trapping cross section in terms of the bubble radius and the bubble velocity. It is found that the dynamic bubble deformations observed in the three-dimensional (3D) particle-in-cell (PIC) simulations influence the trapping process significantly. The bubble elongation reduces the gamma-factor of the bubble, thereby strongly enhancing self-injection. The obtained analytical results are in good agreement with the 3D PIC simulations.


Introduction
Intense laser-plasma and beam-plasma interactions are highly nonlinear phenomena, which, besides being of fundamental interest, attract a great deal of attention owing to a number of important applications. One of the key applications is particle acceleration based on the excitation of the strong plasma wakefield by a laser pulse or by an electron beam. The first acceleration regime is called laser wakefield acceleration (LWFA) [1] and the second one is plasma wakefield acceleration (PWFA) [2,3]. The phase velocity of the plasma wake is close to the speed of light, so the relativistic charged particles loaded in the acceleration phase of the wake can be accelerated during a long period of time to very high energy. In the linear regime of interaction, when the laser intensity is low or the electron beam charge is small, the plasma wake is a linear plasma wave.
If the laser intensity is sufficiently high, then plasma wave breaking can occur. It has been shown by Akhiezer and Polovin [4] that wave breaking occurs in the cold plasma when the electric field of the plasma wake reaches the value E (2γ 0 − 2)mcω p /e, where ω p = 4π n 0 e 2 /m is the (non-relativistic) plasma frequency for the background electron density n 0 , γ 0 = (1 − V 2 /c 2 ) −1/2 , V is the wake phase velocity and m, e and c are the electron mass, the electron charge and the speed of light, respectively. Moreover, the ponderomotive force of the laser pulse pushes out the plasma electrons from the high intensity region, leaving behind the laser pulse the plasma cavity-the bubble, which is almost free from plasma electrons. This is the bubble regime of the laser-plasma interaction [5]. A similar regime, which is usually called the blowout regime, exists for beam-plasma interaction [6]. In the blowout regime, the plasma cavity is formed due to the space charge effect from the electron beam.
Although the bubble propagates with a velocity that is close to the speed of light, the large charge of unshielded ions inside the plasma cavity can trap the cold plasma electrons. Moreover, the electrons are trapped inside the accelerated phase of the bubble plasma field, thereby leading to efficient electron acceleration. Electron self-injection is an important advantage of plasmabased acceleration, which allows us to exclude the beam loading system requiring accurate synchronization and additional space. Recent experiments on LWFA [7]- [10] and PWFA [11] have demonstrated the high efficiency of electron self-injection. The charge and the energy of the electron bunches self-injected and accelerated in the bubble have ranged to a maximum 3 of about nC and 1 GeV, respectively. The self-injection also strongly affects the quality of the resulting bunch of accelerated electrons. The beam quality is often of crucial importance in many applications ranging from inertial confinement fusion to x-ray free-electron lasers.
The theory of one-dimensional (1D) cold relativistic plasma waves where the threshold of wave breaking was considered has been presented by Akhiezer and Polovin [4]. Later, Dawson used the Lagrangian approach to study nonlinear plasma oscillations and to associate the trajectory crossing with plasma wave breaking [12]. The 1D model based on Hamiltonian formalism has been proposed to study electron trapping in a relativistic plasma wave [13,14]. The model predicts the threshold electron momentum required for trapping and the fraction of the trapped electrons. Several optical and plasma techniques have been proposed to enhance electron trapping in plasma wakefield such as the collision of two counter-propagating laser pulses [15,16], density transitions [17,18], ionization in the wake driven by relativistic electron bunch [19], etc. However, electron self-trapping is mainly relevant to the bubble or blowout regimes, when 1D models are no longer valid.
The typical electron density distribution in the bubble regime calculated by 3D particle in cell (PIC) simulation using VLPL code [20] is plotted in figure 1(a). The laser pulse propagates along the x-axis and is linearly polarized along the y-axis. The pulse envelope is a = a 0 exp(−r 2 /r 2 L ) cos(π t/T L ), where T L = 30 fs is the pulse duration, r L = 9 µm is the focused spot size and a 0 ≡ e A/mc 2 = 1.5 is the normalized vector potential, which corresponds to the laser intensity I = 3 × 10 18 W cm −2 . It is seen from figure 1(a) that the bubble is formed when laser pulse passes 900λ. There are typical parts of a bubble: plasma cavity, electron sheath around the cavity, dense plasma at the bubble front, the wave-breaking area at the trailing edge of the bubble and the dense bunch of the relativistic electrons inside the cavity. An example of the trajectory of a trapped electron is shown in figure 1(b). It follows from figure 1 that neither the plasma response structure nor the trajectory of the trapped electrons can be described by the 1D model.
There is a number of analytical models describing the electron dynamics of untrapped electrons [21,29] and dynamics of the ultrarelativistic electrons after self-injection [22]- [25]. However, despite the great interest, there is still little theory for electron trapping including non-1D electron dynamics. The model where the electromagnetic field in the wake was treated as the Coulomb field has been presented in [26]. The electron self-trapping is considered as the trajectory winding around the Coulomb center in the relativistic limit. However, the structure of the electromagnetic field in the bubble strongly differs from the Coulomb field [27]- [29]. The electron self-injection has extensively been studied numerically [28,30,31]. It was found in [31] that a sufficiently large bubble, R > 4, can trap plasma electrons. Here, R is the bubble radius normalized to c/ω p . A different condition for electron self-injection has been proposed in [28]: R > γ 0 , where γ 0 is the bubble gamma factor. Recently the effect of the bubble expansion on the electron self-injection has been studied with numerical simulations [37]. However, the numerical approach has its natural limitations in accuracy and validity conditions of the obtained results and does not provide insight into the self-injection physics. In this paper, we present an analytical model for self-injection and verify it by direct measurements of model parameters in 3D PIC simulations.
The paper is organized as follows. The basic equations and theory propositions are formulated in section 2. In section 3 we consider the piecewise model of the bubble field. The self-trapping condition and the cross-section of the electron trapping are derived. In section 4, we use the perturbation theory to study the self-injection dynamics in the linear approximation for the bubble field. In section 5, we compare the theory predictions with the results of the PIC simulations. In section 6, we discuss the obtained results.

Basic equations
The self-consistent description of electron self-injection is a very complex problem. Here, we present a simplified model that does keep the main features of self-injection physics. The model describes the electron trapping by a bubble. A bubble is a plasma cavity, which is free from plasma electrons and moves in plasma with velocity close to the speed of light. The plasma electrons, which are initially at rest, are scattered by the bubble. Part of the plasma electrons can be trapped in the bubble. To make analytical advance, we adopt the following assumptions: (i) the bubble is a sphere with large radius R c/ω p , which holds for short and intense drivers [29]; (ii) the bubble motion is ultrarelativistic γ 0 = (1 − V 2 /c 2 ) −1/2 1, where γ 0 is the gamma factor of the bubble, V is the bubble velocity; (iii) the electron motion is planar because of axial symmetry of the driver; (iv) the driver field can be neglected because it is weak in the region where the trapping occurs (far behind the laser pulses or electron beams); (v) the field of the trapped electrons is also neglected as the number of trapped electrons is assumed to be small.
To calculate the trajectory of the trapped electron, we have to know the distribution of the electromagnetic field in the bubble. It can be shown [28,29] that the space-time distribution of the electromagnetic field inside the spherical bubble can be approximated as follows: where V is the bubble velocity that is close to the speed of light (V 1). We use dimensionless units, normalizing the time to ω −1 p , the lengths to c/ω p , the velocity to c, the electromagnetic fields to mcω p /|e|, and the electron density, n, to the background plasma density n 0 , respectively. The distribution of the electromagnetic field inside the plasma cavity can also be expressed through vector and scalar potentials [28] where gauge A x = −ϕ is used. Furthermore, we assume that the electron trajectory is planar (z = 0) because of axial symmetry of the wakefield and the driver. The driver field can be neglected because it is weak in the region where trapping occurs (far behind the laser pulse or electron beam). The wake field of the trapped particles is also neglected. The electron motion in the plasma cavity is governed by the Hamiltonian where is the canonical momentum of the electron. We change variables in the Hamiltonian (6) from x and x to ζ = x − V t and ζ = x by a canonical transformation with the generating function S = (x − V t) ζ . The Hamiltonian takes the following form in the new variables: where γ = 1 + p 2 x + p 2 y is the relativistic gamma-factor of the electron and p = vγ = + A is the kinetic momentum of the electron. After canonical transformation the Hamiltonian does not depend on time and it is an integral of motion.
We assume that the electron is initially immobile at the bubble border p(t = 0) = 0, ζ (t = 0) = R 2 − ρ 2 , y(t = 0) = ρ, where ρ is the impact parameter of the electron. It follows from the initial condition that H = 1 + R 2 /4. The Hamiltonian equations are Integrating equations (8)- (11), the electron trajectory can be calculated. We will use the following condition for electron self-injection: the electron is trapped in the bubble if r = ζ 2 + y 2 < R for a long period of time. The appropriate time period will be defined later.

Piecewise model
Before studying the electron dynamics in the bubble field given by equations (1)-(4), we simplify our calculations to distill the self-injection physics. First we consider the piecewise field model. We will show that the model is integrable, which makes analytical calculations of the electron trajectories and the trapping cross-section possible. The model qualitatively reproduces some characteristic features of 3D trapping and shows the differences in trapping processes for the 1D and 3D cases. In the piecewise model, electron motion inside the plasma cavity obeys the following equations: d p y dt = −sign(y) where F is a constant approximating the magnitude of electromagnetic fields in the plasma cavity. The model is based on the piecewise approximation of the bubble field. In each quadrant (a quarter of the plane ζ -y where the coordinates ζ and y do not change sign) the electromagnetic fields are constant and change sign when the particle crosses the quadrant border y = 0 and ζ = 0. The field inside the bubble is generated by the current density and the charge density in the planes ζ = 0 and y = 0. The typical trajectory of the trapped electron in the piecewise model is shown in figure 2(a). It can be seen from the figure that the trapped electron undergoes betatron oscillations inside the bubble. The Hamiltonian in the piecewise model is It does not depend on time explicitly; hence, it is conserved. The electron at the time instant t = 0 is immobile and is situated at the border of the bubble, implying that p

These initial conditions yield
where we assume that y 0 > 0 and ζ 0 > 0. The solution of equations (12)-(15) for the initial conditions corresponds to the oscillating electron motion (see figure 2). We introduce the instants of time t n when the electron crosses the ζ -axis and the instants of time τ n when |y| peaks and p y = 0. We choose numbering such that the time instant t 1 corresponds to the first crossing of the electron trajectory and the ζ -axis, and τ n ∈ [t n , t n + 1].
To derive the trapping condition, we consider now the dynamics of the trapped electrons when the number of betatron oscillations is large, n 1, because this limit is decisive for electron trapping, as will be seen later. For equations (12)- (15), there is an additional integral of motion in each quarter of the plane ζ -y where the electron coordinates do not change sign. This integral can be considered as the second Hamiltonian The Hamiltonian determines the electron motion in each quadrant, and the equations for electron motion can be presented as Hamiltonian equations We introduce h n , that is, the value of h in the time interval t n t t n+1 . For the time interval t ∈ (t n , t n+1 ) from equations (12)- (15) and (19), we obtain the following equation: with the solution where χ = p y Fsign(ζ )/(2h n ) and χ n = χ(t = t n ). It follows from the numerical solution of equations (12)-(15) that | p y / p x | 1 and hence |χ| 1 for n 1. Introducing the function µ = χ n (χ n − χ)/(1 − χ n ) and expanding the righthand side of equation (21) in a series of χ, we obtain The time instances t n and τ n are introduced such that Making use of these relations, we can find As the sign of y(τ n ) is opposite to the sign of h n , we can write To express g n through g n+1 , we consider two instants of time t n+1 − δt and t n+1 + δt, where δt is infinitesimally small Then, we can obtain the relationship between I n+1 and I n from equations (26): The right-hand sides of equations (28) and (30) are much smaller than g n and I n , respectively, for n 1. Therefore, in the limit n 1, we can rewrite the relations for I n and g n in the form of a system of ordinary differential equations dg n dn = I n , Dividing the second equation by the first one, we obtain dI n dg n = I n 3g n .
The solutions of the system are the following: where is a constant and l characterizes the initial conditions of equations (31). It follows from equations (34) and (35) that g n and I n increase as n increases. Under the assumptions p x (τ n ) 1 from the definitions of g n and χ n , we have From equations (26), (33) and (37), we obtain The evaluation of the constant C is presented in the appendix. The piecewise model makes the analytic treatment of the electron trajectory possible in the limit n 1. From equations (26), (33), (29) and (37) and the appendix, the coordinate and the momentum of the trapped electrons as a function of the oscillation number n can be derived: To present the coordinate and the momentum of the trapped electrons as a function of t, we have to find the relationship between t n , τ n and n. It can be done from the definition of t n , τ n defined by equations (25) and (24). However, a simpler way is to use equation (12). In the limit n 1, the longitudinal momentum of the electron becomes much more than the transversal momentum so that the second term in equation (12) can be neglected. Thus, we obtain where V 1 is assumed. The transverse coordinate and momentum of the trapped electron as a function of t can be derived in the limit n 1 from equations (39)-(41) and (42) It follows from the obtained equations that the momentum components of the trapped electron increase with time, whereas the amplitude of the betatron oscillations decreases. The derived expression is in good agreement with the results of numerical integration of equations of motion (12)-(15) (see figure 3). The obtained results can be used to derive the trapping condition. Making use of equations (39)-(41), y(τ n ) can be written as a function of p x (τ n ) Assuming |y(τ n )| R, the distance r = ζ 2 + y 2 between the bubble center and the electron can be approximated as follows: which with equation (16) yields Substituting H from equation (17) and neglecting the term y 2 /(2R) in comparison with R as |y| R and the term 1/ p x in comparison with p x /γ 2 0 (we assume that p x γ 0 , which is valid for wide areas of parameters, as shown later), we have The trapping condition can be formulated as follows. The electron is trapped if r remains less than R while ζ < 0. As the numerical simulation shows, the local maxima of r are reached at instants of time close to the instants τ n (when |y| peaks). Thus, for the local maxima of r , equations (48) and (45) yield If p x is less than some characteristic value p x;m , the local maxima of r increase due to a decrease of |y(τ n )|, and then p x becomes more than p x;m and r (τ n ) decreases due to an increase of p x . Thus, r reaches its maximal value r max when p x = p x;m . Treating the p x (τ n ) as a continuous variable, we obtain The condition r max < R yields the quadratic equation for ζ 0 , which under the assumption ζ 0 R gives the following trapping condition: The initial coordinates of the trapped electron are ζ 0 > 0 and y 0 = R 2 − ζ 2 0 . Thus, the condition (52) implies that plasma electrons are trapped from the thin layer y 0 ∈ ( √ R 2 − 2 , R), where = 0.98R/γ 1/2 0 − 2/F. From ζ 0 > 0 and equation (52), we obtain the conditions for existence of this layer: Therefore the bubble velocity must be less than some threshold value in order that electron trapping occurs. If certain γ 0 , R and F obey equation (53), the trapping cross-section is given by It can be seen from figure 4 that the trapping condition and the trapping cross-section are in good agreement with the numerical solution of equations (12)- (15). The electron trapping in the piecewise model is substantially non-1D. In the case of 1D electron motion, the trapping condition is p x > γ 0 [21]. For typical parameters of the laser-plasma acceleration experiments F R 1 so that p x;m γ 0 (see equation (50)). Thus, the electron has to gain much greater energy for trapping than it follows from the 1D model.

Bubble field model
We now return to the bubble field distribution given by equations (1)-(4). Equations (8) and (9) for the electron motion in the bubble field can be solved numerically. To be more realistic and in accordance with PIC simulations, we include an electron sheath around the plasma cavity, which screens the bubble ion field in the surrounding plasma. We model the electromagnetic fields inside the bubble as follows: where r 2 = ζ 2 + y 2 and d is the width of the electron sheath. This force is close to the one observed in the 3D PIC simulation. The potential of this structure is where Li 2 (z) = 0 z (ln(1 − t)/t)dt is the dilogarithm function [32]. In the limit d → 0, function f (r ) reduces to the step function and equation (57) reduces to equation (5). A similar distribution of the wakefield with d < 1 is observed in the 3D PIC simulation [28,31]. A typical trajectory of a trapped electron is very similar to the one calculated in the piecewise model (see figure 2(a)). The electron with a small impact parameter is not trapped, whereas one with a large impact parameter ρ ≈ R can be trapped in the bubble.
To study the electron dynamics in the bubble field analytically, we assume that the electron is located initially at the bubble border p = 0, y = R and ζ = 0 at t = 0. It follows from the initial condition that where ϕ(r = ∞) = ϕ(r = R) = −R 2 /8. As in the piecewise model, the most critical instant of time when the electron can leave the bubble is t = t m > 0 when dy/dt = p y = 0 and the electron excursion from the ζ -axis reaches its maximum. At this moment, we can write γ − V p x p x /(2γ 2 0 ) + 1/(2 p x ), where p x 1 is assumed, γ 0 = (1 − V 2 ) −1/2 1 is the gamma-factor of the bubble. The square of the distance of the electron from the bubble center (ζ = 0, y = 0) at t = t m is where we use equation (58).
We change the variables p y = R 2 P y , p x = R 2 P x , ζ = X R, y = Y R and t = Rs. As a result, equations (8)-(11) take the form dP x ds = − where P x = P y = 0 and dP y /ds = −1/4, dP The electron leaves the bubble if there exists at least one value of s m when the condition (64) is satisfied. We can take only P x (s 1 ) 1.1 since P x (s 1 ) < P x (s 2 ) < P x (s 3 ) < · · ·. As a result, we come to the condition for electron capture in the bubble [34] γ 0 R which is close to the condition obtained numerically in [28]. Now we estimate the effect of the impact parameter on the trapping condition. If ρ < R, then the electron first moves through the decelerating bubble field where ζ > 0. As a result, the electron gains negative momentum, , at ζ = 0 in contrast to the electron with ρ = R, which starts motion from position ζ = 0 with P = 0. The deceleration leads to reduction of the longitudinal momentum at s = s 1 for the electron with ρ < R as compared to the electron with ρ = R. We estimate this reduction, P x (s 1 , ρ = R) − P x (s 1 , ρ < R), as , that is close to the value calculated by numerical integration of equations (9) and (10). To calculate , we can use equation (58) and assume 1/R R − ρ R. As p y 0 at ζ = 0, y ρ then so that where ν = (1 − ρ 2 /R 2 )/4. Therefore and the condition for electron trapping (65) can be rewritten as follows: Therefore the trapping cross-section σ near the trapping threshold γ 0 2 −1/2 R takes the form It follows from equation (71) that the trapping cross-section decreases as γ 0 /R increases and the trapping stops at the threshold γ 0 2 −1/2 R.

Effect of bubble deformation and bubble field enhancement on self-injection
The developed model abstracts from some important effects observed in numerical simulations: bubble-shape deformation during laser pulse propagation [33] and bubble field enhancement at the bubble back due to electron sheet crossing [31,33]. It follows from 3D PIC simulations that the bubble back typically has a gamma-factor that can be much less than the gamma-factor of the bubble front [28,31]. When the bubble size increases, the effective gamma-factor of the bubble back goes down and vice versa: when the bubble size shrinks, the back of the bubble moves much faster than the head. In fact, the back of the bubble can even become superluminal for a while. The bubble-shape deformations can be induced by the nonlinear dynamics of the laser pulse [21], plasma electron loading in the bubble [35], laser-plasma instabilities [21], etc. For example, laser pulse dynamics can be very complex in the nonlinear regime and is accompanied by self-focusing, self-compression, nonlinear absorption, etching of the pulse front and other effects. The bubble inflation prior to and during self-injection is a typical phenomenon (see figure 1 and, for example, simulation results in [5,31]). We can take into account the bubble-shape deformation effect in our model using the bubble back gamma-factor in equation (65) as γ 0 since the electrons are trapped at the bubble back. Although the difference between the gamma-factors of the bubble back and the bubble front can be very large, we can include this effect in our model because the difference between the velocities of the bubble back and front is very small and much less than the bubble velocity (which is very close to the speed of light). We assume that the rate of the bubble elongation is dR/dt and the gamma-factor of the bubble front is γ fr . Then the gamma-factor of the bubble back can be deduced from the equation where it is assumed that γ fr 1 and γ bb 1. In the limit dR/dt 1/(2γ 2 fr ) (which is relevant to our simulation example) γ bb (2dR/dt) −1/2 . If the bubble velocity V = (V bb + V fr )/2 then the bubble gamma-factor is close to the gamma-factor of the bubble back γ 0 γ bb (2dR/dt) −1/2 γ fr . Then the condition for electron trapping equation (65) is reduced to the form If bubble expansion is taken into account the Hamiltonian (7) is no longer constant. The model with the Hamiltonian evolving because of bubble expansion is studied in reference [37].
It follows from the model that a sufficient condition for electron trapping can be associated with change in the Hamiltonian (7), H . If the bubble is expanding so fast that H of the electron is less than −1, then such an electron is trapped. Unfortunately it is not easy to express H through the laser-plasma parameters or even through the bubble parameters. It follows from this condition and numerical simulations performed in [37] that the threshold value of for electron trapping is about 0.009 for R = 5. This is in an excellent agreement with the threshold value = 0.008 obtained from equation (73). Further investigations are needed to develop a self-consistent theory for self-injection including the nonlinear dynamics of the laser pulse and the bubble. Another effect is that the accelerating and focusing fields can be stronger at the bubble back due to electron sheet crossing [31,33] than a uniform spherical cavity would provide. We can roughly estimate the field enhancement effect introducing the field enhancement factor k inside the bubble so that A x = −ϕ = kr 2 /8. Although the field is strengthened only at the bubble back, we assume that the field intensity increases everywhere inside the bubble. We can use this assumption because the trapped electrons spend most of the time in the bubble back and the additional enhancement of the field in the rest of the bubble does not significantly modify the electron trajectory. It follows from the scalability analysis of equations (60)-(63) that the coefficient on the right-hand side of equation (65) becomes √ k/2. It should be noted that extensive loading of the plasma electron in the bubble can lead to the reverse effect [35]. If the charge of the self-injected electron bunch becomes too large, it can reduce the field at the bubble back.
The considered effects can significantly affect electron self-injection. For the parameters of the simulations performed in [31], the minimal radius of the bubble with self-injection is reduced by a factor of 4 when the bubble field enhancement (∼4, see figure 3 in [31]) and the bubble back gamma-factor (∼9, see figure 1 in [31]) are taken into account in equation (65). The estimated radius is about 1.5 times more than observed in the simulations. The numerical approach is used in [28] to estimate the trapping cross-section. The cross-section was expressed through the electron sheath width, which cannot be derived in our analytical model. However, our model and approach used in [28] give similar results for the cross-section if we take into account the field enhancement (the enhancement factor was about 2, see figure 2(a) in [28]) and use the gamma factor for the bubble back. Further studies are needed to include the abovementioned effects in our model accurately.

Numerical simulations
We carry out a numerical simulation of laser-plasma interaction by using PIC VLPL3D code [20]. We use a simulation box of dimensions (60 × 88 × 88)λ, volume (x, y, z), Cartesian geometry, where λ = 2π/k 0 = 0.820 µm is the laser wavelength. The resolution in the transverse direction is k p0 y = 0.26, the resolution in the longitudinal direction is k p0 x = 0.05, where k p0 = ω p0 /c, k 0 /k p0 = 11.95 and ω 2 p0 ≡ 4π e 2 n 0 /m. We use eight macroparticles per cell. The laser and plasma parameters are the same as for figure 1 except that the polarization is now circular. Figure 5  We modify our code to obtain more detailed information on the trapped particles. The list of particle parameters includes (t, x, y, z, p x , p y , p z , i p , x 0 , y 0 , z 0 ), where t is the current time, i p is the unique label of the particle and x 0 , y 0 , z 0 are the initial coordinates of the particle. The 'initial' coordinates of a particle mean that the immobile particle is located at the laser pulse front. By postprocessing these data, we can follow the trajectory of the trapped electrons and determine the initial coordinates of the trapped electrons. To separate the trapped electrons, we introduce an energy threshold: we assume that an electron was trapped if its maximum energy exceeds 75 MeV. The threshold energy is chosen so that it is more than the maximal value of γ 0 mc 2 (in our case 15 MeV) and less than the maximal value of the accelerated electrons (300 MeV). The distribution of the initial position of the trapped electrons at t = 2000 is shown in figure 5(b). It can be seen from the figure that the volume of the initial position of the trapped electron increases in the x-direction because of the bubble inflation. The volume is axially symmetrical, which is different from the case of the linear polarization [38]. Figure 6 represents the distribution of the gamma-factor of the captured electrons as a function of x 0 for three time instants: t = 1000λ/c, 1500λ/c and 2000λ/c. We can conclude from the figure that the selfinjection and acceleration of electrons occur irregularly. It can be seen from figure 6 that the self-injection starts at x 0 600λ, stops at x 0 700λ and is restarted at x 0 820λ. The selfinjection is again absent for 1380λ < x 0 < 1720λ and restarted at x 0 1720λ. The distribution of the number of electrons trapped in the bubble as a function of their initial position x 0 is shown in figure 7(a). Figure 7(a) again demonstrates that the capture process is non-uniform and the number of trapped electrons changes along the laser path. The gamma-factor of the bubble γ 0 can be found from the dependence of the current longitudinal coordinates of trapped electrons on their initial coordinates. The electrons captured at earlier instants of time have greater current coordinate. Therefore the particles inside the electron bunch are moving from the back part to the front part. The velocity of the particles inside the electron bunch is determined by the slope of the approximation line. We assume that at the instance t = 0, the position of the bubble centre is x = 0. The bubble reaches the electron with the initial coordinate x 0 at the instance t 0 , and the electron starts its motion. As follows from our model, the transversal coordinate of the electron which will be trapped is y ≈ R; therefore, at the instance t 0 the coordinate of the bubble centre is equal to x 0 . Thus, where V (t ) is the bubble velocity. The coordinate of the electron at the instance t is given by where v e (t 0 , t ) is the projection on the x-axis of the velocity of the electron which starts motion at the instance t 0 . Differentiating equation (75) with respect to x 0 yields where the formula v e (t 0 , t 0 ) = 0 was used. It follows from our model (see section 6) that v e (t 0 , t ) depends only on the combination η = (t − t 0 )/R(t 0 ), so dv e dt 0 = dv e dη We assume that the velocity of the electron at the time instance t is much closer to the velocity of light than the bubble velocity: t). Under this assumption from equations (76) and (77), we have where we set v e (t 0 , t) = 1. It follows from the scaling presented in section 6 and numerical simulation of equations (60)-(63) that the characteristic time, in which v e becomes close to unity, is of the order R. Thus, the last term in equation (78) can be estimated as the ratio of R to the characteristic time of changing of the bubble radius. If this ratio is much smaller than 1/γ 2 0 (t 0 ), we obtain dx dx 0 ≈ − 1 2γ 2 0 (t 0 ) .
Therefore the gamma-factor of the bubble can be found from the slope of the curve x(x 0 ). It should be noted that γ 0 (−2dx/dx 0 ) −1/2 retrieved from our simulations largely determines the gamma-factor of the bubble back since the trapped electrons spent most of the time in the bubble back. It is seen from figures 7(a), (b) that the self-injection is strongly correlated with ratio γ 0 . It follows from the simulations (see figure 7(a)) that the number of trapped particles per laser period peaks in the regions 820λ/c < t < 1380λ/c and t > 1720λ/c, where γ 0 is minimal. Vice versa, for 1380λ/c < t < 1720λ/c when γ 0 > 2 −1/2 R the number of the trapped particles becomes negligible. The piecewise model and the bubble field model qualitatively agree with simulation results (see figure 7(b)) while the piecewise model implies a less restrictive condition for electron trapping than bubble field model, even for F = 1.

Conclusions
We study the self-injection physics in the framework of the piecewise and bubble field models. The dynamics of the trapped electron is similar in both models. However, the trapping condition in the piecewise model given by equation (53) is less restrictive than that in the bubble field model given by equation (65). Equation (53) reduces to the form 3.8γ 0 < R 4 for F R and to 3.8γ 0 < R 2 for F 1 that is less restrictive than the condition given by equation (65).
The difference can be caused by the strong attenuation of the bubble field near the axes y = 0 and ζ = 0 in the bubble field model.
There are some important effects that are not taken into account in the proposed models. Among the effects that can essentially modify the self-injection are the bubble-shape deformation during laser pulse propagation and the bubble field enhancement at the bubble back due to electron sheet crossing. 3D PIC simulations [5,28,31] demonstrate that these effects are presented for a wide range of parameters. It is shown in section 5 that these effects can strongly change the threshold for the electron trapping in the bubble. We add these effects in our models by making use of simple estimations and further investigations are needed to include the effects more accurately. Particularly, the Hamiltonian of the trapped is no longer the integral of the motion when the bubble deforms. This should be taken into account for a more accurate description of the electron self-injection in the evolving bubble [37].
The proposed model predicts that self-injection becomes less efficient for a large bubble velocity when the plasma density is low. Moreover, there is a density threshold for selfinjection that is different from the self-injection condition proposed in [31]. The evidence of such threshold has been observed in experiments [15]. The proposed model is in a manner consistent with the 1D model requiring that the plasma wave phase velocity has to be low enough for plasma wave breaking and electron trapping [21]. The model also predicts trapping enhancement due to the bubble expansion as the expansion sharply reduces the gamma factor of the bubble back. The model can also explain the electron self-injection scheme based on downward plasma density transition [18]. Such a transition leads to bubble elongation because of the decrease of plasma density. This sharply reduces the gamma-factor of the bubble back, thereby strongly enhancing self-injection. Recently, the strong impact of the wake phase velocity on electron trapping was also observed experimentally [39].
Moreover, the appropriate plasma profile can be used to produce monoenergetic electron bunches in the bubble regime. If the region of downward plasma density transition follows by a region with constant plasma density, the generation of monoenergetic electron bunches due to phase space rotation is possible that is similar to phase space rotation due to laser pulse widening [37]. In the region of downward plasma density transition a big amount of electrons can be trapped. The energy of electrons in the head of the bunch will be greater than the energy of the electrons in the tail of the bunch. However, in the region of constant plasma density, where electron self-injection can stop, the electrons in the bunch head feel weaker accelerating field than ones in the bunch tail (due to the increase of the accelerating field strength with the increase of the distance from the bubble center). Thus, if the length of constant density region chosen properly, resulting electron energy in the tail of the bunch can be equal to the electron energy in the head, and quasimonoenergetic electron bunch can be produced.
In conclusion, we present the model for electron self-injection in the relativistic plasma wakefield generated by a short laser pulse or by an electron beam. However, accurate measurements of the bubble velocity and the structure of the electromagnetic field in the bubble are needed for careful verification of the results obtained. Much effort is now denoted to theory and experiments to control the emission and the energy spread of the accelerated electron bunch. It follows from the obtained results that self-injection strongly depends on wake phase velocity that can be controlled, for example, by plasma density profiling [40]. Therefore the obtained results can be used to optimize the plasma-based acceleration.
The constant C depends only on the initial coordinates of the electron y 0 , ζ 0 and the parameters of the bubble. To calculate C, first we consider the case of y 0 = R and ζ 0 = 0, in that the electron moves only in a half-plane ζ < 0; hence, the electron motion in a cavity does not depend on ζ and γ 0 (see equations (12), (13) and (15)). It follows from equations (12) and (13) that the electron becomes ultrarelativistic in a time period of about 1/F and reaches the ζ -axis in a time period of about R. Thus, for F R 1, the electron motion can be treated as ultrarelativistic during all the time of motion, and in equations (12), (13) and (15) we can set γ ≈ p 2 x + p 2 y . The resulting system is the following: d p x dt = (1 + V ) where α and β are arbitrary constants, yields the same system for the primed variables, but with F = F/β and the initial condition y 0 = y 0 /α. The substitution (A.4) transforms C as follows: We calculate C numerically from equations (12) The results of the numerical calculation of C for several values of F are shown in figure 2.
We consider now the case ζ 0 = 0. The momentum components of an electron with initial coordinates ζ 0 and y 0 = (R 2 − ζ 2 0 ) 1/2 obey the following equation for ζ > 0 and y > 0: d p x d p y = 2 + 3 p y − 2 p x 2 p y + 1 . (A.7) The solution of this equation with the initial condition p x ( p y = 0) = 0 is p x = 3 4 p 2 y + p y p y + 1/2 .
During the motion of an electron in the first quadrant, p x and p y are negative and their absolute values increase in the course of time. However, it follows from equation (A.8) that p y remains more than −1/2. If Fζ 0 1 the electron becomes ultrarelativistic in a short time interval so that | p x | 1 and p y −1/2; therefore |v y | 1 and v x ≈ −1. It follows from equation (13) that the transversal force exerted on the electron is small until p x becomes equal to zero in the second quadrant (ζ < 0 and y > 0). Thus, we can assume that v y and, hence, transversal displacement of the electron is negligible until the time instant when p x = 0. At this instant of time p y 1. Then | p y | becomes much more than unity in a time interval that is much less than R. Therefore, in the second quadrant, where p x = 0, we can set y ≈ y 0 and p y ≈ 0 (see figure 4). This means that the subsequent electron motion is the same as for ρ = R, and C is also determined by equation (A.6). The results of numerical calculation of C as a function of ζ 0 = 0 for several values of F are shown in figure 2(b). It can be seen from figure 2(b) that C is approximately constant and determined by equation (A.6).