Numerical growth of emittance in simulations of laser-wakefield acceleration

Transverse emittance is a crucial feature of laser-wakeﬁeld accelerators, yet accurately reproducing its value in numerical simulations remains challenging. It is shown here that, when the charge of the bunch exceeds a few tens of picocoulombs, particle-in-cell (PIC) simulations erroneously overestimate the emittance. This is mostly due the interaction of spurious Cherenkov radiation with the bunch, which leads to a steady growth of emittance during the simulation. A new computational scheme is proposed, which is free of spurious Cherenkov radiation. It can be easily implemented in existing PIC codes and leads to a substantial reduction of the emittance growth.


I. INTRODUCTION
Over the past ten years, laser-wakefield accelerators have considerably evolved and are now able to produce quasimonoenergetic electron beams [1][2][3] with energy up to 1 GeV [4]. This fast progress was made possible partly thanks to particle-in-cell (PIC) codes [5] and to their ability to reproduce or predict the properties of the experimental beams (e.g. total charge, mean energy, and energy dispersion).
Among the key properties of the beam, transverse emittance has attracted a growing interest lately, as several applications-including a prospective compact free electron laser [6]-require a very low emittance. Accordingly, considerable research effort is currently invested in decreasing the emittance of the beam, and unprecedently low values have indeed been measured in recent experiments [7][8][9][10]. In this context, an accurate calculation of emittance in PIC codes is essential in order to accompany and guide further research. PIC simulations have certainly been able to predict very low emittances (of the order of 0:1 mm mrad) [9,11], yet these occurrences are restricted to cases where the charge of the bunch is very low ( 1 pC). Here we demonstrate that, in the case of a bunch having a more typical charge [tens of picocoulombs (pC) or more], standard PIC codes tend to largely overestimate transverse emittance. This is due to the emittance spuriously growing from the beginning to the end of the acceleration. While numerical heating [12] is known to produce a similar growth, we find that, in our case, the dominating effect is the interaction of the bunch with numerical Cherenkov radiation.
Numerical Cherenkov radiation is a known effect in the PIC community [13]. It is associated with the finite difference time domain Yee scheme [14] and with the fact that, in this scheme, the numerical velocity of electromagnetic waves in vacuum is lower than the speed of light. As a result, macroparticles may travel faster than these waves and emit unphysical high-frequency Cherenkov radiation. Several methods have been proposed in other contexts, in order to mitigate this effect. These methods include lowpass filters [15][16][17] and modified computational scheme [16][17][18][19], for instance. However, none of these methods is ideally suited to the specificities of laser-wakefield simulations. We therefore propose a new computational scheme that is easily implemented in existing PIC codes, efficiently suppresses Cherenkov radiation, and considerably reduces the spurious growth of emittance.
The paper is organized as follows. In Sec. II, we consider a typical simulation of self-injection and report on the strong unphysical radiation that surrounds the electron bunch and the associated growth of emittance. Section III summarizes the main characteristics of Cherenkov radiation and shows that the radiation observed in the simulations is indeed of this type. Section IV proposes a computational scheme that prevents Cherenkov radiation. Finally, in Sec. V we discuss this new scheme and show that it leads to a reduced growth of emittance.

II. GROWING EMITTANCE IN SIMULATIONS
We begin by considering a typical simulation of laserwakefield acceleration. Parameters are chosen so as to lead to the self-injection of a few hundreds of pC of charge. More specifically, we consider a 1.5 J laser pulse ( ¼ 0:8 m) with 35 fs FWHM duration focused into a plasma with density n ¼ 6 Â 10 À3 n c ¼ 1:0 Â 10 19 cm À3 (where n c ¼ 1:75 Â 10 21 cm À3 is the critical density at this wavelength). The laser is focused so that its waist in the focal plane would have been 17 m in vacuum, corresponding to a 0 ¼ 2 (where a 0 is the peak value of the normalized vector potential). Notice that, due to selffocusing, the actual waist is slightly lower and a 0 is somewhat higher. The simulation was run using the fully 3D PIC code CALDER [20], which uses the Yee scheme and an iterative Poisson solver when solving the Maxwell equations. Calculations are performed in a moving window using 1600 Â 320 Â 320 grid points, with grid spacing Áx ¼ 0:032 m longitudinally and Áy ¼ Áz ¼ 0:25 m transversely, and using a time step cÁt ¼ 0:96Áx that is below the Courant-Friedrichs-Lewy (CFL) limit.
As the laser pulse propagates through the plasma, selfinjection occurs, leading to continuous trapping of background electrons. At a distance of 300 m from the point where injection started, the electron bunch has a total charge of about 250 pC, with a peak in the energy spectrum around 120 MeV containing about 150 pC (see the top panel of Fig. 1).
We observed that the emittance of the bunch steadily increases during the acceleration (see the lower panel of Fig. 1). In Fig. 1 and in the rest of this paper, the emittance calculated is the normalized emittance, which is defined, in the y direction for instance, as y ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where brackets denote an average over the electron bunch and where p y ¼ m _ y. It should be noted that in an fully evacuated plasma bubble [21,22], the transverse fields are linear functions of the radius, and as a result, normalized emittance should remain constant during the acceleration. On the contrary, the growth of emittance observed here is by no means negligible and reaches an average slope of about 6 mm mrad per mm of acceleration.
In the same simulation, we also observed that the bunch is surrounded by a strong high-frequency radiation. Figure 2 is a snapshot of the bubble after 300 m of acceleration and shows this radiation very clearly. Notice that this figure is a representation of E y À cB z instead of E y . The reason for this is that E y À cB z is proportional to the Lorentz force felt by the electrons F y ¼ ÀeðE y À v x B z þ v z B x Þ % ÀeðE y À cB z Þ. Moreover, a plot of E y would be dominated by the strong space-charge electric field of the bunch, whereas its effects are known to be almost canceled by the correspondingṽ ÂB term. Importantly, the same effects (growth of emittance and simultaneous radiation around the bunch) also appeared in simulations of colliding-pulse injection-although not shown here-in which the charge of the bunch was 50 pC and the emittance increased by 2 mm mrad per mm of acceleration. In both cases, the observed growth of emittance has serious consequences for the interpretation of the simulation, and it is of paramount importance to determine whether its origin is numerical or physical. Here we show that the growth of emittance is largely due to the observed radiation, which scatters electrons in the transverse direction by the means of the Lorentz force. This  Shaded regions correspond to zones of high electron density (in this case the sheath of the bubble and the accelerated bunch). While the laser pulse is clearly apparent on the right-hand side of the figure, one also notices higher-frequency radiation around the bunch. effect is indeed suggested by a close examination of the bunch in Fig. 2, which reveals transverse oscillations of its structure having the same spatial frequency as that of the radiation. The same effect will be confirmed in Sec. V, where the radiation is intentionally suppressed and the growth of emittance is observed to be slower.

III. NUMERICAL CHERENKOV RADIATION
It seems very likely that the radiation seen around the bunch is produced by the bunch itself, and a reasonable mechanism for this is the numerical Cherenkov effect. In this section, we summarize the theoretical properties of numerical Cherenkov radiation and we show that the radiation seen here is indeed of this type.
Cherenkov radiation [23]-whether physical or numerical-may appear when the dispersion relation of electromagnetic waves allows for some modes to travel with a phase velocity lower than the speed of light (v < c). In this case, energetic charged particles traveling close to the speed of light can excite some of those modes, thus producing a characteristic radiation. In the case of a charged particle traveling, e.g., along the x axis with speed c ( < 1), the modes that are excited are those that satisfy the relation where v ;x is the phase velocity of the mode along the x axis. In particular, v ;x ¼ !=k x for harmonic waves of the form expðik Ár À i!tÞ.
For electromagnetic waves traveling in vacuum the dispersion relation is ! 2 ¼ c 2k 2 , and therefore v ;x > c. As a result, Eq. (1) cannot be satisfied and no Cherenkov emission occurs. This is shown in Fig. 3 where the plot of v ;x and the plane corresponding to c do not intersect. However, due to spatial and temporal discretization in the Yee scheme, the actual numerical dispersion relation of these waves is and thus some values of ! andk may satisfy Eq. (1) in this case. This is shown in Fig. 4 for modes propagating in the x-y plane. The modes that are predicted to be excited by numerical Cherenkov effect are those for which the two FIG. 3. Physical phase velocity along x (v ;x , as deduced from ! 2 ¼ c 2k 2 ) as a function of k x and k y for electromagnetic modes propagating in the x-y plane (k z ¼ 0). The blue horizontal plane corresponds to v ;x ¼ c with ¼ 1 À 9 Â 10 À6 (which is representative of particles having an energy of 120 MeV).
In order to prove that the radiation observed in Sec. II is due to numerical Cherenkov effect, we verify that it corresponds to modes satisfying Eq. (3). To this end, we perform the spatial Fourier transform of the field E y À cB z in the x-y plane of our simulation and show the result in Fig. 5. As can be seen, the only modes with sizeable amplitude are exactly those that are predicted by Eq. (3) to be excited through numerical Cherenkov effect. We conclude that this radiation is indeed of numerical origin. We will now show that it can be removed by using a modified numerical scheme.

IV. A PROPOSED SCHEME TO AVOID NUMERICAL CHERENKOV RADIATION
A number of solutions have already been proposed in order to avoid numerical Cherenkov effect, albeit in different contexts (including for instance the case of boosted-frame simulations of laser-wakefield acceleration [17]). Yet, none of them is well suited to model laserwakefield acceleration in the laboratory frame.
For instance, some of these solutions [15][16][17] consist of using digital filtering, either in space or time, in order to damp high-frequency radiation. However, any filtering operation in real space (as opposed to Fourier space) that strongly damps high-frequency radiation inevitably also affects the physics at lower frequency. The effect of such a filter is particularly deleterious for physical phenomena having frequencies close to those that are to be filtered. In the case of laser-wakefield acceleration, the laser has a relatively high frequency, which makes it impossible to filter Cherenkov radiation without damping the laser itself-especially if the filter has to be applied at each time step and if the laser propagates over large distances. A notable exception to this rule is the case of boostedframe simulations [24], in which the laser-as seen in a proper relativistic frame-has a low spatial frequency [25], which then allows aggressive filtering.
Another set of solutions consist of modifying the discretization of the Maxwell equations, which then modifies the dispersion relation of electromagnetic modes. In [18,19], the numerical scheme is modified so that the CFL condition allows one to choose cÁt ¼ Áx. In this case, the dispersion relation is exact for modes propagating along the x axis, and Cherenkov radiation can potentially be suppressed. However, it was shown in [17] that using cÁt ¼ Áx triggers the onset of spurious numerical oscillations at the Nyquist frequency (k ¼ =Áx). Again, strong filtering is needed in this case, in order to remove these oscillations. Along the same line of thought, another set of numerical schemes was proposed in [16], and some of them can prevent numerical Cherenkov radiation without requiring cÁt ¼ Áx. These methods could in principle be successfully used in laser-wakefield simulation, yet they were only developed for isotropic grids (Áx ¼ Áy ¼ Áz), which limits their range of application. In this section, we propose a variation of the schemes proposed in [16], which is optimally adapted to simulations of laser-wakefield acceleration. Importantly, the proposed scheme is applicable to anisotropic grids and minimizes the impact of the modified dispersion relation on the propagation of the laser.

A. Description of the proposed scheme
In the proposed scheme, the fieldsẼ andB are defined on the same lattices as in the Yee scheme. The Maxwell equations take the following form: where D t andr symbolize the same discretized leap -frog operators as in the Yee scheme, i.e., for a field F n i;j;k defined on the nodes of the computational lattice, D t Fj nþ1=2 i;j;k ¼ðF nþ1 i;j;k ÀF n i;j;k Þ=Át and r x Fj n iþ1=2;j;k ¼ ðF n iþ1;j;k ÀF n i;j;k Þ=Áx. However,r Ã is a modified operator that depends on a set of coefficients, and which is defined by r Ã x Fj n iþ1=2;j;k ¼ x ðF n iþ1;j;k À F n i;j;k Þ=Áx þ x;y ðF n iþ1;jþ1;k À F n i;jþ1;k Þ=Áx þ x;y ðF n iþ1;jÀ1;k À F n i;jÀ1;k Þ=Áx þ x;z ðF n iþ1;j;kþ1 À F n i;j;kþ1 Þ=Áx þ x;z ðF n iþ1;j;kÀ1 À F n i;j;kÀ1 Þ=Áx þ x ðF n iþ2;j;k À F n iÀ1;j;k Þ=Áx (6) and by similar relations along the directions y and z. The position of the nodes that are used to calculate r Ã x Fj n iþ1=2;j;k are represented in Fig. 6. In order for the operator r Ã x to reduce to @ x in the limit Áx ! 0, the following relation has to be imposed: x ¼ 1 À 2 x;y À 2 x;z À 3 x : Similar relations also have to be imposed in the directions y and z. By injecting functions of the formẼ,B / expðikÁxÀ!tÞ into the discretized Maxwell equations, it is possible to show that the numerical dispersion relation in this scheme is with s t ¼ sinð!Át=2Þ and s u ¼ sinðk u Áu=2Þ for u ¼ x, y, z. Selecting an appropriate scheme now reduces to the choice of the and coefficients. In order to determine the coefficients that efficiently suppress Cherenkov radiation, we first impose that the phase velocity of modes propagating along the x axis is larger than c [v ;x ðk ¼ k xũx Þ ! c, for any k x , whereũ x is the unit vector along x]. From Eq. (8), the phase velocity of such modes is Thus, for a given Át and Áx, v ;x increases when x decreases (even for negative values of x ). In fact, v ;x is larger than c for any k x , when x x;0 , where Notice, incidentally, that x should not be chosen much below x;0 , since CFL-type instabilities start appearing for x < 1=4½1 À Áx 2 =ðcÁtÞ 2 . In the proposed algorithm, we opt for x ¼ x;0 . However, implementing x ¼ x;0 also affects modes having an obliquek, and makes them more prone to CFL instabilities. Adopting proper nonzero coefficients stabilizes these modes and efficiently prevents CFL instabilities. In the end, the proposed scheme corresponds to the following coefficients: As long as cÁt Áx, this scheme is CFL stable. Moreover, it satisfies v ;x ! c for all modes (see Fig. 7), thus preventing Cherenkov radiation for particles traveling along the x direction. We emphasize that this scheme is only efficient in the cases where relativistic particles all travel along the same axis. This is, however, clearly the case in simulations of laser-wakefield acceleration.

B. Spurious oscillations at the Nyquist frequency and velocity of the laser
We point out that, in the proposed scheme, the ratio cÁt=Áx is a free parameter (on which the coefficient x depends). Choosing cÁt=Áx ¼ 1 leads to x ¼ 0, and the scheme reduces to that of [19]. However, as observed in [17], this scheme develops spurious oscillations at the Nyquist frequency, which-unless filtered-rapidly grow to unacceptable levels. We observed that reducing cÁt=Áx in our scheme readily leads to a reduced level of oscillations. In practice, even for cÁt=Áx ¼ 0:96, these oscillations remain at a limited level instead of growing unbound, so that no filtering is required in the simulation.
Another important point to bear in mind is that the proposed scheme artificially increases the velocity of the laser-in the same way as the Yee scheme artificially reduces it. In the limit k laser Áx ( 1 and 1 À cÁt=Áx ( 1, Eq. (9) with x ¼ x;0 reduces to v ;x ðk laser Þ=c ¼ 1 þ where v g;x @! @k x ¼ @ @k x ðk x v ;x Þ is the group velocity of the laser along x.
The laser speedup here is of the same order of magnitude as the laser slowdown in the Yee scheme, since v Yee ;x ðk laser Þ=c ¼ 1 À 1=3 Â ð1 À cÁt=ÁxÞðk laser Áx=2Þ 2 and v Yee g;x ðk laser Þ=c ¼ 1 À ð1 À cÁt=ÁxÞðk laser Áx=2Þ 2 . When simulating laser-wakefield acceleration, it is important to ensure that the numerical alteration of the laser group velocity is negligible compared to the physical alteration caused by the plasma. In other terms, the criterion ð1 À cÁt=ÁxÞðk laser Áx=2Þ 2 ( ðn=n c Þ has to be satisfied. Notice that this holds for both schemes. Not satisfying the above criteria would lead, in the Yee scheme, to the injection of a spuriously high level of charge, while it would lead, in the proposed scheme, to an artificially low charge or no injection at all. Satisfying the above criterion implies that the laser wavelength has to be sufficiently resolved (i.e. k laser Áx sufficiently small) but also that cÁt=Áx should not be too different from 1. There is therefore a trade-off between cÁt=Áx far from 1, which keeps numerical oscillations low and cÁt=Áx close to 1, which reduces the alteration of the velocity of the laser.

V. SUPPRESSION OF CHERENKOV RADIATION
We implemented the scheme proposed in Sec. IV in CALDER and ran a simulation with the same parameters as in Sec. II. In particular, we used cÁt ¼ 0:96Áx. Figure 8 compares the field E y À cB z obtained in those two simulations, after 300 m of acceleration. Consistently with the analysis of Sec. IV, no Cherenkov radiation is emitted in the case of the modified scheme. We emphasize that this is indeed a consequence of the dispersion relation of the proposed scheme, since no filtering was used in the simulation. On the other hand, oscillations at the Nyquist frequency appear behind the accelerated bunch in the case of the modified scheme. Nonetheless, these oscillations seem to have much less impact on the bunch than the Cherenkov radiation does. This can be seen, in particular, in Fig. 9, where the evolution of the transverse emittance is plotted in both cases. Suppressing the Cherenkov radiation allowed us to reduce the growth rate of emittance by roughly a factor of 2. In simulations of controlled injection where the charge of the bunch was lower (' 50 pC), we were able to observe a reduction by a factor of 7. This strong reduction proves that Cherenkov radiation is indeed the main cause for the growth of emittance. It is unclear whether the remaining growth of emittance is the consequence of a physical effect (such as a partially evacuated bubble) or of another numerical artifact. The impact of the numerical scheme on the group velocity of the laser, which was mentioned in Sec. IV B, is also clearly visible in Fig. 8, where the laser pulse has a slight advance in the proposed scheme, as compared to the Yee scheme. In the two simulations, the laser pulse was initialized at the same position, but it propagated then over about 725 m before the snapshot from Fig. 8 was taken. From Eq. (13) and its counterpart for the Yee scheme, the expected difference in laser group velocity between the two schemes is Áv g;x =c ¼ 3ð1 À cÁt=ÁxÞðk laser Áx=2Þ 2 ¼ 1:9 Â 10 À3 . After 725 m of propagation, this leads to an expected advance of 1:4 m in the proposed scheme as compared to the Yee scheme, which is consistent with the observations in Fig. 8.
Apart from the differences mentioned above, our modified scheme did not affect significantly the physics in the simulation. In particular, the energy spectra of the electrons are very similar in the two simulations, as shown in Fig. 10, and the total charge of the bunch differs by less than 5% from one scheme to the other.

VI. CONCLUSION
In this paper, we showed that simulations of laserwakefield acceleration using standard PIC codes tend to largely overestimate the emittance of the bunch. In addition, we proved that this overestimation is, to a large extent, due to numerical Cherenkov radiation interacting with the bunch. By using a modified numerical scheme, we were able to suppress Cherenkov radiation and substantially reduce the calculated emittance.
The main implication of this work is that calculations of emittance, with bunches of a few tens of pC and more, require much care. The low values of emittance, that are now reached in experiments, are indeed extremely sensitive to the numerical artifacts of PIC simulations. Although we removed an important fraction of these artifacts, more work is needed in order to determine wether the obtained emittance is free of numerical errors.
Another important contribution of this paper is the numerical scheme that we propose. This scheme efficiently suppresses Cherenkov radiation and yet develops only a limited level of high-frequency numerical oscillations, as compared to the scheme used in [17]. As a result, valid simulations can be run with our scheme, without the need of low-pass filters. This algorithm, or modified versions of it, could in principle be used in other contexts where the dispersion relation is crucial. This potentially includes simulations of relativistic shock formation [26,27], as well as simulations of laser-wakefield acceleration in the boosted frame [24,28], since Cherenkov-type instabilities can represent a serious issue in both cases.