2. Formulation of general equation for beam of electromagnetic wave propagation in nonlinear media

As is well known, the electric and magnetic fields in dielectric media are described by the macroscopic Maxwell’s equations as follows:

(1)$$\begin{eqnarray}\displaystyle {\rm\nabla}\boldsymbol{\cdot }\vec{D} & = & \displaystyle 0,\end{eqnarray}$$

(2)$$\begin{eqnarray}\displaystyle {\rm\nabla}\boldsymbol{\cdot }\vec{B} & = & \displaystyle 0,\end{eqnarray}$$

(3)$$\begin{eqnarray}\displaystyle {\rm\nabla}\times \vec{E} & = & \displaystyle -\frac{\partial \vec{B}}{\partial t},\end{eqnarray}$$

(4)$$\begin{eqnarray}\displaystyle {\rm\nabla}\times \vec{H} & = & \displaystyle \frac{\partial \vec{D}}{\partial t},\end{eqnarray}$$

$\vec{E}$

$\vec{B}$

$\vec{D}$

$\vec{H}$

$\vec{E}$

$\vec{B}$

(5)$$\begin{eqnarray}\displaystyle \vec{D} & = & \displaystyle {\it\varepsilon}_{0}{\it\varepsilon}\vec{E},\end{eqnarray}$$

(6)$$\begin{eqnarray}\displaystyle \vec{B} & = & \displaystyle {\it\mu}_{0}{\it\mu}\vec{H},\end{eqnarray}$$

${\it\varepsilon}_{0}$

${\it\mu}_{0}$

${\it\varepsilon}$

${\it\mu}$

Following the known procedure we take the curl vector operator of both sides of Equation (3) and the time derivative of both sides of Equation (4) and assume that the magnetic effect of the media is negligible, i.e., ${\it\mu}=1$. Then, using Equations (5) and (6), we can eliminate the magnetic and magnetization fields, obtaining the equation for the electric field

(7)$$\begin{eqnarray}{\rm\nabla}\times \!({\rm\nabla}\times \vec{E})+{\it\varepsilon}_{0}{\it\mu}_{0}\frac{\partial ^{2}({\it\varepsilon}\vec{E})}{\partial t^{2}}=0.\end{eqnarray}$$

Assuming that the material effect on the electric field is slow compared to the period of optical oscillation or the laser pulse duration, i.e., assuming time independence of the permittivity of material, assuming that the material is weakly absorbing, and recalling that the speed of the electromagnetic wave in a vacuum is $c=1/\sqrt{{\it\varepsilon}_{0}{\it\mu}_{0}}$ and the index of refraction of the material is $n=\sqrt{{\it\varepsilon}}$, we rewrite Equation (7) as

(8)$$\begin{eqnarray}{\rm\nabla}\times ({\rm\nabla}\times \vec{E})+\frac{n^{2}}{c^{2}}\frac{\partial ^{2}\vec{E}}{\partial t^{2}}=0.\end{eqnarray}$$

Now, using the identity ${\rm\nabla}\times ({\rm\nabla}\times \vec{A})={\rm\nabla}({\rm\nabla}\boldsymbol{\cdot }\vec{A})-{\rm\Delta}\vec{A}$ we rewrite Equation (8) in the following form:

(9)$$\begin{eqnarray}{\rm\Delta}\vec{E}-\frac{n^{2}}{c^{2}}\frac{\partial ^{2}\vec{E}}{\partial t^{2}}-{\rm\nabla}({\rm\nabla}\boldsymbol{\cdot }\vec{E})=0.\end{eqnarray}$$

From this point the derivations will significantly deviate from the procedure commonly performed in all of the books and journal publications on nonlinear optics since we will be considering the permittivity of the material and, consequently, the index of refraction, as a coordinate dependent function. In contrast to our approach, customary considerations treat permittivity as a constant, zeroing the third term in the left-hand side of Equation (9) and reducing this equation to the commonly known form that contains only two first terms and is called the wave equation. As we show below, neglecting the third term is a significant mistake that leads to an inadequate description of wave propagation in nonlinear media.

For the conditions of inhomogeneity of the properties of electrically neutral media, such as in the case of propagation of short wavelength waves in the ionosphere^{[Reference Budden13]}, Equation (9) can be rewritten in the following form:

(10)$$\begin{eqnarray}{\rm\Delta}\vec{E}-\frac{n^{2}}{c^{2}}\frac{\partial ^{2}\vec{E}}{\partial t^{2}}+{\rm\nabla}\left(\frac{{\rm\nabla}{\it\varepsilon}\boldsymbol{\cdot }\vec{E}}{{\it\varepsilon}}\right)=0,\end{eqnarray}$$

where we use Equations (1) and (5) from which it follows that ${\rm\nabla}{\it\varepsilon}\boldsymbol{\cdot }\vec{E}+{\it\varepsilon}{\rm\nabla}\boldsymbol{\cdot }\vec{E}=0$ and, therefore, ${\rm\nabla}\boldsymbol{\cdot }\vec{E}=-({\rm\nabla}{\it\varepsilon}\boldsymbol{\cdot }\vec{E}/{\it\varepsilon})$. Finally, recalling that ${\it\varepsilon}=n^{2}$ and using the convenient expression for the displacement field $\vec{D}={\it\varepsilon}_{0}{\it\varepsilon}\vec{E}={\it\varepsilon}_{0}(1+{\it\chi})\vec{E}={\it\varepsilon}_{0}\vec{E}+\vec{P}$, where ${\it\chi}$ is the electric susceptibility and $\vec{P}$ is the polarization density, we rewrite Equation (7) in the following form:

(11)$$\begin{eqnarray}{\rm\Delta}\vec{E}-\frac{1}{c^{2}}\frac{\partial ^{2}\vec{E}}{\partial t^{2}}+2{\rm\nabla}\left(\frac{{\rm\nabla}n\boldsymbol{\cdot }\vec{E}}{n}\right)={\it\mu}_{0}\frac{\partial ^{2}\vec{P}_{L}}{\partial t^{2}}+{\it\mu}_{0}\frac{\partial ^{2}\vec{P}_{NL}}{\partial t^{2}},\end{eqnarray}$$

where the polarization density vector is represented by the sum of linear and nonlinear components denoted using the subscripts ‘$L$’ and ‘$NL$’, respectively.

The solution of the propagation equation expressed in the form of either Equation (10) or (11) can be expressed in terms of the slowly varying amplitude function

(12)$$\begin{eqnarray}\vec{E}(\vec{r},t)=\vec{A}(\vec{r})e^{\text{i}(\vec{k}(\vec{r})\boldsymbol{\cdot }\vec{r}-{\it\omega}t)}+\text{c.c.}\end{eqnarray}$$

In this solution the vector amplitudes of the electric field and the wavevector are coordinate dependent, i.e., expression (12) represents the electric field of a nonplanar wave.

Below we will demonstrate that, within paraxial approximation and considering the propagation range in which variation of the spatial profile of laser beam irradiance due to the effect of nonlinear induced refraction is small and, thus, can be considered as perturbation, the propagation Equations (10) or (11) can be straightforwardly modified into an equation that, as proposed in Refs. [Reference Semak and Shneider11, Reference Semak and Shneider12], has a solution represented by the blending of the solution of the Helmholtz equation for propagation of a laser beam in a medium with a uniform and irradiance independent refractive index similar to the one obtained by Kogelnik and Li^{[Reference Kogelnik and Li14]} and a correction term that represents nonlinear field perturbation expressed in terms of paraxial ray optics (the eikonal equation)^{[Reference Born and Wolf15]}.

3. Examination of current theory for laser beam propagation in nonlinear media

Now we will demonstrate that the current formulations of the propagation equation in nonlinear (self-induced inhomogeneity) media are based on two inconsistent and self-contradictory assumptions. First, in all theoretical works known to us it is assumed that the laser beam has a plane wavefront. Second, as mentioned above, the term responsible for refraction due to media inhomogeneity is disregarded, i.e., the media are assumed as homogeneous.

Under the first assumption, i.e., assumption of a plane wavefront, the solution of the propagation equation is expressed in the form of the slowly varying amplitude function

(13)$$\begin{eqnarray}\displaystyle \vec{E}(x,y,z,t) & = & \displaystyle A_{x}(x,y,z)e^{\text{i}(k_{z}z-{\it\omega}t)}\hat{x}\nonumber\\ \displaystyle & & \displaystyle +\,A_{y}(x,y,z)e^{\text{i}(k_{z}z-{\it\omega}t)}{\hat{y}}+\text{c.c.}\end{eqnarray}$$

Note that here, in accordance with the plane wave assumption, the scalar product of vectors $\vec{k}\boldsymbol{\cdot }\vec{r}$ from Equation (12) is substituted with the product of scalars – $k_{z}z$, where $k_{z}$ is the component of the wavevector along the $z$-axis. The solution form of Equation (13) assumes that the $x$- and $y$-components of the wavevector and $z$-component of the electric field are zero, i.e., the beam propagates exactly along the $z$-axis.

According to the second assumption, the term describing refraction on the gradient of the refractive index is neglected. Then, assuming linear polarization and, since the requirement of slowly varying amplitude implies that

(14)$$\begin{eqnarray}\left|\frac{\partial ^{2}A}{\partial z^{2}}\right|\ll \left|k\frac{\partial A}{\partial z}\right|,\end{eqnarray}$$

Equation (10) can be reformulated in a form retaining only the amplitude $A(x,y,z)$ of the electric vector aligned along either the $x$-axis or the $y$-axis,

(15)$$\begin{eqnarray}{\rm\Delta}_{T}A+\left(\frac{{\it\omega}^{2}n^{2}}{c^{2}}-k_{z}^{2}\right)A+2\text{i}k_{z}\frac{\partial A}{\partial z}=0,\end{eqnarray}$$

where ${\rm\Delta}_{T}$ is the transverse Laplace operator.

By definition, the wavevector $k={\it\omega}n/c$, and since a plane wave is assumed, the second term in Equation (15) must be zero. However, all textbooks and scientific articles at this stage of consideration submit that the wavefront deviates from a plane. This, of course, contradicts the initial assumption of a strictly plane front; however, is necessary, as otherwise the self-focusing and self-trapping would not follow from Equation (15). Thus, the second term of propagation Equation (15) in all current models is assumed to be nonzero. We will show below that this manipulation has a devastating consequence.

Another inconsistent assumption of the current theory of electromagnetic wave propagation in nonlinear media is that the third term in complete propagation Equation (11) is negligible and can be omitted, leading to the commonly used Equation (15). As far as we know a justification for such an omission was never provided. It appears that the origin of this assumption can be traced to the original work^{[Reference Chiao, Garmire and Townes16]}. All of the subsequent works, except for our recent publications^{[Reference Semak and Shneider11, Reference Semak and Shneider12]}, followed the path laid^{[Reference Chiao, Garmire and Townes16]} and labored on various mathematical treatments of equations that can be traced to the propagation equation deduced for homogeneous media. Thus, the ‘foundational’ propagation equation in nonlinear optics was deduced while disregarding media inhomogeneity (inherent, induced or self-induced), and has the form of Equation (15) missing the ‘refraction’ term (see for example, Equation 7.2.9 in Ref. [Reference Boyd3]).

At this point one should wonder how an equation with a term missing the refraction due to media inhomogeneity can describe self-focusing. The answer is hidden in the second term of Equation (15). The ‘nonzeroing’ of the second term in Equation (15) is crucial for constructing all nonlinear propagation effects out of this oversimplified equation. Indeed, introducing nonlinearity of the refractive index, $n=n_{0}+n_{2}\langle |E|^{2}\rangle _{t}$, Equation (15) can be modified into the following form^{[Reference Talanov17, Reference Kelley18]}:

(16)$$\begin{eqnarray}\frac{\partial ^{2}E}{\partial x^{2}}+2\text{i}k\frac{\partial E}{\partial z}=-k^{2}\frac{n_{2}}{n_{0}}|E|^{2}E,\end{eqnarray}$$

that is equivalent to the infamous (in the realm of nonlinear wave propagation) nonlinear Schrodinger’s equation^{[Reference Zakharov and Shabat19]}

(17)$$\begin{eqnarray}\frac{\partial ^{2}{\it\psi}}{\partial x^{2}}+\text{i}\frac{\partial {\it\psi}}{\partial z}=-{\it\kappa}|{\it\psi}|^{2}{\it\psi}.\end{eqnarray}$$

4. Revised propagation equation for slowly varying amplitude

Let us revise the foundation of the current theory for the propagation of a beam of an electromagnetic wave in a nonlinear medium. In this consideration we will use the general form of 3D propagation Equation (7) that follows directly from Maxwell’s equations. Following a traditional approach, we will look for the solution of this equation in the following form of wave with varying amplitude:

(18)$$\begin{eqnarray}\displaystyle \vec{E} & = & \displaystyle E_{x}\hat{x}+E_{y}{\hat{y}}+E_{z}\hat{z}=\vec{A}(x,y,z)p(t)e^{\text{i}(k_{x}\hat{x}+k_{y}{\hat{y}}+k_{z}\hat{z}-{\it\omega}t)}\nonumber\\ \displaystyle & = & \displaystyle (A_{x}(x,y,z)\hat{x}+A_{y}(x,y,z){\hat{y}}+A_{z}(x,y,z)\hat{z})p(t)\nonumber\\ \displaystyle & & \displaystyle \times \,e^{\text{i}(k_{x}\hat{x}+k_{y}{\hat{y}}+k_{z}\hat{z}-{\it\omega}t)},\end{eqnarray}$$

where function $p(t)$ is the dimensionless pulse shape such that the time integral of this function from minus to plus infinity equals unity.

Let us now modify propagation Equation (9) assuming that the pulse shape is a slowly varying function compared to the period of wave oscillation and that the dielectric constant, ${\it\varepsilon}$, and therefore the refractive index, $n$, are time independent. Then, recalling that ${\it\varepsilon}=n^{2}$, and substituting solution of Eqaution (18) into Equation (9), neglecting the time derivative of pulse shape $p(t)$ as it is a slow function, and eliminating the exponent of phase, Equation (9) can be transformed into the equation for the amplitude,

(19)$$\begin{eqnarray}\displaystyle & & \displaystyle \left(k_{x}^{2}+k_{y}^{2}+k_{z}^{2}-\frac{{\it\omega}^{2}n^{2}}{c^{2}}\right)A_{x}\hat{x}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\left(k_{x}^{2}+k_{y}^{2}+k_{z}^{2}-\frac{{\it\omega}^{2}n^{2}}{c^{2}}\right)A_{y}{\hat{y}}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\left(k_{x}^{2}+k_{y}^{2}+k_{z}^{2}-\frac{{\it\omega}^{2}n^{2}}{c^{2}}\right)A_{z}\hat{z}\nonumber\\ \displaystyle & & \displaystyle \quad -\,2\text{i}\left(k_{x}\frac{\partial A_{x}}{\partial x}+k_{y}\frac{\partial A_{x}}{\partial y}+k_{z}\frac{\partial A_{x}}{\partial z}\right)\hat{x}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\left(\frac{\partial ^{2}A_{x}}{\partial x^{2}}+\frac{\partial ^{2}A_{x}}{\partial y^{2}}+\frac{\partial ^{2}A_{x}}{\partial z^{2}}\right)\hat{x}\nonumber\\ \displaystyle & & \displaystyle \quad -\,2\text{i}\left(k_{x}\frac{\partial A_{y}}{\partial x}+k_{y}\frac{\partial A_{y}}{\partial y}+k_{z}\frac{\partial A_{y}}{\partial z}\right){\hat{y}}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\left(\frac{\partial ^{2}A_{y}}{\partial x^{2}}+\frac{\partial ^{2}A_{y}}{\partial y^{2}}+\frac{\partial ^{2}A_{y}}{\partial z^{2}}\right){\hat{y}}\nonumber\\ \displaystyle & & \displaystyle \quad -\,2\text{i}\left(k_{x}\frac{\partial A_{z}}{\partial x}+k_{y}\frac{\partial A_{z}}{\partial y}+k_{z}\frac{\partial A_{z}}{\partial z}\right)\hat{z}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\left(\frac{\partial ^{2}A_{z}}{\partial x^{2}}+\frac{\partial ^{2}A_{z}}{\partial y^{2}}+\frac{\partial ^{2}A_{z}}{\partial z^{2}}\right)\hat{z}\nonumber\\ \displaystyle & & \displaystyle \quad -\,2\left[{\rm\nabla}\left(\frac{{\rm\nabla}n}{n}\boldsymbol{\cdot }(\vec{A}e^{\text{i}(k_{x}x+k_{y}y+k_{z}z)})\right)\right]e^{-\text{i}(k_{x}x+k_{y}y+k_{z}z)}\nonumber\\ \displaystyle & & \displaystyle =0.\end{eqnarray}$$

Note that the first three terms in the left-hand side of Equation (19) are null, since by definition $k_{x}^{2}+k_{y}^{2}+k_{z}^{2}-({\it\omega}^{2}n^{2}/c^{2})=0$, and, as in the previous deduction of Equation (10), we used equation ${\rm\nabla}\boldsymbol{\cdot }\vec{E}=-(\vec{E}{\rm\nabla}{\it\varepsilon}/{\it\varepsilon})$ in order to formulate the right-hand side of Equation (19). Thus, the propagation equation in the nonlinear media has the following form:

(20)$$\begin{eqnarray}\displaystyle & & \displaystyle {\rm\Delta}\vec{A}+2\text{i}\left(k_{x}\frac{\partial \vec{A}}{\partial x}+k_{y}\frac{\partial \vec{A}}{\partial y}+k_{z}\frac{\partial \vec{A}}{\partial z}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,2\left[{\rm\nabla}\left(\frac{{\rm\nabla}n}{n}\boldsymbol{\cdot }(\vec{A}e^{\text{i}(k_{x}x+k_{y}y+k_{z}z)})\right)\right]e^{-\text{i}(k_{x}x+k_{y}y+k_{z}z)}\nonumber\\ \displaystyle & & \displaystyle \quad =0.\end{eqnarray}$$

This equation is cardinally different from the currently used Equations (15)–(17) in both aspects of physics and of mathematics. From the point of view of a physicist, Equation (20) straightforwardly shows that beam refraction is produced by self-induced inhomogeneity of the refractive index (that can be a result of the Kerr effect, material ionization etc.). Also, a physicist will find appealing that the presented theory expressed by Equation (20) satisfies the correspondence principle since, as we will demonstrate below, it contains geometrical optics and Equation (20) transforms into the ray-optics equation under the assumption of infinitely small wavelength. In contrast, Equations (15)–(17) do not lead to the geometrical optics in the limiting case of infinitely small wavelength and thus do not satisfy the correspondence principle.

From the point of view of a mathematician new equation of propagation Equation (20) dramatically differs from the current propagation equation in one very peculiar aspect – it has no self-similar solution. In contrast, Equation (16) rewritten as nonlinear Schrodinger’s Equation (17) has a self-similar solution, or so called soliton solution, that serves as the foundation for the prediction of laser beam self-trapping and all the current ‘filament’ theories that predict mind boggling lengths of self-trapped laser ‘filament’.

5. Modification of the revised propagation equation for cases of slowly variable amplitude and paraxial laser beam propagation at distances shorter than or comparable to the Rayleigh length

Let us explore the complete 3D propagation Equation (20) within paraxial beam approximation while considering a propagation range in which perturbation of the spatial profile of laser beam irradiance by the nonlinear induced refraction is negligible. The latter condition is realized within the range of several Rayleigh lengths for a focused laser beam with pulse energy that is below a certain value (see detailed discussion in Refs. [Reference Semak and Shneider11, Reference Semak and Shneider12]).

Here we will demonstrate that the propagation Equation (20) leads, as proposed in Refs. [Reference Semak and Shneider11, Reference Semak and Shneider12], to the solution represented by the blending of the solution of the Helmholtz equation for propagation of a laser beam in a medium with uniform and irradiance independent refractive index similar to the one obtained by Kogelnik and Li^{[Reference Kogelnik and Li14]} and a correction term that represents nonlinear field perturbation expressed in terms of a paraxial ray-optics (eikonal) equation^{[Reference Born and Wolf15]}.

Assuming paraxial beam propagation it is easy to see that one can neglect in propagation Equation (20) the terms containing $x$- and $y$-components of the wavevector and the terms containing the $z$-component of the electric field amplitude as well as their derivatives (see the schematic of the paraxial beam propagation in Figure 1).

Figure 1. Schematic of evolution of electric field amplitude and wavevector during laser beam propagation.

Assuming linear polarization in the $x{-}z$ plane and neglecting smaller terms as described above allows the following simplification of propagation Equation (20):

(21)$$\begin{eqnarray}\displaystyle & & \displaystyle \left[\left(\frac{\partial ^{2}A_{x}}{\partial x^{2}}+\frac{\partial ^{2}A_{x}}{\partial y^{2}}\right)+2\text{i}\left(k_{x}\frac{\partial A_{x}}{\partial x}+k_{z}\frac{\partial A_{x}}{\partial z}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \quad -\,\text{i}\left(\frac{\partial k_{x}}{\partial x}+\frac{\partial k_{z}}{\partial z}\right)A_{x}+2\text{i}\left(x\frac{\partial k_{x}}{\partial x}\frac{\partial A_{x}}{\partial x}+z\frac{\partial k_{z}}{\partial z}\frac{\partial A_{x}}{\partial z}\right)\nonumber\\ \displaystyle & & \displaystyle \quad -\,\text{i}\left(x\frac{\partial ^{2}k_{x}}{\partial x^{2}}+z\frac{\partial ^{2}k_{z}}{\partial z^{2}}\right)\!A_{x}\!+\!2\left(xk_{x}\frac{\partial k_{x}}{\partial x}+zk_{z}\frac{\partial k_{z}}{\partial z}\right)\!A_{x}\nonumber\\ \displaystyle & & \displaystyle \left.\quad +\left(x^{2}\left(\frac{\partial k_{x}}{\partial x}\right)^{2}+z^{2}\left(\frac{\partial k_{z}}{\partial z}\right)^{2}\right)A_{x}\right]\hat{x}\nonumber\\ \displaystyle & & \displaystyle \quad +\left[\left(\frac{\partial ^{2}A_{z}}{\partial x^{2}}+\frac{\partial ^{2}A_{z}}{\partial y^{2}}\right)+2\text{i}\left(k_{x}\frac{\partial A_{z}}{\partial x}+k_{z}\frac{\partial A_{z}}{\partial z}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \quad -\,\text{i}\left(\frac{\partial k_{x}}{\partial x}+\frac{\partial k_{z}}{\partial z}\right)A_{z}+\,2\text{i}\left(x\frac{\partial k_{x}}{\partial x}\frac{\partial A_{z}}{\partial x}+z\frac{\partial k_{z}}{\partial z}\frac{\partial A_{z}}{\partial z}\right)\nonumber\\ \displaystyle & & \displaystyle \quad -\,\text{i}\left(x\frac{\partial ^{2}k_{x}}{\partial x^{2}}+z\frac{\partial ^{2}k_{z}}{\partial z^{2}}\right)\!A_{z}\!+2\left(xk_{x}\frac{\partial k_{x}}{\partial x}+zk_{z}\frac{\partial k_{z}}{\partial z}\!\right)\!A_{z}\nonumber\\ \displaystyle & & \displaystyle \quad \left.+\left(x^{2}\left(\frac{\partial k_{x}}{\partial x}\right)^{2}+z^{2}\left(\frac{\partial k_{z}}{\partial z}\right)^{2}\right)A_{z}\right]\hat{z}\nonumber\\ \displaystyle & & \displaystyle \quad -\,2\left[{\rm\nabla}\left(\frac{{\rm\nabla}n}{n}\boldsymbol{\cdot }(\vec{A}e^{\text{i}(k_{x}x+k_{y}y+k_{z}z)})\right)\right]e^{-\text{i}(k_{x}x+k_{y}y+k_{z}z)}\nonumber\\ \displaystyle & {\approx} & \displaystyle \left[\left(\frac{\partial ^{2}A_{x}}{\partial x^{2}}+\frac{\partial ^{2}A_{x}}{\partial y^{2}}\right)+2\text{i}k_{z}\frac{\partial A_{x}}{\partial z}\right]\hat{x}+\left[2\text{i}k_{z}\frac{\partial A_{z}}{\partial z}\right]\hat{z}\nonumber\\ \displaystyle & & \displaystyle \quad -\,2\left[{\rm\nabla}\left(\frac{{\rm\nabla}n}{n}\boldsymbol{\cdot }(\vec{A}e^{\text{i}(k_{x}x+k_{y}y+k_{z}z)})\right)\right]e^{-\text{i}(k_{x}x+k_{y}y+k_{z}z)}\nonumber\\ \displaystyle & = & \displaystyle 0.\end{eqnarray}$$

Then, from the simplified propagation Equation (21) we can extract two equations: one for the $x$-coordinate

(22)$$\begin{eqnarray}\displaystyle & & \displaystyle \left(\frac{\partial ^{2}A_{x}}{\partial x^{2}}+\frac{\partial ^{2}A_{x}}{\partial y^{2}}\right)+2\text{i}k_{z}\frac{\partial A_{x}}{\partial z}\nonumber\\ \displaystyle & & \displaystyle \quad +\,2\left[{\rm\nabla}_{x}\left(\frac{{\rm\nabla}n}{n}\boldsymbol{\cdot }\vec{A}e^{\text{i}(k_{x}x+k_{y}y+k_{z}z)}\right)\right]e^{-\text{i}(k_{x}x+k_{y}y+k_{z}z)}\nonumber\\ \displaystyle & & \displaystyle =0,\end{eqnarray}$$

and the second for the $z$-coordinate

(23)$$\begin{eqnarray}\displaystyle & & \displaystyle \text{i}k_{z}\frac{\partial A_{z}}{\partial z}+\left[{\rm\nabla}_{z}\left(\frac{{\rm\nabla}n}{n}\boldsymbol{\cdot }(\vec{A}e^{\text{i}(k_{x}x+k_{y}y+k_{z}z)})\!\right)\right]e^{-\text{i}(k_{x}x+k_{y}y+k_{z}z)}\nonumber\\ \displaystyle & & \displaystyle \quad =0.\end{eqnarray}$$

The third term in the left-hand side of Equation (22) for the $x$- (transverse) component of the amplitude of the electric field is small compared to the first two terms, and for the short propagation distances it can be neglected. Then, Equation (22) acquires a form similar to the equation obtained for diffraction dominated laser beam propagation in empty resonators^{[Reference Kogelnik and Li14]}. The solutions of Equation (22) while neglecting the third term can be found in the classical article of Kogelnik and Li^{[Reference Kogelnik and Li14]}, and according to this work, the fundamental mode of the solution is represented by the field that has Gaussian distribution of amplitude in the radial direction with a beam width that changes along the $z$-axis and has a spherical shape of the wavefront with a radius that is also a function of $z$. Of course, for far field propagation the third term in the left-hand side of Equation (22) must be accounted for, since the relatively small deviations of the amplitude of the electric field and the shape of the wavefront should accumulate while propagating at long distances, resulting in significant modification of both the beam intensity distribution and the wavefront shape.

Both terms in the left-hand side of Equation (23) have similar magnitude. Now, let’s demonstrate that the solution of Equation (23) describes perturbation of the wavefront of the ‘carrier’ field given by the approximate solution of Equation (22) described above.

(24)$$\begin{eqnarray}\displaystyle \text{i}k_{z}\frac{\partial A_{z}}{\partial z}\hat{z} & = & \displaystyle -\left[{\rm\nabla}_{z}\left(\frac{{\rm\nabla}n}{n}\boldsymbol{\cdot }(\vec{A}e^{\text{i}(k_{x}x+k_{y}y+k_{z}z)})\right)\right]\nonumber\\ \displaystyle & & \displaystyle \times \,e^{-\text{i}(k_{x}x+k_{y}y+k_{z}z)}\nonumber\\ \displaystyle & = & \displaystyle -e^{-\text{i}(k_{x}x+k_{y}y+k_{z}z)}\frac{\partial }{\partial z}\!\left(\frac{1}{n}\frac{\partial n}{\partial x}A_{x}e^{\text{i}(k_{x}x+k_{y}y+k_{z}z)}\!\right)\!\hat{z}\nonumber\\ \displaystyle & = & \displaystyle -\Bigg[\left(\frac{1}{n}\frac{\partial ^{2}}{\partial z\partial x}-\frac{(\partial n/\partial z)(\partial n/\partial x)}{n^{2}}\right)A_{x}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{n}\frac{\partial n}{\partial x}\left(\frac{\partial A_{x}}{\partial z}+\text{i}k_{z}A_{x}\right)\Bigg]\hat{z}\nonumber\\ \displaystyle & {\approx} & \displaystyle -\text{i}k_{z}\frac{1}{n}\frac{\partial n}{\partial x}A_{x}\hat{z}\end{eqnarray}$$

or

(25)$$\begin{eqnarray}\frac{\partial A_{z}}{\partial z}\approx -\frac{1}{n}\frac{\partial n}{\partial x}A_{x}.\end{eqnarray}$$

Finally, recalling that in paraxial approximation change of the angle between the wavevector and $z$-axis, $d{\it\varphi}$, equals the ratio of the $x$ and $z$ components of the wavevector or the ratio of the $z$- and $x$-components of the electric field amplitude, Equation (25) can be rewritten as

(26)$$\begin{eqnarray}d{\it\varphi}=\frac{-{\it\delta}A_{z}}{A_{x}}\approx \frac{1}{n}\frac{\partial n}{\partial x}dz,\end{eqnarray}$$

which is the eikonal equation^{[Reference Born and Wolf15]}. We use integration of this equation (see, for example, Equation (7) in Ref. [Reference Semak and Shneider12]) in our previous works^{[Reference Semak and Shneider11, Reference Semak and Shneider12]} in order to compute radial distribution of the angle of the wavevector after the propagation of a focused laser beam through the caustic under conditions when focusing due to the Kerr effect and defocusing due to ionization take place (see, for example, Equation (8) in Ref. [Reference Semak and Shneider12]).

At the same time, the local projections of the wavevector are uniquely determined from the system of equations

(27)$$\begin{eqnarray}\displaystyle k_{x}^{2}+k_{z}^{2} & = & \displaystyle \frac{{\it\omega}^{2}n^{2}}{c^{2}}\nonumber\\ \displaystyle \frac{k_{x}}{k_{z}} & = & \displaystyle \text{tg}{\it\varphi},\end{eqnarray}$$

with the solution

(28)$$\begin{eqnarray}\displaystyle k_{z}(x,z) & = & \displaystyle \frac{{\it\omega}n(x,z)}{c(1+\text{tg}{\it\varphi}(x,z))^{1/2}},\quad \nonumber\\ \displaystyle k_{x}(x,z) & = & \displaystyle \frac{{\it\omega}n(x,z)\text{tg}{\it\varphi}(x,z)}{c(1+\text{tg}{\it\varphi}(x,z))^{1/2}}.\end{eqnarray}$$

Finally, complete Equation (20) for electromagnetic wave propagation written, without affecting its generality, for linear polarization in the $x$–$z$ plane (or instead its simplified equivalent Equation (21) deduced for paraxial and near field propagation) in combination with Equations (26)–(28) and with the addition of the equation describing inhomogeneity of the index of refraction (either self-induced due to the Kerr effect and ionization or externally induced due to thermal effect, large scale turbulence, etc. or inherent due to spatially variable material properties) represents a closed system of equations that provides a unique solution describing electromagnetic wave propagation in an inhomogeneous medium.

6. Example of computational simulation of paraxial propagation of a focused laser beam in nonlinear media

The results of the simulation of laser beam propagation through air (normal conditions) assuming Gaussian temporal shape of the laser pulse with duration ${\it\tau}=100\;\text{fs}$ and a Gaussian spatial beam profile with the beam radius on $1/e^{2}$ level in the waist, $w_{0}=50\;{\rm\mu}\text{m}$, are shown in Figure 2, that presents the local angle of the wavevector, ${\it\varphi}_{z_{s}}$, computed for propagating from $z_{0}=-2z_{\text{R}}$ to $z_{\text{s}}=2z_{\text{R}}$ (here Rayleigh length is $z_{\text{R}}=19.62\;\text{mm}$) as a function of dimensionless laser radius, $r/w(z_{s})$ (ratio of radial coordinate to the Gaussian beam radius $w(z_{s})$ at the end of the computational range, $z_{s}$), for different moments of dimensionless time, ${\it\eta}/{\it\tau}$ (negative – before the pulse), and for two pulse energies of 0.2 and 0.4 mJ.

Figure 2. The angle of individual rays at the exit of the zone with high beam intensity where induced refraction is large computed for laser beam propagation through air (normal conditions) assuming Gaussian temporal shape of the laser pulse with duration ${\it\tau}=100\;\text{fs}$ and a Gaussian spatial beam profile with the beam radius on $1/e^{2}$ level in the waist, $w_{0}=50\;{\rm\mu}\text{m}$, shown as a function of dimensionless laser radius, $r/w(z_{s})$, for different moments of dimensionless time, ${\it\eta}/{\it\tau}$, (negative – before the pulse): (a) total radiated energy per pulse $E_{\text{pulse}}=0.2~\text{mJ}$ (maximum irradiance $I_{0}=2.87\times 10^{17}\;\text{W}\;\text{m}^{-2}$); (b) total radiated energy per pulse $E_{\text{pulse}}=0.4~\text{mJ}$ (maximum irradiance $I_{0}=5.75\times 10^{17}\;\text{W}\;\text{m}^{-2}$). All computational results correspond to the laser wavelength ${\it\lambda}=800\;\text{nm}$.

The simulation results show that the contribution to the refractive index due to the Kerr effect (Figure 2a) results in focusing of the central part of the beam. The diffraction divergence of the outer area of the beam is partially compensated and in the locations with dimensionless laser radius $r/w(z_{s})$ larger than approximately unity, the self-induced refraction is negligible and the beam divergence is governed predominantly by diffraction. For dimensionless laser radius $r/w(z_{s})<1$ the nonlinear focusing is dominant. In this region the optical power of the Kerr self-focusing ‘lens’ is increasing with time, reaching its maximum when laser power is at the maximum value (${\it\eta}=0$) and then symmetrically (in time) decreasing. Simultaneously with refraction the nonlinear Kerr produces a shift of the laser frequency numerically described in Ref. [Reference Semak and Shneider12].

In nonlinear media that are ionized by the laser radiation the contribution of the Kerr effect that produces focusing is combined with the defocusing effect due to nonuniform ionization (Figure 2b). The ionization is considered as in Ref. [Reference Semak and Shneider12]. We assume that the main ion produced in air is $\text{O}_{2}^{+}$. For the laser wavelength ${\it\lambda}=800\;\text{nm}$ and typical values of the radiation intensity at which self-focusing occurs, the Keldysh parameter ${\it\gamma}>1$ and the limit of multiphoton ionization are valid. The electron number density is given by the equation

(29)$$\begin{eqnarray}N_{e}(r,z,t)=\int _{-\infty }^{t}{\it\sigma}_{k}[N_{0}-N_{e}(r,z,{\it\xi})][I(r,z,{\it\xi})]^{k}d{\it\xi}\end{eqnarray}$$

where $N_{0}$ is the number density of neutrals, and $I(r,z,t)$ is the laser pulse intensity; ${\it\sigma}_{k}$ is the ionization rate due to absorption of $k$ photons such that $k=\text{int}[I_{i}/\hbar{\it\omega}]+1$, where $I_{i}$ is the potential ionization of gas and ${\it\omega}$ is the laser angular frequency.

The effect of media ionization on the beam propagation results in strong divergence of the near-axis part of the beam (Figure 2b). The maximum angle of divergence is increasing in time as the degree of ionization grows, reaching its maximum at the end of the laser pulse. Simultaneously, the outer part of the beam behaves similarly to the above described dynamics, dominated by the focusing due to the Kerr effect.

It was shown in our previous work^{[Reference Semak and Shneider11]} that the critical laser power for self-focusing depends on the laser intensity, which is determined by the effective radius of a Gaussian beam. This power is substantially less than independent of laser intensity critical power computed in accordance with the current theory^{[Reference Chiao, Garmire and Townes16]}. For example, our theory predicts that for the interaction conditions as in Figure 2 the critical power is $P_{\text{cr}}\approx 0.135\;\text{GW}$. Whereas, for the same conditions, the current theoretical model predicts the almost order of magnitude larger value $P_{\text{cr}}\approx 1.7\;\text{GW}$.

7. Concluding remarks

For obvious reasons, it is usual practice in the majority of nonspecialized educational courses and textbooks to ignore the contribution of the term containing ${\rm\nabla}{\it\varepsilon}$ and describe electromagnetic wave propagation using the equation deduced for uniform media. Unfortunately, without much consideration, this propagation equation was adopted in nonlinear optics and after modifications it took the form of Equations (15) and (16). Then, a relatively recent trend mandated further modification of this simplified equation in order to acquire an appearance similar to the nonlinear Schrodinger’s equation (i.e., Equation (17)). Here we would like to reinforce that omitting the third term in the left-hand side of the above derived complete propagation Equation (9) leads to inadequate description of the physics involved. Indeed, it is self-contradictory to disregard the inhomogeneity of optical properties of media (both induced and inherent) at the stage of deduction of the propagation equation and then to reintroduce into the obtained simplified propagation equation the nonlinear dependence of the index of refraction, $n$, on the laser irradiance. Trivial estimates show a non-negligible contribution of the term containing ${\rm\nabla}{\it\varepsilon}$ in complete Equation (10) for laser beam propagation in nonlinear media. Therefore, ignoring this term leads to false predictions. One of the examples of such false prediction is waveguide like propagation of a laser beam in nonlinear media^{[Reference Chiao, Garmire and Townes16]} that leads to the development of fascinating concepts of optical soliton and laser beam filamentation that recently produced a flurry of extensive theoretical research. It is easy to see that for a laser beam with intensity distribution that is near-Gaussian a possible solution of an incorrectly abbreviated propagation equation, such as Equations (15)–(17), indeed has a self-similar form. The currently accepted interpretation of this self-similarity is that the diffraction divergence and divergence produced due to media ionization is compensated by self-focusing^{[Reference Agrawal2–Reference Chin4]}. One of the results, following on from the self-similarity of the solution of this inadequate propagation Equations (15)–(17), is a captivating (however, contradictory to experimental observations) prediction that kilometers long transmission of the laser beam can be achieved in an atmosphere without the beam diverging^{[Reference Couairon and Mysyrowicz6]}.

Our work demonstrates that the solution of complete Equation (10) that adequately describes laser beam propagation in nonlinear media does not have self-similar form. As a demonstration we solved a complete propagation equation for the conditions when input from the nonlinear refraction can be treated as perturbation of the solution of the linear Helmholtz equation describing propagation of a focused laser beam^{[Reference Semak and Shneider11, Reference Semak and Shneider12]}. This solution demonstrated that laser beam divergence is affected by Kerr self-focusing and plasma defocusing differently in different radial locations of the laser beam and in different times during the laser pulse, i.e., self-similar beam propagation does not occur.

Another inconsistency of the customary approach in which Equation (17) or its derivatives are utilized for modeling of diverging or converging laser beam propagation in nonlinear media is that Equation (17) is derived under the assumption of a plane wavefront. Consequently, it is inapplicable for the description of nonplanar wavefronts. However, in a self-contradictory manner similar to the above illustrated reinsertion of the nonlinear effect into Equation (17), the complex nonplanar wave propagation is described in all current theoretical models using plane wave propagation Equation (17) or its derivatives. At this point it is worth mentioning that the Poynting vector maintains its direction in the approximation of a plane wavefront. In contrast, for a converging or diverging laser beam the local direction of the Poynting vector varies as a function of the distance from the axis. Our previous eikonal-paraxial model^{[Reference Semak and Shneider11, Reference Semak and Shneider12]} correctly reflects this behavior of the Poynting vector.

The effect of inhomogeneity of the dielectric constant on electromagnetic wave propagation is known in several relatively narrowly specialized fields. So far, the practical application of this concept was limited to the theory of radio wave propagation in the ionosphere (see Ref. [Reference Budden13]) and electromagnetic and acoustic wave propagation in stratified media^{[Reference Brekhovskikh20]}, such as radar propagation in atmospheric boundary layers^{[Reference Dockery21]}.

In conclusion we summarize the contribution of this work to the field of nonlinear optics as follows. (1) Realization that the gradient of dielectric constant always provides a non-negligible contribution in the propagation equation of a laser beam with realistic beam profile because the characteristic length of change of irradiance is comparable to the ‘beam size’ (for any reasonable definition of this physical property). (2) Development of the method described in Refs. [Reference Semak and Shneider11, Reference Semak and Shneider12] in which we integrated diffractive and geometrical optics by blending the solution of the linear Maxwell’s equation and a correction term that represents nonlinear field perturbation expressed as solution of a paraxial ray-optics (eikonal) equation that opens an elegant means of numerical computation of the ray trajectories (avoiding singularities) as the focused laser beam propagates in a nonlinear and ionized medium through its caustic (the area near the focal plane that extends several Rayleigh lengths).

The realm of nonlinear optics that deals with laser beam propagation has benefited from a multitude of experimental works and significant experience has been accumulated in solving complex mathematical problems. However, it seems that substantial improvement of the theoretical part of nonlinear optics is needed and revision of the fundamentals of the theoretical model provided in our work can revitalize and substantially advance this field.