Essential advantage in elaboration of plasma accelerators [1] and plasma beam optics [2,3] stimulates the developing of other plasma elements able to reduce the size and cost of conventional devices. One of such instruments is the undulator [4], especially for XFEL generation [5], which needs of 10 s of meters of it to obtain the desired radiation. Adding plasma elements in its construction or even developing a full plasma undulator [6,7] may drastically improve the characteristics of these radiation devices operating at far shorter periods with the same K-number [4]. Both of these cases require thorough numerical investigation to understand the way how to develop, as the analytical approaches [6] are surely not enough, as one needs to carefully comprehend the intricate interactions between the electron beam particles and everything else.

There are several numerical softwares for simulating the characteristics of vacuum undulators [8–10]. However, simulation of complex setups like the plasma interaction with particle beams in the presence of an external fields can be performed only by particle-in-cell methods.

Since its inception in the 1960s [11,12], particle-in-cell (PIC) simulations have become an indispensable tool to study the underlying physics of complex phenomena which cannot be accessed in any other way. PIC permits to model the kinetics of a grand ensemble of particles interacting with each other and with electromagnetic fields in a self-consistent manner. In a PIC code, a finite number of macroparticles, i.e., groups of real particles that conserve the mass-to-charge ratio, sample the system phase space well enough to achieve a good results without the need of a large number of real particles. Such codes offer a great flexibility, thus, virtually any system could be simulated in an “ab initio” fashion. However, two major factors limit the applications of PIC. First, for highly complex systems, a lot of different types of phenomena can occur at the same time (ionization [13], collision [14,15], acceleration [16], etc.), which requires careful planning of the code and testing to make sure that everything occurs as expected. Regarding the code, the different numerical methods have distinct drawbacks (accuracy [17–19], charge conservation [20], border conditions [21], etc.). On top of that, testing that everything is correct needs of the creation of benchmark cases where the result should be understood beforehand, which, in some cases, is a challenge on itself. Second, the sheer number of macroparticles and field points defined in a fine spatial grid (in order to resolve the smallest elements) may demand huge amounts of computational resources for a 6D phase space. In some instances, the phase space can be reduced by just simulating a 2D geometry or taking advantage of symmetries like with the quasi-3D cylindrical codes [22,23]. However, such dimensional reduction can be only used in specific problems and even then, energy conservation, numerical heating [12], or the physics of the problem may be affected, especially for long simulations. So, the use of such tools requires a previous study to ascertain if they are appropriate and to what extent. Even using the dimensional reduction, problems that require long propagation distances (  >cm) are not practical to do due to the usual required resolution and simulation box size. For example, in principle, the simple propagation of a relativistic electron beam along a magnetic planar undulator to generate free electron laser (FEL) radiation [24], becomes a difficult task as it needs meters of propagation with a spatial resolution smaller than the electron beam radiation wavelength. A quick estimation tells us that for a radiation wavelength of 𝜆𝑟=200  nm, a proper simulation will require a spatial resolution of Δ⁢𝑥=𝜆𝑟16=12.5  nm, thus, a meter of propagation needs 8 ×107 steps; however, the bunching and initialization of coherent radiation for FEL [25] requires of typically ≈1–3  m, i.e.,  ≈8–24 ×107 steps. Simulating such number of steps for  ≈108 particles requires of huge resources including simulation time. In addition, the meter size simulation box cannot be calculated directly as it will require a prohibitive amount of RAM (  ≈1 PB) to be stored during the simulation. While the later issue can be solved by using a “moving window” [21] simulation box, the former is not trivial.

In the particular case of FEL generation, one might think that the already existing codes, e.g., GENESIS [8], SIMPLEX [9], ASTRA [10], are enough as they allow to predict the experimental FEL performance even though they are not “ab initio” and use multiple approximations. While this is true if one only cares about such usage, they are not enough to properly study the beam dynamics, and most importantly, the three element interaction previously alluded to, e.g., between the undulator field, electron beam, and a plasma. The addition of another element like a plasma to this process cannot be done with the aforementioned codes and it is of great interest for the future undulators as cm size plasmas show great promise as an electron beam early buncher. Furthermore, even the electron beam propagation inside 10 s of cm long plasmas cannot also be fully done in a reasonable time and set of resources, and such simulations are of capital importance for the research of plasma optics and plasma undulators, which require both long propagation and enough resolution for the radiation emitted.

With the objective of understanding the use of a plasma optic to improve the beam characteristics to increase the efficiency of the FEL radiation generation and exploring the real use possibilities different PIC existing solutions were explored; however, none adapted to such research was discovered. Consequently, this novel research required the development of appropriate simulation tools.

In this work, an extension of the PIC code FPlaser [26,27] to approach such problematic (arbitrary plasma-beam-external field interactions) is presented and used for our research, i.e., in the context of the interaction between a relativistic electron beam, a planar undulator, and a plasma. Then, the first results of the use of a plasma element in order to quickly imprint a regular longitudinal electron beam spatial modulation, i.e., microbunching, are obtained using PIC.

The code is based on the use of a booster reference frame that moves, in the electron beam propagation direction, at a constant relativistic speed, i.e., 𝛾𝑅 >1. In this relativistic reference frame, following the relativity equations, the meter (cm) size undulator (plasma) is reduced to ≈cm ( ≈mm) size. On the other hand, the μ⁢m Gaussian electron beam elongates to only mm size. These changes in distance make suddenly possible to use PIC to calculate the electron beam propagation even for meters without the need of excessive computer resources. However, being in a relativistic reference frame (RRF) instead of a laboratory reference frame (LRF) raises a lot of questions about the way to do and understand the simulation results of the complex three-way interaction that require of careful research.

Previously, such RRF has been tried in the code WARP and WARPX [28–32] for external injection LPA [33], interaction between a ultrarelativistic proton beam and simple uniform magnetic field [34]. However, multiple instabilities related to the backscattering radiation and at the plasma column that can end up introducing large quantities of noise are a recurring problem. These issues limit the 𝛾𝑅 possible for each simulation case. Nevertheless, great efforts are being carried to filter the instabilities manually using different methods and modified pushers [34,35], thus, allowing in some cases high 𝛾𝑅. In addition, a 2D FEL simulation was tried (low current electron beam in a periodic magnetic field) [36] in which no moving window was used and the electron beam was selected so the space charge was not prevalent and no interaction with the window borders appeared. However, while nothing inherently wrong was identified, the simulation could not properly reproduce the SASE (self-amplified spontaneous emission) phenomena and even with a prepared artificial seed the results differed from the expected results by a factor from 2 to 4 times. Some of the issues encountered were due to the choice of 𝛾. Around the same time, and following WARP’s, the OSIRIS code [37] also implemented a RRF for its PIC simulations. In this case, again, it has been mostly oriented for self-injection, injection, and external guiding and LPA simulations [38,39], and to simulate electron beam radiation during betatron oscillations in a ion channel like electric field [40]. Nevertheless, the same instability problems are found. The simulations are limited by the quantity of backward radiation and its wavelength [38,39], and in applicable cases, they are filtered. Moreover, field and plasma current smoothing are used to avoid instabilities in growth [38].

While not reinventing the PIC methods (e.g., leap-frog scheme, Poisson, etc.), these previous works are steps forward in the right direction to generalize PICs methods to phenomena that require of cm and meter sizes. However, as shown in the different previous works, numerous problems of the existing solution limit the utility and applicable circumstances and are not appropriate for every situation. This work is another step forward to expand the PIC RRF to complex systems, including proper ion motion, transformation of beam distribution from LRF to RRF, and simulation of high charge density electron beams. And while not applicable to all cases, the here presented code has been conceived for the accurate simulation of the beam dynamics.

Numerical Cherenkov radiation (NCR), which is caused by the discretization of space and time in simulations deriving in an inaccurate speed of light can cause instabilities and nonphysical phenomena [31,35,41–44]. While in this work a laser is not taken into account, it is important to remark that there is an essential difference between changing the reference frame for a particle beam (like is being done here) and a laser pulse. The RRF of the beam at rest ( 𝛾𝑅 =𝛾𝐵) does not exist for a laser, always moving at the speed of light. In the LRF, the most important parameter for the interaction between a plasma and a laser is 𝛿 =𝑁𝑒𝑁cr, with 𝑁𝑒 the plasma density and 𝑁cr =𝑚⁢𝜔24⁢𝜋⁢𝑒2 the laser critical density. If 𝛿 >1 the laser cannot propagate (overdense plasma) and if 𝛿 <1 the laser moves through the plasma (underdense plasma). Going to an RRF of 𝛾𝑅, the previous quantities are transformed as: 𝑁𝑒,RRF=𝛾𝑅⁢𝑁𝑒,LRF; 𝑚RRF=𝛾𝑅⁢𝑚LRF; 𝜔RRF=𝜔LRF2⁢𝛾𝑅. Therefore, the 𝛿 parameter becomes 𝛿RRF=4⁢𝛾2𝑅⁢𝛿LRF. For LWFA, the typical density of interest is 𝑁𝑒≈1019  cm−3 and for a Ti-Sapphire laser pulse 𝑁cr=1.7×1021  cm−3, thus, in the RRF, the plasma becomes overdense at 𝛾𝑅 ≈10. It is well known from PIC simulations, that, even in the LRF, a laser pulse propagation into plasma with 𝑁𝑒 close to 𝑁cr is unstable resulting in filamentation and soliton formation [45,46]. Meaning that, for a laser pulse the use of a RRF with higher gamma should result in some numerical instabilities related to the incorrect dispersion calculation, for example, NCR. Electron beams are free of this problem, and many instabilities vanish or affect only plasma after relaxation, thus, not the beam dynamics. This work also goes deeper about the discussion regarding the extraction of a full picture of the system in the LRF from RRF data, and how it is limited by the particles own momentum and the time resolution. Nonetheless, as in the mentioned previous works [29,33–35,38,39], the full transformation is not done here.

In this work, first the basis for such PIC code [47] is explained, how the RRF is created and some 2D results to confirm the correctness of the simulation. Then, to present the promising usefulness of such PIC code and the concept and feasibility of the plasma element to improve the beam parameters, thus, increasing the Pierce Parameter [48] (figure of merit for FEL) for later FEL radiation purposes, new results about the utilization of a low-density plasma as a tool to bunch the electron beam in cm lengths are shown; however, a full deep study will be made separately. In Sec. I, the general PIC scheme used and the implications of working in a relativistic frame will be presented. Sec. II deepens in the problematic of understanding the results in a RRF instead of a LRF. Then, in Sec. III, the code is tested with a vacuum and inside undulator beam propagation simulations and the results are compared with the expected results. Sec. IV presents an overview of new results about the bunching imprinted into an electron beam during its propagation in a low-density plasma, obtained thanks to the advantages offered by such code [47]. Before the conclusion, Sec. V swiftly shows the simulation of a complex system of a beam propagating inside a plasma inside an undulator.

The PIC code solves the electromagnetic equations via the Maxwell equations, constituting a system enough for the simulation:

Here 𝜌𝑖, 𝜌𝑒, and 𝐣𝑒, 𝐣𝑖 are the electron and ion charge densities and corresponding current densities. Typically, Eq. (4) is neglected in many problems, i.e., assuming the plasma ions are immovable. Equation (3) is solved by methods of current weighting [49]. In a reference frame moving with velocity 𝐯𝑅, the ion current density is not zero, it becomes 𝐣𝑖=𝐯𝑅⁢𝜌𝑖. Similarly, the electron current is also 𝐣𝑒 =𝐯𝑅⁢𝜌𝑒, and the total current is zero. However, with any charge separation under external forces the effect of plasma currents may become essential. One may exclude artificially 𝑣𝑅 from the electron part of the current calculation. However, the charge separation effects in this case vanish.

The simulation of the FEL process via PIC requires of meters of propagation, i.e., the order of magnitude for the electron beam to achieve bunching allowing the start of coherent radiation amplification is of (gain length [50]). Let’s consider a monoenergetic electron beam of FWHM length, in the LRF and an RRF moving with . Following Lorentz transformations [51]:

with the prime denoting the RRF, the longitudinal frame speed, and the speed of light. One finds that in , the undulator length is reduced by (Fig. 1). However, for a Gaussian initial electron beam, the transformation of the electron distribution using Eqs. (5)and (6) gives the following elongation of the beam:

with the beam longitudinal velocity, the RRF velocity, and the beam size on the LRF. Considering the distribution as done at the initial time of the RRF ( ) and the relation Eq. (7), the beam size becomes in the RRF:

with the electron beam gamma. Therefore, the aforementioned electron beam becomes in such RRF. These values permit to perform a reasonable PIC simulation of the entire process. As it is to be expected, the changes of lengths also affect the resolution needed. Since the beam propagation length in plasma and/or undulator is (with being or ) and plasma propagation length through the beam is the optimal is . However, if computer resources allow it, can be reduced.

FIG. 1.

Considering the FEL simulation case, the smaller feature that has to be resolved is, in the LRF, the radiation wavelength given by [52]:

with 𝜆𝑙 the radiation length, 𝜆𝑢 the undulator period, and 𝑘𝑢 the deflection parameter of the undulator. For a 1 T peak field and a 2 cm period undulator for our example electron beam, one obtains 𝜆𝑙 ≈304  nm. However, as the simulation is done in a RRF, the RRF 𝜆𝑙 adds a dependence in 𝛾𝑅 as follows [53]:

Thus, changing the required resolution on the LRF from nm to μ⁢m for 𝛾𝑅 >10, which is a reasonable spatial resolution for PIC, e.g., in our example 𝜆′𝑙 ≈6  μ⁢m. While it may be tempting to increase 𝛾𝑅 as much as possible, this will cause the electron beam to increase in length, therefore requiring a longer simulation window and number of cells to keep the μ⁢m resolution. A badly chosen simulation 𝛾𝑅 can increase too much the computational resources needed. The interplay between the different element sizes and their implication on the resolution and resources has to be carefully considered case by case. The resolution requirements can be summarized with a condition on the relativistic resolution factor 𝑅RF:

with 𝑁𝑥 the number of simulation cells in the longitudinal (propagation) direction and 𝐿wind the total longitudinal window size in the RRF. Furthermore, using Eq. (10) and simplifying (1+𝐾2𝑢2)≈1 (for undulator 𝐾2𝑢 <1), one can obtain:

The factor 8 of Eq. (14) comes from the fact that at least 8 points evenly distributed should be needed to properly determine a sinusoidal wave. The relativistic resolution factor 𝑅RF shows well the interplay of the changes on the different components of the simulation input parameters due to the relativistic nature of the RRF. Figure 2 presents the changes in 𝑅RF with respect to both 𝛾𝐵 and 𝛾𝑅.

FIG. 2.

Regarding the fields, due to the nature of the Maxwell equations, Eqs. (1)–(4) stay the same. However, the undulator dipolar magnetic field in the LRF is transformed as follows:

with 𝑢 the relative velocity between frames, ⊥ and ∥ the perpendicular and parallel directions with respect to the RRF propagation direction.

When adding a plasma, the same length contraction occurs as with the undulator. However, in the case of the plasma, as a consequence, the density ( ) in the RRF is higher:

Thus, depending on the chosen , if the plasma is too dense and thin a higher longitudinal resolution may be needed. Also to make sure that the plasma effect is well resolved for such cases, the number of macroparticles per cell for the plasma may also be increased.

The change of the plasma density raises the question about the plasma wavelength and wake size. Even though the wake behavior in plasma is out of scope of FEL topic, because the wake and its stability do not affect characteristics of beams, it is interesting to understand possible further plasma evolution. Since the plasma electron density and their mass increase with linearly, the plasma frequency does not depend on the reference frame. Therefore, should be also constant. However, just behind a driver, there is a wave moving with the driver group velocity, , and the density is a running wave in LRF. After the Lorentz transformation, the length of this wave increases times, in RRF. Nevertheless, after the driver passes through plasma, the running wave vanishes and should become independent on , with a peak density distribution inside the wave period depending on . In our test simulations, we observed the increase of wake length in first periods with , while the length of post density perturbation shows the length not increasing with . Nevertheless, the problem, which is not trivial, requires special consideration that is out of the scope of this paper.

All static elements in the LRF (e.g., undulator, plasma) move in the opposite direction to the RRF propagation with a speed of:

For the particle push, the leap-frog scheme [12] is here used, providing a second-order accuracy, a good compromise between exactitude and calculation speed.

The simulation uses a moving window carefully elaborated to work for any 𝛾𝑅, automatically adjusting its movement to be well adequate to the RRF time step. The boundaries of the problem are extremely important to avoid any artificial interference with the undulator and radiation fields and prevent particles from reflecting back into the simulation domain after reaching the edge of the grid. Regarding the fields, here an absorbing boundary condition [21] is used, i.e., the fields are absorbed when at a certain distance from the boundary and the absorption strength is nonlinear and proportional to the distance to the window edge. The integration of this boundary condition is similar to the LRF PIC case as it is independent of the reference frame and of the absorbed field as they still propagate at 𝑐. While the electron beam macroparticles become dummy particles when going beyond the boundaries. For all simulations in this work, the number of particles per cell is 4. The ion current is directly calculated (essential for plasma in the RRF) in this code. In addition, the approximation of using the electron current as an analog of the ions motion was tried. Results of both calculations were different, which implies that the ion current effect due to charge separation was essential. Therefore, the proper direct method was chosen. The particle populations (ion, electron, beam electrons) are split. The plasma ions and electrons use the same weighting, otherwise, the discrepancies cause instabilities.

Regarding the plasma, as mentioned before, in the RRF it moves opposite to the electron beam and window, therefore, it will always enter the window from one side and disappear through the other boundary.

The code itself started as a modified version of the extensively used FPLaser PIC code. However, the multiple challenges to adapt it for RRF made it clear that the changes required were numerous enough to justify the making of another branch of the code. FPLaser has been extensively tested and successfully compared with experiment and other codes for years, e.g., FBPIC [3], REMP [54], OSIRIS [55]. This relativistic code, when using 𝛾𝑅 =1, gives the same results as the nonrelativistic frame FPLaser (which already agrees well with other codes), thus, proving that the electromagnetic field calculation, current calculation, particle weighting, etc., are done correctly.

One of the difficulties for data treatment and comparisons between RRF and LRF that arise from the Lorentz transformations [Eqs. (5) and (6)] of the particles is the time difference. From Eqs. (5) and (6), it can be seen that as the electron evolves in the RRF, its momentum and position change, and therefore, the corresponding LRF 𝑡 to its 𝑡′. While straightforward for a single electron, when considering two electrons, one finds that for an instant 𝑡′, if the electrons have different position and momentum, their respective 𝑡𝑒⁢1 and 𝑡𝑒⁢2 are not necessarily equal, thus, the instant 𝑡′ translates into a span of time (Fig. 3):

FIG. 3.

The quantity of simultaneous electrons in a beam only aggravates this issue. Therefore, it is difficult to judge the exact state in the LRF of the full beam, especially during interaction with other elements, e.g., plasma. Two solutions present themselves to this conundrum.

For a group of particles of similar momentum in the RRF propagation direction, one can select 𝛾𝑅 and Δ⁢𝑋′ so that the Δ⁢𝑇 after transformation to the LRF is small enough to be negligible or at least acceptable:

However, the transformation from RRF to LRF carries always an error in the quantities of around the order of its variation during Δ⁢𝑇.

Due to the time quantization inside the simulation, the evolution of Δ⁢𝑇 corresponding to a step in 𝑡′ depends on both the simulation time step and the system particles momentum in the RRF propagation direction. Therefore, if the momentum in the RRF propagation direction for certain particles is significantly larger than others the Δ⁢𝑇 could vary substantially between time steps. Furthermore, it may be impossible to directly reconstruct a specific time in the LRF even with full information of all time steps during the RRF simulation. With the only direct solution being using a huge time resolution in the RRF, which could counter the simulation speed gains obtained over the LRF. That is why the condition of similar momentum particles was assumed for Eq. (22). Once each group dynamics have been transformed to the LRF, the time evolution in the LRF should be able to be reconstructed for a good resolution simulation. In particular cases, this LRF reconstruction process can be easier and simpler if the change in momentum is gradual and/or periodic, e.g., undulator field. Furthermore, for two different 𝛾𝑅 RRFs, let’s say 10 and 100, the time window produced by an instant in their respective reference frames will differ by a factor 10 (Fig. 3), for a same (Δ⁢𝑋′,Δ⁢𝑝′𝑥). So, even between RRFs the comparison is not necessarily intuitive.

The issue with the previous solution is the amount of postprocessing necessary to do it albeit being technically possible. The second solution is to base the analysis to only invariants. For example, the transverse dynamics are invariant between reference frames, radiated energy, beam charge in vacuum.

In addition, the use of 𝛾𝑅 =𝛾𝐵 difficults the spectral analysis as the beam is static in the RRF and the emission is in 2⁢𝜋. On top of this, for such high 𝛾𝐵 beams, 𝛾𝑅 =𝛾𝐵 reduces considerably the undulator period and increases the FEL wavelength, thus, having a bunching saturation wavelength much larger than an undulator period. Finally, the beam size increases substantially and the plasma is shortened causing sometimes the need for too high resolution. Therefore, for such studies 𝛾𝑅 =𝛾𝐵 does not seem to be appropriate.

These differences between reference frames also affect the understanding of the wakes generated by a driver in an RRF as already alluded previously. For a fixed time in the RRF the plasma particles are distributed in a Δ⁢𝑡LRF. So, when observing a wake in the RRF, one has to keep in mind that it does not correspond to any particular wake that may appear in the LRF, thus, being a different object, a wake in space time and not only space. Again, such counterintuitive subtleties make even more important the need for a careful postprocessing of the data.

Even though the RRF can allow for fast long simulations, each problematic, should have an “appropriate range” where the lower end offers speed ups and easier comparison with LRF and at higher ends, make possible quite long distances without still outweighing the computational gains by the need of a huge window and resolution, e.g., like in our case due to the beam becoming too large or the plasma too short and dense.

Our main motivation for the development of the RRF PIC code is to be able to simulate the electron beam dynamics during propagation inside long plasmas (10s of cms) for their use as prebuncher, optics, and also for the research of plasma undulators. Therefore, the most appropriate benchmarks are the propagation of electron beams in vacuum, in an external field (e.g., undulator) and inside a long plasma. All simulations are done in 2D as the main phenomena that need to be observed (i.e., plasma and planar undulator) would not meaningfully change from a 3D simulation.

The simulations parameters are presented in Table I. In all cases a beam of 30 pC, transverse diameter and length in the LRF has been used. After definition in the LRF the code itself calculates the equivalence to the RRF with the chosen following the prior equations. As previously mentioned, the electron beam length is increased [Eq. (9)]; however, the same total charge is kept constant, thus, making the charge density lower in the RRF. For all simulations, an initial transverse momentum has been applied to the beam. The time and space are normalized by 𝜔 and 𝑐/𝜔, respectively, with 𝑐 the speed of light, and 𝜔=𝜔plasma,cgs=2⁢𝑒⁢√𝜋⁢𝑁𝑒⁢𝛾𝑅𝑚𝑒 (with 𝑚𝑒 the electron mass and 𝑒 the electron charge), except in the case where 𝜔plasma,cgs is lower than 2⁢𝜋⁢𝛾𝑅⁢𝑐𝜆𝑢, in that case the latter is used instead for convenience.

In vacuum the only forces on the electron beam are those created by its own electrons electric field and current. In order to avoid any nonphysical currents due to the sudden calculation of the relativistic electron beam movement in vacuum during the first steps, an initial background of artificial ions perfectly compensating the electron beam charge is used. During a time 𝑇ion, the charge of these ions will be damped linearly until their total disappearance. This allows for a better initialization of the fields generated by the electron beam currents. 𝑇ion is here set as  ≈0.1  ps in the RRF. In addition, a very small transverse divergence (much smaller than that in realistic cases) is added to simulate a realistic case and also to avoid any nonphysical self-focusing of the electron beam due to Cherenkov instabilities [41]. The introduction of a transverse momentum reduces the normalized momentum, i.e., 𝑝𝑥𝑝total, thus, avoiding the numerical error that may allow electrons to go above the speed of light. Figure 4 shows the propagation of an electron beam with initial 𝜎′𝑦,rms=0.5  mrad (root mean squared) transverse divergence in vacuum during  ≈35.5  ps, which corresponds to  ≈0.55  m in the LRF, with the simulation parameters in Table I (V.). Figures 4(a) and 4(b) presents the initial beam density distribution and after 30.5 ps RRF of propagation and Fig. 4(c) the transverse density projection evolution at each time step. It is seen that the beam properly diverges as expected. After  ≈0.55  m the beam size grows up to 𝜎𝑦 =364  μ⁢m, which corresponds to the 𝜎𝑦 =315  μ⁢m due to initial divergence plus the effect of the beam space charge which increases slightly the divergence specially during the start of the simulation (highest charge density), approaching 𝜎′𝑦,rms≈0.75  mrad at 0.55 m. While an exact space charge calculation would require of 3D, due to the relatively low beam charge density of the beams here used, the differences between 2D and 3D are of no importance. Still this shows that a high resolution and accurate simulation of vacuum propagation along  ≈0.55  m taking into account the intrabeam forces, can be done in just 14 real time hours and 220 cores, which is already not possible in normal PIC codes.

FIG. 4.

Figure 5 presents a propagation of the electron beam inside an undulator [parameters at Table I (Und.)]. Special care has to be taken with the start of the undulator field in the simulation. The field cannot start suddenly as the electrons would jump from a zone without field to one with a strong enough field in a single step, causing a kick in the transverse direction. To solve it, a gentle slope is used, which has been found to erase this effect. The limitation comes from the leap frog scheme accuracy. This quick solution allows to prevent the nonphysical kick with minimum effect on the beam dynamics and avoiding the need of unnecessary high resolution. Inside a periodic undulator field the electron beam starts to oscillate [Fig. 5(a)] and eventually, due to its own charge and emitted radiation the density starts to acquire a longitudinal modulation [Fig. 5(b)] [56]. Following the evolution of the transversely integrated beam density [Fig. 5(c)], the slow apparition of the density bunches can be clearly seen. Initially, the electron beam starts its oscillation creating a density wave of wavelength , i.e., equal to the undulator (RRF). After around 8 ps in the RRF, the bunches start to appear [Fig. 5(d)], with an interbunch (i.e., distance between the distribution peaks of two consecutive bunches) distance of corresponding to . Following the 3D phenomenological equations for FEL obtained by Xie [57], for a self-amplified-spontaneous-emission (SASE [58]) FEL, the gain ( ) and saturation lengths ( ) for our configuration are ( in RRF) and ( in RRF), respectively. In the simulation, the beam starts to show clear differentiable bunches after , which agrees with the equations prediction ( or 6 ps in RRF). By the end of the simulation, i.e., , the mean bunch size (i.e., longitudinal size of single bunch inside the electron beam) starts to saturate, thus, it is close enough to the Ming Xie equations results for saturation ( in RRF). The saturation mean bunch size is of around [Fig. 5(d)]. The slight differences are easily explained by the phenomenological character of the equations and the simulation being done in 2D. A proper longitudinal density modulation can be observed in Fig. 5(c).

FIG. 5.

With time, due to the periodical change of longitudinal momentum of the electrons during their transverse oscillation, the beam propagates slightly slower than the simulation window causing its drift to the back of the window [Figs. 5(b) and 5(c)]. However, this can be easily solved by tweaking the window velocity.

Additional simulations were done with the individual elements (plasma, external field, beam) and their combinations for long distances and no numerical instabilities were observed.

As previously mentioned, the propagation of an electron beam inside a cm size plasma can be quite useful for its optics properties [3,59] and also as a prebuncher. In Fig. 6, the results of an 80 ps, in RRF (  ≈24  cm in LRF), simulation of a beam propagation through a low density 20 cm plasma [Table I (Pl.)] are displayed. Once inside the plasma, the head of the Gaussian beam generates a wakefield strong enough to cause a localized focusing on the back of the beam (bunch 1). At the same time, bunch 1 interacts with the wakefield and propagates through the beam toward the head [Fig. 6(c)] as the front part slightly defocuses, i.e., diverges [Figs. 6(a) and 6(b)]. Concurrently, when the bunch 1 achieves enough local density [Fig. 6(c)], it creates a wakefield that focuses a posterior part of the beam (bunch 2) while diverging in a half-moon shape. This bunching process due to the plasma that continues during the entire propagation is well illustrated in Fig. 6(c). The time between the apparition of a bunch and the creation of the next one takes  ≈7.7  ps ( std=0.7  ps), which is  ≈2  cm in the LRF. By the end of the 20 cm of plasma, a total of 11 bunches separated by 47.7±12.5  μ⁢m can be observed inside the beam [Figs. 6(a) and 6(c)], i.e., the back bunches are closer together than the front ones. In addition, even though a fishbone structure appears, transversely, the beam is focused inside the plasma by more than 1.4 times its initial density [Fig. 6(d)]. Once outside the plasma, the beam starts to diverge while keeping its structure.

FIG. 6.

Same parameter simulations have been carried with different 𝛾𝑅, and all of them show the fishbone structure due to the plasma.

The observed beam evolution in RRF is quite clear but to make sure the physics are correct, as it has not been done before, a simulation in the RRF 𝛾𝑅 =3 and its equivalent in FPlaser, i.e., 𝛾𝑅 =1, were carried [Table I ( RRFcomp and LRFcomp)]. Regarding the time difference between calculations, the 𝛾𝑅 =1 case needed close to 48 h to achieve the result shown at Fig. 7, while the 𝛾𝑅 =3 (Fig. 7) case arrived to a similar time position in only around 30 h but with  ≈3 times more resolution and a window 3 and 1.3 times bigger in the longitudinal and transverse directions, respectively. Both cases were performed in the same workstation of 128 cores. Due to the long time required to achieve some mm propagation in the LRF with just enough resolution, the simulation of 10 s of cm of plasma was not performed in LRF. Figure 7 presents the results of 𝛾𝑅 =3 and 𝛾𝑅 =1 at an approximately same LRF time of  ≈18  ps. It can be clearly seen that the previously described physical phenomena occurs in a similar fashion in the 𝛾𝑅 =1 simulation, exhibiting the focused node with an arc spawning from it forming the start of the fishbone structure. One can also see the back part starting to focus again which with time will become a third microbunch. Figure 7(a) displays the RRF result by transforming the longitudinal positions to the LRF compared to the 𝛾𝑅 =1 result. Due to the different longitudinal positions at the RRF, such direct transformation does not correspond to an exact 𝑇LRF but as the beam is short enough and 𝛾𝑅 low enough, it can be used for a comparison. Despite the difference in resolution and reference frame, a good agreement is found between both simulations, and as expected, the transverse evolution is equal even quantitatively between reference frames. A 1D density distribution of the central vertical position of both cases can be seen at Fig. 7(b). Both are quite similar, with some minor differences in the peak density value and rise edge. However, this can easily be due to the difference in resolution (the RRF higher resolution show smaller features) and the slight time difference of the RRF data (Fig. 3).

FIG. 7.

The code is also capable of higher energies and densities. Figure 8 presents the case of a 600 MeV beam in a 3×1016  cm−3 plasma [Table I (Pl.2)]. The initial Gaussian beam [Fig. 8(a)] enters the plasma and as in Fig. 6 starts to be focused but in this case the focusing is strong enough to expel electrons of the wake creating a plasma modulation [Fig. 8(b)]. Once out of the plasma, only the front half of the electron beam continues to propagate with a partial longitudinal density modulation while the rest stays trapped in the strongly modulated plasma [Fig. 8(c)]. For a previously prepared electron beam (longitudinally energy sorted beam, e.g., magnetic chicane [60]), this effect could be used as a low or high energies filter, reducing the beam energy spread.

FIG. 8.

The full high resolution simulation shown in Fig. 6 (Fig. 8) was made in only 50 (30) h with 150 (210) cores and less than 200 GB of RAM, e.g., a single work station, allowing the research of 10 s of cm simulation with a plasma. These simulation already present a great potential for the swift beam bunching using a plasma; however, an in depth exploration to properly understand the use of plasma as prebuncher using this code [47] is in progress and will be presented in another work as it is out of the scope of this article.

As a last demonstration of the capabilities of the code and future plasma prebuncher use case, Fig. 9 presents the case of an electron beam propagating inside a 20 cm plasma, which is at the same time inside a 1 T and 2 cm period undulator [Table I ( Pl.+Und.)]. The electron beam enters the undulator field slightly earlier than the plasma, thus, starting its oscillation [Fig. 9(a)]. Once inside the plasma [Fig. 9(b)], its modulation due to the wakefield starts as in Fig. 6. However, due to the oscillation, the wakefield does move in the transverse direction, which imprints a shift in the plasma wakefield induced bunches [Fig. 9(c)]. Again, an in-depth study of this phenomena is under progress.

FIG. 9.

In this work, a new PIC code that uses the advantages of space time dilation to allow the simulation of complex physical systems for long distances (up to meters) in Cartesian PIC using even a single work station has been presented. Then, using this tool the first steps toward the use of plasma elements to swiftly imprint a longitudinal microbunching to an electron beam for later use on an FEL has been achieved. The nonintuitive aspects of working in a relativistic reference frame (e.g., Gaussian beam elongation, fields transformation, relativistic Doppler) have been discussed. Due to the intrinsic features of the Lorentz transformation [Eqs. (5) and (6)], an instant 𝑡′ in the RRF can result in a time window of Δ⁢𝑡 in the LRF, therefore, making the interpretation of the code raw results nontrivial in several instances. Nevertheless, multiple solutions for the results postprocessing have been discussed; however, the appropriated steps to follow can heavily depend on the explored physics and cannot be all elucidated beforehand. Thanks to the code capacities, we have shown the simulation of an electron beam inside an undulator, with its longitudinal bunching agreeing with the widely used phenomenological FEL equations [57]. It also has been shown that the code is an appropriate tool for the study of plasma optics and prebuncher with excellent results, which needs of 10 s of cm of beam plasma interaction. The creation of bunches inside the electron beam due to the plasma wakefield generated by the beam itself shows a promising future for such device for uses in FEL and more. The insights about this problematic achieved by this code are quite useful and unique and will be properly explored in future works. Finally, the propagation of a beam inside a plasma and undulator has been simulated showing interesting features on the plasma generated bunches, e.g., the bunches shift in the oscillation direction. The usefulness of this PIC tool and capacity to unravel the new physics behind such plasma prebuncher by taking advantage of the relativistic reference frame in PIC has been demonstrated and plenty of new results will be obtained with it in the future allowing for an eventual experimental plasma prebuncher design and more. In addition, a 3D version is ready, extending the relativistic reference frame advantages of this code to simulations for 3D phenomena. We have to note that all these advantages are applicable only for charge particle beams. Simulations for laser pulses or other radiation sources requires special consideration to be done in general, specially its interaction with overdense plasma. In some other cases, the use of RRF may not even be suitable due to the resolution needs or too difficult understanding of the RRF with respect to the LRF.

This work used computational resources of the supercomputer Fugaku provided by RIKEN through the HPCI System Research Project (Project ID: hp240331).

D. O. E. and A. Z. developed the PIC code and the theoretical background. D. O. E. performed the simulations and analyzed and treated the data. The article was written by D. O. E. with discussions and corrections by A. Z. and A. R. additional discussions and support from M. M. and M. T.

The authors declare no competing interests.

The datasets used and/or analyzed during the current study available from the corresponding author on reasonable request.