Modeling classical and quantum radiation from laser-plasma accelerators

The development of models and the “Virtual Detector for Synchrotron Radiation” (vdsr) code that accurately describe the production of synchrotron radiation are described. These models and code are valid in the classical and linear (single-scattering) quantum regimes and are capable of describing radiation produced from laser-plasma accelerators (LPAs) through a variety of mechanisms including betatron radiation, undulator radiation, and Thomson/Compton scattering. Previous models of classical synchrotron radiation, such as those typically used for undulator radiation, are inadequate in describing the radiation spectra from electrons undergoing small numbers of oscillations. This is due to an improper treatment of a mathematical evaluation at the end points of an integration that leads to an unphysical plateau in the radiation spectrum at high frequencies, the magnitude of which increases as the number of oscillation periods decreases. This is important for betatron radiation from LPAs, in which the betatron strength parameter is large but the number of betatron periods is small. The code vdsr allows the radiation to be calculated in this regime by full integration over each electron trajectory, including end-point effects, and this code is used to calculate betatron radiation for cases of experimental interest. Radiation from Thomson scattering and Compton scattering is also studied with vdsr. For Thomson scattering, radiation reaction is included by using the Sokolov method for the calculation of the electron dynamics. For Compton scattering, quantum recoil effects are considered in vdsr by using Monte Carlo methods. The quantum calculation has been benchmarked with the classical calculation in a classical regime.

One of the promising applications of such stable high quality beams from LPAs is an incoherent radiation source in the ultraviolet and x-ray regime for which there are several mechanisms. Electrons accelerated in the wakefield produce betatron radiation if their trajectories are off axis [22]. Experiments have shown the emitted betatron radiation photons can be at the tens of keV level [23]. The LPA accelerated electron beams can also be sent to a magneticbased undulator to set up a synchrotron radiation source. Several groups have already carried out these experiments and nm level wavelength radiation was demonstrated [24,25]. Another short laser pulse (or, alternatively, a reflected pulse by a plasma mirror) can be used to produce high energy photons by scattering with the LPA electron beams [26][27][28]. Dependent on the beam energy, both Thomson or Compton scattering processes can occur. Besides working as a radiation source, the radiation can also be used as a diagnostic for the beam itself [29][30][31]. Because of the unique character of the electron beam in LPA, the beam size and internal structure inside the plasma are important and cannot be directly detected by normal diagnosis techniques. Betatron radiation gives a way to deduce this information.
Theoretical studies on these radiation mechanisms have been carried out for many years [22,32,33]. However, due to the unique character of the LPA and the accelerated beams (e.g., for betatron radiation, a large strength parameter a and small betatron oscillation number N ), conventional analytical results cannot be used directly. Furthermore, real world effects, such as the synchrotron strength parameter varying as a function of space and time (such as the case of realistic laser pulse profiles in Thomson scattering or acceleration effects in betatron radiation) make analytical calculations intractable. Numerical calculations are needed and they have been used recently to calculate the radiation spectrum. Previous numerical methods include photonrecording based calculation within particle-in-cell (PIC) simulations to obtain the betatron radiation spectrum [34][35][36], as well as a postprocessing code to obtain the radiation spectrum [37,38].
Previous formulations and numerical calculations of synchrotron radiation in the classical regime are inadequate to describe the radiation spectrum produced by electrons undergoing a small number of oscillations, as is the case in experiments on betatron radiation from LPAs. This inadequacy is the result of an improper treatment of a mathematical evaluation at the end points of the integration. This results in an unphysically large amount of radiation presented in the spectrum at high frequencies. We show that by properly including end-point effects, the radiation spectrum can be accurately calculated.
In this paper, we describe the development of models and the code ''Virtual Detector for Synchrotron Radiation'' (VDSR) to accurately calculate the synchrotron radiation from LPAs for a variety of physical processes, including undulator radiation, betatron radiation, and Thomson/ Compton scattering. These models and code are valid in the classical regime and in the quantum (linear, singlescattering regime) including the effects of radiation reaction and quantum recoil. Comparisons between the conventional analytical results and VDSR calculation have been made and the inadequacy of the analytic formulas are pointed out. The VDSR code is used to calculate betatron radiation using electron trajectories obtained from PIC simulations, and to calculate the radiation from beam-laser Thomson scattering and Compton scattering. For the last case, a benchmark of the quantum calculation is made by comparing the results with a classical calculation. This paper is organized as follows. In Sec. II, the inadequacy of the conventional formulas describing synchrotron radiation in the classical regime for a small number of particle oscillations is examined. By properly including end-point effects, the radiation spectrum can be accurately calculated. This is quantified by numerical calculations using the VDSR code. In Sec. III, several examples of synchrotron radiation calculations are presented, including betatron radiation from LPAs with 100 MeV level electron beams, Thomson scattering from LPAs with 300 MeV level electron beams, and Compton scattering from LPAs with 10 GeV level electron beams. Included is a benchmarking of the Compton models in VDSR via a comparison to the Thomson models in the classical regime. Section IV presents a summary and conclusion. This paper includes two appendices. The first presents the conventional analytical formulation of synchrotron radiation in the classical regime and the second presents details on the VDSR code.

II. INADEQUACY OF CONVENTIONAL FORMULAS
The conventional formulas that describe synchrotron radiation from a highly relativistic electron in the classical limit are summarized in Appendix A. These formulas are traditionally used to describe radiation from undulator magnets [32], betatron radiation from plasma focusing channels [22], and Thomson scattering when interacting with a counterstreaming laser pulse [33]. Analytical formulas presented were derived under several assumptions, including 2 z0 ) 1, 2 ( 1, and ð1 þ a 2 S =2Þ= 2 z0 ( 1, along with the assumption that the amplitude a S of the force giving rise to the synchrotron motion is constant, where z0 ¼ v z0 =c is the initial axial velocity of the electron, z0 ¼ ð1 À 2 z0 Þ À1=2 , and is the angle of observation with respect to the z axis. As is shown below, we find that the conventional formulas are only valid in the limit where the electron undergoes a large number of synchrotron oscillations, N S ) 1. For undulator radiation, this is typically well satisfied, since undulator magnetics are usually composed of a large number of periods, N u $ 100. However, for betatron radiation in a LPA, this condition (N ) 1) may not be satisfied, and for the parameters of many experiments of interest the number of betatron periods is small, N < 10.
In the limit N < 10, the conventional formulas predict an unphysically large amount of radiation at high frequencies.
To accurately model the radiation from a small number of oscillation periods, it is necessary that the end-point contributions to the classical radiation integrals be properly included, as discussed below.
To correctly describe synchrotron radiation for a small number of oscillation periods, we return to the most general form of the radiation spectrum d 2 I=d!d, which is the energy radiated per frequency, !, per solid angle, . From classical electromagnetic theory, the far field radiation spectrum of a charged particle is [39] whereñ is the unit vector pointing to the detector pixel, andr are the normalized velocity and spatial coordinates of the particle. We now assume that the radiation we are interested in is only produced in the time interval t 1 < t < t 1 þ T, during which the acceleration is nonzero, _ Þ 0. Outside of this interval, we assume the acceleration is zero. In this case we can replace the limits of integration in the above expression with t 1 and t 1 þ T. Note that, since we typically assume an initially highly relativistic electron with initial velocity z0 prior to entering the interaction region t 1 < t < t 1 þ T, this formulation intrinsically neglects the radiation produced by acceleration required to bring the electron up to its initial velocity z0 . However, we can assume that this initial acceleration prior to t 1 is adiabatically slow such that the radiation from t < t 1 is small and lies in a frequency region outside the one of interest. Hence, we assume and integrating by parts yields When the end-point term R 1 is neglected, the resulting expression yields the conventional synchrotron radiation formulas discussed in Appendix A (for the idealized cases in which the integral in R 2 can be performed analytically). For a large number of oscillation periods N ) 1, the second term (R 2 ) in Eq. (4), which is an integration along the particle trajectory, is usually the main contribution to the radiation spectrum. The first term (R 1 ) in Eq. (4) is usually omitted for a long oscillation trajectory, such as in the calculation for the undulator radiation, and this is typically a good approximation.
The end-point term R 1 , however, is important for the calculation of betatron radiation from the electrons oscillating inside a laser wakefield, in which the number of oscillation periods (N ¼ L= ) is often small ($ 1-10). If this term is omitted, there will be always a flattop region on the high frequency part of the spectrum. This is easy to see by the mathematical calculation of the (unphysical) radiation from a nonaccelerated particle during the time interval from t 1 to t 1 þ T. For simplicity, let us consider the radiation from a particle which has a constant velocity of and an initial position ofr 0 . One finds R 2 ¼ Ài! R t 1 þT t 1ñ Â ðñ ÂÞe i!½tÀñÁrðtÞ=c dt ¼ ÀR 1 , and jR 1 j 2 ¼ 2 2 sin 2 ½1 À cosðk ! !Þ=ð1 À cosÞ 2 , where is the angle betweenñ and and k ! ¼ Tð1 À cosÞ. As we see, if the first term R 1 is omitted, the spectrum will always have a flat region with a frequency oscillation of Á! ¼ 2=k ! ¼ 2=Tð1 À cosÞ. The amplitude of this term is negligible in the region of the main radiation frequency when N ) 1. However, for experiments on betatron radiation from a LPA, the effective radiation length of the particle is usually small such that N $ 1-10. In this case, the end-point term R 1 should be included, which is omitted in the conventional theoretical formulation.
To quantitatively examine the contribution of the terms R 1 and R 2 , we utilize the code VDSR that we have recently developed. The code VDSR is a parallel C++ code that can calculate the incoherent synchrotron radiation from electron beams for a wide variety of mechanisms, including undulator radiation, betatron radiation, and Thomson scattering. In VDSR, charged particle trajectories can be calculated by fourth-order Runge-Kutta method or be input from other files such as the output of particle-in-cell simulations. The radiation spectrum is then calculated for each electron by evaluating the expression Eq. (4) numerically, including both terms R 1 and R 2 . From Whittaker-Shannon-Nyquist sampling theory we know to get the high frequency radiation (! r ), the trajectories should at least have a resolution with a time step of Át < 1=! r , where ! r is the maximum radiated photon frequency. For a radiated photon of keV energy level, this means the trajectory step should less than 0.2 nm. Since the betatron radiation is typically emitted over cm scale, the integration steps should be 10 7 level. This is a significantly large computational cost. Fortunately in our case, the trajectories of the particles are relatively smooth within a relatively long scale time (such as a laser period or undulator period) which is much larger than the wavelength of the radiated photons. Over such a scale time, it is reasonable to think of the acceleration force as a constant and the particles as undergoing a uniform acceleration. The trajectories are then interpolated over this time step. Then the sampling frequency is infinitely small in principle, so the calculation can be carried out for keV to MeVand beyond photon radiation. The detailed integration and the limitations of this method have been described by Thomas [38] recently. The code VDSR adopts similar methods for the single particle radiation calculation. Further discussion on the code VDSR is given in Appendix B. Figure 1 shows the on-axis radiation spectrum obtained from the code VDSR for a single electron undergoing idealized betatron motion. The strength parameter of the betatron oscillation is held constant at a ¼ 2:087 and acceleration is neglected. Calculation results with or without the term R 1 for different N are shown. The asymptotic curve from Eq. (52) in Ref. [22] is also shown, which is equivalent to Eq. (A10) of Appendix A. As we see the asymptotic result fits the numerical results well when N is large and term 1 is included. All the spectra contain comb-like peaks and continuous or oscillating parts between peaks. The peaks are corresponding to the resonant and harmonic frequency of the betatron radiation The fine oscillations in each spectrum basically come from two parts: the end-points effect and numerical noise. They show different envelope profiles. The oscillations induced by end-points effect show a plateau envelope profile independent of frequency (see the yellow dashed lines). The numerical oscillation envelope decreases with frequency (see the white dashed line). All the calculated spectra in Fig. 1 are normalized by the maximum value to compare with the asymptotic line, which shows that including term 1 can reduce the high frequency oscillation a lot (not completely, due to numerical integration noise). To get rid of numerical noise, a higher accurate integration method (such as Fresnel integration, more computationally costly) should be used instead of the current Taylor expansion method (see the discussion in Appendix B and Ref. [38]).
Since the oscillation frequency in the plateau region (neglecting R 1 ) is much higher than the radiation frequency separation, the results for the radiation from a beam of electrons will deviate from the real value much worse than the single electron case due to the incoherent addition of the radiation from the collection of beam electrons. Hence, in the spectrum, the ratio of the plateau value to the radiation peak value in the beam case is much higher than the single electron case. The high energy part of the beam radiation spectrum will not be correctly calculated due to the omission of term 1 (R 1 ). On the contrary, if we include this term, the calculated spectrum is more accurate and is limited by the accuracy of the numerical integration. The real signal at the detector also includes the radiation beyond the time region of t 1 and t 1 þ T. However, a detector typically has a frequency window beyond which the radiated photons cannot be detected well. So the neglect of the radiation in those time regions is reasonable once we are sure that the radiation lies outside the relevant frequency region of the detector.

III. EXAMPLES OF VDSR CALCULATIONS
In this section we will give some examples of using the VDSR code. The simulation geometry is shown in Fig. 2. Three kinds of radiation have been calculated, i.e., radiation from betatron oscillation in a LPA, radiation from Thomson scattering, and radiation from Compton scattering by a short intense laser beam. Incoherent radiation from electron beam interacting with an undulator magnet can also be calculated once the external magnetic fields of the undulator are specified.

A. Betatron radiation from a LPA
When the electrons make betatron oscillations in the wakefield, they emit radiation. This can provide a source of broadband x rays based on a LPA [22,23]. The radiation intensity distribution in the code VDSR is calculated according to the following equation: where j represents different electrons. VDSR reads the electron trajectories from PIC simulations and then calculates the radiation from every electron according to their trajectories and sums them incoherently. The radiation spectrum has been simulated for a colliding pulse injection case [13]. The selected electron trajectories are coming from a particle-in-cell simulation using the (VORPAL) code [40]. The driving laser parameters are a 0 ¼ eA=m e c 2 ¼ 2:0, L FWHM ¼ 39:7 fs, and W FWHM ¼ 14:13 m. The plasma is composed of preionized electrons with a density of 1:7 Â 10 18 cm À3 . The injection is caused by a colliding pulse whose intensity is a 1 ¼ 0:3 and duration is L FWHM ' 10 fs. 1000 particles are used to represent the 1 pC accelerated electron beam. The statistical information of these particles is shown in Figs where the average value of hi ¼ 100 has been used. According to betatron radiation theory the peak radiation frequency on axis should be at @! c ¼ ' 0:003 47 eV, the radiation peak is at @! c ' 45 eV when ¼ 0. As we can see, this value is close to the on-axis peak photon energy in our simulation [see Fig. 3(a) at ¼ 0]. The simulation, however, shows the full radiation spectrum and angular dependence, which the theory does not. The angularly integrated spectrum shows the radiation peak is at the photon energy of about 90 eV. This means the higher energy radiation is not exactly on axis. However, most of the radiation is still within 10 mrad $ 1=hi.
The code has also been used as a diagnostic method to detect the beam transverse size inside the wakefield [29]. By using some experimental results such as final electron charge, energy spread, and betatron radiation spectrum, we set different beam sizes in the code and get the radiation on a faraway x-ray CCD camera. The raw radiation calculation has been convolved with experimental factors such as detector filters. Betatron oscillation length and transverse field are specified in the VDSR code based on the experimental and PIC simulation results. The final processed spectrum is then used to fit the diagnostic results in experiments. A large range of beam sizes was tested by the code to get the best fit. The beam emittance is then determined by using the beam radius inferred from the radiation spectrum and divergence from the magnetic spectrometer. The final deduced transverse beam size is 0:1 m for an electron beam of 463 MeV in a plasma with density of 5 Â 10 18 =cm 3 , which cannot be detected by normal diagnostic methods. Further work on beam diagnosis by using undulator radiation analysis is in progress.

B. Thomson scattering of LPA beam with a laser pulse
Thomson or Compton scattering by interaction of a LPA beam with a short intense laser pulse is another way to generate an ultrashort ultrahigh energy radiation source. Generally, Thomson scattering applies to photon emission with energies sufficiently below m e c 2 , such that the recoil of the electron can be neglected. Compton scattering applies to photon emission with energies near and beyond m e c 2 , such that the recoil of the electron is significant. Femtosecond or attosecond x-ray or -ray radiation could be obtained by this method. In the VDSR code, electron dynamics in a short laser pulse has been calculated by solving the Sokolov equation [41] to include the radiation reaction effects. A tightly focused laser pulse has also been implemented for describing the laser field, which includes high order modifications to the paraxial approximation [42].
Here we present Thomson scattering cases calculated using VDSR. The electron beam initially has an energy of 300 MeV, and charge is 160 pC. The beam is an ideal beam with zero transverse emittance at the beginning. Both the energy spread and beam rms radius are changed in different simulations. A laser pulse with 1 J energy and FWHM length of 13 ps is focused to a rms size of 25 m. A laser pulse with wavelength of 1 m propagates oppositely to the electron beam (colliding angle is 180 ). The two beams collide at the laser focus position except in Figs. 4(g) and 4(h), in which the colliding position is 4.968 mm away from the focus position. In the simulation we use 2400 simulation particles to represent the beam, the calculation convergence has been proved by using higher number of simulation particles. The time step for the trajectory calculation is t ¼ 0:02T 0 ' 0:066 fs, where T 0 is the laser period.
The spectrum distribution in the energy-angle space is shown in Figs. 4(a)-4(f) and 4(h). A typical radiation intensity distribution on the virtual CCD camera 4.7 m downstream from the laser focus is shown in Fig. 4(g). In the simulation we find totally about 4:5 Â 10 8 photons are projected to the virtual CCD. As we see the central energy of the emitted photons varies along with the interaction parameters. The rms energy spread of the emitted photons is 0.0165 MeV. When the electron beam has a small radius ( $ 1 m) and low energy spread (0.2%), the frequency at the radiation peak is about 1.55 MeV [see Figs. 4(a) and 4(b)], which is close to the theoretical prediction E r ¼ 4 2 E p = ð1 þ 2 2 þ a 2 L =2Þ ' 1:544 MeV, where E p ' 1:24 eV is the laser photon energy and a L ¼ eA=mc 2 ' 0:4628 is the normalized laser vector potential at the collision position.
Increasing the beam energy spread changes both the central frequency and frequency spread. For example, when the beam energy spread increases from 0.2% to 2%, the central radiation photon energy is at 1.565 MeV and the FWHM spectrum width in ¼ 0 direction is 0.11 MeV. When the beam radius increases from r b ¼ 1:0 m to r b ¼ 10 m, the central photon energy has changed to 1.5975 MeV and energy spread 0.145 MeV. When the colliding position has been changed to 4.968 mm away from the laser focus, the central radiation energy is increased to about 1.7 MeV, which is consistent with the value of E r with a L ¼ 0. This is because of the lower laser intensity at the colliding region. The figures also show the radiation is more concentrated in the forward direction in the electron motion plane [ ¼ 0 o , see

C. Compton scattering of LPA beam with a laser pulse
If the electron beam energy is higher, the radiated photon energy increases further. When the photon energy in the electron rest frame is close to the electron rest energy, the quantum recoil effects for the scattering process (Compton scattering) should be considered, which cannot be included in the classical formula. VDSR uses the Monte Carlo method to deal with the quantum scattering process. The method to simulate this process is described in Appendix B. Here we give two examples. The first one is for parameters where the classical theory applies. We use this case to benchmark the quantum calculation part of VDSR. Later we turn to the quantum regime and benchmark our codes with existing theories.
For the first case, a beam with energy of 100 MeV and 5% energy spread is collided with a laser pulse whose normalized intensity is a ¼ 0:1, wavelength is 0 ¼ 800 nm, duration is T ¼ 30T 0 (T 0 ¼ 0 =c), and the width at the focus is w ¼ 100 0 . The laser is linearly polarized and the electric field is in the direction of ¼ 0 . The classical calculation of the photon distribution is shown in Fig. 5 for the emitted photons in the plane of ¼ 0 (pink solid line) and ¼ 90 (light blue dotted line), respectively. For the quantum calculation, we collect the photon emission angle in two bins. The first one is À5 < < 5 , the emission intensity distribution is shown in Fig. 5 by the light blue solid line. As we see it is close to the classical calculation result. The quantum calculation has a much larger noise due to the low photon emission number. The second set of photon emission angles is 85 < < 95 (represented by the gray line) and it is compared with the classical calculation result ( ¼ 90 ). Similar results are given by the two calculation methods except the higher noise level in the quantum calculation. The averaged values agree very well between the quantum calculation (black solid line) and the classical calculation (green dotted line). Both of the two simulations use the same simulation particle number N ¼ 10 000. The time step of the classical calculation is dt ¼ 0:01T 0 . However, for the quantum calculation it is dt ¼ 1:0 which is only needed to resolve the laser envelope. The total simulation time for the classical calculation is about 4000 CPU hours; the quantum run only costs a few CPU minutes.
In the second case, we calculate the Compton scattering spectrum of a 10 GeV electron beam interaction with a laser pulse of a ¼ 0:1, 0 ¼ 1:0 m, L ¼ 30T 0 , and w ¼ 100 0 . The beam has a total charge of 1 pC, energy spread of 0.1%, and initially zero transverse emittance (p ? ¼ 0:0). The rms beam radius is 1 m. The centers of the beam and laser collide with each other at the laser focus. The quantum calculated radiated intensity distributions along with the photon energy is shown in Fig. 6(a). The different lines use a different number of simulation particles and different amplification factors (detailed in Appendix B). As we can see in a large simulation parameter range, the simulations give similar results except the noise level. The spectrum has a very sharp cut in the high frequency region, and the frequency at the peak is around 1.92 GeV, which agrees well with the quantum theoretical calculation [4 2 E p =ð1þ 2 2 þ4 2 E p =E e Þ'1:92 GeV]. Increasing simulation particle number and the amplification number can give better simulation statistics [see the  The signal is integrated in the direction. Here X ¼ ½ðp þ kÞ 2 À ðm e cÞ 2 =ðm e cÞ 2 ' 0:002 ( 1, quantum recoil effect is negligible and the classical calculation is applicable. In the quantum calculation an integration range for has been used. The averaged line means the average for the ¼ 0 and ¼ 90 . red solid curve in Fig. 6(a)]. Too large of an amplification number will result in multiple scattering for a single simulation particle, and broaden the final radiation spectrum. These simulation parameter effects on the resulting spectrum are discussed in detail in Appendix B. A classical calculation for a single electron scattering with the same laser pulse is also checked. The spectrum is shown in Fig. 6(b). As we see it gives a significantly different result. Both the frequency at the peak and the intensity are different. In the classical calculation, the radiation spectrum is centered at the frequency of 2.37 GeV, which is close to the classical theoretical prediction. We also check the radiation damping (reaction) effect on the radiation spectrum (red line). This makes the peak frequency redshift a little bit but peak energy is still far higher than the quantum results. The difference for the emitted photon redshift in the two kinds of calculation models is due to single photon quantum recoil effect, which is not included in a classical theory. Usually radiation reaction is important when a single electron emits lots of low energy photons, however, quantum recoil is important even if a single photon is emitted by an electron once the photon energy is high enough. In such a case only the quantum calculation can give correct modeling of the photon emission.

IV. SUMMARY AND CONCLUSIONS
We have presented the development of models and the code VDSR that accurately describe synchrotron radiation produced by a variety of methods, including betatron radiation, undulator radiation, and Thomson/Compton scattering. These models are valid in the classical regime and in the linear quantum (single-scattering) regime, including radiation reaction and recoil effects. The code VDSR calculates the synchrotron radiation using the trajectories of a collection of electrons in a beam, and these trajectories can either be calculated internally by VDSR based on electron motion in prescribed laser, plasma, and magnetic fields, or the trajectories can be imported from the output of a PIC simulation.
Previous formulations of synchrotron radiation in the classical regime are inadequate in describing the radiation emitted by a small number of oscillation periods. Specifically, previous models neglect the end-point term R 1 in Eq. (4), as does the conventional analytic formulas describing classical synchrotron radiation presented in Appendix A. With the end-point term R 1 neglected in Eq. (4), an unphysical plateau appears in the radiation spectrum at high frequencies, and the magnitude of this plateau increases as the number of oscillation periods N S decreases. By properly including the end-point term R 1 , as is done in VDSR, this plateau can be eliminated. This is particularly important for modeling betatron radiation, in which the betatron strength parameter is large a > 1 and the number of betatron periods is small N $ 1-10.
Various examples of the radiation spectra calculated from VDSR were presented, including betatron radiation from 100 MeV level electron beams, Thomson scattering from 300 MeV level electron beams, and Compton scattering from 10 GeV level electron beam. The quantum models in VDSR were benchmarked by comparison to calculations in the classical regime. The effect of simulation parameters, such as the simulation particle number and amplification factor, on the final calculated radiation spectrum were discussed. It was shown that by increasing both of these two parameters, a correct and lower noise spectrum can be obtained. Future development of the VDSR will include nonlinear Compton scattering and the modeling of other nonlinear QED phenomena.

APPENDIX A: CONVENTIONAL THEORY OF SYNCHROTRON RADIATION
In the classical limit in which the radiation reaction recoil of the electron can be neglected, the synchrotron radiation spectrum can be calculated using the standard methods of classical electromagnetic theory. In particular, the energy spectrum of the radiation emitted by a single electron on an arbitrary orbitrðtÞ andðtÞ can be calculated from the Lienard-Wiechert potentials [39], Note that this expression has neglected the end-point contributions that arise from the integration by parts used in deriving this expression. This integral can be performed analytically for the case in which the amplitude a S of the force giving rise to the synchrotron motion is constant, e.g., a S cosk S z for the case of undulator magnet, where a S is the normalized vector potential of the magnetic field and k S is the undulator period. In this case, the orbit of the electron can be written in a form similar tox ' ða S =k S z0 Þ sink S ct andz ' " z ct À ða 2 S =8k S 2 z0 Þ sin2k S ct, where the electron is initially assumed to be traveling along the z axis with normalized velocity z0 ; z0 ¼ ð1 À 2 z0 Þ À1=2 , and " z is the average axial velocity inside the medium causing the synchrotron motion (e.g., inside the undulator). Here the energy loss of the electron as it radiates is neglected. With orbits of this form, the radiation spectrum can be calculated with conventional techniques [22,32,33]. In the limits 2 z0 ) 1, 2 ( 1, and ð1 þ a 2 S =2Þ= 2 z0 ( 1, the radiation spectrum can be written as is the resonance function, where ! ¼ ck is the radiation frequency, ! n ¼ ck n is the resonant frequency for the nth harmonic, n is the harmonic number, a s is the synchrotron strength parameter, N S ¼ k S L=2 is the number of betatron periods that the electron undergoes, L ¼ cT is the interaction length, J m are Bessel functions, and f ¼ e 2 =@c ' 1=137 is the fine structure constant. Here is the observation angle with respect to the electron propagation axis and is the azimuthal observation angle. The above formulas are general expressions that can be used to describe synchrotron radiation produced by a variety of processes, such as with undulator magnets, betatron radiation from plasma focusing channels, and Thomson scattering from the interaction with a counterpropagation laser pulse. To describe undulator radiation, betatron radiation, or Thomson scattering, one needs to identify the proper synchrotron wave number k S , synchrotron strength parameter a S , and resonant frequency ! S , as discussed at the end of this Appendix.

Resonance function
The resonance function R n ð!; ! n Þ determines many properties of the radiation spectrum. Provided the number of synchrotron periods is large, N S ) 1, radiation is emitted in a series of harmonics and is confined in a narrow bandwidth about the resonant frequency (! ¼ ! n ) of each harmonic. The intrinsic frequency width Á! n of the spectrum R n about ! n is given by Á! n =! n ¼ 1=nN S . Furthermore, R n ð!; ! n Þ ! Á! n ð! À ! n Þ as N S ! 1. For a single harmonic n, the angular width Á I about the axis of a cone containing radiation with frequencies in a small bandwidth Á! about ! n is given by Â Á! n =! n for Á! Á! n ; Á!=! n for Á! ! Á! n : (A8)

On-axis radiation
Of particular interest is the radiation emitted along the axis, ¼ 0, where only the odd harmonics are finite, i.e., the even harmonics vanish. Setting ¼ 0 in the above expressions gives, for the nth odd harmonic, x ¼ 0, z ¼ n , and

Asymptotic behavior
For values of a 2 S ( 1, the emitted radiation will be narrowly peaked about the fundamental resonant frequency, ! 1 (n ¼ 1). As a S approaches unity, emitted radiation will appear at harmonics of the resonant frequency as well, ! n ¼ n! 1 . When a S ) 1, high harmonic (n ) 1) radiation is generated and the resulting synchrotron radiation spectrum consists of many closely spaced harmonics. Finite variations in the parameter a S and/or the electron beam parameters can broaden the line width and cause the spectrum to overlap. Hence, in the asymptotic limit, i.e., a S ) 1, the gross spectrum appears broadband, and a continuum of radiation is generated which extends out to a critical frequency, ! c , beyond which the radiation intensity diminishes.
Asymptotic properties of the radiation spectrum for a 2 S ) 1 and for large harmonic numbers, n ) 1, can be obtained with standard methods. In particular, the asymptotic spectrum along the axis ¼ 0 is given by where ¼ !=! c and ! c ' ð3a 3 S =4Þ! 1 is the critical frequency, which corresponds to a critical harmonic number of n c ¼ 3a 3 S =4. The function YðÞ ¼ 2 K 2 2=3 ðÞ is maximum at ¼ 1=2 and decreases rapidly for > 1. Half the total power is radiated at frequencies ! < ! c =2 and half at ! > ! c =2.

Undulator radiation
For an electron traveling along the axis of an idealized undulator magnet with a normalized vector potential of the form eA ? =m e c 2 ¼ K cosk u z, the synchrotron strength parameter is the undulator strength a S ¼ K and the resonant frequency is

Betatron radiation
For an electron traveling along the axis of an idealized plasma focusing channel (e.g., a plasma wakefield in which the accelerating field is neglected) with a radial focusing force that is a linear function of radius E r À B $ r, the synchrotron strength parameter is the betatron strength parameter a S ¼ a and the resonant frequency is Here a ¼ z0 k r , where r is the amplitude of the betatron oscillation (transverse electron orbit amplitude) and k is the wave number of the betatron oscillation, which is dependent on the strength of the radial focusing force, k 2 r $ E r À B . In the blowout (cavitation or bubble) regime, the focusing force is maximum and k ¼ k p =ð2 z0 Þ 1=2 , where k p ¼ ð4n e e 2 =m e c 2 Þ 1=2 is the plasma wave number with n e the ambient electron plasma density.

Thomson scattering
Consider an electron traveling along the z axis interacting with a counterpropagating, linearly polarized laser pulse with a normalized vector potential of the form eA ? =m e c 2 ¼ a 0 cosk 0 ðz þ ctÞ. In this case the synchrotron strength parameter is the laser strength parameter a S ¼ a 0 and the resonant frequency is Here ! 0 ¼ ck 0 is the laser frequency. The laser strength parameter a 0 is related to the laser intensity I L and wavelength 0 of a linearly polarized laser pulse by a 0 ' 8:5 Â 10 À10 ½mðI L ½W=cm 2 Þ 1=2 .

APPENDIX B: INTRODUCTION TO VDSR CODE
In this section, we describe the program structure and radiation calculation methods in the VDSR code.
The Virtual Detector for Synchrotron Radiation (VDSR) code is a parallel C++ code. It runs on the NERSC machine and uses thousands of processes well. It can be used as a PIC postprocess code for VORPAL [40] and VLPL [43] in a seamless transition style by using a script file combining the normal PIC run and postprocess VDSR run. The program framework is a parallel C++ code. To accurately simulate radiation from beams containing many particles, the code is parallelized to run on supercomputers such as NERSC and to efficiently use thousands of processors. Radiation calculation methods are implemented for classical and quantum radiation calculation, including radiation reaction.
The program flowchart is shown in Fig. 7. As one can see, except the global frame of the code (the initialization and final output parts), there are two different parts for classical and quantum calculation. However, since the code is object oriented the main structure of the code is uniform. Only the functions of the related classes are different when different models are used. The code mainly consists of the Beam, Particle, Detector, Pixel, and Field classes. Just as the meaning of the class names, each of the classes describes a single object. In the quantum calculation, a new class named Photon is additionally used to describe the radiated photons. The Domain and Control classes deal with the initialization, parallelization, and code efficiency diagnosis. Finally the photon information is projected to the detector object in the form of a radiation intensity distribution d 2 I=d!d. The output of the code includes field, beam, and radiation information. Spatial and temporal evolution of the field can be output, single particle and beam distribution in spatial and temporal space can be outputted. The radiation is outputted for each pixel (determined by , ) in the detector. The output files use a hierarchical data format (HDF5 library is used), which makes the post visualization and data process very convenient. The calculation model is specified in the input file, and all the particles use the same model to calculate the radiation. Work is in progress to allow the code to decide the calculation model for each particle automatically according to the interaction parameters.
Parallelization is made based on the particle tags. This means all the processes have the same virtual detector and only process the particles' trajectories and radiation assigned to them. This allows efficient parallelization because usually the particle number is very large and the detector information is relatively small. The final radiation of all the electrons is then incoherently accumulated to one detector and output. The Detector information specified in the input file includes the detector's spatial position, pixel number and distribution, detector frequency region and resolution.
For classical calculations, radiation is calculated by integration of Eq. (4). This requires knowledge of the electrons' trajectories. Trajectories in specified fields (which may include analytic laser and/or plasma fields) can be calculated by VDSR using a fourth-order Runge-Kutta method. In this case, the user specifies the fields and the initial phase space of the beam. Alternatively, trajectories may be read from a PIC code (such as VORPAL and VLPL). In these PIC codes, electron trajectories are output by a history function, which outputs the phase space information of a user-specified list of particles at every simulation step. For cases with high intensity fields where radiation reactions are important and cannot be neglected, only the first method can be used as the output trajectories from these PIC codes do not presently include the radiation reaction term.
The radiation from a single particle within a single time step is calculated using essentially the same methods used by Thomas [38]. Interpolations have been performed both for velocities and positions as shown in Fig. 7. The end-point effects have been used by adding the modified step integration term, which has been discussed in detail in Sec. II. The neglected terms [ R t 1 À1 Rðt;rÞdtþ R 1 t 1 þT Rðt;rÞdt] can be accounted in the code by selecting a velocity damping at the beginning and at the end within a short time period to make the particle be at rest beyond the simulation time region [44]. Different damping forms can affect the final spectrum when the radiation length is relatively small. For the integration along the particle trajectory, as shown in Ref. [38], the single step integration for particle j at time step k has the form: Àdt=2 e i2!ð1ÀñÁã 1 Þt ð j;k þb 1 tÞ Â e i2!ðÀñÁã 2 t 2 Þ dt; whereã 1 ,ã 2 ,b 1 are the coefficients of the interpolated positions and velocities as shown in Fig. 7. This integration can be changed to the form of standard Fresnel integrals; however, this needs much more computational cost. In VDSR, the integration has used two methods, the more accurate Fresnel integration and the Taylor expansion on the exponential function and truncate at order dt 2 : e i2!ðÀñÁã 2 t 2 Þ ' 1 À i2!ñ Áã 2 t 2 . The latter reduces simulation time greatly and is sufficiently accurate once !ñ Áã 2 t 2 ( 1, which is almost always satisfied in our simulation cases. In the quantum radiation calculation part, currently the code can only process laser scattering interactions. The particles' trajectories are calculated according to their momentum. The momentum evolution is calculated based on momentum conservation during the particle-photon scattering process. This means if the electron has not been scattered by a photon, it will drift without any momentum change even if it is inside a laser pulse. Since our scattering cross section is based on single photon-electron scattering, only linear Compton scattering can be calculated in our model. Our method is essentially the same as the methods used by Sun et al. [45]. Detailed information can be found in Ref. [45] and references therein. Here we describe the calculation process and the effects of simulation particle number and amplification factor, which have not been discussed in detail in Ref. [45]. As shown in Fig. 7 the particle's position (r) from the local laser intensity. The laser pulse has the form ofẼðtÞ / exp½Àðz À ctÞ 2 =L 2 Â sin½ðr; tÞ. The photon vector potential satisfies the distribution f / exp½Àðk À k 0 Þ 2 = 2 k , where 2 k ¼ 4=L 2 and k 0 ¼ ! 0 =c is the central wave number. The direction ofk is determined by the local spatial gradient of the laser phase, r. The total scattering cross section is calculated by where X ¼ ½ðp þ kÞ 2 À ðm e cÞ 2 =ðm e cÞ 2 is a Lorentz invariant which describes the recoil effect. When X ( 1, the recoil effect can be neglected, and t ' 8r 2 e =3 is the classical Thomson scattering cross section. Here p ¼ ðE e =c;pÞ, k ¼ ðE p =c; @kÞ are the standard four vector of particles and photons, respectively. When X is small, it is better to use the asymptotic formula of t to calculate the cross section, otherwise the error due to the limited machine precision of the computer will cause errors. The scattering possibility is calculated by P s ¼ t cð1 À Ák=jkjÞn p dtQ w f, where is the normalized velocity of the simulation particle, Q w is the weight of the simulation particle (the number of the real particles represented by a simulation particle, Q w ¼ N e =N s ), and f is an amplification factor to increase the sampling case. As we will discuss later, both Q w and f affect the simulation results. A random number P r which is uniformly distributed in [0,1] is used to compare with P s and to determine whether scattering happens (P r P s ) or not (P r > P s ). If no scattering happens, the particle makes a drift according to its current velocity and goes to the next time step. If the scattering case happens, a photon will be sampled according to the local laser field information. Both the particle's and photon's four vector are transformed to the scattering frame, in which the particle is at rest and the photon has a vector potential ofk ¼ k zẽz . Then, the scattered photon's momentum will be sampled according to the differential cross section as shown in Ref. [45] by the Monte Carlo method. After that, the scattered electron momenta are calculated according to the momentum conservation law. Later the four vectors of the scattered photon and particle are transformed back to the lab frame. The photon then will be recorded in a photon list and be judged whether to be projected to the virtual detector or not later. The particles' new position will be calculated according to their momenta. Since the total cross section has been multiplied by an amplification factor f, the intensity of the projected photon on the virtual detector should be divided by f.
The effects of the simulation particle number (N psim ) and amplification factor (f) effects on Compton scattering calculation are shown in Table I. In the table, a series of simulations using different simulation number and amplification factor are listed. The total radiated photon energy (E rad ), number of total radiated photons in the simulation (N phrad ), average radiated photon number per simulation particle (N phrad =N psim ), and per real particle numbers are shown (N phrad =N preal ). The last is also calculated by a classical formula and is about N phrad =N preal ' N K ' 0:0066 (see Ref. [22]). The reason to use the amplification factor f is just to increase the efficiency of the code and the number of scattering events, which will give a better accuracy for the Monte Carlo sampling to fit the scattering distribution when it is not practical to run the physical electron number, or when suppression of statistical fluctuation is desired. Since in our study case, the average fraction of real electrons which scatter is about 0.0066 (from classical estimation), if we use the real scattering cross section for the real number of particles, as we can see most of the simulation time is waste without any photon emission. In simulations one can use a simulation particle to represent many real particles, which can increase the simulation particle's cross section. In the case of Fig. 6(a), if we use 0.1 million simulation particles to represent 1 pC charge of electrons (Q w ' 62:5), the average radiation photon number per electron has been increased from 0.004 668 to 0.291 76, which is about 62.5 times (see simulation 6 in Table I). However, in some cases, this is still not enough to sufficiently increase the scattering possibility. To get a highly efficient simulation, amplification factor can be used to increase the scattering fraction TABLE I. Simulation particle number (N psim ) and amplification factor (f) effects on Compton scattering calculation. The beam and laser parameters are the same as the ones in Fig. 6 Table I, the N phrad =N preal has been increased about 10 times when f is 10 times larger). However, the amplification factor f should not be too much larger to result in an artificial multiscattering for a single particle. If that happens the emitted spectrum will be broadened due to the multiscattering for a single particle. The scattered particles lose some energy after the first scatter, and the second emitted photon energy from this particle will be smaller than the first case. This does not represent the real physics picture since the multiscattering is due to the artificial amplification factor. These effects are shown in simulation 7 in Table I, and the corresponding spectrum is shown in Fig. 6(a) by the blue dotted line. From the expression of the scattering possibility, we see that the particle weight Q w affects the simulation results. On one hand, too few simulation particles makes correct representation of the beam in phase space difficult. On the other hand, it will also increase the single simulation particle's scattering possibility and may result in multiscattering, which is shown in simulation 3 in Table I (compare with simulation 6). The radiation spectrum is similar to the blue dotted line in Fig. 6(a) with a broadened line. Increasing the simulation particle number certainly will increase the simulation cost. As we can see from Table I, increasing the particle number does not increase the final emitted simulation photons if f keeps the same, which means the statistic character of the spectrum cannot be improved by increasing the simulation number. This is because of the cancellation between the increase of the particle number and reduction of the scattering cross section. In fact, only the amplification factor affects the number of the scatterings. In the simulation, an appropriate particle number with an appropriate amplification factor can give reasonable results and optimal simulation cost. As shown in Table I, there is a large parameter region to get similar results. Usually increasing particle number and amplification factor simultaneously can lower the spectrum noise and keep the average scattering number less than one [see the red line in Fig. 6(a) and the corresponding simulation 8 in Table I].