Accurate expansion of cylindrical paraxial waves for its straightforward implementation in electromagnetic scattering

TL;DR

This paper develops a Bessel series expansion based on the exact Helmholtz solution in cylindrical coordinates, enabling accurate and straightforward modeling of low-numerical-aperture focal waves in electromagnetic scattering.

Contribution

It introduces a novel Bessel series expansion method for precise description of paraxial waves, facilitating easier implementation in scattering problems.

Findings

Validated with Gaussian beams and focused fields

Achieves accurate wave field evaluation in the paraxial regime

Facilitates implementation in electromagnetic scattering simulations

Abstract

The evaluation of vector wave fields can be accurately performed by means of diffraction integrals, differential equations and also series expansions. In this paper, a Bessel series expansion which basis relies on the exact solution of the Helmholtz equation in cylindrical coordinates is theoretically developed for the straightforward yet accurate description of low-numerical-aperture focal waves. The validity of this approach is confirmed by explicit application to Gaussian beams and apertured focused fields in the paraxial regime. Finally we discuss how our procedure can be favorably implemented in scattering problems.

Equations

$$E(x,z)=\frac{\exp (ikz)}{\sqrt{i\lambda z}}{\int }_{-\infty }^{\infty }E(x_0,0)\exp \left[\frac{ik}{2z}{\left(x-x_0\right)}^2\right]d{x}_0.$$

$$E_F^{+}(x,z)=\frac{\exp (ikr)}{\sqrt{r}}\exp (-i\pi /4)\sqrt{\lambda }{a}_F(\theta ),$$

$$a_F(\theta )=\frac{1}{\lambda }{\int }_{\infty }^{\infty }E(x_0,0)\exp \left(-ik\theta x_0\right)d{x}_0,$$

$$E(x,z)=\exp (ikz){\int }_{-\infty }^{\infty }{a}_F(\theta )\exp \left(-i\frac{k}{2}z{\theta }^2\right)\exp \left(ik\theta x\right)d\theta .$$

$$a_F(\theta )=\sum _{n=-\infty }^{\infty }{a}_n\exp (in\theta ),$$

$$a_n=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{a}_F(\theta )\exp (-in\theta )d\theta .$$

$$a_n=\frac{1}{2\pi }E\left(-\frac{n}{k},0\right),$$

$$H_n^{(1)}(kr)\to \sqrt{\frac{2}{\pi kr}}\exp (ikr)\exp (-in\pi /2)\exp (-i\pi /4).$$

$$E_F^{+}(x,z)=\pi \sum _{n=-\infty }^{\infty }{a}_n\exp (in\pi /2)H_n^{(1)}(kr)\exp (in\theta ),$$

$$E_F^{-}(x,z)=\frac{\exp (-ikr)}{\sqrt{r}}\exp (i\pi /4)\sqrt{\lambda }{a}_F({\theta }^{\prime }),$$

$$H_n^{(2)}(kr)\to \sqrt{\frac{2}{\pi kr}}\exp (-ikr)\exp (in\pi /2)\exp (i\pi /4),$$

$$E_F^{-}(x,z)=\pi \sum _{n=-\infty }^{\infty }{a}_n\exp (-in\pi /2)H_n^{(2)}(kr)\exp (in{\theta }^{\prime }).$$

$$E(r,\theta )=2\pi \sum _{n=-\infty }^{\infty }{i}^n{a}_n{J}_n(kr)\exp (in\theta ).$$

$$E(r,\theta )=E_0{J}_0(kr)+4\pi \sum _{n=1}^{\infty }{i}^n{a}_n{J}_n(kr)\cos (n\theta ),$$

$${\mathbf{E}}_{in}=\hat{z}{E}_0\sum _{n=-\infty }^{+\infty }2\pi a_n{i}^n{J}_n\left(kr\right)\exp \left(in\theta \right),$$

$${\mathbf{H}}_{in}=H_0\sum _{n=-\infty }^{+\infty }2\pi a_n{i}^n\left[\hat{r}n\frac{J_n\left(kr\right)}{kr}+\hat{\theta }i{J}_n^{\prime }\left(kr\right)\right]\exp \left(in\theta \right),$$

$${\mathbf{E}}_{sca}=\hat{z}{E}_0\sum _{n=-\infty }^{+\infty }{b}_n{i}^n{H}_n^{(1)}\left(kr\right)\exp \left(in\theta \right),$$

$${\mathbf{H}}_{sca}=H_0\sum _{n=-\infty }^{+\infty }{b}_n{i}^n\left[\hat{r}n\frac{H_n^{(1)}\left(kr\right)}{kr}+\hat{\theta }i{H_n^{(1)}}^{\prime }\left(kr\right)\right]\exp \left(in\theta \right),$$

$$b_n=-2\pi a_n\frac{J_n(kR)}{H_n^{(1)}(kR)},$$

$$E_{PW}(r,\theta )=\sum _{n=-\infty }^{\infty }{i}^n{J}_n(kr)\exp (in\theta ).$$

$$E_{GB}(x,z)=\sqrt{\frac{q(0)}{q(z)}}\exp \left[ik\left(z+\frac{x^2}{2q(z)}\right)\right],$$

$$E_{GB}(x,z)\to \sqrt{\frac{z_R}{iz}}\exp \left[ik\left(z+\frac{x^2}{2z}\right)\right]\exp \left[-{\left(\frac{x/z}{w_0/z_R}\right)}^2\right].$$

$$E_{AF}$$

$$E_{AF}(x,0)=\frac{\sin \left(k{\theta }_{0}x\right)}{k{\theta }_{0}x}.$$

Full text

††journal: A Journal

1 Introduction

The description of unbounded wave fields in the form of series expansion represents a widely used tool in free space propagation [1, 2, 3]. For instance, the electric field emerging from a laser device is commonly expressed as a combination of either Hermite-Gauss or Laguerre-Gauss basis functions taken from the natural modes existing in mirror cavities [4, 5, 6]. Such procedure, however, seems inappropriate for apertured focal waves where diffraction sidelobes will impose a strong limitation in its efficient characterization by an acceptable truncation of the sequence [7, 8]. In these cases, expansions in terms of Lommel functions demonstrate a suitable implementation [9, 10].

In planar arrangements, two-dimensional (2D) waves satisfying the Helmholtz equation can be described in terms of Bessel wave functions when they are formulated in cylindrical coordinates [11]. This is a method followed for instance in the Mie-Lorenz theory applied to scatterers with circular cross section [12, 13, 14]. Only a small part of the cases analyzed in scattering problems consider non-uniform beams, among other reasons due to the complexity of the series expansions with cylindrical vector wave functions (also spherical vector wave functions) that can be found [15, 16, 17, 18, 19, 20]. Particularly interesting is the paper from Shogo Kozaki which analyzes the case of a Gaussian beam illuminating a cylindrical particle, where it is possible to remarkably simplify the analytical description of the beam, under certain approximations, in terms of Bessel cylindrical waves [21].

A salient feature of these Bessel wave functions is their inherent localization around the chosen origin of coordinates. Therefore, focal waves are generally expected to be favorable candidates to be effectively represented by means of a series expansion using Bessel wave fields.

In this paper, we consider 2D paraxial wave fields which are localized in the vicinity of a given point, which serves for the construction of the Bessel wave-functions basis to be used in a series expansion. The coefficients of such sequence are estimated by analyzing the far field of the paraxial wave. This choice reduces the resultant calculation to a simple Fourier expansion. The validity of this method is finally verified when applied to a Gaussian laser beam and also to an apertured focal wave.

2 Theoretical analysis

Let us consider a 2D wave field $E(x,z)=\exp(ikz)U(x,z)$ propagating in free space and satisfying the paraxial wave equation, $\partial_{x}^{2}U+2ik\partial_{z}U=0$, where $z$ represents the spatial coordinate along the optical axis, $x$ is the transverse coordinate and $k=2\pi/\lambda$ is the wavenumber. Note that a harmonic time variation $\exp(-i\omega t)$ is here assumed, where $\omega$ is the time-domain frequency. If the focal point is set at the origin of coordinates, $(x,z)=(0,0)$, the wave field at any out-of-focus plane $z\neq 0$ can be determined by means of the Fresnel diffraction integral [22]

[TABLE]

The characteristic Fraunhofer pattern can be observed well beyond the focal plane in the limit $|z|\to\infty$. In this case, the quadratic phase term $\exp\left({ik}x_{0}^{2}/{2z}\right)$ in the integral equation (1) is negligible.

When $z>0$, the far field of $E(x,z)$ can be expressed as

[TABLE]

where the angular spectrum

[TABLE]

is simply the Fourier transform of the wave field at the focal plane. The azimuthal angle measured from the optical axis is given as $\theta=x/z$ in the paraxial approximation, whereas the radial coordinate is given by $r=|z+x^{2}/2z|$. In the paraxial regime, the wave function $a(\theta)$ has significant values only when $\theta\ll 1$. Furthermore, the focal field can be evaluated once the angular spectrum $a_{F}(\theta)$ is known as

[TABLE]

This diffraction integral is well known within the Debye theory and has been successfully applied to apertured focal waves, as we will consider below [23].

Now it is clear that a wave expansion can be given by means of a Fourier series of the angular spectrum, that is

[TABLE]

where the Fourier coefficient

[TABLE]

Since the angular spectrum $a(\theta)$ takes values identically zero in the range $|\theta|>\pi/2$ in the semi-space $z>0$, the definite integral given in Eq. (6) can be further simplified by extending the interval of integration to $-\infty<\theta<\infty$. Under such approximation, and substituting Eq. (3) into (6), we finally infer that the Fourier coefficient of order $n$

[TABLE]

depends on the focal field as measured at the off-axis point $x=-n/k$.

To extend the series expansion of the wave field into the near field (note that here we disregard evanescent waves and near field refers to waves in the vicinities of focus), we will use the Bessel wave functions of the first kind and order $n$, $J_{n}(kr)$, and the Bessel functions of the second kind and order $n$, $Y_{n}(kr)$, which are solutions of the 2D Helmholtz wave equation, $\nabla^{2}E+k^{2}E=0$, provided that the angular variation of the wave field is given by $\exp(in\theta)$ [11]. In particular, the Hankel wave function of the first kind $H_{n}^{(1)}(kr)=J_{n}(kr)+iY_{n}(kr)$ has an asymptotic limit far from the origin of coordinates given as

[TABLE]

Therefore, an outgoing cylindrical wave field with a specific angular momentum can be given in terms of a Hankel wave function of the first kind and unique order $n$. As a consequence, a wave field exhibiting an angular spectrum $a_{F}(\theta)$ in the semi-space $z>0$ can be expressed as a combination of Hankel wave functions of different orders $n$ as

[TABLE]

which in addition can be utilized in the near field.

Equation (9) provides the expression of a wave field whose Fraunhofer pattern is given in terms of the wave function $a_{F}(\theta)$, leading to a vanishing far field in $z<0$. However, in such semi-space, the far field of $E(x,z)$ given in Eq. (1) can be expressed as

[TABLE]

where the angle $\theta^{\prime}=x/z$ in the paraxial approximation. With the new angular coordinate $\theta^{\prime}=\theta-\pi$, we again may describe the Fraunhofer pattern observed at $z\to-\infty$ by means of the wave function $a_{F}(\theta^{\prime})$ which takes significant values for $|\theta^{\prime}|\ll 1$. Furthermore, the Fourier series given in Eqs. (5) and (6) are still applicable in this case. Taking into account the asymptotic limit of the Hankel wave function of the second kind $H_{n}^{(2)}(kr)=J_{n}(kr)-iY_{n}(kr)$, written as

[TABLE]

we may infer a series expansion of the wave field in terms of incoming cylindrical waves, resulting in

[TABLE]

The appropriate description of the paraxial wave field $E(x,z)$ includes a far field at $z\to\infty$ in the form of an outgoing cylindrical wave patterned by $a_{F}(\theta)$, and simultaneously at $z\to-\infty$ representing an ingoing cylindrical wave shaped by the same angular spectrum $a_{F}(\theta^{\prime})$. Note that such symmetry around the focal point has been analyzed by Collet and Wolf [24]. As a consequence, the wave field must be computed as $E(x,z)=E_{F}^{+}(x,z)+E_{F}^{-}(x,z)$, which in cylindrical coordinates is expressed as

[TABLE]

This is the main result of our study, together with the fact that the Fourier coefficients as given in Eq. (7) can be achieved in terms of the focal wave field. We conclude that a cylindrical, paraxial wave field can be accurately described by a series expansion sustained by Bessel wave functions which are solutions of the 2D Helmholtz equation, and with expansion coefficients which are determined by the own wave field at specific points of the focal plane. We point out that Eq. (13) accurately provides the wave field even at angles far from the paraxial regime were $\tan\theta=x/z$; of course $r=\sqrt{x^{2}+z^{2}}$ denotes the distance from focus to the observation point.

In centrosymmetric field distributions where $E(-x,z)=E(x,z)$, what occurs simply if the field is symmetric with respect to the origin of coordinates at the focal plane, a reduced expression of the series expansion can be deduced. In this case, $a_{-n}=a_{n}$ as inferred from Eq. (7). By using the following property of Bessel functions of negative order, $J_{-n}(\alpha)=(-1)^{n}J_{n}(\alpha)$, we finally may reduce Eq. (13) to

[TABLE]

where $E_{0}=2\pi a_{0}$ is the in-focus wave field.

3 Implementation in scattering problems

Let us consider a scatterer embodied in a cylindrical region whose axis is set along the $y$ axis and with a radius $R$, as depicted in Fig. 1. In order to analytically estimate the scattered wave field of the body, we may follow the Lorenz-Mie scattering method which is described in detail for instance in Refs. [12, 13, 25]. Here we consider TMy polarized waves in such a way that the electric field is oriented along the $y$ axis. Note that the case of TEy polarized waves where the magnetic field lies on the $y$ axis is straightforwardly treated by using the duality theorem [11]. Using Eq. (13), the electromagnetic field of the incident wave field may be written as

[TABLE]

and $\mathbf{H}_{in}=-i{\nabla\times\mathbf{E}_{in}}/{\omega\mu_{0}}$, which reads as

[TABLE]

where $H_{0}={E_{0}k}/{\omega\mu_{0}}$ is a constant amplitude. Here the prime appearing in $J^{\prime}_{n}\left(\alpha\right)$ denotes derivative with respect to the variable $\alpha$. The scattered electromagnetic field in the environment medium, $r>R$, may be set as [13]

[TABLE]

and

[TABLE]

where the coefficients $b_{n}$ must be determined. The total electric field in the environment medium is simply $\mathbf{E}_{tot}=\mathbf{E}_{in}+\mathbf{E}_{sca}$.

The Lorenz-Mie scattering coefficients $b_{n}$ are determined by means of the proper boundary conditions established at $r=R$, which are continuity of tangential components of the electromagnetic field. For illustrative purposes, let us consider a perfect electric conductor at the boundary. Therefore, the tangential component of $\mathbf{E}_{tot}$, which is oriented along the $y$ axis, must vanishes at $r=R$. This finally yields

[TABLE]

where the coefficients with different index $n$ can be treated separately due to the linear independence of the multipolar components of the wave field.

4 Practical application of the series expansion

The Bessel wave expansion can be applied to well-known patterns to confirm the validity of our approach. First we consider a plane wave propagating along the $z$-axis, where the normalized field $E_{PW}(x,0)=1$ is uniform on the “focal” plane. In this case, the far field wave function $a_{F}(\theta)=\delta(\theta)$, corresponding to the Dirac delta function. Therefore, the Fourier coefficients take the constant value $a_{n}=1/2\pi$ apart from the index $n$, which in addition is consistent with the result attained by using Eq. (7). Finally, the wave field expansion of $E_{PW}(x,z)=\exp(ikz)$ can be set as

[TABLE]

Note that such expression can also be determined by using the Jacobi-Anger identity as commonly carried out in scattering problems in two dimensions [13]. It is noteworthy that precisely a plane wave requiring an infinite number of terms in the series expansion provides, on the other hand, an exact solution to the Helmholtz equation.

4.1 Gaussian beams

A Gaussian beam in two dimensions can be described as an analytical solution to the paraxial wave equation. Provided that the normalized wave field at the plane of the beam waist, that is the focal plane, is given by $E_{GB}(x,0)=\exp\left(-x^{2}/w_{0}^{2}\right)$, where $w_{0}$ is the Gaussian beam width, the field distribution in the $xz$ plane is then evaluated as [3]

[TABLE]

where the complex radius $q(z)=z-iz_{R}$ and the Rayleigh range $z_{R}=\pi w_{0}^{2}/\lambda$. The Fraunhofer pattern can be observed in the far field region $|z|\to\infty$ where

[TABLE]

Therefore, the angular spectrum of the Gaussian beam yields $a_{F}(\theta)=1/(\sqrt{\pi}\theta_{0})\exp\left(-{\theta^{2}}/{\theta_{0}^{2}}\right)$ which also shows a Gaussian distribution with angular width given by the so-called far-field beam angle $\theta_{0}=\lambda/\pi w_{0}$. Finally, the Fourier coefficients of the series expansion are then written as $a_{n}=1/2\pi\exp\left(-n^{2}\theta_{0}^{2}/4\right)$. Note that the resulting wave expansion of the Gaussian beam by applying our procedure is analogous to that found in Ref. [21] following a completely different approach.

Figure 2 shows the field distribution of a Gaussian beam of waist $w_{0}=10\lambda$, which corresponds to a far-field beam angle $\theta_{0}=1.8^{\circ}$, evaluated in the plane of the beam waist $z=0$. For the sake of comparison, we use the series expansion given in Eq. (14) which is truncated by a finite Bessel order, $n\leq n_{\mathrm{max}}$, the latter here called truncation index. When $n_{\mathrm{max}}=10$, the wave field is accurately evaluated only in the vicinities of the focal point, and deviations of the truncated series expansion are evident if $|x|>1.5\lambda$, as shown in Fig. 2(a). This result can be explained by taking into account that the Bessel wave field $J_{n}(kr)$ features a central circular shadow bounded by a main peak set at a distance roughly estimated as $r_{n}=|n|/k$, together with ripples dispersed further from focus [26]. As a consequence, the Bessel expansion of Eq. (14) which is subsequently truncated at $n_{\mathrm{max}}$ provides a rigorous estimation of the wave field at points of the focal region set well inside the circle of radius ${n_{\mathrm{max}}}/k$. For $n_{\mathrm{max}}=10$, $50$, $100$, and $150$ we estimate a boundary of radius $1.6\lambda$, $8.0\lambda$, $15.9\lambda$, and $23.9\lambda$ respectively, which is in good agreement with the results displayed in Figs. 2(a)-(d). Since the focal field of the Gaussian beam is negligible at $|x|>3w_{0}$, we finally conclude that $n_{\mathrm{max}}=3kw_{0}$ (that is $n_{\mathrm{max}}=6/\theta_{0}$) yields a sufficient number of Bessel wave functions to be used in the series expansion. Note that for $w_{0}=10\lambda$ we estimate that a number of elements higher than $n_{\mathrm{max}}=60\pi\approx 188$ in the Bessel expansion will not improve significantly the evaluation of the wave field in the focal plane.

Inasmuch as the Bessel wave functions constituting the series expansion exhibit radial symmetry, it is reasonable to identify equivalent restrictions on the sequence truncation when evaluating the field distribution of the Gaussian beam along the optical axis. The latter is proved in Fig. 3(a), which shows the magnitude of the series expansion at $x=0$ for a Gaussian beam of $w_{0}=10\lambda$, considering a different number of Bessel wave functions used in the truncated sequence. The exact value $|E_{GB}(0,z)|=\left[1+(z/z_{R})^{2}\right]^{-1/4}$ is also represented for the sake of comparison [3]. Due to the slower decay of the Gaussian field distribution along the optical axis, which is characterized by the Rayleigh range ($z_{R}=314\lambda$ is much longer than the beam width), it might occur that an accurate computation of the on-axis field distribution using the Bessel series expansion, at points reaching significant values, would require a higher number $n_{\mathrm{max}}$. To verify the validity of this argument, Fig. 3(b) shows a semi-log plot of the magnitude derived from the series expansion from focus to a point located well beyond the Rayleigh range. Though using $n_{\mathrm{max}}=100$ it is expected a significant deviation of our calculation moving away to a distance further than roughly $16\lambda$, we surprisingly obtain a tolerable estimation up to $3\times 10^{4}\lambda$. The strong sidelobes featuring the Bessel function also play a key role in the determination of the wave field, particularly on the optical axis. We conclude that the truncation index $n_{\mathrm{max}}$ indicating the number of elements used in the truncated series expansion does not strictly determines the region around the focus where the Eq. (14) is applicable but the working precision of our calculation.

The above discussion can include the analysis of the phase distribution of the Gaussian beam along its optical axis, as shown in Fig. 4(a) and (b). Note that the exact phase distribution is given by $-\phi(z)/2$, where $\tan\phi=z/z_{R}$, provided that the fast-evolving term $\exp(ikz)$ is ignored for the sake of clarity [3]. The Gouy phase shift accumulated when passing through the focus from $z\to-\infty$ to $z\to+\infty$ is $-\pi/2$, a fact that is primarily carried out in the range $|z|<z_{R}$. As evidenced in Fig. 4(b), the series expansion (14) always provides the correct Gouy phase shift independently of the applied truncation. However, underestimating the value of $n_{\mathrm{max}}$ will reduce the path wherein such a phase shift is naturally induced.

Figures 5-7 illustrate the convergence of the series expansion (14), by varying the truncation index $n_{\mathrm{max}}$, which is evaluated at different points of the first quadrant in the $xz$ plane. The chosen points have one of these three different radii: $r=\{\lambda,10\lambda,50\lambda\}$, given at three different angular orientations, $\theta=10^{\circ}$, $45^{\circ}$, and $80^{\circ}$. For each one of these nine points we represent the magnitude and phase of the complex wave field separately. When the computation of the field is performed close to focus, at $r=\lambda$, a fast convergence of the series expansion is reached where $n_{\mathrm{max}}=10$, which is a few integers higher than the product $kr$, is enough to accurately determine the complex field. Such behavior does not depend on the azimuthal angle $\theta$, as shown in Figs. 5(a)-(c). By analyzing the points set at $r=10\lambda$, where now the magnitude of the field drops more significantly (yet remains in the same order) at higher values of $\theta$, that is closer to the focal plane than the optical axis, the series expansion on the other hand converges isotropically. In Fig. 6 we recognize that a truncation index $n_{\mathrm{max}}=70$ is appropriate for all angles; note that here $kr\approx 63$.

However, a dramatic variation of the series convergence occurs if we draw attention to points at a distance $r=50\lambda$ far from focus, as shown in Fig. 7. In this case, $n_{\mathrm{max}}=100$ suffices for the point set near to the optical axis at $\theta=10^{\circ}$; on the contrary the truncation index $n_{\mathrm{max}}$ should surpass a value of $300$ when $\theta=80^{\circ}$ in order to achieve an accurate estimation of the wave field. Importantly, the magnitude of the Gaussian field $|E_{GB}|$ reaches a value of $0.48$ in the first case, falling to $|E_{GB}|=2.5\times 10^{-11}$ in the latter, which is a consequence of the natural resistance of paraxial fields to decay along the direction of propagation. Therefore, the requirement establishing a truncation index lightly higher than the parameter $kr$ ($\approx 314$ here) is greatly moderated near the optical axis where the field reaches its highest values by several orders of magnitude.

4.2 Apertured focused fields

Let us consider an apertured focal wave that can be produced for instance by a lens system, where in addition pupil truncation is taken into account. In this case, the angular spectrum might be assumed to be flat, which is expressed as $a_{F}(\theta)=a_{F}(0)$ within the interval $|\theta|<\theta_{0}$; otherwise the angular spectrum of the wave field vanishes. Note that $\theta_{0}$ represents the numerical aperture of the converging wave in the paraxial regime, provided that propagation takes place in free space. The wave field in the focal volume can be evaluated by means of the Debye diffraction integral given in Eq. (4), which yields

[TABLE]

where the normalized spatial coordinate $\zeta=2\theta_{0}^{2}z/\lambda$ and $\tilde{F}^{*}(\cdot)$ is the complex conjugate of the complex Fresnel integral function [3]. Also we used $a_{F}(0)=1/2\theta_{0}$ to obtain a normalized in-focus wave field $E_{0}=1$. Specifically at the focal plane $z=0$, the wave field reduces to

[TABLE]

The first zero of the wave field which is found at $x_{1}$ around the focal point satisfies $k\theta_{0}x_{1}=\pi$ and determines the limit of resolution following the Rayleigh criterion. If we conserve the definition of the beam width previously used for Gaussian beams ($w_{0}=\lambda/\pi\theta_{0}$) as measured at the focal plane, we may write $x_{1}=(\pi/2)w_{0}\approx 1.57w_{0}$. Therefore, here $w_{0}$ can also be considered as the beam width of the apertured focused wave. Finally, the Fourier coefficients to be employed in the series expansion (14) now read as $a_{n}={\sin\left(n\theta_{0}\right)}/{\left(2\pi n\theta_{0}\right)}$.

The focal field of the apertured converging wave is represented in Fig. 8 for different numerical apertures $\theta_{0}$, together with its corresponding Bessel series expansion which is truncated at $n_{\mathrm{max}}=70$. The truncated series expansion proves an accurate evaluation of the wave field from focus to a distance $~{}11\lambda$, which on the other hand is roughly the estimator $n_{\mathrm{max}}/k$. When the numerical aperture is low, it may happen that the region of validity of the series expansion even cannot reach the central lobe of the sinc pattern, as illustrated in Fig. 8(a) and (b) for $\theta_{0}=1^{\circ}$ and $\theta_{0}=2^{\circ}$, respectively. This peak of highest intensity will be conveniently reproduced by the series expansion (14) provided that $n_{\mathrm{max}}>kx_{1}$, occurring if $\theta_{0}>\pi/n_{\mathrm{max}}$ ($\theta_{0}>2.6^{\circ}$ for $n_{\mathrm{max}}=70$). In Fig. 8(c) we consider the increased numerical aperture $\theta_{0}=5^{\circ}$, which associated Bessel expansion covers not only the central lobe but the first sidelobe accurately. These sidelobes are determined by the zeros of the sine function found at $x_{m}=m/k\theta_{0}$, thus a given truncation index ensures an accurate reproduction of the wave field containing the central peak and $m-1$ sidelobes determined by $m=n_{\mathrm{max}}\theta_{0}/\pi$. Specifically, for $\theta=5^{\circ}$ we appropriately represent up to $m-1\approx 1$ sidelobes, whereas such quantity is tripled for $\theta=10^{\circ}$, as shown in Figs. 8(c) and (d).

The analysis of the on-axis field distribution is illustrated in Fig. 9. The complex wave field depends on the normalized spatial coordinate $\zeta$ exclusively, leading to a magnification of the pattern in direct proportion to $1/\theta_{0}^{2}$, in opposition to the transverse pattern whose magnification depends on $1/\theta_{0}$. When the series expansion is truncated at $n_{\mathrm{max}}=150$, accurate estimations are derived within the interval ranged from focus to roughly $20\lambda$ (note that $n_{\mathrm{max}}/k=23.87\lambda$). Out of that range fast oscillations in phase and magnitude make the series expansion deviates from the exact on-axis field distribution, as evidenced in Fig. 9(a)-(d). However, the amplitude of such oscillations are notably reduced when the numerical aperture $\theta_{0}$ increases. Nevertheless, an adequate field decay at much longer distances as $1/\sqrt{z}$ is reproduced by the series expansion even at the lowest numerical aperture shown in Fig. 9(a). The exact phase distribution now exhibits a strong linear variation, conjoined with the central peak of intensity, along with some maxima and minima around the asymptotes $\pm\pi/4$. The series expansion provides a good estimation of the slope in the central linear region even in the case of $\theta_{0}=1^{\circ}$ where $n_{\mathrm{max}}=150$ remains suboptimal to accurately mimic the whole central peak of the field pattern. The asymptotes are also well established in the latter case. By increasing the value of the angle $\theta_{0}$ we can in addition follow appropriately the profile of the continuous side lobes that appear when we move further from focus. When $\theta_{0}=10^{\circ}$, the central peak and the first sidelobe of the intensity profile are conveniently determined by the series expansion, as well as the focal linear variation of the phase distribution and its three closest lateral extrema, as shown in Fig. 9(e) and (f).

5 Conclusions

Based on the far field characteristics of paraxial beams, we are able to represent the full wave field in a series expansion in terms of Bessel wave functions which are solutions of the Helmholtz equation in two dimensions. The coefficients of the series expansion in fact stem from a Fourier expansion of the wave-field angular spectrum, and they can be analytically determined by using the pattern in the focal plane wherein the field is mostly localized. This representation of the wave field is adequate for scattering problems where a scatterer is localized in a circular section around the origin of coordinates, as described thoroughly.

Application of the Bessel expansion to a uniform plane wave leads us to one expression which, alternatively, may be inferred from the well-known Jacobi-Anger identity, demonstrating the validity of our approach. Other paraxial beams have been analyzed in detail. Firstly a Gaussian beam has been critically examined, particularly in comparison with a practical truncation given by the index $n_{\mathrm{max}}$ of the corresponding series expansion. The wave field is accurately evaluated at points whose distance from the origin of coordinates is shorter that $n_{\mathrm{max}}/k$. Provided that such origin of coordinates is established at the maximum of intensity, then a truncation index $n_{\mathrm{max}}=3kw_{0}$ (where $w_{0}$ is the beam waist) should be enough for the series expansion to provide a good estimation of the wave field in the whole space.

By analyzing apertured wave fields we further introduce the main aspects of the series truncation. For instance, intensity sidelobes existing in the vicinity of the focal peak are critically reproduced under the appropriate truncation index. Due to the specific characteristics of the paraxial beams, such a circumstance is much more evident transversally to its direction of propagation. Nevertheless, both the inverse root squared decay of the field with distance from focus, and the $\pi/2$ phase shift undergone on axis for general wave fields is always implicit in the Bessel-based series expansion. Finally, this procedure can be applied to any kind of two-dimensional paraxial beam, in particular when it will interact with particles set in its trajectory. Such scatterers are not restricted in size, ranging from nanowires to large cylinders, which in addition can be made of anisotropic and even multilayered materials [27, 28].

Acknowledgments

This work was supported by the Spanish Ministry of Economy and Competitiveness (MINECO) (TEC2014-53727-C2-1-R).