Solutions to Maxwell's Equations using Spheroidal Coordinates

Analytical solutions to the wave equation in spheroidal coordinates in the short wavelength limit are considered. The asymptotic solutions for the radial function are significantly simplified, allowing scalar spheroidal wave functions to be defined in a form which is directly reminiscent of the Laguerre-Gaussian solutions to the paraxial wave equation in optics. Expressions for the Cartesian derivatives of the scalar spheroidal wave functions are derived, leading to a new set of vector solutions to Maxwell's equations. The results are an ideal starting point for calculations of corrections to the paraxial approximation.


Introduction
Solutions to the wave equation in spheroidal coordinates have application to a wide range of problems in physics [1]. In particular, in the short wavelength limit, consideration of the wave equation in oblate spheroidal coordinates leads directly to the well known Gauss-Laguerre solutions to the paraxial approximation of the wave equation in optics [2]. The paraxial approximation is an extremely versatile tool for studying beams of coherent radiation (e.g. laser beams), and successfully describes phenomena such as propagation of beams through lens systems, focus size, wave front curvature and phase shifts in the focus of beams, as well as eigenmodes and eigenfrequencies of spherical mirror resonators [3]. Despite its great success, the paraxial approximation is only an approximation. Consideration of corrections beyond the paraxial approximation is relevant not only in order to establish a bound on its validity but is in fact necessary to describe strongly focused beams [4] and to explain the fine structure of the frequency spectrum of high finesse resonators [5].
Corrections to the paraxial approximation can be calculated by reconsidering the term which is neglected when deriving the paraxial approximation from the wave equation [6]. As an alternative, corrections can be effectively considered using exact solutions to the wave equation in spheroidal coordinates. This has the advantage of removing ambiguities in the definition of higher-order terms, leads to simpler expressions, and is in general a more natural framework for exact consideration of beamlike solutions to the wave equation.
The wave equation is separable in spheroidal coordinates, allowing solutions to be written as the product of a so called radial function depending only on the ξ coordinate, a so called angle function depending only on the η coordinate and the function e imφ depending only on the azimuthal φ coordinate [1]. For short wavelengths, solutions for the radial and angle functions in the form of asymptotic expansions have been known for a long time [7,8,9]. The close relationship between spheroidal coordinates and Laguerre-Gaussian beams suggests that the asymptotic solutions are in some way related to the Gauss-Laguerre solutions to the paraxial wave equation. While this relationship is directly obvious for the asymptotic solutions for the angle functions, the opposite is true for the previously known solutions for the radial functions. In fact, the previously known asymptotic solutions for the radial functions consist of multiple sums with a huge number of terms.
Using solutions to the wave equation for exact calculations in optics requires taking the transverse vector character of the electromagnetic fields into account. Full vector solutions to Maxwell's equations were already considered by Ref. [10]. These vector solutions are derived by applying vector operators to the scalar spheroidal solutions and are expressed as derivatives in the spheroidal coordinates of the scalar solutions. As for the asymptotic solutions for the radial functions, the resulting expressions are quite complicated.
In this paper, we derive a new expression for the asymptotic expansion of the radial functions which to lowest order is directly reminiscent of solutions of the paraxial approximation. The new expression is significantly simpler, making its analytic application practicable. In addition, we derive expressions for the Cartesian derivatives of the spheroidal wave functions which allow us to define a significantly simpler set of vector spheroidal wave functions. The component of the field transverse to the direction of propagation is equal to a single scalar spheroidal wave function. This creates a direct correspondence between our vector solutions and Gauss-Laguerre solutions to the paraxial approximation.
After a brief discussion on notation used in this paper in section 1.1, we introduce oblate spheroidal coordinates in section 2. The wave equation in spheroidal coordinates as well as separation of coordinates is discussed in section 2.1. Section 3 is a review of the asymptotic expansion of the angle functions. The derivation is equivalent to the one given e.g. in Ref. [1]. The radial functions are considered in section 4. Starting with the previously known asymptotic expansion allows us to derive a new expression for the two lowest-order terms. Taking the lowest-order term as an ansatz for the radial function leads to a differential equation which allows higher-order terms to be calculated iteratively. In section 5 we define scalar spheroidal wave functions and demonstrate their equivalence in the short wavelength limit to the Gauss-Laguerre solutions to the paraxial wave equation. This is followed by a brief discussion concerning the symmetry properties of the scalar spheroidal wave functions upon reflection through a plane containing the symmetry axis of the spheroidal coordinate system. We proceed in section 6 by deriving expressions for the Cartesian derivatives of the scalar spheroidal wave functions. These expressions allow us to define a new set of vector spheroidal wave functions in section 7 which are transverse vector solutions to the wave equation. section 7 concludes with expressions for the curl of the vector spheroidal wave functions.

Notation
The notation used in this paper is generally consistent with the notation in Ref. [1], with the exceptions noted here. Previously, the symbol c has been used to quantify the scaling of the spheroidal coordinate system relative to the wavelength. Since it seems problematic to use the same symbol as for the speed of light in a theory with strong application to optics, we usec instead.
Additionally, Ref. [1] generally uses m and n to label the spheroidal functions and associated parameters, using the labels m and ν with (m, n) = (m, 2ν + m) only when Laguerre functions are directly involved. In the short wavelength limit, both sets of labels can be used. The advantage of the m, ν labeling is that it is identical to the labeling of the Laguerre polynomials, which are essential for the expansion of the angle functions in the short wavelength limit. As a result, we use the labels m and ν exclusively.
Outside the short wavelength limit, not considered in this paper, use of the labels m and n is necessary. This is due to the fact that when the short wavelength limit does not apply, the labels m, ν do not uniquely label the set of solutions to the wave equation whereas the labels m, n do. Each pair of values m, ν corresponds to two solutions for the angle function and two solutions for the radial function with opposite parity about the plane ξ = 0. In the short wavelength limit, the two distinct angle functions converge to the same function asymptotically, and the two radial functions can be considered as the real and imaginary part of a single complex valued function, so that m, ν effectively labels all solutions uniquely.
In contrast to Ref. [1], we attempt to consistently place the labels m and ν as subscripts, generally placing other labels as superscripts. The only exception is the Laguerre functions, where we use the standard notation L (m) ν . Since we are only interested in oblate spheroidal coordinates and real arguments η and ξ, we write the angle and radial functions respectively as S mν (η) and R mν (ξ) instead of as S mn (−ic, η) and R mn (−ic, iξ) as in Ref. [1], thereby ignoring the relationship between solutions in oblate spheroidal coordinates and in prolate spheroidal coordinates. Finally, the indices for the sums over Laguerre functions as e.g. in Eq. (12) range from 0 to ∞ rather than from −ν to ∞ as in Ref. [1]. This allows us to treat all solutions for the angle functions on an equal footing, without choosing a specific one in advance by fixing ν ahead of time. The two-dimensional elliptic coordinate system is defined from the set of all ellipses and all hyperbolas with a common set of two focal points [1]. We denote the separation of the two focal points by d. Oblate spheroidal coordinates are derived from elliptic coordinates by rotating the elliptical coordinate system about the perpendicular bisector of the focal points [1]. The focal points thereby sweep out a circle of radius d/2. Oblate spheroidal coordinates are commonly mapped to Cartesian (x, y, z) coordinates by placing this circle in the x − y plane with center at the origin. The coordinates are often labeled η, ξ and φ with the transformation to Cartesian coordinates given by [1]

Spheroidal coordinates
Oblate spheroidal coordinates reduced to the x − z plane are shown in Fig. (1). The coordinates r and θ defined by are useful to describe aspects of the coordinate system. Note that near the z-axis, ξ quantifies distance along the z-axis and θ quantifies distance from the z-axis. To lowest order in θ, z = d 2 ξ and r = d 2 1 + ξ 2 θ.

Wave equation in spheroidal coordinates
The wave equation, written in oblate spheroidal coordinates, is given by [1] ∂ ∂ξ This equation is separable. Therefore ψ can be written as a product of three functions depending only on η, ξ and φ, respectively. Due to cylindrical symmetry, the function depending on φ can be written as e imφ for integer values of m, The label ν is used to denote the various solutions of Eq. (4). The so called radial function R mν (ξ) and angle function S mν (η) satisfy the two differential equations as can be seen from comparison with Eq. (4). Here, λ mν is the separation constant.

Angle functions for largec
In the limitc → ∞, a solution to Eq. (7) expressed as a series in 1 c can be found as follows. The ansatz allows Eq. (7) to be rewritten as We have dropped the index ν in order to refer to a general solution with arbitrary λ m(ν) . Note that by introducing the factor (1−η 2 ) m 2 in Eq. (8), we break the symmetry between positive and negative m. This symmetry nonetheless persists in the original Eq. (7) so that the solutions must in the end be symmetric under a transformation from m to −m. We address this issue in more detail in section 5.
Forc → ∞, the terms on the second line of Eq. (9) inside the bracket vanish. What is left is simply the associated Laguerre differential equation. A possible set of approximate solutions to Eq. (9) is therefore the set of associated Laguerre polynomials L (m) and λ m given by The form of Eq. (10) suggests expanding the exact solution s m (x) in terms of associated Laguerre polynomials as The A s m are expansion coefficients. Inserting this expression into Eq. (9) and using appropriate recursion relations for Laguerre polynomials allows the left hand side of Eq. (9) to be rewritten as a linear combination of Laguerre polynomials with fixed m. Note that a number of useful recursion relations for Laguerre Polynomials are listed in Appendix A. Since the Laguerre polynomials with fixed m are linearly independent, the coefficient of each term in this sum must be zero, leading to the following relations between the A s m , These relations can be written in matrix form as Here, M 0 is a diagonal matrix, M 1 is a tridiagonal matrix, I is the identity matrix and A is the vector of coefficients A s m according to We use the Kronecker delta function, δ s,t . Eq. (14) transforms the problem of determining the coefficients in Eq. (12) into an eigenvalue problem. λ ′ m is an eigenvalue ofcM 0 + M 1 and A is the corresponding eigenvector.
In the limitc → ∞, M 1 can be neglected compared tocM 0 . Since M 0 is diagonal, finding an approximate solution to Eq. (14) is therefore trivial, this solution being given by Eqs. (10) and (11). Improved approximate solutions can be found via perturbation theory [11]. According to convention we reintroduce the label ν to denote the solution whose eigenvector is given to lowest order in 1 c by A s mν = δ s,ν . We expand λ ′ mν and A s mν in powers of 1 c as Inserting these expressions into Eq. (14) and matching terms of equal order in 1 c leads to recursive expressions for β k mν and A s,k mν , Since M 1 is tridiagonal, A s,k mν is zero for |s − ν| > k. As a result, each sum over s or t contains only a finite number of nonzero terms.
Eqs. (17) lead to the expressions for β k mν and A s,k mν for k ≤ 4 given explicitly in Ref. [1]. Higher-order values can easily be obtained numerically from Eqs. (17). Note that the expression for β mn 4 in Ref. [1] contains an incorrect "−" sign.

Radial functions for largec
Following a similar procedure as the one in section 3, it is also possible to find an asymptotic expansion for the radial function, given by [1] Note that Ref.
[1] defines a total of four different radial functions R For integer values of m, the denominator and numerator in Eq. (19) both equal zero, in which case U (m) ν (z) is defined as a limit in m. For large |z|, an asymptotic expansion for U (m) ν (z) exists, given by [1,12] Unlike the asymptotic expansion for the angle function, which is quite transparent since it contains a single Laguerre polynomial to lowest order and additional Laguerre polynomials as higher-order corrections, the asymptotic expansion for the radial function is rather opaque since the terms s = 0 through s = ν all contain terms of equal and lowest order in 1 c as can be seen from Eqs. (18) and (20). We attempt to rewrite Eq. (18) as a sum over terms of equal order in 1 c . To this end we insert the asymptotic expansions for U (m) ν (z) and A s mν , Eqs. (20) and (16), into Eq. (18). This leads to sums over r, s and k. We define p = r + s + k − ν and observe that the total order in 1 c of a term in the triple sum is p + 1. Replacing the sum over k by a sum over p, one obtains The symbol ⌊.⌋ denotes rounding down to the nearest integer. As we will demonstrate, the sums over r and s can be performed resulting in relatively compact analytic expressions for all values of p. We start by considering the two lowest-order cases, p = 0 and p = 1. To proceed, we need analytic expressions for A s,k mν for s + k ≤ ν + 1. Using Eqs. (17), it can be shown by induction that and For p = 0, one obtains The first line is Eq. (21) with fixed p = 0. To obtain the second line, define z = 1 + iξ and observe that the first line is essentially the binomial expansion of (1 − z 2 ) ν . The case p = 1 is already slightly more complicated, The first line is again Eq. (21) but with fixed p = 1. The second line is obtained in a similar manner as before, by observing that for r = 0 the first line is essentially the binomial expansion of z 3ν+m+1 d dz z −(4ν+2m) 1 − z 2 ν−1 , and for r = 1 the first line is Using Eqs. (24) and (25), we obtain a preliminary result for the summation of terms of equal order in 1 c for the asymptotic expansion of R mν (ξ) through first order in 1 c . In writing the preliminary result, we must confront the issue that the higherorder terms are to some degree arbitrary. R mν (ξ) is defined by its differential equation only up to some constant factor. In particular, we can multiply R mν (ξ) by a constant which depends onc. The asymptotic expansion in 1 c of such a product contains different higher-order terms. Specifically, a constant a times the zeroth order term can be added to the p th order term as long as a times the nth order term is simultaneously added to the (p + n)th order term, resulting in a different asymptotic expansion. This issue is raised by the presence of A 0 mν as part of the normalization constant in Eq. (18), which depends on 1 c through first order according to Replacing A 0 mν according to Eq. (26) leads to This is a new and significantly simpler expression for the asymptotic expansion of the radial function through first order in 1 c . Continuing the previous procedure for larger values of p becomes increasingly cumbersome, not least because analytic expressions for A s,k mν for larger values of s + k − ν become increasingly difficult to find. The similar form of the results for p = 0 and p = 1 suggests an alternative approach. We factor out the ξ dependent terms in front of the bracket in Eq. (27), Inserting this expression into the differential equation for R mν (ξ) leads to a differential equation for r mν (ξ), We can write r mν (ξ) as an asymptotic expansion in 1 c , Eq. (29) can be transformed into a recursive expression for the r The β k mν are the expansion coefficients of λ mν according to Eq. (16). Through p = 2 one obtains The r
We briefly demonstrate the equality to lowest order in 1 c of Eq. (33) to the Gauss-Laguerre solutions to the paraxial wave equation. According to Eq. (2) we have 1 − η ∼ θ 2 /2. Due to the presence of the term e −c(1−η) , ψ mν takes on finite values only for θ = O(c −1/2 ). We can therefore make the following substitutions in Eq. (33).
This is the equation for a Gauss-Laguerre beam with confocal parameter d [3]. While the differential equations for the spheroidal functions are independent of the sign of m, the approach used to find the asymptotic series for the spheroidal functions introduces a dependence on the sign of m. As a result, it is by no means obvious that there is some symmetry between positive and negative m in Eq. (33). Demonstrating this symmetry is in fact nontrivial. The fact that the Laguerre polynomials transform from positive to negative m according to [13] L (−m) Surprisingly the transformation from positive to negative m involves a factor which depends onc.

From scalar to vector wave functions
The scalar spheroidal wave functions defined by Eq. (33) are solutions to the wave equation. Nonetheless, they are not suitable for exact calculations in optics, since solutions to Maxwell's equations are vector functions. Let E be the electric field of a solution to Maxwell's equations. In order to specify E, it is useful to choose a vector basis and to specify the components of E in this basis. Since our wave functions are defined in spheroidal coordinates, it seems reasonable to express E in terms of the spheroidal unit vectors, E = E ξêξ + E ηêη + E φêφ . Doing so has the major disadvantage that the individual components E ξ , E η and E φ are not solutions to the wave equation and therefore have no obvious relationship to the scalar spheroidal wave functions. We avoid this problem by using Cartesian unit vectors and write E = E xx + E yŷ + E zẑ . A sufficient condition for E to be the electric field of a solution to Maxwell's equations in free space is that each of the components E x , E y and E z satisfies the scalar wave equation and that ∇ · E = 0 [14]. As a result, we can construct solutions for E from our scalar spheroidal wave functions by writing each of E x , E y and E z as a linear combination of the ψ mν such that ∇ · E = 0 is satisfied. Since ∇ · E = ∂Ex ∂x + ∂Ey ∂y + ∂Ez ∂z , satisifying ∇ · E = 0 requires knowledge of ∂ψmν ∂x , ∂ψmν ∂y and ∂ψmν ∂z . Since Cartesian derivatives commute with the wave equation, ∂ψmν ∂x , ∂ψmν ∂y and ∂ψmν ∂z are again solutions to the wave equation, and can therefore be written as linear combinations of the ψ mν . We try to determine the coefficients of these linear combinations.
We begin by deriving expressions for ∂ ∂x , ∂ ∂y and ∂ ∂z in spheroidal coordinates. The scale factors in spheroidal coordinates, defined by h u i = ∂x ∂u i with u i ∈ {ξ, η, φ}, are [1] h ξ = d 2 The spheroidal unit vectors, given byê u i = 1 hu i ∂x ∂u i and expressed in terms of their Cartesian components, arê The gradient of a function expressed in spheroidal coordinates is Obtaining expressions for ∂ ∂x , ∂ ∂y and ∂ ∂z is simply a matter of combining Eqs. (38) and (39), and From here on, one can obtain expressions for the Cartesian derivatives of the ψ mν essentially by brute force. Eqs. (40) and (41) are applied to Eq. (33). The result is compared to successive orders in 1 c to Eq. (33). One obtains, We use J as the total angular momentum, equal to the sum of the orbital angular momentum and the spin angular momentum. So far we have ignored the fact that our solutions to the wave equation are traveling waves which can travel in both the +ξ and −ξ direction. All solutions so far have consistently represented waves traveling in the +ξ direction when an e −iωt time dependence is assumed. This is acknowledged in Eqs. (45) and (46) by the + superscript label in E + Jσν . We can also define vector spheroidal wave functions E − Jσν traveling in the −ξ direction. These are given by the same expressions, Eqs. (45) and (46), except thatẑ is replaced by −ẑ and the scalar spheroidal wave functions are evaluated at −ξ instead of ξ.
For many problems in physics, knowledge of both the electric as well as the magnetic field of an electromagnetic beam is needed [14]. Obtaining the magnetic field from the electric field requires taking the curl of the electric field and vice versa. The curl of the vector spheroidal wave functions can be obtained rather elegantly using Eqs. (42-44) to evaluate the necessary derivatives. We present the result here as an example for the usefulness of these equations.

Outlook
The results obtained here are an excellent starting point for calculating corrections to the paraxial approximation. Due to the equivalence to lowest order in 1 c of the Laguerre-Gaussian solutions to the paraxial approximation and the spheroidal wave functions considered here, solutions to physical problems based on the paraxial approximation are effectively lowest-order solutions in 1 c using spheroidal coordinates. Higher-order terms in 1 c can then be included as perturbations to obtain results to any desired degree of accuracy, providedc is large enough that the resulting series converge. This procedure can be used in particular to calculate corrections to the eigenfrequencies and eigenmodes of a Fabry-Perot resonator. The resulting corrections can be observed experimentally in a high-finesse resonator [5].
A similar procedure as the one used here to write the radial functions in a new form could possibly be applied for solutions to the wave equation in other coordinate systems, in particular for elliptic coordinates. Such solutions should reduce to the Hermite-Gaussian solutions to the paraxial wave equation in the short wavelength limit. Applying such solutions to calculate eigenfrequencies and eigenmodes of two-dimensional resonators would allow comparison with results obtained previously for two-dimensional resonators using other methods [15,16,17].