Algorithm for calculating spectral intensity due to charged particles in arbitrary motion

An algorithm for calculating the spectral intensity of radiation due to the coherent addition of many particles with arbitrary trajectories is described. Direct numerical integration of the Lienard-Wiechert potentials, in the far-field, for extremely high photon energies and many particles is made computationally feasible by a mixed analytic and numerical method. Exact integrals of spectral intensity are made between discretely sampled trajectories, by assuming the space-time four-vector is a quadratic function of proper time. The integral Fourier transform of the trajectory with respect to time, the modulus squared of which comprises the spectral intensity, can then be formed by piecewise summation of exact integrals between discrete points. Because of this, the calculation is not restricted by discrete sampling bandwidth theory, and hence for smooth trajectories, time-steps many orders larger than the inverse of the frequency of interest can be taken.


Introduction
Radiation from synchrotrons is a well developed field, and a number of numerical methods for calculating the radiation spectrum exist [1,2,3,4,5]. Laser driven sources of radiation have also triggered interest in measurements [6,7,8,9,10,11,12,13,14] and numerical calculations [15,16,17,18] of radiation, particularly in the area of laser-plasma particle acceleration. As sources of high-energy particle beams and radiation, laser-plasma based techniques may be used for a large range of future applications. Ultrafast X-ray sources would be useful in, for example, time resolved diffraction, medical imaging, spectroscopy and microscopy of transient physical, chemical, or biological phenomena. Laser wakefield acceleration [19] of high energy electrons beams has recently successfully demonstrated the production of GeV peak energy electron beams [20], and has become a highly cited field of research. Such beams could be used directly for radiotherapy or radiographic imaging, or alternatively can be converted into fs duration, high brightness sources of x-rays. Numerical calculation of the x-ray spectrum is therefore of interest. Scattering of laser pulses from relativistic electron beams is another area which may require well characterized angularly resolved spectra from a realistic bunch, for comparison with experiment.
In general, the trajectories of particles in laser driven experiments are quite complicated; wakefields, electron beam interactions with intense laser fields, electron orbits in laser generated channels and in laser-solid interactions all represent sources of radiation. The radiation fields can be explicitly calculated by a fast Fourier transform method or by finite differencing of the Liénard-Wiechert fields in time, but these are computationally intensive processes due to the constraints of the Whittaker-Shannon-Nyquist sampling theorem [21,22]. This essentially states that for discretely sampled data, frequencies higher than half the sampling frequency are aliased to frequencies lower than half the sampling frequency. Since the spectral power as a function of frequency emitted is essentially equivalent to a Fourier transform, the highest resolvable frequency in the spectrum is constrained. However, for relativistic particles radiation can be produced at much higher frequencies than the actual time-scale corresponding to the change in momentum of the particles, because the radiation co-propagates with the particle.
Here, an algorithm is developed which uses a combination of numerical and analytic methods to integrate the Liénard-Wiechert fields, for arbitrary trajectories of many charged particles, to frequencies greatly exceeding the Nyquist frequency, ν N . This means that the sampling rate to accurately reproduce a particular spectrum can be many orders of magnitude lower than 1/ν N , and therefore faster to solve numerically. This is particularly relevant to calculations of radiation from betatron oscillations in laser wakefield accelerators and beam-laser interactions, but the technique is applicable to numerous other areas of physics.

The numerical method
The spectral intensity of radiation emitted by a number N P of accelerating point charges, with the jth particle at position r j and with normalized velocity β j = v j /c, can be expressed, in the far-field, as: where the unit vectorŝ is in the direction of observation, at a distance far compared with the scale of the emission region. This can be written alternatively in terms of proper time, τ : where the four-wave-vector κ α = ω{1,ŝ/c}, and v j is the momentum part of the jth particle's four-velocity defined as: One of the advantages of using proper time rather than 'laboratory' time for numerical calculations is that for a uniform step finite differencing scheme, the time resolution is effectively adaptive; as the particle gains inertia and is therefore accelerated at a decreased rate for a similar force, the laboratory time step-size increases. To numerically integrate the equations of motions for charged particles, both x α and v have to be recorded at a number of discrete points. To then perform the spectral integration numerically, a 'zeroth order' model would be to reduce equation 2 to a summation over finite differenced points: The problem with this method is that frequencies higher than half the sampling rate are aliased to lower frequencies, and therefore an upper limit is put on the maximum frequency that can be effectively resolved to ∼ 1/∆τ [21,22]. For attempts to simulate high energy photons from laser interactions, this can be computationally prohibitive. However, the motion of the particles which lead to such high-energy photons generally consists of changes on timescales much larger than the radiation frequencies produced. Hence, a different approach is taken, which is to assume that the motion of a particle between time-steps can be approximated by an interpolating function. The spectral integral can be expressed as a summation over analytic integrals between each time-step: Here the notation used is that Greek character sub/super-scripts represents components of 4-vector quantities, subscript Roman characters denote either particle number, j, or (proper) time-step, n, and bold font represents 3-vector quantities. Proper time is discretized into steps of size ∆τ . An exact analytic solution can then be employed to calculate the sub-integral between timesteps, using an interpolating function for v j and x α j between discrete time steps. A quadratic interpolation for x α can accurately model both linear and harmonic accelerations, and a linear interpolation for the velocity four-vector is consistent. Thus, the four-velocity and four-displacement at a time τ − τ n are approximated by: where x α 0n = x α (τ n ), v 0j,n = v j (τ n ), and v 1j,n , x α 1n and x α 2n are interpolation coefficients. Each integral over time in equation 5 can then be broken up into a series of analytic integrals between time steps from τ = n∆τ to τ = (n + 1)∆τ : where, with quadratic interpolation: and the space and velocity fourvectors are assumed real so thatĨ * j,n (ω) = I j,n (−ω). A change of variables in the integral leads to: 1j,n τ e ±iκαx α 2j,n τ 2 dτ . (9) If each particle trajectory is accurately described by the interpolation function between grid-points, then the total integral is exactly solved. This is the key difference which allows calculation of the spectrum to far beyond the Nyquist frequency corresponding to the time step between known position and moment values. The integral with respect to τ , I j,n (ω) = exp(−iκ α x α 0j,n )Ĩ j,n (ω), can be split into real and imaginary parts: It can be shown that these integrals have the solutions: where C(x) = x 0 cos(πt 2 /2)dt and S(x) = x 0 sin(πt 2 /2)dt are the Fresnel integrals, and: Here, the Fresnel integrals are solved numerically by using the power series and continued fraction expressions in Numerical Recipes in C++ [23].
Although the benefit gained in using less time steps by this method far outweighs the cost of calculating these functions, this is a computationally expensive process, and if the use of Fresnel integrals can be avoided it would be beneficial. If the second exponent in the integral is small, it is appropriate to Taylor expand the exponential function and truncate at order τ 2 , hence: The integrals when x α 2j,n τ 2 is small are: where: where sinc(x) is the unnormalized cardinal sine function, sinc(x) = sin(x)/x. The Taylor expanded solution for the integral I j,n is significantly faster to solve than the Fresnel integral expressions, so a condition statement can switch between solutions depending on the size of χ 2j,n ∆τ 2 . For practical purposes, the expansion solution is necessary because the exact solution, when numerically solved, diverges for small χ 2j,n ∆τ 2 due to floating point truncation error, the accuracy of the trigonometric functions and the inverse χ 2j,n ∆τ 2 relationship. Some threshold value, T , for χ 2j,n ∆τ 2 can be used to choose between the two models. This must be between 1 and 0, as the Taylor expansion is only valid for χ 2j,n ∆τ 2 < 1. For values less than the threshold, χ 2j,n ∆τ 2 < T , the Taylor expansion solution is used and for values greater than threshold, χ 2j,n ∆τ 2 > T , the exact solution is used. In addition, long double precision floating point numbers are necessary for sufficiently accurate calculation of the functions. A good threshold value is T = 10 −3 , as a balance between the accuracy of the solution and the speed of the algorithm. Figure 1 shows calculations of integral I j,n for dτ = 1, v 0j,n = 1, v 1j,n = 2, χ 1j,n = 1 and varying χ 2j,n . (a) Calculation of integral using Fresnel integral form. (b) Calculation using Taylor series expansion form of integral. (c) Calculation using finite differenced integral with 10 7 steps and Simpson's rule. It is clear in this specific example that the Fresnel integral solution is correct for all but the smallest values of χ 2j,n . For these values, since the equations 12 and 13 involve the subtraction of large terms and division by χ 2j,n , the truncation errors are amplified to significant values. The Taylor expanded solution does not contain a term 1/χ 2j,n , and therefore approaches the correct result for χ 2j,n 1.
!" The radiated spectral intensity is then given by: The radiation would normally be considered in a spherical polar coordinate system, {r, θ, φ}, where θ is the azimuth and φ is the polar angle. For sim-plicity, a cartesian coordinate system, x, y, z, can be chosen such that the radiation is calculated at an angle θ with respect to the z axis in the y − z plane. To observe radiation in a particular direction in φ, the coordinate system can be rotated about the z axis. Defining x, y, z components of combinations of integrals I j,n : the radiated spectral intensity is: (29)

Synchrotron radiation
A test for the algorithm for relativistic motion is reproducing the well known continuum spectrum from a single ultra-relativistic electron rotating in a magnetic field without losing energy. Mathematically this means no radiation damping term in the force equation, physically this can be accelerating electric fields that compensate for the energy losses. An analytic expression for the spectrum is the sum of modified Bessel functions of the second kind, K ν (x): where ξ = ωρ/3cγ 3 (1 + θ 2 γ 2 ) 3/2 and ρ is the radius of curvature. Figure 2 shows calculated synchrotron spectra (left) and the percentage error relative to the analytic solution (right) for different numerical integration techniques. The relative percentage error is defined as the difference between the analytic and numerical solutions divided by the maximum value of the analytic solution. The electron Lorentz factor is γ = 1000, and the magnetic field strength, B, is such that eB/m e ω 0 = 1. The time-step is ω 0 ∆τ = π × 10 −4 , and the maximum frequency calculated is 10 7 ω 0 . Note that this means that the maximum frequency resolved is π × 10 3 times the inverse of the time-step size. (a-b) show the calculated spectrum and error for finite differenced numerical integration, (c-d) show the spectrum and error for the summation of exact integrals using first-order interpolation of position only and (e-f) show the spectrum and error for the summation of exact integrals using second-order interpolation of position. In the left-hand images, the red dashed line corresponds to the analytic solution, and the blue solid line corresponds to the numerically calculated solution. It is clear from (a-b) that using the simple finite differencing method (equivalent to a discrete Fourier transform) that the resulting spectrum is of no use whatsoever with this particular time-step. The error grows with frequency to be 12 orders of magnitude larger than the maximum of the spectral intensity that the method is trying to reproduce. Decreasing the time-step size eventually yields an error smaller than 100%, but to produce accurate spectra, ω 0 ∆τ 1. This is a very limiting factor in calculations of this kind.
By performing a first-order interpolation, as in figure 2 (c-d), which corresponds to using the real part of the Taylor expanded solution only, equation 21, the calculated spectrum is clearly now representative of the analytic spectrum, albeit with an error of up to 10%. The simplicity of this algorithm (the only non algebraic operation involving the calculation of a single sinusoid for each time-step) makes it fast and easy to implement. However, beyond 10 7 ω 0 , the error in the numerical solution starts to grow. Using second-order interpolation, (e-f), yields an error that is at most more than an order of magnitude smaller than the maximum error in the first-order interpolation calculation. Importantly, it is free of spurious oscillations which could lead to misinterpretation of more complex spectra. In addition, the error remains less than the maximum error shown in figure 2(f) up to 5 × 10 9 ω 0 , which corresponds to a frequency of ω > 10 6 /∆τ . Despite being a more complex algorithm, the second-order method greatly speeds up a spectral calculation due to a larger time-step being allowed. Provided the motion of a particle can be well-described by piecewise quadratic functions (i.e. cubic or higher terms would be small), the method will produce accurate spectra.
In figure 3, angularly resolved polarization components of synchrotron spectra for the same parameters are shown. In this case, the frequency scale is shown as an energy scale corresponding to a magnetic field of B 0 = 13 T, which is similar to the parameters of a small synchrotron. The two polarized components of the radiation are u ⊥ = µ 0 e 2 c/16π 3 ω 2 ( (S x ) 2 + (S x ) 2 ) and u = µ 0 e 2 c/16π 3 ω 2 ( (S y ) cos θ − (S z ) sin θ) 2 +( (S y ) cos θ − (S z ) sin θ) 2 , i.e. the components perpendicular and parallel to the plane of observation. Such calculations can be performed in a matter of minutes on a single processor (the figure shows 500 frequency bins by 500 angular bins resolution and took 30 minutes to calculate on a single processor of a 2 × 1.4 GHz iMac, but a calculation of 100 frequency bins by 100 angular bins resolution took ∼1 minute, including the calculation of the particle trajectory.)

Linear and non-linear Thomson scattering
To test non-relativistic motion extending to moderately relativistic motion, Thomson scattering of a laser pulse from a single electron can be analyzed. In non-relativistic Thomson scattering from a linearly polarized plane wave, an analytic solution for the radiated power is possible. The field strength corresponding to the interaction is parameterized by the normalized peak vector potential a 0 = eE 0 /mcω 0 , where E 0 is the peak electric field strength and ω 0 is the angular frequency of the laser. Returning to equation Figure 4: Polar plot of radiation emitted by an electron oscillating in a linearly polarized plane electromagnetic wave, with low electric field strength (having a normalized vector potential of a 0 = 0.01) with the field vector along the x axis, calculated numerically. The vertical axis is the frequency normalized to the laser frequency, and the horizontal plane gives the radiated intensity. The angular shape of the structure is a set of sin 2 Θ lobes, in agreement with the analytic expression. 1 for a single particle, and inserting a time varying velocity β = a 0 sin(ω 0 t), and using a coordinate system so that angle Θ is measured between the polarization axis and the observation direction, the spectral power is: To first-order in a 0 , assuming a 0 1, this can be integrated to give: Here, the Dirac delta function δ(x) is a representation of (k/π)sinc 2 (x/k) in the limit that k → ∞. The radiated energy is therefore distributed at a single frequency with an angular structure consisting of a sin 2 Θ shape. In figure 4, the angular distribution of spectral power is shown from calculations using the quadratic algorithm with a time-step of ω 0 ∆τ = π/25, with a plane electromagnetic wave polarized along the x-axis with a field strength of a 0 = 0.01. Time is integrated to include 20 wave-periods. The angular shape of the structure is a set of sin 2 Θ lobes, in agreement with the analytic expression. Although the spectrum as a function of energy is a peaked distribution centered at ω = ω 0 , it is not a Dirac delta function. However, the analytic solution is for an infinite summation of wave-periods, whereas for obvious reasons a finite summation is calculated numerically. In terms of frequency, the calculation yields a sinc 2 (ω − ω 0 ) shape, characteristic of a finite window, which approaches the Dirac delta as the size of the window is increased.
As the field strength increases in the interaction, the electron motion becomes more complicated due to the magnetic field, and tends towards a figure of eight motion. This means that additional harmonics of the motion appear in the spectrum. For the case of exact back-scatter from an electron with an initial Lorentz factor of γ 0 , the motion of the electron means that it experiences a Doppler shifted plane wave, and therefore the period of oscillation is modified. In reference [24] an analytic expression for the shifted fundamental frequency, ω 1 , was given as: where β 0 is the initial velocity of the particle, and the power per unit solid angle in the mth harmonic of ω 1 is: where J ν (x) is a Bessel function of the first kind and ζ = a 2 0 /(2a 2 0 + 4). In figure 5, the power in harmonics of the fundamental frequency, normalized to the first harmonic, as a function of harmonic m, due to non-linear Thomson backscattering from a 5 MeV electron beam colliding with 10 periods of a linearly polarized plane wave with a normalized vector potential of a 0 = 1 is shown. The temporal resolution is such that ω 0 ∆τ = π/50, which represents only 100 sample points to calculate the spectrum from. Also shown is the (normalized) power from the analytic solution. Again the analytic solution models an infinite plane wave solution, and hence the harmonics are discrete. In the numerical solution, a finite number of periods is calculated, and therefore the harmonics have a sinc 2 (ω − ω m ) shape. The scaling of the amplitude of the peaks and the position of the harmonics agrees well between the two solutions. For the numerical solution, the spectrum was calculated as a function of frequency and then divided by the analytically calculated frequency ω 1 to give the horizontal axis. The 12th harmonic corresponds to a frequency of 784ω 0 , which for an 800 nm laser interaction corresponds to a photon energy of 1.2 keV. Figure 6 shows the angularly resolved plot of the same calculations as figure 5, and can be compared with figure 2(b) in reference [25], which shows a similar plot calculated from an analytic solution. Note that in their figure, only the first three harmonics were calculated. In the top right of the meshplot in figure 6, the next 3 harmonics can be observed. The finite width of the pulse, in both the numerical calculations here and the analytic calculations of reference [25], results in the sinc 2 (ω − ω m ) shape of the spectral peaks.
It should be noted that the radiation reaction force model of reference [26] has also been included to particle motions for other laser-particle interactions, and the energy loss by the particle compared with the energy in the calculated spectra, with good agreement. The question arises as to what resolution is required to accurately reproduce spectra using this method. Firstly, the particle trajectories themselves have to be reproduced accurately. Secondly, the quadratic interpolation must be an accurate representation of the function κ α x α . Since the calculation is an integral, slight discontinuities at grid points are not as important as the accuracy of the function between grid-points. This means that a spline interpolation method is not necessarily better, and has been found to be worse, than an interpolation that only relates to local grid-points. If the coefficient for a cubic interpolation term is small, χ 3n τ 3 1, then the spectrum should be accurately modeled. The result of this is that the algorithm is not particularly efficient for calculations of highly relativistic particles ( GeV) performing large radius orbits such as in a classic synchrotron, although it will yield correct results given sufficient resolution (the required resolution for such orbits scales as γ, because the trigonometric nature of the motion means a cubic term exists and κ α x α scales with γ). For such calculations other methods may be preferred. However, for laser or plasma based accelerator interactions, for example, where typically the electron energies are < 10 GeV but the oscillation frequencies are very fast, 10 13 − 10 16 Hz, the algorithm works well.