Numerical Simulations of the Nonlinear Quantum Vacuum in the Heisenberg-Euler Weak-Field Expansion

The Heisenberg-Euler theory of the quantum vacuum supplements Maxwell's theory of electromagnetism with nonlinear light-light interactions. These originate in vacuum fluctuations, a key prediction of quantum theory, and can be triggered by high-intensity laser pulses, causing a variety of intriguing phenomena. A highly accurate numerical scheme for solving the nonlinear equations due to the leading orders of the Heisenberg-Euler weak-field expansion is presented. The algorithm possesses an almost linear vacuum dispersion relation even for comparably small wavelengths and incorporates a nonphysical modes filter. The implemented solver is tested in one spatial dimension against a set of known analytical results for vacuum birefringence and harmonic generation. More complex scenarios for harmonic generation are demonstrated in two and three spatial dimensions.


Highlights
• A universal 13 th order of accuracy numerical scheme with nonphysical modes filter • Universality with respect to pulse configurations in contrast to analytical treatments • Scalability on distributed computing systems • Inclusion of up to six-photon processes in the Heisenberg-Euler weak-field expansion Email addresses: and.lindner@physik.uni-muenchen.de(Andreas Lindner), b.oelmez@physik.uni-muenchen.de(Baris Ölmez), hartmut.ruhl@physik.uni-muenchen.de(Hartmut Ruhl)

Introduction
Virtual electron-positron pairs are omnipresent in the vacuum of quantum electrodynamics (QED) and can mediate effective photon-photon interactions.Taking those into account in the low-energy limit below the Compton scale of the electron, nonlinear terms in the electromagnetic field extend Maxwell's theory, invalidating the classical superposition principle in the vacuum.Radiation on the other hand can influence the behavior of the particle-antiparticle dipoles in the vacuum.Therefore, the vacuum of QED with its quantum fluctuations acts as a nonlinear, polarizable medium with quantifiable effects.The corresponding low-energy effective theory was devised by Euler, Kockel, Heisenberg, and Weisskopf [1][2][3].Schwinger incorporated the theory into the larger QED picture [4].For reviews, see [5][6][7][8][9][10][11][12][13][14].
A promising approach to their detection is the investigation of the asymptotic dynamics of probe photons after passing through a strong-field region in the form of a high-intensity laser pump pulse.Probe photons can in this way indirectly sense the applied pump field via the quantum fluctuations which couple to both probe and pump fields [40].As a result of the nonlinear interaction with the power pulses, signal photons of the probe field scatter in such a way that their dynamics or polarization distinguishes them from the main pulses.The properties and dynamics of the quantum vacuum, in turn, are hereby encoded into the signal photons.This class of quantum vacuum experiment, where one electromagnetic field drives the nonlinear effect while the other carries its signature, is called all-optical.
Experiments at the high-intensity frontier promise unprecedented detail in the study of the quantum vacuum in the near future.The hope is that the rapid and steady advancements in ultra-intense laser physics [79,80] combined with theoretically optimized specifications will soon facilitate the experimental discovery [7,10,12].
Almost all analytical approaches to compute nonlinear vacuum effects, as found in the literature cited above, have shortcomings.As a consequence, these attempts are still limited to simple scenarios, which impedes experimental verification.Approximations limit the accuracy of predictions and the precision the theory can be tested with.Precision tests require accurate theoretical predictions for arbitrary laser field configurations [81].
The present paper introduces a numerical algorithm for the solution of the Heisenberg-Euler equations in weak-field approximation, where "weak" means below an ultra-high critical field strength of more than 10 18 V m −1 , to be detailed below.The work is based on the numerical scheme outlined in one spatial dimension in [66][67][68] and extended to full three dimensions in [82,83].The accuracy of the numerical scheme relies on the assumptions that the fields vary on scales larger than the Compton scale of the electron and the field strengths involved are below the critical field for pair creation in the vacuum.In addition, a viable numerical algorithm is restricted by the computational load it generates, particularly in the realm of high frequencies.The advantage of a numerical approach, however, is that it can principally solve almost any interaction scenario for nearly any arrangement of laser pulses.
The numerical predictions of the presented solver are benchmarked with analytical predictions for birefringence as calculated in [37,38] and the generation of higher harmonics as predicted in [66].

Comparison to other approaches
There exist other approaches to simulate all-optical QED vacuum effects.The algorithm outlined in [84] and the solver presented in [81] are based on the vacuum emission picture [85].In the vacuum emission picture the fields are split into strong background and weak signal fields and the Heisenberg-Euler Lagrangian is expanded to linear order in the signal field strength.The strong electromagnetic fields driving the nonlinear effects are propagated by the linear Maxwell equations in the absence of vacuum nonlinearities.This is done by means of a Maxwell solver detailed in [86].
The vacuum subjected to the high-intensity fields constitutes a source term for an outgoing photon via a scattering amplitude using the Heisenberg-Euler interaction Lagrangian.The picture can be interpreted as describing laser-stimulated signal photon emission from the vacuum.In certain configurations excellent signal to background separation can be achieved.
However, there is no feedback on the effect-generating electromagnetic fields themselves by the response of the quantum vacuum.In other words, nonlinearities induced by laser fields, and thus also back-reactions into the latter, are neglected.Hence, there is no notion of pump and probe laser fields in contrast to the scenarios considered throughout the present work.
The vacuum emission picture is capable of making predictions for asymptotic states of ultra-short emission wavelengths, while the solver outlined in the present paper is limited by the affordable grid resolution.Numerous analytical works are based on the vacuum emission picture [38,40,47,[87][88][89][90][91][92][93].In addition, the solver based on the vacuum emission picture has already been employed in some studies of optical signatures of the quantum vacuum [27,94,95].
Besides vacuum birefringence, the approach has also been employed to photon-photon scattering as well as merging and splitting processes.The emission process is not restricted to cubic order in the background field [57,58].Hence, multi-photon emission processes can also be incorporated in the description.
In [96] a generalized Yee scheme of second order accuracy in space and time is devised.The accuracy order of the numerical scheme discussed in the present paper is arbitrary and implementations for the second to the thirteenth order are given.The Yee scheme approach in [96] requires interpolations in space and time to compute the vacuum nonlinearities.Also, the dispersion relations of the numerical scheme in [96] and in the present paper differ from each other.The dispersion relation outlined in [96] has an imaginary part that can lead to the amplification of nonphysical modes and requires a high grid resolution for a given wavelength.The dispersion relations outlined in the present paper have imaginary parts that always damp nonphysical modes and at high integration order can afford lower grid resolution.
The solvers in [96] and [81] include only the four-photon order of the weak-field expansion of the Heisenberg-Euler interaction.Some all-optical vacuum nonlinear effects are tiny.Hence, high-order, high-precision numerical solvers may bear advantages.An implementation for distributed computing systems of the scheme discussed in the present paper is available and allows costly simulations in full three spatial dimensions plus time.High-resolution grids, however, are subject to the curse of dimensionality.Some applications require very high resolutions in order to model high-frequency waves.These easily go beyond the limits of any computing system and thus extrapolation techniques have to be employed.

Outline
The paper begins with a recapitulation of the weak-field expansion of the Heisenberg-Euler effective theory in Section 2. In Section 3, the nonlinear modifications to Maxwell's equations due to the Heisenberg-Euler theory are detailed.In Section 4 the numerical scheme that solves the nonlinear Maxwell equations is outlined.Dispersion effects on the lattice and the scaling behavior of the code implementation are discussed in Sections 5 and 6.
In Sections 7-10 the capabilities of the Heisenberg-Euler solver is demonstrated by solving various scenarios of nonlinear vacuum phenomena and by successfully benchmarking with analytical results, where they exist.In Section 7 the phase velocity reduction of a probe pulse propagating through a strong electromagnetic background is successfully compared to analytical predictions for sufficiently large background field strengths.In Section 8 the phenomenon of vacuum birefringence is investigated with the help of simulations for probe-pump setups.The theoretical prediction of the polarization flipping probability and its scaling with various parameters is verified.An extrapolation of the results to wavelengths in the x-ray regime is performed.
In Section 9 the capability of simulations to predict higher harmonics generated by vacuum nonlinearities is demonstrated.It is shown that for specific settings analytical results from [66] are accurately reproduced in the simulations presented in this work.Examples of harmonic generation in 2D and 3D simulations are presented in Section 10.The latter are hard, if not impossible, to obtain by analytical means.The results are used for comparison with the previous 1D simulations.Higher-dimensional simulations enable the investigation of complicated scenarios that cannot be solved analytically in the future.
A summary and an outlook to future research employing the Heisenberg-Euler solver are given in Section 11.

Heisenberg-Euler weak-field expansion
The Heisenberg-Euler Lagrangian can be expressed as [97] where and e, m e are the charge and mass of the electron.The quantities F and F * denote the electromagnetic field strength tensor and its dual and E cr is the above mentioned critical field strength, which is the scale where the magnitude of the vacuum nonlinearities becomes dominant.
The Lagrangian for the quantum vacuum is then given by the classical Maxwell Lagrangian and the quantum corrections due to vacuum polarization captured in the Heisenberg-Euler Lagrangian ( The Maxwell Lagrangian is given by with the vacuum permeability To derive L HE it has been assumed that the constant background approximation is valid [98,99] as a consequence of the fact that the variations of the field strengths are small on the scale of the Compton wavelength [100] The weak-field expansion of (1) up to fourth order in F and G is given by with the vacuum permittivity and the fine-structure constant It should be noted that the Heisenberg-Euler Lagrangian solely depends on the electromagnetic field invariants F and G , and subsequently (9a) is of the order O ((E/E cr ) 4 ).
The effective interaction L HE can be graphically represented by a double-lined loop as in Figure 1.As a consequence, (9a) represents the four-photon contributions, (9b) the six-photon contributions, and (9c) the eight-photon contributions to the closed loop as illustrated in Figure 1.Since the field strengths are constrained to below the critical field strength, c|F µν | < E cr , higher order terms can be neglected in the low-intensity regime.
It has to be mentioned that pair production is exponentially suppressed in the Heisenberg-Euler Lagrangian.Therefore, the properties of the nonlinear vacuum can be described exclusively by radiation fields.

Modified Maxwell equations
To obtain the nonlinear modifications to the linear Maxwell equations an explicit rescaling of the electric and magnetic fields ⃗ E and ⃗ B to the critical field strengths is instructive such that the field strengths are given in units of E cr ,

≈ + +
Figure 1: Graphical representation of the Heisenberg-Euler Lagrangian in weak-field expansion.The double-lined loop to the left represents L HE and corresponds to the irradiated electron-positron loop, denoting the coupling of the external field to the fermion loop to all orders in the field strength.The diagrams to the right are the nonlinear four-, six-, and eight-photon contributions due to the weak-field expansion of L HE , corresponding to the terms in (9a), (9b), and (9c).
From the Lagrangian of the quantum vacuum (5) the equations of motion are obtained with the help of the Euler-Lagrange equations with respect to the normalized fields with the definition The classical Maxwell-Ampère circuital law is compared to (13) in such a way that photon-photon interactions are included, which implies the appearance of macroscopic nonlinear polarization and magnetization terms given by The modified Maxwell-Ampère law (15) with the terms in ( 16) is merged with the Maxwell-Faraday law of induction which persists in the nonlinear vacuum, such that a single partial differential equation that describes the whole dynamics of the system is formulated.The curl of ( ⃗ B − ⃗ M) can be written as With the help of the vector the Maxwell equations ( 15) and ( 17) can then be expressed as with the 6 × 6 identity matrix 1 6 and the matrices A and Z j are given by where 0 3 is the 3 × 3 zero matrix and with matrix elements Equation ( 20) contains the full dynamics of electromagnetic fields in the weak-field approximation of the Heisenberg-Euler model.For illustrative purposes, first the linear vacuum is considered by setting ⃗ P = ⃗ M = 0, leading to To deduce modes with specific propagation directions, which becomes important in the construction of the numerical scheme, Z lin j is diagonalized using the matrices so that for Z lin j the expressions are obtained.Now, Equation (24) becomes where Equations ( 27) and (28) imply backward propagation in x-direction of the field components (E y − B z ) and (E z +B y ), forward propagation of the field components (E y +B z ) and (E z −B y ), and non-propagation of the first and fourth component.The components for the other propagation directions can be read off analogously.

From PDE to ODE
For demonstration purposes only the x-direction is considered in the following.In order to transform the partial differential equation into an ordinary one, finite difference approximations are introduced.The spatial derivative ∂ x is replaced by finite differences where The S ν are stencil matrices and ∆ x is the spatial resolution.The stencil matrices S ν are derived with the help of a Taylor expansion of a generic function g(x + ν∆ x ) around ν∆ x = 0.The Taylor expansion up to accuracy order six is given by (ν∆ x ) 5 g (5) (x) + 1 720 (ν∆ x ) 6 g (6) (x) + O (ν∆ x ) 7 . ( For the finite differences approximations of derivatives, forward and backward steps can be made use of.For the values ν ∈ {−3, ..., 3} the expansion in matrix form reads x g ′′ (x) ∆ 3 x g ′′′ (x) ∆ 4 x g (4) (x) ∆ 5 x g (5) (x) ∆ 6 x g (6) x) and can be extended to larger ranges of ν and inverted to obtain approximations for the derivatives taking into account more distant points.Recall that only the first derivative is required for the merged equation of motion (20).
To illustrate how the finite difference approximation of the first derivative is calculated, ν ∈ {0, 1} is considered for simplicity.It is obtained which leads to the first derivative of g given by The derivative g ′ ν∈{0,1} in (34) is called first order forward discrete derivative.In the same way the first order backward discrete derivative g ′ ν∈{−1,0} is obtained with the help of where In order to obtain specific dispersion properties and numerical stability of the discretized version of Equation (27) it is necessary to introduce forward and backward finite differences for D x , as discussed further below.With Equation (30) the full expression for the spatial derivatives for distinct propagation directions is formulated Setting aside non-propagation of modes, the first three elements of the rotated field vector R x ⃗ f repre- sent forward-propagating electromagnetic modes, while the last three components of the latter represent backward-propagating ones, as discussed at the end of Section 3.This understanding is adopted for the discussion of dispersion relations in Section 5.
For the nonlinear vacuum the system of ODEs is obtained from the merged equation of motion (20) and Equation (37) to be The inversion of (1 6 + A) is an expensive operation.An approximation in terms of a truncated geometric series is not satisfying since also larger nonlinear corrections ought to be taken into account.The problem of inverting the matrix can be reduced from a 6×6-dimensional problem to a 3×3-dimensional one by making use of the block form in (21) such that [83] ( with This 3×3 matrix can be inverted explicitly.An external numerical solver introduced in Subsection 6.1 is responsible for the time integration of the ODE system.

Simple examples
To investigate dispersion effects in the finite differences approximation an elementary example is worked out in the following.For the values ν ∈ {−1, 0, 1} the expansion of g(x) in matrix form is given by It has to be noted that Equation (41) is not restrictive to one specific kind of finite difference method.
To obtain an expression for the first derivative the matrix in (41) has to be inverted.This leads to Assuming that g(x) is one of the two polarizations of an electromagnetic mode that propagates in positive (+) or negative (−) x-direction in 1D, the corresponding equation of motion reads Making use of (42), Equation (43) becomes for symmetric forward and backward differentiation where the spatial derivative of g is accurate up to second order in ∆ x .Equation ( 44) can be solved analytically by a plane wave ansatz where k + > 0 and k − < 0 to obtain However, the derivative in (43) can also be approximated for a forward-propagating mode g + by a first order forward finite difference.In this case, Equation ( 44) can be written as where the solution of (48) with ω ∈ C reads Figure 2: Dispersion relations obtained with loworder finite differences for a simple plane wave.
Here c = 1.Top left: for a second order symmetric forward and backward finite difference scheme.
Top right: for a first order forward finite difference for the forward-propagating mode and a first order backward finite difference for the backwardpropagating mode.Bottom: for a first order forward finite difference for the backward-propagating mode and a first order backward finite difference for the forward-propagating mode.The linear vacuum is given by the black line as reference with the vacuum speed of light c set to unity.Note that all stencils result in a symmetric dispersion relation.The phase velocity is the same for both directions.
Further, the derivative in (43) for the backward-propagating part can be approximated by a first order backward finite difference.In this case, Equation ( 44) can be approximated for the backwardpropagating case as with the solution The dispersion relations connect a simulated wave with a specified wavelength to its phase velocity on the lattice.While in the case of symmetric differentiation there is no imaginary part of ω, it does show up for biased differentiation.The derived dispersion relations are shown in the top row of Figure 2. The Nyquist frequency f Ny = ∆ −1 x /2, corresponding to k • ∆ x = π, marks the point where Re(ω) = 0 for k ̸ = 0 in all cases.

The scheme at fourth order
The imaginary part amplifies the wave and thus makes this scheme unstable.By differentiating biased against the propagation direction, however, the imaginary part switches signs and has a damping effect.The calculations are straightforward and the result is shown at the bottom of Figure 2.This behavior can be tuned with higher order schemes in order to get a more vacuum-like real part and to defer the imaginary part such that it becomes relevant only for short wavelengths and damps only those.As the real part of ω eventually decreases in any scheme as it reaches the Nyquist frequency, the damping of those high-frequency waves serves as anti-aliasing effect.This is shown below.Now, the unbalanced coefficients for the first order discrete derivative in fourth order accuracy in ∆ x , g ′ ν∈{−3,1} and g ′ ν∈{−1,3} in the above notation, are given with the help of Equation ( 32).This results in To make the connection to Equation ( 27), g is replaced by the components of R x ⃗ f .As explained above, R x ⃗ f has components propagating in different directions that are differentiated with biases in the re- spective opposite direction.This means that the first three components of R x ⃗ f are differentiated in the same way as g in the first line of ( 56) and the last three components of R x ⃗ f are differentiated as g in the second line of (56).Following these instructions for Equation (30) yields with the fourth order stencils given by 0 = diag(−5/6, −5/6, −5/6, 5/6, 5/6, 5/6) .
In consequence of the symmetry and redundancy, it is instructive to rewrite the stencil matrices in terms of its components applied to forward-and backward-propagating modes.The fourth order stencil matrices can be expressed as where the components applied to forward-and backward-propagating modes s 4 f and s 4 b are given by and s 4 f /b = 0 for out-of-range values of ν.Due to the obvious symmetry, the stencil matrix can also be written as for the whole range ν = −3, ..., 3.

Dispersion relations at order four and thirteen
The solver has implementations of the scheme up to order thirteen.In the present work the currently maximal available accuracy is used.By going seven steps of ∆ x forward and six steps backward in the forward-biased differentiation and vice versa in the backward-biased case, the unbalance is kept very small.The components of the low-bias stencils up to order thirteen are listed in Appendix A.
For a single plane wave the Heisenberg-Euler Lagrangian reduces to the Maxwell Lagrangian, as in this case F = G = 0 holds.The vacuum is thus the linear one.Inserting a plane wave, without loss of generality propagating on the x-axis, into the propagation equation ( 38) yields Since the stencil matrices are diagonal, this can be written as which, when multiplied with R x from the left, gives Inserting the stencil matrices expressed in terms of their components applied to forward-and backwardpropagating modes s f and s b results in It can be seen that the equation above holds irrespective of the axis of propagation since R x can be canceled.Hence, the equation becomes There are the distinct cases of forward-and backward-propagating modes, and non-propagating ones.
It is obtained for • a forward-propagating mode (implying k > 0): • a backward-propagating mode (implying k < 0): • a non-propagating mode (implying k = 0) : The stencil components to all orders up to order thirteen are listed in Appendix A and plots of the dispersion relations are given.The results for order four and thirteen are visualized in Figure 3, the latter being the ones used for simulations of this work.The minimal resolvable distance ∆ is the spacing between the lattice points, the resolution of the lattice.It is the decisive parameter to tune the dispersion effects for a given wavelength and is given by the ratio of physical length and lattice points.The real parts of ω start in the vicinity of the black line and deviate from it for shorter wavelengths or smaller grid resolutions.There is less deviation at higher orders.The imaginary part of ω starts with values close to zero for large wavelengths/high grid resolutions and decreases until the Nyquist  frequency is reached.The negative imaginary part in this scheme has the positive effect to promptly annihilate those modes that deviate strongly from the vacuum dispersion relation.At higher orders, the damping effect of the imaginary part is deferred to shorter wavelengths and is overall smaller.It can further be seen that for higher orders the real part stays closer to the linear vacuum for smaller wavelengths.Higher orders also defer and decrease the rise of the imaginary part.A relatively strong bias as in the fourth order scheme (56) causes a superluminal phase-velocity in a given k • ∆ range, a critical regime where the dispersion relation deviates remarkably from the linear vacuum.The well balanced order thirteen scheme has a smaller critical regime and is well-behaved for a larger range of wavelengths.
For small wavelengths the curve showing the real part of ω falls off and the imaginary part of ω causes a damping.As a result, there are two k's for each Re(ω).At the wavelength corresponding to the Nyquist frequency, k • ∆ = π, Re(ω) becomes zero as in the simple scenarios above.The highk values cause nonphysical standing waves, as Re(ω)/k → 0, which would cause a self-heating of the system.The damping takes care of their annihilation.The superluminosity which can be seen in the fourth order finite difference differentiation is thus an acceptable add-on to the favor of a damping of nonphysical modes.With increasing accuracy and balancing as in the order thirteen scheme there is no superluminal range left.

Damping of modes
In the numerical investigation it is demonstrated that the damping has noticeable effects already at an earlier stage than a first look at the above plot would suggest.Note that the Nyquist frequency, on the other hand, f Ny = ∆ −1 /2 (corresponding to k • ∆ = π), is not the limiting factor for wave modeling in this scheme.It has to be stayed well below the Nyquist limit for accurately time-evolved waves.Note also that the simulation of a high-frequency wave does not overshoot the Nyquist limit, as it would be sampled as a wave with lower frequency, see Figure 4.The relevant frequency scale of the dispersion relation in Figure 3 thus ranges only to the Nyquist limit and the scheme is generally stable for any frequency.
It has to be noted that energy is not conserved on the grid as a consequence of the amplitude damping during the propagation.The imaginary part of ω and therefore the damping effect might be small if the grid resolution is high, but there is a trade-off since high grid resolutions are computationally more expensive.
To visualize the effects of the dispersion relation at varying wavelengths the propagation of a plane wave in one direction of a two-dimensional square grid with side length 80 µm divided into 1024×1024 points is investigated.The resolution is therefore given by ∆ = 80 µm/1024, ∆ −1 = 128 × 10 5 m −1 .The grid resolution is given by ∆ = 80 µm/1024 in the propagation direction.In the free vacuum the corresponding frequency is to be associated to half the Nyquist frequency.32 periods have passed after 10 µm/c and 64 periods after 20 µm/c.The damping and also the barely subluminal phase velocity are, after a long enough propagation time, noticeable even at this wavelength.
The initial peak position has clearly shifted after 20 µm/c.Right: a plane wave with λ = 156.25 nm at different time steps with the same grid settings.In the free vacuum the corresponding frequency is to be associated to the Nyquist frequency.The wave is standing and the damping is strong, resulting in a rapid annihilation.

Domain of process 1
Domain of process 2 From Figure 5 it can be deduced that the waves are well-behaved and quite well-modeled for k • ∆ ≲ 0.5.From Figure 6, it can be inferred that already at half the Nyquist frequency the damping is non-negligible for relevant time scales.The Nyquist frequency in this scenario is given by f Ny = 64 × 10 5 m −1 (λ = 1.5625 × 10 −7 m).A video demonstrating the evolution of a plane wave on the lattice with a wavelength corresponding to half the Nyquist frequency can be found in the Mendeley Data repository [101].The conclusion is that for a proper modeling and to avert damping effects, one is obliged to adapt the grid resolutions to the lowest wavelength such that ∆ ≲ 1/12 λ.While this relation is not a hard limit and can be relaxed in many cases without detriments, it poses a safe rule of thumb for long-time simulations.A sufficiently fine grid resolution is pertinent for a clean analysis of the polarization flipping and harmonic generation effects in the present work.

Solving the ODE and processing the data
For the numerical solution of the nonlinear system of ODEs for ⃗ f in (38) the CVODE solver is used, which is part of the SUNDIALS family of solvers [102][103][104][105].The implicit Adams method (Adams-Moulton formula) in conjunction with a fixed-point iteration is utilized.
The numerical time integration error is controllable with CVODE by setting relative and absolute integrator tolerances per user-defined step.The solver adapts its internal time step sizes according to the system's dynamics, guided by the tolerances.Larger time steps are performed in quiet regions and shorter steps in highly dynamic regions.Errors per step accumulate to a global error.
Integrator error tolerances are set to below 10 −12 for the present work.Since CVODE integrates with high accuracy, the size of a time step can be defined as is convenient and in this work varies in the range of 1 fs to 6 fs, or 0.3 µm/c to 2 µm/c .The computation time to achieve a given numerical integration error threshold is approximately independent of the stencil order, but the correct physical solution is rather approached with a higher stencil order.
Output data are written at each user-defined time step optionally into convenient comma-separated values (CSV) files or a compact and fast binary file.A CSV output file contains six columns for the field components E x , E y , E z , B x , B y , and B z ; each column holding the grid values of the corresponding field component.A binary output file aligns all six field components for every lattice point.The postprocessing of the data for this work is done with the help of Python scripts, Jupyter notebooks employing the SciPy library, and with the help of Mathematica and Paraview [106-109].

Parallelization of the algorithm
For faster computation and scalability the code is parallelized.To this end a Cartesian domain decomposition of the lattice is performed, allowing individual compute cores to process patches of the lattice.The number of these sub-lattices is determined by the user who may subdivide each dimension of the lattice into a number of equally sized parts.The finite differences scheme used to discretize derivatives requires values from neighboring sub-lattices when it is applied at the boundaries of a patch.For this purpose, ghost cells are placed at the boundaries and updated via message passing making use of MPI [110,111], see Figure 7.
Since the sub-lattices form a Cartesian grid, the communication scheme is conveniently implemented in an MPI virtual Cartesian topology.Output to CSV format is written to one file per process (patch), whereas output in binary form is written to one single file per step, making use of MPI IO.The lattice is sliced into sub-lattices in one dimension (blue line and dots).The runtime speedup decreases where one memory domain (socket) is fully occupied with sixteen cores.The code for this benchmark configuration is memory-bound as this typical behavior at the socket saturation shows.The intra-socket scaling is not linear since the memory as a shared resource on the socket does not scale along with additionally used cores.Slicing only in one dimension leads to a communication overhead when the sub-lattices become too narrow.This is remedied in this scenario by an equal slicing in every dimension such that the patches are cubic (red line and dots).The scaling across nodes is then again optimal.Each setting ran twice and is averaged in order to take into account minor runtime variations.More tests prove that cubic patches form the most efficient decomposition.
Since efficient parallelization is paramount for expensive 3D simulations, strong and weak scaling tests proving the basic parallelization capability are conducted.For the tests no output data are generated, yet the MPI IO variant is highly efficient.The employed high-performance computing system contains sixteen cores per memory domain (here a socket) and two sockets per node.The results of the scaling tests are visualized in Figures 8 and 9.
Runtimes and bottlenecks vary strongly with the particular setting for the problem under consideration.Therefore, it is hard to give universal benchmarks and provide general scaling properties.Bottlenecks have been investigated using the Intel ® oneAPI HPC Toolkit [112].For simulation configurations as used in the present work the code is compute-bound in 1D, memory-bound in 2D, and MPI-bound in 3D.
Remarkable speedups in the latter case of full three-dimensional simulations can thus be achieved through a reduction of the communication overhead by lowering the order of the numerical scheme.This is because the required depth of the ghost layers decreases, e.g., from seven to three from order thirteen to order 4, as can be seen from the stencils in Section 5. To further ease the messaging pressure and memory requirements at large scales, additional parallelism by multi-threading with the help of OpenMP [113,114] is used on top of the multi-processing.As mentioned above, the stencil order is on the other hand not decisive for computation speed alone.

Phase velocity in a strong background
A probe plane wave is propagated along the x-axis, through a linearly polarized strong electromagnetic background with field strength A b .The polarization of latter breaks the isotropy of space, giving rise to different refractive indices [29,30,72,115] for a probe polarization orthogonal (+) and parallel (-) to the background polarization.Since Equation ( 72) only takes into account the four-photon interaction contribution, i.e., it neglects all but the first nonlinear term in the weak-field expansion (9), the results are verified turning off six-photon processes in the simulations.Figure 9: Weak scaling test for a simulation such as to obtain 3D results shown in Section 10.While the strong scaling test uses a fixed overall lattice, a weak scaling test is obtained by keeping the patch size constant and thereby increase the size of the overall lattice with the number of patches and compute cores.Hence, the problem size is increased along with the number of parallel workers.For such 3D simulations the code becomes more and more communication-bound.Runtime variations and scaling issues are mainly attributable to data transfer via the node interconnect and resource contention thereon.The weak scaling can be considered satisfying with the only large increase in runtime observed at the last run.This one occupies 50 of a total of 153 nodes on the cluster.These simulations ran only once each, in order to save resources.

Lattice Points 1000
Background ⃗ A (0, 3, 0) × 10 −6 E cr up to (0,0.9,0)The resulting phase velocity change v nli from the vacuum speed of light, given by can be extracted with the help of a Fourier analysis of a time-propagated wave.When the passed time of propagation is chosen to be an integer multiple of λ, it is obtained with l ∈ N [83] The nonlinear phase velocity contribution can then be extracted from the phase after a spatial Fourier transformation evaluated at λ −1 .This results in The spatial Fourier transformation in Equation ( 75) can be replaced by a Fast Fourier Transformation in the analysis.
To analyze the phase velocity variation numerically, the background field strength is varied.In a second step the relative polarization of the waves is changed from parallel to orthogonal.The configurations are given in Table 1.The total simulation time is chosen to be 200 µm/c, conveniently divided into 100 steps of 2 µm -the chosen wavelength of the probe wave.
It can be seen that the nonlinear interactions give note to themselves in a reduction of the phase velocity of the simulated waves in Figure 10.For sufficiently large background field strengths the numerical values are in very good agreement with the analytical predictions.

Polarization flipping -vacuum birefringence
Polarization, or helicity, flipping of a fraction of photons in a probe pulse propagating through a strong background pump pulse is a result of vacuum birefringence.The origin of the effect is again the breaking of the isotropy of space by the polarization of the strong background.The refractive indices from above,  for parallel and 4.5% for orthogonal relative polarization of the probe.This difference originates in the probe wave not being a perfect "probe" as in the idealized theoretical scenario and thus is contributing with its own polarization of A p = 3×10 −6 E cr .The values for background field strengths larger than 10 −4 E cr have a mean absolute percentage error of 1.8% for parallel relative polarization and 1.2% for orthogonal relative polarization of the probe.
Figure 11: Qualitative depiction of the electric fields in a coaxial probe-pump experiment for the measurement of vacuum birefringence.The probe (green) traverses the counter-propagating pump (blue), experiencing a polarization rotation as a result of different refractive indices.The originally linearly polarized probe is afterwards marginally elliptical (turquoise).The effect is depicted greatly exaggerated for visibility.This effect is attributable the fact that the isotropy of space for the charged particleantiparticle fluctuations in the vacuum is broken by the polarization of the strong background pulse.The coupling of these fluctuating particles in turn to the probe pulse results in different refractive indices for the polarization modes of the probe, as given in Equation ( 76).
generate a difference in optical path length for the probe pulse polarization components parallel and orthogonal to the pump polarization, which results in birefringence.On a microscopic level, a portion of the probe pulse's quanta flip their polarization.Macroscopically, the overall polarization experiences a tiny rotation.
A typical probe-pump scenario devised for the observation of helicity flips is sketched in Figure 11.A probe pulse propagates through a strong low-frequency pump field in which spatial isotropy is broken.While propagating through a pump field a fraction of probe pulse photons flips their polarization by 90°which results in a tiny ellipticity of an initially linearly polarized probe pulse.A corresponding simulation configuration showing one polarization direction is shown in Figure 12.
In order to backtest the numerical solver, firstly, in Subsection 8.1, parametric checks of the flipping probability as derived in [37] are performed.The settings given in Table 2 are chosen to this end with Gaussian pulses given by Equation ( 78) below.With that result being verified, secondly, in Subsection 8.2, the parametric scaling properties are made use of to extrapolate results to wavelenghts in the x-ray regime.This is compared to a calculation for a realistic polarization flipping scenario calculated in [38].The parameters for that setup are listed in Table 3.
In both cases the normalized vectors parallel and orthogonal to the initial probe polarization are Note that in actual simulations (parameters in Table 2) the probe field strength and wavelength are significantly smaller.Adjustments are made here for better visibility.
given by For the simulations 1D Gaussian pulses are used in the form where the vector ⃗ A comprises amplitude and polarization, ⃗ x 0 is the center of the pulse, τ its width, λ the wavelength, and k the unit propagation direction vector.
Pulses are implemented in space without explicit time dependence.With spatial derivatives via the finite difference scheme the ODE in time (38) is formulated, which is solved by CVODE for the time evolution.For convenience, the normalized vector k indicating the propagation direction is stated in the parameter tables.
Since analytical estimates for polarization flipping such as in [38] make use of a photon picture, while the numerical simulations presented in the present paper propagate coherent modes, the mapping is used.The energies in the respective polarization directions are proportional to the electric field strength projections squared, All other factors in (79) cancel out.It has to be mentioned that Equation (79) implies that the frequencies of the signal photons equal that of the probe pulse.However, as to be shown in Section 9, the nonlinear interaction results in a small fraction of photons with altered frequency.

Vacuum birefringence -parametrical dependencies
For a probe pulse coaxially counter-propagating to a plane wave background field, the polarization flipping probability, taking into account again only up to four-photon interactions, in the low-energy approximation, is given by [37] where σ is the initial angle between the probe and pump polarizations, λ p the wavelength of the probe pulse, and A b the amplitude of the background pulse.The propagation direction of the probe is assumed to be perpendicular to the pump polarization and the probe field strength to be negligible compared to the pump.The probability directly translates to the flip ratio, Formula ( 81) yields all the parametric dependencies for the probability of polarization flips and indirectly excludes other parameters.There is a strong dependence on the optical path of the pump pulse  81).The probe wavelength, the pump field strength, and their relative polarizations are varied to obtain the parametric scaling results of Figure 14. and on the probe wavelength.Notably, the ratio is independent of the shapes of the pulses.Limitations of the above formula for focused background pulses are discussed in the following subsection.In 1D simulations the background can be modeled as a Gaussian pulse.To investigate the scaling properties of the numerical solver the settings in Table 2 are used, where only those parameters affecting ( 81) are actually relevant.A time-resolved flipping process for those parameters is depicted in Figure 13.

Grid
The results of the parametric scaling tests are visualized in Figure 14.There is perfect agreement between the 1D simulation results and formula (81).With these scaling properties being verified in the algorithm, in the following subsection the analytical result in case (a) of [38], where a small probe pulse traverses a strong pump field, is compared to an extrapolation of simulation results.
Neglecting the signals for σ = 0, π/2, where P flip = 0 analytically, that cannot be respected in a relative error calculation with the true values as baseline, the mean absolute percentage errors for each scaling test are below 0.1%.

Vacuum birefringence -extrapolation to an analytical value in the x-ray regime
Simulations of birefringence effects are computationally expensive in higher dimensions for the small wavelengths and field strengths targeted in experiments.Making use of the scaling properties in Equation ( 81), the phenomenon of birefringence can still be predicted for the parameters accessible in planned near future experiments by simulating numerically feasible, quasi-1D setups with consecutive extrapolation.
To this end this approach is used to reproduce an analytical result for a coaxial probe-pump setup with Gaussian laser pulses in a realistic scenario by considering case (a) of [38].In this case the radius of the probe pulse is taken to be much smaller than the waist of the pump beam, such that the probe does not sense the transverse structure of the pump.This scenario thus amounts to a 1D case.Settings adapted to the scenario described in case (a) in [38] are given in Table 3.The calculation in [38] is performed in the vacuum emission picture.Some parameters devised for experimental verification impair the numerical approach.First, the pump field strength is too low to extract the flipping process from the numerical noise.To combat this, the field strength scaling properties can be used to simulate with a larger background amplitude.Second, the probe wavelength of λ p = 96 pm is in the x-ray regime to amplify the effect, c.f. Equation (81), and is therefore problematically small for modeling on a discrete grid.An extremely fine grid would be necessary to model that pulse.To evade computations that expensive, the wavelength scaling properties can be made use of.The resulting extrapolation is thus a combination of two scaling methods in the way shown in the bottom right of Figure 14 and described in that caption.
Figure 16 shows that the extrapolated flipping ratio of 2.72 × 10 −12 is almost twice as high as the flipping ratio of 1.39 × 10 −12 obtained in [38].This is attributable to the neglect of a longitudinally localizing term with the Rayleigh length in the pulse form in Equation ( 78), c.f. the higher-dimensional Gaussian pulse in Equation ( 100).This yields a further suppression of about a factor of two.The value obtained at the corresponding probe frequency via formula ( 81) is 2.73 × 10 −12 and agrees with  2 and simultaneously an extrapolation to frequencies in the x-ray regime is performed.This last way of combining scaling properties is useful in order to extrapolate to relevant probe frequency regimes while keeping the numerical accuracy high and the computational load low.
Table 3: Parameters for probe and pump beam adapted to [38].The pump field strength is obtained as the square root of the ratio of intensity to critical intensity.The pump pulse duration is 30 fs (2τ in the 1/e 2 criterion).3.One time step corresponds to 1.5 fs.The adaptations are: The pump field strength is magnified by a factor of 100 to 34 ×10 −3 E cr in order to reduce numerical noise.The probe wavelength is enlarged to a computationally acceptable value of 25 nm.The distance between the vertical red lines is the total interaction time t I = 24 fs.The horizontal red line denotes the asymptotic flip ratio.

N /N Simulated Analytical (VEP) Extrapolation
Figure 16: Extrapolation of polarization flipping ratios to the x-ray regime and comparison to a result obtained via the vacuum emission picture in [38] (red dot).The blue line is an extrapolation of simulation results (blue dots) with various probe wavelengths.
Table 4: Initial settings to observe harmonic generation in 1D simulations, see Figure 17.

Grid Length 300 µm
Lattice Points 4000 Figure 17: Visualization of the pulse configuration to detect higher harmonics with a zero-frequency background pulse, see Table 4.
the simulations.Accordingly, in future higher-dimensional simulations a reconciliation of the results should be achieved.

Harmonic generation
To further crosscheck the solver, the prominent probe-pump scenario for the detection of nonlinear vacuum signatures shown in Figure 17 is considered again with two head-on colliding pulses, a strong background pulse and a weaker probe pulse.For this analysis, the former pulse is assumed to have zero frequency.The initial settings are listed in Table 4.
Approximate analytical results for this scenario are derived in [66,68].The effective vertices for four-and six-photon scattering in Figure 18 (a) can produce outgoing photons with higher frequency by photon merging.For example, in Figure 18 (b) two probe photons and a zero-frequency background photon merge into an outgoing photon with frequency 2ω p .The possible contributions of two-wave scattering that result from the first orders of the weak-field expansion are listed in Figure 19.
In the present case of ω b = 0, there are for four-photon processes • the scattering of a background and a probe photon contributing with one photon to the zeroth harmonic (ω r = 0), also called dc component, and with one to the first harmonic (ω r = ω p ), also called fundamental harmonic;  • two background photons and one probe photon merging to produce a photon of the fundamental harmonic (ω r = ω p ); and • two probe photons and one background photon merging to produce a photon of the second harmonic (ω r = 2ω p ).
For six-photon processes it is obtained • the sheer scattering of background and probe photons contributing to the dc component and the fundamental harmonic; • two background and two probe photons merging and producing one photon contributing to the second harmonic and one contributing to the dc component; • two background and two probe photons merging and producing two photons contributing to the fundamental harmonic; • two background photons and three probe photons merging and producing a photon contributing to the third harmonic (ω r = 3ω p ); and • the merging of three background and two probe photons producing a photon contributing to the second harmonic.
A visualization of the contributions at selected points in time is provided with Figure 20.These simulations are strictly 1D and the use of plane waves leads to strong constraints.It can be seen that the asymptotic contribution to the second harmonic is only attributable to six-photon processes.That is because wave-mixing -resulting pulses that are combined of photons of both pulses -is not allowed asymptotically.The reason behind this is energy conservation, since for coaxial pulses and a photon resulting of wave-mixing it is found where n p and n b are the numbers of the contributing photons and k is the unit propagation direction vector of the probe pulse.These states are thus only visible in the overlap position.The six-photon process, on the other hand, can produce second harmonics without wave-mixing, see the second point for six-photon processes above.
As discussed in the context of birefringence in Section 8, the 1D case corresponds to a simplified handling of the experimentally relevant scenario of counter-propagating pulses.where the pulses are directly overlapping, at the final state they have separated again -the asymptotic field is left.It can be deduced that the third harmonic and the asymptotic part of the second harmonic are solely ascribable to six-photon processes.
For the highest generated harmonic, the short-lived third harmonic, the rule-of-thumb resolution limit for the grid defined at the end of Section 5 is slightly exceeded.This comes without noticeable accuracy problems as is shown in the next subsections.

Harmonic generation -analytical results
Analytical methods in [66,68] contain a derivation of iterative solutions to the nonlinear equations of motion for zero-frequency backgrounds.With the probe (p) and background (b) pulses as timedependent 1D Gaussian pulses with the parameters of Table 4 the pulses are obtained to be, compared to (78), with k µ j x jµ = ω j t − ⃗ k j ⃗ x j .The shifted coordinates of probe and background field read Writing a combined electric field as the inhomogeneous wave equation in 1D can be written as The source term T[ ⃗ E] comprises the nonlinear Heisenberg-Euler interactions.Decomposing the electric field into where ⃗ E (0) solves the free vacuum wave equation (∂ 2 t /c 2 −∂ 2 x ) ⃗ E (0) = 0, an iterative procedure is obtained [66,68] With the polarization for both fields given by ⃗ ϵ = (0, 1, 0) and defining the shorthands it is then obtained for the solution to the nonlinear wave equation ( 89) at the first iterative order for [66,68] • the dc component: and the asymptotic field • the fundamental harmonic: and the asymptotic field • the second harmonic: the overlap field and the asymptotic field • the third harmonic: and the asymptotic field ⃗ E (1)  3,a = 0 .
Using the values from Table 4 the solid lines in Figure 22 are obtained.Those are employed to scrutinize the correctness of the simulation results.
A Mathematica [108] analysis of the analytical results for harmonic generation can be found in [116].Animations of the arising and evolution of the various harmonics are provided in the Mendeley Data repository [101] (thumbnails and description in Figure 21).

Harmonic generation -simulation results
To put the focus on the nonlinear contributions to the generation of harmonics, each setting is simulated three times with varying combinations of interactions included: once including only the nonlinear effects of four-and six-photon contributions; once including only six-photon processes; once excluding all nonlinear interactions, i.e., keeping only the linear vacuum.By subtracting the dynamics in the linear vacuum from the full dynamics, the higher order processes of the weak-field expansion are extracted.Furthermore, he simulations of six-photon processes permit to isolate their sole contribution, as is shown in Figure 20.It is a fruitful feature of the simulation code that the contributions of fourand six-photon diagrams can be turned on and off to make them separately visible.
In order to extract the amplitudes of the various arising harmonics and their time evolution, their respective frequencies have to be filtered in Fourier space and then transformed back to position space [83].The amplitude is then obtained as the maximum norm.Its time evolution can be observed, in this case with a chosen resolution of a time step of 1 µm/c.The results for the zeroth harmonic, first harmonic, second, and third harmonic are displayed in Figure 22.
As can be seen in Figure 22, there is good agreement between the analytical approximation and the simulation results.Small systematic errors are unavoidable as a consequence of the back and forth Fourier transformation and slicing of frequency ranges.The mean absolute percentage errors of the simulation results are calculated to be less than 1% in the regions where the amplitudes are non-vanishing.
Asymptotic states are constrained by energy-momentum conservation, while the overlap state has a richer spectrum.The overlap spectrum becomes even more pronounced and versatile when the pulses collide at a non-zero angle in higher dimensions.First such results are demonstrated in Section 10, where also the zero-frequency background restriction is relaxed.Signals, which are degenerate in the case of a non-zero-frequency pulse, split up in that case.
Simulations in higher dimensions provide a powerful means to analyze varying collision configurations.Situations that pose no further difficulty to the numerical code are considerably hard to cope with analytically.Employing simulations of the solver it is possible to track harmonic frequencies in time and space for scenarios of arbitrary pulse parameters -with the only restrictions posed by the applicability of the Heisenberg-Euler weak-field expansion and computational feasibility.

Higher-dimensional simulations
For simulations in 2D a special adaptation of 3D Gaussian pulses is used to model the diffusion behavior.The pulse is assumed to propagate along the z-axis.The widening of the beam with respect to the longitudinal coordinate z is given by the waist with w 0 the waist of the beam at position z 0 , where the amplitude is 1/e of the initial value.The crosssectional area in 1D is a point, in 2D is a line, and in 3D is an area.The field intensity scales with w 0 /w(z) and with the surface area ∼ z 2 in 3D.Lateral dispersion has to be taken into account for the 2D Gaussian pulses, where the surface scales as ∼ z.Hence, the factor w 0 /w(z) appearing as prefactor in the Gaussian pulses gets a square root in the lower-dimensional case.The pulse can thus be written as It is further defined • the parameter A determining the peak pulse amplitude and the polarization ⃗ ϵ ; • the distance to the propagation axis (here taken to be z) r = x 2 + y 2 ; • the wavenumber k = 2π/λ ; • the pulse width in z-direction τ z (pulse duration) and the envelope center z τ ; and • the Rayleigh length z R = πw 2 0 /λ as the longitudinal distance from z 0 at which the waist has increased by a factor of 2 , which is contained in • the Gouy phase ζ(z) = arctan(z/z R ) , and The settings used for 2D simulations, which are confined to the x-y-plane for propagation, are listed in Table 5.
For coaxial pulse collisions there is a correspondence to the 1D scenario with respect to the generated harmonics.Comparing a head-on collision in 1D (Figure 20) and 2D (Figure 23), a similar frequency spectrum is found at the different states in time.The contributions of the different orders in the weak-field expansion are demonstrated again (Figure 24).The reasoning is the same as in Section 9.The main difference stems from the fact that the two conceptually equal pulses have the same non-zero frequency ω p .There is still a degeneracy of signals present by virtue of the equal frequencies.
A feature which is not present in 1D simulations is a lateral broadening of the pulses in frequency space, which arises during the interaction and remains.The presence of a lateral beam profile of the pulses is a prerequisite to invoke this outcome.This results in outgoing signal photons with transverse momentum components, effectively yielding a diffraction effect [35].Diffraction spreading opens up the opportunity to detect signal photons off the beam axis with a background free measurement.The scattering of polarization-flipped signal photons outside the forward cone of a probe beam may thus constitute an essential key ingredient for the detection of vacuum birefringence [38].

Lattice Points 1024×1024
Pulse 1 ⃗ ϵ (0,0,1)  Correspondingly, from the position space point of view, the transversal momenta might imply a slight focusing of the pulses.This can be explained with the lensing effect of a power pulse which creates a refractive index influencing the propagation speed and direction, as detailed in Section 7.With a lower refractive index at the outer waist regions of at least one of the pulses, light passing through the strong-field zone experiences a phase velocity change comparable to the one in a convex lens.Focusing of light by light is an interesting topic to be further investigated with the help of adequate simulations with tailored pulse parameters.
Going beyond coaxial pulses by varying the collision angle and thereby lifting the degeneracy of frequencies of the harmonics, geometry effects with rich spectra in 2D simulations can be observed, see Figures 25-28 for perpendicularly propagating and colliding pulses as well as for a collision angle of 135°.Most signals vanish again in the asymptotic state.These off-axis contributions occur due to the field spatio-temporal inhomogeneities [96].
The asymptotic harmonics that can be seen in the right frames of the frequency plots are on account of the self-interaction of the 2D Gaussian pulses [65,96].Time-resolving the processes reveals that these harmonics arise immediately as the dynamics begin, directly after the initial configuration shown at the lower left of the figure, and thus already before the pulses overlap.This can be seen in the corresponding simulation videos in the Mendeley Data repository [101].Both the four-and six-photon processes contribute to the asymptotic signals with | ⃗ k| = ω p /c and | ⃗ k| = 3ω p /c. Six-photon processes contribute to all harmonics, in the overlap as well as in the asymptotic state.In the overlap states various merging processes become directly visible.A more precise analysis shows that the four-photon spectrum in the overlap state is even richer.Some signals, however, are of the order of accumulating numerical errors.Such errors appear, e.g., at the corners of the frequency plots for four-photon and six-photon processes, irrespective of the pulse alignment.
The signals of asymptotic harmonics generated by four-photon interactions are aligned on the initial propagation axes, while six-photon processes seem to generate a tiny off-axis twist.Note that the initial symmetry of the two-pulse system is thereby conserved.
Giving the pulses different polarization directions, the harmonics can be tagged.The frequency space of the simulations visualized in Figures 29 and 30 show again collision angles of 90°and 135°, but the pulse propagating from the left is now polarized along the E y -direction.Hence, the shown polarization components can be identified with one of the two pulses, and thus also the harmonics.
Ultimately, in order to take into account all geometry effects, simulations have to be conducted in full three spatial dimensions.A demonstration of configurations similar to those in 2D as shown in Figure 25, in 3D yields results visualized in Figure 31.Pulses with frequencies differing by a factor two colliding collinearly produce the results of Figure 32.By virtue of the differing frequencies expressly rich harmonics spectra can be observed.The corresponding 1D case with a breakdown of the frequencies is demonstrated in the examples provided in the code repository [117].Colliding the two pulses at an angle of 135°generates the harmonics shown in Figure 33.It can be seen that the weakest generated signals have about the magnitude of numerical artifacts.Figure 33: Rich harmonics spectrum in 3D.This spectrum is generated by the two Gaussian pulses of Figure 32 colliding at an angle of 135°.There are numerical remnants present on the coordinate axes, where the zero values in frequency space are located and small numerical errors accumulate.These can also be observed in Figure 31.

Conclusion and outlook
Within the limits of the Heisenberg-Euler weak-field approximation the numerical approach represents an efficient solver for the complete dynamical response of the nonlinear vacuum for complicated pulse setups in up to three dimensions.
Every simulation of the Heisenberg-Euler model incorporates the complete nonlinear physics in the weak-field approximation, whereas analytical calculations normally have to concentrate on single, isolated effects, leaving others aside, and hence often miss the complete picture of the interaction.By taking into account the whole dynamics of the nonlinear vacuum, the solver captures in particular back-reactions to the radiation fields.Withal, simulations permit to describe the temporal evolution of nonlinear vacuum processes, which is beyond the reach of many analytical approaches.Simulations permit to time-resolve the nonlinear vacuum processes, which is beyond the reach of many analytical approaches.
The validity of the presented numerical scheme for solving the modified Maxwell equations in the Heisenberg-Euler weak-field expansion relies on two basic assumptions: I) field strengths below the critical values E cr and B cr , and II) wavelengths larger than the Compton length of the electron.The laser pulses are considered on a macroscopic level implying that the individual photons the pulses consist of are not resolved.
A validated solver of the Heisenberg-Euler dynamics is very useful for strong-field QED research going on at present.The current C++ implementation of the numerical scheme discussed in this paper allows distinct simulations of the linear Maxwell vacuum, the four-photon processes, the six-photon processes, and any combination of the latter processes.
Of paramount significance is the dispersion relation, lying at the heart of the algorithm, which ensures stability throughout the frequency spectrum and moreover creates an imaginary part that annihilates nonphysical modes.Furthermore, the linear vacuum-like behavior of the dispersion relation for a large frequency range is an essential ingredient.
A good agreement with analytical results is achieved by making use of high discretization orders of the numerical scheme and high expansion orders of the effective Heisenberg-Euler model.The computational cost scales strongly with the number of lattice points but weakly with the discretization orders of the scheme and the expansion orders of the Heisenberg-Euler model.The impact of the discretization order on the computational cost, however, does become relevant in 3D with increasing MPI communication required in order to transfer spatial derivative data.Making use of high discretization orders and their favorable dispersion relations facilitates the use of comparably small lattices to accurately model the involved waves.Nonetheless, for high-frequency pulses the required grid resolutions can become unfeasible.
There are facilities constructed with the purpose to detect vacuum birefringence [118,119].The goal is to provide a tool for experimentalists to second their setups with simulation data.While vacuum birefringence effects in 2D are simulated in a parallel project, even the employed cluster computing system does not provide enough computing power to for such simulations in 3D as a consequence of the small probe wavelengths required to enhance the effect, c.f. Equation (81).
Ideas are being developed in order to overcome the obstacle of extremely large 3D grids.One promising, ongoing, and important project is on multi-scale simulation capability.This research is intended to pave the way for a dynamical grid, adapting its resolution regionally on the fly and on demand in order to reduce the computational load, while at the same time maximizing the accuracy in important spatial regions.An adaptive grid is particularly suited for the prominent low-frequency pump pulse and high-frequency probe pulse setup to detect vacuum nonlinearities, where the computational load can in principle be dramatically reduced.

4 πFigure 3 :
Figure3: Dispersion relations of the numerical scheme for minimally biased finite differences at order four (left) and order thirteen (right).The black lines represent the real vacuum dispersion relation with c = 1.∆ is the grid spacing, the physical distance between lattice points.For better visibility of the symmetric plots only the values for k ≥ 0 are shown.
Effect of overshooting the Nyquist frequency.Left: two cosine functions that, when periodically evaluated at multiples of the distance ∆ = 0.4π, have the same values.Right: simulations to verify the behavior on the grid.On a 1D line with point spacing ∆ = 0.1µm two waves with wavelengths 5∆ and 5/4∆, respectively, are simulated, corresponding to the analytical scenario on the left.It can be seen that the discretization makes no difference between the two.The simulated waves are shown in the initial state (top right) and after they have propagated the distance of 90 periods (bottom right).The damping effect initiated by the imaginary part of ω in the dispersion relation is the same for both waves.

Figure 5 :Figure 6 :
Figure 5: Numerical tests of the dispersion for low frequencies.Left: numerical results of dispersion relation tests.Simulationresults (blue dots) are in agreement with the free vacuum (black line) for wavelengths from 20 µm to 1 µm.Right: the plane wave corresponding to the rightmost point of the left figure with λ = 1 µm at different time steps.The wave is not that accurately modeled anymore but still does not lose energy in the given amount of time -the amplitude stays the same.Even after 32 periods the wave is perfectly overlapping with its initial state.

Figure 7 :
Figure 7: Ghost-cell exchange.The blue boxes represent sub-lattices (patches), batches of the interdependent grid data, each being processed on single cores.The processes (compute cores) exchange their boundary cells' values with the neighboring processes.The boundary regions, whose required depths depend on the order of the finite differences scheme, are indicated in red.This exchange is performed sequentially for each dimension.Here it is shown for one boundary region between two processes.

Figure 8 :
Figure8: Strong scaling test for a relatively cheap simulation in 3D.The overall lattice keeps its size and is split into smaller and smaller patches to be processed by single cores.The runtimes nearly halve for a doubling of cores for up to eight cores.The lattice is sliced into sub-lattices in one dimension (blue line and dots).The runtime speedup decreases where one memory domain (socket) is fully occupied with sixteen cores.The code for this benchmark configuration is memory-bound as this typical behavior at the socket saturation shows.The intra-socket scaling is not linear since the memory as a shared resource on the socket does not scale along with additionally used cores.Slicing only in one dimension leads to a communication overhead when the sub-lattices become too narrow.This is remedied in this scenario by an equal slicing in every dimension such that the patches are cubic (red line and dots).The scaling across nodes is then again optimal.Each setting ran twice and is averaged in order to take into account minor runtime variations.More tests prove that cubic patches form the most efficient decomposition.

Figure 10 :
Figure10: Nonlinear contribution to the phase velocity slow-down in a background with varying field strength.The simulation results converge to a value of v nli ≈ 2 × 10 −12 for small background field strengths.This is the phase velocity reduction caused by numerical errors, which are getting larger than the physical effect for background field strengths A b < 10 −4 E cr .Deviations from the analytical expectation are higher for lower background field strengths.The error at A b = 3 × 10 −4 E cr is still 6.1% for parallel and 4.5% for orthogonal relative polarization of the probe.This difference originates in the probe wave not being a perfect "probe" as in the idealized theoretical scenario and thus is contributing with its own polarization of A p = 3×10 −6 E cr .The values for background field strengths larger than 10 −4 E cr have a mean absolute percentage error of 1.8% for parallel relative polarization and 1.2% for orthogonal relative polarization of the probe.

Figure 12 :
Figure 12: Sketch of the pulse configuration for the simulation of polarization flipping as shown in Figure 11.On the left is the weaker probe pulse which propagates to the right; on the right is the strong pump pulse propagating leftwards.Note that in actual simulations (parameters in Table2) the probe field strength and wavelength are significantly smaller.Adjustments are made here for better visibility.

]Figure 13 :
Figure 13: Time evolution of the polarization flipping ratio for the parameters presented in Table2.A simulation time of 30 µm divided into 100 steps, so 1 fs per step is used.The settings used here provide an interaction time of 25 fs, indicated by the red vertical lines.The red horizontal line corresponds to the asymptotic relative flip ratio.

2 σ 5 . × 10 - 10 1. × 10 - 9 1. 5 × 10 - 9 2. × 10 - 9 2. 5 × 10 - 9 N ⟂ /N 4 . 9 5. × 10 - 9 1. × 10 - 8 NFigure 14 :
Figure 14: Scaling of the polarization flipping probability with variations of the parameters given in Table2.The solid lines are the analytical curves obtained with the help of(81).The red dots are simulation results.Top left: varying the background field strength.Top right: varying the relative polarization angle of probe and pump.Bottom left: varying the probe wavelength.Bottom right: combining the scaling of the pump field strength and probe wavelength; A b is scaled down by a factor of 100 compared to the parameter in Table2and simultaneously an extrapolation to frequencies in the x-ray regime is performed.This last way of combining scaling properties is useful in order to extrapolate to relevant probe frequency regimes while keeping the numerical accuracy high and the computational load low.

5 N /N [10 9 ]Figure 15 :
Figure 15: Time evolution of the polarization flipping probability for an adaptation of the parameters presented in Table3.One time step corresponds to 1.5 fs.The adaptations are: The pump field strength is magnified by a factor of 100 to 34 ×10 −3 E cr in order to reduce numerical noise.The probe wavelength is enlarged to a computationally acceptable value of 25 nm.The distance between the vertical red lines is the total interaction time t I = 24 fs.The horizontal red line denotes the asymptotic flip ratio.
(a) Effective vertices for four-and six-photon scattering.(b) Example of high-harmonic generation with a zero-frequency background.

Figure 18 :Figure 19 :
Figure 18: Feynman diagrams for harmonic generation with a probe (subscript p) and background (subscript b) wave.In (b) there is an implicit time axis from left to right.

Figure 20 :
Figure20: Log-scale making the higher harmonics visible in frequency space.Top: a full simulation in the linear vacuum supplemented by four-and six-photon nonlinear interactions.Bottom: after subtraction of the linear vacuum.The initial state contains only the main signals of the pulses with no nonlinear interaction yet present.The overlap state denotes the time step where the pulses are directly overlapping, at the final state they have separated again -the asymptotic field is left.It can be deduced that the third harmonic and the asymptotic part of the second harmonic are solely ascribable to six-photon processes.

Figure 21 :Figure 22 :
Figure 21: Thumbnail of animations available in the Mendeley Data repository[101] showing the time evolution of the various harmonics in 1D as a consequence of a pulse collision with the specifications of Table4with nonlinear four-photon and six-photon interactions.

Figure 23 :
Figure 23: Coaxially colliding pulses with the same polarization.The left plots show the initial state, those in the middle the overlap state, and the right ones the final state.Top: position space.Bottom: frequency space.

Figure 24 :
Figure 24: Frequency space of coaxially colliding pulses with the same polarization.The left plots show the initial state, those in the middle the overlap state, and the right ones the final state.Top: only four-photon diagrams included.Bottom: only six-photon diagrams included.

Figure 25 :
Figure 25: Perpendicularly colliding pulses with equal polarization.The left plots show the initial state, those in the middle the overlap state, and the right ones the final state.Top: position space.Bottom: frequency space.

Figure 26 :
Figure 26: Frequency space of perpendicularly colliding pulses with same polarization.The left plots show the initial state, those in the middle the overlap state, and the right ones the final state.Top: only four-photon diagrams included.Bottom: only six-photon diagrams included.

Figure 27 :
Figure 27: Pulses with the same polarization colliding at an angle of 135°.The left plots show the initial state, those in the middle the overlap state, and the right ones the final state.Top: position space.Bottom: frequency space.

Figure 28 :
Figure 28: Frequency space of pulses with the same polarization colliding at an angle of 135°.The left plots show the initial state, those in the middle the overlap state, and the right ones the final state.Top: only four-photon diagrams included.Bottom: only six-photon diagrams included.

Figure 29 :
Figure 29: Frequency space of perpendicularly colliding pulses with orthogonal relative polarization.The left plots show the initial state, those in the middle the overlap state, and the right ones the final state.Top: the E z component is shown.Bottom: the B z component is shown.

Figure 30 :
Figure 30: Frequency space of pulses colliding at an angle of 135°with orthogonal relative polarization.The left plots show the initial state, those in the middle the overlap state, and the right ones the final state.Top: the E z component is shown.Bottom: the B z component is shown.

Figure 31 :
Figure 31: 3D simulation of two perpendicularly colliding Gaussian pulses.The visualizations on the left show the overlap state, the visualizations on right the final state.Top: position space.Bottom: frequency space.

Figure 32 :
Figure 32: 3D simulation of two coaxially colliding Gaussian pulses with different frequencies.Left: initial pulse configuration.Right: harmonics spectrum at the overlap position.

Table 1 :
Settings for phase velocity variation tests.

Table 2 :
Parameters for probe and pump beams chosen to test the parametric dependencies in Equation (

Table 5 :
Settings for 2D simulations with two conceptually equal Gaussian pulses.The wavelength is obtained via λ = πw 2 0 /z R .The Rayleigh length and waist are chosen such that the wavelength equals one micrometer.