Vacuum high-harmonic generation and electromagnetic shock

Abstract

When one takes into account the presence of virtual charged states in the quantum vacuum, a nonlinear self-interaction can arise in the propagation of electromagnetic fields. This self-interaction is often referred to as ‘real photon–photon scattering’. When the centre-of-mass energy of colliding photons is much lower than the rest energy of an electron–positron pair, this quantum effect can be included in the classical field equations of motion as a vacuum current and charge density using the Heisenberg–Euler Lagrangian. Using analytical and numerical methods for subcritical fields, the intrinsic solution to Maxwell’s equations has been found for counterpropagating probe and pump plane waves in the presence of vacuum four- and six-wave mixing. In the corresponding all-order solution for the scattered probe, a route to vacuum high-harmonic generation is identified in which a long phase length can compensate for the weakness of interacting fields. The resulting shocks in the probe carrier wave and envelope are studied for different parameter regimes and polarisation set-ups. In this special issue, we study two additional set-ups: that of a slowly varying single-cycle background to highlight the effect of an oscillating background on the probe harmonic spectrum, and that of a few-cycle probe to highlight the smoothing of the harmonic peaks produced by a wider spectrum of probe photons. We also correct sign errors in an earlier publication.

1 Introduction

That real photons scatter off one another due to their interaction on mutual virtual pairs has been predicted to occur for over eighty years (Sauter 1931; Halpern 1934), with the first low-energy calculations performed by Euler & Kochel (1935), Heisenberg & Euler (1936) and, independently, Weisskopf (1936). Recent advances in laser technology (Danson, Hillier, Hopps & Neely 2015) have generated much interest in measuring this effect using high-power laser pulses (reviews can be found in Marklund & Shukla (2006), Di Piazza et al. (2012) and King & Heinzl (2016)), in which the high flux of photons compensates for the very small cross-section (Berestetskii, Lifshitz & Pitaevskii 1982). In such scenarios where the centre-of-mass energy of the colliding photons is much less than the electron rest energy, an effective interaction based on the Heisenberg–Euler Lagrangian in which fermion dynamics have been integrated out, is a good approximation (King 2010). Most scenarios suggested to measure real photon–photon scattering involve a stronger, polarising ‘background’ pump pulse in which a weaker ‘probe’ pulse is scattered. The leading-order signal of real photon–photon scattering in such a laser-based set-up is well described by a four-wave mixing process, which provides signals of modified probe angular dispersion (Di Piazza, Hatsagortsyan & Keitel 2006; Heinzl et al. 2006; King, Di Piazza & Keitel 2010b ; Kryuchkyan & Hatsagortsyan 2011; Monden & Kodama 2011; King & Keitel 2012), polarisation (King, Di Piazza & Keitel 2010a ; Dinu et al. 2014a , b ) and frequency (Lundström et al. 2006; King & Keitel 2012). Four-wave mixing has been the subject of several experiments (Bernard et al. 2000; Zavattini et al. 2012; Cadène et al. 2014) and is a planned experiment at the HIBEF (Helmholtz International Beamline for Extreme Fields) on the European XFEL facility at DESY (Schlenvoigt et al. 2016).

In contrast, of particular interest in this project was the production of higher frequencies through real vacuum photon–photon scattering. Such ‘vacuum high-harmonic generation’ has been studied for a range of field strengths and frequencies by considering single events of $2n$ -photon scattering (Di Piazza, Hatsagortsyan & Keitel 2005; Fedotov & Narozhny 2006), where field strengths of the order of the critical field were calculated as required for the process to be observable. In the current study (King, Böhl & Ruhl 2014; Böhl, King & Ruhl 2015), an alternative route to producing higher harmonics in a vacuum by using arbitrary numbers of four- and six-photon scattering, was calculated. This can be achieved by solving the wave equation in the presence of the polarised vacuum current. Such a calculation is similar to studies of a plane wave propagating in a constant magnetic background, relevant in astrophysical contexts (Bialynicka-Birula 1981; Zheleznyakov & Fabrikant 1982; Rozanov 1993; Heyl & Hernquist 1998, 1999) but appears to have been performed for the first time for a plane-wave background.

In this special issue, we build upon results in Böhl et al. (2015) by considering set-ups more relevant for collisions of intense laser pulses. First, we consider the head-on collision of an oscillating probe pulse with a single-cycle background to exhibit the change in the probe’s harmonic spectrum when the background, however slowly varying, oscillates an integer number of cycles. Second, we simulate a few-cycle probe pulse colliding head on with a constant background to show how, even though the harmonic peaks become less well defined, the shock-wave effect persists and the power-law behaviour is identical.

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 1936) 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. 2008), 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 2006)). 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 2002). 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. (2015)):

(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. 2015).

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. (2015)):

(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. (2015) 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 2004), 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. (2015), 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 2007). We note that this all-order solution resembles the Fubini solution (Fubini-Ghiron 1935) for the propagation of lossless finite-amplitude planar acoustic waves in nonlinear media (Rossing 2007).

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 1980), matrix inversion to convert the discretised system into a system of ordinary differential equations (ODEs) and the ODE solver CVODE (Hindmarsh et al. 2005).

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 1991), 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. (2015). 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 2012). 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. 2015). 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).

7 Discussion

It was found that when the nonlinear parameter ${\it\nu}_{2}=64{\it\alpha}{\mathcal{E}}_{s}^{3}{\mathcal{E}}_{p}{\it\omega}_{p}{\it\tau}_{s}/105{\rm\pi}$ tended towards unity, the high-harmonic spectrum produced in an oscillating probe wave counterpropagating with a much slower varying background moved from an exponential suppression of higher harmonics to just a power-law suppression. A similar transition is found for increasing field strength in high-harmonic generation in laser-irradiated foils. The ‘oscillating mirror model’ of laser–foil interactions has been used to predict different power-law exponents ${\it\gamma}$ , for example ${\it\gamma}=-5/2$ (Gordienko et al. 2004) and ${\it\gamma}=-8/3$ (Baeva, Gordienko & Pukhov 2006) have been postulated, and experiments on solid targets have recorded intensity-dependent power-law exponents, for example in Norreys et al. (1996) of $-5.50<{\it\gamma}<-3.38$ . For vacuum six-wave mixing we found power-law exponents $-4<{\it\gamma}<-2.4$ for $0.9<{\it\nu}_{2}<1.4$ . Also similar is the generation of even harmonics with a polarisation perpendicular to the probe and odd harmonics with a polarisation parallel to the probe. In the current study, this was observed when the initial probe and strong field polarisations were perpendicular, whereas it is observed in laser-irradiated foil experiments for harmonics generated with $s$ - (perpendicular to plane of incidence) and $p$ - (parallel to plane of incidence) polarisation (Lichters, Meyer-ter-Vehn & Pukhov 1996).

The presented work has made an analogy of photon–photon scattering in counterpropagating plane waves with the interaction of a probe and a ‘vacuum plasma’, and found a shock-wave instability for this idealised case. The effect on the probe due to six-photon scattering in the parallel-polarised background was found to be describable in terms of a probe-dependent vacuum refractive index, which when added to the usual four-photon refractive index gives:

(7.1) $$\begin{eqnarray}\displaystyle n_{vac}^{\Vert }=1+\frac{{\it\alpha}\,E_{s}^{2}}{{\rm\pi}}\left[\frac{8}{45}+\frac{64}{105}\,E_{p}E_{s}\right]. & & \displaystyle\end{eqnarray}$$

Other than four-photon scattering, no dispersive mechanism, such as pair creation from vacuum (Schwinger 1951) or photon-seeded pair creation (Nikishov & Ritus 1964) have been included, nor have transverse degrees of freedom. Therefore to what extent the polarised vacuum can indeed be used as a lens to focus ultra-intense laser pulses to the highest intensities remains an uncertain but tantalising question.