Generation of high order harmonics in Heisenberg–Euler electrodynamics

Present and forthcoming high power laser developments [1–3] open a vast area in the study of nonlinear physics phenomena related to the electromagnetic field interaction with matter and vacuum [4–12]. Vacuum polarization effects leading to high order harmonics generation in vacuum in the course of the collision of the extreme intensity laser pulses is one of the more important directions in these investigations [13–16]. This problem has attracted substantial attention [17–29] because it sheds light on the dynamical properties of quantum electrodynamics (QED) vacuum in strong electromagnetic fields. The generation of high order harmonics plays an important role in the steepening of a nonlinear electromagnetic wave, in the intersecting of strong laser beams [30] and in the formation of the relativistic electromagnetic solitons in the QED vacuum [31, 32].

In spite of a number of publications [19, 22–27] devoted to the detailed theoretical analysis of the harmonic generation, the theory of this process is far from complete. There are several theoretical questions that need to be more thoroughly clarified. They concern the role of the various small parameters governing the problem that should be ordered and considered simultaneously, taking into account a specific form of the Heisenberg–Euler Lagrangian [33–36] describing the QED vacuum polarization effects.

The present paper is devoted to the theoretical consideration and clarification of the high order harmonics generation mechanism within the framework of the Heisenberg–Euler electrodynamics. We consider here the case where the high order harmonics are generated by two crossing plane electromagnetic waves in vacuum under the conditions where only one of the two Poincaré invariants of the electromagnetic field do not vanish. We assume that the finite length electromagnetic pulses are infinite in the two transverse directions. Applying a Lorentz transformation we can find a boosted frame of reference where the crossing electromagnetic pulses appear in the form of two counter-propagating waves (e.g. see [38]). Thus the problem becomes one-dimensional, non-stationary and is described by two independent variables. We consider the case of electromagnetic waves with wavelengths much larger than the electron Compton scattering length and with amplitudes substantially lower than the Schwinger field, which corresponds to the conditions of validity of the Heisenberg–Euler approach [33–36]. In this context we recall that the problem of extending the Heisenberg–Euler electrodynamics has attracted much interest from the quantum field theory community. One of the ways along which such a generalization can be performed is by taking into account the two-loop correction to the ground state energy of the virtual electron-positron sea in an almost constant external electromagnetic field, see for example [37] and references therein. Here we prefer to consider only the lowest non-vanishing order in terms of powers of fine-structure constant α, assuming the amplitudes of in-coming waves in terms of the Schwinger field as given. This approach allows us to restrict our investigation to the original Heisenberg–Euler Lagrangian and to simplify considerably our results in the framework of the Heisenberg–Euler electrodynamics approach.

The paper is organized as follows. In section 2 we present the nonlinear wave equations within the framework of the Heisenberg–Euler electrodynamics which are used throughout the paper. Section 3 describes the electromagnetic configuration and the Dirac light-cone coordinates convenient for analysing the problem under consideration. Section 4 is devoted to the derivation of a convenient form of nonlinear wave equation. The formulation of the scattering problem and the introduction of a perturbation theory which is used for obtaining final results are in section 5. In section 6 we write the expressions giving the intensity of the high order harmonics. Section 7 is devoted to a general case corresponding to different polarizations. The discussion and the summary of the results obtained are presented in section 8.

The analysis of the electromagnetic wave interaction in the QED vacuum is based on the Heisenberg–Euler Lagrangian density, ${\mathcal{L}}_{\text{HE}}$ [33–36]. The sum of the Lagrangians,

$$\begin{equation}\mathcal{L}={\mathcal{L}}_{0}+{\mathcal{L}}_{\text{HE}},\end{equation} \tag{ 1 }$$

describes the electromagnetic field in the long-wavelength limit. Here

$$\begin{equation}{\mathcal{L}}_{0}=-\frac{{m}^{4}}{16\pi \alpha }{F}_{\mu \nu }{F}^{\mu \nu }\end{equation} \tag{ 2 }$$

is the Lagrangian of classical electrodynamics with the electromagnetic field tensor F_μν defined in terms of the four-vector potential A_μ as [39]

$$\begin{equation}{F}_{\mu \nu }={\partial }_{\mu }{A}_{\nu }-{\partial }_{\mu }{A}_{\nu }.\end{equation} \tag{ 3 }$$

Here and in the following we use dimensionless variables with ℏ = c = 1, m_e and e are the electron mass and elementary electric charge, and α = e²/(4π) ≈ 1/137 is the fine structure constant. Thus the electromagnetic fields are measured in units of ${m}_{\text{e}}^{2}/e$, i.e. are normalized on the QED critical field ${E}_{\text{S}}={m}_{\text{e}}^{2}{c}^{3}/e\hslash =1.32\times 1{0}^{18}\enspace \text{V}\enspace {\text{m}}^{-1}$.

In the Heisenberg–Euler theory, the radiative corrections are described by the ${\mathcal{L}}_{\text{HE}}$ term on the right-hand side of equation (1). It can be written as [36]

$$\begin{align}{\mathcal{L}}_{\text{HE}}& =\frac{{m}_{\text{e}}^{4}}{8{\pi }^{2}}\mathcal{M}(\mathfrak{e},\mathfrak{b})=\frac{{m}^{4}}{8{\pi }^{2}}{\int }_{0}^{\infty }\frac{\mathrm{exp}(-\eta )}{{\eta }^{3}}\left[-(\eta \mathfrak{e}\enspace \mathrm{cot}\enspace \eta \mathfrak{e})(\eta \mathfrak{b}\enspace \mathrm{coth}\enspace \eta \mathfrak{b})+1-\frac{{\eta }^{2}}{3}({\mathfrak{e}}^{2}-{\mathfrak{b}}^{2})\right]\mathrm{d}\eta .\end{align} \tag{ 4 }$$

The invariant fields $\mathfrak{e}$ and $\mathfrak{b}$ are expressed in terms of the Poincaré invariants

$$\begin{equation}\mathfrak{F}=\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad \mathfrak{G}=\frac{1}{4}{F}_{\mu \nu }{\tilde{F}}^{\mu \nu }\end{equation} \tag{ 5 }$$

as

$$\begin{equation}\mathfrak{e}=\sqrt{\sqrt{{\mathfrak{F}}^{2}+{\mathfrak{G}}^{2}}-\mathfrak{F}}\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad \mathfrak{b}=\sqrt{\sqrt{{\mathfrak{F}}^{2}+{\mathfrak{G}}^{2}}+\mathfrak{F}},\end{equation} \tag{ 6 }$$

respectively. The dual tensor ${\tilde{F}}^{\mu \nu }$ is defined by

$$\begin{equation}{\tilde{F}}^{\mu \nu }=\frac{1}{2}{\varepsilon }^{\mu \nu \rho \sigma }{F}_{\rho \sigma },\end{equation} \tag{ 7 }$$

where ɛ^μνρσ is the Levi-Civita symbol in four dimensions (ɛ⁰¹²³ = −ɛ₀₁₂₃ = −1). In the 3D notation the Poincaré invariants are

$$\begin{equation}\mathfrak{F}=\frac{1}{2}\left({\mathbf{B}}^{2}-{\mathbf{E}}^{2}\right),\qquad \mathfrak{G}=\mathbf{B}\cdot \mathbf{E}.\end{equation} \tag{ 8 }$$

As explained in reference [36], the Heisenberg–Euler Lagrangian in the form given by equation (4) should be used for obtaining an asymptotic series over the field invariant $\mathfrak{e}$ and $\mathfrak{b}$ assuming that they are small. In this limit, the function $\mathcal{M}(\mathfrak{e},\mathfrak{b})$ in equation (4) can be expanded for small arguments as:

$$\begin{align}\mathcal{M}(\mathfrak{e},\mathfrak{b})& =\frac{{\Gamma}(2)}{45}\left[({\mathfrak{e}}^{4}+{\mathfrak{b}}^{4})+5{\mathfrak{e}}^{2}\enspace {\mathfrak{b}}^{2}\right]-\frac{{\Gamma}(4)}{945}\left[2({\mathfrak{b}}^{6}-{\mathfrak{e}}^{6})+7{\mathfrak{e}}^{2}\enspace {\mathfrak{b}}^{2}({\mathfrak{b}}^{2}-{\mathfrak{e}}^{2})\right]\\ & \quad +\frac{{\Gamma}(6)}{7\times 4{5}^{2}}\left[3\left({\mathfrak{e}}^{8}+{\mathfrak{b}}^{8}\right)+10{\mathfrak{e}}^{2}{\mathfrak{b}}^{2}\left({\mathfrak{e}}^{4}+{\mathfrak{b}}^{4}\right)-7{\mathfrak{e}}^{4}{\mathfrak{b}}^{4}\right]+\dots \enspace .\end{align} \tag{ 9 }$$

Here Γ(x) is the Euler gamma function [40]. The properties of the expansion of $\mathcal{M}(\mathfrak{e},\mathfrak{b})$ in series of powers of the fields $\mathfrak{e}$ and $\mathfrak{b}$ are discussed in reference [41]. The expression (9) yields for the Lagrangian ${\mathcal{L}}_{\text{HE}}$ in the weak field approximation (see also reference [43])

$$\begin{align}{\mathcal{L}}_{\text{HE}}=\kappa \mathcal{M}& =\frac{\kappa }{45}\left(4{\mathfrak{F}}^{2}+7{\mathfrak{G}}^{2}\right)-\frac{4\kappa }{315}\mathfrak{F}\left(8{\mathfrak{F}}^{2}+13{\mathfrak{G}}^{2}\right)+\frac{8\kappa }{945}\left(48{\mathfrak{F}}^{4}+88{\mathfrak{F}}^{2}{\mathfrak{G}}^{2}+19{\mathfrak{G}}^{4}\right)+\dots \end{align} \tag{ 10 }$$

with the constant κ = m⁴/8π². In the Lagrangian given by equation (10) the first, second and third terms on the right-hand side correspond to four-, six- and eight-photon interactions, respectively.

As explained in the introduction we consider two counter propagating waves along the x-axis. The vector potential corresponding to the electromagnetic field can be presented in the form

$$\begin{equation}\mathbf{A}={A}_{z}(t,x){\mathbf{e}}_{z}+{A}_{y}(t,x){\mathbf{e}}_{y}.\end{equation} \tag{ 11 }$$

Below we use Dirac's light cone coordinates x⁻ and x⁺ defined as (see e.g. reference [45])

$$\begin{equation}{x}^{+}=\frac{x+t}{\sqrt{2}},\qquad {x}^{-}=\frac{t-x}{\sqrt{2}}.\end{equation} \tag{ 12 }$$

The Lorentz transform of the light-cone variables under a boost along x with velocity β is given by

$$\begin{equation}{x}^{\prime +}=\frac{{x}^{\prime }+{t}^{\prime }}{\sqrt{2}}={\text{e}}^{-\eta }\frac{x+t}{\sqrt{2}}={\text{e}}^{-\eta }{x}^{+},\end{equation} \tag{ 13 }$$

$$\begin{equation}{x}^{\prime -}=\frac{{t}^{\prime }-{x}^{\prime }}{\sqrt{2}}={\text{e}}^{+\eta }\frac{t-x}{\sqrt{2}}={\text{e}}^{+\eta }{x}^{-},\end{equation} \tag{ 14 }$$

with a prime denoting variables in the boosted frame of reference and η equal to

$$\begin{equation}\eta =\mathrm{ln}\sqrt{\frac{1+\beta }{1-\beta }}.\end{equation} \tag{ 15 }$$

The following shorthand stands for the partial derivatives:

$$\begin{equation}{\partial }_{+}f={\left(\frac{\partial f}{\partial {x}^{+}}\right)}_{{x}^{-}},\qquad {\partial }_{-}f={\left(\frac{\partial f}{\partial {x}^{-}}\right)}_{{x}^{+}}.\end{equation} \tag{ 16 }$$

Then the derivatives are transformed as:

$$\begin{equation}{\left({\partial }_{-}\right)}^{\prime }={\text{e}}^{-\eta }{\partial }_{-}\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad {\left({\partial }_{+}\right)}^{\prime }={\text{e}}^{+\eta }{\partial }_{+},\end{equation} \tag{ 17 }$$

respectively.

For simplicity we consider in this section interacting electromagnetic waves of the same polarization. The general polarization case is considered in section 7. In the case of the same polarizations, the invariant $\mathfrak{G}$ vanishes identically. The coordinate system can be chosen so that A_y = 0. We will use a notation a(x, t) for A_z (x, t).

Calculating the derivatives of a(x, t) with respect to the coordinates x⁺ and x⁻,

$$\begin{equation}u={\partial }_{-}a\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad w={\partial }_{+}a,\end{equation} \tag{ 18 }$$

we find the relationships between u and w and the electric, e_z = −∂_t a, and magnetic, b_y = −∂_x a, field components. They read

$$\begin{equation}u=-\frac{{e}_{z}-{b}_{y}}{\sqrt{2}}\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad w=-\frac{{e}_{z}+{b}_{y}}{\sqrt{2}}.\end{equation} \tag{ 19 }$$

The field variables u and w are not independent of each other. Due to commutation of partials, ∂₋₊ a = ∂₊₋ a, the field variables u and w obey the equation

$$\begin{equation}{\partial }_{-}w={\partial }_{+}u.\end{equation} \tag{ 20 }$$

The Lorentz transformation of the fields u and w

$$\begin{equation}{w}^{\prime }=-\frac{{e}_{z}^{\prime }+{b}_{y}^{\prime }}{\sqrt{2}}=-{\text{e}}^{+\eta }\frac{{e}_{z}+{b}_{y}}{\sqrt{2}}={\text{e}}^{+\eta }w,\end{equation} \tag{ 21 }$$

$$\begin{equation}{u}^{\prime }=-\frac{{e}_{z}^{\prime }-{b}_{y}^{\prime }}{\sqrt{2}}=-{\text{e}}^{-\eta }\frac{{e}_{z}-{b}_{y}}{\sqrt{2}}={\text{e}}^{-\eta }u\end{equation} \tag{ 22 }$$

shows that the w and u are transformed as the − and + components of a contravariant four-vector, respectively. The field product $uw=({b}_{y}^{2}-{e}_{z}^{2})/2$,

$$\begin{equation}{u}^{\prime }{w}^{\prime }=uw,\end{equation} \tag{ 23 }$$

is Lorentz invariant and is equal to the first Poincaré invariant $\mathfrak{F}$.

In the case of interacting electromagnetic pulses with the same polarization, the Lagrangian given by equation (1) can be written in terms of the field variables u and w, defined by equation (18), as

$$\begin{equation}\mathcal{L}=-\frac{{m}^{4}}{4\pi \alpha }[wu+\alpha Q(wu)]=-\frac{{m}^{4}}{4\pi \alpha }\bar{\mathcal{L}}(u,w).\end{equation} \tag{ 24 }$$

Here

$$\begin{equation}\bar{\mathcal{L}}=wu+\alpha Q(wu)\end{equation} \tag{ 25 }$$

is the normalized Lagrangian, where it is assumed that uw ≪ 1 and

$$\begin{equation}Q(uw)=-\frac{1}{2\pi }\mathcal{M}\left(0,\sqrt{2uw}\right).\end{equation} \tag{ 26 }$$

The function Q(ζ) can be represented in the form of the series

$$\begin{equation}\quad Q(\zeta )=\sum\limits _{m=2}^{\infty }{b}_{m}{\zeta }^{m}\end{equation} \tag{ 27 }$$

with coefficients b_m

$$\begin{equation}{b}_{m}=\frac{{2}^{3(m-1)}{B}_{2m}}{\pi m(2m-1)(m-1)}\end{equation} \tag{ 28 }$$

proportional to the Bernoulli numbers [40], B_n . These coefficients can be obtained from the general expansion of $\mathcal{M}(\mathfrak{e},\mathfrak{b})$ in series of powers of $\mathfrak{e}$, and $\mathfrak{b}$, (see reference [41]) setting $\mathfrak{e}=0$ and $\mathfrak{b}=\sqrt{2\zeta }$, where ζ is the argument of the function Q(ζ) (see equations (26) and (27)). It represents actually the field invariant $\mathfrak{F}$ in equation (8). Several leading order coefficients b_m in the expansion (27) are presented in table 1. Using the asymptotic dependence of the Bernoulli numbers at m → ∞ (see reference [40] (section 24.11))

$$\begin{equation}{B}_{2m}\simeq {(-1)}^{m+1}4\sqrt{\pi m}{\left(\frac{m}{\pi e}\right)}^{2m},\end{equation} \tag{ 29 }$$

we obtain the asymptotic expression for the coefficients b_m ,

$$\begin{equation}{b}_{m}\simeq \frac{{(-1)}^{m+1}}{4\sqrt{\pi {m}^{5}}}{\left(\frac{{2}^{3/2}m}{\pi e}\right)}^{2m}.\end{equation} \tag{ 30 }$$

We see that the coefficients of the series for Q(ζ) in equation (27) grow faster than exponentially for large m. Thus this series can be considered only as an asymptotic series and should be truncated by taking into account only a finite number m₀ of its terms. Nevertheless, this finite number m₀ is rather large for the relevant values of |uw| ≪ 1. The truncation number m₀ can be estimated as a maximum of m for which the terms in equation (27) are still decreasing. This criterion gives

$$\begin{equation}{m}_{0}\sim \frac{\pi e}{\sqrt{8\vert uw\vert }}.\end{equation} \tag{ 31 }$$

Thus the condition m ≲ m₀ is very weak allowing the use of a substantially large number of the terms in equation (27). This problem was considered from different points of view also in references [41, 42].

Table 1. Several coefficients in the series of Q(ζ), equation (27).

m	bm	m	bm
2	$\quad -\frac{2}{45\pi }$	6	$\quad -\frac{5660\enspace 672}{225\enspace 225\pi }$
3	$\quad \frac{16}{315\pi }$	7	$\quad \frac{65\enspace 536}{117\pi }$
4	$\quad -\frac{64}{315\pi }$	8	$\quad -\frac{474\enspace 087\enspace 424}{26\enspace 775\pi }$
5	$\quad \frac{512}{297\pi }$	9	$\quad \frac{45\enspace 997\enspace 883\enspace 392}{61\enspace 047\pi }$

The Euler–Lagrange equation corresponding to the Lagrangian (24) and the equation (20) gives the field equations. They can be written as a current conservation equation

$$\begin{equation}{\partial }_{+}\left(\frac{\partial \bar{\mathcal{L}}}{\partial w}\right)+{\partial }_{-}\left(\frac{\partial \bar{\mathcal{L}}}{\partial u}\right)=0.\end{equation} \tag{ 32 }$$

We can define the canonical momenta

$$\begin{equation}{\bar{{\Pi}}}^{+}=\frac{\partial \bar{\mathcal{L}}}{\partial w}=u[1+\alpha {Q}^{\prime }(wu)],\end{equation} \tag{ 33 }$$

$$\begin{equation}{\bar{{\Pi}}}^{-}=\frac{\partial \bar{\mathcal{L}}}{\partial u}=w[1+\alpha {Q}^{\prime }(wu)],\end{equation} \tag{ 34 }$$

where Q'(ζ) denotes differentiation of Q(ζ) with respect to its single argument. Then equation (32) takes the form

$$\begin{equation}{\partial }_{+}{\bar{{\Pi}}}^{+}+{\partial }_{-}{\bar{{\Pi}}}^{-}=0.\end{equation} \tag{ 35 }$$

Using the above obtained relationships we can find the explicit form of the field equation. They can be written as

$$\begin{equation}{\partial }_{+}[u+u\alpha {Q}^{\prime }(wu)]+{\partial }_{-}[w+w\alpha {Q}^{\prime }(wu)]=0,\end{equation} \tag{ 36 }$$

$$\begin{equation}{\partial }_{-}w={\partial }_{+}u.\end{equation} \tag{ 37 }$$

We can use equation (37) to exclude either ∂₋ w or ∂₊ u from equation (36). As a result we obtain a set of two equivalent symmetric field equations:

$$\begin{equation}{\partial }_{+}\left[u+\frac{\alpha }{2}u{Q}^{\prime }(uw)\right]=-\frac{\alpha }{2}{\partial }_{-}\left[w{Q}^{\prime }(uw)\right],\end{equation} \tag{ 38 }$$

$$\begin{equation}{\partial }_{-}\left[w+\frac{\alpha }{2}w{Q}^{\prime }(uw)\right]=-\frac{\alpha }{2}{\partial }_{+}\left[u{Q}^{\prime }(uw)\right].\end{equation} \tag{ 39 }$$

This form of the field equations is convenient for applying a perturbation approach because in the region where w = 0 equation (38) takes the form

$$\begin{equation}{\partial }_{+}u=0,\quad \text{hence}\enspace u=u({x}^{-}),\end{equation} \tag{ 40 }$$

whereas in the region where u = 0 equation (39) becomes

$$\begin{equation}{\partial }_{-}w=0,\quad \text{hence}\enspace w=w({x}^{+}).\end{equation} \tag{ 41 }$$

The momenta introduced above can be represented in the same regions as

$$\begin{equation}{\bar{{\Pi}}}^{+}=u=u({x}^{-}),\qquad {\bar{{\Pi}}}^{-}=0,\end{equation} \tag{ 42 }$$

$$\begin{equation}{\bar{{\Pi}}}^{-}=w=w({x}^{+}),\qquad {\bar{{\Pi}}}^{+}=0,\end{equation} \tag{ 43 }$$

respectively.

It is also useful to introduce the energy momentum tensor

$$\begin{equation}{\bar{T}}_{\lambda }^{\mu }={a}_{\lambda }\frac{\partial \bar{\mathcal{L}}}{\partial {a}_{\mu }}-{\delta }_{\lambda }^{\mu }\bar{\mathcal{L}}\quad (\lambda \enspace ,\mu =\pm ),\end{equation} \tag{ 44 }$$

where a_λ = ∂_λ a and δ_μν is the Kronecker delta. The energy momentum tensor is useful to express the conservation of energy and momentum of the field

$$\begin{equation}{\partial }_{\mu }{T}_{\lambda }^{\mu }=0.\end{equation} \tag{ 45 }$$

In terms of the variables u and w the tensor ${\bar{T}}_{\lambda }^{\mu }$ components can be written as

$$\begin{equation}{\bar{T}}_{+}^{+}=w\frac{\partial \bar{\mathcal{L}}}{\partial w}-\bar{\mathcal{L}}=\alpha \left[uw{Q}^{\prime }(uw)-Q(uw)\right],\end{equation} \tag{ 46 }$$

$$\begin{equation}{\bar{T}}_{-}^{+}=u\frac{\partial \bar{\mathcal{L}}}{\partial w}={u}^{2}\left[1+\alpha {Q}^{\prime }(uw)\right],\end{equation} \tag{ 47 }$$

$$\begin{equation}{\bar{T}}_{+}^{-}=w\frac{\partial \bar{\mathcal{L}}}{\partial u}={w}^{2}\left[1+\alpha {Q}^{\prime }(uw)\right],\end{equation} \tag{ 48 }$$

$$\begin{equation}{\bar{T}}_{-}^{-}=u\frac{\partial \bar{\mathcal{L}}}{\partial u}-\bar{\mathcal{L}}=\alpha \left[uw{Q}^{\prime }(uw)-Q(uw)\right].\end{equation} \tag{ 49 }$$

We see that the diagonal components, ${\bar{T}}_{+}^{+}={\bar{T}}_{-}^{-}$, are equal to each other due the tensor symmetry T_tx = T_xt , see also reference [44].

We formulate the scattering problem for the field equations (36) and (37) or (38) and (39).

In the case of non-interacting waves, the functions u and w do not change since they are equal to u = u(x⁻) and w = w(x⁺). Let us consider the case where the functions u = u(x⁻) and w = w(x⁺) describe the two pulses of finite length. The length is L₊ for the pulse propagating in the x⁺ direction and it equals L₋ in the case of the pulse propagating in the in the x⁻ direction. We assume that the amplitude of the u(x⁻) function is approximately equal to u₀ at |x⁻| ≲ L₋, and it becomes exponentially small as compared with u₀ outside this region. We assume also that there are no spatial structures in the u(x⁻) pulse with scale-length much less than L₋. We set analogous assumptions for the w(x⁺) function with the obvious substitution: − → +. As was noticed above, the pulse amplitudes are substantially small, u₀ w₀ ≪ 1. We will consider the case where the nonlinear effects, leading to harmonics generation of different orders, can be treated perturbatively. Thus, we should exclude too small and too large values of the ratio L₋/L₊. This point will be considered below. Here we present only a simplified condition for the applicability of the perturbation approach used below:

$$\begin{equation}\alpha {w}_{0}^{2}\ll \frac{{L}_{-}}{{L}_{+}}\ll \frac{1}{\alpha {u}_{0}^{2}}.\end{equation} \tag{ 50 }$$

This strong inequality is Lorentz-invariant, that means that all its terms have the same dependence on the boost parameter η given by equation (15). This criterion will be made more precise below when we will further specify the parameters of the pulses.

Let us now take into account the interaction between the waves u and w. It takes place only in the following 'rectangular' domain, Ω, on the (x⁺, x⁻)-plane, where both the following inequalities are valid:

$$\begin{equation}\vert {x}^{-}\vert \lesssim {L}_{-},\quad \text{and}\quad \vert {x}^{+}\vert \lesssim {L}_{+}.\end{equation} \tag{ 51 }$$

In this 'rectangular' region Ω, both amplitudes u and w depend on both the independent coordinates: x⁺ and x⁻. Outside this region, when |x⁻| ≫ L₋ or |x⁺| ≫ L₊, the waves do not interact with each other, and w depends only on the coordinate x⁺, whereas u depends only on the coordinate x⁻. We may distinguish incoming and outgoing waves in the latter combined region:

$$\begin{equation}w(\vert {x}^{+}\vert \lesssim {L}_{+}\enspace ,{x}^{-}\gg {L}_{-})={w}^{\text{in}}({x}^{+}),\end{equation} \tag{ 52 }$$

$$\begin{equation}u(-{x}^{+}\gg {L}_{+}\enspace ,\vert {x}^{-}\vert \lesssim {L}_{-})={u}^{\text{in}}({x}^{-}),\end{equation} \tag{ 53 }$$

$$\begin{equation}w(\vert {x}^{+}\vert \lesssim {L}_{+}\enspace ,-{x}^{-}\gg {L}_{-})={w}^{\text{out}}({x}^{+}),\end{equation} \tag{ 54 }$$

$$\begin{equation}u({x}^{+}\gg {L}_{+}\enspace ,\vert {x}^{-}\vert \lesssim {L}_{-})={u}^{\text{out}}({x}^{-}).\end{equation} \tag{ 55 }$$

Outside all these regions, including the Ω region, either the field u or the field w, or both are exponentially small or vanish.

The scattering problem can be formulated now as follows: to find the outgoing waves w^out(x⁺) and u^out(x⁻), knowing the incoming waves wⁱⁿ(x⁺) and uⁱⁿ(x⁻).

5.1. Perturbative solution of 1st order

Equations (38) and (39) can be rewritten in the following integral form that takes into account the initial conditions (52) and (53):

$$\begin{equation}u+\frac{\alpha }{2}u{Q}^{\prime }(uw)={u}^{\text{in}}-\frac{\alpha }{2}\underset{-\infty }{\overset{{x}^{+}}{\int }}{\partial }_{-}[w{Q}^{\prime }(uw)]\enspace \mathrm{d}{x}^{+},\end{equation} \tag{ 56 }$$

$$\begin{equation}w+\frac{\alpha }{2}w{Q}^{\prime }(uw)={w}^{\text{in}}-\frac{\alpha }{2}\underset{+\infty }{\overset{{x}^{-}}{\int }}{\partial }_{+}[u{Q}^{\prime }(uw)]\enspace \mathrm{d}{x}^{-}.\end{equation} \tag{ 57 }$$

In view of equations (54) and (55) we can formally send x^± to ±∞ in equations (56) and (57) and obtain the following scattering relationships:

$$\begin{equation}{u}^{\text{out}}={u}^{\text{in}}-\frac{\alpha }{2}\underset{-\infty }{\overset{+\infty }{\int }}{\partial }_{-}[w{Q}^{\prime }(uw)]\enspace \mathrm{d}{x}^{+},\end{equation} \tag{ 58 }$$

$$\begin{equation}{w}^{\text{out}}={w}^{\text{in}}+\frac{\alpha }{2}\underset{-\infty }{\overset{\infty }{\int }}{\partial }_{+}[u{Q}^{\prime }(uw)]\enspace \mathrm{d}{x}^{-}.\end{equation} \tag{ 59 }$$

For sufficiently small uw we are able to treat equations (58) and (59) perturbatively, and to substitute into the integrands uⁱⁿ and wⁱⁿ instead of u and w, respectively. Then we obtain:

$$\begin{equation}{u}^{\text{out}}={u}^{\text{in}}-\frac{\alpha }{2}\left[\underset{-\infty }{\overset{+\infty }{\int }}{({w}^{\text{in}})}^{2}{Q}^{{\prime\prime}}\left({u}^{\text{in}}{w}^{\text{in}}\right)\mathrm{d}{x}^{+}\right]{\partial }_{-}{u}^{\text{in}},\end{equation} \tag{ 60 }$$

$$\begin{equation}{w}^{\text{out}}={w}^{\text{in}}+\frac{\alpha }{2}\left[\underset{-\infty }{\overset{\infty }{\int }}{({u}^{\text{in}})}^{2}{Q}^{{\prime\prime}}\left({u}^{\text{in}}{w}^{\text{in}}\right)\mathrm{d}{x}^{-}\right]{\partial }_{+}{w}^{\text{in}}.\end{equation} \tag{ 61 }$$

For this perturbation procedure to be relevant, the second terms in the right hand sides of these equations must be much less than the first terms. If we take for this estimation only leading term in Q(ζ → 0) ∝ ζ², then we obtain the strong inequality (50).

The leading order of Q(ζ) ∝ ζ², caused by four-photon interaction, leads to the following delays

$$\begin{equation}\frac{\alpha }{2}{Q}^{{\prime\prime}}\left(0\right)\underset{-\infty }{\overset{+\infty }{\int }}{({w}^{\text{in}})}^{2}\enspace \mathrm{d}{x}^{+},\end{equation} \tag{ 62 }$$

and

$$\begin{equation}\frac{\alpha }{2}{Q}^{{\prime\prime}}\left(0\right)\underset{-\infty }{\overset{\infty }{\int }}{({u}^{\text{in}})}^{2}\enspace \mathrm{d}{x}^{-}.\end{equation} \tag{ 63 }$$

of the u^out and w^out pulses and with respect to uⁱⁿ and wⁱⁿ, respectively.

This corresponds to the phase shift between the interacting waves discussed in references [44, 46]. The strong inequality (50) ensures that the delays are much shorter than L₋ and L₊, respectively. These delays do not correspond to harmonic generation, although they determine the applicability of our approach. The next, subleading orders of Q(ζ → 0) cause generation of the harmonics.

Our approach uses the following set of small parameters: α, u₀ w₀, $\alpha {w}_{0}^{2}{L}_{+}/{L}_{-}$, and $\alpha {u}_{0}^{2}{L}_{-}/{L}_{+}$. Here we address the question whether the next order corrections to equations (60) and (61), as they would appear after the next iterations of equations (56) and (57), or the subleading terms in Q are more important for harmonic generation. To answer this question we use the structure of $\bar{\mathcal{L}}$, given by equations (24), (26), (27), (56) and (57) as the equations of our perturbation theory starting from zero order: u = uⁱⁿ and w = wⁱⁿ. The combined result of the iterations of equations (56) and (57) and of retaining the higher order terms in the expansion of the function Q(u, w) can be written symbolically as

$$\begin{equation}{u}^{\text{out}}{\left.={u}^{\text{in}}+\sum\limits _{\ell =1}^{\infty }\sum\limits _{m=0}^{\infty }{\hat{U}}_{\ell m}w{(\alpha uw)}^{\ell }{(uw)}^{m}\right\vert }^{\text{in}},\end{equation} \tag{ 64 }$$

$$\begin{equation}{w}^{\text{out}}{\left.={w}^{\text{in}}+\sum\limits _{\ell =1}^{\infty }\sum\limits _{m=0}^{\infty }{\hat{W}}_{\ell m}u{(\alpha uw)}^{\ell }{(uw)}^{m}\right\vert }^{\text{in}},\end{equation} \tag{ 65 }$$

where ${\hat{U}}_{\ell m}$ and ${\hat{W}}_{\ell m}$ are linear integro-differential operators that do not contain any small or large parameters. They act implicitly on all arguments x^± to the right of them separately. Here we use the expansion (27) of Q(ζ).

Due to the conservation of the kinematic momentum, the harmonic with the number n of the wave u will be generated by the terms proportional to uⁿ , uⁿ⁺², etc. The largest term of this type corresponds to the term in equation (64) with (ℓ, m) = (1, n − 1). All other terms contributing to the same n-harmonic will have higher orders of α or/and higher orders of uw. Thus the leading order contribution to the amplitude of the n-harmonic comes from the term of the form w(αuw)uⁿ⁻¹ wⁿ⁻¹ in agreement with equation (60). The same conclusion can be drawn for the w-wave.

We may conclude that the higher order terms in the expansion of the function Q(uw) determine the leading order contribution to harmonic generation, and hence the expressions (60) and (61) are sufficient to calculate the intensities of the harmonics in the leading order relatively to the small parameters.

6.1. General considerations

We consider here two plane waves of the same polarization, when the incident spectrum corresponding to the wave uⁱⁿ is narrow enough. In this case we may choose uⁱⁿ in the form:

$$\begin{equation}{u}^{\text{in}}({x}^{-})=\text{Re}\left\{{U}_{0}({x}^{-}){\text{e}}^{\text{i}{\omega }_{-}{x}^{-}}\right\},\end{equation} \tag{ 66 }$$

where U₀ is a smooth function of its argument with maximum value ${\bar{U}}_{0}$, and decaying at least exponentially outside the region |x⁻| ≲ L₋. The carrier frequency of the wave is ω₋. We assume that ω₋ L₋ is significantly larger than one. We assume that U₀ is almost constant at the spatial-time scales of the order of 1/ω₋. The spectral relative half width of the wave is about

$$\begin{equation}\frac{{\Delta}{\omega }_{-}}{{\omega }_{-}}\sim \frac{1}{{\omega }_{-}{L}_{-}}.\end{equation} \tag{ 67 }$$

Instantaneous intensity of this wave (averaged over its carrier period) is equal to

$$\begin{equation}{I}_{1}=\frac{{m}^{4}}{8\pi \alpha }{U}_{0}{U}_{0}^{{\ast}}.\end{equation} \tag{ 68 }$$

Here the asterisk '*' denotes the complex conjugation. This expression is in accordance to the expressions for the energy momentum tensor (46)–(49) in terms of u and w. We will skip analogous comments below.

It is obvious that the spectral width of the n-harmonic, ${\Delta}{\omega }_{-}^{(n)}$ grows with n. We calculate ${\Delta}{\omega }_{-}^{(n)}$ below in this section. We are able to separate the harmonics with adjacent numbers, when ${\Delta}{\omega }_{-}^{(n)}$ is less than, for example, half spectral distance between the harmonics, ω₋:

$$\begin{equation}2{\Delta}{\omega }_{-}^{(n)}< {\omega }_{-}.\end{equation} \tag{ 69 }$$

We will not be interested here in the detailed form of the harmonic spectral line for separated harmonics, but only in the total intensity in the harmonic of order n, integrated over its spectral line.

As we explain above, the leading term of the n-harmonic wave field, ${u}_{(n)}^{\text{out}}({x}^{-})$, is determined by the following expression, obtained from equation (58):

$$\begin{equation}-\alpha (n+1){b}_{n+1}{L}_{+}{\partial }_{-}{({u}^{\text{in}})}^{n}\left\langle {({w}^{\text{in}})}^{n+1}\right\rangle ,\end{equation} \tag{ 70 }$$

where

$$\begin{equation}\left\langle {({w}^{\text{in}})}^{n+1}\right\rangle =\frac{1}{2{L}_{+}}\underset{-\infty }{\overset{\infty }{\int }}{\left[{w}^{\text{in}}\left({x}^{+}\right)\right]}^{n+1}\enspace \mathrm{d}{x}^{+},\end{equation} \tag{ 71 }$$

and

$$\begin{align}{\partial }_{-}{({u}^{\text{in}})}^{n}& ={\partial }_{-}{\left[\frac{1}{2}\left({U}_{0}{\text{e}}^{\text{i}{\omega }_{-}{x}^{-}}+{U}_{0}^{{\ast}}{\text{e}}^{-\text{i}{\omega }_{-}{x}^{-}}\right)\right]}^{n}\\ & =\frac{1}{{2}^{n-1}}{\partial }_{-}\enspace \mathrm{R}\mathrm{e}\left({U}_{0}^{n}{\text{e}}^{\text{in}{\omega }_{-}{x}^{-}}\right)+\dots \enspace .\end{align} \tag{ 72 }$$

Only the first term in equation (72) contributes to ${u}_{(n)}^{\text{out}}({x}^{-})$, the other terms contribute to the n − 2-, n − 4- and lower harmonics. Hence we obtain:

$$\begin{equation}{u}_{(n)}^{\text{out}}({x}^{-})=\alpha \frac{n(n+1){b}_{n+1}}{{2}^{n-1}}{\omega }^{-}{L}_{+}\left\langle {({w}^{\text{in}})}^{n+1}\right\rangle \mathrm{I}\mathrm{m}\left({U}_{0}^{n}{\text{e}}^{\text{i}n{\omega }_{-}{x}^{-}}\right).\end{equation} \tag{ 73 }$$

Then the leading term of the 'instantaneous' intensity ${I}_{n}^{\text{out}}({x}^{-})$ of the n-harmonic averaged over its carrier period can be expressed as:

$$\begin{equation}{I}_{n}^{\text{out}}({x}^{-})=\frac{{m}^{4}\alpha }{8\pi }\frac{{n}^{2}{(n+1)}^{2}{b}_{n+1}^{2}}{{2}^{2(n-1)}}{\left({\omega }_{-}{L}_{+}\right)}^{2}{\left\langle {({w}^{\text{in}})}^{n+1}\right\rangle }^{2}\vert {U}_{0}({x}^{-}){\vert }^{2n}.\end{equation} \tag{ 74 }$$

We recall that this intensity is integrated over the spectral line of the n-harmonic. The normalized intensity becomes

$$\begin{equation}\frac{{I}_{n}^{\text{out}}}{{I}_{1}^{\text{in}}}={\alpha }^{2}\frac{{n}^{2}{(n+1)}^{2}{b}_{n+1}^{2}}{{2}^{2(n-1)}}{\left({\omega }_{-}{L}_{+}\right)}^{2}{\left\langle {({w}^{\text{in}})}^{n+1}\right\rangle }^{2}\vert {U}_{0}({x}^{-}){\vert }^{2n-2}.\end{equation} \tag{ 75 }$$

This expression is invariant relative to boosting. We may rewrite the latter expression by order of magnitude:

$$\begin{equation}\frac{{I}_{n}^{\text{out}}}{{I}_{1}^{\text{in}}}\propto {\left({\omega }_{-}{L}_{-}\right)}^{2}{\left[\alpha {w}_{0}^{2}\frac{{L}_{+}}{{L}_{-}}\right]}^{2}{({w}_{0}{u}_{0})}^{2n-2},\end{equation} \tag{ 76 }$$

where a coefficient depending only on n has been dropped. The last two terms in equation (76) are small parameters ensuring the validity of our approach. It breaks completely even on a qualitative level, when $\alpha {w}_{0}^{3}{u}_{0}{L}_{+}{\omega }_{-}$ becomes approximately equal to unity. In this case, secular effects in evolution of the wave shape become very important. The theory presented in references [30, 31] treated the case when the latter parameter is approximately equal to or much higher than unity.

Equation (74) is written for a general profile of wⁱⁿ(x⁺) that does not contain any additional spatial parameter besides L₊. However, when this profile has the form similar to equation (66) we find that

$$\begin{equation}{w}^{\text{in}}({x}^{+})=\mathrm{R}\mathrm{e}\left\{{W}_{0}({x}^{+}){\text{e}}^{\text{i}{\omega }_{+}{x}^{+}}\right\},\end{equation} \tag{ 77 }$$

with ω₊ L₊ ≫ 1, then $\left\langle {({w}^{\text{in}})}^{n+1}\right\rangle$ becomes considerably less than ${W}_{0}^{n+1}$ for even values of n. It means that in this case harmonics with even numbers will be significantly suppressed with respect to the estimation (76).

The spectrum of the n-harmonic is determined by the multiplier ${U}_{0}{({x}^{-})}^{n}\mathrm{exp}(\text{i}n{\omega }_{-}{x}^{-})$ in the expression (73). If we approximate U₀(x⁻) in the region close to its maximum as

$$\begin{equation}{U}_{0}({x}^{-})\approx {\bar{U}}_{0}\left[1-\frac{1}{2}{\left(\frac{{x}^{-}}{{L}_{-}}\right)}^{2}\right],\end{equation} \tag{ 78 }$$

where ${\bar{U}}_{0}$ is a constant equal to an amplitude of the uⁱⁿ-wave, then for sufficiently high n the spectrum of the n-harmonic will be determined by the following x⁻ dependent factor:

$$\begin{equation}\mathrm{exp}\left[-\frac{n}{2}{\left(\frac{{x}^{-}}{{L}_{-}}\right)}^{2}+\text{i}n{\omega }_{-}{x}^{-}\right].\end{equation} \tag{ 79 }$$

Its Fourier transform is proportional to

$$\begin{equation}\mathrm{exp}\left[-\frac{{(\omega -n{\omega }_{-})}^{2}{L}^{-\enspace 2}}{2n}\right].\end{equation} \tag{ 80 }$$

Hence the spectral width is

$$\begin{equation}{\Delta}{\omega }_{-}^{(n)}\simeq \frac{\sqrt{n}}{{L}_{-}}.\end{equation} \tag{ 81 }$$

In this case, the condition (69) for harmonic spectral separation is fulfilled, when

$$\begin{equation}n< {n}_{\mathrm{max}}\sim \frac{{\left({\omega }_{-}{L}_{-}\right)}^{2}}{4}.\end{equation} \tag{ 82 }$$

The latter inequality is non-restrictive. It allows us to skip consideration of the opposite case, when n is considerably is above n_max.

6.2. Formulas for the n-harmonic generation in physical units

We use in this section the presentations (66) and (77) for the in-waves, with the profiles given by

$$\begin{equation}{U}_{0}({x}^{-})={\bar{U}}_{0}\enspace \mathrm{exp}\left[-\frac{{({x}^{-})}^{2}}{2{({L}_{-})}^{2}}\right],\end{equation} \tag{ 83 }$$

$$\begin{equation}{W}_{0}({x}^{+})={\bar{W}}_{0}\enspace \mathrm{exp}\left[-\frac{{({x}^{+})}^{2}}{2{({L}_{+})}^{2}}\right].\end{equation} \tag{ 84 }$$

We introduce a relationship between the intensities of two incoming waves, averaged over carrier periods, and the parameters |U₀| and |W₀| as follows:

$$\begin{equation}\vert {\bar{U}}_{0}{\vert }^{2}=\frac{{I}_{u}}{5.8095\times 1{0}^{28}{\text{W}\;\text{cm}}^{-2}},\end{equation} \tag{ 85 }$$

$$\begin{equation}\vert {\bar{W}}_{0}{\vert }^{2}=\frac{{I}_{w}}{5.8095\times 1{0}^{28}{\text{Wcm}}^{-2}},\end{equation} \tag{ 86 }$$

where I_u and I_w are maximal intensities of the u and w in-waves, respectively.

The results described by equations (74) and (75) contain the dimensionless parameter $\left\langle {({w}^{\text{in}})}^{n+1}\right\rangle$ that depends on the w-pulse profile. We present the parameters below in the particular case where W₀ is a real number:

$$\begin{equation}\left\langle {({w}^{\text{in}})}^{n+1}\right\rangle ={W}_{0}^{n+1}{g}_{n}({\omega }_{+}{L}_{+}).\end{equation} \tag{ 87 }$$

We do not consider here phase envelope effects that are important for a short enough w-pulse. We set actually a certain phase to have maximal yields for the even harmonics of the u-wave. The coefficient g_n (ω₊ L₊) in equation (87) depends only on the harmonic number and on the parameter ω₊ L₊ describing the in-pulse of the w-wave. The plots of g_n versus ω₊ L₊ are shown in figures 1 and 2 separately for odd and even n. For even n the coefficient g_n becomes equal to zero at ω₊ L₊ → ∞, whereas for odd n they tend in this limit to constant values:

$$\begin{equation}{g}_{2\ell -1}(\infty )=\frac{\sqrt{\pi }\enspace (2\ell )!}{{2}^{2\ell +1}\sqrt{\ell }{(\ell !)}^{2}}.\end{equation} \tag{ 88 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The ratio of the energy of the n-harmonic, ${\mathcal{E}}_{n}^{\text{out}}$, in the outgoing u-wave to the initial incoming energy of this wave is equal to

$$\begin{equation}\frac{{\mathcal{E}}_{n}^{\text{out}}}{{\mathcal{E}}_{1}^{\text{in}}}={f}_{n}{g}_{n}^{2}\enspace {\left[\alpha {\omega }_{-}{L}_{+}\vert {\bar{W}}_{0}{\vert }^{2}\right]}^{2}{\left\vert {\bar{U}}_{0}{\bar{W}}_{0}\right\vert }^{2(n-1)}.\end{equation} \tag{ 89 }$$

where

$$\begin{equation}{f}_{n}=\frac{{n}^{2}{(n+1)}^{2}{b}_{n+1}^{2}}{{2}^{2n-2}\sqrt{n}}.\end{equation} \tag{ 90 }$$

Table 2 presents the values of f_n and g_n (∞) for several n.

Table 2. Coefficients f_n and g_n (∞) in equation (89).

n	fn	gn (∞)	n	fn	gn (∞)
2	0.06568	0	9	2.842 × 1014	0.08162
3	0.8580	0.1599	10	4.837 × 1017	0
4	37.15	0	11	1.203 × 1021	0.07017
5	3.973 × 103	0.1212	12	4.227 × 1024	0
6	8.826 × 105	0	13	2.043 × 1028	0.06153
7	3.629 × 108	0.09753	14	1.326 × 1032	0
8	2.540 × 1011	0	15	1.133 × 1036	0.05479

The dimensionless and Lorentz invariant factor in the brackets in equation (89) is one of the small parameters of the theory. It is used instead of the simplified form $\alpha {w}_{0}^{2}{L}_{+}/{L}_{-}$ introduced above to provide a more accurate description in the limit ω₋ L₋ ≫ 1.

The limit ω₊ L₊ → 0 in equation (89) corresponds to the case when the w-wave becomes actually a constant amplitude cross electromagnetic field. This fact is in accordance with the results of references [30, 31], where the generation of the second harmonic due to six-photon interaction plays a main role in the steepening of nonlinear waves on the constant background of a counter-propagating cross electromagnetic field.

7.1. Definitions and governing equations

We use here the same frame of reference that was introduced in section 3. Taking into account the general expression (11) we may define the fields

$$\begin{equation}{u}_{i}={\partial }_{-}{A}_{i};\qquad {w}_{i}={\partial }_{+}{A}_{i},\end{equation} \tag{ 91 }$$

where, and below in this section, i, j, ... span the set {y, z}. The electric and magnetic field components can be expressed as

$$\begin{equation}{E}_{i}=\frac{{u}_{i}-{w}_{i}}{\sqrt{2}};\qquad {B}_{i}=-{e}_{ij}\frac{{u}_{j}+{w}_{j}}{\sqrt{2}}.\end{equation} \tag{ 92 }$$

Here e_ij is the skew-symmetric two-matrix with e_zy = 1. The field invariants f and g can be written as

$$\begin{equation}f={u}_{i}{w}_{i}=({B}^{2}-{E}^{2})/2;\end{equation} \tag{ 93 }$$

$$\begin{equation}g={e}_{ij}{u}_{i}{w}_{j}=\left(\mathbf{B}\cdot \mathbf{E}\right)/2.\end{equation} \tag{ 94 }$$

The Lagrangian for this problem can be written in a way similar to the case given by equation (24). It reads

$$\begin{equation}\mathcal{L}=-\frac{{m}^{4}}{4\pi \alpha }\left[f+\alpha Q(f,g)\right]=-\frac{{m}_{\text{e}}^{4}}{4\pi \alpha }\bar{\mathcal{L}},\end{equation} \tag{ 95 }$$

where

$$\begin{equation}\bar{\mathcal{L}}=f+\alpha Q(f,g),\end{equation} \tag{ 96 }$$

and

$$\begin{align}Q(f,g)& =-\frac{1}{2\pi }\mathcal{M}\left(\sqrt{\sqrt{{f}^{2}+{g}^{2}}-f},\sqrt{\sqrt{{f}^{2}+{g}^{2}}+f}\right)\\ & =\sum\limits _{m+n\geqslant 2}{b}_{m,\enspace n}\enspace {f}^{m}{g}^{n}.\end{align} \tag{ 97 }$$

Although $\bar{\mathcal{L}}$ and Q are now functions of four field components they can be presented in the form of functions of two arguments which are two field invariants, f and g. We note that the Lagrangian $\bar{\mathcal{L}}$ and Q, introduced in section 4, that describe the interaction of two electromagnetic waves of the same polarization are functions of the two field components u and w, which were combined in one field invariant uw.

The field invariant g is a pseudo-scalar, hence $\mathcal{L}$ should be an even function of it. As a result all coefficients b_m,n in the sum given by equation (97) with odd n vanish. It is possible to obtain explicit analytic expressions for the coefficients b_m,n . However the final expressions are too cumbersome for their presentation in this paper. We will present below their particular values, appearing in the problem for perpendicular polarizations.

In a way similar to that used in section 4, the normalized Lagrangian (96) and (97) gives the equations of the Heisenberg–Euler electrodynamics

$$\begin{equation}{\partial }_{+}\left(\frac{\partial \bar{\mathcal{L}}}{\partial {w}_{i}}\right)+{\partial }_{-}\left(\frac{\partial \bar{\mathcal{L}}}{\partial {u}_{i}}\right)=0.\end{equation} \tag{ 98 }$$

They can be written as

$$\begin{equation}{\partial }_{+}\left[{u}_{i}{\bar{\mathcal{L}}}_{f}-{e}_{ij}{u}_{j}{\bar{\mathcal{L}}}_{g}\right]+{\partial }_{-}\left[{w}_{i}{\bar{\mathcal{L}}}_{f}+{e}_{ij}{w}_{j}{\bar{\mathcal{L}}}_{g}\right]=0.\end{equation} \tag{ 99 }$$

As a result, using equation (96), we obtain

$$\begin{equation}{\partial }_{+}\left[{u}_{i}+\alpha \left({\delta }_{ij}{Q}_{f}-{e}_{ij}{Q}_{g}\right){u}_{j}\right]+{\partial }_{-}\left[{w}_{i}+\alpha \left({\delta }_{ij}{Q}_{f}+{e}_{ij}{Q}_{g}\right){w}_{j}\right]=0.\end{equation} \tag{ 100 }$$

To get a closed system of the field equations in terms of u_i and w_i (instead of A_i ) we may use the following direct consequence of equation (91)

$$\begin{equation}{\partial }_{-}{w}_{i}={\partial }_{+}{u}_{i}.\end{equation} \tag{ 101 }$$

The following system of the field equations, that is equivalent to equations (100) and (101), is more suitable for the perturbation theory in terms of α:

$$\begin{equation}{\partial }_{+}{u}_{i}+\frac{1}{2}\alpha \left\{{\partial }_{+}\left[\left({\delta }_{ij}{Q}_{f}-{e}_{ij}{Q}_{g}\right)\enspace {u}_{j}\right]+{\partial }_{-}\left[\left({\delta }_{ij}{Q}_{f}+{e}_{ij}{Q}_{g}\right)\enspace {w}_{j}\right]\right\}=0,\end{equation} \tag{ 102 }$$

and

$$\begin{equation}{\partial }_{-}{w}_{i}+\frac{1}{2}\alpha \left\{{\partial }_{-}\left[\left({\delta }_{ij}{Q}_{f}+{e}_{ij}{Q}_{g}\right){w}_{j}\right]+{\partial }_{+}\left[\left({\delta }_{ij}{Q}_{f}-{e}_{ij}{Q}_{g}\right){u}_{j}\right]\right\}=0.\end{equation} \tag{ 103 }$$

Integrating equations (102) and (103) along the lines x^∓ = const. respectively, we obtain the following exact relationships between in and out states:

$$\begin{align}{u}_{i}^{\text{out}}({x}^{-})& ={u}_{i}^{\text{in}}({x}^{-})-\frac{\alpha }{2}\underset{-\infty }{\overset{\infty }{\int }}{\partial }_{-}\left[\left({\delta }_{ij}{Q}_{f}+{e}_{ij}{Q}_{g}\right)\enspace {w}_{j}\right]\enspace \mathrm{d}{x}^{+},\end{align} \tag{ 104 }$$

and

$$\begin{equation}{w}_{i}^{\text{out}}({x}^{+})={w}_{i}^{\text{in}}({x}^{+})+\frac{\alpha }{2}\underset{-\infty }{\overset{\infty }{\int }}{\partial }_{+}\left[\left({\delta }_{ij}{Q}_{f}-{e}_{ij}{Q}_{g}\right)\enspace {u}_{j}\right]\enspace \mathrm{d}{x}^{-}.\end{equation} \tag{ 105 }$$

7.2. Perturbative approach

Here we consider the case of small |f| and |g|. We substitute in the expressions on the right hand sides of equations (104) and (105) the unperturbed u and w i.e. uⁱⁿ(x⁻) and wⁱⁿ(x⁺). Assuming strong inequalities similar to those given by equation (50), we obtain for the out-wave

$$\begin{equation}{u}_{i}^{\text{out}}({x}^{-})={u}_{i}^{\text{in}}({x}^{-})-\frac{\alpha }{2}\left({\partial }_{-}{u}_{k}^{\text{in}}\right)\underset{-\infty }{\overset{\infty }{\int }}{w}_{j}^{\text{in}}{w}_{\ell }^{\text{in}}{\left[{\delta }_{ij}{\delta }_{k\ell }{Q}_{ff}+\left({\delta }_{ij}{e}_{k\ell }+{e}_{ij}{\delta }_{k\ell }\right){Q}_{fg}+{e}_{ij}{e}_{k\ell }{Q}_{gg}\right]}^{\text{in}}\enspace \mathrm{d}{x}^{+}.\end{equation} \tag{ 106 }$$

The subscript ⁱⁿ at the bracket means that f and g entering it are calculated with ${u}_{i}^{\text{in}}$ and ${w}_{i}^{\text{in}}$. Here the ${u}_{i}^{\text{in}}$ waves do not depend on the x⁺ coordinate. Similar expressions can be written for ${w}_{i}^{\text{out}}-{w}_{i}^{\text{in}}$ waves.

7.3. Harmonic generation by two counter-propagating waves with perpendicular polarizations

To demonstrate how polarization can affect harmonic generation we consider here the case where the incoming waves have only u_z and w_y components, i.e. where ${u}_{y}^{\text{in}}={w}_{z}^{\text{in}}=0$. In this case equation (106) becomes:

$$\begin{equation}{u}_{z}^{\text{out}}({x}^{-})={u}_{z}^{\text{in}}({x}^{-})-\frac{\alpha }{2}\enspace \left({\partial }_{-}{u}_{z}^{\text{in}}\right)\underset{-\infty }{\overset{\infty }{\int }}{Q}_{gg}(0,{u}_{z}^{\text{in}}{w}_{y}^{\text{in}}){\left({w}_{y}^{\text{in}}\right)}^{2}\enspace \mathrm{d}{x}^{+}.\end{equation} \tag{ 107 }$$

Note that f = 0 and g = u_z w_y exactly in this case. Hence we need to know only the coefficients b_0,n in the expansion (97). They are equal to:

$$\begin{equation}{b}_{0,2\ell }=-\frac{{2}^{4\ell -1}{\Gamma}(4\ell -2)}{\pi }\sum\limits _{\sigma =0}^{2\ell }{(-1)}^{\sigma +1}\frac{{B}_{4\ell -2\sigma }{B}_{2\sigma }}{(4\ell -2\sigma )!\enspace (2\sigma )!}\end{equation} \tag{ 108 }$$

for ℓ = 1, 2, .... These coefficients can be obtained from the general expansion [41] of $\mathcal{M}(\mathfrak{e},\mathfrak{b})$ in series of powers of $\mathfrak{e}$, and $\mathfrak{b}$, setting $\mathfrak{e}=\mathfrak{b}=g$. Several leading coefficients b_0,n are presented in table 3. We have at ℓ → ∞:

$$\begin{equation}{b}_{0,2\ell \gg 1}\simeq -\frac{8}{{(4\ell )}^{5/2}\sqrt{\pi }}{\left(\frac{4\ell }{e\pi }\right)}^{4\ell }.\end{equation} \tag{ 109 }$$

Table 3. Coefficients b_0,n .

	n
	2	4	6	8
b0,n	$-\frac{7}{90\pi }$	$-\frac{76}{945\pi }$	$-\frac{185\enspace 984}{75\enspace 075\pi }$	$-\frac{7007\enspace 744}{16\enspace 065\pi }$

We see that expression (107) is quite similar to what we have in section 5.1, see for comparison equation (60). The differences are only in the subscripts y and z and in the appearance of this different coefficients in the present expansion of Q_gg in comparison with the coefficients in the expansion of Q_ff (in our present terminology). These correspondences allow us to use our previous expressions (75) and (89) to obtain new expressions for the intensities of the high order harmonics for the present combinations of polarizations.

In this way we obtain the following results. The normalized intensity of the n-harmonic (integrated over the separate spectral line of the harmonic) becomes (n ⩾ 3):

$$\begin{equation}\frac{{I}_{z,n}^{\text{out}}}{{I}_{z,1}^{\text{in}}}={\alpha }^{2}\frac{{n}^{2}{(n+1)}^{2}{b}_{0,n+1}^{2}}{{2}^{2(n-1)}}{\left({\omega }_{-}{L}_{+}\right)}^{2}{\left\langle {({w}_{y}^{\text{in}})}^{n+1}\right\rangle }^{2}\vert {U}_{z,0}({x}^{-}){\vert }^{2n-2}.\end{equation} \tag{ 110 }$$

Here the subscripts y and z indicate the direction of the electric fields of the u- and w-waves respectively. The subscripts 1 and n indicate the number of the harmonic in incoming and outgoing u_z -waves. The parameters ω₋ and U_{z, 0}(x⁻) describe the incoming u_z -wave in accordance with definition (66). We imply the obvious insertion of the subscript z into that expression. The parameters L₊ and $\left\langle {({w}_{y}^{\text{in}})}^{n+1}\right\rangle$, describing the incoming w_y -wave are defined by equation (71) with the obvious insertion of the y subscript. Here, L₊ is the finite length of the w_y -wave in the x⁺ direction.

For the total relative energy of the n-harmonic in the outgoing u_z -wave we now have

$$\begin{equation}\frac{{\mathcal{E}}_{z,n}^{\text{out}}}{{\mathcal{E}}_{z,1}^{\text{in}}}={\bar{f}}_{n}\enspace {g}_{z,\enspace n}^{2}\enspace {\left[\alpha {\omega }_{-}{L}_{+}\vert {\bar{W}}_{y,\enspace 0}{\vert }^{2}\right]}^{2}{\left\vert {\bar{U}}_{z,\enspace 0}{\bar{W}}_{y,\enspace 0}\right\vert }^{2(n-1)}.\end{equation} \tag{ 111 }$$

where

$$\begin{equation}{\bar{f}}_{n}=\frac{{n}^{2}{(n+1)}^{2}{b}_{0,\enspace n+1}^{2}}{{2}^{2n-2}\sqrt{n}},\end{equation} \tag{ 112 }$$

and ${\bar{U}}_{z,\enspace 0}$, ${\bar{W}}_{y,\enspace 0}$, and g_{z, n} (ω₊ L₊) are defined analogously to equations (83), (84) and (87). Plots for the functions g_z,n (ω₊ L₊) are presented in figure 1, whereas table 4 presents the values of ${\bar{f}}_{n}$ for several n. The limits g_z,n (∞) are presented in table 2. We may recall here that only the odd harmonics is generated for such a combination of polarizations at least in this approximation, whereas the suppression of the even harmonics for parallel polarization was only approximate even in our present order. We see also that the numerical coefficients f_n are considerably larger than the numerical coefficient ${\bar{f}}_{n}$ with the same n: compare tables 2 and 4. It means that when the polarizations of the two waves are parallel the high harmonic generation, is much more effective than in the case of perpendicular polarizations.

Table 4. Coefficients ${\bar{f}}_{n}$ in equation (111).

n	${\bar{f}}_{n}$	n	${\bar{f}}_{n}$
3	0.005 321	11	8.476 × 1013
5	1.331	13	8.802 × 1019
7	7.061 × 103	15	2.999 × 1026
9	3.305 × 108	17	2.871 × 1033

We analyse the high order harmonics generation during the interaction in vacuum of two strong plane electromagnetic waves. The harmonics generation occurs inside the region of the waves intersection. The results obtained are presented in the frame of reference where the waves are counter-propagating.

The theory presented in our paper takes into account the nonlinear polarization of the QED vacuum caused by the existence of a sea of virtual electron–positron pairs. The theory is developed within the framework of the Heisenberg–Euler electrodynamics implying relatively low amplitude and long-wavelength interacting waves.

The regime of high order harmonics generation is described by using the perturbation theory giving the expressions for the high order harmonics to leading order in terms of the small parameters. We show that the leading contributions to the number n harmonic are obtained by expanding the Heisenberg–Euler Lagrangian in power series of the electromagnetic field amplitudes and not by iterations of the nonlinear contribution to the Heisenberg–Euler Lagrangian. In terms of the QED theory we can say that the amplitude of the nth harmonic is determined by the 2(n + 1)-photon Feynman diagram with one electron loop, but not by the combination of several diagrams with fewer photon lines attached to the electron loops. Although this statement is based on the consideration of two interacting plane waves we may expect that it is valid in a general case. To see this, one can evaluate the impact of nonlinear currents entering the electromagnetic field equations of the Heisenberg–Euler electrodynamics, using only unperturbed electric Eⁱⁿ(x, t) and magnetic Bⁱⁿ(x, t) fields. As a result, the harmonic intensities can be calculated as the electromagnetic radiation emitted by these nonlinear currents. The unperturbed Eⁱⁿ(x, t) and Bⁱⁿ(x, t) fields can be obtained by solving the linear Maxwell equations. The emitted electromagnetic fields, E^out(x, t) and B^out(x, t), can also be found by solving the linear Maxwell equation with an external source produced by the nonlinear currents. This approach gives the leading order intensity of all harmonics for any polarization of the in-waves (see also references [22, 24, 26, 47]).

From the expressions obtained above it follows that to generate one quantum of the 5th number harmonic in the collision of two counter-propagating laser beams focused in the 1 μm focus spot a laser intensity approximately equal to 5 × 10²⁶ W cm⁻² is needed. Here we assume a 1 μm wavelength laser beam with 30 fs duration. To have the same photon amount for the 3rd harmonic an intensity approximately equal to 3.3 × 10²⁵ W cm⁻² is required. This estimate corresponds to interacting waves with parallel polarizations. The interaction of perpendicular polarized pulses is significantly less efficient in generating high order harmonics.

This work was supported by the project High Field Initiative (CZ.02.1.01/0.0/0.0/15 003/0000449) from the European Regional Development Fund.

No new data were created or analysed in this study.