Dispersion properties, nonlinear waves and birefringence in classical nonlinear electrodynamics

One of the most fundamental and important features of classical electrodynamics is the fact that macroscopically it is a linear theory. This property of Maxwell electrodynamics described in every standard textbook finds remarkable validation in physics experiments and relevant applications almost daily. Put otherwise in a formal language, the linearity of classical electrodynamics sounds like this: any superposition of two or more electromagnetic fields obeying Maxwell's equations satisfies the latter as well [1].

At the subatomic level, however, deviations from the principle of linear superposition can be expected. The standard simplification commonly used in theoretical physics consists of the representation of a charged particle as a localized distribution of charge. Unfortunately, this leads to the infinite growth of its electromagnetic energy with the decrease of the localization dimensions, thus approaching a point-like distribution. To avoid infinite self-energies of point particles, it is natural to speculate that a particular field strength saturation (upper bound) exists. One well-known example of such nonlinear theories developed in the past is the theory of Born and Infeld [2]. Classical electrodynamics cannot describe the interaction between two electromagnetic waves, whereas in quantum electrodynamics, such a scattering of light by light is possible. The two incident plane waves with wave vectors k₁ and k₂ do not merely add coherently, as predicted by linear superposition, but interact and (with a small probability) transform into two different plane waves with corresponding wave vectors k₃ and k₄. These results were first obtained by Euler and Kockel [3] in 1935 and further elaborated by Heisenberg and Euler [4] in 1936. Interestingly enough in both Born-Infeld and Heisenberg-Euler theories of nonlinear electrodynamics in the weak-field limit, the electric polarization ${\bf{P}}$ and the magnetization M vectors are of third-order (up to leading terms) in the electric E and the magnetic B field. This implies that the polarized quantum vacuum can be regarded as a nonlinear medium of Kerr-type.

In the present article, we shall investigate the properties of the nonlinear Heisenberg-Euler model in the weak-field limit from a purely classical electrodynamics point of view. The reason for this is very simple: in the small field approximation, the nonlinear electrodynamics of Heisenberg-Euler and Born-Infeld are surprisingly similar, coinciding up to scaling constants (compare equations (2) and (3) versus equations (5) and (6)) measuring the nonlinearity of the corresponding theories. For the purposes of this work, the origin of these scaling constants, whether quantum or classical, is absolutely irrelevant and is a matter of initial choice and/or preference. Important is the fact that in both cases the media (be it quantum polarized vacuum, or purely classical one of Born-Infeld type) under consideration is of a similar type. This implies that after appropriate and obvious reshaping of the underlying parameters, the results obtained here remain valid for both the Heisenberg-Euler and the Born-Infeld models of nonlinear electrodynamics.

It seems that the manifestation of altered dispersion properties of photons in a constant, or weakly varying electromagnetic background fields have been first studied in detail by Toll [5]. Continuous studies of the linear dispersion properties of photons as well as the vacuum birefringence and dichroism in similar background conditions have been carried out in the early 1960s [6–10], and more recently [11–15]. In order to take into account external field variations for space and time, the weak-field Heisenberg-Euler Lagrangian must be modified by including a term containing derivatives of the electromagnetic field tensor [16, 17]. These corrections result in an additional linear term in both the polarization P and the magnetization M vectors, which is proportional to $\square {\bf{E}}$ and $\square {\bf{B}}$, respectively. The additional linear term would yield a small correction to the linear dispersion properties, which gives us the right to safely ignore it in what follows. The systematic consideration and discussion of dispersion effects in the case of strong electromagnetic fields have been performed in [18].

Nonlinear contributions to the scattering cross section of photons colliding elastically in the presence of an external, static electric field or an electromagnetic wave with the same frequency as the colliding photons have been obtained in the mid 1960s [19, 20]. The analysis of a number of applications, in which nonlinear collective effects among photons may play a dominant role has been initiated by Rozanov [17, 21]. Predominant here refers to processes occurring at sufficiently strong electromagnetic fields in which nonlinear terms in the polarization and magnetization vectors begin to play a significant role. Rozanov considered the perpendicular propagation of high-intensity laser pulses traveling in externally applied electric and magnetic cross fields. By choosing an initial pumping laser pulse polarized in the direction of the background fields, Rozanov obtained a nonlinear Schrödinger equation for the slowly varying envelope of the perturbation [17]. Soljacic and Segev [22] studied the case of a 'beam' resulting from the superposition of two carrier plane waves, modified by a slowly varying envelope in the horizontal direction. They derived a nonlinear Schrödinger equation governing the dynamics of the 'beam' envelope due to the crossing laser waves. Brodin et al [23] considered the propagation of a single TE-mode with weakly modulated amplitude between two parallel conducting planes, and as a result, came up with a similar nonlinear Schrödinger equation for the slowly varying mode envelope.

Compared with the usual properties of various media widely used in nonlinear optics, the vacuum is characterized by both the nonlinear electric polarization and the nonlinear magnetization simultaneously. In the weak-field limit, these nonlinearities give rise to three and four-wave mixing, which has even been argued to be one of the prime candidates to demonstrate the manifestation of nonlinear quantum electrodynamics effects in laboratory conditions. There exist extensive literature on the subject, which is impossible and unnecessary to mention here, so we restrict ourselves to a few selected quotes [21, 24–28].

In the present article, we adopt a different approach, as compared to other works dedicated to nonlinear wave phenomena in polarized vacuum existing in the literature. It is a simple and rather intuitive one, and is widely used in finite and infinite degree-of-freedom dynamical systems. The perturbation method we utilize here is quite straightforward, and as a starting point, it uses the availability of an exact solution to the initial nonlinear wave equations describing the properties of the quantum vacuum. The perturbation expansion about the exact solution of the underlying equations is then built up to third-order, followed by the renormalization of secular terms Thus, the relevant dynamics is being split into two parts: the first one describing the fast wave oscillations, and a second one exhibiting the nonlinear behavior of specific wave amplitudes on much slower spatial and time scales.

The article is organized as follows. In the next section, we present some basics of the weak-field nonlinear electrodynamics of Heisenberg-Euler type. In section 3, we study the free (from external electromagnetic field) quantum vacuum. Small perturbation about the exact plane wave solution of the nonlinear wave equation for the electromagnetic vector potential has been analysed. Employing the Renormalization Group (RG) method [29, 30], we derive a nonlinear wave equation for the (relatively) slowly varying amplitude of the first-order perturbation. An intriguing property of this equation is that it governs the propagation of electromagnetic waves with a reduced speed of light. In section 4, we consider the case where an external magnetic field is applied. The latter introduces spatial anisotropy with all the ensuing consequences. Similar to the case of polarized vacuum without external field, the underlying wave equation for the electromagnetic vector potential possesses an exact solution in the form of a plane wave with constant amplitude. For the analysis carried out in the present article, we choose the circularly polarized exact solution and repeat the RG procedure in a similar manner as has been done in section 3. As a result, we obtain a system of coupled nonlinear wave equations for the two slowly varying amplitudes of the first-order perturbation, which describe the two polarization states. In section 5, we present the stationary (traveling) wave solutions of the nonlinear wave equation, thus obtained for the case of polarized vacuum without external field. The slowly varying wave amplitude behavior is shown to be similar to that of a cnoidal wave, characterized by a near-periodic swell in shallow water. Finally, in section 6 we draw some conclusions.

In Quantum Electrodynamics, photon-photon scattering is known to be a second-order effect in terms of the fine-structure constant α_fs. For constant or weakly varying fields it can be formulated in standard notation by using the Lagrangian (Euler-Heisenberg Lagrangian density [31]) or the Hamiltonian approach [32] in classical field theory. The effects of the quantum electrodynamics vacuum polarization and magnetization introduce (third-order) nonlinearity, such that vacuum itself can be considered a nonlinear medium with appropriate constitutive relations

$$\begin{eqnarray}&&{\bf{D}}={\epsilon }_{0}{\bf{E}}+{\bf{P}},\qquad \qquad {\bf{H}}=\displaystyle \frac{1}{{\mu }_{0}}{\bf{B}}-{\bf{M}}.\end{eqnarray} \tag{ 1 }$$

Here ₀ and μ₀ are the electric permittivity and the magnetic permeability of vacuum $\left({\epsilon }_{0}{\mu }_{0}=1/{c}^{2}\right)$, respectively, D and H are the electric displacement vector and the intensity of the magnetic field, respectively, while P is the electric polarization and M is the magnetization. The latter two quantities are third-order in E and B

$$\begin{eqnarray}&&{\bf{P}}=2\kappa {\epsilon }_{0}^{2}\left[2\left({E}^{2}-{c}^{2}{B}^{2}\right){\bf{E}}+7{c}^{2}\left({\bf{E}}\cdot {\bf{B}}\right){\bf{B}}\right],\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\bf{M}}=2\kappa {\epsilon }_{0}^{2}{c}^{2}\left[-2\left({E}^{2}-{c}^{2}{B}^{2}\right){\bf{B}}+7\left({\bf{E}}\cdot {\bf{B}}\right){\bf{E}}\right],\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&\kappa =\displaystyle \frac{2{\alpha }_{{fs}}^{2}{{\hslash }}^{3}}{45{m}_{e}^{4}{c}^{5}}=\displaystyle \frac{{\alpha }_{{fs}}}{90\pi }\displaystyle \frac{1}{{\epsilon }_{0}{E}_{S}^{2}}\approx \displaystyle \frac{1}{3\times {10}^{29}\ {\rm{J}}\,{{\rm{m}}}^{-3}}.\end{eqnarray} \tag{ 4 }$$

Here m_e is the electron rest mass, and ${E}_{S}={m}_{e}^{2}{c}^{3}/(e{\hslash })\approx 1.32\times {10}^{18}\ {\rm{V}}\,{{\rm{m}}}^{-1}$, is the Schwinger limit above which the electromagnetic field is expected to become nonlinear. Equations (2) and (3) indicate that the nonlinear corrections take the same form as in nonlinear optics, where the material properties of optical fibers for example, give rise to cubic nonlinear terms in Maxwell's equations, so-called Kerr effect. The essential difference however, is that quantum vacuum polarization leads to nonlinearities in both the electric polarization and the magnetization. Note that the effective self-interaction term is proportional to the fine structure constant squared, which implies that field strengths must reach values close to the Schwinger limit until these effects become important.

For comparison, let us adduce the corresponding expressions for the electric polarization and the magnetization in the weak field approximation ensuing from the classical Born-Infeld model [33]. They are

$$\begin{eqnarray}&&{{\bf{P}}}_{{BI}}=\displaystyle \frac{{\epsilon }_{0}}{2{b}^{2}}\left[\left({E}^{2}-{c}^{2}{B}^{2}\right){\bf{E}}+2{c}^{2}\left({\bf{E}}\cdot {\bf{B}}\right){\bf{B}}\right],\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\bf{M}}}_{{BI}}=\displaystyle \frac{{\epsilon }_{0}{c}^{2}}{2{b}^{2}}\left[-\left({E}^{2}-{c}^{2}{B}^{2}\right){\bf{B}}+2\left({\bf{E}}\cdot {\bf{B}}\right){\bf{E}}\right],\end{eqnarray} \tag{ 6 }$$

where the parameter b has the dimension of an electric field strength (Born and Infeld called it the absolute field) estimated initially as b ∼ 3 · 10²⁰ ${\rm{V}}\,{{\rm{m}}}^{-1}$.

The macroscopic Maxwell equations taking into account the vacuum polarization and magnetization can be written as

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot {\bf{D}}={\varrho }_{f},\qquad \qquad {\boldsymbol{\nabla }}\cdot {\bf{B}}=0,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\times {\bf{E}}=-{\partial }_{t}{\bf{B}},\qquad \qquad {\boldsymbol{\nabla }}\times {\bf{H}}-{\partial }_{t}{\bf{D}}={{\bf{J}}}_{f},\end{eqnarray} \tag{ 8 }$$

where ϱ_f and J_f is the free charge and current densities, respectively. Manipulating equations (7) and (8) in an obvious manner with due account of equation (1), we obtain

$$\begin{eqnarray}&&\square {\bf{E}}={\mu }_{0}\left[{\partial }_{t}^{2}{\bf{P}}-{c}^{2}{\boldsymbol{\nabla }}\left({\boldsymbol{\nabla }}\cdot {\bf{P}}\right)+{\partial }_{t}\left({\boldsymbol{\nabla }}\times {\bf{M}}\right)\right]+{\mu }_{0}\left({\partial }_{t}{{\bf{J}}}_{f}+{c}^{2}{\boldsymbol{\nabla }}{\varrho }_{f}\right),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\square {\bf{B}}=-{\mu }_{0}\left[-{{\boldsymbol{\nabla }}}^{2}{\bf{M}}+{\boldsymbol{\nabla }}\left({\boldsymbol{\nabla }}\cdot {\bf{M}}\right)+{\partial }_{t}\left({\boldsymbol{\nabla }}\times {\bf{P}}\right)\right]-{\mu }_{0}{\boldsymbol{\nabla }}\times {{\bf{J}}}_{f},\end{eqnarray} \tag{ 10 }$$

where $\square ={{\boldsymbol{\nabla }}}^{2}-\left(1/{c}^{2}\right){\partial }_{t}^{2}$ is the well-known d'Alembert operator. These are the nonlinear generalizations of the classical wave equations for the electric and the magnetic field, often serving as a starting point in the description of a variety of nonlinear phenomena in quantum vacuum.

In what follows, however, more convenient will be the alternative formulation of nonlinear electrodynamics in the weak-field limit. Since the second of equations (7) and the first of equations (8) remain unaltered as compared to the microscopic electrodynamics, they allow the introduction of the scalar Φ and the vector ${\bf{A}}$ potentials defined in the standard manner

$$\begin{eqnarray}&&{\bf{E}}=-{\boldsymbol{\nabla }}{\rm{\Phi }}-{\partial }_{t}{\bf{A}},\qquad \qquad {\bf{B}}={\boldsymbol{\nabla }}\times {\bf{A}}.\end{eqnarray} \tag{ 11 }$$

Manipulating the remaining two of the Maxwell equations and taking into account the definition (1) for the electric displacement vector and the intensity of the magnetic field, we obtain

$$\begin{eqnarray}&&\square {\bf{A}}-{\boldsymbol{\nabla }}\left(\displaystyle \frac{1}{{c}^{2}}{\partial }_{t}{\rm{\Phi }}+{\boldsymbol{\nabla }}\cdot {\bf{A}}\right)=-{\mu }_{0}\left({\partial }_{t}{\bf{P}}+{\boldsymbol{\nabla }}\times {\bf{M}}\right)-{\mu }_{0}{{\bf{J}}}_{f},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}^{2}{\rm{\Phi }}+{\partial }_{t}\left({\boldsymbol{\nabla }}\cdot {\bf{A}}\right)=\displaystyle \frac{1}{{\epsilon }_{0}}{\boldsymbol{\nabla }}\cdot {\bf{P}}-\displaystyle \frac{{\varrho }_{f}}{{\epsilon }_{0}}.\end{eqnarray} \tag{ 13 }$$

In addition to the nonlinear wave equations above, one can use the Lorentz $\left[\left(1/{c}^{2}\right){\partial }_{t}{\rm{\Phi }}+{\boldsymbol{\nabla }}\cdot {\bf{A}}=0\right]$ or the Coulomb $\left({\boldsymbol{\nabla }}\cdot {\bf{A}}=0\right)$ gauge condition.

Equations (12) and (13) together with equations (2) and (3) describe the nonlinear wave properties of polarized vacuum in the weak-field limit. The case, where free electric charges ϱ_f and currents J_f are absent will be the starting point for our subsequent analysis. In other words, throughout the subsequent exposition only the case of zero free electric charges ϱ_f = 0 and currents J_f = 0 will be considered.

We assume a field configuration as follows

$$\begin{eqnarray}&&{\bf{E}}=\left(0,\ -{\partial }_{t}A,\ 0\right),\qquad {\bf{B}}=\left(-{\partial }_{z}A,\ 0,\ {\partial }_{x}A\right),\end{eqnarray} \tag{ 14 }$$

which is determined by a vector potential

$$\begin{eqnarray}&&{\bf{A}}\left(x,z;t\right)=\left[0,\ A\left(x,z;t\right),\ 0\right].\end{eqnarray} \tag{ 15 }$$

Since E · B = 0 in this configuration, the second terms in equations (2) and (3) depending on the scalar product of the electric field E and the magnetic field B vanish. As far as components are concerned, obviously, the vacuum polarization P and the vacuum magnetization M follow the pattern of E and B

$$\begin{eqnarray}&&{P}_{y}=-4\kappa {\epsilon }_{0}^{2}{c}^{2}{\mathscr{F}}{\partial }_{t}A,\qquad {\bf{M}}=-4\kappa {\epsilon }_{0}^{2}{c}^{4}{\mathscr{F}}{\boldsymbol{\nabla }}\times {\bf{A}},\end{eqnarray} \tag{ 16 }$$

respectively. Here

$$\begin{eqnarray}&&{\mathscr{F}}=\displaystyle \frac{1}{{c}^{2}}{\left({\partial }_{t}A\right)}^{2}-{\left({\boldsymbol{\nabla }}A\right)}^{2}.\end{eqnarray} \tag{ 17 }$$

It is convenient to use the Coulomb gauge, in which the scalar potential Φ vanishes identically.

It can be verified in a straightforward manner that the nonlinear wave equation (12) for the vector potential can be rewritten as

$$\begin{eqnarray}&&\square A=4\kappa {\epsilon }_{0}{\partial }_{t}\left({\mathscr{F}}{\partial }_{t}A\right)-4\kappa {\epsilon }_{0}{c}^{2}\left({\boldsymbol{\nabla }}A\cdot {\boldsymbol{\nabla }}{\mathscr{F}}+{\mathscr{F}}{{\boldsymbol{\nabla }}}^{2}A\right).\end{eqnarray} \tag{ 18 }$$

Note that a plane wave of the form

$$\begin{eqnarray}&&{A}_{0}={\mathscr{C}}{{\rm{e}}}^{i{\varphi }_{0}}+{{\mathscr{C}}}^{* }{{\rm{e}}}^{-i{\varphi }_{0}},\end{eqnarray} \tag{ 19 }$$

where ${\mathscr{C}}$ is a constant complex amplitude, ${\varphi }_{0}={{\bf{k}}}_{0}\cdot {\bf{x}}-{\omega }_{0}t$ is the wave phase, and the wave frequency ω₀ and the wave number k₀ satisfy the dispersion relation

$$\begin{eqnarray}&&{\omega }_{0}^{2}={c}^{2}\left({k}_{0x}^{2}+{k}_{0z}^{2}\right)={c}^{2}{k}_{0}^{2},\qquad \qquad {k}_{0y}=0,\end{eqnarray} \tag{ 20 }$$

is an exact solution of equation (18). This property is not surprising—it immediately follows from the important relation ${{\mathscr{F}}}_{0}={\mathscr{F}}\left({A}_{0}\right)=0$, which can be checked by direct substitution. Put another way, the plane wave minimizes the effective Heisenberg-Euler action.

Following the standard procedure of the renormalization group method, we represent $A\left(x,z;t\right)$ as a perturbation expansion

$$\begin{eqnarray}&&A\left(x,z;t\right)=\sum _{n=0}^{\infty }{\epsilon }^{n}{A}_{n}\left(x,z;t\right),\end{eqnarray} \tag{ 21 }$$

in the formal small parameter . The next step consists in expanding equation (18) in the small parameter , and obtaining its naive perturbation solution order by order.

3.1. First order

The first-order vector potential A₁ obeys the equation

$$\begin{eqnarray}&&{\widehat{{\bf{L}}}}_{1}{A}_{1}=-\alpha \left[{{\rm{e}}}^{2i\left({\varphi }_{0}+\delta \right)}+{{\rm{e}}}^{-2i\left({\varphi }_{0}+\delta \right)}\right]{\widehat{{\bf{D}}}}_{1}^{2}{A}_{1},\end{eqnarray} \tag{ 22 }$$

where

$$\begin{eqnarray}&&{\widehat{{\bf{D}}}}_{1}=\displaystyle \frac{{\omega }_{0}}{c}{\partial }_{t}+c{{\bf{k}}}_{0}\cdot {\boldsymbol{\nabla }},\qquad \quad {\widehat{{\bf{L}}}}_{1}=\square -2\alpha {\widehat{{\bf{D}}}}_{1}^{2},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\alpha =8\kappa {\epsilon }_{0}{\left|{\mathscr{C}}\right|}^{2},\end{eqnarray} \tag{ 24 }$$

and δ is the phase of the constant complex amplitude ${\mathscr{C}}$. The right-hand-side of equation (22) is fast oscillating, and can be averaged away. Thus, we obtain the first-order solution

$$\begin{eqnarray}&&{A}_{1}={\mathscr{A}}{{\rm{e}}}^{i{\varphi }_{1}}+{{\mathscr{A}}}^{* }{{\rm{e}}}^{-i{\varphi }_{1}}.\end{eqnarray} \tag{ 25 }$$

Here, ${\mathscr{A}}$ is an additional constant (for the present moment in time) complex amplitude and ${\varphi }_{1}={\bf{k}}\cdot {\bf{x}}-\omega t$ is the wave phase with k_y=0. In addition, the wave frequency

$$\begin{eqnarray}&&\omega \left({\bf{k}}\right)=\displaystyle \frac{1}{1+2\alpha {\omega }_{0}^{2}}\left[2\alpha {c}^{2}{\omega }_{0}{{\bf{k}}}_{0}\cdot {\bf{k}}+\sqrt{{c}^{2}{k}^{2}\left(1+2\alpha {\omega }_{0}^{2}\right)-2\alpha {c}^{4}{\left({{\bf{k}}}_{0}\cdot {\bf{k}}\right)}^{2}}\right].\end{eqnarray} \tag{ 26 }$$

is a solution to the dispersion relation

$$\begin{eqnarray}&&{\mathscr{D}}\left({\bf{k}},\omega \right)=\displaystyle \frac{{\omega }^{2}}{{c}^{2}}-{k}^{2}+2\alpha {\left(\displaystyle \frac{{\omega }_{0}\omega }{c}-c{{\bf{k}}}_{0}\cdot {\bf{k}}\right)}^{2}=0.\end{eqnarray} \tag{ 27 }$$

An important comment is now in order. Apart from the solution (25), the first-order equation (22) possesses a solution of the form (19), which has been disregarded, since it is already contained in the exact plane wave solution of the initial nonlinear wave equation. Furthermore, equation (22) can be regarded as a wave analog of the Mathieu equation. Perturbation analysis shows that the fast oscillating term proportional to $\cos 2\left(\varphi +\delta \right)$ does not give rise to additional higher-order contribution in terms of the polarization parameter κ. In other words, the neglection (or equivalently, the averaging procedure) of the fast oscillating term is justified up to (at least) second-order in perturbation theory as applied to the Mathieu equation (22). Further technical details can be found in appendix.

3.2. Second order

The second-order perturbation equation for the vector potential A₂ reads as

$$\begin{eqnarray}&&\square {A}_{2}=4\kappa {\epsilon }_{0}{\partial }_{t}\left({{\mathscr{F}}}_{2}{\partial }_{t}{A}_{0}+{{\mathscr{F}}}_{1}{\partial }_{t}{A}_{1}\right)-4\kappa {\epsilon }_{0}{c}^{2}\left({\boldsymbol{\nabla }}{A}_{0}\cdot {\boldsymbol{\nabla }}{{\mathscr{F}}}_{2}+{\boldsymbol{\nabla }}{A}_{1}\cdot {\boldsymbol{\nabla }}{{\mathscr{F}}}_{1}+{{\mathscr{F}}}_{2}{{\boldsymbol{\nabla }}}^{2}{A}_{0}+{{\mathscr{F}}}_{1}{{\boldsymbol{\nabla }}}^{2}{A}_{1}\right).\end{eqnarray} \tag{ 28 }$$

Substituting equations (19) and (25) into the right-hand-side of equation (28), we rewrite the latter in a detailed form as follows

$$\begin{eqnarray}&&\begin{array}{l}{\widehat{{\bf{L}}}}_{1}{A}_{2}+2\alpha \cos 2\left({\varphi }_{0}+\delta \right){\widehat{{\bf{D}}}}_{1}^{2}{A}_{2}\\ \quad =\,8\kappa {\epsilon }_{0}{c}^{2}{\square }_{0k}\left[\left({\square }_{0k}+3{\square }_{k}\right){{\mathscr{CA}}}^{2}{{\rm{e}}}^{i\left(2{\varphi }_{1}+{\varphi }_{0}\right)}+\left({\square }_{0k}-3{\square }_{k}\right){{\mathscr{C}}}^{* }{{\mathscr{A}}}^{2}{{\rm{e}}}^{i\left(2{\varphi }_{1}-{\varphi }_{0}\right)}+c.c.\right],\end{array}\end{eqnarray} \tag{ 29 }$$

where

$$\begin{eqnarray}&&{\square }_{k}=\displaystyle \frac{{\omega }^{2}}{{c}^{2}}-{k}^{2},\qquad \qquad {\square }_{0k}=\displaystyle \frac{{\omega }_{0}\omega }{{c}^{2}}-{{\bf{k}}}_{0}\cdot {\bf{k}}.\end{eqnarray} \tag{ 30 }$$

The solution of the second-order perturbation equation (29) can be written as

$$\begin{eqnarray}&&{A}_{2}=2\kappa {\epsilon }_{0}{c}^{2}\left[\left({\square }_{0k}+3{\square }_{k}\right){{\mathscr{CA}}}^{2}{{\rm{e}}}^{i\left(2{\varphi }_{1}+{\varphi }_{0}\right)}-\left({\square }_{0k}-3{\square }_{k}\right){{\mathscr{C}}}^{* }{{\mathscr{A}}}^{2}{{\rm{e}}}^{i\left(2{\varphi }_{1}-{\varphi }_{0}\right)}+c.c.\right].\end{eqnarray} \tag{ 31 }$$

Note that the second-order vector potential A₂ includes terms with combination phases of the form ${{\rm{e}}}^{i\left(2{\varphi }_{1}\pm {\varphi }_{0}\right)}$. They describe the three wave mixing and are giving rise to secular terms as will be shown below.

3.3. Third order. Derivation of the nonlinear amplitude equation

The third-order perturbation equation for vector potential A₃ reads as

$$\begin{eqnarray}\begin{array}{rcl}\square {A}_{3} & = & 4\kappa {\epsilon }_{0}{\partial }_{t}\left({{\mathscr{F}}}_{3}{\partial }_{t}{A}_{0}+{{\mathscr{F}}}_{2}{\partial }_{t}{A}_{1}+{{\mathscr{F}}}_{1}{\partial }_{t}{A}_{2}\right)\\ & & -4\kappa {\epsilon }_{0}{c}^{2}\left({\boldsymbol{\nabla }}{A}_{0}\cdot {\boldsymbol{\nabla }}{{\mathscr{F}}}_{3}+{\boldsymbol{\nabla }}{A}_{1}\cdot {\boldsymbol{\nabla }}{{\mathscr{F}}}_{2}+{\boldsymbol{\nabla }}{A}_{2}\cdot {\boldsymbol{\nabla }}{{\mathscr{F}}}_{1}+{{\mathscr{F}}}_{3}{{\boldsymbol{\nabla }}}^{2}{A}_{0}+{{\mathscr{F}}}_{2}{{\boldsymbol{\nabla }}}^{2}{A}_{1}+{{\mathscr{F}}}_{1}{{\boldsymbol{\nabla }}}^{2}{A}_{2}\right).\end{array}\end{eqnarray} \tag{ 32 }$$

Unlike the second-order, where the perturbation equation possesses a unique regular solution, there are two types of terms on the right-hand-side of the third-order equation (32). The first type comprises a collection of secular (resonant) terms, which follow the pattern of the basic wave mode in linear approximation (proportional to ${{\rm{e}}}^{\pm i{\varphi }_{1}}$). Such terms, which appear in higher orders as well, would provide a divergent counterpart in the naive perturbation solution. Thus, they must be renormalized by an elegant procedure described below. The rest of the terms contribute to the regular solution of the third-order perturbation equations, involving higher harmonics and/or higher order harmonic combinations of the modes proportional to ${{\rm{e}}}^{i{\varphi }_{0}}$ and ${{\rm{e}}}^{i{\varphi }_{1}}$. Omitting straightforwardly reproducible calculation's details, we write down the resonant part of equation (32)

$$\begin{eqnarray}&&{\widehat{{\bf{L}}}}_{1}{A}_{3}+2\alpha \cos 2\left({\varphi }_{0}+\delta \right){\widehat{{\bf{D}}}}_{1}^{2}{A}_{3}=-36\kappa {\epsilon }_{0}{c}^{2}{\square }_{k}^{2}{\left|{\mathscr{A}}\right|}^{2}{\mathscr{A}}{{\rm{e}}}^{i{\varphi }_{1}}+c.c.,\end{eqnarray} \tag{ 33 }$$

The above equation possesses an exact solution

$$\begin{eqnarray}&&{A}_{3}\left(x,z;t\right)={{\mathscr{P}}}_{3}\left(x,z;t\right){{\rm{e}}}^{i{\varphi }_{1}},\end{eqnarray} \tag{ 34 }$$

where the amplitude ${{\mathscr{P}}}_{3}\left(x,z;t\right)$ satisfies the equation

$$\begin{eqnarray}&&\left[{\widehat{{\bf{L}}}}_{1}+i\left({\partial }_{\omega }{\mathscr{D}}\right){\partial }_{t}-i\left({{\boldsymbol{\nabla }}}_{k}{\mathscr{D}}\right)\cdot {\boldsymbol{\nabla }}\right]{{\mathscr{P}}}_{3}=-36\kappa {\epsilon }_{0}{c}^{2}{\square }_{k}^{2}{\left|{\mathscr{A}}\right|}^{2}{\mathscr{A}}.\end{eqnarray} \tag{ 35 }$$

Here ${\mathscr{D}}\left({\bf{k}},\omega \right)$ is the dispersion function defined by equation (27).

We follow an elegant approach, known as the proto RG operator scheme [34–36], which has been proposed in the early 2000s to free as much as possible the standard RG theoretical reduction from the necessity of explicit (in the majority of cases, rather cumbersome) calculation of secular terms [see e.g. 29 and 30]. Here we present only a brief description of the renormalization procedure; the interested reader can consult the references cited above, where further details and other interesting applications of the proto RG operator approach can be found. Since, according to the first-order perturbative solution (25) the complex amplitude ${\mathscr{A}}$ (regarded as a free parameter) is constant, the solution of equation (35) will contain secular terms, which must be resummed appropriately. The usual procedure to do so is first to define the free parameters x_R, z_R and τ and the renormalized amplitude ${{\mathscr{A}}}_{R}$ according to the relation ${\mathscr{A}}={{\mathscr{ZA}}}_{R}$, where ${\mathscr{Z}}$ is the renormalization constant. Omitting zero-order and second-order solutions, which are regular and thus independent of the free parameters just introduced, we can write

$$\begin{eqnarray}&&A\left(x,z;t\right)=\epsilon {{\mathscr{A}}}_{R}\left({x}_{R},{z}_{R};\tau \right){{\rm{e}}}^{i{\varphi }_{1}}+{\epsilon }^{3}\left[{{\mathscr{P}}}_{3}\left(x,z;t\right)-{{\mathscr{P}}}_{3}\left({x}_{R},{z}_{R};\tau \right)\right]{{\rm{e}}}^{i{\varphi }_{1}}+...+c.c..\end{eqnarray} \tag{ 36 }$$

Let $\widehat{{\boldsymbol{\Sigma }}}\left({x}_{R},{z}_{R};\tau \right)$ denotes the operator in the square brackets on the left-hand-side of equation (35) with $x,z$ and t being replaced by x_R, z_R and τ. Provided that $A\left(x,z;t\right)$ should not depend on the renormalization parameters, we obtain

$$\begin{eqnarray}&&0=\widehat{{\boldsymbol{\Sigma }}}\left({x}_{R},{z}_{R};\tau \right)A=\epsilon \widehat{{\boldsymbol{\Sigma }}}\left({x}_{R},{z}_{R};\tau \right)\left[{{\mathscr{A}}}_{R}-{\epsilon }^{2}{{\mathscr{P}}}_{3}\left({x}_{R},{z}_{R};\tau | {{\mathscr{A}}}_{R}\right)\right]{{\rm{e}}}^{i{\varphi }_{1}}.\end{eqnarray} \tag{ 37 }$$

Replacing back the original independent variables and making use of equation (35), we finally arrive at the sought for nonlinear amplitude equation

$$\begin{eqnarray}&&\left[{\widehat{{\bf{L}}}}_{1}+i\left({\partial }_{\omega }{\mathscr{D}}\right){\partial }_{t}-i\left({{\boldsymbol{\nabla }}}_{k}{\mathscr{D}}\right)\cdot {\boldsymbol{\nabla }}\right]{\mathscr{A}}=-36\kappa {\epsilon }_{0}{c}^{2}{\square }_{k}^{2}{\left|{\mathscr{A}}\right|}^{2}{\mathscr{A}}.\end{eqnarray} \tag{ 38 }$$

It governs the relatively slow dynamics of the wave envelope ${\mathscr{A}}$ and describes the formation of nonlinear waves and coherent structures. As it will be shown in section 5 the equation just derived is a nonlinear wave equation with a characteristic reduced light velocity.

In this section, we shall analyze the case, where a constant external magnetic field ${{\bf{B}}}_{0}=\left({B}_{0},0,0\right)$ is applied along the x-axis. Since the applied magnetic field introduces a special type of spatial anisotropy, we assume that the induced electromagnetic waves propagate in the axial x-direction only. The electromagnetic field configuration can be specified as follows

$$\begin{eqnarray}&&{\bf{E}}=\left(0,\ -{\partial }_{t}{A}_{y},\ -{\partial }_{t}{A}_{z}\right),\qquad {\bf{B}}=\left({B}_{0},\ -{\partial }_{x}{A}_{z},\ {\partial }_{x}{A}_{y}\right),\end{eqnarray} \tag{ 39 }$$

which is defined by an electromagnetic vector potential of the form

$$\begin{eqnarray}&&{\bf{A}}\left(x;t\right)=\left[0,\ {A}_{y}\left(x;t\right),\ {A}_{z}\left(x;t\right)\right].\end{eqnarray} \tag{ 40 }$$

Since for plane waves this configuration yields ${\bf{E}}\cdot {\bf{B}}=0$, the second terms in equations (2) and (3) depending on the scalar product of the electric field E and the magnetic field B vanish correspondingly.

It will prove convenient for the subsequent exposition to introduce a new field variable according to the relation

$$\begin{eqnarray}&&{\mathscr{A}}\left(x;t\right)={A}_{y}\left(x;t\right)+{{iA}}_{z}\left(x;t\right).\end{eqnarray} \tag{ 41 }$$

Obviously,

$$\begin{eqnarray}&&{E}^{2}={\partial }_{t}{\mathscr{A}}\left({\partial }_{t}{{\mathscr{A}}}^{* }\right),\qquad {B}^{2}={B}_{0}^{2}+{\partial }_{x}{\mathscr{A}}\left({\partial }_{x}{{\mathscr{A}}}^{* }\right).\end{eqnarray} \tag{ 42 }$$

As far as components are concerned, similarly to the preceding section the vacuum polarization P and the vacuum magnetization M follow the vectorial pattern of E and B

$$\begin{eqnarray}&&{P}_{x}=0,\quad \qquad {P}_{y}=-4\kappa {\epsilon }_{0}^{2}{c}^{2}{\mathscr{G}}{\partial }_{t}{A}_{y},\quad \qquad {P}_{z}=-4\kappa {\epsilon }_{0}^{2}{c}^{2}{\mathscr{G}}{\partial }_{t}{A}_{z},\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{M}_{x}=-4\kappa {\epsilon }_{0}^{2}{c}^{4}{\mathscr{G}}{B}_{0},\quad \qquad {M}_{y}=4\kappa {\epsilon }_{0}^{2}{c}^{4}{\mathscr{G}}{\partial }_{x}{A}_{z},\quad \qquad {M}_{z}=-4\kappa {\epsilon }_{0}^{2}{c}^{4}{\mathscr{G}}{\partial }_{x}{A}_{y},\end{eqnarray} \tag{ 44 }$$

respectively. Here

$$\begin{eqnarray}&&{\mathscr{G}}=\displaystyle \frac{1}{{c}^{2}}{\partial }_{t}{\mathscr{A}}\left({\partial }_{t}{{\mathscr{A}}}^{* }\right)-{\partial }_{x}{\mathscr{A}}\left({\partial }_{x}{{\mathscr{A}}}^{* }\right)-{B}_{0}^{2}.\end{eqnarray} \tag{ 45 }$$

As it was done in the case of a polarized quantum vacuum without externally applied fields, we again use the Coulomb gauge, in which the scalar potential Φ vanishes identically.

It can be verified in a straightforward manner that the nonlinear wave equation (12) for the vector complex potential can be rewritten as

$$\begin{eqnarray}&&\square {\mathscr{A}}=4\kappa {\epsilon }_{0}{\partial }_{t}\left({\mathscr{G}}{\partial }_{t}{\mathscr{A}}\right)-4\kappa {\epsilon }_{0}{c}^{2}{\partial }_{x}\left({\mathscr{G}}{\partial }_{x}{\mathscr{A}}\right).\end{eqnarray} \tag{ 46 }$$

Again, it can be checked by direct substitution that a plane wave of the form

$$\begin{eqnarray}&&{{\mathscr{A}}}_{0}={\mathscr{E}}\left({{\rm{e}}}^{i{\varphi }_{0}}+{\mathscr{P}}{{\rm{e}}}^{-i{\varphi }_{0}}\right),\end{eqnarray} \tag{ 47 }$$

where ${\mathscr{E}}$ is a constant complex amplitude, ${\varphi }_{0}={k}_{0}x-{\omega }_{0}t$ is the wave phase, and the wave frequency ω₀ and the wave number k₀ satisfy the dispersion relation ${\omega }_{0}^{2}={c}^{2}{k}_{0}^{2}$, is an exact solution of equation (46). In addition, ${\mathscr{P}}$ is the polarization parameter, taking into account various plane wave polarizations. The assertion that the plane wave (47) is an exact solution immediately follows from the important property ${{\mathscr{G}}}_{0}={\mathscr{G}}\left({{\mathscr{A}}}_{0}\right)=-{B}_{0}^{2}$. For the sake of simplicity we consider here the case of circular polarization, that is ${\mathscr{P}}=0$. Since the direction of propagation is collinear with E × B, the above simplification implies that the driving wave propagates along the x-axis and is polarized in the yz-plane.

Similar to the preceding section, we represent ${\mathscr{A}}$ as a perturbation expansion

$$\begin{eqnarray}&&{\mathscr{A}}\left(x;t\right)=\sum _{n=0}^{\infty }{\epsilon }^{n}{{\mathscr{A}}}_{n}\left(x;t\right),\end{eqnarray} \tag{ 48 }$$

in the formal small parameter . The first-order vector potential ${{\mathscr{A}}}_{1}$ obeys now the equation

$$\begin{eqnarray}&&{\widehat{{\bf{L}}}}_{M}{{\mathscr{A}}}_{1}=-\beta {{\rm{e}}}^{2i\left({\varphi }_{0}+\delta \right)}{\widehat{{\bf{D}}}}_{M}^{2}{{\mathscr{A}}}_{1}^{* },\end{eqnarray} \tag{ 49 }$$

where

$$\begin{eqnarray}&&{\widehat{{\bf{D}}}}_{M}=\displaystyle \frac{{\omega }_{0}}{c}{\partial }_{t}+{{ck}}_{0}{\partial }_{x},\qquad \quad {\widehat{{\bf{L}}}}_{M}={\rm{\Gamma }}\square -\beta {\widehat{{\bf{D}}}}_{M}^{2},\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&\beta =4\kappa {\epsilon }_{0}{\left|{\mathscr{E}}\right|}^{2},\qquad \qquad {\rm{\Gamma }}=1-4\kappa {\epsilon }_{0}{c}^{2}{B}_{0}^{2},\end{eqnarray} \tag{ 51 }$$

and δ is the phase of the constant complex amplitude ${\mathscr{E}}$. Applying again the averaging procedure described in the preceding section, we neglect the fast oscillating term on the right-hand-side of equation (49), and obtain the first-order solution

$$\begin{eqnarray}&&{{\mathscr{A}}}_{1}={{\mathscr{B}}}_{+}{{\rm{e}}}^{i{\varphi }_{1}}+{{\mathscr{B}}}_{-}{{\rm{e}}}^{-i{\varphi }_{1}}.\end{eqnarray} \tag{ 52 }$$

Here, ${{\mathscr{B}}}_{\pm }$ are supplementary constant complex amplitudes and ${\varphi }_{1}={kx}-\omega t$ is the wave phase. In addition, $\omega \left(k\right)$ is the wave frequency

$$\begin{eqnarray}&&\omega \left(k\right)=\displaystyle \frac{1}{{\rm{\Gamma }}+\beta {\omega }_{0}^{2}}\left[\beta {c}^{2}{\omega }_{0}{k}_{0}k+\sqrt{{\rm{\Gamma }}{c}^{2}{k}^{2}\left({\rm{\Gamma }}+\beta {\omega }_{0}^{2}\right)-{\rm{\Gamma }}\beta {c}^{4}{k}_{0}^{2}{k}^{2}}\right],\end{eqnarray} \tag{ 53 }$$

which is a solution to the dispersion relation

$$\begin{eqnarray}&&{{\mathscr{D}}}_{M}\left(k,\omega \right)={\rm{\Gamma }}{\square }_{k}+\beta {c}^{2}{\square }_{0k}^{2}=0.\end{eqnarray} \tag{ 54 }$$

Since the main considerations and explicit calculations in higher orders are confined in the same mainstream of particulars as described in detail in section 3, we will skip them and will present only the final result. The sought for nonlinear amplitude equations for the slowly varying wave amplitudes ${{\mathscr{B}}}_{\pm }$ are

$$\begin{eqnarray}&&\begin{array}{l}\left[{\widehat{{\bf{L}}}}_{M}+i\left({\partial }_{\omega }{{\mathscr{D}}}_{M}\right){\partial }_{t}-i\left({\partial }_{k}{{\mathscr{D}}}_{M}\right){\partial }_{x}\right]{{\mathscr{B}}}_{+}\\ \quad =\,-24\kappa {\epsilon }_{0}{c}^{2}{\square }_{k}^{2}{\left|{{\mathscr{B}}}_{-}\right|}^{2}{{\mathscr{B}}}_{+}+2\kappa {\epsilon }_{0}{c}^{2}\displaystyle \frac{{\square }_{k}}{{\square }_{0k}}\left(4{\square }_{k}^{2}-6{\square }_{k}{\square }_{0k}+{\square }_{0k}^{2}\right){\left|{{\mathscr{B}}}_{+}\right|}^{2}{{\mathscr{B}}}_{+},\end{array}\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&\begin{array}{l}\left[{\widehat{{\bf{L}}}}_{M}-i\left({\partial }_{\omega }{{\mathscr{D}}}_{M}\right){\partial }_{t}+i\left({\partial }_{k}{{\mathscr{D}}}_{M}\right){\partial }_{x}\right]{{\mathscr{B}}}_{-}\\ \quad =\,-24\kappa {\epsilon }_{0}{c}^{2}{\square }_{k}^{2}{\left|{{\mathscr{B}}}_{+}\right|}^{2}{{\mathscr{B}}}_{-}-2\kappa {\epsilon }_{0}{c}^{2}\displaystyle \frac{{\square }_{k}}{{\square }_{0k}}\left(4{\square }_{k}^{2}+6{\square }_{k}{\square }_{0k}+{\square }_{0k}^{2}\right){\left|{{\mathscr{B}}}_{-}\right|}^{2}{{\mathscr{B}}}_{-}.\end{array}\end{eqnarray} \tag{ 56 }$$

The above equations describe the nonlinear optical anisotropy of the polarized quantum vacuum and the presence of two intrinsic polarization states with different indices of refraction (nonlinear birefringence of the vacuum in an externally applied magnetic field). In other words, coherent interaction of photons in the pumping plane wave causes each of the induced waves to acquire an effective nonlinear refractive index. Obviously, the value of this index depends on the intensity of the pumping wave, as well as on the amplitude of the waves characterized by the two polarization states. After traversing a certain distance, a circularly polarized impacting wave is transformed into an elliptically polarized wave of the most general form. Polarization phenomena accompanying the interaction of two laser beams, which are due to scattering of light by light have been considered in the past [37].

It is instructive to transform equations (38), (55) and (56) in a more convenient and familiar form. Let us deal with equation (38) first. Taking into account the translational invariance of our problem, we introduce new spatial variables according to the relation

$$\begin{eqnarray}&&{\boldsymbol{\xi }}={\bf{x}}-{\bf{v}}t,\qquad {\bf{v}}=\displaystyle \frac{2\alpha {\omega }_{0}^{2}}{1+2\alpha {\omega }_{0}^{2}}c{\widehat{{\bf{k}}}}_{0},\end{eqnarray} \tag{ 57 }$$

where ${\widehat{{\bf{k}}}}_{0}$ is the unit wave vector. This means that the new coordinate system is firmly fixed with an observer moving with a speed v (which for realistic impact wave intensities is a small number), as defined by the second of the above relations. In the new spatial variables defined above, the nonlinear amplitude equation (38) can be rewritten as

$$\begin{eqnarray}&&i\left({\partial }_{\omega }{\mathscr{D}}\right)\left[{\partial }_{t}+\left({{\bf{v}}}_{g}-{\bf{v}}\right)\cdot {{\boldsymbol{\nabla }}}_{{\boldsymbol{\xi }}}\right]{\mathscr{A}}+\left[{{\boldsymbol{\nabla }}}_{{\boldsymbol{\xi }}}^{2}-\displaystyle \frac{2\alpha {\omega }_{0}^{2}}{1+2\alpha {\omega }_{0}^{2}}{\left({\widehat{{\bf{k}}}}_{0}\cdot {{\boldsymbol{\nabla }}}_{{\boldsymbol{\xi }}}\right)}^{2}-\displaystyle \frac{1}{{c}^{2}}\left(1+2\alpha {\omega }_{0}^{2}\right){\partial }_{t}^{2}\right]{\mathscr{A}}=-\sigma {\left|{\mathscr{A}}\right|}^{2}{\mathscr{A}}.\end{eqnarray} \tag{ 58 }$$

Here ${{\bf{v}}}_{g}=-\left({{\boldsymbol{\nabla }}}_{k}{\mathscr{D}}\right){\left({\partial }_{\omega }{\mathscr{D}}\right)}^{-1}$ is the wave group velocity [38] and

$$\begin{eqnarray}&&\sigma =36\kappa {\epsilon }_{0}{c}^{2}{\square }_{k}^{2},\end{eqnarray} \tag{ 59 }$$

is the nonlinear coupling coefficient.

Let us now assume that the vector potential (15) depends on one spatial coordinate (say the x-coordinate) only. Thus, equation (58) acquires the form

$$\begin{eqnarray}&&i{{\mathscr{D}}}_{\omega }\left({\partial }_{t}+{\rm{\Delta }}v{\partial }_{\xi }\right){\mathscr{A}}+\displaystyle \frac{1}{1+2\alpha {\omega }_{0}^{2}}{\widehat{\square }}_{p}{\mathscr{A}}=-\sigma {\left|{\mathscr{A}}\right|}^{2}{\mathscr{A}},\end{eqnarray} \tag{ 60 }$$

where

$$\begin{eqnarray}&&{\widehat{\square }}_{p}={{\boldsymbol{\nabla }}}_{{\boldsymbol{\xi }}}^{2}-\displaystyle \frac{1}{{c}_{1}^{2}}{\partial }_{t}^{2},\qquad \quad {c}_{1}=\displaystyle \frac{c}{1+2\alpha {\omega }_{0}^{2}}.\end{eqnarray} \tag{ 61 }$$

In addition, ${{\mathscr{D}}}_{\omega }$ is a short-hand notation of the derivative of the dispersion function (27) with respect to ω and Δ v stands for the difference between the group velocity and the speed v defined by equation (57) $\left({\rm{\Delta }}v={v}_{g}-v\right)$.

Here is the place to make a small diversion from the main exposition with an important comment. The nonlinear wave equation (60) governs the evolution of electromagnetic waves with a reduced speed of light c₁. For sufficiently small values of the parameter α, the relative light velocity reduction ${\rm{\Delta }}c/c=2\alpha {\omega }_{0}^{2}{\left(1+2\alpha {\omega }_{0}^{2}\right)}^{-1}$ is roughly proportional to the intensity of the pumping plane wave (cf equation (24)). For example, at the maximum achievable value of the laser parameter a₀ = 170 at ELI-NP, this decrease is of the order of ${\rm{\Delta }}c/c\sim 1.045\cdot {10}^{-10}$. In the interstellar space, where the electromagnetic wave intensity (corresponding to an equivalent laser parameter) can be of several orders of magnitude greater, the corresponding reduction in the speed of light is of the order of ${\rm{\Delta }}c/c\sim 0.01-1.0 \%$. This implies that the effect mentioned above is expected to become more profound at a very high intensity of the impacting electromagnetic fields, and could be confirmed experimentally by monitoring and collecting of relevant data in a supernova explosion. A similar effect, although in a different context, was recently reported [39, 40].

In order to find a stationary wave solution, we will take advantage of the universal ansatz

$$\begin{eqnarray}&&{\mathscr{A}}\left(\xi ;t\right)=B\left(\zeta \right){{\rm{e}}}^{i\psi \left(\xi ;t\right)},\qquad \quad \zeta =\lambda \xi -{ut},\end{eqnarray} \tag{ 62 }$$

where λ is a constant, and u is a constant phase velocity of the traveling wave, both to be specified in the sequel. Substituting the above expression for ${\mathscr{A}}$ into equation (60) and separating real and imaginary parts, we obtain two coupled differential equations relating the amplitude B and the phase ψ. These can be written as

$$\begin{eqnarray}&&B{\widehat{\square }}_{p}\psi +2\left(\lambda {\psi }_{\xi }+\displaystyle \frac{u\dot{\psi }}{{c}_{1}^{2}}\right){B}^{{\prime} }-{{\mathscr{D}}}_{\omega }\left(1+2\alpha {\omega }_{0}^{2}\right)\left(u-\lambda {\rm{\Delta }}v\right){B}^{{\prime} }=0,\end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray}&&\left({\lambda }^{2}-\displaystyle \frac{{u}^{2}}{{c}_{1}^{2}}\right){B}^{{\prime\prime} }-\left[{\psi }_{\xi }^{2}-\displaystyle \frac{{\dot{\psi }}^{2}}{{c}_{1}^{2}}+{{\mathscr{D}}}_{\omega }\left(1+2\alpha {\omega }_{0}^{2}\right)\left(\dot{\psi }+{\rm{\Delta }}v{\psi }_{\xi }\right)\right]B=-\sigma \left(1+2\alpha {\omega }_{0}^{2}\right){\left|B\right|}^{2}B,\end{eqnarray} \tag{ 64 }$$

where the prime indicates differentiation with respect to the new variable ζ, the over-dot implies differentiation with respect to time t and the ξ-subscript denotes differentiation with respect to the ξ-variable.

Let us now explore the travelling wave solutions possessed by the nonlinear wave equation (60). First, we consider the simple case, where $\psi \left(\xi ,t\right)=\mu \zeta$. Moreover, μ is an arbitrary constant and without loss of generality, we can also assume that λ = 1. The phase equation (63) simply yields a relation between the known constants and the ones characterizing the travelling wave ansatz, which need to be determined additionally. Thus, we have

$$\begin{eqnarray}&&2\mu -{{\mathscr{D}}}_{\omega }{\gamma }_{u}^{2}\left(1+2\alpha {\omega }_{0}^{2}\right)\left(u-{\rm{\Delta }}v\right)=0,\end{eqnarray} \tag{ 65 }$$

where

$$\begin{eqnarray}&&{\gamma }_{u}=\displaystyle \frac{1}{\sqrt{1-\tfrac{{u}^{2}}{{c}_{1}^{2}}}},\end{eqnarray} \tag{ 66 }$$

while the amplitude equation (64) can be rewritten as

$$\begin{eqnarray}&&{B}^{{\prime\prime} }+{\mu }^{2}B=-{\sigma }_{1}{B}^{3},\qquad \quad {\sigma }_{1}=\sigma {\gamma }_{u}^{2}\left(1+2\alpha {\omega }_{0}^{2}\right).\end{eqnarray} \tag{ 67 }$$

The above equation is the well-known Duffing equation, which in the case of positive σ₁ (which is the case, since σ > 0) has a simple solution expressed in terms of Jacobi cosine function. Assuming that ${B}^{{\prime} }\left(\zeta =0\right)=0$, the exact solution of equation (67) can be expressed as

$$\begin{eqnarray}&&B={{\mathscr{B}}}_{0}\mathrm{cn}\left[\nu \left(\zeta +{\psi }_{0}\right),{k}_{e}\right].\end{eqnarray} \tag{ 68 }$$

Here $\mathrm{cn}\left(w\right)$ denotes the elliptic Jacobi cosine function. In addition, ${{\mathscr{B}}}_{0}=B\left(\zeta =0\right)$, the frequency ν and the elliptic modulus k_e are expressed as follows

$$\begin{eqnarray}&&\nu =\sqrt{{\mu }^{2}+{\sigma }_{1}{{\mathscr{B}}}_{0}^{2}},\qquad {k}_{e}=\displaystyle \frac{{{\mathscr{B}}}_{0}}{\sqrt{2\left({{\mathscr{B}}}_{0}^{2}+\tfrac{{\mu }^{2}}{{\sigma }_{1}}\right)}}.\end{eqnarray} \tag{ 69 }$$

and ψ₀ is an arbitrary phase. Note that ν ≈ μ and for that reason the arbitrary parameter μ can be interpreted as a harmonic number. Figures 1 and 2 show the typical behaviour of the nonlinear quasi-periodic wave for harmonic numbers μ = 1 and μ = 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, we consider the general case of travelling wave solutions to the nonlinear wave equation, for which $\psi =\psi \left(\zeta \right)$. The phase equation (63) can be written as

$$\begin{eqnarray}&&{\psi }^{{\prime\prime} }B+2{\psi }^{{\prime} }{B}^{{\prime} }-\gamma {B}^{{\prime} }=0,\end{eqnarray} \tag{ 70 }$$

where

$$\begin{eqnarray}&&\gamma ={{\mathscr{D}}}_{\omega }{\gamma }_{u}^{2}\left(1+2\alpha {\omega }_{0}^{2}\right)\left(u-{\rm{\Delta }}v\right).\end{eqnarray} \tag{ 71 }$$

It possesses a first integral of the form

$$\begin{eqnarray}&&C={B}^{2}\left({\psi }^{{\prime} }-\displaystyle \frac{\gamma }{2}\right),\end{eqnarray} \tag{ 72 }$$

which after being substituted into the amplitude equation (64) yields the Ermakov-Pinney equation

$$\begin{eqnarray}&&{B}^{{\prime\prime} }+\displaystyle \frac{{\gamma }^{2}}{4}B+{\sigma }_{1}{B}^{3}=\displaystyle \frac{{C}^{2}}{{B}^{3}}.\end{eqnarray} \tag{ 73 }$$

Multiplying both sides of equation (73) by ${B}^{{\prime} }$ and integrating once, we obtain the first integral of the Ermakov-Pinney equation

$$\begin{eqnarray}&&W={B}^{{\prime} 2}+\displaystyle \frac{{C}^{2}}{{B}^{2}}+\displaystyle \frac{{\gamma }^{2}}{4}{B}^{2}+\displaystyle \frac{{\sigma }_{1}}{2}{B}^{4}.\end{eqnarray} \tag{ 74 }$$

Similar to the case worked out above, we assume that ${B}^{{\prime} }\left(\zeta =0\right)=0$. Then, it can be verified by direct substitution that the Ermakov-Pinney equation (73) possesses an exact solution in the form of a cnoidal wave

$$\begin{eqnarray}&&B=\sqrt{\displaystyle \frac{{X}_{2}-{k}_{e}^{2}{X}_{3}{\mathrm{sn}}^{2}\left[\nu \left(\zeta +{\psi }_{0}\right),{k}_{e}\right]}{1-{k}_{e}^{2}{\mathrm{sn}}^{2}\left[\nu \left(\zeta +{\psi }_{0}\right),{k}_{e}\right]}}.\end{eqnarray} \tag{ 75 }$$

The frequency ν and the elliptic modulus k_e are given by the expressions

$$\begin{eqnarray}&&\nu =\sqrt{\displaystyle \frac{{\sigma }_{1}\left({{\mathscr{B}}}_{0}^{2}-{X}_{3}\right)}{2}},\qquad {k}_{e}=\sqrt{\displaystyle \frac{{{\mathscr{B}}}_{0}^{2}-{X}_{2}}{{{\mathscr{B}}}_{0}^{2}-{X}_{3}}}.\end{eqnarray} \tag{ 76 }$$

The constants X₂ and X₃ entering the equations above can be expressed as follows

$$\begin{eqnarray}&&{X}_{\mathrm{2,3}}=-\displaystyle \frac{1}{2}\left[{{\mathscr{B}}}_{0}^{2}+\displaystyle \frac{{\gamma }^{2}}{2{\sigma }_{1}}\mp \sqrt{{\left({{\mathscr{B}}}_{0}^{2}+\displaystyle \frac{{\gamma }^{2}}{2{\sigma }_{1}}\right)}^{2}+\displaystyle \frac{8{C}^{2}}{{\sigma }_{1}{{\mathscr{B}}}_{0}^{2}}}\right].\end{eqnarray} \tag{ 77 }$$

The analogy with the surface gravity waves on shallow water is quite interesting and impressive. Despite the fact that cnoidal waves have been proposed by Korteweg and de-Vries in 1895 and later modified by Benjamin, Bona and Mahony in 1972 in a different context [41, 42], it is intriguing that they can play an important role in wave propagation in polarized quantum vacuum. Another intriguing feature of the cnoidal waves is that their spatio-temporal structure does not change much in terms of spatial and temporal scales with increasing the intensity of the pumping wave. Changes only the corresponding amplitude, which is to be expected. This important characteristics is shown in figures 3 and 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Let us now turn back to equations (55) and (56) and perform the change of variables defined by equation (57). This time we define the translational velocity as

$$\begin{eqnarray}&&{\bf{v}}=\displaystyle \frac{\beta {\omega }_{0}^{2}}{{\rm{\Gamma }}+\beta {\omega }_{0}^{2}}c{\widehat{{\bf{k}}}}_{0}.\end{eqnarray} \tag{ 78 }$$

In the new variables the coupled nonlinear equations (55) and (56) can be written as

$$\begin{eqnarray}&&\begin{array}{l}i{{\mathscr{D}}}_{M\omega }\left({\partial }_{t}+{\rm{\Delta }}{v}_{M}{\partial }_{\xi }\right){{\mathscr{B}}}_{+}+\displaystyle \frac{{{\rm{\Gamma }}}^{2}}{{\rm{\Gamma }}+\beta {\omega }_{0}^{2}}{\widehat{\square }}_{M}{{\mathscr{B}}}_{+}=-24\kappa {\epsilon }_{0}{c}^{2}{\square }_{k}^{2}{\left|{{\mathscr{B}}}_{-}\right|}^{2}{{\mathscr{B}}}_{+}\\ \quad +\ 2\kappa {\epsilon }_{0}{c}^{2}\displaystyle \frac{{\square }_{k}}{{\square }_{0k}}\left(4{\square }_{k}^{2}-6{\square }_{k}{\square }_{0k}+{\square }_{0k}^{2}\right){\left|{{\mathscr{B}}}_{+}\right|}^{2}{{\mathscr{B}}}_{+},\end{array}\end{eqnarray} \tag{ 79 }$$

$$\begin{eqnarray}&&\begin{array}{l}-i{{\mathscr{D}}}_{M\omega }\left({\partial }_{t}+{\rm{\Delta }}{v}_{M}{\partial }_{\xi }\right){{\mathscr{B}}}_{-}+\displaystyle \frac{{{\rm{\Gamma }}}^{2}}{{\rm{\Gamma }}+\beta {\omega }_{0}^{2}}{\widehat{\square }}_{M}{{\mathscr{B}}}_{-}\\ \quad =\,-24\kappa {\epsilon }_{0}{c}^{2}{\square }_{k}^{2}{\left|{{\mathscr{B}}}_{+}\right|}^{2}{{\mathscr{B}}}_{-}-2\kappa {\epsilon }_{0}{c}^{2}\displaystyle \frac{{\square }_{k}}{{\square }_{0k}}\left(4{\square }_{k}^{2}+6{\square }_{k}{\square }_{0k}+{\square }_{0k}^{2}\right){\left|{{\mathscr{B}}}_{-}\right|}^{2}{{\mathscr{B}}}_{-},\end{array}\end{eqnarray} \tag{ 80 }$$

where now

$$\begin{eqnarray}&&{\widehat{\square }}_{M}={{\boldsymbol{\nabla }}}_{{\boldsymbol{\xi }}}^{2}-\displaystyle \frac{1}{{c}_{M}^{2}}{\partial }_{t}^{2},\qquad \qquad {c}_{M}=\displaystyle \frac{{\rm{\Gamma }}c}{{\rm{\Gamma }}+\beta {\omega }_{0}^{2}},\end{eqnarray} \tag{ 81 }$$

and the other notations $\left({{\mathscr{D}}}_{M\omega },\ {{\bf{v}}}_{{gM}}\ \mathrm{and}\ {\rm{\Delta }}{{\bf{v}}}_{M}\right)$ are obvious and self-explanatory, in analogy to the case of polarized vacuum without external field, which has been just considered. Equations (79) and (80) possess a simple solution in the form of a plane wave $\left[{{\mathscr{B}}}_{\pm }={B}_{\pm }{{\rm{e}}}^{i\left(K\xi -{{\rm{\Omega }}}_{\pm }t\right)}\right]$ with constant complex amplitudes B_± . The two frequencies ${{\rm{\Omega }}}_{\pm }$ are solutions to the dispersion equations

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Gamma }}}^{2}}{{\rm{\Gamma }}+\beta {\omega }_{0}^{2}}{\square }_{\pm }\pm {{\mathscr{D}}}_{M\omega }\left({{\rm{\Omega }}}_{\pm }-{\rm{\Delta }}{v}_{M}K\right)={{\rm{\Sigma }}}_{\pm },\end{eqnarray} \tag{ 82 }$$

where

$$\begin{eqnarray}&&{\square }_{\pm }=\displaystyle \frac{{{\rm{\Omega }}}_{\pm }^{2}}{{c}_{M}^{2}}-{K}^{2},\end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\pm }=-24\kappa {\epsilon }_{0}{c}^{2}{\square }_{k}^{2}{\left|{B}_{\mp }\right|}^{2}\pm 2\kappa {\epsilon }_{0}{c}^{2}\displaystyle \frac{{\square }_{k}}{{\square }_{0k}}\left(4{\square }_{k}^{2}\mp 6{\square }_{k}{\square }_{0k}+{\square }_{0k}^{2}\right){\left|{B}_{\pm }\right|}^{2}.\end{eqnarray} \tag{ 84 }$$

Since in the general nontrivial case both ${{\rm{\Sigma }}}_{\pm }$ cannot be equal to zero, we would obtain ${{\rm{\Omega }}}_{-}\left(-K\right)\ne -{{\rm{\Omega }}}_{+}\left(K\right)$, that is ${B}_{-}\ne {B}_{+}^{* }$. This implies that the two wave modes ${{\mathscr{B}}}_{\pm }$ are independent and they propagate at different frequencies. It is well-known that birefringence makes use of the refractive index of the medium, defined as ${n}_{\pm }={c}_{M}K/{{\rm{\Omega }}}_{\pm }$. Birefringence Δ n is defined as the difference between the refractive indices characterizing the two polarizations, or

$$\begin{eqnarray}&&{\rm{\Delta }}n={c}_{M}K\left(\displaystyle \frac{1}{{{\rm{\Omega }}}_{+}}-\displaystyle \frac{1}{{{\rm{\Omega }}}_{-}}\right).\end{eqnarray} \tag{ 85 }$$

The effect described above by taking into account that polarized vacuum in external magnetic field can be considered as a nonlinear medium is usually called nonlinear birefringence of vacuum.

Starting from the basics, we have studied the implications of quantum corrections to classical electrodynamics and the propagation of electromagnetic waves and pulses. Used essentially is the property of the polarized vacuum, that the nonlinear wave equation for the electromagnetic vector potential possesses an exact solution in the form of a plane wave in both the case of free from external fields vacuum, as well as when a constant magnetic field is applied.

The initial nonlinear wave equation is solved perturbatively about the known exact plane wave solution. Following the elegant approach of the proto RG operator, we derive a nonlinear wave equation with a nonzero convective part (containing first-order spatial derivative with respect to the dependent variable) for the (relatively) slowly varying amplitude of the first-order perturbation. This equation governs the propagation of electromagnetic waves with a reduced speed of light. The reduction is roughly proportional to the intensity of the pumping plane wave.

An external background magnetic field introduces spatial anisotropy coercing electromagnetic pulses to travel in the direction of the applied field. In analogy to the case of polarized vacuum without external field, the initial wave equation for the electromagnetic vector potential has been analysed in detail by following the same guidelines. In order to avoid insignificant complications in the specific calculations, we consider here the case of a circularly polarized initial pumping plane wave. The analysis of the general case of an elliptically polarized exact solution is also possible, but with it, the technical details conceal the physical nature of the problem. Repeating the proto RG procedure in a similar manner as has been done in the case of a quantum vacuum without externally applied magnetic field, we obtain a system of coupled nonlinear wave equations for the two slowly varying amplitudes of the first-order perturbation solution, which describe the two polarization states.

The stationary (traveling) wave solutions of the proposed nonlinear wave equation with nonzero convective part for the case of a polarized vacuum without external fields have been obtained. The slowly varying wave amplitude behaviour is shown to be similar to that of a cnoidal wave, known to describe surface gravity waves in shallow water. It has been demonstrated that the two wave modes describing the two polarization states are independent, and they propagate at different wave frequencies. This effect is usually called nonlinear birefringence.

It is worthwhile to be pointed out that our considerations are based on the leading effective interaction of classical electromagnetic fields, which give rise to nonlinear terms in the polarization P and the magnetization M vectors proportional to the cube of E and B (involving the well-known electromagnetic invariants). Specifically, in one-loop Heisenberg-Euler theory this effective interaction scales as ${\alpha }_{{fs}}^{2}$ in terms of the fine structure constant. Sub leading corrections arising from higher-loop contributions and scaling as ${\alpha }_{{fs}}^{3}$ give rise to quintic (in E and B) terms in the corresponding expressions for P and M (see e.g. the recent [43]). Limiting oneself to the effective interaction of second order in the fine structure constant, terms of the order ${\alpha }_{{fs}}^{3}$ are obviously not accounted for from the outset. Correspondingly, all the results derived within this truncation can be reliably and consistently accounted for in the fundamental physics only up to terms of this order in α_fs.

We wish to express our gratitude to Prof SV Bulanov and Dr H Kadlecova for many enlightening discussions and suggestions related to the topics touched upon in the article. We are especially indebted to the unknown reviewer for their interesting comments and valuable suggestions aimed at improving the quality of the exposition in the present article.

This work has been supported by Extreme Light Infrastructure—Nuclear Physics (ELI-NP) Phase II, an innovative project co-financed by the Romanian Government and the European Union through the European Regional Development Fund.

We consider the rapidly oscillating term on the right-hand-side of equation (22) as perturbation and rewrite the latter accordingly

$$\begin{eqnarray}&&{\widehat{{\bf{L}}}}_{1}{A}_{1}=-\epsilon \alpha \left[{{\rm{e}}}^{2i\left({\varphi }_{0}+\delta \right)}+{{\rm{e}}}^{-2i\left({\varphi }_{0}+\delta \right)}\right]{\widehat{{\bf{D}}}}_{1}^{2}{A}_{1},\end{eqnarray} \tag{ A1 }$$

where, once again, is a formal small parameter introduced to measure successive orders of magnitude in the perturbation expansion. Similar to equation (21), we expand A₁ according to the expression

$$\begin{eqnarray}&&{A}_{1}\left(x,z;t\right)=\sum _{n=0}^{\infty }{\epsilon }^{n}{A}_{1}^{(n)}\left(x,z;t\right).\end{eqnarray} \tag{ A2 }$$

Obviously, the zero-order term ${A}_{1}^{(0)}$ is just the expression on the right-hand-side of equation (25).

The first-order perturbation equation acquires the form

$$\begin{eqnarray}&&{\widehat{{\bf{L}}}}_{1}{A}_{1}^{(1)}=-\displaystyle \frac{{\square }_{k}}{2}\left[{\mathscr{A}}{{\rm{e}}}^{i\left({\varphi }_{1}+2{\varphi }_{0}+2\delta \right)}+{\mathscr{A}}{{\rm{e}}}^{i\left({\varphi }_{1}-2{\varphi }_{0}-2\delta \right)}+c.c.\right].\end{eqnarray} \tag{ A3 }$$

Its solution can be written as

$$\begin{eqnarray}&&{A}_{1}^{(1)}=-\displaystyle \frac{{\square }_{k}{\mathscr{A}}}{8{\square }_{0k}}\left[{{\rm{e}}}^{i\left({\varphi }_{1}+2{\varphi }_{0}+2\delta \right)}-{{\rm{e}}}^{i\left({\varphi }_{1}-2{\varphi }_{0}-2\delta \right)}\right]+c.c..\end{eqnarray} \tag{ A4 }$$

Finally, the second-order perturbation equation can be written as

$$\begin{eqnarray}&&{\widehat{{\bf{L}}}}_{1}{A}_{1}^{(2)}=\displaystyle \frac{{\square }_{k}^{2}}{16{\square }_{0k}}\left[{{\rm{e}}}^{2i\left({\varphi }_{0}+\delta \right)}+{{\rm{e}}}^{-2i\left({\varphi }_{0}+\delta \right)}\right]\left[{\mathscr{A}}{{\rm{e}}}^{i\left({\varphi }_{1}+2{\varphi }_{0}+2\delta \right)}-{\mathscr{A}}{{\rm{e}}}^{i\left({\varphi }_{1}-2{\varphi }_{0}-2\delta \right)}+c.c.\right].\end{eqnarray} \tag{ A5 }$$

It is straightforward to verify that resonant terms are missing on the right-hand-side of the above equation. Its solution contains only regular terms proportional to ${{\rm{e}}}^{i\left({\varphi }_{1}\pm 4{\varphi }_{0}\pm 4\delta \right)}$ and their complex conjugated counterpart. This implies that the dispersion relation (27) remains unaltered by the fast oscillating term in question, which completes the proof concerning the justification of the averaging procedure pursued in section 3.1.