###### Abstract

We construct solutions of the paraxial and Helmholtz equations which are polynomials in their spatial variables. These are derived explicitly using the angular spectrum method and generating functions. Paraxial polynomials have the form of homogeneous Hermite and Laguerre polynomials in Cartesian and cylindrical coordinates respectively, analogous to heat polynomials for the diffusion equation. Nonparaxial polynomials are found by substituting monomials in the propagation variable with reverse Bessel polynomials. These explicit analytic forms give insight into the mathematical structure of paraxially and nonparaxially propagating beams, especially in regards to the divergence of nonparaxial analogs to familiar paraxial beams.

Paraxial and nonparaxial polynomial beams and the analytic approach to propagation

Mark R Dennis Jörg B Götte Robert P King Michael A Morgan and Miguel A Alonso

OCIS codes:

A major aim of optical physics is to find exact mathematical solutions of electromagnetic wave equations, which can be identified with observations of the electromagnetic field. We explore here a set of solutions of the scalar paraxial and nonparaxial (Helmholtz) equations arising naturally from Taylor series expansions, which are polynomials in and and (Cartesian coordinates), or and (cylindrical coordinates).

These polynomial beams arise out of Taylor expansions of optical fields in free space, and describe the optical amplitude behavior local to the focal plane and optical axis, although they diverge at infinity. They are therefore a powerful tool in understanding the local structure of optical fields, without concern for global properties (which is the domain of Fourier analysis). Examples of such analysis by polynomials and Taylor series includes the interaction of optical vortices [1, 2, 3, 4] and optical superoscillation [5]. Polynomial beams also illuminate well-known problems in constructing nonparaxial counterparts to simple paraxial beams such as Gaussians [6, 7, 8]. Paraxial polynomials, which are more straightforward than nonparaxial, have recently been considered in studies of this relationship [9, 10]. We will derive explicit forms of polynomial solutions to the paraxial equation and reduced Helmholtz equation ( solves the 3D Helmholtz equation when solves this latter equation). These forms give further insight into the difference between paraxial and exact propagation.

Optical fields with a known analytic form for some constant may readily be represented by power series in the spatial variables. The properties of beams based on the special functions of mathematical physics, such as Gaussian, Bessel and Airy beams, have analytic forms of this type. They are usually specified by an initial amplitude distribution when given by a power series in or with a natural inverse length associated with the initial field. is for a 2D plane wave, for a Bessel beam and for a Gaussian beam. Of course, power series expansions must converge with respect to their variables to make physical sense.

For instance, a paraxially propagating cylindrical Gaussian beam has the form with waist and Rayleigh distance In the waist plane this function is expanded as an exponential power series in However, for fixed the function can also be expanded in around and the appearance of in the denominator of the exponent implies that for there is an essential singularity at and the power series in does not converge beyond the Rayleigh distance. By comparison, similar series in of paraxial, propagation-invariant beams, such as plane waves and Bessel beams converge for all (infinite radius of convergence). The nonparaxial counterparts of many paraxial beams (i.e. with the same initial amplitude distribution), such as do not converge for any

We define paraxial polynomials to be solutions of the paraxial equation which are monomials in the initial plane, such as or In Cartesian coordinates, we write the paraxial polynomial with The Fourier transform (in ) of is the derivative -function so we write the paraxially propagating polynomial in terms of an angular spectrum integral:

(1) | |||||

(2) | |||||

(3) |

where we apply, in the second line, integration by parts times, and in the third line, the substitution The expression inside square brackets in (3) is the generating function of the Hermite polynomials[11], from which it immediately follows that

(4) |

that is, the paraxial polynomial which is when is a homogeneous polynomial of order with the coefficients of Hermite polynomials, each of whose terms is made up of powers of and to give the same dimension of [Length]. Since only includes odd or even powers of depending on whether is odd or even, only even powers of occur in the paraxial polynomial, which is then genuinely a polynomial in and (no fractional powers). Generalizing this argument, Cartesian paraxial polynomials with initial data are simply the products

Paraxial polynomials in cylindrical coordinates with vortices of order on axis (with ), are very easily derived by an identical argument to that above, with analytic initial data

(5) |

which is a homogeneous, associated Laguerre polynomial [11] in and Since by itself solves the paraxial equation (and also the reduced Helmholtz equation). The occurrence of Hermite and Laguerre polynomials here is somewhat analogous to the appearance of these polynomials with complex arguments in elegant Hermite-Gaussian and Laguerre-Gaussian beams [12].

Paraxial polynomials have been derived by different methods and studied before, for example in [3, 4, 9, 10]. The paraxial equation may be thought of as the heat equation with imaginary time and in this guise the paraxial polynomials are known as ‘heat polynomials’ [13]. Any analytic solution of the heat or paraxial equation, with a specified initial field distribution at expressed as a power series, can be constructed directly by substituting the appropriate power of and or with the appropriate paraxial polynomial, by linearity of the differential equations.

As described above, the amplitude distribution of any beam involves a power series in or with an inverse length; usually we think of this as an expansion in the spatial variable or about 0. However, when we replace or with the polynomials (4), (5), the propagated function is formally a power series expansion in about 0. In fact, this is how paraxial polynomials occur in free space paraxial beams: functions such as the Bessel beam or Gaussian beam expanded about must term-by-term satisfy the wave equation, since does not appear in the paraxial equation. Each term in must therefore be a paraxial polynomial of order thus, all analytic paraxial beams can be expressed as appropriate sums of paraxial polynomials, and paraxial beams are generating functions for paraxial polynomials. The nonparaxial series expansions considered in [6] are effectively in and we consider these briefly below.

We may also think of the propagation of paraxial polynomials as the appropriate initial monomial times a Gaussian of asymptotically large width (limiting to a plane wave). In this limit, approaches infinity, and the beam close to propagates like the polynomial. The Fourier transform of such a wide Gaussian limits to the derivative -function in the derivation (1) above. This approach was used to create fields with desired geometric properties – knotted optical vortices – using polynomial methods [4], which were then embedded analytically in more physical Gaussian beams without affecting their local knot topology.

Nonparaxial polynomials, which solve the reduced Helmholtz equation given above, may be constructed by modifying the angular spectrum method. The nonparaxial propagator, instead of involving in the angular spectrum integral (1) (or its cylindrical counterpart), involves to which the paraxial propagator is the Fresnel approximation. Now,

(6) |

can be thought of as a generating function for monomials by expanding in Substituting with therefore means that the powers in in the paraxial polynomials are replaced by the polynomials generated by

In fact, there are two related polynomial sequences, the reverse Bessel polynomials [14, 15]

(7) |

where is a modified Bessel function of half-integer order [11]. In the literature, is usually called the ‘reverse Bessel polynomial’, and clearly The ‘Bessel polynomials’ are an orthogonal polynomial sequence, but themselves are not. has the explicit form and for

(8) |

The reverse Bessel polynomials are given by the exponential generating function [14, 15]

(9) |

Thus, by elementary substitution, the nonparaxial exponential in the angular spectrum can be expanded

(10) |

Thus by (6) and (10), nonparaxial polynomials can be found from the corresponding paraxial polynomial by the substitution of all powers of by polynomials,

(11) |

This statement is the main result of this Letter.

Explicitly, the nonparaxial polynomials corresponding to (4) and (5) are

(12) | |||||

(13) |

from the appropriate forms of the Hermite and associated Laguerre polynomials, and the paraxial polynomials are recovered asymptotically as Nonparaxial polynomials can be generated by expanding explicit solutions of the reduced Helmholtz equation around This was done in [1, 2, 3] for the plane wave and the Bessel beam although the simple forms (12), (13) were not realized.

The replacement (11) is equivalent to the approach of Wünche [7], in which operators and in and were constructed to find the nonparaxial counterpart to a paraxial beam. It can be shown that so substitution (11) is equivalent to the operation of The action of gives a second family of nonparaxial polynomials. Since these polynomials are the same as (12) and (13) but with replacing Rather than satisfying the Dirichlet condition on the beam’s amplitude at this second family satisfies a Neumann-like condition when (Ref. [7] Eq. (3.30)). Similar operators involving Bessel polynomials were considered for the heat and paraxial equations in [16].

A general nonparaxial beam is therefore represented by a triple power series. On rearranging, this gives a series in , whose leading-order term is the paraxial beam, and later terms are higher-order corrections [6, 17]. The known failure of many such nonparaxial power series to converge for can be seen explicitly using nonparaxial polynomials. Each monomial in the paraxial series expansion becomes a reverse Bessel polynomial, whose coefficients increase factorially as powers of decrease. For instance, the coefficient for of a cylindrical beam with initial expansion on nonparaxial propagation on-axis (), from (13) and (8) is a series To converge, must approach 0 quickly enough to defeat this factorial growth, which is not the case for Gaussian and Airy beams, although it is for Bessel beams.

The analytic approach suggests another interpretation about this divergence of nonparaxial beams: the spectra of nonparaxial beams with divergent power series expansions contain evanescent components which are neither forward nor backwards propagating, as their wavevector component in is imaginary. Evanescent plane waves solving the reduced Helmholtz equation have the form with Such waves, however, cannot be represented by a power series expansion about as this is on the opposite side of the branch point singularity at This suggests that nonparaxial power series expansions never formally converge for beams whose spectra include evanescent waves, such as Gaussian [18] and Airy beams [19]. In such a situation, it is appropriate to asymptotically resum the divergent tail of the series in and such an approach has been proposed [8, 10]. Further work will reveal the relationship between nonparaxiality, evanescence, and divergent series, and we believe the polynomial solutions constructed here will be a useful tool in these investigations.

## References

- [1] J. F. Nye, J. Opt. Soc. Am. A 15, 1132–1138 (1998).
- [2] M. V. Berry, J. Mod. Opt. 45, 1845–1858 (1998).
- [3] M. V. Berry and M.R. Dennis, J. Phys. A 34, 8877–8888 (2001).
- [4] M. R. Dennis, R. P. King, B. Jack, K. O’Holleran and M. J. Padgett, Nature Phys. 6, 118–121 (2010).
- [5] Y. Aharonov, F. Colombo, I. Sabadini, D. C. Struppa and J. Tollaksen, J. Phys. A 44, 365304 (2011).
- [6] M. Lax, W. H. Louisell, W. B. McKnight, Phys. Rev. A 11, 1365–1370 (1975).
- [7] A. Wünsche, J. Opt. Soc. Am. A 9, 765–774 (1992).
- [8] R. Borghi and M. Santarsiero, Opt. Lett. 28, 774–776 (2003).
- [9] A. Torre, J. Opt. 13, 015701 (2011).
- [10] R. Borghi, F. Gori, G. Guattari and M. Santarsiero, Opt. Lett. 36, 963–965 (2011).
- [11] National Institute of Standards and Technology, Digital Library of Mathematical Functions (NIST 2010), http://dlmf.nist.gov/.
- [12] A. E. Siegman, J. Opt. Soc. Am. 63, 1093–1094 (1973).
- [13] D. V. Widder, The heat equation (Academic Press 1976).
- [14] E. Grosswald, Bessel polynomials (Springer 1979).
- [15] S. Roman, The umbral calculus (Dover 2005).
- [16] L. R. Bragg and J. W. Dettman, Rocky Mountain J. Math. 25, 887–917 (1995).
- [17] Q. Cao and X. Deng, J. Opt. Soc. Am. A 15, 1144–1148 (1998).
- [18] C. J. R. Sheppard, J. Opt. Soc. Am. A 18, 1579–1587 (2001).
- [19] A. V. Novitsky and D. V. Novitsky, Opt. Lett. 34, 3430–3432 (2009).