Ion-channel laser growth rate and beam quality requirements

In this paper, we determine the growth rate of the exponential radiation amplification in the ion-channel laser, where a relativistic electron beam wiggles in a focusing ion channel that can be created in a wakefield accelerator. For the first time the radiation diffraction, which can limit the amplification, is taken into account. The electron beam quality requirements to obtain this amplification are also presented. It is shown that both the beam energy and wiggler parameter spreads should be limited. Two-dimensional and three-dimensional particle-in-cell simulations of the self-consistent ion-channel laser confirm our theoretical predictions.

The ion-channel laser (ICL), initially proposed by Whittum, Sessler & Dawson (1990), relies on the injection of a relativistic electron beam in an ion channel (IC) to create a coherent and highly amplified radiation source. Such an IC can be produced in a plasma-based wakefield accelerator in the blowout or bubble regime (Faure et al. 2004;Geddes et al. 2004;Mangles et al. 2004): while propagating in a plasma, a laser pulse or a particle beam pushes the electrons off-axis and leaves an IC in its wake. The fields generated in the IC provide a focusing force for the relativistic electrons on-axis. The resulting wiggling motion of the electron along the IC axis then leads to the emission of the so-called betatron radiation (Esarey et al. 2002;Rousse et al. 2004). For appropriate conditions, betatron radiation can interact with the electron beam and bunch it at the radiation wavelength, allowing for the exponential amplification of the emitted radiation, as in a conventional free-electron laser (FEL). One of the most important advantages of the ICL are the strong fields generated in the plasma, which can lead to amplification in the UV to X-ray range, with very high brightness within much shorter distances than those obtained in the † Email address for correspondence: xavier.davoine@cea.fr X. Davoine, F. Fiúza, R. A. Fonseca, W. B. Mori and L. O. Silva conventional FEL sources. Previous works analysed the ICL gain length and the associated Pierce parameter assessments (Chen, Katsouleas & Dawson 1990;Whittum et al. 1990;Whittum 1992;Liu, Tripathi & Kumar 2007;Ersfeld et al. 2014).
In order to take full advantage of this scheme, it is critical to correctly estimate the gain length and to understand the requirements in terms of the beam quality to obtain a high gain, since the focusing structure is easily determined solely by the plasma density and the radius of the blowout/bubble region. In this paper, we present a detailed analysis of the beam requirements when the Pierce parameter ρ is much smaller than 1, as required for FEL-like amplification. In an ICL, the wiggler parameter K depends on the electron properties, so it can be different for each electron. Therefore, we show that both the beam energy spread and beam wiggler parameter spread should be limited and satisfy: (1.1a,b) Multi-dimensional particle-in-cell (PIC) simulations of ICL are performed to confirm that if those conditions are fulfilled then a good amplification is observed. As the spread limitations are a function of the Pierce parameter, this parameter should be carefully calculated. However, two important effects were neglected in most of the previous works (Chen et al. 1990;Whittum et al. 1990;Whittum 1992;Liu et al. 2007): (i) the radiation diffraction and (ii) the Pierce parameter dependence on the wiggler parameter K. These effects are included in our theoretical calculation of the Pierce parameter and the associated gain length, and are confirmed by PIC simulations in Lorentz boosted frames.
2. Theory 2.1. Radiation emission As a first step, we analyse the motion of an electron in an IC whose boundary is described by a radius, r b (ξ ) which depends on the variable, ξ = z − ct, with z the longitudinal coordinate (corresponding to the beam propagation direction), t the time and c the speed of light in vacuum. In general the motion of a particle moving near the speed of light in an IC can be described in terms to the so-called wake potential ψ ≡ (e/m e c 2 )(φ − cA z ) where φ and A z are the scalar potential and axial component of the vector potential generated by the IC, m e and e are respectively the electron mass and charge. The accelerating and focusing fields are obtained from (∂/∂ξ )ψ and (∂/∂r)ψ where we assume azimuthal symmetry and where r is the radial position. Inside the IC, (∂/∂r)ψ is given by (Lu et al. 2006) −k 2 p r/2 where k p ≡ ω p /c and ω p ≡ (n e e 2 / 0 m e ) 1/2 is the plasma frequency, with n e the plasma density and 0 the permittivity of free space. Note that these expressions are valid when the IC is created by long (negligible accelerating fields) or short pulse particle beams or lasers (large accelerating fields) and if there are large surface currents in the IC (as there is in the highly nonlinear channels). Therefore, in all these cases, the focusing force is m e c 2 k 2 p r/2 as was used by Esarey et al. (2002). With this focusing force, the Lorentz factor γ , the transverse radial position r and the transverse radial momentum p r of an electron with an initial longitudinal momentum p 0 (all momentum quantities are normalized to m e c), a maximum radius of oscillation r 0 and no azimuthal momentum, are given by γ = γ 0 + r 2 0 k 2 p sin 2 (θ r )/4, r = r 0 cos(θ r ) and p r = K sin(θ r ) with γ 0 = (1 + p 2 0 ) 1/2 , K = r 0 k p (γ 0 /2) 1/2 and θ r = −Kct/r 0 γ 0 + θ r0 = −ω β t + θ r0 , where θ r0 is the initial angle and ω β = ω p /(2γ 0 ) 1/2 is the betatron frequency. Here K γ 0 has been assumed. This assumption is made throughout the paper. Hereafter, the second-order terms proportional to γ −2 0 are neglected. The electrons wiggling in the focusing potential generate a betatron radiation with a fundamental wavelength λ 1 = 2πc/ω 1 with ω 1 = 4γ 2 0 ω β /(2 + K 2 ). The interaction between the electron beam and the radiation can lead to the amplification of the radiation. In order to get micro-bunching, the spread in the radiation wavelength must be limited. In an ICL, the K parameter depends on r 0 and γ 0 which can be different for each electron, so the radiation wavelength spread can be induced by both the beam energy spread and the K spread. A good approximation of the limiting spread can be found by assuming that λ 1 /λ 1 < ρ should be satisfied, much in the same way as for FELs (Huang & Kim 2007). Knowing that λ 1 = 2πc(2 + K 2 )(2γ 0 ) −3/2 /ω p , we find that the energy spread and K spread must then approximately satisfy the conditions given by (1.1).

ICL Pierce parameter and gain length
To further explore the optimal parameters for the ICL it is fundamental to determine the Pierce parameter. To start with, we analyse the bunching mechanism, which is a consequence of the energy exchange between the electrons and the radiation. We first consider an electron propagating in the z direction and a co-propagating electromagnetic (EM) wave. This wave is polarized in the x direction and characterized by its vector potential normalized to m e c/e: where A 1 and Ψ 1 are respectively the wave amplitude and phase. We assume that the electron oscillates in the (x, z) plane. We then define φ, the electron phase in the EM wave, and η, the relative electron energy as: with γ η the electron Lorentz factor after its interaction with the wave and z the longitudinal position of the electron averaged over one betatron oscillation. As shown in appendix A, the interaction with the wave leads to the following equations of motion for the electron in the (φ, η) phase space: where [JJ] = J 0 (K 2 /(4 + 2K 2 )) − J 1 (K 2 /(4 + 2K 2 )), with J 0 and J 1 the Bessel functions. Equation (2.5) indicates that a beam of electrons is bunched by the EM wave at the phase φ = −Ψ 1 + π/2 + 2mπ, with m an integer, which leads to a bunching at the position r = r 0 sin(k 1 z − ω 1 t + Ψ 1 ). Therefore, due to the correlation between the radial and longitudinal position, the electron beam gets a continuous and oscillating shape after the bunching, with a period equal to λ 1 . This is different from a conventional FEL, in which a succession of separated bunches is obtained.
Knowing the equations of motion, the amplification growth rate can be derived from the Vlasov and paraxial equations, as has been described by Huang & Kim (2007) for the conventional FEL case. As explained in appendix B, this method can be adapted to the ICL case by taking into account equations (2.4) and (2.5). The equivalent of the Pierce parameter ρ 1D and the gain length of the radiation power L 1D GP in the one-dimensional limit (radiation diffraction is neglected) for the ICL case is then given by: with I the beam current and I A ∼ 17 kA the Alfvèn current. We note that those results have been obtained assuming that K ρ 1/2 1D and ρ 1D 1. Using ρ 1D ∼ 1 may also lead to amplification, but the analytics have to be redone for this case. In addition, to avoid the damping of the bunching due to plasma oscillation in the beam, the gain length should be smaller than the longitudinal plasma oscillation wavelength characterized by its wavenumber (Rosenzweig et al. 1997): with n b = 4 0 m e Iω 2 β γ 2 0 /(I A K 2 e 2 ) the beam density. The condition L 1D GP 1/k pb then leads to: if ρ 1D ≈ (I/(2I A γ 0 )) 1/3 is assumed (limit obtained with K = 0). Note that with K > 0 or even K 1, then L 1D GP k pb is different from the approximation (2/K)(I/I A ) 1/6 γ −1/6 0 by only 30 % at maximum. The condition given by (2.9) can be difficult to fulfil for high current and low K cases.

Radiation diffraction effects
If the electrons have similar γ 0 and K values, then the beam transverse size is limited to 2r 0 = 2K(2/γ 0 ) 1/2 /k p . In an ICL, the radiation is emitted with a waist close to r 0 , so the associated Rayleigh length is Z r ∼ r 2 0 k 1 /2 L 1D GP , since: with ρ 1D 1. As a result, the radiation diffraction can reduce or even stop the amplification and it should not be neglected. This is a major difference to conventional FEL, where this limitation is not present.
As explained in appendix C, taking into account the diffraction can lead to the following solution for the Pierce parameter and the power gain length: (2.11) where Γ depends on a parameter µ, with Γ and µ given by: (2.14) with B(z, r) the amplitude of a Gaussian beam characterized by its waist W 0 , its wavelength λ 1 and B(0, 0) = 1, the focal plane being in z = 0. ρ and L GP correspond to the two-or three-dimensional (3-D) solution, depending on whether B is the solution of respectively the 2-D or 3-D paraxial wave equation. In two dimensions, B is thus given by: and in three dimensions, by: with: As explained in appendix C, a good approximation for the waist is W 0 = 3r 0 /4 in two dimensions and W 0 = 3r 0 /(4 √ 2) in three dimensions. The solution of the coupled equations (2.13)-(2.14) can be found iteratively: we start from Γ = 1 (1D limit), then (2.14) and (2.13) can be solved iteratively until a converged solution is obtained. Alternatively, as shown in appendix D, an analytical solution can be found if Z r L 1D GP is assumed. This approximated solution Γ is given in two or three dimensions by: (2.20) where ζ = (Z r e i(π/6) )/(L 1D

GP
√ 3), γ e ≈ 0.577 is the Euler-Mascheroni constant and LambertW is the Lambert-W function (also called the product logarithm).

Theory validation
In order to validate the theoretical conditions given in (1.1), we have performed 2-D simulations with the PIC code Osiris 2.0 (Fonseca et al. 2002). PIC codes are well suited to correctly and self-consistently model the radiation emission, diffraction, particle bunching and radiation amplification, as the full set of Maxwell's equations is solved. As the typical IC size is much larger than the radiation wavelength λ 1 , the IC formation is not self-consistently calculated in our simulation, allowing for a considerable reduction of the simulation size. We initialize our simulations with a preformed field profile that matches the IC focusing fields. A simulation technique that uses a Lorentz boosted frame (Vay 2007;Martins et al. 2010) is used in order to considerably speed up the calculations, by performing simulations in the beam frame instead of the laboratory frame. In this new frame, ω β = ω 1 , so the required number of time steps is reduced by a factor of 4γ 2 0 /(2 + K 2 ). For instance, a speed up of three orders of magnitude is obtained with γ 0 = 50 and K = 1. Moreover, running the ICL simulations in the beam frame prevents the numerical noise due to the numerical Cerenkov radiation (Godfrey 1974). The numerical noise can often perturb the bunching and artificially reduce or even stop the amplification. Perfectly matched layer (PML) absorbing boundary conditions (Vay 2000) are used on the transverse side of the box, and periodic boundaries are used in the longitudinal direction. In the boosted frame, the box length was chosen between 2λ r (for the shortest 3-D simulations) and 40λ r (for most of the 2-D simulations). The box transverse size was typically equal to 40r 0 . The longitudinal and transverse cell sizes used are typically dz = dr = λ r /50. In the following 2-D simulations, the total current I is meaningless due to the lack of the third dimension and only current density j can be properly used as an input for the simulation. However, for the sake of comparison with the real 3-D case, we will still introduce in two dimensions the beam current I defined as in three dimensions by I = πr 2 0 j. The self-consistent field amplitude in the simulation box is initially equal to 0, so the initial self-forces are neglected. This assumption is consistent with the fact that, in a FEL, the beam self-fields can be neglected as long as ρ 1 (Huang & Kim 2007). The condition L GP k pb 1 is also often fulfilled in the simulations discussed in the following.
All the physical values used as inputs for the simulation (initial particle momentums and positions, beam density, external focusing field) are first converted from their laboratory-frame values to the corresponding values in the boosted frame. Thanks to the periodic boundaries in the longitudinal direction, the radiated field amplitude is usually homogeneous in the longitudinal direction. Assuming that the field propagates at c, the average radiated power in the boosted frame can then be deduced from the total radiated field energy integrated in the simulation box and divided by L/c, with L the simulation box length. This power value is then converted to its corresponding value in the laboratory frame, still assuming that the radiated field propagates at c.
In figure 1, the simulation results for a beam characterized by γ 0 = 50, K = 1 and a current I = 0.8 kA injected in the IC field are presented. The beam parameters are chosen such that the computational costs of the simulations are reduced but the main physical features are captured. In the simulations, γ and K are initialized within a Gaussian distribution and the electrons are initialized with a random angle θ r0 . If γ = K = 0, the Pierce parameter and power gain length determined in the 1-D limit or in two dimensions are given by respectively ρ 1D = 0.082, L 1D GP = 8.4 c/ω β , ρ 2D = 0.048 and L 2D GP = 13.4 c/ω β . The corresponding approximated Pierce parameter and FIGURE 1. Evolution of the radiation growth as a function of the energy spread (a) and K spread (b). Two-dimensional simulations with γ 0 = 50, K = 1 and I = 0.8 kA. (a) The green, light blue, dark blue, red and purple curves correspond to respectively E/E = 0, The green, light blue, dark blue, red and purple curves correspond to respectively K/K = 0, K/K = 0.02, K/K = 0.04, K/K = 0.08 and K/K = 0.12. The dotted black and dotted red lines correspond to respectively the theoretical growth rate in the 1-D limit and in two dimensions. The γ and K spreads correspond to root-mean-square values. power gain length given by (2.19) are ρ 2D ap = 0.05 and L 2D GP,ap = 12.7 c/ω β . The 1-D and 2-D theoretical growth rates are also represented in figure 1. We can observe a very good agreement between the 2-D theoretical growth rate and the simulation results. Note that L 2D GP k pb = 0.82 in this case so the longitudinal plasma oscillation in the beam can be neglected, which is confirmed by the good agreement between the theory and the simulation. In the simulation, the initial noise produced by the macro-particles is amplified up to the saturation level. This is reached when the particles are fully bunched. However, with a high γ or K spread, the growth rate is reduced or even stopped. We observe that the change between a maximal and reduced growth rate matches the theoretical limits given by K/K = 0.072 and γ /γ = 0.032 with ρ 2D = 0.048.

Towards more realistic beams
The condition K/K 1 can be parameterized by different complex configurations of the electron distribution in the transverse phase space. For example, in the 2-D case, the electrons can be distributed over a ring in the transverse phase space. This ring is parameterized by r = r 0 cos(θ r ) and p r = K sin(θ r ). We propose more realistic distributions, with a spot shape instead of a ring shape. In a first configuration, K ∼ 1 and K/K 3ρ/2 1 are used, so (1.1) is satisfied, but the electrons are only distributed over a ring fraction, with an initial angle θ r0 that satisfies |θ r0 | < θ r,max . If θ r,max π, the initial beam transverse size is much smaller than r 0 and the beam corresponds to an off-axis injected beam oscillating in the IC. In that case, the beam shape in the transverse phase space is close to a spot with an initial transverse size and transverse momentum spread approximately equal to respectively r 0 K/K and Kθ r,max . In a second configuration, we choose K ∼ ρ 1/2 1 and K/K ∼ 1, which still satisfies the K spread condition in (1.1). In that case, as K/K ∼ 1, the maximum radial momentum p rm of a given electron roughly satisfies K − K p rm K + K so |p rm | 2K. We also have |r 0m | 2r 0 with r 0m the maximum radial position of a FIGURE 2. (a) In blue: radiated power for an off-axis beam, initialized with |θ r0 | < π/16, K = 1, K/K = 0.01 and γ /γ = 0.005. In green: radiated power in the reference case, with |θ r0 | < π and K/K = γ /γ = 0. (b) In blue: radiated power for an on-axis beam, which is initialized with |θ r0 | < π, K = 0.1, K/K = 0.3 and γ /γ = 0.002. In green: radiated power in the reference case, with K/K = γ /γ = 0. The dotted black and dotted red lines correspond to respectively the theoretical growth rate in the 1-D limit and in two dimensions.
given electron. Therefore, the spread around the ring is such that the beam distribution in the transverse phase space becomes a spot. As r 0 ∝ K, using K 1 corresponds to a narrow on-axis injected beam.
The two configurations are highlighted by 2-D simulations. In the first case, an offaxis beam with γ 0 = 50, K = 1 and I = 0.27 kA (L 2D GP k pb = 0.73) is injected with |θ r0 | < π/16. The corresponding Pierce parameter is ρ 2D = 0.031 (ρ 2D ap = 0.032) and the beam is initialized with K/K = 0.01 and γ /γ = 0.005. In the second case, an on-axis beam with γ 0 = 50, K = 0.1 and I = 42 A (L 2D GP k pb = 11.6) is injected with |θ r0 | < π. The corresponding Pierce parameter is ρ 2D = 6.4 × 10 −3 (ρ 2D ap = 6.4 × 10 −3 ) and the beam is initialized with K/K = 0.3 and γ /γ = 0.002. Reference simulations have been performed for both cases, using |θ r0 | < π and K/K = γ /γ = 0. The evolution of the amplified radiation power for these different simulations is presented in figure 2. In both cases, we observe that the use of more realistic beams, with a finite spot in the transverse phase space and an energy spread, can still lead to exponential radiation amplification, even if the growth rate and final power are lower than in the reference simulations, for the idealized scenarios. The discrepancy between the 2-D theoretical growth rate and the idealized simulation result in the on-axis case can be explained by two reasons: (i) the use of K = 1.25ρ 1/2 whereas our theoretical model is valid in the limit K ρ 1/2 , and (ii) the fact that L 2D GP k pb = 11.6 > 1 so the longitudinal plasma oscillation can significantly damp the bunching and reduce the growth rate. Despite these facts, it is interesting to see that an exponential growth of the radiation is still obtained in the simulation.

Three-dimensional case
We have also performed 3-D simulations to confirm the 3-D theoretical results. The electrons are initialized with a radial momentum but no azimuthal momentum, so the electrons still oscillate in a plane and do not gain helical trajectories. The uniform distribution of the electrons along the azimuthal angle leads to the initialization of a cylindrical beam. The results are shown in figure 3. In the simulation with γ 0 = 50, FIGURE 3. (a) Radiation growth with γ 0 = 50, K = 1, I = 8 A, K/K = γ /γ = 0 in a 3-D simulation (blue) and given by the 1-D theory (dotted black) and 3-D theory (dotted red). (b) Radiation growth with γ 0 = 50, K = 1, I = 0.8 kA, K/K = 0, γ /γ = 0.2 % in a 3-D simulation (blue) and given by the 1-D theory (dotted black) and 3-D theory (dotted red). (c) Shape of the electron beam at saturation in the 3-D simulation with I = 0.8 kA: a helical bunching is observed (iso-surface of the electron density). K = 1 and I = 8 A, the corresponding Pierce parameter and power gain length obtained in the 1-D limit or in three dimensions are given by respectively ρ 1D = 0.018, L 1D GP = 39 c/ω β , ρ 3D = 2.7 × 10 −3 and L 3D GP = 226 c/ω β (L 3D GP k pb = 1.38). The approximated values using (2.20) are ρ 3D ap = 2.8 × 10 −3 and L 3D GP,ap = 217 c/ω β . A good agreement between simulation and theory is found (figure 3a). Since the initial noise in the simulation is too low to start the amplification mechanism in the 3-D simulations, we have injected a circularly polarized seed in the IC. The seed wavelength is λ 1 , like the expected amplified radiation. As the seed diffracts, most of its energy gets out from the simulation box from the side. This explains the power dip at the beginning of the simulation at t ∼ 500 ω −1 β . The amplification is initiated and the saturation level is reached at the end of the simulation. A similar case has been run with I = 0.8 kA and γ /γ = 0.2 % (figure 3b), leading to ρ 1D = 0.082, L 1D GP = 8.4 c/ω β , ρ 3D = 0.0225, L 3D GP = 26.9 c/ω β and L 3D GP k pb = 1.65 (ρ 3D ap = 0.0235 and L 3D GP,ap = 25.5 c/ω β ). The presence of the small energy spread and the fact that L 3D GP k pb is farther away from 1 may explain the difference between the theoretical prediction and the simulated growth in this case. Nevertheless, as expected, the bunch shape is helical at saturation (figure 3c). This result is consistent with a circularly polarized seed.

Discussion
Our results show that an amplification of several orders of magnitude of the radiated power can be achieved, even if the Rayleigh length of the generated radiation is shorter than the gain length in the typical ICL configurations. The diffraction is responsible for the growth rate reduction. The gain length is 1.6 times larger than the 1-D limit in the 2-D case presented in figure 1, and respectively 5.8 and 3.2 times in the 3-D cases with I = 8 A and I = 0.8 kA. We have also confirmed in our ab initio simulations that the amplified radiation wavelength and the oscillating shape period of the bunching are λ 1 , matching the theory. Odd harmonics have also been observed. Yet, as expected with K = 1, their amplitudes are much smaller than the fundamental harmonic amplitude. This demonstrates that PIC simulations in the beam frame might be an efficient tool to study the self-consistent dynamics of harmonics and its feedback to the growth rate in scenarios where K > 1 in an ICL or in FELs.
Numerical applications of our analytical results show that the most stringent conditions will be (i) to inject an electron beam with a very low emittance at a precise radius in the IC and (ii) to generate a stable IC over a long enough distance while the beam longitudinal acceleration remains negligible. For instance, if we consider a laser or particle beam driving a wakefield in a plasma with a density n e = 5 × 10 17 cm −3 , the injection of a 25 MeV, 0.8 kA beam with K = 1 can generate a source with a wavelength λ 1 = 146 nm. The associated gain length is L 3D GP = 2 mm, so the radiated power can be multiplied by 1000 after 13.7 mm if an ideal electron beam is considered. To get this amplification, the electron beam should have a relative energy spread lower than 1.5 % and be injected at r 0 = 1.5 µm off-axis with a normalized transverse emittance N < 0.02 mm mrad. The production of shorter wavelengths can be achieved by increasing the electron energy. For instance, a 250 MeV beam (the other parameters are kept constant) would produce photons with λ 1 = 4.6 nm. However the higher energy induces a longer gain length (1.8 cm in this case) and a more stringent limit for the energy spread (0.5 %). Increasing the beam current can help to reduce these constraints. Using a 250 MeV beam, 10 kA beam with K = 1 leads to λ 1 = 4.6 nm and L 3D GP = 5.6 mm, so the radiated power can be multiplied by 1000 after 3.8 cm. The effect of the longitudinal plasma oscillation should stay limited in this case as L 3D GP k pb = 1.66. The relative energy spread limit is then 1.7 % and the beam should be injected at r 0 = 0.48 µm with N < 0.01 mm mrad. As all the lengths scale with 1/ω p (e.g. see (2.12)), the radial position of the beam injection (r 0 ) and the transverse emittance limit can become larger with a lower density (∝ n −1/2 e ), but the radiated wavelength and the gain length also increase with the same factor. Using n e = 1.25 × 10 17 cm −3 still with the same 250 MeV and 10 kA beam is less restrictive as the beam need to be injected at r 0 = 0.96 µm with N < 0.02 mm mrad. However, λ 1 is elongated to 9.2 nm and L 3D GP = 11.3 mm. The other parameters (energy spread limit, L 3D GP k pb ) stay unchanged. Even if the emittance value is still far from the best values obtained in a wakefield accelerator, optimization or mix of new injection schemes in a laser wakefield accelerator, such as optical (Faure et al. 2006;Davoine et al. 2009), ionization (McGuffey et al. 2010Pak et al. 2010) or magnetic (Vieira et al. 2011) injection, might help to improve emittance and control the off-axis injection. It is also important to note that the accelerating field, present in a typical IC can affect the amplification process, as the electron energy will change in time. This effect can be reduced, for instance, by injecting the electron beam close to the centre (longitudinally) of the bubble in a wakefield, where the accelerating field is zero. To keep a stable IC structure and a negligible accelerating field over several millimetres or centimetres, which can be challenging, a particle beam driver instead of a laser driver can be used as its propagation in a plasma is usually more stable. Even if these issues should be handled before demonstrating the possibility of generating an ICL, the constraints on the IC generation and on the beam injection techniques are outside of the scope of this work, as this paper focuses on the constraints on the beam parameters and on the derivation of the correct ICL growth rate.

Summary
In this paper, we have determined analytically the amplification growth rate of an ICL while taking into account the diffraction effect. The required conditions on the electron beam quality in order to observe ion channel lasing have also been presented. It is shown that it is not necessary to use a guiding structure for the radiation as was considered in previous work on ICLs: the radiation defocusing reduces the growth rate but does not stop the amplification. Two-and three-dimensional PIC simulations, which are the first fully relativistic electromagnetic 3-D simulations of ICL, have confirmed our analytical findings, illustrating the possibility of achieving high-gain radiation amplification in ICL. Despite the still needed efforts to experimentally reach a sufficient beam quality and generate the required stable IC, these results pave the way for the generation of high brilliance coherent radiation in compact plasma structures.
In this appendix, we show that the interaction between an electron following a betatron motion and an EM defined by (2.1) leads to the motion equations in the (φ, η) phase space given by (2.4) and (2.5).
In an ICL, K is a function of γ . Therefore, we first need to determine how K evolves, as well as γ 0 , p 0 or r 0 , when the EM wave exchanges energy with the electron. This result is first presented and the description of the bunching process, leading to (2.4) and (2.5), is addressed in the second part of this appendix.
A.1. Influence of an EM wave on the betatron oscillation parameters As mentioned, K can evolve and is now a function of time. We defined K 0 and the time-dependent longitudinal momentum p z and maximum radius r m such that K 0 = K(t = 0), p z (t = 0) = p 0 and r m (t = 0) = r 0 . We still consider that γ 0 = (1 + p 2 0 ) 1/2 . In the following, we use the notation:Ẋ = dX d(ω 1 t) .
(A 1) In the presence of an EM wave defined by (2.1), the energy and momentum change of an electron following betatron motion in the (x, z) plane is given by: where α = A 1 sin(k 1 z − ω 1 t + Ψ 1 ) and β r and β z are the normalized transverse and longitudinal electron velocities. In the above equation, we have used the following identity to define the ion-channel focusing field, which is normalized to m e cω 1 /e: By using r = r m cos(θ r ), p r = K sin(θ r ) and K = r m k p (γ /2) 1/2 , we can show that: r m k p = −γ p 2 r + 2γ p rṗr + γ 2 rṙk 2 p r m k p γ 2 .
(A 9) AsK = 0 andṙ m = 0 when α = 0, we can simplify the equation and get: We now consider the average of the derivatives over one betatron period, and we assume that the change of γ , K and r m is small during one betatron period (γ ω β γ ). The averaged derivatives are then given by: We introduce the parameter ν defined as: Then, we obtain:K We also assume that all parameter evolutions are small during the whole interaction (e.g. γ (t) − γ 0 γ 0 for all t). This leads to: We introduce η = t 0 ν dτ . According to our last assumption we have η 1 and K 0 η 1/2 . By neglecting the terms proportional to η 2 , we then obtain: where we rename γ , K and r m by respectively γ η , K η and r η for convenience. Hereafter, all the items with η as a subscript are functions of γ η , K η or r η , and if η is not mentioned, it means that the value is taken at η = 0. For example, ω βη = ω p /(2γ η ) 1/2 , and from this equation we find: We also deduce from ω β /ω 1 = (2 + K 2 0 )/(4γ 2 0 ) that: A.2. Electron motion in the presence of an EM wave We can now rewrite (A 12) as: Here, only the dominant term is relevant and the terms proportional to η can be neglected. This is also true for the terms in the phase k 1 z − ω 1 t + Ψ 1 . Thus, we determine k 1 z while neglecting the terms proportional to η: where z 0 and θ r0 are the initial position and phase. We define a new phase φ as By using (A 27), we can note that φ is a constant of time when we neglect the terms proportional to η, so φ = φ. We then find that: where J 0 and J 1 are the Bessel functions, we can then write: To get the derivative of φ with time, we need to rewrite (A 28) without neglecting the terms proportional to η. Provided that: we eventually obtain the equation of motion in the (φ, η) phase space:

Coupling between the Maxwell and motion equations
In this appendix, we follow the method used in Huang & Kim (2007) to calculate the one-dimensional and ideal growth rate for FEL. Here, this method has been adapted to the ICL context.
To start with, we consider the presence of an EM wave polarized along the x direction, propagating along the z direction and characterized by its normalized vector potential A x = A 0 (x, y, ξ , τ ), where ξ = k 1 z − ω 1 t and τ = ω 1 t. We introduce the wave amplitude A ν (x, y, τ ) in the frequency domain through: where c.c. is the complex conjugate. The Maxwell equations for A ν give: where ∇ 2 ⊥ is the transverse Laplacian normalized to k 2 1 . j x is the transverse current density along the x direction, and it is normalized to eω 1 k 2 1 . I A = ec/r e is the Alfvèn current, with r e = e 2 /(4π 0 m e c 2 ) the classical electron radius. By using the slowly varying envelope approximation (|∂ 2 A ν /∂τ 2 | 2ν|∂A ν /∂τ |), we get: Calculation of the transverse current The normalized transverse current density j xn of the particle n, which follows betatron motion in the (x, z) plane is given by: where K n , γ n and θ rn are the parameters of the electron at the time τ , δ is the Dirac function and x n and ξ n are respectively the transverse position of the electron and its position over the ξ direction at the time τ . x n and ξ n are given by: x n = r m,n cos(θ rn ) (B 6) where r m,n is the maximum radius of oscillation of the electron n at time τ . Equation (B 7) is obtained by using (A 28). j νn is then given by: j νn (x, τ ) = − K n 2iπγ n e −i( νθ rn +νφ n ) δ[k 1 x − k 1 r m,n cos(θ rn )]e iν(K 2 n /(4+2K 2 n )) sin(2θ rn ) (1 − e −2iθ rn ).

(B 8)
Where ν = ν − 1. We define the function G(x, ν, K, γ , θ r ) as: where r m is a function of γ and K, since we have the following identity: The current created by the electron n is then given by: j νn (x, τ ) = − K n 2iπγ n e −i( νθ rn +νφ n ) G(x, ν, K n , γ n , θ rn ).

(B 11)
The electron distribution at the time τ in the phase space can be parameterized by the 4 parameters φ, η, K and θ r . Therefore, the distribution function F is given by: where I is the longitudinal beam current (absolute value so I > 0), N is the number of electrons and 2πeω 1 /I is a normalizing factor. The total current j ν (x, τ ) is then given by: where, according to the definition of η, we have γ = γ 0 (1 + η).
The normalizing factor of F has been chosen so that if we consider a beam distribution with the parameters φ, η, K and θ r which are not correlated, then F can be written as: where D 1 dφ = Lk 1 (B 15) with L = Nec/I the beam length.

B.2. One-dimensional approximation
We now consider that the EM wave is a plane wave, so the term ∇ 2 ⊥ and the transverse position can be neglected. Equation (B 4) then becomes: where j ν (τ ) is the current averaged over the beam transverse size and is given by: G 2 (ν, K, γ , θ r ) = e iν(K 2 /(4+2K 2 )) sin(2θ r ) (1 − e −2iθ r ), (B 21) with S = πr 2 0 the beam transverse size. Since only the radiations with a wavelength close to λ 1 are generated and amplified, we assume that ν 1, so the term νθ rn + νφ n evolves slowly and can be considered as constant over one betatron period. We also get G 2 = [JJ]. The average over one betatron period of the current is then: where we have also considered a small energy and K spread, so K ∼ K 0 and γ ∼ γ 0 . We finally get: where we have also assumed that exp[−i νθ r ] ∼ exp[i ν(ω β /ω 1 )τ ]. Indeed, as we have the three identitiesθ r = −ω βη /ω 1 , ω βη ∼ ω β and ν 1, then the difference between νθ r and − ν(τ ω β /ω 1 + θ r0 ) is small, even if τ 1. Moreover, according to equations (B 6), (B 7), it is possible to choose θ r0 so that θ r0 ∈ [0, 2π], so νθ r0 1.
The Vlasov equation is defined byḞ = 0. We thus get: Moreover, with an EM wave described by (B 1), (A 37) becomes: νA ν e i( νθ r +νφ) dν + c.c., (B 25) where the dependence of [JJ] as a function of ν has been neglected. We introduce the following scaled variables to simplify our equation: We thus obtain: By defining ρ (the equivalent of the Pierce parameter in FEL theory) as: we find: Moreover, the Vlasov equation becomes: where X = dX/dτ . B.3. Calculation of the growth rate To calculate the growth rate, we need to solve the coupled equations (B 34) and (B 35).
Equation (B 35) can be linearized in the small signal regime before saturation when the scaled radiation field is small, i.e.: a ν e iνφ d ν + c.c. =η 1. (B 36) Let us split f in two parts: where f 0 is the distribution function averaged over φ and f 1 contains the noise fluctuation and the modulation induced by the bunching. The average over φ of (B 35) leads to: The small signal regime also implies that f 1 f 0 . We can then assume that the second term on the left-hand side of (B 38) can be neglected, which leads to: The corresponding equation for f 1 is therefore: To solve this equation, we consider the trajectory of an electron, which is parameterized by φ (0) ,η (0) , K (0) and θ (0) r . According to the Vlasov equation (B 35), we have: where φ (0) ,η (0) , K (0) and θ (0) r are here given at the time s. Thanks to equation (B 40), we can write: So: ) where φ,η, K and θ r are the values of φ (0) ,η (0) , K (0) and θ (0) r at timeτ . Moreover, we have: So we find that: In the small signal regime, provided thatη 1, we can assume thatη (0) (s) ∼η. From (A 22), we can deduce that K = (1 + K 2 0 )ρη /(2K 2 0 ) so K 1 and K (0) ∼ K. Based on (B 39) and on the definition of θ (0) r , we can deduce that f 0 (θ (0) r , s) is a constant if we assume that K = 0. Therefore, f 0 (θ (0) r , s) = f 0 (θ r ,τ ), which leads to: As f 0 does not depend on φ, we have the following result if we assume that the electron beam is very long in comparison to the fundamental radiation wavelength λ 1 : Then, (B 34) becomes: By using (B 45), (B 46) and (B 48), we obtain: This equation shows that each frequency component of the radiation field is independently amplified. The right-hand side of (B 49) corresponds to the initial fluctuation and is the source term that creates the initial radiation in the absence of seed.
To determine the growth rate, we only consider the homogeneous part of (B 49). We seek a solution in which a ν is proportional to exp(−iμτ ), whereμ is the complex growth rate. Then, we have a ν (s) = a ν (τ ) exp [−iμ(s −τ )]. This leads to: We first calculate the integral over the time s. Then, we assume that η, K and θ r are not correlated at the time τ . Thanks to (B 14)-(B 18), the integration over K and θ r leads to: Here we have also assumed that |exp[iμτ ]| 1, as exp[−iμτ ] is supposed to grow exponentially with time. After integrating by part overη, we obtain: In the limit where there is no energy spread ( f 0 (η) = δ(η)), this equation becomes: At the optimal frequency ( ν = 0), we obtain: The solution with the largest imaginary part is associated with the largest growth rate. Thus, we only consider the following solution: Thanks to equation (B 27), we finally find that the field amplitude is proportional to: In the following, the parameter ρ given in the 1-D approximation by (B 33) will be referred as ρ 1D . The 1-D gain time for the field amplitude is then: The associated power or intensity gain length is then: Appendix C. Transverse effect: influence of the Rayleigh length In appendix B, we have assumed that the transverse variation of A ν can be neglected, as we have used ∇ 2 ⊥ = 0 to simplify (B 4). However, if we consider that the electron beam creates a radiation beam with a waist close to the electron beam radius r 0 , then the associated Rayleigh length is Z r ∼ r 2 0 k 1 /2. This length is much shorter than the gain length, since: so Z r /(L 1D GP ) 1 as ρ 1. Therefore, the intensity of the emitted radiation is strongly reduced after one gain length, which reduces the growth rate. To take into account this phenomenon, we assume that the current j ν generates a Gaussian beam with a waist W 0 = w 0 r 0 , where w 0 is a free parameter. It implies that the source j ν also gets a Gaussian transverse distribution. We therefore assume that j ν (τ , x) = j ν (τ , 0) exp(−x 2 /W 2 0 ) in two dimensions. W 0 (and thus w 0 ) should be chosen so that the function j ν (τ , 0) exp(−x 2 /W 2 0 ) + c.c. provides the best fit of the real j ν + c.c. given by (B 13). Assuming ν ≈ 1, the transverse distribution of the average current j ν + c.c. generated by one particle over a betatron period only comes from the radial dependence of the function G(x, ν = 1, K, γ , θ r ). The use of (B 9) then leads to: ∝ e i(K 2 /(4+2K 2 )) sin(2 arccos(x/r 0 )) (1 − e −2i arccos(x/r 0 ) ) + c.c.
where G(x/r 0 , K) corresponds to the transverse distribution of the current. As can be seen in figure 4, G can be well approximated by a Gaussian with root-mean-square value σ = 3/4. In this whole study, we have then use w 0 = 3/4 and thus W 0 = 3r 0 /4 in two dimensions. For the 3-D case, that fact that the electrons oscillate along different transverse direction reduces this transverse size and we have always use w 0 = 3/(4 √ 2). These are not analytically determined values. Nevertheless, these approximations allow us to obtain a good agreement with the simulation results.
In the following, for the sake of simplicity, we consider only the 2-D case, so ∇ 2 ⊥ = ∂ 2 /∂(k 1 x) 2 . (B 4) then becomes: We then use a method based on Green's functions to solve this equation, with G ν (τ , x) the solution of the equation with j ν (τ , 0) = δ(τ ): We introduce the function B which describes a Gaussian beam with a waist W 0 : B is a solution of the 2-D paraxial wave equation, so: The function G ν (τ , x) is then given by: where H is the Heaviside function. A ν (τ , x) is then given by: The limits of the integral can be changed from (−∞, +∞) to (0, τ ) because we consider that nothing happens when τ < 0 (i.e. j ν = A ν,τ = 0 if τ < 0), and thanks to the presence of the function H, we have G ν (τ −τ , x) = 0 ifτ > τ . As in appendix B, we seek a solution where the current and field are proportional to exp(−iµτ ), with µ the complex growth rate. The current thus satisfies j ν (τ , 0) = j ν (0, 0) exp(−iµτ ). By defining we obtain We can also deduce from (C 6) that: 2iν ∂ ∂τ + ∂ 2 ∂(xk 1 ) 2 A ν (τ , x) = 2iνJ ν0 e −iµτ B(0, x).

(C 19)
To go further, we now assume that the term τ 0 e iµτ B(τ , x) that appears in (C 18) becomes constant after some time (after few gain times). Indeed, as we have supposed that A ν is exponentially growing, then e iµτ is exponentially decreasing and the integral stays constant if τ 1/Im(µ). This assumption has been verified numerically: the integral reaches a nearly constant value after few gain times. We can then write: (C 20) So: If we define the function Γ (x), which is constant with time, as follows: then we obtain from (C 19): ∂τ .
(C 30) To be more consistent, we can rewrite (C 22) as: The final solution can be found by an iterative method. We first start from the 1-D result Γ = 1. With this value, we can then solve (C 29) and (C 31). Finally, by solving iteratively those two equations, the result found after few loops converges to the solution of those two coupled equations.

Appendix D. Analytical solution for the 2-D and 3-D gain length and Pierce parameter
In two and three dimensions, the Pierce parameter and gain length are respectively given by (2.11) and (2.12). This implies that we first calculate the value of Γ by solving the coupled equations (2.13) and (2.14), which is equivalent to solving the following equation: with ζ = (Z r e i(π/6) )/(L 1D