{'fraction_non_alphanumeric': 0.03781590112525364, 'fraction_numerical': 0.012227580572904314, 'mean_word_length': 4.521315291721228, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Radiation spectra of laser-driven quantum relativistic electrons

A procedure to calculate the radiation spectrum emitted by an arbitrarily prepared Dirac wave packet is developed. It is based on the Dirac charge current and classical electrodynamic theory. Apart from giving absolute intensity values, it is exact in terms of relativistic retardation effects and angular dependence. We employ a laser driven free electron to demonstrate the advantages of our method as compared to traditional ones that merely rely on the Fourier transform of the dipole operator's expectation value. Classical reference calculations confirm the results obtained for the low-frequency part of the spectrum, especially in terms of the observed red-shifts, which clearly deviate from non-relativistic calculations. In the high-frequency part of the spectrum, we note appreciable deviations to the purely classical calculations which may be linked to quantum averaging effects.


Introduction
The appearance of harmonics in the light scattered at free electrons is known both in theory [1][2][3][4] as well as in experiments [5][6][7] for a rather long time now. Atomic systems [8][9][10][11][12][13] in highly intense laser fields have turned out to produce especially high orders, the so-called High Harmonics [14][15][16]. In this case, the experimentally observed spectra [17][18][19] are usually explained in terms of the so-called recollision model [20]. Therefore, in the context of quantum mechanical simulations such as those shown in [21][22][23], apart from the obvious observables such as the expectation value of the position or momentum operator and the spatial probability density distribution, there is a lot of interest in calculating the radiation spectrum as well.
The numerical methods that were described in our previous work [22] provide the Dirac wave function as a function of time and position, where both variables are discretized and limited to certain regions. There exists a wellestablished procedure [24] that derives the emission spectrum from such a given data set (e.g. [12,25]). However, this one was developed for the nonrelativistic Schrödinger theory, and contains further approximations. The question arises whether, in the highly relativistic regime, where non-relativistic quantum theory breaks down, it is still advisable to employ a partly nonrelativistic method for the calculation of the quantum emission spectrum.
In this paper, we present a procedure that takes into account all relativistic spectral effects. It is based on the Dirac charge current and the classical theory of electromagnetism. Apart from giving absolute intensity values, it is exact in terms of relativistic retardation effects and angular dependence. To take advantage of the new scheme, one needs to obtain the time-resolved Dirac wave function of some physically interesting scenario in the first place. We made use of our advanced two-dimensional Dirac split-operator code [22] and considered the motion of a free electron in a laser field [26] to put the new scheme to the test. Classical calculations, which should not deviate much from the quantum calculations for this particular case as far as the lower part of the spectrum is concerned, serve as a reference and prove the correctness of the new scheme. However, significant deviations from the established method [24] are revealed. In addition to that, the higher order harmonics are averaged out. The delocalized quantum charge distribution is likely to be responsible for this effect. This remainder of this paper is organized as follows: Section 2 briefly reviews the formulas required to evaluate the radiation spectrum of a classical point charge and their non-relativistic limit. Section 3 briefly discusses the established standard methods, then turns to the derivation of our new method. In the end, the non-relativistic limit of the latter is evaluated to show that it matches the former. This section also discusses the solution of some numerical problems that are involved. Section 4 discusses the results obtained with our method and relates them to those obtained with the standard methods and classical formulas. We draw our conclusions in section 5 and provide the classical equations of motion in appendix A and the Dirac equation in appendix B, as well as some remarks and references concerning the way that we solved both problems numerically.

Classical radiation spectrum
The physical quantity of interest here is the differential amount of energy (d 2 W ) emitted into a certain interval of frequency (dω) and solid angle (dΩ), i.e. the quantity d 2 W dωdΩ .

Spectrum of a charged classical point particle in relativistic motion
From [27] we obtain a formula for the far-field radiation of a classical charged point particle which is ready for application: Here, y(t) represents the particle's trajectory and x the observer's position. The particle's charge is labelled q, and c is the speed of light. The integral dt(. . .) is replaced by a sum ∆t(. . .) which is carried out for all values ω that are of interest by using the discrete trajectory data for y,˙ y,¨ y, that is available from the numerical integration of the classical equation of motion (A.5).

Non-relativistic limit
If one omits relativistic corrections (| β| ≪ 1) and accounts (with respect to the retardation R) only for particle displacements in the direction towards the observer that are very small (x · y ≈ 0) as compared to the distance of the fixed observer | x|, then R is constant, ≈ const. (2) and instead of the integral in equation (1), the result is just the Fourier transform of the acceleration component a ⊥ that is perpendicular to the momentary direction of observation n: Using˙ β ⊥

Quantum mechanical treatment
The usual way (e.g. [12,25]) to calculate the emission spectrum is the following: Starting from the classical non-relativistic formula (4) the spectrum for a quantum mechanical system is evaluated using the corresponding quantity that is based on the acceleration expectation value a : For a more precise derivation and in addition to that a detailed discussion of the role of the "coherent" part (see above) and the "incoherent" part (not discussed here) of the spectrum, see [28,29]. The next step [24] is the application of the Ehrenfest theorem, i.e. the replacement which is easy to evaluate if the magnetic field B is neglected (the so-called dipole approximation). However, it is the magnetic field effects and the forward drift caused by them which have to be taken care of in this work with an emphasis on relativistic velocities. This precludes equation (4) from being the basis of any highly-relativistic theory. In a weakly-relativistic work such as [30], the starting point remains equation (4), however, the relativistically correct equation (A.3) is employed instead of equation (6). Utilising several essential approximations, the authors obtain an expression that can be evaluated numerically. Those approximations are no longer valid in the more strongly relativistic cases considered here, where even for classical calculations, the substitution of equation (1) by (4) can already lead to noticeable frequency shifts in the spectrum.
Unfortunately, in a quantum mechanical context, equation (4) is not easily replaced by the relativistically correct variant (1) and a calculation of the required expectation values in the framework of the Dirac theory. On the one hand there are fundamental reasons, such as the difficulty to define a reasonable velocity operator in the Dirac theory 1 , on the other hand one needs to take into account certain practical considerations 2 . Therefore we have dropped all corresponding considerations in favour of the following far more promising idea.
We start by calculating the emission spectrum of a spatially distributed classical charge current, instead of the one of a point charge. The following derivation differs slightly from those found in the literature [27] in that the one and only approximation made is the far-field expression for the electric field of a charge current. Otherwise it is exact. The transition to quantum mechanics is carried out at the very end of the calculation by substituting a suitable expression for the charge current.

Fully relativistic spectrum of a spatially distributed charge current
If we call the magnitude of the Poynting-vector (energy per area and unit time) | S|, its direction n, and the electric and magnetic field E and B, then we obtain We assume the radiation detector to be located at position r = rr and to be activated during the time interval [t 0 , t 1 ]. Let W(t) be a function 3 that describes the detection interval: If we designate energy by W , power by P and angular frequency by ω, then the energy per unit solid angle Ω collected by the detector at position r in the Here, we insert the Fourier transform of the electric field and its inverse: Since and therefore one finally arrives at The integrand of (13) yields The electric field itself is given by the time derivative of a certain vector potential A f as obtained in the usual far field approximation for a classical charge current  (t ′ , r ′ ) that is essentially localized at the origin of the coordinate system: Here,r is a constant unit vector that points towards the observer who is situated at a fixed distance r, i.e. at position r = rr. If we insert (16) into (15), convert the result to the frequency domain using (11) and, while keeping in mind (8), insert everything into (14), we finally obtain, neglecting contributions that vanish as O 1 r for large r In the above equation and some more to follow the symbol "⊗" is just a multiplicative line continuation marker, which we have to distinguish from the cross product "×". Note that so far, the calculation is an entirely classical one.
We are aiming at replacing the classical current by the quantum mechanical one. Since our numerical Dirac code delivers ψ and hence  at fixed, equally spaced times, whereas equation (17) requires very different times depending on position r ′ , we have to rewrite the above as follows: We substitute which, if we restrict ourselves to the special case h = 1, means This is still difficult to evaluate because of the variable bounds of integration. We expand them as needed and compensate the error in the integrand: Here, the time interval t 0 . . . t 1 is the detection interval, and t ′ 0 . . . t ′ 1 are the times for which we have to supply the current data. Thus, the latter mark the start and end of our numerical simulation. In the exponent, ωr c is large and constant, and can therefore be omitted. In fact, it should be omitted for numerical reasons: The functions sine and cosine do not work well for large arguments. The inversion of (22) and (23) tells us the detection time limits: Let V r 0 = supp(  (r 0 , r )) be the support of the current at the time r 0 and V = ∀r 0 V r 0 the sum of all these. Then, in (21), it is possible to replace the spatial integral over R 3 by another one over just V , and to calculate the times t 0 and t 1 in (25) and (26) as well (use ∀ r ′ ∈ V instead of R 3 ) -at least in principle. As we do not want to run our code twice -once to exactly evaluate V , the second time to actually calculate the spectrum (to save all currents from the first run would definitely exceed any storage capacity) -we simply choose V based upon classical estimates of the particle trajectory as a spherical volume (which facilitates arbitrary observation directions) of radius R, which is large enough to enclose all V r 0 . This is shown in section 1. If the current is not localized at the origin, as in the above derivations, but at the middle r K of the spherical volume V , then all positions vectors r, r ′ have to be moved accordingly: If r− r K | r− r K | points from the mid point r K towards the observer r, then one obtains from (25) and (26), taking into account (27) Changes according to (27) have to be applied to all equations (21) to (24), of course. The following explanation is meant to improve the readers intuitive understanding of the above purely mathematical description: If we designate the edge of volume V that lies furthest away from a given, fixed observer by the label a) and the other, closest edge by b), as shown in section 1, then obviously the time that the radiation needs to travel from a) to the observer is longer than from b). The simulation starts at r ′ 0 = ct ′ 0 . One would be inclined to activate the detector at t 0 = t ′ 0 + [travel time b) to the detector] and to close it at t 1 = t ′ 1 + [travel time from a) to detector] in order to avoid letting any radiation go by undetected. However, for consistency, then one would also have to take into account any radiation that reaches the detector at time t 0 as defined before, but which originates from position a) at a time before the beginning of our actual calculation, and about which therefore, no information is available. Thus, our conservative approach is to open the detector a little bit later. The choice t 0 = t ′ 0 +[travel time from a) to the detector] throws away some radiation of origin b), but now everything that is detected is properly taken care of. Similarly, at time t 1 as defined above there would be a lack of any information about radiation that is coming from b) and superimposing the radiation that is still arriving from a). Therefore, the conservative choice t 1 = t ′ 1 + [travel time from b) to the detector] has to be made here, too. Considering equation (28), one easily notices that the detection interval can become rather short for large values of R. Thus, a spatially extended motion can only be analyzed if it is slow, or more exactly speaking, if t ′ 0 and t ′ 1 are widely separated in time. At least for a particle initially at rest, t ′ 0 can be moved into the past as much as required without altering R, whenever this is not the case.

Application to quantum mechanics
In relativistic quantum mechanics, the probability current density is given in the following way [32] and can be easily decomposed into the vectorial current density and the scalar density For the charge current  in the context of Dirac theory, we choose This choice is less obvious than it seems. For a extensive discussion of this matter, see [33]. The time derivative ∂ ∂r 0  = 1 c ∂ ∂t  is derived in the following way: Note that (34) is an immediate result from the Dirac equation (B.1). The momentum operator (35) is transformed by means of Fourier transforms into a simple multiplication operator.

Reduction to two dimensions
We have to reduce (21) by means of (B.3) to the two-dimensional case: If the observer is located inside the plane, i.e.r · r ′ =r · r ′ ⊥ then equations (21), (22), (23) and (24) are transformed into To evaluate the above, we need and as a part of the above Finally, insert equation (41) (apart from the factor δ(0), which has to be omitted, cf. appendix B) into (40) and the result into (36). Together with (37), (38) and (39) this would give (for the case h = 1) the final formula used in our code. Because of its size, we do not print it again.

Numerical limitations
Tests showed that the derivative of the current as shown before in equations (40) and (41), although correct in theory, cannot be numerically calculated correctly this way. Apparently, the reason is to be found in equation (34): This equation is only true, if ψ is a exact solution of the Dirac equation. Small errors in, for example, the phase of the wave function can have a large influence on the result. However, the current (32) itself is independent of the phase. Therefore, a less elegant, but numerically superior method is a finite difference scheme like the following: A higher order scheme would be desirable, however this would require a rather impracticable multiple storage of the wave function. Fortunately, ∆t has to be chosen only small enough as to let ω max = π ∆t cover the whole spectral range that is to be expected. It can therefore be much larger than the step size of the numerical propagation of ψ.

Non-relativistic limit
Just to gain some trust into the above derivation, we now look at its nonrelativistic and spin-less limit. It can be obtained by ignoring small retardation effectsr · r ′ ≈ 0 and using the Schrödinger charge current 4 (which depends on the Schrödinger wave function ψ t ∈ C) Similar to the classical theory's limit in section 2.2, we insertr · r ′ ≈ 0 and (43) into (20) and, after some work, arrive at: Here r ′ 0 = ct ′ and p designates the momentum operator whose expectation value is defined, for the case of perfect conservation of normalization (ψ t , ψ t ) = 1 ∀t, as follows: At this point, we can employ the Ehrenfest theorem and, for the case of a vanishing magnetic field B ≈ 0, we write to arrive at a simple formula: Calling a Def.
= q m E and a ⊥(t ′ ) Def.
=r × a ψ ′ t ×r we obtain If we expand the integration interval to infinity, as indicated above, and replace a ⊥ ↔ a ⊥ , then we recover the classical result (4) in (50). It returns the spectrum as proportional to the Fourier transform of the perpendicular component of the mean acceleration. 5 If however, we neglect the retardation r in the integration limits of equation (50), then we obtain the result from [24]. 6 4 Spectra of free electrons in a laser field In order to appreciate the work carried out in section 3, we have made a comparative study of five different methods. In this study, the physical parameters are primarily chosen in order to magnify the differences in the spectra and to speed up calculations.
The scenario is the following one: A free electronic wave packet is propagated over 20 cycles, two of which are dedicated to the field amplitude's sin 2 -shaped turn-on and turn-off phase, respectively. In order to see a reasonable number of harmonics in the spectrum of an otherwise free electron [1], such a large number of cycles is necessary. A free electron is favourable to the numerics: There are, at least in polarization direction, no large canonical momentum components, and the wave packet is not torn apart by scattering in addition to its unavoidable natural spreading. Therefore, it is possible to use small grids of rather low spatial resolution. The large number of cycles implies a high frequency in order to keep the overall propagation time short. Note that the time resolution for Dirac calculations is almost independent of the problem 7 and only given by ∆t E ≈ mc 2 . Therefore the frequency was chosen as ω = 7.120 a.u. (corresponding to 6.4 nm). 8 In order to stay in the strongly relativistic regime -it is only here that we can expect appreciable differences -the intensity has to be set to the challenging value E 0 = 300 a.u. (I = 3.15 × 10 21 W/cm 2 ). The following calculations were carried out: 5 Note that most authors ignore the correct prefactor and give the spectrum in "arbitrary units" only. 6 The authors of [24] employ the length gauge E(t) = − ∇A 0 with A 0 = − E(t) · x and the dipole approximation A = 0. 7 Extremely high kinetic energies in the MeV regime are an exception to this rule, see [23]. They require an even higher temporal resolution since E = γmc 2 and γ > 1. 8 This particular frequency is the highest one that is indicated for the future FEL of the so-called ,,TESLA-Test-Facility (Phase 2)" at DESY in Hamburg [34].
• the classically-relativistic trajectory of a point particle according to the procedure described in section A. • the quantum mechanical propagation of a Gaussian wave packet according to the Dirac equation as shown in [22]. • the quantum mechanical propagation of a Gaussian wave packet according to the Schrödinger equation. 9 From these data sets, the following spectra were produced: • Based on the classical trajectory data, the following was calculated: · the relativistically correct spectrum according to equation (1). This will be labelled "classical theory (exact)" in the figures to follow. · the non-relativistically approximated spectrum according to equation (4), labelled "classical theory (approx.)". In both cases we used the same (relativistically correct) trajectory data set, i.e. the distinction relativistic vs. non-relativistic only refers to the way that the spectrum is evaluated.
• From the single data set according to the Dirac theory, two different kinds of spectra were extracted: · the fully relativistic one ("Dirac") according to the method presented in subsections 3.1 to 3.4: It contains no non-relativistic approximations and takes into account all retardation and spin effects. · a half-relativistic one ("Dirac-Schrödinger"): the expectation value of the electric field was evaluated from the density distribution, and then, according to the procedure in section 3.5, the spectrum was calculated. In other words, the relativistically correct Dirac probability density was fed into a non-relativistic procedure, which was originally developed for the Schrödinger probability density. Spin effects are taken care of for the propagation, but not for the spectrum calculation, i.e. for the charge current. • From Schrödinger data set ("Schrödinger"), the spectrum was evaluated according to the method in section 3.5. This represents an entirely nonrelativistic version, with no spin effects at all.
Additionally, the observation angle ϑ is varied from 0 • to 90 • in steps of 30 • . The "Dirac-Schrödinger" method allows us to differentiate between effects caused by different wave function propagation methods and those that originate from the different treatment of the emitted radiation.

Comparison of the results
We point out the following features of the spectrum that was calculated using the Dirac charge current in section 2 and the magnification thereof in section 3: For ϑ = 0 • , there are only harmonics of first, third and fifth order. Third and fifth order are very weak and hard to distinguish from the oscillating underground. No red shift is found. At ϑ = 30 • , there are more harmonic orders: Up to the seventh order both even and odd harmonics can be seen, and a small red shift is found. This effect is enlarged for ϑ = 60 • . For ϑ = 90 • we find the largest red shift and harmonics of up to ninth order. The intensity of the fundamental line decreases to a level even below the one of the second to fourth order harmonics. The underground intensity also decreases. The underlying reasons are best understood by looking at the similar, but more intuitive classical picture. Looking at the spectrum in section 4 and its magnified version section 5, which are based upon the classical trajectory and the corresponding relativistically correct radiation formula, we notice that for ϑ = 0 • there are only odd harmonics, but both even and odd ones up to about 19th order for all other angles. For order above 20, the peaks disappear in the noise. There is no red shift at ϑ = 0 • , but an increasing one for larger angles. At ϑ = 90 • , the intensity of the fundamental line drops below the one of the second to fifth order harmonic. The similarities between this classical calculation and the Dirac based one are obvious. The main difference is that due to the smeared out charge distribution in the quantum mechanical Dirac calculation, some sharp structures of the classical spectrum are averaged out, especially its higher frequency part.
The growth of the red shift with increasing angle can be understood qualitatively quite easily (cf. [9]): In the rest frame of the electron, the laser frequency is red shifted because of the forward drift of the electron in the laser field, and the same is true for the oscillation frequency of the electron in polarization direction due to the laser's electric field component. If the electron drift is directed exactly towards the observer (ϑ = 0 • ), the latter sees this actually red-shifted frequency in the electron's radiation, which however is Dopplershifted back to exactly its original position in the frequency domain. If the observation angle is larger (ϑ > 0 • ), the original red shift stays the same, but the compensating Doppler blue-shift decreases with increasing angle, leaving behind an increasing net red shift.
The appearance of even and odd harmonics can be understood by the means of a symmetry consideration. One needs to consider the approximative formula (4), i.e. the magnitude squared of the Fourier transform of the acceler-ation component a ⊥ perpendicular to the direction of observation. Reference [4] provides the classical analytical solution for the particle velocity β(η) in terms of the laser phase η. If we insert a simple harmonic cw-laser field in those formulas and take the time derivative, we obtain for the acceleration something like where u(η) is a rather complicated function, with odd symmetry with respect to both η = 0 und η = π 2 , whereas g(η) features even symmetry with respect to η = 0 and odd symmetry with respect to η = π 2 . In the above, σ 1 and σ 2 are unit vectors in propagation and polarization direction, respectively. The properties of the Fourier transform then automatically lead to even harmonics for the σ 1 -component of acceleration and odd ones for its σ 2 -component. The latter ones are exclusively visible for the observation direction designated by ϑ = 0 • , and only the former ones for ϑ = 90 • . At any angle in between, one would expect the appropriate mixture of both even and odd harmonic orders. The above analytical discussion explains the complete disappearance of even orders in the numerical calculations for ϑ = 0 • just as well as the suppression of the fundamental (which is odd by itself) at ϑ = 90 • . All deviations from the mathematically ideal behaviour have to be attributed to the relativistic corrections, the retardation r −r · r ′ = 0 and last but not least the usage of a time-limited, pulse-shaped laser pulse.
Summing up the above, we can say that the classical and the quantum mechanical description are in good agreement and both deliver well understood spectra. One might therefore ask immediately for the motivation of the computationally expensive Dirac calculations. In fact, this particular scenario was set up mainly because we wanted to verify the correctness of the procedure described in subsections 3.1 to 3.4 using a system from which one would not expect large quantum effects, so any deviations from the well-known classical behaviour would easily reveals errors in the method or the code. The agreement between classical and quantum theory is likely to disappear however for laser-driven bound systems. The parameter regime usually referred to as "tunnelling regime", for which a description of the "High Harmonic Generation" process according to the recollision model [20] is applicable, would be of high interest. Unfortunately, this regime is characterized by low frequencies and therefore long computational runtimes, which renders it impossible to study with our Dirac code, at least when simultaneously employing field strengths that are large enough to render relativistic and spin-related effects visible. Low intensities would be feasible, but not very interesting, as there are already a number of completely sufficient Schrödinger [12] and Pauli [25,30] calculations for this regime.
However, we note some interesting discrepancies when we compare the method from subsections 3.1 to 3.4 and the ordinary procedure as shown in section 3.5. The latter, or more exactly equation (48), can, although originally deduced from Schrödinger theory, also be evaluated using the Dirac probability density: This is achieved by Fourier transforming the perpendicular acceleration component a ⊥ as in equation (50), where the acceleration is taken as the expectation value of q m E( x, t) with respect to the Dirac wave function. As an example, this is shown in section 6 for an observation angle ϑ = 30 • and enlarged in section 7 (second graph from top to bottom in both figures) side by side with the exact Dirac procedure, the also exact classical method according to equation (1) and its non-relativistic approximation (4), as well as the non-relativistic Schrödinger calculation 10 , which is also based upon equation (50), but is derived from a non-relativistic Schrödinger wave function. As one can see, the three not entirely relativistic methods ("Dirac-Schrödinger", "classical theory (approx.)", and "Schrödinger") significantly over-estimate the red shift 11 , because in these calculations, the Doppler effect due to the electron's drift in propagation direction is not included. The even harmonics, which should be visible for any non-zero observation angle, are missing in both spectra according to equation (50) ("Dirac-Schrödinger" and "Schrödinger"). This is because only the perpendicular component of the electric field is taken into account, weighted either with the Dirac or the Schrödinger probability density. As a result, the spectrum features a simple cos(ϑ)-dependence and even vanishes for ϑ → 90 • . There is no change in structure with varying angle, it has always the shape of the version for ϑ = 0 • , although it is only valid for this particular angle.
Apart from the methodical differences laid out above, we also see differences between the two classical calculations and all three quantum mechanical ones. Especially the higher order harmonic peaks that are visible in the classical curves are averaged out in the case of the quantum spectra.
10 Note: The Schrödinger calculation only seems to exhibit far less noise than all the other calculations. The reason is that it was carried out using a much coarser frequency resolution (about one fifth). Apart from that, its prefactor is not known. To indicate this, the vertical scale in the corresponding plots was omitted. 11 Note: These three methods even predict a red shift at ϑ = 0 • , which is not shown here.

Conclusions
From the above, it can be concluded that the procedure laid out in detail in subsections 3.1 to 3.4 is superior to the established techniques in terms of angle dependence of especially the even harmonic orders, and in terms of the predicted red shifts, whenever the parameters are chosen in a way that leads to a significant influence of retardation and magnetic field effects. In principal, even spin effects 12 could be analyzed this way, since the method is based on the Dirac theory. Although originally developed for the usage in combination with Dirac wave functions, the method could be employed in the context of the Schrödinger or Pauli theory as well 13 , if one abandoned the Ehrenfest theorem and evaluated equation (44) directly, and it would at least correct the angle dependence. Employing equation (20) would even lead to correct red shifts, and the only remaining quantitative deviations would be due to the underlying non-relativistic wave packet propagation. The computational cost for the spectrum alone is not essentially higher than for the more simple schemes, therefore its general usage can be recommended. However, the parameter regime in which its advantages are best visible are already very demanding as far as the wave packet propagation is concerned.
For the free electron scenario, classical results are very similar to the ones obtained with our new method. We see no true quantum effects apart from the disappearance of the higher order harmonics. This may be attributed to averaging effects by the spatially extended Dirac charge current. More significant effects are very likely to be observed once that atomic systems are considered. It must be noted, though, that our main goal was to devise a relativistically correct scheme, as opposed to the non-relativistic existing ones. Our results show that the relativistic corrections are significant for the parameters employed here.

A Classical equations of motion
In classical relativistic mechanics, momentum change equals the applied force, where the force f is the well-known Lorentz-force: Equation (A.1), which represents a second-order differential equation with respect to the three-dimensional position variable x can be easily transformed [10] into a first-order system of differential equations in terms of the sixdimensional combined variable ( x, v): Formally, the above means where G ( x, v,t) represents the six-dimensional function that is obtained from the right hand side of equations (A.3) and (A.4) after inserting the force (A.2). The first-order system (A.5) is then easily numerically integrated by means of a standard fourth-order Runge-Kutta algorithm [36].

B The Dirac equation
The Dirac equation in 3 + 1 dimensions reads as follows: In the above, is Planck's constant, q and m represent the particle's charge and rest mass, c is the speed of light, α i and β are the usual Dirac matrices (i = 1, 2, 3) [32]. We use atomic units, so for an electron = m = −q = 1 a.u., and c = 137.036 a.u. Suppose we have a potential that is invariant under translations in one particular direction that is designated by n then we can use the separation ansatz 3) in the Dirac equation (B.1). The result is: (B.4) In general, there will not be a single momentum p 0 , but rather a momentum distribution giving specific weights for arbitrary p ∈ R and the solution is a wave packet instead of (B.3). We restrict ourselves to the simple choice p 0 = 0 to eliminate this constant. Other values would not provide more physical insight, and true wave packet calculations would cause the computing times to grow into the regime of true three-dimensional calculations. Note however, that the ansatz (B.3) inevitably leads to a meaningless factor δ(0) in bilinear expressions of ψ, such as expectation values or the integrated current density.
We finally obtain an effective Dirac equation in two dimensions: The above equation is integrated [37] using the split-operator scheme [38] presented in our earlier work [22]. The dynamical grid techniques described in there are especially well-suited for the free electron problem discussed here, and dramatically reduce computation times.
A vector potential that describes a laser beam of electric field amplitude E 0 and frequency ω, which is linearly polarized in directionâ and propagates in direction of the wave vector k, such as is invariant in direction n Def.
=k ×â. As can be easily seen, in this context "twodimensional" does not mean that there is no magnetic field B (which points in the x -direction). Also note that the above two-dimensional treatment is exact for the special case of a free electron in a laser field that we considered in this paper, because the vector potential (B.6) features the required symmetry. Fig. 1. Sketch of the radiation scenario: The three black-shaded areas that carry time (t ′ ) labels represent the electron at the beginning (t ′ = t ′ 0 ) and the end (t ′ = t ′ 1 ) of the simulation, as well as at an arbitrary point in time (t ′ ∈ [t ′ 0 , t ′ 1 ]) in between, each enclosed inside the dashed boundaries of the co-moving and -growing numerical grid. The oscillatory curve is the trajectory of the particle's centre. The observer (or detector), which is activated during the time interval t 0 . . . t 1 , is located under an angle ϑ with respect to the laser propagation direction at a distance | r − r k | as measured from the middle of the volume V (large grey-shaded circle), which is a sphere of radius R and origin r K . At the same time, the radius R represents half the diagonal line of the solid rectangle, which is the enclosing envelope of all occurring numerical grids. The two parallel lines perpendicular to the line of observation labelled "a)" und "b)" mark two specific distances to the observer mentioned in the text.      The plot shows the spectrum of a free electron in interaction with a 20 cycle laser pulse (including two cycles for turn-on and turn-off, respectively) with an amplitude of E 0 = 300 a.u. (I = 3.15 × 10 21 W/cm 2 ) and a frequency of ω = 7.12 a.u. (λ = 6.4 nm). The fully relativistically calculated spectrum due to two-dimensional motion of a classical point particle is shown for several different angles of observation ϑ with respect to the propagation direction.    Fig. 6. Comparison of several methods to evaluate the radiation spectrum: The top graph is based on the method from subsections 3.1 to 3.4 and shows the fully relativistic spectrum due to a two-dimensional Dirac charge current. The second graph from top to bottom displays the spectrum obtained by Fourier transforming the perpendicular acceleration component a ⊥ , where the acceleration is taken as the expectation value of q m E( x, t) with respect to the Dirac wave function, similar to the way shown in section 3.5 for a Schrödinger wave function. The third graph shows the fully relativistic spectrum calculated for the classical trajectory according to section 2.1, and the fourth graph the one that can be obtained from the same trajectory data in the non-relativistic limit according to section 2.2. The bottom graph displays the result according to the method in section 3.5 for a non-relativistic Schrödinger calculation. No vertical scale is given for this last calculation, which is a contribution of Andreas Staudt, because the prefactor is not known. The differences of the distinct methods in terms of red shift and harmonic orders are obvious. All methods that contain non-relativistic elements (Dirac-Schrödinger, approx. classical theory and Schrödinger) feature significant red shifts at the observation angle displayed here (30 • ). For the exact classical and the Dirac-calculation noticeable red shifts only occur for much larger angles. Apart from that, even harmonic orders are missing in both the Dirac-Schrödinger and the Schrödinger calculation. The latter only seems to have less numerical noise than all the others, but this is due to a coarser frequency resolution.