Simple model of the slingshot effect

We present a detailed quantitative description of the recently proposed"slingshot effect"[Fiore, Fedele, De angelis 2014]. Namely, we determine a broad range of conditions under which the impact of a very short and intense laser pulse normally onto a low-density plasma (or matter to be locally completely ionized into a plasma by the pulse) causes the expulsion of a bunch of surface electrons in the direction opposite to the one of propagation of the pulse, and the detailed, ready-for-experiments features of the expelled electrons (energy spectrum, collimation, etc). The effect is due to the combined actions of the ponderomotive force and the huge longitudinal field arising from charge separation. Our predictions are based on estimating 3D corrections to a simple, yet powerful plane magnetohydrodynamic (MHD) model where the equations to be solved are reduced to a system of Hamilton equations in one dimension (or a collection of) which become autonomous after the pulse has overcome the electrons. Experimental tests seem to be at hand. If confirmed by the latter, the effect would provide a new extraction and acceleration mechanism for electrons, alternative to traditional radio-frequency-based or Laser-Wake-Field ones.

We present a detailed quantitative description of the recently proposed "slingshot effect". Namely, we determine a broad range of conditions under which the impact of a very short and intense laser pulse normally onto a low-density plasma (or matter locally completely ionized into a plasma by the pulse) causes the expulsion of a bunch of surface electrons in the direction opposite to the one of propagation of the pulse, and the detailed, ready-for-experiments features of the expelled electrons (energy spectrum, collimation, etc). The effect is due to the combined actions of the ponderomotive force and the huge longitudinal field arising from charge separation. Our predictions are based on estimating 3D corrections to a simple, yet powerful plane 2-fluid magnetohydrodynamic (MHD) model where the equations to be solved are reduced to a system of Hamilton equations in one dimension (or a collection of) which become autonomous after the pulse has overcome the electrons. Experimental tests seem to be at hand. If confirmed by the latter, the effect would provide a new extraction and acceleration mechanism for electrons, alternative to traditional radio-frequency-based or Laser-Wake-Field ones.

I. INTRODUCTION AND SET-UP
Laser-driven Plasma-based Acceleration (LPA) mechanisms were first conceived by Tajima and Dawson in 1979 [1] and have been intensively studied since then. In particular, after the rapid development [2,3] of chirped pulse amplification laser technology -making available compact sources of intense, high-power, ultrashort laser pulses -the Laser Wake Field Acceleration (LWFA) mechanism [1,4,5] allows to generate extremely high acceleration gradients (>1GV/cm) by plasma waves involving huge charge density variations. Since 2004 experiments have shown that LWFA in the socalled bubble (or blowout) regime can produce electron bunches of high quality (i.e. very good collimation and small energy spread), energies of up to hundreds of MeVs [6][7][8] or more recently even GeVs [9,10]. This allows a revolution in acceleration techniques of charged particles, with a host of potential applications in research (particle physics, materials science, structural biology, etc.) as well as applications in medicine, optycs, etc.
In the LWFA and its variations the laser pulse travelling in the plasma leaves a wakefield of plasma waves behind; a bunch of electrons (either externally [11] or self injected [12]) can be accelerated "surfing" one of these plasma waves and exit the plasma sample just behind the pulse, in the same direction of propagation of the latter (forward expulsion). In Ref. [13] a new LPA mechanism, named slingshot effect, has been proposed, in which a bunch of electrons is expected to be accelerated and expelled backwards from a low-density plasma sample shortly after the impact of a suitable ultra-short and ultra-intense laser pulse in the form of a pancake normally onto the plasma (see fig. 1). The surface electrons (i.e. plasma electrons in a thin layer just beyond the vacuum-plasma interface) first are all displaced for-ward (with respect to the ions) by the ponderomotive force F p := −e( v c × B) z generated by the pulse, leaving a layer of ions completely depleted of electrons (here is the average over a period of the laser carrier wave, E, B are the electric and magnetic fields, v is the electron velocity, c is the speed of light,ẑ is the direction of propagation of the laser pulse); F p is positive (negative) while the modulating amplitude s of the pulse respectively grows (decreases). These electrons are then pulled back by the longitudinal electric force F z e = −eE z exerted by the ions and the other electrons, and leave the plasma. [In the meanwhile the pulse proceeds deeper into the plasma, generating a wakefield.] Tuning the electron density in the range where the plasma oscillation period T H [28] is about twice the pulse duration τ , we can make these electrons invert their motion when they are reached by the maximum of s , so that the negative part of F p (due to the subsequent decrease of s ) adds to F z e in accelerating them backwards; thus the total work W = τ 0 dt F p v z done by the ponderomotive force is maximal [29]. Provided the laser spot size R is sufficiently small a significant part of the expelled electrons will have enough energy to win the attraction by ions and escape to infinity.
Very short τ 's and huge nonlinearities make approximation schemes based on Fourier analysis and related methods (slowly varying amplitude approximation, frequency-dependent refractive indices,...) unconvenient. On the contrary, in the relevant space-time region a MHD description of the impact is self-consistent, simple and predictive (collisions are negligible, and recourse to kinetic theory is not needed). Here we develop and improve the 2-fluid MHD approach introduced in [13,14] and apply it to determine a broad range of conditions enabling the effect, as well as detailed quantitative predictions about it (a brief summary is given in [15,16]). In section II we study the plane problem (R = ∞) and show that for sufficiently low density and small times (after the impact) we can neglect the radiative corrections [backreaction of the plasma on the electromagnetic (EM) field (3)] and determine the motion of the surface electrons in the bulk by (numerically) solving a single system of two coupled first order ordinary differential equations of Hamiltonian form, if the initial density n 0 is step-shaped, or a collection of such systems, otherwise; the role of 'time' is played by the light-like coordinate ξ = ct − z. The rough model of [13] considered only step-shaped n 0 and was based on neglecting: F z e during the forward motion, F p during the backward motion of the electrons; the estimates could be considered reliable only for very low, unrealistic n 0 . Here n 0 needs no longer to be so low, nor step-shaped, as in [13], because we take F z e , F p in due account during the whole motion of the electrons. In section III we heuristically modify the potential energy outside the bulk to account for finite R and determine a R-range such that the motion of the surface electrons coming from some inner cylinder ρ 2 ≤ r 2 < R 2 be (by causality) well approximated by the solution of the correspondingly modified Hamilton equations; we then find which electrons indeed escape to infinity and estimate in detail their final energy spectrum, collimation, total number, charge and energy. To be specific, in section IV we specialize predictions to potential experiments at the FLAME facility (LNF, Frascati) or the ILIL laboratory (INO-CNR, Pisa). We welcome 3D simulations and experiments checking these predictions; the experimental conditions are at hand in many laboratories. In section V we discuss the results, the conditions for their validity and draw the conclusions.
As a context remark, we recall that relatively simple 2-fluid magnetohydrodynamic models can be used also to describe the complicated physics of the impact of very intense and short laser pulses on overdense solid targets. If the density gradient of the target is sufficiently steep, the massive displacement of electrons (induced by the ponderomotive force) with respect to ions (named snowplow in [18,19]) produces a longitudinal electric force which may accelerate also protons or other light ions, either backward or forward, by the socalled Skin-Layer Ponderomotive Acceleration [17] or Relativistically Induced Transparency Acceleration [18,19] mechanisms.

The 2-fluid magnetohydrodynamic framework
The set-up is as follows. We assume that the plasma is initially neutral, unmagnetized and at rest with electron (and proton) density equal to zero in the region z < 0. We describe the plasma as consisting of a static background fluid of ions (the motion of ions can be neglected during the short time interval in which the effect occurs) and a fully relativistic collisionless fluid of electrons, with the "plasma + EM field" system fulfilling the Lorentz-Maxwell and the continuity equations. We show a posteriori that such a MHD treatment is self-consistent in the spacetime region of interest. We denote as x e (t, X) the position at time t of the electrons' fluid element initially located at X ≡ (X, Y, Z), and for each fixed t as X e (t, x) the inverse map [x ≡ (x, y, z)]. For brevity, we refer: to such a fluid element as to the "X electrons"; to the fluid elements with arbitrary X, Y and specified Z, or with X in a specified region Ω, respectively as the "Z electrons" or the "Ω electrons". We denote as m, n e , v e the electron mass, Eulerian density and velocity and often use the dimensionless fields β e ≡ v e /c, u e ≡ p e /mc = β e / 1−β 2 e , γ e ≡ 1/ 1−β 2 e = 1+u 2 e . The equations of motion are in CGS units (d/dt ≡ ∂ t +v e ·∇ x is the electrons' material derivative) and the initial conditions are p e (0, X) = 0, x e (0, X) = X for Z ≥ 0. The Lagrangian fields depend on t, X, rather than on t, x, and are distinguished by a tilde, e.g. n e (t,X) = n e [t,x e (t,X)]. The continuity equation dn e /dt + n e ∇ x · v e = 0 follows from the local conservation of the number of electrons, which amounts toñ We assume that n 0 is independent of X, Y and, as said, vanishes if Z < 0; also as a warm-up to more general Z-dependence, we start by studying the case that it is constant in the region Z ≥ 0: n 0 (Z) = n 0 θ(Z), where θ is the Heaviside step function. We consider a purely transverse EM pulse in the form of a pancake with cylindrical symmetry around the z-axis, propagating in the positiveẑ direction and hitting the plasma surface z = 0 at t = 0. We schematize the pulse as a free plane pulse multiplied by a "cutoff" function χ R (ρ) which is approximately equal to 1 for ρ ≡ x 2 +y 2 ≤ R and rapidly goes to zero for ρ > R (with some finite radius R, see fig. 1

II. PLANE WAVE IDEALIZATION
In the plane problem (R = ∞) the invertibility of x e : X → x for all fixed t amounts to z e (t, Z) being strictly increasing with respect to Z for all t. Eq. (2) becomes ñ e ∂ze ∂Z (t,Z) = n 0 (Z) ⇔ n e (t,z) = n 0 [Z e (t,z)] ∂Ze ∂z (t,z).
(4) Regarding ions as immobile, the Maxwell equations imply [14] that the longitudinal component of the electric field is related to N (Z) ≡ Z 0 dZ n 0 (Z ) (the number of electrons per unit surface in the layer 0 ≤ Z ≤ Z) by We partially fix the gauge [14] imposing that the transverse (with respect toẑ) vector potential itself is independent of x, y, and hence is the physical observ- As known, the transverse component of the Lorentz equation (1) 1 implies p ⊥ e − e c A ⊥ =const on the trajectory of each electron; this is zero at t = 0, hence p ⊥ e = mcu ⊥ e = eA ⊥ /c. Hence u ⊥ e is determined in terms of A ⊥ . As in [14], we introduce the positive-definite field which we name electron s-factor. u z e , γ e , β ⊥ e , β z e are recovered from u ⊥ e , s e through the formulae (44) of [14]: Remarkably, all of (7) are rational functions of u ⊥ e , s e (no square roots appear). Moreover, fast oscillations of u ⊥ e affect γ e , u z e but not s e [see the comments after (15)]. For these reasons it is convenient to use u ⊥ e , s e instead of u ⊥ e , u z e as independent unknowns. The evolution equation of s e (difference of the ones of γ e , u z e ; the former is the scalar product of (1) 1 with p e /γ e m 2 c 2 ) reads The Maxwell equation for A ⊥ takes the form ( Using the Green function of the D'Alembertian ∂ 2 0 −∂ 2 z , abbreviating x ≡ (t, z), these equations can be equivalently reformulated as the integral equation (42) of [14] A The past, future causal cones D x , T , the supports of A ⊥ , n 0 (z), and their intersections are shown in fig. 2. For t < 0 D x ∩T is empty, and the right-hand side of (9) 1 is zero, as it must be. Below we shall analyze the consequences of neglecting it also for small t, and determine the range of validity of such an approximation.

II.1. Motion of the electrons
F z e (t,Z) ≡ F z e [z e (t,Z),Z] is the longitudinal electric force acting on the Z electrons at time t; it is conservative, as it depends on t only through z e (t,Z). The approxi- , and the last term of (8) vanishes. Replacing (5) in the Lagrangian version of (8), we find for each Z ≥ 0 the equationγ e ∂ tse = −s e F z e /mc. The initial condition is s e (0, Z) ≡ 1. The other equation to be solved is (1) 2 with the initial condition x e (0, X) = X. By (7) one is thus led to the Cauchy problems (parametrized by Z ≥ 0) x ⊥ e (t,X) is obtained from the solutions of (11-12) using (1), (7): For all fixed Z the map t →ξ(t,Z) ≡ ct − z e (t,Z) is invertible, because the speed of electrons is always smaller than c. We can simplify (11) by the change of variables (t,Z) → ξ, Z , making the argument of v an independent variable. Denoting the dependence on ξ, Z by a caret [e.g.ŝ(ξ,Z) =s e (t,Z)] and introducing the displacement from the initial position∆(ξ,Z) ≡ẑ e (ξ,Z)−Z, we find ∂ ξ = (γ e /cs e )∂ t , and (11) becomeŝ (the prime means differentiation with respect to ξ). For ξ ≤ 0 v(ξ) ≡ 0,∆,ŝ remain constant, and we can replace the initial conditions∆(−Z,Z) = 0,ŝ(−Z,Z) = 1 bŷ An alternative derivation of (14-15) with a deeper insight on the role of the s-factor is given in [20]. In the zero density limit N (Z) ≡ 0,ŝ ≡ 1, (14-15) is integrable, and all unknowns are determined explicitly from ⊥ [14,21]).
Both decrease for ξ > ξ r (Z), untilŝ becomes so small, and the right-hand side of (14) 1 so large, that first∆, and thenŝ−1, are forced to abruptly grow again to positive values. This preventsŝ to vanish anywhere, consistently with (6). In ξ-intervals where v(ξ) ≡ v c ≡const, H is conserved, and all trajectories P (ξ; Z) in phase space (paths) are level curves H(∆, s; Z) = h(Z), above the line s = 0, integrable by quadrature [22]. For Z = 0 the paths are unbounded with∆(ξ, 0) → −∞ as ξ → ∞. For Z > 0 the paths are cycles around the only critical point The function ξ ex (Z) is strictly increasing if ∂ Zẑe > 0. For any family P (ξ; Z) of solutions of (14-15) let (note thatŶ z =∆). The so definedû,γ,x e are the solutions -expressed as functions of ξ, X -of all equations and initial conditions [31]. Note thatx e ,t can be obtained also solving the system of functional equations [by (19) the second is actually equivalent to the zcomponent of the third] with respect to t, x. Clearlŷ Ξ(ξ,Z) is strictly increasing and invertible with respect to ξ for all fixed Z. Solving (20) with respect to ξ, x (resp. ξ, X) as functions of t, X (resp. of t, x) and replacing the results inû,γ,ŝ, ... one obtains the solutions in the Lagrangian (resp. Eulerian) description: in particular one finds (generalizing [14]) Indeed, it is straightforward to check that z e (t, Z),s e (t, Z) is the solution of (11)(12) and From (17), (18), (19) 5 , the times of maximal penetration and of expulsion of the Z electrons arē Deriving (21) and the identity y ≡Ξ Ξ −1 (y,Z),Z we obtain a few useful relations, e.g.
By (23), ∂ Zẑe ≡ 1+∂ Z∆ > 0 is thus a necessary and sufficient condition for the invertibility of the maps z e : Z → z, x e : X → x (at fixed t), justifying the hydrodynamic description of the plasma adopted so far and the presence of the inverse function Z e (t,z) in (21). Finally, from (4), (23) we find also We can test the range of validity of the approximation A ⊥ (t, z) = α ⊥ (ct − z) by showing that the latter makes the modulus of the right-hand side of (9) much smaller than defined below), or equivalently [multiplying by e/mc 2 and using (24)] ; actually, it suffices to check this inequality on the worldlines of the expelled electrons.

II.2. Auxiliary problem: constant initial density
As a simplest illustration of the approach, and for later application to a step-shaped initial density, we first consider the case that n 0 (Z) = n 0 .
Then F z e is the force of a harmonic oscillator (with equilibrium at z e = Z) F z e (z e ,Z) = −4πn 0 e 2 [z e − Z] = −4πn 0 e 2 ∆; the Z-dependence disappears completely in (14)(15), which reduces to the auxiliary Cauchy problem where M ≡ 4πe 2 n 0 /mc 2 . The potential energy in (16) takes the form U(∆, Z) ≡ M ∆ 2 /2. Problem (26), and hence also its solution ∆(ξ), s(ξ) , the value of the energy as a function of ξ and the functions defined in (19), are Z-independent. It follows ∂ Z∆ ≡ 0 and by (23) the automatic invertibility of z e (t,Z); moreover, the inverse function Z e (t,z) has the closed form [here Ξ(ξ) ≡ ξ+∆(ξ)], what makes the solutions (21) of the system of functional equations (20), as well as those of (1), completely explicit in terms of Ξ and the inverse Ξ −1 only. As a consequence, all Eulerian fields depend on t, z only through ct−z (i.e. evolve as travelling-waves). In fig.  3-left we plot some solution of (26). If v(ξ) ≡ v c ≡const all paths P (ξ;Z) are cycles around C ( fig. 3-right), corresponding to periodic solutions. Within the bulk electron trajectories for slowly modulated laser pulse like the ones considered in section IV are tipically as plotted in fig. 8; in average they have no transverse drift, but a longitudinal forward/backward one. Fig. 4 shows a couple of corresponding charge density plots.

III. 3-DIMENSIONAL EFFECTS
We now discuss the effects of the finiteness of R. For brevity, for any nonnegative r, L we shall denote as C r the infinite cylinder of equation ρ ≤ r, as C L r the cylinder of equations ρ ≤ r, 0 ≤ z ≤ L. The ponderomotive force of the pulse will boost forward (as in fig. 8) only the small-Z electrons located within (or nearby) C R . These Forward Boosted Electrons (FBE) will be thus completely expelled out of a cylinder which will reach its maximal extension C ζ R around the timet(0) of maximal longitudinal penetration ζ ≡∆[ξ(0),0] of the Z = 0 electrons. The displaced charges modify E. By causality (see appendix A), for x near the z axis E(t, x) is the same as in the plane wave case for t t (0)+R/c, and smaller afterwards. We choose n 0 , R so that they fulfill (28) and condition (25) for all x = (t, x) such that t t (0)+R/c; heret ≡t(0), t ex ≡ t ex (0) are the times of maximal penetration and of expulsion from the bulk of the Z = 0 electrons [see (22)]. [As n 0 grows from zero the righthand side of (25) does as well, whereast(0), t ex (0)−t(0) decrease]. In appendix A we show that conditions (28) respectively ensure that these FBE, at least within an inner cylinder ρ ≤ r ≤ R: i) move approximately as in section II until their expulsion; ii) are expelled before Lateral Electrons (LE), which are initially located outside the surface of C R and are attracted towards the z-axis (see fig. 1.3), obstruct their way out. For the validity of our model we must a posteriori check also that the expelled C r electrons remain in C R , i.e. their transverse oscillations ∆x e are R.
In the plane model the expelled Z > 0 electrons cannot escape to z → −∞ because are decelerated by the constant electric force F z e (t, Z) = 4πe 2 N (Z) > 0, see (10). The real electric force F zr e > 0 acting on the C r electrons after expulsion is generated by charges localized in C R ; hence F zr e ∝ 1/z 2 e as z e → −∞, and the escape of expelled electrons is no more excluded a priori. Moreover, since F z e (t,0) = 0, it should be also F zr e (t,0) = 0, allowing the escape of the Z = 0 electrons; by continuity there will  fig. 7 after about 25 fs (up) and 70 fs (down), i.e. resp. before and after the maximal penetration timet(0) = 51fs; in the latter picture the electrons travelling backwards make light yellow-blue striped the region between the yellow and the white-blue striped ones. exist some positive Z M ≤ Z b such that the C Z M r electrons escape to infinity. We stick to estimate F zr e on the z-axis electrons; we assume that after the pulse has overcome them, they move along the z-axis. Actually this will be justified below ifû ⊥ (l) 0, which in turn holds if, as usual, l λ [see the comments after (33)]. In fig. 5 a) we schematically depict the charge distribution expected shortly after the expulsion. The light blue area is occupied only by the X ∈ C Z M r electrons. The orange area is positively charged due to an excess of ions. For any Z-electrons moving along the z-axis consider the surfaces S 0 , S 1 , S 2 occupied at time t by the X ∈ C r electrons respectively having Z = 0, Z, Z 2 (Z), where Z 2 (Z) is defined by the condition N (Z 2 ) = 2 N (Z), which ensures that the electron charges contained between S 0 , S 1 and S 1 , S 2 are equal (in the figure S 0 , S 1 , S 2 are respectively represented by the left border of the blue area, the dashed line and the solid line). The longitudinal electric force F zr e acting at time t on this Z-electron is nonnegative and can be decomposed and bound as follows [13]: Here E z − (t, Z) stands for the part of the longitudinal electric field generated by the electrons between S 0 , S 2 ; since those between S 0 , S 1 have by construction the same charge as those between S 1 , S 2 , but are more dispersed, generated by the ions and the remaining electrons (at the right of S 2 ) will be smaller than the force F z er generated by the charge distribution of fig. 5 b), where the remaining electrons are located farther from (0, 0, z e ) (in their initial positions X , not in the ones at t) and hence generate a smaller repulsive force. This explains the second inequality in the equation. In appendix B we show that for z e ≡ Z +∆ ≤ 0 Commendably, F z er is conservative, nonnegative and goes to zero as ∆ → −∞, while it reduces to zero for Z = 0 and to 4πe 2 N(Z) as r → ∞, as F z e in (10) 3 ; it becomes a function of t (resp. ξ) through ∆(t, Z) [resp.∆(ξ,Z)] only. We therefore modify the dynamics outside the bulk replacing F z e by F z er , or equivalently U by U r in (16), where U r is continuous and equals U for z e ≡ Z +∆ ≥ 0, and the potential energy (B2) associated to F z er for z e ≡ Z+∆ ≤ 0; there U r is a decreasing function of ∆ with finite left asymptotes (B3). We will thus underestimate the final energy of the electrons, because F z er is larger than the real electric force F zr e decelerating the electrons outside the bulk; this makes our estimates safer. In fig. 6 we plot suitably rescaled U and U r for n 0 (Z) = n 0 θ(Z). After the pulse is passed we can compute γ as a function of ∆, Z using energy conservation mc 2 γ +U r (∆,Z) =const. For the expelled electrons the final relativistic factor γ f (Z) ≡ γ e (∆ = −∞,Z) is the decreasing function (B4). The maximum of γ f (Z) is γ M ≡ γ f (0). Let Z M ≤ Z b be the Z fulfilling γ f (Z) = 1. The estimated total number N e , electric charge (in absolute value) Q, and kinetic energy E of the X ∈ C Z M r escaped electrons are thus The number of escaped X ∈ C Z M r electrons with Z ≤ Z ≤ Z +dZ is estimated as πr 2 n 0 (Z)dZ, that with relativistic factor between γ and γ + dγ is estimated as dN = πr 2 [ n 0 (Z)/|dγ f /dZ|] Z=Ẑ(γ) dγ, whereẐ(γ) is the inverse of γ f (Z) (a strictly decreasing function, see appendix B). Hence the fraction of escaped electrons with final relativistic factor between γ and γ+dγ is estimated determines the associated energy spectrum. As α ⊥ (ξ) = α ⊥ (l) if ξ ≥ l, by (7) the final transverse deviation of the escaped electrons will be where u ⊥ f ≡û ⊥ (l). This is an increasing function of Z, because γ f (Z) is decreasing. If λ l then u ⊥ f 0 (see next section), and (33) is negligible unless Z Z M .

III.1. Step-shaped initial density
If n 0 (Z) = n 0 θ(Z) then N (Z) = n 0 θ(Z)Z, and for Z ≥ 0 Since the first expression is the same as in the case n 0 (Z) = n 0 , the motion of the Z-electron will be as in subsection II.2 until ξ = ξ ex (Z). The second expression goes to the constant force 4πn 0 e 2 Z as r → ∞, as expected. The motion for ξ > ξ ex (Z) will be studied in detail in [22]; we plot the graphs of a typical solution (until the expulsion) in fig. 7 and a few corresponding electron trajectories in Fig. 8. We can readily understand that it will be ∂ Zẑe (ξ,Z)> 0 for all ξ and 0 ≤ Z ≤ Z M , since this holds for ξ ≤ ξ ex (Z) [by the comments following (26)], and both ξ ex (Z) and the decelerating force F z er (∆,Z) (outside the bulk) increase with Z, while the speed of exit from the bulk decreases with Z, whence the distance between

IV. NUMERICAL RESULTS
We assume for simplicity that the pulse is a slowly modulated sinusoidal function linearly polarized in the x direction: ⊥ (ξ) = s (ξ)x cos kξ, the modulating amplitude s (ξ) ≥ 0 is nonzero only for 0 < ξ < l, and slowly varies on the scale of the period λ ≡ 2π/k l, i.e. λ| s | | s | on the support of s . Integrating by parts we find α ⊥ (ξ) =x s (ξ)(sin kξ)/k +O(1/k 2 ) [20] and, in terms of the rescaled amplitude w(ξ) ≡ e s (ξ)/kmc 2 , where a b means a = b+O(1/k 2 ). Note that, as s (ξ) = 0 for ξ ≥ l, this implies u ⊥ f =û ⊥ (l) 0, as anticipated. If we approximate as χ R (ρ) ≡ θ(R−ρ) the cutoff function in (3), the average pulse intensity on its support is I = c E/πR 2 l. Here E is the EM energy carried by the pulse, High power lasers produce pulses where λ ∼ 1µm and s is approximately gaussian, s (ξ) ∝ exp −(ξ −ξ 0 ) 2 /2σ ; σ is related to the fwhm (full width at half maximum) l of 2 s by σ = l 2 /4 ln 2. If initially matter is composed of atoms then s (ct−z) can be considered zero where it is under the ionization threshold, because the pulse has not converted matter into a plasma yet. Hence we adopt as a modulating amplitude s (ξ) the cut-off Gaussian where U i is the first ionization potential (for Helium U i 24eV ); the formula for b 2 g follows replacing the Ansatz (39) 1 in (38) [neglecting the tails left out by the cutoff θ(ξ)θ(l − ξ)]. Numerical computations are easier if we adopt [14] as s (ξ) the following cut-off polynomial: b p , l p are determined by the requirement to lead to the same fwhm and E: b 2 p = 5040 E/R 2 l p and l p = 5l /2. We now present the results of extensive numerical simulations based on the experimental parameters available already now at the FLAME facility [23] or in the near future at the ILIL facility [32]: l 7.5µm (implying l p = 18.75µm), λ 0.8µm (implying kl p = 2πl p /λ 147), E = 5J, and R tunable by focalization in the range 10 −4 ÷ 1 cm. We model the electron density: first as the step-shaped one n 0 (Z) = n 0 θ(Z) (this allows analytical derivation of more results); then as a function smoothly increasing from zero to the asymptotic value n 0 , with substantial variation in the interval 0 ≤ Z ≤ L ≡ 20µm (as motivated by experiments, see section V), more precisely n 0 (Z) = n 0 θ(Z) tanh(Z/L). We have numerically solved the corresponding systems (14)(15) and proceeded as in section III, for R = 16, 15, 8, 4, 2, 1µm [resp. leading to average intensities I/10 19 (W/cm 2 ) 1, 1.1, 4, 16, 64, 255], n 0 in the range 10 17 cm −3 ≤ n 0 ≤ 3 × 10 20 cm −3 and Z ≤ Z M ; all results follow from these solutions.
In fig. 9-left we plot the maximal final relativistic factor γ M of the expelled electrons as a function of n 0 , with the above values of I and n 0 (Z) = n 0 θ(Z); each graph stops where n 0 becomes too large for conditions (25), (28) 1 , or (29) to be fulfilled and is red where condition (28) 2 is no more fulfilled. The latter prevents collisions with the LE and becomes superfluous if the target is a solid cylinder of radius R (since then there are no LE) [33]; the I = 64, 255 × 10 19 W/cm 2 graphs are plot green for densities corresponding to the lightest solids (aerogels) available today. As expected [13]: 1) as n 0 → 0 γ M − 1 ∝ n 0 I 2 ; 2) each graph γ M (n 0 ; I) has a unique maximum γ M M (I) ≡ γ M (n 0M ; I) at n 0M ∼n 0 , wheren 0 is the density makingξ(0) = l/2, namely such that the Z = 0 electrons reach the maximal penetration ζ =∆(l/2, 0) when they are reached by the pulse maximum. The dependence of γ M on n 0 is anyway rather slow. The striking γ M M (I) ∝ I behaviour shown in fig.  9 up-center hints at scaling laws and will be discussed elsewhere. In figures 10 we plot sample spectra ν(γ) for I/10 19 (W/cm 2 ) 1, 4, 16, 64 and n 0 compatible with (25), (28), (29). In table I we report our main predictions for the same I (equivalently, R) and n 0 . The final energies of the expelled electrons range from few to about 15 MeV. The spectra (energy distributions) are rather flat for the step-shaped densities, albeit they become more peaked near γ M as n 0 grows; if n 0 (Z) grows smoothly from zero to about the asymptotic value n 0 in the interval 0 ≤ Z ≤ L ∼ 20µm, they can be made much better  Table I. merical study of the mapẑ e (·, ξ) : Z → z shows that this holds true also in the continuous density case. Numerical computations show that (25) is fulfilled at least on the Z ≤ Z M electrons' worldlines, even with the highest densities considered here (see e.g. fig. 9 right). Finally, the data in table I show that (28), (29) are fulfilled.
If we choose s (ξ) as the cut-off gaussian, instead of the cut-off polynomial, convergence of numerical computations is slower, but the outcomes do not differ significantly. Sample computations show that choices of other continuous n 0 (Z) lead to similar results, provided the function n 0 (Z) is increasing and significantly approaches the asymptotic value n 0 in the interval 0 ≤ Z ≤ L ∼ 20µm.

V. DISCUSSION, FINAL REMARKS, CONCLUSIONS
These results show that indeed the slingshot effect is a promising acceleration mechanism of electrons, in that it extracts from the targets highly collimated bunches of electrons with spectra which can be made peaked around the maximum energies by adjusting R, n 0 ; with laser pulses of a few joules and duration of few tens of femtoseconds (as available today in many laboratories) we find that the latter range up to about ten MeV (it would increase with more energetic pulses). The spectra (distributions of electrons as functions of the final relativistic factor γ f ), their dependence on the electron density and pulse intensity, the collimation and the backward pulse energy E 5J, wavelength λ 0.8µm, fwhm l 7.5µm, spot radius R 1 ÷ 16µm p PoP14 p PoP14 p g cp cg cp cg cp cg cp cg pulse spot radius R (µm) 15  In the 'p,g' columns the initial electron densities are stepshaped, n0(Z) = n0 θ(Z), and the amplitudes are resp. of the gaussian, polynomial forms (39), (40); in the 'PoP14' columns we report results computed in [13] with poorer approximation. In the 'cg,cp' columns the initial electron densities are the continuous ones n0(Z) = n0 θ(Z) tanh(Z/L) with L = 20µm, and the amplitudes are resp. of the forms (39), (40).
direction of expulsion in principle allow to discriminate the slingshot effect from the LWFA or other acceleration mechanisms. In table I and fig. 10 we have reported detailed quantitative predictions of the main features of the effect for some possible choices of parameters in experiments at the present FLAME, the future upgraded ILIL facilities, or similar laboratories. Low density gases or aerogels (the lightest solids available today) are targets with appropriate electron densities. The steepest z-oriented density gradient of a gas sample isolated in vacuum is attained just outside a nozzle expelling a supersonic gas jet in the xy plane; across the lateral border of the jet the density may vary from about zero to almost the asymptotic value n 0 in about L ∼ 20µm [23]. Hence if we choose a supersonic helium jet as the laser pulse target the initial electron density is reasonably approximated by the choice n 0 (Z) = n 0 θ(Z) tanh(Z/L), and the predictions of table I, fig. 10 (a-d) are reliable. By the way, the values of n 0 considered in table I are considerably higher than in typical LWFA experiments.
Step-shaped n 0 (Z) are unrealistic approximations of densities of gas samples, but reasonable ones of solids (for which ∆Z λ), provided n 0 exceeds 48×10 18 cm −3 , which is the electron density of aerographene (the lightest aerogel so far: mass density=0.00016 g/cm 3 ). Silica areogels, with a wide range of densities from 0.7 to 0.001g/cm 3 , electron densities of the order of 10 20 /cm −3 and porosity from 50 nm down to 2 nm in diameter (i.e. much smaller than λ) have been produced and extensively studied [24,25]. Therefore the results of the last two 'p,g' columns of table I [and the corresponding spectra, figure 10 (e)] are applicable to aerogels, while those of the first four are presently only of academic interest.
The quantitative predictions of our model are based on a rather rigorous plane-wave, 2-fluid magnetohydrodynamic model [34] and simple, but heuristic approxima-tions for the 3D corrections, which certainly affect their liability. We welcome numerical 3D simulations (particlein-cell ones, etc.) to improve the latter. Experimental tests are easily feasible with the equipments presently available in many laboratories. We welcome experiments testing the effect. approximate as a potential leading to (3). It is easy to show that D + (D 0 R ) is the union of three regions, resp. of equations: a. z ≥ ct ≥ 0; b. ct ≥ z ≥ 0 and ρ+ √ c 2 t 2 −z 2 ≤ R; c. t ≥ 0, z < 0 and ρ+ct ≤ R (see fig. 11). In D + (D 0 R ) the two solutions coincide, in particular a "real" electron worldline x e (t, X) remains equal to the plane solution worldline x e (t, X) as long as x e (t, X) ∈ D + (D 0 R ). By continuity, we expect that the two solutions remain close to each other also in a neighbourhood of D + (D 0 R ). This is confirmed by estimates [13] involving the retarded electromagnetic potential (in the Lorentz gauge ∂ · A = 0) i.e. the general solution of the Maxwell equation A µ = 4πj µ with a current j µ (t, x) vanishing for t < 0; here t r (t,x−x ) ≡ t−|x−x |/c, A µ f (t,x) fulfills A µ f = 0 (determining the t → −∞ behaviour), and E = −1 c ∂ t A−∇ A 0 , B = ∇ × A. Since the formation of C ζ R is completed at t =t(0), and the 'information' [encoded in (A1)] about the finite radius of C ζ R takes a time R/c to go from the lateral surface ρ = R to the z-axis, then if eq. (28) 1 is fulfilled the X = 0 electrons (red worldline in fig. 11) move approximately as in section II until the expulsion. Similarly, the Z 0, ρ r electrons (yellow worldlines in fig.  11) move approximately as in section II untilt+(R−r)/c, i.e get the main backward boost (acceleration is maximal aroundt). Eq. (28) 2 is equivalent to t ex l/c; ⇒ r R; or 0 < (t ex −l/c)v ρ a < R ⇒ r R − (t ex −l/c)v ρ a > 0. If the left-hand side of the first line is fulfilled the surface electrons are expelled while the laser pulse is still entering the bulk and thus producing an outward force that keeps the LE out of C ζ R . Otherwise, the left-hand side of the second line ensures that the distance inward travelled by the most dangerous LE (the Z = 0 ones) after the pulse has completely entered the bulk is less than R; v ρ a stands for the average ρ-component of the velocity of these LE. By geometric reasons v ρ a < v z a ≡ average z-component of the X = 0 electrons velocity in their backward trip within the bulk; our rough estimate v ρ a v z a /2 = ζ/(t ex −t)2 gives (28) 2 . Eq. (28) is thus explained.

Appendix B: Finite R energies
Using cylindrical coordinates (y, ρ, ϕ) for X , one easily finds that for z e ≡ Z+∆ ≤ 0 the electric force generated by the static charge distribution of fig. 5 b), the associated potential energy mc 2 U r and the left asymptotes of U r are The last equality holds only ifẑ e (l,Z) ≥ 0, i.e. l ≤ ξ ex (Z); the right-hand side is the electrons' energy when expelled one in vacuum; the backward acceleration takes place afterwards and is due only to F z e , hence the final energy is smaller. Whereas if τ TH -which was the standard situation in laboratories until a couple of decades agothen Fpv z oscillates many times about 0, W 0, and the backward acceleration is washed out.
[30] Albeit the pump (3) violates the Maxwell equations (due to the ρ-dependence), we adopt it as for our purposes it is essentially equivalent to one that fulfills the Maxwell equations and at the time of impact coincides with it in the main part of its support, while rapidly decaying out-side (this and similar schematizations, e.g. the paraxial one, are currently used in the literature).
[33] We thank L. Gizzi for this remark.
[34] There is no need of a recourse to kinetic theory taking collisions into account, e.g. by BGK [26] equations or effective linear inheritance relations [27].