2 Modified electromagnetic wave propagation

Photons can scatter off one another by a mutual coupling to the virtual fermionic states of the vacuum (see figure 1). When the centre-of-mass energy is much less than the electron rest energy $mc^{2}$ , the scattering can be well described by an effective description in which the fermion dynamics have been integrated out (it is exact for constant fields). The corresponding Lagrangian was derived by Heisenberg and Euler (Heisenberg & Euler Reference Heisenberg and Euler1936) and takes the form:

(2.1) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{HE}=-\frac{m^{4}}{8{\rm\pi}^{2}}\int _{0}^{\infty }\text{d}s\frac{\text{e}^{-s}}{s^{3}}\left[s^{2}ab\,\text{cot}\,as\,\text{coth}\,bs-1+\frac{s^{2}}{3}(a^{2}-b^{2})\right], & & \displaystyle\end{eqnarray}$$

(unless they occur explicitly, $\hslash =c=1$ ), where the secular invariants $a$ and $b$ are given by:

(2.2a,b ) $$\begin{eqnarray}\displaystyle a=\left[\sqrt{{\mathcal{F}}^{2}+{\mathcal{G}}^{2}}+{\mathcal{F}}\right]^{1/2};\quad b=\left[\sqrt{{\mathcal{F}}^{2}+{\mathcal{G}}^{2}}-{\mathcal{F}}\right]^{1/2}, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{F}}=(E^{2}-B^{2})/2$ and ${\mathcal{G}}=\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}$ and to simplify expressions, we express electric and magnetic fields in units of the critical field $E_{cr}=m^{2}c^{3}/e\hslash =1.3\times 10^{16}~\text{V\,cm}^{-1}$ . The critical or ‘Sauter’ field $E_{cr}$ gives the field strength scale (field strengths occur normalised by this quantity) for vacuum polarisation effects. Although $E_{cr}$ is several orders higher than the current record for the strongest field produced in a laser (Yanovsky et al. Reference Yanovsky2008), the field scattered due to vacuum effects also depends on the phase length of the probe, which can compensate for the small ratio of $E$ in laser experiments. This aspect is exploited in the reported work.

Applying the Euler–Lagrange equations to ${\mathcal{L}}={\mathcal{L}}_{M}+{\mathcal{L}}_{HE}$ where ${\mathcal{L}}_{M}=m^{4}(E^{2}-B^{2})/8{\rm\pi}{\it\alpha}$ is the purely classical term, modifies the inhomogeneous Maxwell equations:

(2.3a,b ) $$\begin{eqnarray}\displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{E}={\it\rho}_{vac}[\boldsymbol{E},\boldsymbol{B}];\quad \boldsymbol{{\rm\nabla}}\wedge \boldsymbol{B}=\boldsymbol{J}_{vac}[\boldsymbol{E},\boldsymbol{B}]+\frac{\partial \boldsymbol{E}}{\partial t}, & & \displaystyle\end{eqnarray}$$

with the vacuum charge density ${\it\rho}_{vac}$ and vacuum current $\boldsymbol{J}_{vac}$ :

(2.4a,b ) $$\begin{eqnarray}\displaystyle {\it\rho}_{vac}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{P};\quad \boldsymbol{J}_{vac}=\boldsymbol{{\rm\nabla}}\wedge \boldsymbol{M}+\partial _{t}\boldsymbol{P}, & & \displaystyle\end{eqnarray}$$

where the vacuum polarisation $\boldsymbol{P}=\partial {\mathcal{L}}_{HE}/\partial \boldsymbol{E}$ and vacuum magnetisation $\boldsymbol{M}=\partial {\mathcal{L}}_{HE}/\partial \boldsymbol{B}$ . Equations (2.3) demonstrate the neoclassical approach (modifying classical equations of motion to subsume the effect of quantum phenomena on smaller length scales (Delphenich Reference Delphenich2006)). This approach is valid as photon energies are below the (effective) threshold for pair creation, clearly a quantum effect beyond the classical framework. It has the advantage that established methods from classical electrodynamics can be employed to solve for the resultant field. Standard approaches calculating observables for laser-based experiments neglect the dependency of vacuum polarisation on the scattered field itself, and a highlight of the project was to take this into account to all orders in the corresponding Dyson-like iterated series.

Figure 1. For weak fields, the full interaction with the charged electron–positron states of the vacuum (left) can be expanded in terms of four-, six- and eight-wave mixing diagrams (right).

For application to terrestrial experiments, it suffices to take a weak field ( $E,B\ll 1$ ) expansion of the Heisenberg–Euler Lagrangian equation (2.1) (depicted in figure 1):

(2.5) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{HE} & = & \displaystyle \frac{m^{4}}{{\it\alpha}}\mathop{\sum }_{i=1}^{\infty }{\mathcal{L}}_{i},\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{1} & = & \displaystyle \frac{{\it\mu}_{1}}{4{\rm\pi}}[(E^{2}-B^{2})^{2}+7(\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B})^{2}],\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{2} & = & \displaystyle \frac{{\it\mu}_{2}}{4{\rm\pi}}(E^{2}-B^{2})[2(E^{2}-B^{2})^{2}+13(\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B})^{2}],\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle {\mathcal{L}}_{3} & = & \displaystyle \frac{{\it\mu}_{3}}{4{\rm\pi}}[3(E^{2}-B^{2})^{4}+22(E^{2}-B^{2})^{2}(\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B})^{2}+19(\boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B})^{4}],\end{eqnarray}$$

where ${\it\mu}_{1}={\it\alpha}/90{\rm\pi}$ , ${\it\mu}_{2}={\it\alpha}/315{\rm\pi}$ , ${\it\mu}_{3}=4{\it\alpha}/945{\rm\pi}$ and ${\it\alpha}\approx 1/137$ is the fine-structure constant. The term ${\mathcal{L}}_{n}$ corresponds to $2(n+1)$ -photon scattering, or equivalently in this neoclassical approach, vacuum $2(n+1)$ -wave mixing.

To investigate the long propagation behaviour of electromagnetic fields when photon–photon scattering is present, the modified Maxwell equations (2.3) were solved for the case of counterpropagating probe and strong plane waves (depicted in figure 2):

(2.9a,b ) $$\begin{eqnarray}\displaystyle \boldsymbol{E}_{p}^{(0)}({\it\varphi}_{p})={\bf\varepsilon}_{p}\,{\mathcal{E}}_{p}\text{e}^{-({\it\varphi}_{p}/{\it\Phi}_{p})^{2}}\cos {\it\varphi}_{p}\quad \boldsymbol{E}_{s}^{(0)}({\it\varphi}_{s})={\bf\varepsilon}_{s}\,{\mathcal{E}}_{s}\,\text{Rect}\left(\frac{{\it\varphi}_{s}}{{\it\Phi}_{s}}\right), & & \displaystyle\end{eqnarray}$$

where $\text{Rect}({\it\varphi}/{\it\Phi})={\it\theta}({\it\varphi}+{\it\Phi}/2)-{\it\theta}({\it\varphi}-{\it\Phi}/2)$ is the rectangular function, ${\it\theta}(\cdot )$ the Heaviside function, ${\it\varphi}_{p}=k_{p}x={\it\omega}_{p}x^{-}$ , ${\it\varphi}_{s}=k_{s}x={\it\omega}_{s}x^{+}$ , $x^{\pm }=t\pm z$ , ${\it\Phi}_{p}={\it\omega}_{p}{\it\tau}_{p}$ , ${\it\Phi}_{s}={\it\omega}_{s}{\it\tau}_{s}$ with the probe and strong-field wave vectors $\boldsymbol{k}_{p,s}$ and polarisation vectors ${\bf\varepsilon}_{p,s}$ obeying ${\bf\varepsilon}_{p}\boldsymbol{\cdot }{\bf\varepsilon}_{p}=1$ , ${\bf\varepsilon}_{s}\boldsymbol{\cdot }{\bf\varepsilon}_{s}=1$ , $\boldsymbol{k}_{p}\boldsymbol{\cdot }{\bf\varepsilon}_{p}=0$ , $\boldsymbol{k}_{s}\boldsymbol{\cdot }{\bf\varepsilon}_{s}=0$ . Initially, the probe pulse is assumed to be much weaker than the strong background ${\mathcal{E}}_{p}\ll {\mathcal{E}}_{s}$ , although this is just for convenience. We also assume ${\mathcal{E}}_{s}\ll 1$ . Initially, the probe and strong fields are well separated: $\lim _{t\rightarrow -\infty }{\mathcal{F}},{\mathcal{G}}=0$ . The configuration equation (2.9) was chosen to represent a tractable case when the probe is much more rapidly oscillating than the background. The index $(0)$ represents the classical vacuum solution to the wave equation: $[\partial _{t}^{2}-\partial _{z}^{2}]\boldsymbol{E}_{p,s}^{(0)}=\mathbf{0}$ and the total electric field is $\boldsymbol{E}=\boldsymbol{E}_{p}+\boldsymbol{E}_{s}$ . We concentrate on the forward-scattered probe wave and search for plane-wave solutions to (2.3) ( ${\it\omega}_{p}\boldsymbol{B}_{p}=\boldsymbol{k}_{p}\wedge \boldsymbol{E}_{p}$ ).

Figure 2. The calculational and simulational set-up studied. The oscillating probe pulse moves from left to right and the strong field from right to left.

3 Multiple four-photon scattering

To illustrate the technique used, we begin with the familiar example of four-photon scattering. Since for counterpropagating plane waves ${\it\rho}_{vac}=0$ , the probe field can be solved for with the ansatz:

(3.1) $$\begin{eqnarray}\displaystyle \boldsymbol{E}_{p}^{(n+1)}({\it\varphi}_{p},{\it\varphi}_{s})=\boldsymbol{E}_{p}^{(0)}({\it\varphi}_{p})-\int _{-\infty }^{{\it\varphi}_{s}}\frac{\text{d}y}{2}~\boldsymbol{J}_{vac}[\boldsymbol{E}_{p}^{(n)}({\it\varphi}_{p},y),\boldsymbol{E}_{s}^{(0)}(y)]. & & \displaystyle\end{eqnarray}$$

Equation (3.1) is derived by solving iteratively the wave equation for the probe field $[\partial _{t}^{2}-\partial _{z}^{2}]\boldsymbol{E}_{p}=-\partial _{t}\boldsymbol{J}_{vac}[\boldsymbol{E}_{p},\boldsymbol{E}_{s}]$ using the retarded Green’s function $G_{R}(t,z)={\it\theta}(t){\it\theta}(t-|z|)/2$ (Mahan Reference Mahan2002). If $\boldsymbol{J}_{vac}$ is derived from ${\mathcal{L}}_{1}$ in (2.6), then $\boldsymbol{E}_{p}^{(n)}$ is accurate up to order $O({\mathcal{E}}_{s})^{2n}$ , which acts as a perturbative parameter.

A diagrammatic illustration of this iteration is shown in figure 3, where strong-field legs have been suppressed and the virtual pair has been shrunk to an effective interaction vertex. Suppose we ignore the envelope of the probe and choose $\boldsymbol{E}_{p}^{(0)}={\bf\varepsilon}_{p}\,{\mathcal{E}}_{p}\,\cos {\it\varphi}_{p}$ . If the probe and strong-field polarisations are parallel then the asymptotic form ( ${\it\varphi}_{s}\rightarrow \infty$ ) of the probe that would reach detectors, for the first few iterations is (the replacement ${\it\upsilon}_{1}\rightarrow -{\it\upsilon}_{1}$ corrects a sign error in Böhl et al. (Reference Böhl, King and Ruhl2015)):

(3.2) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\boldsymbol{E}_{p}^{(0)}={\bf\varepsilon}_{p}\,{\mathcal{E}}_{p}\,\cos {\it\varphi}_{p},\quad \boldsymbol{E}_{p}^{(1)}={\bf\varepsilon}_{p}\,{\mathcal{E}}_{p}\left[\cos {\it\varphi}_{p}+{\it\upsilon}_{1}\,\sin {\it\varphi}_{p}\right],\\ \boldsymbol{E}_{p}^{(2)}={\bf\varepsilon}_{p}\,{\mathcal{E}}_{p}\left[\cos {\it\varphi}_{p}+{\it\upsilon}_{1}\,\sin {\it\varphi}_{p}-\displaystyle \frac{{\it\upsilon}_{1}^{2}}{2}\,\cos {\it\varphi}_{p}\right],\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where ${\it\upsilon}_{1}={\it\upsilon}_{1}^{\Vert }=(n_{vac}^{\Vert }-1){\it\Phi}$ , $n_{vac}^{\Vert }=1+8{\it\alpha}{\mathcal{E}}_{s}^{2}/45{\rm\pi}$ and ${\it\Phi}={\it\omega}_{p}{\it\tau}_{s}$ . The full solution exponentiates and becomes:

(3.3) $$\begin{eqnarray}\displaystyle \boldsymbol{E}_{p}({\it\varphi}_{p})=\lim _{{\it\varphi}_{s}\rightarrow \infty }\lim _{n\rightarrow \infty }\,\boldsymbol{E}_{p}^{(n)}({\it\varphi}_{p},{\it\varphi}_{s})=\boldsymbol{E}_{p}^{(0)}({\it\varphi}_{p}-{\it\upsilon}_{1}). & & \displaystyle\end{eqnarray}$$

This all-order solution was recently demonstrated using the full polarisation operator (Meuren et al. Reference Meuren, Hatsagortsyan, Keitel and Di Piazza2015).

Figure 3. Arbitrary numbers of four-photon scattering are considered in the propagation of the probe in the strong background.

4 Multiple six-photon scattering

To calculate the change in the probe wave due to only six-wave mixing, (3.1) can be used with $\boldsymbol{J}_{vac}$ derived from ${\mathcal{L}}_{2}$ in (2.7). From the six photon legs of the hexagon diagram in figure 1, one is for the scattered field and three for the strong field, leaving two for the probe. Suppressing the low-frequency strong legs, gives the three-pronged diagram in figure 4.

Due to two photon legs coming from the probe, which is iterated, many more topological structures are produced in this series rather than the simple chain produced for four-photon scattering in figure 3. For higher orders, the series becomes complicated with the $n$ th iteration containing terms from the first to the $(2^{n}-1)$ th perturbative order. For parallel probe and strong field polarisation, the first few iterations are (the replacement ${\it\nu}_{2}\rightarrow -{\it\nu}_{2}$ corrects a sign error in Böhl et al. (Reference Böhl, King and Ruhl2015)):

(4.1) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}\displaystyle \boldsymbol{E}_{p}^{(0)}\hspace{2.0pt}=\hspace{2.0pt}{\bf\varepsilon}_{p}\,{\mathcal{E}}_{p}\,\cos {\it\varphi}_{p}\quad \boldsymbol{E}_{p}^{(1)}={\bf\varepsilon}_{p}\,{\mathcal{E}}_{p}\left[\cos {\it\varphi}_{p}+\frac{{\it\nu}_{2}}{2}\,\sin 2{\it\varphi}_{p}\right]\\ \boldsymbol{E}_{p}^{(2)}\hspace{2.0pt}=\hspace{2.0pt}\displaystyle {\bf\varepsilon}_{p}\,{\mathcal{E}}_{p}\left[\left(1-\frac{1}{2}\left(\frac{{\it\nu}_{2}}{2}\right)^{2}\right)\cos {\it\varphi}_{p}\right.\\ \displaystyle \quad \hspace{24.0pt}+\left.\frac{{\it\nu}_{2}}{2}\,\sin 2{\it\varphi}_{p}-\frac{3}{2}\left(\frac{{\it\nu}_{2}}{2}\right)^{2}\cos 3{\it\varphi}_{p}-\hspace{-2.0pt}\frac{2}{3}\left(\frac{{\it\nu}_{2}}{2}\right)^{3}\sin 4{\it\varphi}_{p}\right],\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where ${\it\nu}_{2}=(n_{vac,2}^{\Vert }-1){\it\Phi}$ was referred to as the ‘shock parameter’ in Böhl et al. (Reference Böhl, King and Ruhl2015) and $n_{vac,2}^{\Vert }=1+192{\it\mu}_{2}{\mathcal{E}}_{s}^{3}{\mathcal{E}}_{p}$ defines a refractive index for six-photon scattering. If the probe has an envelope, this multiplies the parameter. So for (2.9), the nonlinear parameter depends linearly on the probe amplitude and therefore also the envelope and accordingly becomes ${\it\nu}_{2}\exp [-({\it\varphi}_{p}/{\it\Phi}_{p})^{2}]$ . Unlike the familiar ‘nonlinear vacuum refractive index’ in the context of four-photon scattering, when six-photon scattering is considered, the refractive index becomes probe-dependent. This implies the polarised vacuum in this set-up displays a self-focusing or Benjamin–Weir instability (Moloney & Newell Reference Moloney and Newell2004), in which the probe can be focussed to very high intensities before dissipative effects – such as pair creation – break the pulse apart. Due to the similarity with classical plasma instabilities, we liken the polarised vacuum to a ‘vacuum plasma’. The full solution to Maxwell’s equations for the probe scattering due to six-photon scattering is then (the replacement ${\it\nu}_{2}\rightarrow -{\it\nu}_{2}$ leads to a slightly different formula to that published in Böhl et al. (Reference Böhl, King and Ruhl2015), but this does not affect the numerical solution nor the plots presented):

(4.2) $$\begin{eqnarray}\displaystyle \boldsymbol{E}_{p}({\it\varphi}_{p})=2{\bf\varepsilon}_{p}\,{\mathcal{E}}_{p}\mathop{\sum }_{j=1}^{\infty }(-1)^{j+1}\left[\frac{\text{J}_{2j}(2j{\it\nu}_{2})}{2j{\it\nu}_{2}}\sin 2j{\it\varphi}_{p}+\frac{\text{J}_{2j-1}[(2j-1){\it\nu}_{2}]}{(2j-1){\it\nu}_{2}}\cos (2j-1){\it\varphi}_{p}\right], & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $\text{J}_{l}(\cdot )$ is the $l$ th-order Bessel function of the first kind (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2007). We note that this all-order solution resembles the Fubini solution (Fubini-Ghiron Reference Fubini-Ghiron1935) for the propagation of lossless finite-amplitude planar acoustic waves in nonlinear media (Rossing Reference Rossing2007).

Figure 4. Arbitrary numbers of six-photon scattering are considered in the propagation of the probe in the strong background. For higher orders of interaction, the number of possible topological structures becomes increasingly complicated.

5 Numerical simulation of field propagation in the nonlinear polarised vacuum

The numerical method applied is based on the ‘pseudocharacteristic method of lines’ (PCMOL) (Carver Reference Carver1980), matrix inversion to convert the discretised system into a system of ordinary differential equations (ODEs) and the ODE solver CVODE (Hindmarsh et al. Reference Hindmarsh2005).

For simplicity, let us first look at Maxwell equations without quantum corrections. For the set-up of two counterpropagating plane-wave pulses with only transverse polarisations, the linear dynamic Maxwell equation (2.3) with vanishing $\boldsymbol{J}_{vac}=0$ together with the homogeneous equation $\partial _{t}\boldsymbol{B}+\boldsymbol{{\rm\nabla}}\wedge \boldsymbol{E}=0$ reduce to

(5.1) $$\begin{eqnarray}\displaystyle \mathbf{1}_{4}\,\partial _{t}\,\boldsymbol{f}+\unicode[STIX]{x1D64C}\,\partial _{z}\,\boldsymbol{f}=0, & & \displaystyle\end{eqnarray}$$

with the vector $\boldsymbol{f}=(E_{x},E_{y},B_{x},B_{y})^{\text{T}}$ , where $\mathbf{1}_{4}$ is the identity matrix in four dimensions and the coefficient matrix $\unicode[STIX]{x1D64C}=\text{adiag}(1,-1,-1,1)$ is anti-diagonal. Within the PCMOL, the system is transformed to a new basis $\boldsymbol{u}:=\unicode[STIX]{x1D64E}\,\boldsymbol{f}$ in which the new coefficient matrix ${\it\bf\Lambda}=\unicode[STIX]{x1D64E}\unicode[STIX]{x1D64C}\unicode[STIX]{x1D64E}^{-1}=\text{diag}(-1,-1,1,1)$ is diagonal with real eigenvalues. The transformation matrix is given by:

(5.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64E}=\frac{1}{\sqrt{2}}\left(\begin{array}{@{}cccc@{}}-1 & 0 & 0 & 1\\ 0 & 1 & 1 & 0\\ 1 & 0 & 0 & 1\\ 0 & -1 & 1 & 0\end{array}\right). & & \displaystyle\end{eqnarray}$$

In this new basis, the system is decoupled into four independent advection equations:

(5.3) $$\begin{eqnarray}\displaystyle \mathbf{1}_{4}\partial _{t}\left(\begin{array}{@{}c@{}}u_{1}\\ u_{2}\\ u_{3}\\ u_{4}\end{array}\right)+\left(\begin{array}{@{}cccc@{}}-1 & & & \\ & -1 & & \\ & & 1 & \\ & & & 1\end{array}\right)\partial _{z}\left(\begin{array}{@{}c@{}}u_{1}\\ u_{2}\\ u_{3}\\ u_{4}\end{array}\right)=0. & & \displaystyle\end{eqnarray}$$

The new components of $\boldsymbol{u}$ are related to the $E$ - and $B$ fields via

(5.4) $$\begin{eqnarray}\displaystyle \boldsymbol{u}:=\unicode[STIX]{x1D64E}\,\boldsymbol{f}=\frac{1}{\sqrt{2}}\left(\begin{array}{@{}c@{}}B_{y}-E_{x}\\ E_{y}+B_{x}\\ E_{x}+B_{y}\\ B_{x}-E_{y}\end{array}\right). & & \displaystyle\end{eqnarray}$$

The sign of the eigenvalues ${\it\lambda}_{i}$ of $\boldsymbol{{\it\Lambda}}$ determine the direction in which the component $u_{i}$ propagates. For a positive sign, the component $u_{i}$ travels to the right, for ${\it\lambda}_{i}<0$ it travels to the left. Now the system is discretised in space by introducing a co-located grid with spacing ${\rm\Delta}z=L/(N-1)$ where $L$ is the length of the simulation box and $N$ is the number of spatial points. Each grid point is thereby assigned the field values $\boldsymbol{u}^{l}=(u_{1}^{l},u_{2}^{l},u_{3}^{l},u_{4}^{l})$ where $u_{i}^{l}:=u_{i}({\rm\Delta}z\cdot l)$ and $0<l\leqslant N$ is the index on the grid. Following (Schiesser Reference Schiesser1991), the spatial derivatives are approximated with biased five-point-stencils as, for a fixed grid method, they show good behaviour with minimal oscillations where high gradients in the solution are present and thereby introduce only minimal numerical diffusion. Due to the discretisation, we are left with a $4N$ -dimensional system of ordinary differential equations, $\tilde{\boldsymbol{u}}^{\prime }(t)=g[\tilde{\boldsymbol{u}}]$ with $g[\tilde{\boldsymbol{u}}]=-\tilde{\boldsymbol{{\it\Lambda}}}\,\unicode[STIX]{x1D63F}\,\tilde{\boldsymbol{u}}$ . The prime $^{\prime }$ denotes the time derivative, $\unicode[STIX]{x1D63F}$ is the matrix representing the finite differencing and $\tilde{\boldsymbol{{\it\Lambda}}}=\mathbf{1}_{N}\otimes \boldsymbol{{\it\Lambda}}$ where $\otimes$ is the Kronecker product. We use the tilde ${\sim}$ to denote quantities on the grid. The system is then integrated using CVODE. The initial conditions are set-up in the $\tilde{\boldsymbol{f}}$ -basis, integrated in $\tilde{\boldsymbol{u}}$ and then transformed back to $\tilde{\boldsymbol{f}}$ . The number of grid points is thereby chosen sufficiently high to minimise the effect of numerical diffusion and properly resolve steep gradients in the solution. Since the $\text{Rect}$ -function in (2.9) actually shows an infinite slope, it is approximated in the simulation by some mirrored Fermi–Dirac-function $\text{FD}(y)$ which is given by

(5.5) $$\begin{eqnarray}\displaystyle \text{Rect}(y)\approx \text{FD}(y)=\frac{1}{1+\exp \left(\displaystyle \frac{|y|-{\it\omega}_{s}z_{m}}{{\it\omega}_{s}z_{b}}\right)}. & & \displaystyle\end{eqnarray}$$

The parameters $z_{b}$ and $z_{m}$ can be thought of as the ‘temperature’ and ‘chemical potential’ which control the width and steepness of the pulse. The simulational set-up is shown in figure 2.

For the nonlinear case, $\boldsymbol{J}_{vac}$ contains nonlinear corrections with spatial and temporal derivatives of polynomials in the fields $f_{i}$ . Using the chain rule, one can always write the derivative of monomials like $\partial _{t}(f_{i})^{n}$ of the field components as some (field dependent) coefficient times the linear derivative, e.g. $\partial _{t}(f_{i})^{n}=n(f_{i})^{n-1}\partial _{t}\,f_{i}$ . In this way, the nonlinear corrections to (5.1) can be written as:

(5.6) $$\begin{eqnarray}\displaystyle (\mathbf{1}_{4}+\unicode[STIX]{x1D653})\partial _{t}\boldsymbol{f}+(\unicode[STIX]{x1D64C}+\unicode[STIX]{x1D654})\partial _{z}\boldsymbol{f}=0. & & \displaystyle\end{eqnarray}$$

For the detailed form of $\unicode[STIX]{x1D653}$ and $\unicode[STIX]{x1D654}$ the reader is referred to Böhl et al. (Reference Böhl, King and Ruhl2015). In the discretised version, the matrices $\unicode[STIX]{x1D653}$ and $\unicode[STIX]{x1D654}$ also become $4N\times 4N$ dimensional matrices. Since both only depend on local field values, they are of block diagonal form:

(5.7a,b ) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D653}}=\bigoplus _{l=1}^{N}\unicode[STIX]{x1D653}^{l},\quad \tilde{\unicode[STIX]{x1D654}}=\bigoplus _{l=1}^{N}\unicode[STIX]{x1D654}^{l}, & & \displaystyle\end{eqnarray}$$

where $\bigoplus$ is the direct sum and $\unicode[STIX]{x1D653}^{l}$ and $\unicode[STIX]{x1D654}^{l}$ are the matrices $\unicode[STIX]{x1D653}$ and $\unicode[STIX]{x1D654}$ evaluated at grid point $l$ . It remains to invert the matrix $(\mathbf{1}_{4N}+\tilde{\unicode[STIX]{x1D653}})$ to bring the discretised version of (5.6) into the form of a set of ODEs. Since the inverse of a block diagonal matrix is also block diagonal, the problem reduces to $N$ times the inversion of the matrix $\mathbf{1}_{4}+\unicode[STIX]{x1D653}^{l}$ . To further reduce the computational cost, we rewrite $\unicode[STIX]{x1D653}^{l}$ as $\unicode[STIX]{x1D653}^{l}=\unicode[STIX]{x1D642}\,\unicode[STIX]{x1D643}^{l}$ with

(5.8a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D642}=\left(\begin{array}{@{}cc@{}}1 & 0\\ 0 & 1\\ 0 & 0\\ 0 & 0\end{array}\right),\quad \unicode[STIX]{x1D643}^{l}=\left(\begin{array}{@{}cccc@{}}x_{11}^{l} & x_{12}^{l} & x_{13}^{l} & x_{14}^{l}\\ x_{21}^{l} & x_{22}^{l} & x_{23}^{l} & x_{24}^{l}\\ \end{array}\right).\end{eqnarray}$$

The inversion can then be further reduced to

(5.9) $$\begin{eqnarray}\displaystyle (\mathbf{1}_{4}+\unicode[STIX]{x1D653}^{l})^{-1}=\mathbf{1}_{4}-\unicode[STIX]{x1D642}(\mathbf{1}_{2}+\unicode[STIX]{x1D643}^{l}\unicode[STIX]{x1D642})^{-1}\unicode[STIX]{x1D643}^{l}, & & \displaystyle\end{eqnarray}$$

where we used the Woodbury formula (Golub & Van Loan Reference Golub and Van Loan2012). The inverse of the $2\times 2$ matrix is calculated via a LU-factorisation (Lower Upper factorisation) for each grid point at each evaluation of the now nonlinear function $g[\tilde{\boldsymbol{u}}]$ in $\tilde{\boldsymbol{u}}^{\prime }(t)=g[\tilde{\boldsymbol{u}}]$ with

(5.10) $$\begin{eqnarray}\displaystyle \tilde{\boldsymbol{u}}^{\prime }(t)=-\tilde{\unicode[STIX]{x1D64E}}(\mathbf{1}_{4N}+\tilde{\unicode[STIX]{x1D653}})^{-1}(\tilde{\unicode[STIX]{x1D64C}}+\tilde{\unicode[STIX]{x1D654}})\tilde{\unicode[STIX]{x1D64E}}^{-1}\unicode[STIX]{x1D63F}\tilde{\boldsymbol{u}}, & & \displaystyle\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D64E}}=\mathbf{1}_{N}\otimes \unicode[STIX]{x1D64E}$ and $\tilde{\unicode[STIX]{x1D64C}}=\mathbf{1}_{N}\otimes \unicode[STIX]{x1D64C}$ . Since the nonlinear corrections do not change the signs of the eigenvalues ${\it\lambda}_{i}$ for the parameters considered here, we use the same biased differencing $\unicode[STIX]{x1D63F}$ as in the linear case. The signals are analysed using a spatial Fourier transform in Wolfram Mathematica (Wolfram Research 2012) under the assumption ${\it\omega}=|\boldsymbol{k}|$ .

Figure 5. Harmonic spectra generated by six-photon scattering of a probe, polarised parallel to the background, for shock parameter ${\it\nu}_{2}=0.05,0.6,1$ (left-to-right). The squares show the leading-order perturbative term for that harmonic, the smooth solid brown line is the all-order analytical solution equation (4.2) and the jagged solid red line is from numerical simulation.

Figure 6. (a) For parallel probe and strong polarisations, a shock wave accompanies the generation of higher harmonics as ${\it\nu}_{2}\rightarrow 1$ . (b) For perpendicular probe and strong polarisations, different shock waves are generated in the field parallel ( $E_{p}^{\Vert }$ ) and perpendicular ( $E_{p}^{\bot }$ ) to the probe. This shock can change form depending on the ratio of four- to six-photon scattering. Here ${\it\upsilon}_{1}=5$ , ${\it\nu}_{2}=1$ and the shock resembles that from a nonlinear Kerr medium.

6 Vacuum high-harmonic generation

In general, both four- and six-photon scattering will be present. The extent of four-photon scattering can be quantified with ${\it\upsilon}_{1}$ and six-photon scattering with ${\it\nu}_{2}$ . Since we are interested in the case ${\mathcal{E}}_{p}\ll {\mathcal{E}}_{s}\ll 1$ , we must consider when ${\it\upsilon}_{1}$ is larger than ${\it\nu}_{2}$ . For the case of parallel probe and strong polarisations, the produced harmonic spectrum is similar to that in the limit ${\it\upsilon}_{1}\rightarrow 0$ given by (4.2). The numerical and analytical solutions are plotted for ${\it\nu}_{2}=0.05,0.6,1$ in figure 5. As ${\it\nu}_{2}\rightarrow 1$ , higher harmonics move from being exponentially to only power-law suppressed and the shock wave in figure 6(a) is generated. The amplitude of the field of the $n$ th harmonic relative to the initial probe scales as ${\sim}({\mathcal{E}}_{p}{\mathcal{E}}_{s}^{3}{\it\Phi})^{n-1}$ for multiple six-photon scattering whereas it scales as ${\sim}{\mathcal{E}}_{p}^{n-1}{\mathcal{E}}_{s}^{n}{\it\Phi}$ for single $2n$ -photon scattering. If the amplitude of the $n$ th harmonic is not to be vanishingly small for large $n$ , the latter route requires ${\mathcal{E}}_{p}$ and ${\mathcal{E}}_{s}$ being close to $1$ for large $n$ , whereas the former route of repeated six-photon scattering allows for the possibility of smaller field strengths being compensated for by a long phase length ${\it\Phi}$ .

For the case of perpendicular probe and strong polarisation, a shock wave also accompanies high-harmonic generation, but both are suppressed compared to the parallel case. This is due to the problem of phase matching. Since the vacuum is birefringent, a probe photon polarised parallel to the background will propagate with a different phase to a probe photon polarised perpendicular to the background. To reach higher harmonics, the probe photon scatters multiple times. The higher the harmonic, the more channels there are to reach that harmonic. These channels involve the photon being in intermediate even (polarisation parallel to the background) and odd (polarisation perpendicular to the background) harmonics to differing degrees. This leads to contributions for higher harmonics from different channels being added incoherently and hence less efficiently than in the case of parallel probe and background polarisations.

The effect of dispersion is most clear when probe and background polarisations are perpendicular. For the non-dispersive case ( ${\it\upsilon}_{1}\approx 0$ ), the shock generated in the odd harmonics is square-wave-like and in the even harmonics is sawtooth like (Böhl et al. Reference Böhl, King and Ruhl2015). However, for the dispersive case of ${\it\upsilon}_{1}=5$ , when ${\it\nu}_{2}=1$ , the shock in figure 6(b) was observed in the numerical solution. This is reminiscent of shocks generated in a Kerr medium, which contains a cubic nonlinearity in the polarisation $P_{i}={\it\chi}_{ij}^{(1)}E_{j}+{\it\chi}_{ijkl}^{(3)}E_{j}E_{k}E_{l}$ . In the strongly dispersive case ${\it\upsilon}_{1}\gg {\it\nu}_{2}$ , the nature of the scattered field depends more sensitively on the ratio of ${\it\upsilon}_{1}$ to ${\it\nu}_{2}$ . The probe pulse envelope has frequency components of the order of ${\sim}1/{\it\tau}_{p}$ , so when the lag induced by four-wave mixing for these frequencies becomes of the order of or larger than unity i.e. when ${\it\upsilon}_{1}\gtrsim {\it\omega}_{p}{\it\tau}_{p}$ , we find the pulse envelope itself is deformed.

The set-up in figure 2 is useful for studying the physics of vacuum high-harmonic generation as the probe propagates through a half-cycle of the more slowly varying background. However, in any laser-based experiment, one should take in account at least a complete cycle of the background field. For this special issue, we include the set-up in figure 7 of a single sinusoidal cycle in an otherwise square-wave background and compare the harmonics generated with the constant background case of figure 2. Since the $n$ th occurrence of high-harmonic generation through six-photon scattering of the probe scales as ${\sim}(E_{s}^{3}E_{p}^{2})^{n}$ , one would expect even harmonics to be heavily suppressed due to integrating over an odd power of the oscillating background $E_{s}$ . This is displayed in the spectrum in figure 7, which shows also an overall suppression compared with the square-wave case. $n=2$ , which contributes to the fundamental and third harmonic of the probe, is the leading-order term that contains a constant, which does not disappear when the background is integrated over. Since the corresponding constant from the $E_{s}^{6}$ term ${\sim}\sin ^{6}{\it\varphi}$ is equal to $5/16$ , we choose to compare the spectrum of ${\it\nu}_{2}=1$ for a square-wave background to ${\it\nu}_{2}=(16/5)^{1/2}$ (since ${\it\nu}_{2}\sim E_{s}^{3}$ and for $n=2$ is for $E_{p}\sim {\it\nu}_{2}^{2}$ ) for the single-cycle sine wave background. The result in figure 7 shows (i) the harmonic spectrum is considerably suppressed for this set-up; (ii) the third and fifth harmonics are less suppressed than other harmonics and (iii) for this parameter, several high harmonics are no longer exponentially suppressed with respect to one another. This demonstrates that even though ${\it\nu}_{2}$ was increased to compensate for the integration over the sinusoidal shape of the new background, harmonic generation is greatly suppressed. The constant that arises due to this integration is $[1+(-1)^{n}](3n)!/[(3n/2)!]^{2}/2^{3n+1}$ . When the increased value of ${\it\nu}_{2}$ is included, this coefficient for the $n$ th occurrence of six-photon scattering increases with $n$ . This suggests that the background pulse shape plays a more complicated role than a simple power scaling.

Figure 7. (a) The modified set-up of a right-moving oscillating probe with envelope and a left-moving background with a single sinusoidal cycle. (b) A comparison of the harmonics produced in a constant background (figure 2) for ${\it\nu}_{2}=1$ and in an oscillating background for ${\it\nu}_{2}=(16/5)^{1/2}$ . ‘Box’ and ‘Hexagon’ refer to the four- and six-photon-scattering shown in figure 1.

Even though the background now has a frequency ${\it\omega}_{s}$ which is not much smaller than the probe frequency ${\it\omega}_{p}$ , we do not observe the phenomenon of vacuum frequency mixing. Although the merging of $n_{p}$ probe photons and $n_{s}$ background photons would generate a range of frequencies, for example: $n_{p}{\it\omega}_{p}+n_{s}{\it\omega}_{s}$ by frequency mixing, it can be shown by energy–momentum conservation for our head-on set-up, that all mixed frequencies with $n_{s}\neq 0$ cannot propagate.

Another source of effects on the probe harmonic spectrum is the probe pulse duration. For this special issue, we have simulated the collision of a few-cycle probe ( ${\it\Phi}_{p}=5$ in (2.9)) with a constant background. The original probe, the generated shockwave and its corresponding harmonic spectrum are plotted in figure 8. The result of a shorter probe pulse and hence a wider initial probe spectrum can be seen by contrasting the generated spectrum in figure 8 with that for a long-cycle probe in figure 5. Although the plot in figure 8 only includes six-photon scattering, the result is similar (but phase shifted) when four-photon scattering is taken into account. The effect of the wider probe spectrum is to smooth out the harmonic structure, but to maintain the same power-law exponent. This is straightforward to understand. For a long pulse the spectrum is quasi-monochromatic so the generated frequencies are well defined around the $n$ th harmonic $n{\it\omega}_{p}$ . For a short pulse, probe photons have a range of frequencies so the merging of $n$ photons gives a signal at $\sum _{l=1}^{n}{\it\omega}_{p,l}$ , where ${\it\omega}_{p,l}$ are the individual frequencies taken from the probe spectrum. That the gradient of the spectrum is very similar to that for a long pulse is also what our original analysis would have predicted, so is valid here too.

Figure 8. (a) The shock wave generated when a few cycle ( ${\it\Phi}_{p}=5$ ) probe pulse with counterpropagates through a constant background. (b) The corresponding harmonic spectrum for ${\it\nu}_{2}=1$ (blue thick line), the harmonic spectrum for the first set-up with a long pulse ( ${\it\Phi}_{p}\gg 1$ ) when ${\it\nu}_{2}=1$ (black thin line) and original probe spectrum (red dashed line).