New basis for rotationally symmetric nonparaxial fields in terms of spherical waves with complex foci

A scheme for computing rotationally-symmetric nonparaxial monochromatic scalar fields is proposed, based on a new orthonormal basis of solutions of the Helmholtz equation given by combinations of spherical waves focused at imaginary points. These basis fields are found through a mapping of the angular spectra of the multipolar basis over the sphere of directions. The convergence of the basis can be optimized by an appropriate choice of the location of the imaginary focus. The new scheme is tested for the case of converging spherical waves of different numerical apertures, with and without aberrations. © 2006 Optical Society of America OCIS codes: (260.1960) Diffraction theory; (260.0260) Physical optics; (000.3860) Mathematical methods in physics. References and links 1. E. Heyman and L.B. Felsen, “Complex source pulse-beam fields,” J. Opt. Soc. Am. A 6, 806-817 (1989). 2. C.J.R. Sheppard and S. Saghafi, “Beam modes beyond the paraxial approximation: A scalar treatment,” Phys. Rev. A 57, 2971-2979 (1998). 3. C.J.R. Sheppard, “High-aperture beams,” J. Opt. Soc. Am. A 18, 1579-1587 (2001). 4. M.A. Alonso, R. Borghi, and M. Santarsiero, “Joint spatial-directional localization features of wave fields focused at a complex point,” J. Opt. Soc. Am. A 23, 933-939 (2006). 5. B. Tehan Landesman and H.H. Barrett, “Gaussian amplitude functions that are exact solutions to the scalar Helmholtz equation,” J. Opt. Soc. Am. A 5, 1610-1619 (1988). 6. Z. Ulanowski and I.K.Ludlow, “Scalar field of nonparaxial Gaussian beams,” Opt. Lett. 25, 1792-1794 (2000). 7. J.J. Wen and M.A. Breazeale, “A diffraction beam field expressed as the superposition of Gaussian beams,” J. Acoust. Soc. Am. 83, 1752-1756 (1988). 8. D. Ding and Y. Zhang, “Notes on the Gaussian beam expansion,” J. Acoust. Soc. Am. 116, 1401-1405 (2004). 9. Y. Li, “Light beams with flat-top profiles,” Opt. Lett. 27, 1007-1009 (2002). 10. J. D. Jackson, Classical Electrodynamics (Wiley , New York, 1998), Chapters 3 and 9. #68236 $15.00 USD Received 6 March 2006; revised 11 July 2006; accepted 12 July 2006 (C) 2006 OSA 24 July 2006 / Vol. 14, No. 15 / OPTICS EXPRESS 6894 11. V. Bargmann, “Irreducible unitary representations of the Lorentz group,” Ann. Math. 48, 568-640 (1947). 12. A. E. Siegman and E. A. Sziklas, “Mode calculations in unstable resonators with flowing saturable gain. 1: Hermite-Gaussian expansion,” Appl. Opt. 13, 2775-2792 (1974). 13. R. Borghi, F. Gori, and M. Santarsiero, “Optimization of Laguerre-Gauss truncated series,” Opt. Commun. 125, 197-203 (1996). 14. M. Born and E. Wolf, Principles of Optics, 7th ed. (Cambridge University Press, Cambridge, 1999), Sec. 13.2.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.

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.

Multipolar fields
The spherical harmonics of order l and index m are defined by [10] Y l,±m (u) = (±1) m 2l + 1 4π where θ and φ are the polar and azimuthal angle, respectively, and l (•) 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 4π 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:

Complex focus fields
As can be seen from Eq. (1), shifting the spatial coordinate by r 0 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 0 is proportional to exp(−iu • r 0 ) and the corresponding field is proportional to sin |r − r 0 |/|r − r 0 |.Suppose now that the focus is at the imaginary point r 0 = (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 (sin r)/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., U 0,0 (r; q) ≈ 8π q exp(iz) 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).

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 twodimensional 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 cos Θ(θ 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)| 2 is uniform, and the radial lines are equally spaced.As q increases, the radial lines drift towards the positive zdirection, where |A 0,0 (u; q)| 2 peaks.

z x q=0 z x q=1
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)| 2 .For increasing q, these radial lines bunch up around the positive z axis.
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.
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 (m) l (•) 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 2 l + 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 max denotes the truncation order and N 2 = π a 2 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 circ function of Log 10 (ε) 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/θ 2 max .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]

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.

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 A l,0 (u; q) = q (2l + 1) 2π sinh2q exp(q cos θ ) P l exp(2q cos θ ) − cosh2q sinh2q . ( Because P l is a polynomial, the expression in Eq. ( 24) is a sum of terms proportional to exp[(2n + 1) q cos θ ].Therefore, the corresponding fields take the form where U 0,0 is the closed form given in Eq. ( 11), and with coef[Pol(μ), μ n ] 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 orthonormal 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 c l (q) = 2π 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.

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 where Ω 0 = {θ ∈ [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/θ 2 max .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 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.
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.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.

Example: Aberrated field
Let us consider now a smooth "flat-top" angular spectrum given by the expression for w 0 = 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, 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).
we also include some amount of spherical aberration.The resulting angular spectrum is given by where the part of the phase proportional to θ 4 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  32) (green curve), and inversemapped spectrum A (θ ) for q = 3 (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 = 3, for the aberrated angular spectrum in Eq. ( 33) with s = 15. aberration.

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.

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

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

Fig. 7 .
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.

Fig. 8 .
Fig.8.a) Magnitude of the angular spectrum A(θ ) in Eq. (32) (green curve), and inversemapped spectrum A (θ ) for q = 3 (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 = 3, for the aberrated angular spectrum in Eq. (33) with s = 15.

Fig. 9 .
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.