Quantum Electrodynamics vacuum polarization solver

The self-consistent modeling of vacuum polarization due to virtual electron-positron fluctuations is of relevance for many near term experiments associated with high intensity radiation sources and represents a milestone in describing scenarios of extreme energy density. We present a generalized finite-difference time-domain solver that can incorporate the modifications to Maxwell's equations due to vacuum polarization. Our multidimensional solver reproduced in one-dimensional configurations the results for which an analytic treatment is possible, yielding vacuum harmonic generation and birefringence. The solver has also been tested for two-dimensional scenarios where finite laser beam spot sizes must be taken into account. We employ this solver to explore different types of laser configurations that can be relevant for future planned experiments aiming to detect quantum vacuum dynamics at ultra-high electromagnetic field intensities.


Introduction
The prospects offered by ultra-intense laser sources in the infra-red (IR) or X-ray central wavelengths [1] have triggered a renewed interest in Quantum Electrodynamics (QED) and its impact on quantum processes at a macroscopic scale, namely how such phenomena can affect well studied interactions in the fields of plasma and laser dynamics. The most relevant QED processes in strong fields and high intensity laser interactions have been explored in several reviews [2,3,4]. Among these effects, the second order QED process of photon-photon scattering mediated by the vacuum fluctuation of virtual electron-positron pairs has been a topic of renewed interest motivated by several exotic consequences [3,5,6,7] that originate directly from the original Heisenberg & Euler Lagrangian [8]. However, many of these effects, such as the virtual polarization of the vacuum, remain to be experimentally observed with the use of ultra-high intensity laser physics. With expected peak intensities up to 10 23 − 10 24 Wcm −2 to be delivered by large-scale facilities such as the Extreme Light Infrastructure (ELI) [9], the VULCAN 10 PW project [10], or the HERCULES laser upgrade [11], the regime where these virtual fluctuations can be detected in the laboratory is close to being within reach. In particular, experiments are now being planned to study the quantum dynamics of the vacuum [12,13] by combining ultra-intense optical lasers with X-ray lasers [14]. The increasing consensus regarding the importance of quantum dynamics in the collective effects of many extreme laser plasma interactions has motivated the development of novel numerical tools that couple the multiple scales associated with the problem. Numerical codes that simulate quantum radiation reaction [15,16,17] and pair production effects [18,19,20,21,22], have already made important predictions in extreme energy density scenarios [22,23,24,25]. We propose here an algorithm, different from the one suggested by Domenech & Ruhl [26], that includes the effect of vacuum polarization in the spatio-temporal evolution of the electromagnetic field in multi-dimensions and allows for a broad set of initial conditions. Without dwelling deeply in the details of each approaches, both have advantages that deserve to be highlighted. Our approach is computationally less expensive and can be implemented in the PIC loop. On the other hand, the algorithm proposed by Domenech & Ruhl [26] is based on a gridless method which removes any interpolation issues and the precision of the numerical scheme is of fourth order. These vacuum quantum effects can be integrated via an effective nonlinear permeability and permittivity, and therefore to use a semi-classical approach. The effects of the quantum vacuum can be important and appreciable, not only in scenarios involving high intensity electromagnetic radiation, but also in extreme astrophysical environments surrounded by near critical Schwinger magnetic fields (e.g. neutron stars) where the propagation of electromagnetic waves is modified [27,28].
The electron-positron pair vacuum fluctuations were first taken into account by Heisenberg and Euler (HE) who calculated the first full Lagrangian to all orders [8].
In the low field E E s , low frequency ω ω c limit of the electromagnetic (EM) fields, the leading corrections of the standard Maxwell (M) Lagrangian density [8] can be written as where ω c = m e c 2 / is the Compton frequency, E sch = m 2 e c 3 /e the Schwinger critical field, and the EM invariants F = (E 2 − B 2 ) /2, G = E · B, E and B the electric field and magnetic field respectively. The nonlinearity coupling parameter is with α = e 2 / c. The parameter weights the relative importance of the quantum corrections compared to the classical fields and vanishes in the limit → 0. From the Euler-Lagrange equations for the electromagnetic fields, we obtain a set of modified QED Maxwell's equations [6] ∇ · D = 0 (3a) with and The nonlinear vacuum polarization, P , and magnetization, M read We would like to stress that we are using CGS units throughout the article. This semi-classical formulation treats the vacuum as a nonlinear medium, when the EM invariants F and G do not vanish. The algorithm aims to solve the nonlinear set of corrected Maxwell's equations leveraging the smallness of the vacuum non linearity. This algorithm described in this article is second order accurate in time and space. Due to its design, a key feature of the algorithm is that it can be seemingly incorporated into massively parallel fully relativistic electromagnetic particle-in-cell codes such as OSIRIS [29]. This could also allow for studying self-consistently scenarios where charged particles are also present in the system, and to be explored in future publications. This paper is organized as follows. In section 2, we describe the numerical algorithm, a generalization of the Yee algorithm, that solves Eqs.(3a-6) in multi-dimensions. Section 3 is devoted to the induced vacuum birefringence on an electromagnetic probe beam while crossing a region of high field. Several configurations of high field regions are considered such as a static field or realistic optical laser pumps. Harmonic generation constitutes the focus of section 4 in one-dimensional and two-dimensional geometry. Finally in section 5 we state the conclusions.

Description
The standard second-order finite-difference time domain (FDTD) method to solve Maxwell's equations is the Yee Algorithm [30]. The Yee scheme solves simultaneously for both electric and magnetic fields by solving Faraday's and Ampère's law, respectively. The explicit linear dependence of Maxwell's equations on the fields allows the field solver to be centered both in space and time (leapfrog scheme), thus providing a robust, second order accurate scheme without the need to solve for simultaneous equations or perform matrix inversions [31]. Moreover, the efficiency and simplicity of the Yee scheme allow an easy incorporation into numerically parallel codes.
To solve the QED Maxwell's equations, a modified Yee scheme was developed to address the two main difficulties which arise from the nonlinear terms. Firstly, all fields must be evaluated at all grid positions as opposed to spatially staggered fields. This permits to accurately evaluate quantities such as the EM invariants F and G, since a given component of the nonlinear polarization and magnetization vectors now fully couples all other field components as can be understood from Eqs. (6)(7). This is a significant obstacle regarding the essence of the Yee scheme since the algorithm may no longer be correctly spatially centered. Secondly, the temporal derivative of the nonlinear polarization term in Ampère's equations prevents each electric field component to be advanced straightforwardly as it requires the knowledge of future quantities. This is easily understood through the discretization of the modified Ampère's law in one dimension where the indices n and i denote the temporal and spatial positions, respectively. Usually, one would isolate the electric field term of temporal index n + 1 to advance this field in time. However, to calculate this component one must know the polarization at time step n + 1, which is a nonlinear function of all other fields at the new time step. The latter difficulty served as the main motivation to develop the modified Yee scheme proposed here. The numerical scheme is illustrated in figure 1 for a time step ∆t: -we begin by advancing the fields using the standard Yee scheme (i.e. without accounting for the polarization and magnetization of the vacuum). This setup allows us to obtain predicted quantities for the values of the fields at the new time. This approach is based on the standard technique of the predictor-corrector method, where the linear Maxwell's equations are solved as the zeroth order solution to the fields; -the predicted field values are then interpolated at all spatial grid points using a spline interpolation method thus allowing to calculate quantities such as the EM invariants and respective polarization and magnetization of the vacuum, to lowest order; -the polarization and magnetization are then used to advance the electric field via the modified Ampère's law; -the convergence loop re-injects this new electric field value back into the polarization and magnetization source terms Eqs. (6)(7) to refine these quantities and re-calculate the electric field iteratively. This loop is reiterated until the electric field converges to a value within the desired accuracy; It must be emphasized that this method is only valid as long as the effects of the polarization and magnetization of the medium are small compared to the non-perturbed propagation of the fields given as solutions to Maxwell's equations in classical vacuum. This condition is automatically satisfied for realistic values of electromagnetic fields available in current, or near future, technology. In this regime, the QED theory is valid since the Schwinger field, around which spontaneous pair creation (Schwinger effect) becomes non-negligible, corresponds to an electric field of E s ∼ 10 18 V/m, whereas ambitious laser facilities aim to push available intensities to the 10 23 − 10 24 W/cm 2 (E ∼ 10 −3 E sch ) range. The order of the ξ parameter in Eqs. (6)(7) which is on the order 10 −6 E −2 sch clearly helps to ensure the validity of the method. Therefore, this scheme highly benefits from the fact that the nonlinear QED corrections of the vacuum are perturbative in nature. The convergence loop can be seen as a Born-like series since for every re-insertion of the fields back into the nonlinear source term, there is a gain in accuracy of one order in the expansion parameter to the result. The algorithm proposed here solves Ampère's law by treating the nonlinear corrections as a source term, in an iterative manner, where S NL = ∇×M +∂ t P . From this discussion and equation (9), we can conclude that this generalization of the Yee scheme can be extended beyond the framework of QED corrections to the vacuum as it is valid to solve Maxwell's equations in any nonlinear medium provided that the polarization and magnetization are given and that their order is such that they can be treated as a perturbation. This possible generalization enhances the range of applicability of our algorithm. Furthermore, the inclusion of a current in the algorithm (J = 0 in Ampère's law) can be done, both within a PIC framework or for a macroscopic field dependent current by including the current term in the initial standard Yee scheme loop where the predictor quantities are computed. This is another key feature regarding the ability to couple our proposed generalized Yee solver to the PIC framework, of paramount importance to model scenarios where charged particles (even in small numbers) are present.

Interpolation of the fields
The algorithm requires that all fields are calculated at the same spatial positions. When considering the spatial interpolation of the self-consistent fields given by the Yee Algorithm we found a clear asymmetry between interpolating the electric field at the magnetic field position or vice-versa in terms of the precision of the EM invariant E 2 − B 2 for both cases. Since a plane wave is an exact solution of the QED Maxwell's equations, the invariants calculated in the simulation should be identically zero [32]. Figure 2 shows the distribution of the EM fields within a two-dimensional Yee grid cell. We found that all the standard interpolation schemes yield invariants with much greater precision at the lower left corner of the cell compared to the other positions. This difference in precision was found to be of two orders of magnitude when tested for a plane wave in 1D, which can affect the stability and precision of the code. The reason for this artefact is due to the way that the fields are initialized within the simulation domain. In particular, the fact that the electric and magnetic fields must be initialized with a shift both in space and time, creates an asymmetry between interpolating a field to the corresponding position of the other field, even if this interpolation is done in a centered manner. The solution we have adopted to address this problem is to calculate all the fields at the cell corner where the invariants are known to be of higher precision. For instance, the B z component at the left corner of the cell becomes where I is an interpolation function. Once all fields are calculated at the (i, j) positions, we can compute the invariants at these positions and then re-interpolate these invariants directly to the other grid cell points in a similar fashion. The correct calculation of the EM invariants is necessary in order to evaluate the nonlinear polarization and magnetization of the vacuum via equations (6)-(7).

Numerical Stability
The method adopted to study the numerical stability of the QED polarization solver follows the standard mode analysis [31]. With the linear Yee scheme, the onedimensional numerical dispersion relation for a plane wave propagating on a grid with spatial and temporal resolution ∆x and ∆t respectively is [31] ω = 1 ∆t arccos 1 + c∆t ∆x 2 (cos(k∆x) − 1) .
2D grid A notable case is when ∆t = ∆x/c for which equation (10) reduces to the EM dispersion relation for a plane wave in vacuum, ω = kc. To study the stability of the new set of QED-corrected Maxwell's equations using this method, a self-consistent numerical dispersion relation was derived. Due to the nonlinearity of the equations, the new numerical dispersion relation can be written as where E 0 is the amplitude of the wave and F NL is a nonlinear function of ω, k, and the spatial and temporal steps. In the classical limit ξ → 0 the RHS goes to zero and the dispersion relation reduces to equation (10). A numerical plane wave propagating via our QED solver will therefore obey equation (11). For a numerical plane wave the EM invariant E · B is identically zero, whereas the invariant E 2 − B 2 will not vanish identically due to finite spatial resolution and the fact that the fields must be interpolated in space to evaluate the invariants, as already discussed above. Therefore, the amplitude of this EM invariant depends on the interpolation method and grid resolution. We calculate this dependence by evaluating E 2 − B 2 at a given grid point, taking into account that a correct centering in space implies that one of the fields must be interpolated to the position of the other (in this case the B field, using linear interpolation). This yields This expression was compared to the results extracted from one-dimensional simulations. Figure 3 shows a comparison between equation (12) and simulations with several seeded k modes and with ξE 2 0 = 10 −4 , c∆t = 0.98∆x and k∆x = π/100. The simulation points agree with the trend presented by the theoretical curve. This result shows that equation (12) provides an upper bound to the interpolation error when seeding a particular k mode. In particular, the results show that for higher wave numbers, up to the resolution limit, the order of magnitude of the invariant amplitude increases, tending towards unity. One shall therefore limit the simulations to low k modes to ensure the smallness of the invariants.
The stability of the QED Yee solver, i.e., the nonlinear dispersion relation, equation (11) was solved using three methods: a numerical solution, an analytical solution through the linearization of the system via the ansatz ω = ω 0 + δω with δω ω 0 , and finally by estimating the growth rate of the maximum mode allowed by the grid resolution i.e. k∆x = π. The results are shown in figure 4 for simulations performed with a grid resolution of ∆x = 0.0314, ∆t = 0.98∆x and ξE 2 0 = 10 −4 . One can verify in figure 4 that the analytic solution is in excellent agreement with the numerical integration. The maximum growth rate is given by where ε = 1 − c∆t/∆x. The maximum growth rate predicted theoretically by equation (13) serves, therefore, as an accurate rule-of-thumb criterion to understand how unstable a given simulation setup may be. Finally we took the solution of the perturbative expansion and studied the limit for small k values, which yields lim k∆x→0 Im(δω) = 1 4 Equation (14) suggests that the smallest k modes will be the least unstable as not only does the growth rate scales with the small quantity ξE 2 0 , but also due to the power law applied to the small value of k∆x. This is an important result since, in principle, the low k modes, for a given grid resolution, are those that will be seeded for a simulation setup. These theoretical predictions were compared with one-dimensional simulations  (11) Maximum growth rate, equation (13) Numerical solution of equation (11) 0 20 40 60 80 100 by extracting the growth rate of a given k mode in the simulation domain, the results are shown in figure 5. The growth rate was extracted from the Fourier spectrum of the simulations for different values of ξE 2 0 . Figure 5 shows a good agreement between the maximum growth rates extracted and equation (13). Our theoretical analysis shows that the growth rate of the most unstable mode scales linearly with ξE 2 0 . Furthermore, we performed simulations under the same conditions, varying only the seeded k mode and verified that this does not affect the growth rate of the most unstable high k modes. Instead, it is the amplitude of the seeded mode that affects the growth rate of the higher k modes by nonlinear coupling. It is possible to derive a criterion for the time at which the seeded field starts to be strongly deteriorated by the growing numerical noise, by assuming this blow-up occurs once the amplitude of the fastest growing k modes δẼ (initially this amplitude is at the numerical noise level, and can be measured from the initial spectrum of the fields in the simulation), become of the order of the initial seed amplitude. This criterion yields, For realistic values of ξE 2 0 , this time is far greater than any simulation setup one may wish to perform.

Vacuum Birefringence
A thorough benchmark of the functionality and robustness of the algorithm may only be gained by comparing simulation results with analytical results in 1D simplified cases. One-dimensional scenarios provide excellent opportunities to test the code against analytical predictions. The two cases we present here are the vacuum birefringence in the presence of a strong static field and counter propagating plane waves. Whilst the first case is well studied in the literature [33,34], the second case requires a finer analytical work, yielding nevertheless the well known result of generation of higher harmonics due to the nonlinear interaction as shown in [35,36,37] for different setups and physical regimes that will be addressed in the next section. The true physical value of the parameter in normalized units of the simulation is ξ ∼ 10 −17 for an optical frequency. Simulations were performed with artificially increased values of the ξ parameter to better illustrate the method proposed here. This does not alter the physical relevance of the results. Rather, this is simply a re-scaling of a constant to highlight the effects more clearly. Finally, in all the results presented in this section, the units were normalized to the characteristic laser frequency, ω 0 and wave number used in the simulations, k 0 . The normalizations are thus t → ω 0 t and x → k 0 x. These normalizations of space and time define the normalizations used for the fields, i.e: E → eE/mcω 0 and B → eB/mcω 0 .

Vacuum birefringence with a static field
The birefringence of the vacuum is a thoroughly studied setup of great experimental interest to explore the properties of the quantum vacuum [38,39]. A one-dimensional wave packet traveling in the presence of a strong static field will experience a modified refractive index of the vacuum due to the HE corrections. To obtain an approximate analytical expression, one assumes that the strong background field remains unperturbed by the nonlinearities. This motivates the following ansatz for the solution of modified Maxwell's equations where E p and E s represent the electromagnetic pulse and the static fields, respectively and E p E s . Inserting equation (16) into the QED Maxwell's equations and keeping only the dominant terms in the polarization and magnetization, one obtains the following refractive indices: where the parallel and perpendicular directions refer to the direction of the probe polarization when compared to the static field. In the case of a constant externally imposed magnetic field, the expressions of the indices are swapped [34] (the value of the perpendicular index takes the value of the parallel index and vice-versa). Notice that the product ξE 2 s appears as a relevant quantity. This is a recurring property of several setups. It must be ensured that this product is a small quantity, both for the validity of the theoretical framework but also from the point of view of the algorithm. This quantity controls whether the corrections to the unperturbed fields are small or not, a crucial feature for the stability of the algorithm as already discussed. It is worth mentioning that the effective vacuum refractive indices only depend here on the externally imposed field to the order considered. Bialynicka-Birula [40] showed that if one considers the HE Lagrangian to all orders, a non collimated wave packet will experience, in addition to the external field effects, self-interaction. The dependencies of the index of refraction on the wave field will result in higher harmonics generation, which will eventually alter the shape of the wave pulse along the propagation.
The simulation setup consists of a strong static electric field of 10 −3 E s aligned along the y direction and a Gaussian EM pulse propagating in the x direction and polarized in the y − z plane. The central wavelength of the EM Gaussian pulse is 1 µm and its duration is 5.6 fs. Figure 6 shows two simulations for the same pulse after propagating once through a periodic box. In one case the propagation is in the classical vacuum, whereas the QED solver is used in the other. Qualitatively, the difference in propagation distance and the reduced electric field amplitude is consistent with the theory of a pulse traveling in a refractive medium. To test the accuracy of the algorithm this same setup was run for different values of the product ξE 2 s for both the parallel and perpendicular setup. The difference in phase velocity between the two pulses allows to extract directly the quantum vacuum refractive indices and to compare with the analytical predictions of equations (17) and (18). The results are shown in figure 7(a) and figure 7(b) where an excellent agreement between simulations and theory is found.

Vacuum birefringence: Optical pump and X-ray probe
The effects of the quantum vacuum on the propagation of light waves requires a very strong static electrical field that is unlikely to be produced in the laboratory. However it has been suggested a long time ago, for example by Brezin and Itzykson [34], that high intensity oscillatory fields with frequency small compared with the wave frequency might play the role of an external field. This is nowadays the aim of various experiments based on the interaction between a counter-propagating ultra-intense optical pulse and a low-amplitude X-ray probe pulse [4,13]. It should be emphasized that vacuum birefringence differs for finite times when considering an eternally-constant background (as in section 3.1) or an adiabatically-evolved, quasi-static background [41,42]. The latter configuration has motivated the numerical work presented in this subsection We start with the wave equation for the corrected QED Maxwell's equations [43,35] given by: where It is evident that the probe pulse propagation will be modified whenever the source term S is nonzero. This is verified when the EM invariants are nonzero and the strong field intensity approaches the critical value. If one assumes that the strong field remains unperturbed by the nonlinearities, then we can search for solutions of the HE Maxwell's where p and 0 stand for the probe and strong field quantities, respectively. Plugging Eq. (21) in Eqs.(6-7) and considering just the dominant terms, we can obtain approximate expressions for the resultant probe field after the interaction. This last step can be achieved by solving equation (19) using the Green's function method [35,43]. The increasing availability of high-intensity lasers renders these QED effects closer to being detectable in large-scale facilities. The prospect of coupling a 1 PW optical laser with an XFEL laser pulse will allow measuring the ellipticity induced by the QED vacuum non-linearities [13,35]. The setup shown in figure 8 comprises the most promising configuration for detecting the vacuum birefringence.
Optical pulse Figure 8. Counter-propagating setup proposed in [13] to probe the quantum vacuum.
For this setup, the one-dimensional approximate solutions for the probe (being a plane wave), considering just an Gaussian profile for the pump, are where, for θ p = π/4, with σ 0 being the optical pulse length, k p the wave number of the probe, and θ p the angle between the electric field polarization and the y axis. The expressions given by Eqs.(23a) already show us several important features. The first-order probe field components have a π/2 phase-shift relative to the zeroth-order components. This makes it explicit that QED induced components are contributing to a resulting ellipsoidal polarization [41,42]. The induced components scale with ξE 2 0 , which shows that this effect is only measurable for strong fields approaching the Schwinger critical field value. Dealing with an optical laser pulse of finite duration, the effective interaction distance is now proportional to its pulse duration parameter, σ 0 . The accumulated phase after the interaction can be written as [4] ∆φ = The ellipticity b is related to ∆φ [4,13], and reads, for a static field, where ∆φ = 12πk p dξE 2 0 and where d is the distance of interaction. For a probe interacting with an optical pump one obtains The static field result is recovered if we take into account: (i) a factor of 2 that appears from considering an optical pulse with no longitudinal oscillations; (ii) the limit σ 0 → ∞.
Up to now, we have only considered longitudinal effects. However, finite laser pulses possess transverse Gaussian profiles. We study how these Gaussian profiles applied to the strong and probe pulses modify the amplitude of the ellipticity as a function of several constraints applicable to realistic configurations: probe polarization angle, laser misalignment, timing and spatial jitter. We quantify these modifications through the variation of the ellipticity due to having a more realistic optical pulse duration and pulse diffraction. The standard simulation setup is the one presented in figure 8, where both laser pulses are focused at the same point in time and space and counter-propagate co-axially. The XFEL (X-ray pulse) is represented by a 3 fs pulse with 10 18 W/cm 2 and wavelength λ p = 10 nm. The optical pump pulse has the same duration with 10 23 W/cm 2 and λ 0 = 1 µm. We consider both pulses to have the same spot size W p = W 0 = 30 [c/ω 0 ] = 4.78 µm at the focal place x = 0. For a better illustration of the vacuum nonlinearity effect, we have increased artificially the coupling parameter to ξ = 4 × 10 −13 .

Probe polarization angle
Equation (26) shows that the ellipticity is proportional to sin (2θ p ), thus yielding its maximum amplitude when the X-ray probe polarization has a π/4 angle with respect to the strong optical pulse polarization. In addition, considering a two-dimensional paraxial Gaussian transverse profile, the ellipticity will vary along the y direction due to the transverse dependence of the fields with y 0 and y p being the center of the Gaussian profile along y for both laser pulses and z being the longitudinal distance from the focus of each pulse (z rp,0 = k p,0 W 2 p,0 /2 is the Rayleigh length). We consider here the special case where both pulses interact in the same point in time and space (z = 0). Figure 9 shows the transverse variation of the ellipticity for different probe polarization angles. On the other hand, the dependency of the ellipticity with the angle is verified in figure 10. The ellipticity value shown has been normalized to the best case scenario, corresponding to a transverse position y = y 0 = y p and a probe polarization angle θ p = π/4.

Laser misalignment
An important effect to take into account is the laser misalignment in the case where both pulses are not propagating along the same axis. In this case, the transverse profiles of the pulses do not overlap fully and we define as the misalignment parameter (also known as the transverse impact parameter). As equation (26) consists in the product of two transverse profiles, it is then possible to rewrite the dependency of the ellipticity on the transverse coordinate as It is important to notice that we have written equation (29) as an effective transverse profile, hence, the expression forb contains only the central amplitude of both pulses. Figure 11 shows the transverse variation of the ellipticity for different values of misalignment, while figure 12 shows how the maximum for each case vary with the misalignment parameter. As previously, the ellipticity amplitude is normalized by the best case scenario, corresponding to a coaxial counter-propagation (m = 0) and, thus, y = y 0 = y p . The probe polarization angle is considered to be θ p = π/4. These results are consistent with the ones obtained in section III-D of [41]. In this paper, the authors show that the full exponent for the ellipticity amplitude should decay with −2ρ 2 /(1 +ω), where ρ = r/W 0 andω = 2W 2 p /W 2 0 . If we identify m = r, following the notation of [41], we verify that it is equivalent to the two-dimensional R decay expression in equations (29) and (30d).

Timing and spatial jitter
Potential fluctuations in the alignment of the lasers can lead not only to laser misalignment but also to jitter [13]. The probe pulse can be off focus with the optical pulse both in space (spatial jitter) and in time (timing jitter), as shown in figure 13. This effect has a direct impact on the integrated amplitude of the optical pulse as seen by the XFEL probe pulse during the overlap phase. Consequently, the ellipticity can be written as in equations (29)-(30d), with W 2 0 replaced by W 2 0 1 + (z/z r0 ) 2 . The small corrections were not taken into account in equation (27a) since the XFEL Rayleigh length is much larger than the interaction length (approximately the optical pulse width). Figure 14 shows the transverse variation of the ellipticity for different spatial jitter configurations, while figure 15 shows how the maxima for each case vary with the distance to the focus of the optical pulse. In this last figure, we also superimpose results for timing jitter configurations, verifying that the results follow the same trend. The ellipticity amplitude is as usual normalized by the best case scenario, corresponding to a coaxial counter-propagation (m = 0), probe polarization angle θ p = π/4, and nonexistent jitter (z = 0). In section III-C of [42], the authors explore these effects in a three-dimensional configuration. This implies that the Gaussian paraxial fields will decay faster, hence, yielding that the ellipticity should decay with 1/(1 + τ 2 ), where τ = z/z r0 , instead of our square root decay. Our results and our numerical solutions are consistent with the predictions of [42] further validating our algorithm.

Counter-propagating plane waves
Two counter-propagating plane waves polarized in the same direction, with the same frequency ω 0 and amplitude would normally result in an electromagnetic standing wave in the classical vacuum. However, the vacuum nonlinearities lead to the well-known phenomena of harmonic generation [37,44]. This simple setup represents an ideal benchmark for the numerical algorithm as analytic results can be obtained. To address the problem, it is convenient to decompose the field into a Born series of partial waves as performed by Bohl et al. [35]  with ω 0 = k 0 c. The electric and magnetic fields are aligned with the y and z axis, respectively. Starting from the modified Maxwell's equations and inserting Eqs. (31)(32) as the expressions for the fields, we arrive at the wave equation for the first-order correction to the electric field E (1) where is the d'Alembert operator and the source term These are just the one-dimensional versiosn of equations (19)- (20). Inserting the zero order field in the source term, i.e. taking P, M = f (E (0) , B (0) ), we find This source term only accounts for the unperturbed fields being inserted into the nonlinear polarization and magnetization. The formal solution of this equation is given by the convolution between the source term and the Green's function of the one-dimensional wave operator [45], where the Green's function is The modified electric field reads We notice that the relative amplitude between E (1) and the unperturbed field amplitude is proportional to ξE 2 0 showing again that this perturbative treatment is valid as long as ξE 2 0 1. More specifically, the corrected field exhibits a secular growth term modulated by an oscillating term. This term is dominant for t 1 and should be interpreted as a phase shift due to the induced birefringence of one wave to the other. The total field for t 1 reads E E (0) + 8ξE 3 0 t sin(t) cos(x). Using the trigonometric identity a cos(x) + b sin(x) = √ a 2 + b 2 cos(x − arctan(b/a)), one can write the total corrected field as where n −1 = 1 + 4ξE 2 0 and n the modified refractive index induced by the interaction of the two waves. Taking the spatial Fourier transform of E (1) , we verify that the fundamental mode k = k 0 is corrected by a secular term and the appearance of an harmonic at 3k 0 . Defining the Fourier transform of E(x, t) asẼ(k, t), we obtaiñ The third harmonics correspond to two waves, cos[3(x+t)] and cos[3(x−t)], propagating with the initial zeroth-order field. The next order correction to the field E (2) , reveals a correction to the k 0 mode growing as t 2 , a secular 3k 0 harmonic and an oscillating 5k 0 harmonic. Repeating this process to higher orders, we can show that this nonlinear interaction generates odd higher harmonics from vacuum with the relative amplitude between these harmonics obeying the ordering Nonetheless, it should be emphasized that a rigorous treatment of higher harmonics (beyond the first-order correction) should take into account additional terms in the expansion of the Heisenberg and Euler Lagrangian for weak fields (E < E sch ). As shown by Bohl et al. [35], purely four-photon scattering (first-order term of the Euler Heisenberg Lagrangian) allows the generation of higher harmonics in the counter propagating setup. However the contribution from this twice-iterated process scales as (ξE 2 0 ) 2 E 0 , which when compared to the leading contribution to the fifth harmonic from six photon scattering is suppressed by a factor of ξE 2 0 . Our analytical predictions were compared with the results of the QED solver using a field amplitude of E 0 = 0.025E sch , λ 0 = 1 µm plane waves and ξ = 10 −9 , such that the higher harmonics can be accurately resolved above the numerical noise. The spatial Fourier transform of the fields is shown in figure 16 for two simulations, with and without the self-consistent inclusion of Heisenberg-Euler corrections. We observe that when the non-linearities are present, the odd higher harmonics are generated with a relative amplitude that matches the ordering given in equation (41). To compare the simulation results with equations (39)-(40), we subtracted the classical vacuum electric field to remove the zeroth order standing wave contribution, and performed the Fourier transform of this subtracted field. Finally, we tracked the temporal evolution of the amplitude of the k 0 mode in Fourier space and compared it with equation (39). Figure 17 shows the temporal evolution of E 1 (k 0 ). The simulation shows an excellent agreement with the theoretical predictions for many laser cycles, ensuring that the algorithm is robust. Despite the one-dimensionality of this example, a setup of counterpropagating beams is of great interest for planned experiments at extreme high intensity laser facilities, as outlined in [43].

Interaction of paraxial beams: counter propagating setup
In order to illustrate the generation of harmonics in multi-dimensions, two setups were investigated: the counter propagation of two Gaussian pulses interacting at the focal point, and the perpendicular interaction of two Gaussian pulses focused in the same region of space. For these setups, a consistent analytical treatment becomes cumbersome especially due to the self-consistent treatment of both the transverse and longitudinal components of the pulses. A quantum parameter of ξ = 10 −7 was used for the sake of showing the prominent features of the harmonics of very small amplitudes.
In the first setup, two λ = 1 µm laser beams with a normalized vector potential a 0 = eE 0 /mcω 0 = 50 (which corresponds to laser electric field at the focus of E 0 10 −4 E s ) and duration of 25 fs interacted in the presence of vacuum non-linearities. Both beams have a focal spot of W 0 = 4 µm. Figure 18 shows the transverse electric field of the laser beams before interaction at t = 0, figure 19-(a) the spatial Fourier transform of the beams with ξ = 0 (classical limit) and in figure 19-(b) the Fourier transform of the electric field after the interaction (asymptotic state) including the HE corrections. As shown in figure 19-(b), after the interaction odd higher harmonics are also generated as in the 1D case, with relative amplitudes consistent with Eq. (41). However, in this case, the harmonics generated have the same Gaussian behavior as the unperturbed pulses and attain a greater spread in Fourier space after the interaction. After the pulses have spatially overlapped, the harmonics propagate and leave an imprint of the nonlinear interaction, that co-propagates with the original beams.

Interaction of paraxial beams: 90 • setup
The second setup, shown in Fig 20, comprises two optical Gaussian pulses interacting at 90 • . The two laser beams possess the same parameters: a 0 = 100, a wavelength of 1 µm, a focal spot W 0 = 4µm, a duration of 25 fs. The advantage of this setup is the vast amount of harmonics generated during the interaction of the two pulses. Before discussing the results of the simulations, the reader can develop a valuable intuition of the generated harmonics by carefully computing the electromagnetic invariants and the associated vacuum polarizations for paraxial beams. The theory of paraxial electromagnetic fields has been developed by Davis [46] with a simple method that allows to find a formal solution of a light beam propagating in classical vacuum. The formal solution is based on an expansion in powers of a small parameter s = W 0 /l r = 1/kW 0 where W 0 is the beam waist and l r = kW 2 0 the Rayleigh or diffraction length. For the special case of a two-dimensional beam, varying as e iω 0 t propagating in the x direction and polarized in the y direction, the non vanishing components up to the second order in s are where ω 0 = k 0 c, x = ζl r , y = W 0 η and where σ is the typical duration of the beam and a 0 = eA 0 /mc 2 is the Lorentz invariant parameter which measures the magnitude of the field. Whereas Gaussian paraxial beams are usually described up to the first order in s [47,44], the inclusion of the second-order terms for the transverse components are essential to calculate accurately the electromagnetic invariants which consist of a series of terms proportional to (k 0 A 0 ) 2 , s(k 0 A 0 ) 2 and (sk 0 A 0 ) 2 . In our setup, we refer to the index 1 for the pulse propagating in the x direction with wavenumber k x and to the index 2 for the pulse propagating in the y direction with wavenumber k y . The electromagnetic fields for the pulse 2 can be found by rotating by 90 • the fields described for the pulse 1. The non vanishing invariant is In order to highlight the harmonics generated during the interaction of the two pulses, we use the simplified notation (n, m) ≡ e i(nkxx+mkyy) with n, m ∈ Z. In Fourier space, the invariant is symmetric with respect to the k x and k y axes and we can thus just consider for the sake of simplicity one quadrant of the k space. Keeping only the terms for which k x , k y > 0, the polarization P = 4ξFE calculated at first-order reads The full expression for the invariants F x , F y , F z and the polarisation coefficientP x ij ,P y ij can be found in the appendix Appendix A. The time at which the interaction is strongest occurs at the full overlap of the two pulses, ω 0 t = 100 for our configuration. The Fourier transform of the polarisation P x at that time is shown in figure 21(a). The harmonics predicted theoretically in equation (46a) can be readily identified as well as new harmonics such asP x 23 orP x 32 that result from highest order coupling. The first order harmonics of largest amplitude areP x 10 ,P x 12 (four others are just the symmetric harmonics with respect to the k x and k y axis) and are proportional to ξ(k 0 A 0 ) 3 whilst P x 01 ,P x 21 scale as ξs(k 0 A 0 ) 3 . The harmonics of lowest amplitudeP x 30 do not arise from the interaction of the two pulses but from their self-interaction and are thus proportional to ξ(sk 0 A 0 ) 3 . A more precise comparison between the simulation and the first order theoretical harmonics of the polarisation P x is shown on figure 21(b1-b2-b3-c1). One notes the very good agreement both for the respective amplitude of the harmonics and their shapes in Fourier space. The theoretical calculation of the first order fields can be carried out by convoluting the 2D Green function of the wave propagator with the linearised source term [7,37,43] comprised of partial derivatives of the first order polarisation and magnetisation. We have plotted in figure 22 the temporal evolution of the Fourier transform of the electric field E x . At t = 0, the two pulses are fully separated which implies that no interaction has started yet. As a result, only the central wavelength of each pulse is visible, k x for the longitudinal field of pulse 1 and k y corresponding to the transverse component of pulse 2. When the two pulses fully overlap, one identifies several harmonics that are identical to the ones we have described for the polarization, which is a direct consequence of the linear property of the wave operator (the polarisation being the source term of the wave equation). Nonetheless, the relative amplitude between the harmonics of the electric field differ from the one we have observed for the polarisation. This is somehow obvious since we are showing here the total electric field E x which is the sum of all corrections due to the polarisation and magnetisation of the vacuum. Finally, at time ω 0 t = 200, the two pulses have left the zone of interaction, thus ceasing to feed the nonlinear interaction between them. The remainder contributions stem from the self-interaction alone [36]. As also identified by other authors [35,26], off-axis contributions dominate the spectral region during the overlap and occur due to the field spatio-temporal inhomogeneities. These contributions rapidly fade away with the separation of both pulses. The harmonics that persist after the two pulses have fully separated stem from the self-interaction of each pulse. The amplitude of these harmonics appears to fade away contrary to the 1D case for which the amplitude remains constant. In a 1D simulation, the amplitude of the harmonics generated during the interaction of two pulses does not decrease as they move out of the interaction zone, as seen before whereas in 2D the amplitude of the signal goes down as 1/r (and 1/r 2 in 3D).
The Fourier spectra obtained in these two setups show that the harmonics generated in either case are distinct, thus allowing to clearly distinguish both cases. Future work will include the analytical study of the relative intensity and spectral width of the generated harmonics and their possible relation with other beam parameters. Namely, it is of great interest to understand how the production of these higher harmonics from vacuum may be optimized in terms of the duration of the pulses as these results can provide signatures of experimental relevance. A future setup to explore will also include the interaction of two laser beams at an arbitrary angle 0 < θ < π 2 to model realistic experimental conditions. If this angular dependence of the interaction is well understood, one could in principle determine how well aligned two ultra-intense beams are by examinating at the Fourier spectrum after a vacuum interaction. Finally, the theoretical predictions, made in the case of two intense focused beams overlapping, on photon merging/splitting [47] and four wave mixing [44] could also be verified with an extension of this present code in 3D dimensions, which is computational very demanding but still feasible. A straightforward extension to the code could also be the addition of higher order terms in the Euler Heisenberg Lagrangian. This would allow us to explore the shock formation and asymptotic field generated via six-, eight-or higher wave mixing processes [26,35].

Conclusions
A numerically stable and robust generalized Yee scheme to solve the nonlinear set of QED Maxwell's equations was developed and incorporated in a standard PIC loop. This work represents an important step towards modeling plasma dynamics in extreme scenarios when QED processes significantly alter the collective behavior of the system. Furthermore, the algorithm is fully generalizable to include higher order corrections (such as six-photon scattering or higher order terms). These terms are to be included in the future as they have been shown to be necessary to fully simulate certain scenarios [35]. Our numerical model can be used to design planned experiments, leveraging on ultra-intense laser facilities able to deliver intensities of 10 23 − 10 24 W/cm 2 , to verify for the first time the dynamics of the quantum vacuum below the Schwinger limit. The simulations confirm predicted optical phenomena such as vacuum birefringence and high harmonics generation in one-dimensional setups with an excellent accuracy. The algorithm was also extended for two-dimensional scenarios where two setups of interacting Gaussian beams were studied. The results highlight the importance of transverse beam effects and hint that the generation of higher harmonics from quantum vacuum can be achieved via this interaction. The spectrum of the harmonics could provide a direct measurement of important beam properties such as the peak intensity and alignment. This algorithm may also be used to test two and three dimensional setups that have been proposed in the literature (where transverse and finite spot size effects are taken into account under certain approximations), thus complementing the results of previous theoretical works [41,42,43]. Finally our algorithm contributes to the generalization of the Yee scheme, one of the most successful and commonly used algorithms in computational physics, to scenarios where nonlinear polarization and magnetization can impact EM propagation. This