Introduction

X-ray radiation sources have important applications in biology, medicine, industry and fundamental science because of their capability of resolving the structure and dynamics of matter on the molecular and atomic scales1,2,3. For the past twenty years, the third-generation synchrotron facilities have been the workhorse X-ray sources4. By the virtue of an electron beam instability called microbunching5, much brighter, shorter, and fully coherent XUV and X-ray pulses can be generated in the devices called free-electron lasers (FELs)5,6,7. The excellent FEL radiation properties make it possible to time-resolve molecular structural dynamics and obtain high-resolution images. However, synchrotron and FEL facilities are typically large, expensive, and oversubscribed for the users.

It has been demonstrated that schemes based on laser-plasma interaction are able to produce ultrashort and bright radiation8,9,10,11,12,13,14,15,16. In a typical scenario, electrons accelerated by a laser-excited wakefield17,18, oscillate either in the focusing fields of the wake or inside an (external) undulator, emitting bright radiation with a fundamental wavelength of \(\lambda ={\lambda }_{u}\mathrm{(1}+{K}^{2}\mathrm{/2)/(2}{\gamma }_{0}^{2})\) on axis, where λ u is the undulator wavelength, γ 0 is the Lorentz factor of the electron, and K is the so-called undulator strength parameter. Compared to a conventional magnetic undulator12,13, the period of a plasma-based undulator is short and can be less than a millimeter. Therefore, a laser-created plasma undulator together with a laser-plasma electron accelerator (LPA) make it possible to construct an economical and compact incoherent XUV or X-ray source for university laboratories, hospitals, and even commercial applications15,19. Plasma undulators can be realized by ultra-intense laser plasma interactions in the bubble regime8,9,10, a laser pulse propagating in plasma perpendicularly to the electron beam propagation direction20, a laser pulse interaction with a nanowire array21, or using laser pulse offset injection in a matched plasma channel22,23,24. However, it is still an open question whether these plasma undulators can be used as an FEL. In general, not all undulators are suitable for the onset of the FEL instability. For example, in magnetic undulators matching and field uniformity are critical. In the case of plasma undulators one of the major challenges is the large radiation spread caused by varying values of undulator strength K throughout the beam or by strong focusing and hence large electron beam divergence inside the wakefield, while for the FEL process to develop a very narrow bandwidth is required8,15. Another challenge is that the phase slippage between the electrons and the wakefield limits the length of a plasma undulator. Furthermore, it is not generally the case in plasma undulators that the electron trajectories are independent of the injection positions. We show for the first time that a solution exists for a plasma undulator that meets the stringent conditions required for the FEL lasing.

Here we propose a phase-locked plasma undulator created by the wakefields of a combination of high-order laser modes propagating in a parabolic plasma channel. We demonstrate that the undulator fields can be made uniform along the transverse direction by choosing appropriate intensities of the laser modes. This enables all the electrons to oscillate with the same strength parameter K on the order of unity for a few tens of undulator periods. The plasma density in the channel can be tapered to lock the phase between the electrons and the wakefield, which significantly increases the number of electron oscillations, i.e. the total undulator length. As a result, X-ray radiation with high brightness and narrow bandwidth is generated when a high-energy electron beam is injected into the present plasma undulator. The beam loading limit indicates that the current of the beam can reach approximately 0.3 kilo-Ampere for typical laser-plasma parameters. These properties imply that such a plasma undulator is a miniaturised electron device naturally matching the extremely compact scale of a plasma accelerator, similar in significance to the recent breakthrough development of plasma lenses25,26,27, and plasma accelerator staging28,29, and may have great potential in incoherent XUV and X-ray sources or future compact FELs.

Results

Principles of the plasma channel undulator

Laser pulse guiding for distances much larger than the Rayleigh length of the laser pulse are necessary in LPAs for achieving GeV level of electron energy. Typically a plasma channel that has a parabolic (or close to parabolic) transverse density distribution is used in LPA experiments for laser guiding30. For laser power below the critical power18, laser pulses with transverse profile given by Hermite-Gaussian (or Laguerre-Gaussian) modes will propagate inside the channel without changing their transverse shape given that their spotsize is matched to the channel radius. Plasma undulator can be generated by propagating a mixture of different Hermite-Gaussian laser modes in a matched plasma channel, as shown by the schematic in Fig. 1. An oscillatory behaviour of intensity envelope appears when the modes with the same polarization co-propagate in the plasma channel, since the phase velocities for different modes are dependent on the mode numbers. For example, the wavelength of the oscillations of the total intensity of a mixture of two modes |m 0, p 0〉 and |m 1, p 1〉 is 2πZ R /|(m 0 + p 0) − (m 1 + p 1)|31, where the quantum ket notation |m, p〉 indicates a Hermite-Gaussian mode with an order number m in x direction and an order number p in y direction, and \({Z}_{R}=\pi {w}_{0}^{2}/{\lambda }_{L}\) is the Rayleigh length with w 0 the laser spot size and λ L the wavelength of the laser pulses, which is assumed to be the same for all the modes. It is also important to mention that to obtain the desired undulator fields one has to create an asymmetric initial transverse intensity distribution of the modes mixture, which is achieved by mixing even and odd mode numbers. The wakefield generated by the two modes will also oscillate while propagating in the channel, which provides an additional control of the focusing field. Consider a two-dimensional case, then the transverse wakefield generated by two y-polarized modes m 0 and m 1 can be expressed as (see Methods for the calculation of the wakefield)

$${E}_{x}=\frac{{E}_{0}C}{{k}_{p}}[\frac{4{\alpha }_{\beta }x}{{w}_{0}^{2}}+\frac{4{\alpha }_{u}}{{w}_{0}}\,\cos \,\frac{z({m}_{0}-{m}_{1})}{{Z}_{R}}]\,\sin ({k}_{p}\zeta ),$$
(1)

with

$${\alpha }_{\beta }=\frac{-{a}_{{m}_{0}}^{2}(2{m}_{0}+1)}{{m}_{0}{!2}^{{m}_{0}}}{[\frac{{m}_{0}!}{(\frac{{m}_{0}}{2})!}]}^{2}+\frac{4{a}_{{m}_{1}}^{2}}{{m}_{1}{!2}^{{m}_{1}}}{[\frac{{m}_{1}!}{(\frac{{m}_{1}-1}{2})!}]}^{2},$$
(2)

and

$${\alpha }_{u}=\frac{{a}_{{m}_{0}}{a}_{{m}_{1}}\sqrt{{m}_{0}!{m}_{1}!}}{\sqrt{{2}^{{m}_{0}}{2}^{{m}_{1}-1}}}\frac{{(-\mathrm{1)}}^{({m}_{0}+{m}_{1}-\mathrm{1)/2}}}{\frac{{m}_{0}}{2}!\frac{{m}_{1}-1}{2}!},$$
(3)

where ζ = z − ct, \(C=(\sqrt{\pi }{k}_{p}L\mathrm{/4)}\exp (-{k}_{p}^{2}{L}^{2}\mathrm{/4)}\), E 0 = m e c 2 k p /e, L is the laser pulse longitudinal length, k p is the plasma wave number, e is the electron charge, m e is the electron mass and c is the light speed in vacuum. As can be seen from equation (1), the transverse wakefield is separated into two parts: the betatron part which is linearly proportional to the transverse coordinate x, and the undulator part which is a simple cosine function of the longitudinal coordinate z. The strengths of these two parts α β and α u are determined by both the laser modes intensities and order numbers. Electrons with Lorentz factor γ 0 injected in such a wakefield will experience two kinds of oscillations: betatron oscillations with a wave number \({k}_{\beta }=\sqrt{\mathrm{4|}{a}_{\beta }|C/{\gamma }_{0}{w}_{0}^{2}}\) and undulator oscillations with a wave number k u  = |m 0 − m 1|/Z R . Such a plasma wakefield structure serves as an undulator (or wiggler) and leads to electron oscillations and generation of synchrotron radiation. In contrast to the plasma undulator generated by an off-axis injected laser pulse22, in which case it is difficult to suppress the strong focusing field (large α β ), in the plasma undulator created by high-order modes one can eliminate the betatron part by choosing appropriate laser intensities a m and a n to satisfy α β  = 0 in equation (2). Physically, the focusing field disappears because the sum of the intensities of the two modes keeps constant near the axis, and thus the transverse gradient of the potential reduces to zero. As a result, electrons injected into such a particular wakefield will oscillate only with the undulator frequency k u . This significantly enhances the undulator radiation, which is more interesting because its energy can be located in the soft or even hard X-ray range and its radiation bandwidth is very narrow. It is also worth noting that the undulator strength parameter here K = 4πCα u w 0/λ L , is independent of the electron transverse positions, unlike the typical betatron case8,9,10,14. This property is helpful for decreasing the undulator radiation bandwidth, which should be smaller than the Pierce parameter ρ for the FEL application5,8,15.

Figure 1
figure 1

Schematic of the plasma channel undulator. As an example, three y-polarized Hermite-Gaussian modes, |0, 0〉, |0, 1〉, and |1, 0〉, propagate in a matched parabolic plasma channel. In both the (x, z) plane and the (y, z) plane, the profile of the total intensity oscillates with a wavelength proportional to the Rayleigh length of the pulse, and can be controlled by choosing the mode numbers. The ellipticity of the laser pulse oscillations can be controlled by the phase difference between the modes |0, 1〉 and |1, 0〉. The focusing fields of the induced plasma wakefield serve as an undulator. Electrons injected into the wakefield will experience undulator oscillations and emit bright radiation.

Electron dynamics and radiation

The transverse electric field experienced by the electrons and the trajectories of test electrons from 2D PIC simulations are presented in Fig. 2(a). In the simulation, the fundamental Gaussian mode (a 0 = 0.14) and the first-order Hermite-Gaussian mode (a 1 = 0.1) co-propagate in a matched plasma channel. Detailed parameters are given in the Methods. One can see that the transverse field is periodic with time and almost uniform in transverse coordinate x in a range of 2λ L , analogous to the distribution of the magnetic field in a magnetic undulator32. The trajectories of the test electrons are almost the same, although they are injected into the undulator with different initial transverse positions, which is completely different from the case of electrons undergoing betatron oscillations in a plasma wakefield8,9,10. As a result, the on-axis radiation spectra of the test electrons, shown in Fig. 2(b) are almost identical. The frequency of the first harmonic of the on-axis radiation is located at \(\omega =2{\gamma }_{0}^{2}{\omega }_{u}\mathrm{/(1}+{K}^{2}\mathrm{/2)}\sim 1900\,{\omega }_{L}\), where K = 0.44 is the strength parameter and ω u = c/Z R is the undulator frequency, corresponding to approximately millimeter undulator period. For λ L = 1 μm, the radiation wavelength is 0.5 nm, which is in the soft x-ray range. Theoretically, higher photon energies, even reaching the hard X-ray region, can be achieved by co-propagating even higher order modes together with the fundamental Gaussian mode. The undulator strength parameter K can be increased by changing the intensities of the laser pulses. In the above discussions the injected electrons are assumed to be locked at a certain phase. However, such condition can not be satisfied in most cases, because the injected high-energy electrons always run faster than the wakefield, which is referred to as dephasing in laser wakefield acceleration18. For an electron initially injected into the phase k p ζ = −5π/2, its phase will slip forward to k p ζ = −2π after a dephasing length L D . For a highly relativistic electron \({\upsilon }_{z}\simeq c\) traveling in a plasma channel, the dephasing length is given by L D = λ p /[4(1 − υ p /c)], where \({\upsilon }_{p}=c\sqrt{1-{\omega }_{p}^{2}/{\omega }_{L}^{2}-\mathrm{(4}{c}^{2})/({\omega }_{L}^{2}{w}_{0}^{2})}\) is approximately the phase velocity of the plasma wave, ω L is the laser frequency, ω p is the plasma frequency and λ p is the plasma wavelength. In the simulation presented in Fig. 2, the on-axis plasma density is n p0 = 0.001 n c and electron initial energy γ 0 = 1000, leading to the dephasing length L D ≈ 5000λ L . Dephasing sets a limitation to the number of undulator periods and, hence, reduces the radiation brightness. Moreover, varying of the electron energy broadens the radiation spectrum. This could be a limitation should this undulator be used for a compact FEL. In order to avoid dephasing and obtain radiation with high brightness and narrow-band spectrum, it is necessary to lock the phase of the electrons in the wakefield. In this work, phase-locking by longitudinally tapering the plasma channel is used (see Methods)33. As the electrons slip forward with respect to the driver laser pulses, the plasma density is increased, reducing the plasma wavelength and maintaining the phase of the electrons inside the plasma wave bucket.

Figure 2
figure 2

Distribution of the undulator field, test electron trajectories and their on-axis radiation spectrum. (a) The distribution of the transverse electric field normalized by mc 2 k p /e in the space (x, ct) for a fixed position k p ζ = −5π/2. The wakefield is created by two linearly polarized modes m 0 = 0 and m 1 = 1, with intensities a 0 = 0.14 and a 1 = 0.1, spot radius w 0 = 7λ L , and duration τ = 15T L , in a plasma channel with an on-axis density of 0.001n c , where \({n}_{c}={m}_{e}{\omega }_{L}^{2}\mathrm{/(4}\pi {e}^{2})\) is the critical plasma density with ω L the laser frequency. The blue lines represent the trajectories of the test electrons. (b) The corresponding on-axis radiation spectra for the test electrons with different initial transverse positions. Only the first harmonic is considered here.

We now investigate the dynamics of an electron beam propagating in the phase-locked plasma undulator. Parameters of the electron beam can be found in the Methods. The evolution of the beam is shown in Fig. 3(a). One can see that the beam keeps a constant radius while oscillating in the undulator with the wavelength λ u . Because the plasma channel has been tapered to lock the phase of the beam in the wakefield, the electron beam can stably propagate in the undulator for approximately 20000 λ L , corresponding to the number of undulator periods N u = 20. The radiation spectrum emitted by the beam is shown with the red solid line in Fig. 3(b). The spectrum bandwidth is as narrow as 6%, which is very close to the theoretical bandwidth \(\sqrt{{(1/{N}_{u})}^{2}+{(2{\rm{\Delta }}\gamma /{\gamma }_{0})}^{2}+{({\gamma }_{0}{\rm{\Delta }}{\theta })}^{4}/16}\) 34, where \({\rm{\Delta }}\gamma \simeq 2.35\,{\sigma }_{\gamma }\) is the FWHM electron beam energy spread, and \(\Delta \theta \simeq 2.35\,{\sigma }_{\theta }\) is the FWHM electron beam angular spread. Also minor betatron radiation located at ω = 400ω L can be seen in the spectrum, which comes from the electrons far from the axis where α β ≠ 0. The radiation of an electron beam in an untapered plasma channel is also presented with the blue dash-dot line in Fig. 3(b). The spectrum exhibits a two-peak structure with a larger bandwidth and a much lower brightness. This is because the un-phase-locked electrons oscillate in two different plasma wave buckets, then escape from the undulator with an angle and become lost finally. As we stated before, higher frequency radiation can be obtained by using higher order mode mixed with the fundamental mode. The green dashed line in Fig. 3(b) shows the radiation spectrum of the same electron beam in an undulator generated by modes m 0 = 0 and m 1 = 3. One can see that the central frequency of the radiation has increased to 5750 ω L , which is approaching the hard X-ray region. The bandwidth becomes even narrower because the electron beam experiences more undulator periods during the same propagation length.

Figure 3
figure 3

Trajectory and radiation spectrum of the electron beam propagating through the plasma undulator(s). (a) The trajectory of the beam in the plasma undulator. The electron density of the beam has been normalized. The parameters of the laser pulses and plasma are the same as in Fig. 2 except that the plasma density has been tapered to lock the phase of the electrons inside the plasma wakefield. (b) The on-axis radiation spectra of the electron beam with taper (red solid line) and without taper (blue dash-dot line) using modes m 0 = 0 and m 1 = 1. The green dashed line corresponds to the case of the modes m 0 = 0 and m 1 = 3 with taper. (c) The on-axis radiation spectra from a single undulator (red line) and two phase-matched undulators (magenta line).

It should be mentioned that tapering does not remove the mode slippage that appears due to different group velocities of the laser modes, and can become a main factor limiting the total undulator length. One can estimate the mode slippage length as the length it takes for two modes with numbers m 0 and m 1 to be separated by the longitudinal length L of each of the modes; hence, the mode slippage length is \({L}_{m,d}\simeq L{({k}_{L}{w}_{0})}^{2}/(2|{m}_{1}-{m}_{0}|)\). For the parameters of our simulations, this estimate gives the mode slippage length roughly equal to 14 mm, close to the 20 mm slippage length obtained from the numerical simulations. One solution for increasing the mode slippage length and extending the overlap of the laser modes is to employ longer pulses albeit with higher intensities. For the applications in FELs, the required number of periods is on the order of 1/ρ (\(\rho \sim 0.001\)), which is still much larger than that of the plasma undulator. A multiple-stage scheme could be considered to extend the number of periods. In Fig. 4(a) a schematic for a staged plasma undulator system is presented. The electron beam from the first undulator is transported to the second undulator by a discharge capillary, which acts as an active plasma lens25,28. An azimuthal focusing magnetic field is produced when an axial discharge current is introduced in the gas-filled capillary. It was reported that the field gradient could be larger than 3000 T/m25, enabling cm-scale focal length for GeV-level beam energies. The second undulator can be excited by the second pair of high-order laser modes reflected by a tape-based plasma mirror28. The trajectory of the electron beam in the two-stage plasma undulator obtained from PIC simulations is presented in Fig. 4(b). Each stage has 18 periods. Such a two-stage undulator also works in a 3D geometry, as shown by the typical trajectories of one electron in a linear undulator (Fig. 4(c)) and in a circular undulator (Fig. 4(d)). An important point that should be mentioned here is that in order to preserve the phase of the electrons with respect to the radiation field, one needs to carefully design the distance between the two undulators. For an electron with a gamma factor of γ 0, the phase shift after a length L in the radiation field is \(\delta {\phi }=L\mathrm{/(2}{\gamma }_{0}^{2})\). For proper phase matching δφ should be an integral multiple of the radiation wavelength λ. For our case here, γ 0 = 1000 and λ = 0.53 nm, the required distance is L = 1.06η mm, where η is an integer. In the simulation we choose L = 1.06 mm and get the on-axis radiation for an electron beam, as shown in Fig. 3(c). The radiation brightness in a two-stage phase-matched undulator is almost 4 times higher than in a single undulator. Also the bandwidth decreases from 6% to 3%. Such a multi-stage scheme provides the possibility to realize an FEL by employing many segments of plasma undulator.

Figure 4
figure 4

Staging of two plasma undulators using a plasma lens. (a) Scheme of two-stage undulator system. The electron beam from the first undulator is transported to the second undulator by a plasma lens. The second plasma undulator is created by the second pair of high-order laser modes which are injected via a plasma-mirror tape. The distance of the two undulators is carefully designed to preserve the optical phase of the electrons. (b) The trajectory of an electron beam in the two-stage plasma undulator. Each undulator has 18 periods. (c) A typical trajectory of an electron in a 3D two-stage linear undulator. (d) A typical trajectory of an electron in a 3D two-stage circular undulator.

Radiation polarization control

The polarization control of the X-ray is motivated by the applications in studying the dynamics of magnetization35 and polarization-dependent X-ray absorption spectroscopy36. For the present plasma undulator, the flexibility of the radiation polarization can be achieved by controlling the laser pulse phases. Consider three Hermite-Gaussian modes in a three dimensional geometry, |0, 0〉, |0, 1〉 and |1, 0〉, propagating in a matched plasma channel. The mixture of modes |0, 0〉 and |1, 0〉 makes the laser pulse oscillate in the x direction, while |0, 0〉 and |0, 1〉 leads to oscillations in the y direction. Different phases of these two orthogonal oscillations, determined by the phase difference of the two modes |0, 1〉 and |1, 0〉, and assuming that the amplitudes of these two modes are the same, control the polarization of the laser pulse centroid oscillation, and, hence, the polarization of the wakefield, electron oscillations, and radiation. Figure 5 shows the 3D PIC simulation results of the spatial distribution and polarization distribution of the radiation for the circular and linear polarization cases. Schematically the circular polarization case is presented in Fig. 1. When φ 01φ 10 = ±π/2, the far field of the radiation concentrates in a circular region with a radius less than one milliradian. Here φ 01 and φ 10 indicate the phase of the mode |0, 1〉 and |1, 0〉, respectively. According to Fig. 5(b), on-axis the radiation is almost perfectly circularly polarized, while it becomes elliptically polarized when the angle θ between the observer and the z axis increases. Finally, at θ = 1/γ 0, the radiation is linearly polarized.

Figure 5
figure 5

Distributions of the radiation intensity and polarization. The spatial distribution (a) and polarization distribution (b) of circularly polarized radiation generated when the phase difference between the modes |0,1〉 and |1,0〉 is π/2. (c) and (d) are the spatial distribution and polarization distribution, respectively, for the linearly polarized radiation, generated when the phase difference between the two modes is 0. In figures (a) and (c), the unit of the polar angle θ is milliradian, while the unit of the azimuthal angle ϕ is degree. In figures (b) and (d) the color-coded image represents the normalized energy-angular spectrum, while the contour lines represent the radiation ellipticity, with ellipticity equal to 1 corresponding to circular polarization and 0 corresponding to linear polarization.

When φ 01φ 10 = 0 or π, the radiation is linearly polarized, as shown in Fig. 5(d). The spatial distribution also demonstrates a linear structure in Fig. 5(c). In general, fine tunability of the distribution of the radiation intensity and polarization can be obtained by controlling the phases between the modes.

Beam loading limit and radiation brightness

A high current electron beam is beneficial for the FEL process15,32. For the plasma undulator, one needs to consider the beam loading limit since the electron beam will generate its own wakefield when propagating in the plasma37,38. This undesirable wakefield will affect the wakefield generated by the laser pulses, and can even destroy the undulator if it is strong enough. The beam-driven wakefield is zero at the head of the bunch and increases toward the tail. To avoid the breakdown of the undulator, the amplitude of the wakefield generated by the bunch should be smaller than that generated by the laser pulses at the bunch tail. For an electron bunch with Heaviside step-function profiles in radial and axial directions, the beam loading density limit, i.e., when the beam-driven wake amplitude equals the amplitude of the laser-driven transverse wakefield, can be written as (see Methods)

$${n}_{bl}=\frac{4\,C{a}_{u}}{{k}_{p}^{2}{w}_{0}{r}_{b}\mathrm{(1}-\,\cos \,{k}_{p}{L}_{b}){K}_{1}({k}_{p}{r}_{b}){I}_{1}({k}_{p}{r}_{b})}{n}_{p0},$$
(4)

where I 1 and K 1 are the first-order modified Bessel function of the first and second kind, n p0 is the on-axis plasma density, r b is the transverse radius of the bunch and L b is the longitudinal length of the bunch. As an example, in a plasma undulator created by a fundamental Gaussian mode and a first-order Hermite-Gaussian mode with a 1 = 0.1, a 0 = \(\sqrt{2{a}_{1}}\), C = 0.38, λ L = 1 μm, w 0 = 7 μm, r b = 0.5 μm, L b = 1 μm, n p0 = 0.001n c , the limit of the electron density n bl is 16n p0. Figure 6(a) plots different radiation spectra for different bunch densities. When the bunch density approaches the density limit, the undulator radiation becomes weaker as the undulator is suppressed. To keep a high performance of the plasma undulator, the applicable electron bunch density can be chosen as around 1% of n bl , as shown in Fig. 6(a). Theoretically, by the choice of the undulator parameters and electron bunch shaping, the current can approach kA limit. For instance, with the same parameters stated above but higher laser intensities of a 1 = 1.0 and a 0 = \(\sqrt{2}\), the 1% of the beam loading limit n bl reaches 0.008n c , which corresponds to an electron number of 0.6 × 107 and a current of 0.3 kA.

Figure 6
figure 6

Beam loading effect and brightness scaling. (a) The on-axis radiation spectra for different bunch densities: n b0 = 0.01 n p0 (red solid line), n b0 = 0.1 n p0 (blue dashed line), n b0 = 0.2 n p0 (black dotted line), n b0 = 1 n p0 (green dash-dot line). (b) The dependence of radiation peak brightness on laser intensity. In the calculation the electron bunch density has been chosen as 1% of the beam loading limit n bl .

From the simulation, one can obtain the peak energy radiated per unit solid angle per unit frequency interval d 2 I/dωdΩ = 10−29 J/Hz for one electron. Then the brightness of the radiation emitted by the electron beam in the simulation can be estimated as 1019 photons/(s mrad2 mm2 0.1%BW), while the electron bunch density here is 0.01 n p0 (see Methods for the calculation of the electron radiation). When the bunch density increases to 0.16n p0 (1% of n bl ), the radiation brightness increases to 7 × 1019 photons/(s mrad2 mm2 0.1% BW). The dependence of the brightness on the laser intensity is presented in Fig. 6(b), where the bunch density is chosen as 1% of n bl . With laser pulses with amplitudes a 1 = 1 and a 0 = \(\sqrt{2}\), the peak brightness of the undulator radiation is about 1022 photons/(s mrad2 mm2 0.1% BW).

Discussion

It is interesting to compare the main parameters of a permanent magnetic undulator, a plasma undulator and an RF undulator (see Table 1). Plasma undulator has a notably shorter period (about 1 mm or even less by using higher modes) as compared to the other two, whose period is typically more than 10 mm. The undulator strength parameter K of the plasma undulator can be tuned in a broad range by changing the laser normalized intensity a 0 of the laser pulse. The main disadvantage of a plasma undulator is the number of periods per segment which is around 20, while the number is on the order of 100 for a magnetic or an RF undulator. Thus a multi-stage scheme (Fig. 4) is very important for the FEL instability to develop. It should be mentioned that when we calculate the shortest wavelength of the radiation from a plasma undulator, we use the highest energy of electrons from an LPA at present30.

Table 1 Comparison between a magnet undulator, a plasma undulator and an RF undulator.

It is also very interesting to calculate the FEL parameters for the plasma undulator. The ideal 1D power gain length can be calculated by L g0 = λ u /(4π \(\sqrt{3{\rho }}\)), where \(\rho =1.78\times {10}^{-5}{A}_{u}^{\mathrm{2/3}}{\lambda }_{u}^{\mathrm{2/3}}[{\rm{cm}}]{n}_{e}^{\mathrm{1/3}}[{{\rm{cm}}}^{-3}]/{\gamma }_{0}\) is the Pierce parameter, A u = K[J 0(χ) − J 1(χ)]/\(\sqrt{2}\) with χ = K 2/(4 + 2K 2), K is the strength parameter and J 0, J 1 are the Bessel functions of the first kind. Using the parameters in the simulation, λ u = 0.97 mm, γ 0 = 1000, K = 0.44, n e = 0.001 × 1021 cm−3, we get ρ = 1.7 × 10−3 and L g0 = 27 mm ≈ 27λ u , which is less than the length of 2 segments of the present plasma undulator. It should be noted that many effects such as electron beam energy spread and emittance, space charge, finite bunch length and radiation diffraction will increase the gain length up to L g ≈ 2L g0 32. The power of the radiation P exp(z/L g ) becomes saturated after a saturation length L s ≈ 20L g . Therefore, the length of the present plasma undulator required to reach saturation is roughly L s ≈ 1 m, corresponding to N s ≈ 1000 undulator periods or 50 plasma undulator segments. We can then estimate the peak brightness of FEL radiation after a saturation length as 1022 × exp((L s L u )/L g ) ≈ 1030 photons/(s mrad2 mm2 0.1% BW), where L u is the length of one undulator.

The stability of a plasma undulator system is an important issue which should be considered in practical applications. Since a plasma undulator is essentially a wakefield generated by lasers, the stability of a plasma undulator is hence strongly dependent on the laser and plasma stability. The period of the plasma undulator is determined by the laser spot size while the strength parameter K is mainly determined by the laser power and plasma density. Note that the fluctuation of K caused by the plasma density error is much smaller than that caused by the laser error when the laser duration is around the optimized value, which can be deduced from the expression of C parameter in equation (1). Thus the fluctuation of the radiation fundamental wavelength can be expressed as

$$\frac{\delta \lambda }{\lambda }=\sqrt{{(2\frac{\delta {w}_{0}}{{w}_{0}})}^{2}+{(\frac{2{K}^{2}}{1+{K}^{2}/2}\frac{\delta P}{P})}^{2}},$$
(5)

where w 0 is the radius of the laser spot and P is the laser power. For a standalone undulator, its parameters are quite stable because δw 0 ≈ 0 and δP ≈ 0 when laser propagating in the underdense plasma. However, for a multiple-stage scheme, the matching between separate undulators becomes important. Technologically, each laser pulse driving an undulator should be generated by beam splitters from a single laser pulse, thus ensuring stability. For the FEL instability to develop, the fundamental wavelength has to be tuned with an accuracy given by δλ/λρ. Thus stable and precise control of beam splitter mirrors and focus mirrors is important for the FEL instability to occur, while shot-to-shot laser fluctuations are less critical as they only change the final radiation wavelength.

In conclusion, we have demonstrated a plasma undulator excited by high-order laser modes in a matched plasma channel, in which high-quality electron beams can make undulator oscillations with a few tens of cycles and emit bright X-ray radiation with a narrow bandwith. The advantages of the plasma undulator rely on matching the intensity relationship between the modes to suppress the betatron oscillation, and tapering the density of the plasma channel to lock the phase of the electron beam in the undulator. The polarization of the radiation can be controlled by changing the phase difference of the modes. The beam loading limit indicates that the tolerated beam charge can lead to currents as high as 0.3 kA, theoretically approaching the kA level. Such a plasma undulator, together with a laser-plasma accelerator, may open the way to realize an extremely compact FEL.

Methods

Calculation of the wakefield by two modes

The Hermite-Gaussian laser beams are guided in a matched plasma channel with a parabolic transverse profile as

$${n}_{p}(r)={n}_{p0}+{\rm{\Delta }}{n}_{c}\frac{{r}^{2}}{{w}_{0}^{2}},$$
(6)

where r is the transverse coordinate, n p0 is the on-axis electron density, \({\rm{\Delta }}{n}_{c}={(\pi {r}_{e}{w}_{0}^{2})}^{-1}\) is the critical channel depth, r e = e 2/m e c 2 is the classical electron radius18. In 2D, a guided Hermite-Gaussian mode can be expressed as31

$${a}_{(m)}(x,z)=\frac{{a}_{m}}{{(m{!2}^{m})}^{1/2}}{H}_{m}(\frac{\sqrt{2}x}{{w}_{0}})\exp [\frac{-{x}^{2}}{{w}_{0}^{2}}+i{\theta }_{m}],$$
(7)

where H m is the Hermite polynomial of order m and \({\theta }_{m}=(-z\mathrm{/2}\,{k}_{L})[{k}_{p}^{2}+\mathrm{2(2}\,m+\mathrm{1)/}{w}_{0}^{2}]\) is the phase, k L is the laser wave number and k p is the plasma wave number. From the phase θ m one can see that different modes propagate with different phase velocities. The phase speed for a specific mode m is \({v}_{ph}/c\simeq 1+[{k}_{p}^{2}+\mathrm{2(2}m+\mathrm{1)/}{w}_{0}^{2}\mathrm{]/2}{k}_{L}^{2}\). Consider two modes with the same linear polarization (y direction) and the same frequency, the total normalized intensity reads

$$I={a}_{\perp }{a}_{\perp }^{\ast }={|{a}_{y,{m}_{0}}|}^{2}+|{a}_{y,{m}_{1}}{|}^{2}+2|{a}_{y,{m}_{0}}||{a}_{y,{m}_{1}}|\,\cos [\frac{({m}_{0}-{m}_{1})z}{{Z}_{R}}].$$
(8)

One can see that the total intensity oscillates with a wavelength of 2πZ R /|m 0m 1|. The wakefield generated by the two modes will also oscillate while propagating in the matched channel. Assuming two modes with mode numbers m 0 (even) and m 1 (odd) and with the same Gaussian profiles in all dimensions (i.e., \(I\propto \exp [-{\zeta }^{2}\mathrm{/(2}{L}^{2})]\,\exp \,\,(-{r}^{2}/{w}_{0}^{2})\)) and \({a}_{\perp }{a}_{\perp }^{\ast } < 1\), the laser-excited, normalized wake potential is18

$${\varphi }(\zeta ,x,z)=-C{a}_{\perp }{a}_{\perp }^{\ast }\,\sin \,({k}_{p}\zeta ),$$
(9)

and

$$\vec{E}/{E}_{0}=-\nabla {\varphi }/{k}_{p},$$
(10)

yielding equation (1). Note that the second and higher order terms are omitted here, which implies that the area of interest is close to the axis (\({x}^{2}\ll {w}_{0}^{2}\)).

Particle-in-cell simulations

The simulations have been performed on JURECA39 at Jülich Supercomputing Centre using 2D/3D PIC code LAPINE40. Two y-polarized laser modes, the fundamental Gaussian mode and the first-order Hermite-Gaussian mode, propagate into a plasma channel from the left boundary. The normalized peak intensities are a 0 = 0.85 × 10−9 \(\sqrt{{I}_{L}({W/\mathrm{cm}}^{2}){l}_{L}({\mu }{\rm{m}})}\) = 0.14 for the fundamental Gaussian mode and a 1 = 0.1 for the first-order Hermite-Gaussian mode, respectively. The spot radius w 0 = 7λ L , duration FWHM (full width at half maximum) τ = 15 T L , laser wavelength λ L = 1 μm and laser period T L = 3.33 fs are the same for the both. The transverse profile of the plasma channel is designed according to equation (6), in order to guide the laser pulses. The plasma density along the propagation axis is n p0 = 0.001n c , where \({n}_{c}={m}_{e}{\omega }_{L}^{2}\mathrm{/(4}\pi {e}^{2})=1.1\times {10}^{21}{\lambda }_{L}(\mu {\rm{m}})/{\rm{c}}{m}^{3}\) is the critical density. Test electrons with γ 0 = 1000 are injected into the plasma channel at a longitudinal position with phase k p ζ = −5/2 π. For the electron beam, the transverse density distribution is \({n}_{b}(x)={n}_{b0}\,\exp (-{x}^{2}/2{\sigma }_{x}^{2})\) with σ x = 0.5 μm and n b0 = 0.01 n p0, and the density distribution in the momentum space is \({n}_{b}({p}_{x}/{m}_{e}c)={n}_{b0}\,\exp [-{({p}_{x}/{m}_{e}c)}^{2}\mathrm{/2}{\sigma }_{px}^{2}]\) with σ px = 0.05, corresponding to normalized emittance ε n = 0.025 μm. The rms energy spread of the electron beam is σ γ /γ 0 = 1% with γ 0 = 1000. The density distribution of the electron beam in longitudinal direction is uniform with a length of L b = 1 μm. In the 2D simulations, a moving window consisting of 600(x) × 1800(z) grids with 9 electrons and 9 protons per cell is used to follow the long-distance propagation. The resolutions are dx = 0.1 λ L and dz = 0.05 λ L . In the 3D simulations, the laser pulse mixture propagates only 2700 laser periods because of the limitation due the computing time. The size of the simulation box is 30λ L (x) × 30λ L (y) × 90λ L (z) corresponding to grids 150(x) × 150(y) × 900(z), with 8 macro-particles per cell for the background plasma and 64 macro-particles per cell for the injected electrons.

In the plasma channel, to maintain a constant phase in the transverse field the plasma density has been tapered as33

$$\frac{d\sqrt{{n}_{p0}/{n}_{c}}}{d(z/{\lambda }_{L})}=2\pi \frac{{n}_{p0}/{n}_{c}}{|{\psi }_{0}|}[\frac{1}{2{\pi }^{2}{({w}_{0}/{\lambda }_{L})}^{2}}-\frac{{n}_{p0}/{n}_{c}}{4}],$$
(11)

where ψ 0 = k p ζ(z = 0, t = 0) is the initial phase.

Calculation of the wakefield by an electron beam

According to the theory of beam-driven plasma wakefield37,38, the transverse field excited by an electron beam with a number density n b (r, ζ) = n b0 ψ(r)f(ζ) can be expressed as

$${W}_{\perp }(r,\zeta )=4\,\pi e{k}_{p}{n}_{b0}{\int }_{0}^{\zeta }d\zeta ^{\prime} {\int }_{0}^{\infty }r^{\prime} dr^{\prime} {\partial }_{r^{\prime} }\psi (r^{\prime} ){I}_{1}({k}_{p}{r}_{ < }){K}_{1}({k}_{p}{r}_{ > })f(\zeta ^{\prime} )\sin \,{k}_{p}(\zeta ^{\prime} -\zeta ),$$
(12)

where I 1 and K 1 are the first-order modified Bessel function of the first and second kind, ψ(r) and f(ζ) are the density profiles in radial and axial directions, respectively, r < = min(r, r ') and r > = max(r, r '). For a electron bunch with Heaviside step-function profiles in both radial and axial directions, \(\psi (r)= {\mathcal H} ({r}_{b}-r),f(\zeta )= {\mathcal H} (-\zeta ) {\mathcal H} (\zeta +{L}_{b})\), r b is the radius of the electron bunch and L b is the longitudinal length of the electron bunch, the radial wakefield in the body of the bunch (rr b , −L b ζ ≤ 0) can be written as

$${W}_{\perp }(r,\zeta )={E}_{0}{k}_{p}{r}_{b}\frac{{n}_{b0}}{{n}_{p0}}\mathrm{(1}-\,\cos \,{k}_{p}\zeta ){I}_{1}({k}_{p}r){K}_{1}({k}_{p}{r}_{b}\mathrm{).}$$
(13)

Let the amplitude of W equal to the amplitude of the transverse wakefield in equation (1), one can obtain the beam loading limit in equation (4).

Calculation of the electron radiation

The radiation is calculated using the trajectories from PIC simulations. We computed the radiation spectra in the far field by41

$$\frac{{d}^{2}I}{d\omega \,d{\rm{\Omega }}}=\frac{{e}^{2}}{16\,{\pi }^{3}{\varepsilon }_{0}c}{|{\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}\frac{\overrightarrow{n}\times [(\overrightarrow{n}-\overrightarrow{\beta })\times \dot{\overrightarrow{\beta }}]}{{(1-\overrightarrow{\beta }\overrightarrow{n})}^{2}}{e}^{i\omega (t-\overrightarrow{n}\overrightarrow{r}/c)}dt|}^{2}.$$
(14)

The radiation spectra of the electron beam is calculated by incoherent adding all the spectra of the 1100 macro-particles. The estimation of the radiation peak brightness [in unit of photons/(s mrad2 mm2 0.1% BW)] is as follows. We first multiply the term d 2 N ph /dω dΩ of one electron by the total electron number \({n}_{b}\pi {r}_{b}^{2}{L}_{b}\) in the bunch, then we divide it by the bunch duration L b /c and the sectional area \(\pi {r}_{b}^{2}\), and then multiply it by 0.1% BW and 10−6 due to conversion from square radians to squared milliradians.

Data availability

The data that support the findings of this study are available from the corresponding authors upon request.