1. Introduction

The description of nonparaxial monochromatic scalar fields requires the use of efficient computational schemes. Customarily, the field is computed by expressing it as a superposition of closed-form solutions of the Helmholtz equation, like plane waves, spherical waves, multipoles, Bessel beams, etc. These superposition schemes often require the evaluation of oscillatory integrals or slowly-converging series for the computation of the field at each point. For special cases, there are bases that are particularly convenient. For example, in the paraxial limit, one can expand a field in terms of (Hermite or Laguerre) Gaussian beams. In the highly nonparaxial case, on the other hand, the multipolar basis gives rapidly-convergent solutions. The goal of this paper is to develop a new scheme that is efficient between these two limiting situations.

The scheme proposed here is valid for free fields with rotational symmetry, and is based on a different class of closed-form solutions of the Helmholtz equation. These solutions are built in terms of what we call here complex focus fields (CFF’s) [1, 2, 3, 4, 5, 6], which correspond to spherical waves focused at a complex point. The CFF’s have rotational symmetry around the direction of the imaginary part of their focal position and reduce to paraxial Gaussian beams when this imaginary part is much larger than the wavelength. The elements of the new basis are the result of transforming the angular spectra of the multipolar basis according to a simple geometrical mapping of the sphere of directions. This mapping is chosen such that the monopole transforms into a CFF. The expansion coefficients for a general field are calculated through a simple integral involving the field’s angular spectrum. The magnitude of the imaginary focus of the CFF is then a free parameter that can be chosen to optimize the convergence of the expansion. In particular, in the limit when this parameter is large, one obtains an orthonormal basis for expanding paraxial beams, which can be written in terms of Legendre polynomials with Gaussian arguments. The resulting paraxial scheme is reminiscent of the one proposed by Wen and Breazeale [7, 8], where a paraxial acoustical beam with rotational symmetry is expanded in terms of coaxial Gaussian beams of different spot sizes, whose coefficients are obtained through a numerical optimization procedure. A similar expansion was proposed by Li [9] for describing a paraxial flat-top optical beam. For fields that deviate significantly from the paraxial regime, a smaller value can be chosen for the above-mentioned parameter, so that the new scheme gives rapidly-convergent results. In fact, in the limit of highly nonparaxial fields, this parameter can be set to zero, so that the new basis reduces to the multipolar one.

2. Real and complex multipoles

2.1. Debye representation of scalar fields

A general monochromatic scalar field in three dimensions with no evanescent components can be written as a sum of propagating plane waves of the form

where k is the wave number, A(·) is the angular spectrum of the wave field, and the integral is over the sphere of directions of the unit vector u. The angular spectrum, which is a complex distribution over the unit sphere of directions, fully characterizes the wave field. For simplicity, we set k = 1 in all subsequent formulas, i.e., all distances are in units of reduced wavelengths.

2.2. Multipolar fields

The spherical harmonics of order l and index m are defined by [10]

where θ and ϕ are the polar and azimuthal angle, respectively, and $P_l^{\left(m\right)}$(·) is the associated Legendre function with integer indices l and m [10]. The spherical harmonics constitute a set of orthonormal functions across the unit sphere of directions, since

The angular spectrum of any free monochromatic wave field A(u) can then be written as

where a_l,m is found through the inner product

It turns out that the field corresponding to an angular spectrum Y_l,m (u) is the scalar multipole

where r = |r| and j_l (·) denotes the lth-order spherical Bessel function [10]. In particular, a perfect spherical wave (or monopole) is characterized by the following angular spectrum and field distributions:

2.3. Complex focus fields

As can be seen from Eq. (1), shifting the spatial coordinate by r ₀ is equivalent to multiplying the angular spectrum by a phase factor, i.e.,

This transformation is analogous to the well known shift-phase property in Fourier analysis. Let us consider the spherical wave in Eq. (7). From Eq. (8), the spectrum of a spherical wave focused at a point r ₀ is proportional to exp(-iu · r ₀) and the corresponding field is proportional to sin |r - r ₀|/|r - r ₀|.

Suppose now that the focus is at the imaginary point r ₀ = (0,0, iq), where q ≥ 0. Accordingly, the new angular spectrum, say A _0,0(u;q), is

where we used u_z = cosθ, and a constant factor was introduced in order for the following normalization condition to be fulfilled:

The corresponding field, say U _0,0(r;q) is given by

These fields are the CFF’s mentioned in the introduction. Note that the choice of branch for the complex square root in Eq. (11) is not important, as the Taylor expansion around r = 0 of (sinr)/r contains only even powers. For q > 0, the spherical symmetry is replaced by rotational symmetry around the z-axis, which becomes the main direction of propagation. In fact, when q ≫ 1, A _0,0(u) becomes a Gaussian in the polar angle, i.e.,

and the field U _0,0 becomes a paraxial Gaussian beam, i.e.,

Higher-order multipoles can also be displaced to the same imaginary location, giving them a directionality that is controlled by q. However, this operation of displacement to an imaginary point destroys the orthogonality of the basis. In the next section we propose an alternative way to transform the multipolar basis such that orthogonality is preserved and that the lowest-order element is the CFF in Eq. (11).

3. Mapping of the unit sphere

The basic idea of this work is to use a geometrical mapping of the unit sphere onto itself in order to transform the spherical harmonics Y _l,m (u) onto a new orthonormal set of angular spectra, say A_l,m (u;q). Mappings on the sphere are not unusual in physics. Think, for example, of the Bargmann mapping of velocity directions associated with a relativistic boost [11], or the mapping of ray directions between the two foci of an ellipsoidal reflector.

Consider two spherical coordinate reference frames, say Ω = (θ,ϕ) and Ω′ = (θ′, ϕ′). The mapping is defined by the following pair of equations:

Let us assume that the ϕ -dependence is preserved, so that Φ¯(ϕ) = ϕ and the function Θ¯(·) depends only on θ. This function can be found through an “energy conservation” approach. In particular, we require that Θ¯ maps a uniform distribution onto the distribution of Eq. (9). For this purpose, consider two infinitesimal solid angles in the two spaces, say dΩ′ and dΩ, respectively, and let the “energy”, i.e., the square modulus of the respective angular spectra, contained within them be equal, that is

By using Eqs. (7) and (9), Eq. (15) leads to the following differential equation:

Equation (16) can be easily solved. (Note that the corresponding equation for the two-dimensional case has no simple solution.) We must also impose, however, the condition that θ′ = 0 is mapped onto θ=0. This yields an expression for θ′ = Θ¯ given by

The effect of this geometrical transformation is illustrated in Fig. 1, which shows meridional sections of the sphere of directions. The squared modulus of A _0,0(u;q), as a function of θ, corresponds to the length of the blue radial lines. For q = 0, |A _0,0(u;q)|² is uniform, and the radial lines are equally spaced. As q increases, the radial lines drift towards the positive z-direction, where |A _0,0(u;q)|² peaks.

 

Fig. 1. Mapping of the angular spectra over a meridional section of the sphere of directions. The lengths of the blue radial line segments are proportional to |A _0,0(u;q)|². For increasing q, these radial lines bunch up around the positive z axis. [Media 1]

Download Full Size | PDF

We now introduce the complete basis of spectra A_l,m (u;q) as the result of applying the mapping defined above to the spherical harmonics Y_l,m (u) of Eq. (2), i.e.,

Figure 2 shows A _l,0, for l = 0,1,2,3,4, as a function of θ/π, for q running from 0 to 4. Because the mapping was defined through an energy-conservation approach, it preserves the inner product between two arbitrary angular spectra. Therefore, the functions A _l,m (u;q) are orthonormal on the unit sphere for any q. Furthermore, in the limit q → 0, they reduce to Y _l,m(u). The opposite limit is discussed in the next section.

4. Paraxial limit: Legendre-Gauss basis

Since for large values of q we expect the angular spectrum to be strongly localized around θ ≈ 0, we can use the approximation cos Θ¯(θ)≈-1 + 2η ² with η = exp(-q θ ²/2). The angular spectrum then reduces to

 

Fig. 2. Plots of A _l,0, for l = 0,1,2,3,4, as functions of θ/π, for q running from 0 to 4. [Media 2]

Download Full Size | PDF

which can be written as

where the orthonormal set of functions B_l,m (·,·) are referred to here as the Legendre-Gauss basis:

For m = 0, $P_l^{\mathit{\left(}m\mathit{\right)}}$ (·) reduces to a polynomial of degree l. Thus B _l,0 is a finite sum of Gaussian functions, so that the angular spectrum in Eq. (20) gives rise to a paraxial field composed of a discrete set of coaxial fundamental paraxial Gaussian beams having different spot-sizes.

Figure 3 shows plots of B _l,0’s, for l = 0,1,2,3,4, as functions of τ. Notice that, unlike the Laguerre-Gaussian functions whose support is proportional to 2l + 1, the Legendre-Gaussian functions have roughly a fixed support. They are therefore not a convenient basis for expanding functions with larger supports. To show this, consider the orthogonal expansion of the function circ(τ/a) in terms of Legendre-Gaussian functions. This function, which is independent of ϕ, takes the value 1 for τ < a and 0 otherwise. Because of its rotational symmetry, only the Legendre-Gaussian functions B _l,0 are needed in the expansion. Due to the orthonormality of the basis, the rms error, say ε, caused by the truncation of the series can be calculated in terms of the expansion coefficients, say c _l , in the following form:

where l _maxdenotes the truncation order and N ² = π a ² equals the total energy of the circ function. Figure 4 shows ε as a function of a for different values of the truncation order l _max. Notice that the value of a that minimizes the error for a fixed truncation order is roughly 1.3, independently of l _max. Suppose now that we have a field whose angular spectrum is a circ function of width θ _max. From Eq. (20) and Fig. 4 we easily see that, in order to optimize the convergence, q must be chosen such that θ _max √q ≈ 1.3, so that q _opt ≈ √3/${\theta }_{\text{max}}^2$. It is worth noting that even the minimized errors in Fig. 4 are relatively large, and decay slowly with increasing l _max. This is due to the discontinuity of the circ function. Other bases, like the Hermite-Gaussian and the Laguerre-Gaussian ones, also converge slowly when there are discontinuities. [12, 13]

 

Fig. 3. Plots of B _l,0’s, as functions of τ, for several values of l.

Download Full Size | PDF

 

Fig. 4. Behavior of the rms error ε, as a function of a, for different truncation orders l _max.

Download Full Size | PDF

5. Application to the computation of focused axis-symmetric fields

In the previous section we proved that any paraxial field with axial symmetry can be expressed in terms of a superposition of fundamental Gaussian beams. The aim of this section is to extend this superposition scheme to the nonparaxial regime.

5.1. Orthonormal expansion

Let us start by investigating whether the field distributions U_l,m (r;q), whose angular spectra are the functions A_l,m (u;q), can be evaluated in closed form. From Eq. (1) these fields are given by

The evaluation of these integrals is not trivial except for the case m = 0, when the angular spectrum reduces to

Because P _l is a polynomial, the expression in Eq. (24) is a sum of terms proportional to exp[(2n+1)qcosθ]. Therefore, the corresponding fields take the form

where U _0,0 is the closed form given in Eq. (11), and

with coef[Pol(μ),μⁿ ] denoting the coefficient of the nth power of μ in the polynomial expression Pol(μ).

In the case of axis-symmetric fields, the functions A_l (θ;q) = A _l,0(u;q) constitute an ortho-normal complete basis over the sphere. Any axis-symmetric field can then be calculated as a superposition of U _l,0’s. For an angular spectrum A(θ), the expansion coefficients c_l ’s are found as

Rather than directly using Eq. (27), one can evaluate the coefficients by using the inverse map, say θ = Θ(θ′), on A(θ). From Eq. (17) we find cos θ = cos[Θ(θ′)], where

From Eqs. (27) and (28), and taking Eq. (16) into account, we get

where

Once coefficients c_l are found, the field can be written as

Equation (31) gives the desired expansion of the field as a series involving simple closed form expressions. The parameter q can be chosen to make the series as rapidly convergent as possible.

5.2. Example: Debye focused fields

As a first example, consider the field produced by a spherical monochromatic wave emerging from a circular aperture and converging towards the axial focal point. In the Debye approximation, where the exit pupil is considered to be infinitely far away, the field can be written as U(r) = ∫_Ω0 exp(iu · r)dΩ, where Ω₀ = {θ ∊ [0, θ _max];ϕ ∊ [0,2π)}, and sin θ _max is the ratio between the radius of the aperture and the focal distance [14]. The angular spectrum is simply given by A(u)=circ(θ/θ _max).

As pointed out in the previous section, in order to improve the numerical convergence of the expansion in Eq. (31), the parameter q must be chosen appropriately. We found that for the paraxial case the truncation error is roughly minimized by choosing q = q _opt = √3/${\theta }_{\text{max}}^2$. This choice is also suitable for the nonparaxial case when θ _max = π/4 as can be seen from Fig. 5, which shows the behavior (on the right) of the truncation error ε (for l _max ranging from 1 to 30) obtained by using the multipolar expansion q = 0 (orange dots), and the expansion in the functions A_l (black dots) for q ranging from 0 to 10. On the left, the green curve represents the angular spectrum A(θ), the black curve is the inverse mapped spectrum A(θ) for the specified q, and the blue curve is A corresponding to q _opt. Note that, as expected, the truncation error ε is smallest when q≈q _opt. Let us now explore the effect of changing the numerical aperture, i.e., θ _max. As can be seen from Fig. 6a, the inverse-mapped spectrum A(θ) corresponding to q _opt (blue curve) is nearly independent of θ _max. Therefore, the expansion coefficients c_l and the truncation error ε (shown as black dots in Figs. 6b and c, respectively) are, for this choice of q, also roughly independent of the numerical aperture. For comparison, Figs. 6b,c also show (as orange dots) the expansion coefficients and truncation error for the multipolar basis (q = 0). Finally, Fig. 7 shows contour plots of the amplitude of the corresponding fields for θ _max ∊ [0, π/2]. These plots were evaluated by using q = q _opt and l _max = 15.

 

Fig. 5. a) Angular spectrum A(θ) (green curve), and inverse-mapped spectrum A(θ) for the specified q (black curve) and for q = q _opt (blue curve). b) Truncation error ε, for l _max ranging from 1 to 30, corresponding to the multipolar expansion (orange dots), and to the expansion in the functions A _l (black dots) for q running from 0 to 10. [Media 3]

Download Full Size | PDF

As in the paraxial case, the truncation errors for both the multipolar basis and the new scheme decay slowly with l _max due to the abrupt edge of the angular spectrum. In the next section, we present an example with no discontinuities.

5.3. Example: Aberrated field

Let us consider now a smooth “flat-top” angular spectrum given by the expression

for w ₀ = 2. This spectrum is shown as a green line in Fig. 8a. Also shown as a blue line is the corresponding inversely-mapped spectrum for q = 3. To make this example more interesting, we also include some amount of spherical aberration. The resulting angular spectrum is given by

 

Fig. 6. a) Angular spectrum A(θ) (green curve), and inverse-mapped spectrum A(θ) for q = q _opt (blue curve), for θ _max running from 0 to π/2. b) Expansion coefficients c_l for the multipolar basis (orange dots), and for the basis of functions A_l (black dots). c) Truncation error ε, for l _max ranging from 1 to 30, corresponding to the multipolar basis (orange dots), and to the basis of functions A_l (black dots). [Media 4]

Download Full Size | PDF

 

Fig. 7. Contour plots of the field amplitude (normalized to the focal amplitude) at the focal region of an apertured spherical wave with half-angle θ _max ∊ [0,π/2], calculated by using Eq. (31) with q = q _opt and l _max = 15. [Media 5]

Download Full Size | PDF

where the part of the phase proportional to θ ⁴ is responsible for the spherical aberration, and the part proportional to cosθ is a shift of magnitude 2s of the origin in the z-direction (i.e. a defocus) included to improve the convergence of the expansion. Figure 8b shows, for s = 15, the truncation error as a function of l _max for the new basis with q = 3 (black dots), as well as the corresponding error for the multipolar expansion (orange dots). Notice that, due to the smooth edges of the spectrum, the error decays much faster than for the previous example. Finally, Fig. 9 shows the field amplitudes at the focal region for several amounts of spherical aberration.

 

Fig. 8. a) Magnitude of the angular spectrum A(θ) in Eq. (32) (green curve), and inverse-mapped spectrum A (θ) for q = 3 (blue curve). b) Truncation errore, for l _max ranging from 1 to 30, corresponding to the multipolar expansion (orange dots), and to the expansion in the functions A_l (black dots) for q = 3, for the aberrated angular spectrum in Eq. (33) with s = 15.

Download Full Size | PDF

 

Fig. 9. Contour plots of the field amplitude (normalized to the unaberrated focal amplitude) at the focal region of the field corresponding to the angular spectrum in Eq. (33) for s between 0 and 30, calculated by using Eq. (31) with q = 3 and l _max = 15. [Media 6]

Download Full Size | PDF

6. Concluding remarks

By using a geometric mapping of the sphere of directions, a new orthonormal basis of free nonparaxial monochromatic fields was introduced. This basis allows any axis-symmetric field to be computed in a numerically efficient fashion. The new scheme can be regarded as a bridge between the multipolar expansion (suitable for highly nonparaxial fields) and a paraxial expansion in terms of Gaussian beams. The convenience of this scheme is illustrated by the numerical study in Sections 5.2 and 5.3. For comparison, the field corresponding to the angular spectrum in Eq. (33) with s = 15 was also computed as an integral over θ of weighted Bessel beams, yielding a plot that is indistinguishable from the one in Fig. 9 (and therefore not shown) but requiring over two orders of magnitude more computation time. The generalization of this scheme for the case of non-rotationally symmetric fields, both scalar and electromagnetic, will be considered in the future.

MAA acknowledges support from the National Science Foundation through the Career Award number PHY-0449708.