Nonparaxial Cartesian and azimuthally symmetric waves with concentrated wavevector and frequency spectra

TL;DR

This paper introduces a theoretical framework for analyzing nonparaxial waves with concentrated spectra, enabling the description of complex azimuthally symmetric fields beyond traditional paraxial limits.

Contribution

It extends the paraxial formalism to handle superpositions of nonparaxial waves, especially azimuthally symmetric Bessel beam superpositions with large cone angles.

Findings

Describes nonparaxial Bessel-Gauss beams with enhanced curvature

Demonstrates superpositions produce complex transverse patterns

Provides analytical tools for azimuthally symmetric wave analysis

Abstract

In this paper, we develop a theoretical analysis to efficiently handle superpositions of waves with concentrated wavevector and frequency spectra, allowing an easy analytical description of fields with interesting transverse profiles. First, we analyze an extension of the paraxial formalism that is more suitable for superposing these types of waves, as it does not rely on the use of coordinate rotations combined with paraxial assumptions. Second, and most importantly, we leverage the obtained results to describe azimuthally symmetric waves composed of superpositions of zero-order Bessel beams with close cone angles that can be as large as desired, unlike in the paraxial formalism. Throughout the paper, examples are presented, such as Airy beams with enhanced curvatures, nonparaxial Bessel-Gauss beams and Circular Parabolic-Gaussian beams (which are based on the Cartesian Parabolic-Gaussian beams), and experimental data illustrates interesting transverse patterns achieved by superpositions of beams propagating in different directions.

Equations

$$\Psi (x,y,z,t)={\int }_{-\infty }^{\infty }\text{d}\omega {\int }_{-\infty }^{\infty }\text{d}{k}_x{\int }_{-\infty }^{\infty }\text{d}{k}_y\tilde{S}(k_x,k_y,\omega )e^{-i\omega t}{e}^{i{k}_{x}x}{e}^{i{k}_{y}y}{e}^{i\sqrt{k^2(\omega )-(k_x^2+k_y^2)}z}$$

$$k_z=$$

$$\sqrt{k^2(\omega )-k_{{\perp }_0}^2}-\frac{k_{x_0}}{\sqrt{k^2(\omega )-k_{{\perp }_0}^2}}{k}_x^{\prime }-\frac{k_{y_0}}{\sqrt{k^2(\omega )-k_{{\perp }_0}^2}}{k}_y^{\prime }$$

$$-\frac{k^2(\omega )}{2{\left[k^2(\omega )-k_{{\perp }_0}^2\right]}^{\frac{3}{2}}}{k_x^{\prime }}^2-\frac{k^2(\omega )}{2{\left[k^2(\omega )-k_{{\perp }_0}^2\right]}^{\frac{3}{2}}}{k_y^{\prime }}^2-\frac{k_{x_0}{k}_{y_0}}{{\left[k^2(\omega )-k_{{\perp }_0}^2\right]}^{\frac{3}{2}}}{k}_x^{\prime }{k}_y^{\prime }$$

$$\sqrt{k^2(\omega )-k_{{\perp }_0}^2}\approx \sqrt{k_0^2-k_{{\perp }_0}^2}+\frac{u}{c}\frac{k_0}{\sqrt{k_0^2-k_{{\perp }_0}^2}}$$

$$a_x\stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}\frac{k_{x_0}}{\sqrt{k_0^2-k_{{\perp }_0}^2}}\text{,}{a}_y\stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}\frac{k_{y_0}}{\sqrt{k_0^2-k_{{\perp }_0}^2}}$$

$$b\stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}\frac{1}{2}\frac{k_0^2}{{\left(k_0^2-k_{{\perp }_0}^2\right)}^{\frac{3}{2}}}$$

$$d\stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}\frac{k_{x_0}{k}_{y_0}}{{\left(k_0^2-k_{{\perp }_0}^2\right)}^{\frac{3}{2}}}$$

$$e\stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}\frac{k_0}{\sqrt{k_0^2-k_{{\perp }_0}^2}}$$

$$\Psi (x,y,z,t)\approx e^{-i{\omega }_{0}t}{e}^{i{k}_{x_0}x}{e}^{i{k}_{y_0}y}{e}^{i\sqrt{k_0^2-k_{{\perp }_0}^2}z}A(x,y,z,t)$$

$$A(x,y,z,t)\stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}{\int }_{-\infty }^{+\infty }\text{d}u{\int }_{-\infty }^{+\infty }\text{d}{k}_x^{\prime }{\int }_{-\infty }^{+\infty }\text{d}{k}_y^{\prime }S(k_x^{\prime },k_y^{\prime },u)e^{-iut}{e}^{i{k}_x^{\prime }x}{e}^{i{k}_y^{\prime }y}{e}^{i{k}_z^{\prime }z}$$

$$k_z^{\prime }=-a_x{k}_x^{\prime }-a_y{k}_y^{\prime }-d{k}_x^{\prime }{k}_y^{\prime }-b(k_x^{\prime 2}+k_y^{\prime 2})+\frac{e}{c}u$$

$$\left|\frac{\partial A}{\partial x_i}\right|\ll \left|k_{x_{i_0}}A\right|\text{,}\left|\frac{{\partial }^{2}A}{\partial x_i^2}\right|\ll \left|k_{x_{i_0}}\frac{\partial A}{\partial x_i}\right|$$

$$\left|\frac{\partial A}{\partial t}\right|\ll |{\omega }_{0}A|\text{,}\left|\frac{{\partial }^{2}A}{\partial t^2}\right|\ll \left|{\omega }_0\frac{\partial A}{\partial t}\right|$$

$$\frac{1}{i}\frac{\partial }{\partial x_i}⟷k_{x_i}^{\prime }\text{,}-\frac{{\partial }^2}{\partial x_i^2}⟷k_{x_i}^{\prime 2}\text{,}-\frac{1}{i}\frac{\partial }{\partial t}⟷u$$

$$\frac{\partial A}{\partial z}=-\frac{e}{c}\frac{\partial A}{\partial t}-\vec{a}\cdot {\vec{\nabla }}_{\perp }A+ib{\nabla }_{\perp }^{2}A+id\frac{{\partial }^{2}A}{\partial x\partial y}$$

$$\Psi (x,z,t)\approx e^{-i{\omega }_{0}t}{e}^{i{k}_{x_0}x}{e}^{i\sqrt{k_0^2-k_{{\perp }_0}^2}z}A(x,z,t)$$

$$A(x,z,t)={\int }_{-\infty }^{+\infty }\text{d}u{\int }_{-\infty }^{+\infty }\text{d}{k}_x^{\prime }S(k_x^{\prime },u)e^{-iut}{e}^{i{k}_x^{\prime }x}{e}^{i{k}_z^{\prime }z}$$

$$k_z^{\prime }=-a_x{k}_x^{\prime }-b{k}_x^{\prime 2}+\frac{e}{c}u$$

$$\frac{\partial A}{\partial z}=-\frac{e}{c}\frac{\partial A}{\partial t}-a_x\frac{\partial A}{\partial x}+ib\frac{{\partial }^{2}A}{\partial x^2}$$

$$\gamma \stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}\frac{k_{x_0}}{k_0}\text{,}\alpha \stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}{k}_{x_0}{x}_0$$

$$s\stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}\frac{x}{x_0}\text{,}\xi \stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}\frac{z}{k_0{x}_0^2}\frac{1}{(1-{\gamma }^2{)}^{\frac{3}{2}}}\text{,}\tau \stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}\frac{c}{k_0{x}_0^2}t$$

$$\frac{{\partial }^{2}A}{\partial s^2}+2i\frac{\partial A}{\partial \xi }+2i\alpha (1-{\gamma }^2)\frac{\partial A}{\partial s}+2i(1-{\gamma }^2)\frac{\partial A}{\partial \tau }=0$$

$$\Psi (s,\xi ,\tau )\approx$$

$$\times e^{-i\frac{{\alpha }^2}{{\omega }_0{\gamma }^2}u\tau }{e}^{i{\tilde{k}}_x^{\prime }s}{e}^{-i\left[\alpha (1-{\gamma }^2){\tilde{k}}_x^{\prime }+\frac{1}{2}{\tilde{k}}_x^{\prime 2}-\frac{{\alpha }^2}{{\gamma }^2}(1-{\gamma }^2)\frac{u}{{\omega }_0}\right]\xi }$$

$$\overline{x}$$

$$\overline{y}$$

$$\overline{t}$$

$$\overline{z}$$

$$\frac{\partial A}{\partial \overline{z}}=i\mathcal{L}[A]\text{,}\mathcal{L}\stackrel{\text{makebox}[0.0pt]\text{mbox}def}{=}b{\nabla }_{\perp }^{\prime 2}+d\frac{{\partial }^2}{\partial \overline{x}\partial \overline{y}}$$

$$\zeta =a_1\overline{x}+a_2\overline{y}\text{,}\eta =a_3\overline{x}+a_4\overline{y}$$

$$\zeta$$

$$\eta$$

$$\zeta =\frac{1}{\sqrt{2(4b^2-d^2)}}[2b\overline{x}-(d+\sqrt{4b^2-d^2})\overline{y}]$$

$$\eta =\frac{1}{\sqrt{2(4b^2-d^2)}}[2b\overline{x}-(d-\sqrt{4b^2-d^2})\overline{y}]$$

$$\frac{\partial A}{\partial \overline{z}}=ib\left[\frac{{\partial }^{2}A}{\partial {\zeta }^2}+\frac{{\partial }^{2}A}{\partial {\eta }^2}\right]$$

$$S(k_x^{\prime },k_y^{\prime },u)=\frac{1}{(2\pi {)}^3}{\int }_{-\infty }^{+\infty }\text{d}{x}^{\prime }{\int }_{-\infty }^{+\infty }\text{d}{y}^{\prime }{\int }_{-\infty }^{+\infty }\text{d}{t}^{\prime }A(x^{\prime },y^{\prime },0,t^{\prime })e^{-i{k}_x^{\prime }{x}^{\prime }}{e}^{-i{k}_y^{\prime }{y}^{\prime }}{e}^{iu{t}^{\prime }}$$

$$A(x,y,z,t)=\frac{1}{4\pi ibz\sqrt{1-\frac{d^2}{4b^2}}}{\int }_{-\infty }^{+\infty }\text{d}{x}^{\prime }{\int }_{-\infty }^{+\infty }\text{d}{y}^{\prime }A\left(x^{\prime },y^{\prime },0,t-\frac{e}{c}z\right)$$

Full text

1 Introduction

Appropriate superpositions of waves propagating in different directions may create interesting and useful interference patterns in the transverse plane that keep their shape up to a certain propagation distance, which is the range within which the waves interfere significantly. Consequently, it is higher when the spot sizes of the superposed waves are bigger, which implies that their wavevector spectra are concentrated around their propagation directions. Since these can make large angles with respect to the longitudinal axis, the paraxial wave equation may not describe the propagation of the superposition properly.

In principle, each wave could be analyzed separately using the paraxial formalism in a system of rotated coordinates whose z-axis coincide with its propagating direction and the application of the inverse rotation would express the wave as a function of the desired non-rotated coordinates. However, this procedure is not straightforward and results in complicated expressions, as the rotations mix the transverse and longitudinal non-rotated coordinates.

To allow a more suitable and analytically simple way of handling superpositions of waves with concentrated spectra, we analyze and apply an extension of the paraxial formalism that deals with combinations of plane waves whose wavevectors are close to any desired direction without relying on coordinate rotations. The resulting expressions do not mix transverse and longitudinal coordinates while still incorporate the inclination of the plane waves very clearly through parameters. Moreover, if desired, the method can be straightforwardly extended to include higher-order corrections, so it is not limited to the accuracy of the paraxial approximation.

As the main contribution of this work, we then leverage these results to develop a novel analytical treatment for azimuthally symmetric waves with concentrated spectra, that is, waves composed of superpositions of zero-order Bessel beams with close cone angles. However, in contrast to the paraxial assumption of small angles with respect to the longitudinal axis, this formalism encompasses cone angles as large as desired and provides simple expressions for this class of nonparaxial waves. In particular, similarities with the unidimensional222We will say that a wave is n-dimensional when it is a function of n transverse spatial coordinates, in addition to the longitudinal coordinate $z$ and the time $t$. Cartesian waves here analyzed allow us to build azimuthally symmetric versions of them, which can be used to derive new types of waves. As examples, we provide an analytical expression for nonparaxial Bessel-Gauss beams and for azimuthally symmetric versions of Parabolic-Gaussian beams (which we refer to as Circular Parabolic-Gaussian beams).

Throughout this work, we also present experimental demonstrations and comparisons to numerically-calculated Rayleigh-Sommerfeld diffraction integrals to corroborate the validity and accuracy of the theoretical results.

2 Integral formulation and wave equation

2.1 Bidimensional case

Let $\Psi(x,y,z,t)$ be an arbitrary solution of the homogeneous scalar wave equation $(\nabla^{2}-\partial_{ct}^{2})\Psi=0$. It can be represented by a superposition of plane waves with different wavevectors $\vec{k}=k_{x}\hat{x}+k_{y}\hat{y}+k_{z}\hat{z}$ and angular frequencies $\omega$, given that $|\vec{k}|=k(\omega)=\omega/c$, where $c$ is the speed of light in the medium. If we consider only waves that propagate in the $+\hat{z}$ direction, $k_{z}=\sqrt{k^{2}(\omega)-(k_{x}^{2}+k_{y}^{2})}$ and $\Psi(x,y,z,t)$ can be written as

[TABLE]

where $\tilde{S}(k_{x},k_{y},\omega)$ is the spectrum, which determines the amplitude and phase of each plane wave. The limits of the integrals in Eq. (1) allow imaginary values of $k_{z}$, accounting for the presence of evanescent waves. If only propagating waves are desired, the integration should be performed in the region $k_{x}^{2}+k_{y}^{2}\leq\omega^{2}/c^{2}$.

If the spectrum $\tilde{S}(k_{x},k_{y},\omega)$ is concentrated around $(k_{x},k_{y},\omega)=(k_{x_{0}},k_{y_{0}},\omega_{0})$, the plane waves have wavevectors that vary only slightly around $\vec{k}_{0}=k_{x_{0}}\hat{x}+k_{y_{0}}\hat{y}+\sqrt{k_{0}^{2}-(k_{x_{0}}^{2}+k_{y_{0}}^{2})}\hat{z}$, where $k_{0}\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}k(\omega_{0})=\omega_{0}/c$. Mathematically, this means that $|\tilde{S}(k_{x},k_{y},\omega)|$ is significant only in the region in which $|k_{x}^{\prime}/k_{x_{0}}|\ll 1$, $|k_{y}^{\prime}/k_{y_{0}}|\ll 1$ and $|u/\omega_{0}|\ll 1$, with $k_{x}^{\prime}\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}k_{x}-k_{x_{0}}$, $k_{y}^{\prime}\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}k_{y}-k_{y_{0}}$ and $u\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}\omega-\omega_{0}$ being defined as the deviation variables 333Additionally, we may assume that $k_{x_{0}}<k_{0}$, $k_{y_{0}}<k_{0}$ and $k_{x_{0}}^{2}+k_{y_{0}}^{2}<k_{0}^{2}$, so that all the waves are propagating. In this case, even though we keep the limits of the integrals in $k_{x}$ and $k_{y}$ from $-\infty$ to $\infty$, the concentration of the spectrum implies that there are no significant evanescent waves in the superposition.. In this region, $k_{z}$ can be approximated to second order in $k_{x}^{\prime}$ and $k_{y}^{\prime}$ using a Taylor series expansion444If desired, higher-order corrections can be easily included to improve accuracy, even though the second-order approximation is usually enough.:

[TABLE]

where $k_{{\perp}_{0}}^{2}\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}k_{y_{0}}^{2}+k_{x_{0}}^{2}$. For the variation in $\omega$, we perform a first order Taylor series approximation only in the most significant term, which is the first one on the right-hand side of Eq. (2), resulting in

[TABLE]

In all the other terms, we take $k(\omega)\approx k_{0}$. Assigning names to the constants of Eq. (1) and Eq. (2):

[TABLE]

and performing the integrals of Eq. (1) on the deviation variables, we get

[TABLE]

[TABLE]

where $S(k^{\prime}_{x},k^{\prime}_{y},u)\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}\tilde{S}(k^{\prime}_{x}+k_{x_{0}},k^{\prime}_{y}+k_{y_{0}},u+\omega_{0})$ and

[TABLE]

is the deviation of $k_{z}$ with respect to $k_{z_{0}}\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}\sqrt{k_{0}^{2}-k_{{\perp_{0}}}^{2}}$.

Eq. (5) shows that $\Psi(x,y,z,t)$ is approximately a plane wave with wavevector $\vec{k}_{0}$ and frequency $\omega_{0}$ modulated by an envelope $A(x,y,z,t)$ which, according to Eq. (6), is composed of a superposition of plane waves with wavevectors $\vec{k}^{\prime}=k^{\prime}_{x}\hat{x}+k^{\prime}_{y}\hat{y}+k^{\prime}_{z}\hat{z}$ and angular frequencies $u$. Since $|k^{\prime}_{x}|$, $|k^{\prime}_{y}|$, $|k^{\prime}_{z}|$ and $|u|$ are small compared to $|k_{x_{0}}|$, $|k_{y_{0}}|$, $|k_{z_{0}}|$ and $\omega_{0}$, respectively, $A(x,y,z,t)$ has slow variations in $x$, $y$, $z$ and $t$. Mathematically, this means

[TABLE]

where $x_{i}$ stands for any spatial coordinate.

The differential equation satisfied by $A(x,y,z,t)$ is uniquely defined by the relation among the frequency and wavevector components of its composing plane waves, given by Eq. (7). Each differentiation of $A(x,y,z,t)$ with respect to a variable implies a multiplication of its spectrum by a factor proportional to the corresponding reciprocal variable, resulting in the correspondences:

[TABLE]

where $x_{i}$ stands for any spatial coordinate. Therefore, Eq. (7) is equivalent to

[TABLE]

where $\vec{a}\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}a_{x}\hat{x}+a_{y}\hat{y}$ and $\vec{\nabla}_{\perp}\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}\hat{x}\partial_{x}+\hat{y}\partial_{y}$. The coefficient of the mixed derivative in Eq. (10) implies that the equation has separable solutions of the form $\Psi(x,y,z,t)=f(x,z)g(y,z)h(t)$ only when $d=0$, that is, when at least one of $k_{x_{0}}$ and $k_{y_{0}}$ is zero. If they are both zero, $\vec{a}=\vec{0}$, $b=(2k_{0})^{-1}$, $d=0$ and $e=1$ and Eq. (10) reduces to the well-known paraxial wave equation $\frac{\partial A}{\partial z}=-\frac{1}{c}\frac{\partial A}{\partial t}+\frac{i}{2k_{0}}\nabla_{\perp}^{2}A$, since in this case the superposed plane waves have wavevectors almost parallel to the z-direction.

2.2 Unidimensional case

If $\Psi$ is a function of only one of the transverse coordinates, say $\Psi(x,z,t)$, the same kind of analysis provides results similar to those of sec. 2.1:

[TABLE]

with $k_{{\perp}_{0}}=k_{x_{0}}$. In this case, it is possible to write Eq. (14) in a dimensionless form by defining new constants and dimensionless variables:

[TABLE]

where $x_{0}$ is a constant. Then, the differential equation for $A(s,\xi,\tau)$ becomes

[TABLE]

When $k_{x_{0}}=0$ it correctly reduces to the dimensionless unidimensional paraxial wave equation $\frac{\partial^{2}A}{\partial s^{2}}+2i\frac{\partial A}{\partial\xi}+2i\frac{\partial A}{\partial\tau}=0$.

Using these constants and dimensionless variables and defining the dimensionless wavevector $\tilde{k}_{x}^{\prime}\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}k_{x}^{\prime}x_{0}$, the wave of Eq. (11) can be written as

[TABLE]

where $\mathcal{S}(\tilde{k}^{\prime}_{x},u)\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}S(\tilde{k}^{\prime}_{x}/x_{0},u)/x_{0}$ is a concentrated spectrum around $(\tilde{k}^{\prime}_{x},u)=(0,0)$.

2.3 Reduction to the paraxial wave equation

By means of some changes of variables, it is possible to reduce the general Eq. (10) to the paraxial wave equation, so we can leverage all we know about this equation to Eq. (10). The first order derivatives in $x$, $y$ and $t$ can be eliminated if we introduce the shifted variables $\overline{x}$, $\overline{y}$ and $\overline{t}$ according to

[TABLE]

Applying the chain rule and defining $\vec{\nabla}^{\prime}_{\perp}\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}\hat{x}{\partial}_{\overline{x}}+\hat{y}{\partial}_{\overline{y}}$, Eq. (10) can be rewritten as

[TABLE]

Note that the shifts in the transverse coordinates express the inclination of $\vec{k}_{0}$ with respect to the z-direction and the shift in the temporal variable represents the wave’s travel time to arrive at the longitudinal position $z$. Eq. (20) shows that any solution of Eq. (10) must be of the form $A(x-a_{x}z,y-a_{y}z,z,t-ez/c)$.

With an additional change of variables of the type

[TABLE]

it is possible to eliminate the mixed derivative of the operator $\mathcal{L}$. Two alternatives that reduce it to $\mathcal{L}=b\left[\frac{\partial^{2}}{\partial\zeta^{2}}+\frac{\partial^{2}}{\partial\eta^{2}}\right]$ are

[TABLE]

and

[TABLE]

so that Eq. (20) becomes

[TABLE]

which is like a paraxial wave equation in the variables $\zeta$ and $\eta$. The factor $b$ was not changed to $(2k_{0})^{-1}$ because, as already pointed out in sec. 2.1, $b=(2k_{0})^{-1}$ if $k_{x_{0}}=k_{y_{0}}=0$. However, if desired, an additional scaling in both $\zeta$ and $\eta$ can reduce Eq. (24) precisely to the paraxial wave equation in the new scaled variables, which are $\overline{\zeta}=\zeta/\sqrt{2bk_{0}}$ and $\overline{\eta}=\eta/\sqrt{2bk_{0}}$. In the unidimensional case, Eq. (14) may be reduced to $\partial_{\overline{z}}A=ib\partial_{\overline{x}\overline{x}}A$ and $A$ also obeys the paraxial wave equation in the scaled variable $\chi=\overline{x}/\sqrt{2bk_{0}}$.

Therefore, any solution of the paraxial wave equation can be extended to the case being considered, which is natural, since this approach should be analogous to a rotation of coordinates. Note, however, that the former is simpler: while the rotation expresses the $z$ coordinate as a combination of both transverse and longitudinal rotated coordinates, the presented method keep the transverse and longitudinal coordinates separated.

The choice of Eq. (22) for $\zeta$ and $\eta$, although asymmetrical in $\overline{x}$ and $\overline{y}$, reduces to $\zeta=\overline{x}$ and $\eta=\overline{y}$ when the “coupling factor” $d$ is zero, showing that the lack of the mixed derivative makes it possible to find a separable solution for $A$ in the variables $\overline{x}$ and $\overline{y}$. The choice of Eq. (23), on the other hand, has a better symmetry in $\overline{x}$ and $\overline{y}$, but reduces to $\zeta=(\overline{x}-\overline{y})/\sqrt{2}$ and $\eta=(\overline{x}+\overline{y})/\sqrt{2}$ when $d=0$.

Also, note that if a beam with envelope $A(x,y,z)$ satisfies Eq. (10), the pulse with envelope $\overline{A}(x,y,z,t)=A(x,y,z)B(t-ez/c)$ is also a solution of that equation for an arbitrary differentiable function $B(t-ez/c)$.

3 Diffraction integral

For waves with concentrated spectrum around $(k_{x},k_{y},\omega)=(k_{x_{0}},k_{y_{0}},\omega_{0})$ it is possible to derive a generalization of the Fresnel diffraction integral.

3.1 Bidimensional case

From Eq. (6), it is evident that $A(x,y,0,t)$ is a Fourier transform of $S(k^{\prime}_{x},k^{\prime}_{y},u)$. Therefore, by inverting the transform, we have

[TABLE]

Substituting Eq. (25) in Eq. (6), the integrals in $k^{\prime}_{x}$, $k^{\prime}_{y}$ and $u$ can be solved by completing squares and using the known results $\int_{-\infty}^{+\infty}e^{-ix^{2}}\text{d}x=(1-i)\sqrt{\frac{\pi}{2}}$ and $\int_{-\infty}^{+\infty}e^{iux}\text{d}u=2\pi\delta(x)$, where $\delta(\cdot)$ is the Dirac delta function. This leads to the diffraction integral:

[TABLE]

For beams, that is, when $A=A(x,y,z)$ is not a function of time, Eq. (26) is still valid if we remove the $t$ dependencies. As it should be, for $k_{x_{0}}=k_{y_{0}}=0$ it reduces to the Fresnel diffraction integral. An alternative procedure to derive Eq. (26) is to write the Fresnel diffraction integral for the transformed variables $\overline{\zeta}$ and $\overline{\eta}$ and arrange the terms in $x-a_{x}z$ and $y-a_{y}z$.

3.2 Unidimensional case

An analogous procedure starting with Eq. (12) results in

[TABLE]

4 Helmholtz-Gauss beams

Since $A(\zeta,\eta,z)$ satisfies Eq. (24), which is exactly a paraxial wave equation with $(2k_{0})^{-1}$ replaced by $b$, there exists a corresponding Helmholtz-Gauss beam with envelope $A_{HG}(\zeta,\eta,z)$ given by [7]

[TABLE]

where $\mu=1+i4q^{2}bz$ 555Note that $\mu=1+i2q^{2}z/k_{0}$ for the paraxial wave equation and we have to replace $(2k_{0})^{-1}$ with $b$ to be consistent with Eq. (24)..

Using the transformation of Eq. (22), the envelope is expressed in the original Cartesian coordinates by

[TABLE]

Note that, as $\zeta$ and $\eta$ are linear combinations of $x$, $y$ and $z$, $A\left(\frac{\zeta}{\mu},\frac{\eta}{\mu},\frac{z}{\mu}\right)$ can be expressed as $A\left(\frac{x}{\mu},\frac{y}{\mu},\frac{z}{\mu}\right)$. At $z=0$, we have

[TABLE]

which is a Gaussian apodization with an extra mixed-product term proportional to $d$.

5 Nonparaxial linearly polarized electromagnetic field with concentrated spectra

Since any Cartesian component of an electric field $\vec{E}$ obeys the scalar wave equation, a scalar wave with concentrated spectra can be assigned to one of them and the other field components can be obtained via Maxwell’s equations.

If we choose a linearly polarized electric field in the x-direction $\vec{E}=E_{x}\hat{x}+E_{z}\hat{z}$ with

[TABLE]

the $E_{z}$ component is given by Gauss’s law ($\vec{\nabla}\cdot\vec{E}=0$), which implies

[TABLE]

Integrating by parts when $A(z,y,z,t)$ is in the integrand and using its slow variation (Eq. (8a)) to write $\int(\partial A/\partial x_{i})e^{ik_{z_{0}}z}\,\text{d}z\approx(\partial A/\partial x_{i})\int e^{ik_{z_{0}}z}\,\text{d}z$ (where $x_{i}$ can be any spatial coordinate), we find that

[TABLE]

From Faraday’s law ($\vec{\nabla}\times\vec{E}=-\partial\vec{B}/{\partial t}$), $\vec{B}$ can be calculated by integration in $t$. Once more, we will use integration by parts along with the slow variation of $A(x,y,z,t)$ (Eq. (8b)) to write $\int(\partial A/\partial t)e^{-i\omega_{0}t}\,\text{d}t\approx(\partial A/\partial t)\int e^{-i\omega_{0}t}\,\text{d}t$. Due to Eq. (8b), it is also possible to neglect the second-order derivatives. Therefore, the results are

[TABLE]

6 Examples of Cartesian solutions

6.1 Gaussian pulse with Gaussian bidimensional spatial profile

If we choose a separable Gaussian spectrum in $k^{\prime}_{x}$, $k^{\prime}_{y}$ and $u$:

[TABLE]

where $C$ is an arbitrary constant and $r_{x}$, $r_{y}$ and $r_{t}$ are related to the spot sizes in $x$, $y$ and $t$, respectively, and use the integral result $\int_{-\infty}^{+\infty}\exp(-p^{2}x^{2}\pm qx)dx=\frac{\sqrt{\pi}}{p}\exp\left(\frac{q^{2}}{4p^{2}}\right)$, valid for $\real(p^{2})>0$ [6], Eq. (5) and Eq. (6) give

[TABLE]

6.2 Unidimensional Gaussian beam (GB)

For a single frequency ($\omega_{0}$) and a Gaussian spectrum in $k^{\prime}_{x}$:

[TABLE]

the resultant wave is a unidimensional Gaussian beam (GB):

[TABLE]

6.3 Nonparaxial Cartesian Parabolic-Gaussian (PG) beams

The solutions of the paraxial wave equation known as Cartesian Parabolic-Gaussian (PG) beams were introduced by Bandres et al. in [3] and further analyzed in [8]. Their envelope ${}_{n}A_{\nu}(x,z)$ is described by

[TABLE]

with

[TABLE]

where $z_{R}=kr_{0}^{2}/2$ is the Rayleigh distance, $r_{0}$ is the waist of the Gaussian envelope, $n$ is $1/2$ for even beams and $3/2$ for odd beams, $\nu\in\mathbb{R}$ is the order of the parabolic cylinder function ${{}_{n}P_{\nu}}(x)$, ${{}_{1}F_{1}}(a;b;x)$ is the confluent hypergeometric function and $h\in\mathbb{R}$ is a real parameter. At $z=0$, it becomes

[TABLE]

To obtain the nonparaxial version of this envelope, is is only necessary to replace $x$ by $(x-a_{x}z)/\sqrt{2bk_{0}}$, as explained in sec. 2.3.

Although the previous nonparaxial examples are just inclined versions of the usual paraxial beams, they are useful for describing superpositions of beams (as commented in sec. 7) and for generating new expressions for azimuthally symmetric waves (as shown in sec. 9 and exemplified in sec. 10). However, it is also possible to use the presented formalism to obtain waves with enhanced characteristics with respect to their paraxial counterparts, as illustrated in the next example.

6.4 Nonparaxial Airy Beam

The spectrum

[TABLE]

with free real parameters $\beta$, $\vartheta$ and $\sigma$ results in a nonparaxial Airy beam. $\sigma$ causes a Gaussian apodization of the spectrum, which is manifested in space as an exponential apodization of the Airy spatial profile at $z=0$. As in the paraxial case, the spectrum $\mathcal{S}(\tilde{k}^{\prime}_{x},u)$ is not concentrated unless the beam is apodized.

Using the integral result $\int_{-\infty}^{+\infty}\exp\left[i\left(\frac{t^{3}}{3}+vt^{2}+wt\right)\right]dt=2\pi\exp\left[i(\frac{2}{3}v^{3}-vw)\right]A_{i}(w-v^{2})$ [11], the expression of the beam becomes

[TABLE]

Denoting $\delta\approx 1.01879$, the parabolic trajectory of the beam’s position of maximum amplitude is

[TABLE]

since $x=\delta$ is where $Ai(x)$ has its maximum value. Therefore, $\beta$ controls the linear coefficient and $\vartheta$, the independent coefficient. If we choose $\sigma=0$, $\beta=\alpha(1-\gamma^{2})$ and $\vartheta=-\delta+\beta^{2}$, the trajectory equation becomes simply

[TABLE]

or, as a function of $x$ and $z$,

[TABLE]

While $x_{0}$ controls the spot size of the beam, $k_{x_{0}}$ can be used to choose a desired curvature for its trajectory. To illustrate this, Fig. 1 shows the variation of the beam’s curvature when $k_{x_{0}}$ is varied for $x_{0}=50\,\mu\text{m}$.

7 Symmetric superposition of Cartesian waves

As explained in the introduction, the basic goal of the analysis so far is to provide a simple analytic way to describe the propagation of a superposition of waves traveling in different directions, thus composing interesting interference patterns in the transverse plane. Although these waves obey different approximate equations, since their dissimilar $k_{x_{0}}$ and $k_{y_{0}}$ give rise to different constants in Eq. (10), they are all approximate solutions of the homogeneous scalar wave equation, which implies that their expressions can be legitimately combined.

Let’s take as an example the superposition of two equal unidimensional beams with the same initial envelope $A(x,z=0)$ traveling in opposite directions so that the resulting wavevector spectrum has peaks concentrated both at $+k_{x_{0}}$ and $-k_{x_{0}}$ 666Note that, according to sec. 2.3, $A$ is a function of $\chi$ and $z$ and, at $z=0$, $\chi=x/\sqrt{2bk_{0}}$. Since $b$ is the same for $+k_{x_{0}}$ and $-k_{x_{0}}$, the initial envelope is equal for the two beams.. Applying the diffraction integral (Eq. (27)) for each beam and combining the results, we have

[TABLE]

where $a_{x}$ is related to $+k_{x_{0}}$. After expanding the squares and rearranging the terms, the result is

[TABLE]

At the initial plane $z=0$, the simple superposition of the two patterns gives

[TABLE]

showing that the slowly-varying envelope is modulated by the fast-oscillating function $\cos(k_{x_{0}}x)$. Therefore, the spot radius $\Delta x_{0}$ of the resulting beam is determined by a quarter of this cosine period, that is, $\Delta x_{0}=\pi/(2k_{x_{0}})$.

In any superposition, the resulting field depth $Z$ can be estimated by the distance within which the waves interfere significantly. In this simple case, we can apply a reasoning analogous to that of the upcoming sec. 9.4 to obtain $Z\approx w_{0}k_{0}/(|k_{x_{0}}|\sqrt{2bk_{0}})$, where $w_{0}$ is a characteristic width of the envelope $A(x,z=0)$.

8 Experimental results

To demonstrate that the proposed theoretical expressions can be used to easily describe a superposition of beams with concentrated spectra propagating in different directions, we experimentally generated two superpositions of Gaussian beams that result in illustrative interference patterns. The experimental setup is shown in Fig. 2 and used a reflective Spatial Light Modulator (SLM) to generate the desired patterns.

A He-Ne laser ($\lambda\approx 632.8\,\text{nm}$) is expanded and collimated by a system of lenses ($L_{1}$ and $L_{2}$) and then sent to the SLM, which is a LC-R1080 model from Holoeye Photonics and possesses a display matrix of 1920x1200 with $8.1\,\mu\text{m}$ pixel size. The device was used in amplitude modulation mode, with the polarizer and analyzer angles ($0^{\circ}$ and $90^{\circ}$, respectively) measured with respect to the SLM axis. The reflected beam is then subject to a 4f spatial filtering system. The SLM is positioned at the input plane (at a focal distance $f$ behind the lens $L_{3}$) and the spatial filtering mask (SF, a band-pass circular pupil) is placed at the Fourier plane to select the first diffraction order containing the information of the desired complex field pattern. After another Fourier transform, the field is obtained at the Fourier plane of lens $L_{4}$ and its intensity profile for different propagation distances $z$ is measured by a CCD camera (DMK 41BU02.H from The Imaging Source, with a 1280x960 display and $4.65\,\mu\text{m}$ pixel size).

Due to the size of the SLM’s pixel, it cannot resolve well the fast amplitude and phase variations of beams with $k_{x_{0}}/k_{0}$ and/or $k_{y_{0}}/k_{0}$ ratios that are not small, so we had to choose close-to-zero values for them. Even though these beams are not significantly nonparaxial, the results illustrate how interesting interference patterns can be produced. With suitable equipment, it is straightforward to generate the same kinds of beams in a nonparaxial regime.

8.1 Superposition of two Gaussian beams

The first example consists in a superposition of two Gaussian beams with same amplitudes and spot sizes $r_{x}=r_{y}=300/\sqrt{2}\,\mu\text{m}$, but with different values of $k_{x_{0}}$ and $k_{y_{0}}$. One of them has $(k_{x_{0}},k_{y_{0}})=(+0.003k_{0},0k_{0})$ and the other has $(k_{x_{0}},k_{y_{0}})=(-0.003k_{0},0k_{0})$, so that they propagate in opposite directions in the $x$-axis but do not shift in the $y$-direction. Fig. 3 shows what we should expect for a $y=0$ cut of the beam’s intensity profile: a sinusoidal pattern enveloped by a Gaussian function. Therefore, it can be viewed as a kind of “cartesian analog” of a Bessel-Gauss (BG) beam, which consists in a superposition of GBs over a cone and has a Bessel function profile enveloped by a Gaussian function.

The propagation distance $Z$ for which the sinusoidal pattern is kept can be estimated by the expression in sec. 7, resulting in $Z\approx r_{x}k_{0}/(|k_{x_{0}}|\sqrt{2bk_{0}})\approx 70.71\,\text{mm}$, which agrees with Fig. 3. Fig. 4 shows the good agreement between the theoretical and the measured intensity patterns for $y=0$ cuts at different propagation distances and also presents the measured patterns for each of these positions.

8.2 Superposition of four Gaussian beams

The second example consists in a superposition of four Gaussian beams with same amplitudes and spot sizes $r_{x}=r_{y}=300/\sqrt{2}\,\mu\text{m}$, but with different values of $k_{x_{0}}$ and $k_{y_{0}}$. Two of them are the same as the ones in sec. 8.1 and the other two have $(k_{x_{0}},k_{y_{0}})=(0k,+0.003k)$ and $(k_{x_{0}},k_{y_{0}})=(0k,-0.003k)$, so that, in addition to the two beams being shifted in the x-direction, we now have two more beams that are shifted along the y-axis in opposite directions. Therefore, the interference pattern will be a little more sophisticated, with more maxima and minima than the one in sec. 8.1.

Fig. 5 shows the good agreement between the theoretical and measured intensity patterns for $y=0$ cuts at different propagation distances and also presents the measured patterns for each of these positions.

9 Nonparaxial azimuthally symmetric waves with concentrated spectra

Following a procedure similar to the one in sec. 2, we now present how the previous unidimensional results can be leveraged to describe nonparaxial azimuthally symmetric waves with concentrated spectra.

9.1 Superposition of Bessel beams with close cone angles

Any azimuthally symmetric wave $\Psi(\rho,z,t)$ can be expressed as a superposition of Bessel beams of order zero, that is

[TABLE]

with the spectrum $\tilde{S}(k_{\rho},\omega)$ defining the amplitudes of each beam. Note that the upper limit of $k_{\rho}$ is not limited to $\omega/c$, so that evanescent waves are allowed.

For an arbitrary spectrum, the integrals in Eq. (53) generally do not have a closed analytical representation, but when $\tilde{S}(k_{\rho},\omega)$ is concentrated around $(k_{\rho},\omega)=(k_{\rho_{0}},\omega_{0})$ some simplifications can be applied using an approximate expression for $J_{0}(\cdot)$. As in sec. 2.1, a concentrated spectrum means that $|\tilde{S}(k_{\rho},\omega)|$ is significant only inside the region $|k_{\rho}^{\prime}/k_{\rho_{0}}|\ll 1$ and $|u/\omega_{0}|\ll 1$, where $k_{\rho}^{\prime}\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}k_{\rho}-k_{\rho_{0}}$ and $u\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}\omega-\omega_{0}$. In other words, all the superposed Bessel beams have almost the same cone angle. Note, however, that it can be as large as desired, unlike in the paraxial approximation, which restricts the angle to small values.

Similarly to what was done in sec. 2, the hypothesis of concentrated spectrum allows us to accurately approximate $k_{z}=\sqrt{k^{2}-k_{\rho}^{2}}$ to second order in $k_{\rho}^{\prime}$ and to first order in $u$:

[TABLE]

In this expression, the constants are defined as:

[TABLE]

Changing the integration variable in Eq. (53) from $k_{\rho}$ to $k_{\rho}^{\prime}$, the result is

[TABLE]

with $S(k_{\rho}^{\prime},u)\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny def}}}}{{=}}}\tilde{S}(k_{\rho}^{\prime}+k_{\rho_{0}},u+\omega_{0})$ and

[TABLE]

For $k_{\rho_{0}}\neq 0$, the spectrum is concentrated around $k_{\rho}^{\prime}=0$ and its amplitude is negligible in the region $\absolutevalue{k_{\rho}^{\prime}}\geq k_{\rho_{0}}$. Therefore, the lower limit of the integral in $k_{\rho}^{\prime}$ in Eq. (56) can be extended to $-\infty$ without loss of accuracy.

9.2 An approximate expression of $J_{0}(x)$

Now, we will make use of an unusual, although very accurate, approximate expression for $J_{0}(x)$. It was first introduced in [9], but the reasoning behind its derivation was not presented. Therefore, we show here how this approximation is obtained.

The asymptotic approximation for $J_{0}(x)$ is [2]

[TABLE]

Although it is rigorously valid only when the argument is large, Fig. 6 shows that this asymptotic approximation is precise even for relatively small values of $x$, failing only due to the divergence when $x\rightarrow 0$. Therefore, one could attempt to improve the accuracy of Eq. (58) for small arguments by introducing a function in the denominator of the square root to eliminate the divergence, that is,

[TABLE]

where $g(x)$ is chosen in such a way that $F(x)\approx J_{0}(x)$ for all $x$. To this end, $g(x)$ must have the following properties:

• •

$g(0)=1$, so that $F(0)=J_{0}(0)=1$.

• •

$g^{\prime}(0)=-(\pi-2)$, so that $F^{\prime}(0)=J_{0}^{\prime}(0)=0$.

• •

$g(x)\rightarrow 0$ sufficiently fast when $x$ increases, so that $F(x)\approx J_{0}(x)$ for $x\neq 0$.

The simplest function that satisfies all these conditions is $g(x)=e^{-(\pi-2)x}$ and we can hence accurately approximate $J_{0}(x)$ with the expression

[TABLE]

for any $x\geq 0$, even small. To allow negative arguments, $x$ in the right-hand side of Eq. (60) should be replaced by $|x|$, so that the expression is even, as is $J_{0}(x)$. Fig. 6 displays a comparison between the exact $J_{0}(x)$, the asymptotic approximation and the unusual approximation, showing that the latter represents $J_{0}(x)$ very accurately.

It is worth noting that the approximation in Eq. (60) for $J_{0}(x)$ separates the slowly-decaying envelope $\sqrt{\frac{2}{\pi x+e^{-(\pi-2)x}}}$ from the fast oscillations expressed by $\cos\left(x-\frac{\pi}{4}\right)$.

9.3 Similarities with the unidimensional Cartesian case

Extending the lower limit of the integral in Eq. (56) to $-\infty$ and using the approximation for the Bessel function, the expression becomes

[TABLE]

Since by hypothesis $|k_{\rho}^{\prime}|\ll k_{\rho_{0}}$ where $|S(k_{\rho}^{\prime},u)|$ is significant, the argument $(k_{\rho}^{\prime}+k_{\rho_{0}})\rho\approx k_{\rho_{0}}\rho$ of the Bessel function is always positive in this region, so there is no need to include the absolute value for the approximate expression to be valid. The effect of the small variation $k_{\rho}^{\prime}$ is significant in the fast oscillations of the cosine, but can be neglected in the envelope, which remains always approximately equal to $\sqrt{\frac{2}{\pi k_{\rho_{0}}\rho+e^{-(\pi-2)k_{\rho_{0}}\rho}}}$ and can be taken outside the integral. Expanding the cosine into complex exponentials, we finally get

[TABLE]

where

[TABLE]

has the same functional form of a unidimensional nonparaxial Cartesian envelope (Eq. (12)) with $x$ replaced by $\rho$. The differential equation satisfied by $A(\rho,z,t)$ is then

[TABLE]

and the unidimensional transformations of sec. 2.3 are applicable if $x$ is replaced by $\rho$.

Therefore, the azimuthally symmetric version of a nonparaxial Cartesian wave can be obtained by simply applying the known solution $A(x,z,t)$ to Eq. (62) with $x$ replaced by $\rho$ and $k_{x_{0}}$, by $k_{\rho_{0}}$. The minus sign in $A(-\rho,z,t)$ can be viewed alternativelly as a change in the sign of $a$, implying that in Eq. (62) there are waves traveling in opposite directions.

Over the axis ($\rho=0$), the wave is

[TABLE]

that is, $|\Psi(0,z,t)|=|A(0,z,t)|$, thus allowing the easy prediction of the on-axis intensity. On the other hand, at $z=0$, the result is

[TABLE]

If $A(\rho,0,t)$ has a definite parity in $\rho$, the numerator of Eq. (66) becomes $2\cos(k_{\rho_{0}}\rho-\pi/4)A(\rho,0,t)$ if $A$ is even and $2i\sin(k_{\rho_{0}}\rho-\pi/4)A(\rho,0,t)$ if $A$ is odd, showing that the slowly-varying envelope is modulated by a fast-oscillating sine or cosine factor with period $2\pi/k_{\rho_{0}}$, as in sec. 7. In these cases, we may expect to have waves with spot radii of the order of $\Delta\rho_{0}\sim\pi/(2k_{\rho_{0}})$.

9.4 Approximate field depth

Since the wave in Eq. (62) can be qualitatively viewed as a superposition of waves with envelope $A(\rho,z,t)$ over a cone of half-angle $\theta=\arcsin\left(k_{\rho_{0}}/k\right)$, we can estimate its field depth $Z$ as the distance within which their interference is significant. From [5], we know that this interpretation, besides being qualitatively appealing, is also quantitatively accurate for determining the field depth of paraxial waves, so we can use the unidimensional transformations of sec. 2.3 to generalize this result to the nonparaxial case considered here.

Let $w_{0}$ be a characteristic width of $A(\rho,z,t)$ at $z=0$ which embraces its main features, such as its spot radius, and let’s assume its change is negligible when propagating the distance $Z$. This is reasonable because $A(\rho,z,t)$ is a slowly-varying envelope by hypothesis. In the paraxial case, $Z$ can be estimated by $Z\approx w_{0}/\sin\theta=w_{0}k_{0}/k_{\rho_{0}}$ [5]. Since the expression for the nonparaxial envelope $A(\rho,z,t)$ is obtained from an initial paraxial expression by replacing $\rho$ by $(\rho-a_{x}z)/\sqrt{2bk_{0}}$ (see sec. 2.3), if the characteristic width $w_{0}$ of $A(\rho,z,t)$ is converted back to a paraxial counterpart of this envelope, the result is a characteristic width of $w_{0}/\sqrt{2bk_{0}}$. Therefore, combining this converted length with the paraxial approximation for $Z$, we get

[TABLE]

which is an accurate estimate for the field depth of any azimuthally symmetric wave with concentrated spectra, as illustrated in the examples of sec. 10.

10 Examples of nonparaxial azimuthally symmetric beams

10.1 Bessel-Gauss (BG) beam

A nonparaxial Bessel-Gauss (BG) beam is obtained by choosing a Gaussian spectrum concentrated around $k_{\rho_{0}}$, as the one in Eq. (37) with $k_{x}^{\prime}$ replaced by $k_{\rho}^{\prime}$ and $r_{x}$, by $r_{0}$. Using the Cartesian expression of Eq. (38) in Eq. (62), the result is

[TABLE]

At $z=0$, the beam’s pattern is a Bessel function $J_{0}(k_{\rho_{0}}\rho)$ (in the form of Eq. (60)) with spot size $\Delta\rho_{0}\approx 2.4/k_{\rho_{0}}$ apodized by a Gaussian function of waist $r_{0}$.

For practical situations, it is desired to have $r_{0}\gg\Delta\rho_{0}$, so that a large number of the Bessel’s rings are unattenuated and the BG beam behaves like a Bessel beam for longer distances, resulting in a higher field depth, which can be estimated by making $w_{0}=r_{0}$ in Eq. (67). Note that a large $r_{0}$ implies a concentrated spectrum, which is consistent with our assumptions.

For comparison, the usual expression for a BG beam, valid for paraxial regime, is [5]

[TABLE]

with $Q=1/r_{0}^{2}-ik_{0}/(2z)$.

Fig. 7 contrasts Eq. (69) and Eq. (68) for two situations: paraxial ($k_{\rho_{0}}=0.001k_{0}$) and highly nonparaxial ($k_{\rho_{0}}=0.8k_{0}$). In both cases, the medium is air, $\lambda=632.8\,\text{nm}$ and $r_{0}=1500\,\mu\text{m}$. The first line of Fig. 7 depicts the Gaussian spectra, showing that they are very concentrated. The following two lines show the intensities $|\Psi(\rho,z,t)|^{2}$ predicted by the paraxial expression (Eq. (69), second line) and the generalized nonparaxial expression (Eq. (68), third line) 777Note that the intensities were normalized by their peak at $z=0$.. It is clear that the two expressions agree in the paraxial regime, but in the nonparaxial case Eq. (69) overestimates the field depth predicted by Eq. (68). The last line shows the intensities at $\rho=0$ predicted by both expressions and compares them with the exact intensity calculated using the Rayleigh-Sommerfeld diffraction integral for the initial pattern $\Psi(\rho,0,t)=e^{-i\omega_{0}t}J_{0}(k_{\rho_{0}}\rho)\exp(-\rho^{2}/r_{0}^{2})$. The three results coincide for the paraxial regime, but only Eq. (68) agrees with the diffraction integral for the nonparaxial situation. According to Eq. (67), the estimated field depths are $Z\approx 1.5\,\text{m}$ and $Z\approx 0.87\,\text{mm}$ for the paraxial and nonparaxial examples, respectively, which accurately match the diffraction integral results.

It is worth mentioning that Eq. (68), while valid for the cases of practical interest in which $r_{0}\gg\Delta\rho_{0}$, is much simpler than other nonparaxial extensions of the BG beams found in the literature, which rely on nonparaxial correction terms [4], series of functions [1] or integral representations [10].

10.2 Circular Parabolic-Gaussian (CPG) beam

If we apply the nonparaxial version of the envelope in Eq. (39) to Eq. (62), we get an azimuthally symmetric version of the PG beams, which we will refer to as Circular Parabolic-Gaussian (CPG) beams.

To illustrate, we will consider two examples. The first is an odd beam ($n=3/2$) with $\nu=h=1$ and $r_{0}=200\,\mu\text{m}$ while the second is an even beam ($n=1/2$) with $\nu=-3/2$, $h=0.5$ and $r_{0}=200\,\mu\text{m}$. Both are analyzed in two regimes: paraxial ($k_{\rho_{0}}=0.01k_{0}$) and nonparaxial ($k_{\rho_{0}}=0.5k_{0}$), with $\lambda=632.8\,\text{nm}$.

Since the envelopes of PG beams have definite parity in $\rho$ at $z=0$, the initial patterns of CPG beams are Gaussian-apodized parabolic cylinder functions with an additional $1/\sqrt{2[\pi k_{\rho_{0}}\rho+e^{-(\pi-2)k_{\rho_{0}}\rho}]}$ decay factor modulated by a fast oscillating sinusoidal pattern with period $2\pi/k_{\rho_{0}}$, as explained in sec. 9.3. Fig. 8 shows the initial patterns of the CPG beams for the paraxial case, in which $\Delta\rho_{0}\sim 16\,\mu\text{m}$. In the nonparaxial regime, the patterns are similar, but with much faster oscillations ($\Delta\rho_{0}\sim 0.32\,\mu\text{m}$).

Fig. 9 presents the resulting spectra and the intensity profiles of the odd beam. The first line depicts the spectra, calculated from the inverse Fourier transform of the initial profile (analogously to Eq. (25)). Since the PG beams are solutions of the paraxial wave equation, their spectra are indeed concentrated and, in the case of their CPG counterparts, they are centered at $k_{\rho}=k_{\rho_{0}}$. The second line shows the resulting intensity profiles, normalized by their peak intensities at $z=0$. Although the odd symmetry makes the on-axis amplitude zero at $z=0$, the sidelobes shown in Fig. 8 interfere constructively after some distance and create a very high-intensity peak over the axis, due to the concentration of the energy they carry in a small area. The effect is even more intense in the nonparaxial regime, in which the spot size is much smaller. The last line compares the on-axis intensity predicted by Eq. (65) and the numerical calculation of Rayleigh-Sommerfeld diffraction integral with the initial pattern given by Eq. (66). It is clear that in both paraxial and nonparaxial regimes the analytical expression is accurate.

Fig. 10 presents the same kinds of results of Fig. 9 but for the even beam. We see again that the spectra are concentrated and that the on-axis intensities predicted by the analytical expressions are accurate. However, the even symmetry makes the peak intensities appear at $(\rho,z)=(0,0)$. Therefore, due to the normalization, we do not see higher-intensity on-axis peaks as in Fig. 9.

The field depths of the examples presented can be estimated by Eq. (67) without recurring to the plots of intensity profiles. We could use $w_{0}=r_{0}$ for a reasonable estimate, but better results are obtained by choosing $w_{0}$ based on the envelopes of Fig. 8. For the odd beam, $w_{0}=0.35\,\text{mm}$ embraces all its significant lobes, while the same happens for the even beam for $w_{0}=0.3\,\text{mm}$. Using these values in Eq. (67), we get $Z\approx 35\,\text{mm}$ and $Z\approx 0.56\,\text{mm}$ for the odd beam under paraxial and nonparaxial regimes, respectively, while $Z\approx 30\,\text{mm}$ and $Z\approx 0.48\,\text{mm}$ for the even beam under paraxial and nonparaxial regimes, respectively.

11 Conclusions

In this work, we developed a theoretical analysis to efficiently describe the propagation of superpositions of waves with concentrated wavevector and frequency spectra, thus allowing a simple analytical description of fields with interesting transverse profiles. Starting with an extension of the paraxial formalism, we show how it can be applied to easily handle combinations of waves traveling in different directions without recurring to coordinate rotations. Moreover, and more importantly, we used a similar procedure to create a novel analytical description of azimuthally symmetric waves with concentrated spectra that can leverage all the previously presented results for unidimensional Cartesian waves, allowing us to build azimuthally symmetric versions of them and, therefore, to possibly derive new types of waves. Since these waves are composed of superpositions of zero-order Bessel beams with close cone angles that can be as large as desired, unlike in the paraxial formalism, we showed that it can also provide simple expressions for nonparaxial versions of known beams, such as the Bessel-Gauss beam. Finally, the validity and accuracy of the theory was corroborated by experimental demonstrations and comparisons to numerically-calculated Rayleigh-Sommerfeld diffraction integrals.

12 Funding information

São Paulo Research Foundation (FAPESP) (2015/26444-8); National Council for Scientific and Technological Development (CNPq) (304718/2016-5).

13 Acknowledgments

The authors thank Pedro Paulo Justino da Silva Arantes for his assistance with the numerical calculations of the diffraction integrals.