In many applications it is natural to seek to extract a characteristic scale for a function’s variations by reference to a frequency spectrum. Although the moments of a spectrum appear to promise simple options to make such a connection, standard Fourier methods fail to yield finite moments when the function’s domain is itself finite. We investigate a family of Fourier-like bases with rapidly decaying spectra that yield well-defined moments for such cases. These bases are derived by considering classes of functions for which a normalised mean square derivative is stationary. They are shown to provide precisely the type of spectrum needed to complete a recent investigation of mid-spatial frequency structure on optical surfaces [K. Liang, Opt. Express 27, 3390-3408 (2019)].
1. Introduction
Fourier methods have a well-earned status in the toolbox of science and technology. The standard methods include Fourier transforms –both continuous and discrete– as well as Fourier series. The associated spectra, typically in terms of spatial or temporal frequencies, have widespread applications. Oftentimes, however, generalisations or variations of the standard methods give a vital boost to utility in certain contexts. These variations include such things as time-frequency analysis [1], wavelets [2], and linear canonical transformations [3]. Further, since the Fourier series is just one special case of Sturm-Liouville methods [4], other instances are close cousins with similar orthogonality-related properties and potential advantages.
It turns out that Fourier methods often encounter significant challenges in situations where a function is known over only a finite domain. Standard methods can sometimes be applied in these cases by taking arbitrary steps such as replicating the function to make it complete and periodic or by assuming the function to be zero outside the finite domain (for the Fourier series or transform, respectively). In both these cases, however, even when the finite domain can be replicated to effectively tesselate a complete function that is periodic in each independent variable, it is difficult to avoid spectral artefacts associated with the edges of the original domain. The methods described below can be considered as potentially useful tools in these frequently encountered situations.
A relevant context for this discussion is the study of the surface shape of as-built optical surfaces with complex nominal shapes including aspheric and freeform optics [5–7]. The tool marks left by the associated fabrication steps are generally referred to as mid-spatial frequency structure (MSF). In order to (i) characterise, (ii) model the optical impact of, and (iii) set tolerances for MSF, it is natural to seek the benefits that typically come from an intuitive spectral decomposition. For most of these purposes, an orthogonal polynomial basis has been shown to be well matched [8]. The second of these steps, however, is typically performed by using a perturbation approximation and we recently showed that the accuracy of that approximation can be estimated in terms of the mean squared Laplacian of the MSF [9]. In a simple one-dimensional model, it is the mean squared second derivative that is required, of course. Although we claimed that such entities could be related to a fourth moment of the associated Fourier spectrum, we did not explicitly show how such a spectrum can be computed in general. In fact, we did not appreciate the subtleties involved and, once they were understood, came to the view that the required methods and ideas are valuable in themselves and potentially of more general interest.
To have a finite fourth moment, a spectrum must be sufficiently well localised. Figure 1(a) shows for an analytic function that orthogonal (Legendre) polynomials yield spectra that decay exponentially [10], and hence possess finite spectral moments of any order. On the other hand, the standard Fourier basis decays slowly and is not localised enough for the calculation of a fourth – or even second – spectral moment. It is shown below that a class of rapidly decaying Fourier-like (RDF) decompositions has increasingly compact spectra, as demonstrated in Fig. 1, and can provide a useful connection between these two limits. In this way, we finally resolve the nature of the spectra referred to in [9] that gives access to an intuitive characteristic scale of a function’s variations. Interestingly, some of these RDF bases are useful in a range of other applications [11–15].
2. RDF bases in one dimension
The derivation of the basis elements with the desired properties proceeds by requiring a certain normalized mean square derivative over the finite domain $a = [-L/2, L/2]$ to take a stationary value. In particular, we look for real functions $\textrm {V}(x)$ such that the ratio
(1)$$\frac{\int_a [ \textrm{V}^{(N)}(x) ]^{2}\, {\textrm{d}}x}{\int_a \textrm{V}^{2}(x)\, {\textrm{d}}x },$$
is stationary, where $N \ge 1$ is a positive integer and $\textrm {V}^{(N)}(x)$ denotes the $N$-th derivative. This can be achieved by introducing $\alpha$, a Lagrange multiplier, in order to construct the following functional: (2)$$\mathcal{L}(\textrm{V},\alpha) = \int_a [ \textrm{V}^{(N)}(x) ]^{2}\, {\textrm{d}}x - \alpha^{2N} \int_a \textrm{V}^{2}(x)\, {\textrm{d}}x.$$
Upon slightly perturbing $\textrm {V}$, namely $\textrm {V}(x) \rightarrow \textrm {V}(x) + \varepsilon (x)$, there should be no change to $\mathcal {L}$ at first order [or equivalently to the value of the expression in Eq. (1)]. By making this substitution and expanding, we find (3)$$\mathcal{L}(\textrm{V}+\varepsilon,\alpha) = \! \mathcal{L}(\textrm{V},\alpha) \!+ \! 2 \!\left[ \int_a \! \textrm{V}^{(N)}(x) \varepsilon^{(N)}(x) \, {\textrm{d}}x - \alpha^{2N} \!\!\int_a \!\textrm{V}(x) \varepsilon(x) \, {\textrm{d}}x \right] \! + \! O(\varepsilon^{2}).$$
Stationarity follows upon choosing $\textrm {V}$ so that the quantity within square brackets in Eq. (3) vanishes for any $\varepsilon (x)$. It is possible to rewrite this term by performing integration by parts $N$ times on the first integral in order to shift all of the derivatives from $\varepsilon$ to $\textrm {V}$. The resulting expression contains terms evaluated at the endpoints plus an integral in $x$ whose integrand has a global factor of $\varepsilon (x)$. Requiring the result to vanish for arbitrary $\varepsilon (x)$ leads to the following $2N$-th order Sturm-Liouville differential (eigenvalue) system: (4)$$\textrm{V}^{(2N)}(x) = (-\alpha^{2})^{N} \textrm{V}(x), \quad \textrm{V}^{(j)} \!\left( - \frac{L}{2} \right) = \textrm{V}^{(j)} \!\left( \frac{L}{2} \right) = 0 \quad \textrm{for} \quad j = N, \dots, 2N-1.$$
The solutions to Eq. (4) naturally split into even and odd forms given by (5a)$$ {{}^{N} \!} \textrm{V}_{\textrm{e}, i}(x) = {{}^{N} \!}\Gamma_{\textrm{e}, i} \sum_{\ell=0}^{N-1} a_\ell \cos\left[ ({-}1)^{\ell/N} \alpha_{\textrm{e}, i} x \right] ,$$
(5b)$$ {{}^{N} \!}\textrm{V}_{\textrm{o}, i}(x) = {{}^{N} \!}\Gamma_{\textrm{o}, i} \sum_{\ell=0}^{N-1} b_\ell \sin \left[ ({-}1)^{\ell/N} \alpha_{\textrm{o}, i} x \right] ,$$
where (without loss of generality) $a_0 = b_0 = 1$, and the remaining values of $a_\ell , b_\ell$ as well as the eigenvalues $\alpha _{\textrm {e}, i}$ and $\alpha _{\textrm {o}, i}$ are determined by the boundary conditions in Eq. (4). Finally, ${{}^{N} \!}\Gamma _{\textrm {e}, i}$ and ${{}^{N} \!}\Gamma _{\textrm {o}, i}$ are normalization constants chosen such that (6)$$\int_a \left[ {{}^{N} \!}\textrm{V}_{\textrm{e}, i}(x) \right]^{2} \, {\textrm{d}}x =1 \quad \textrm{and} \quad \int_a \left[ {{}^{N} \!}\textrm{V}_{\textrm{o}, i}(x) \right]^{2} \, {\textrm{d}}x =1.$$
Although it is not immediately obvious, it turns out that the RDF elements in Eqs. (5) are real-valued. Because this is a Sturm-Liouville problem, it is guaranteed that there are an infinite number of eigenvalues (found to be roughly proportional to $i$).In addition to the solutions in Eqs. (5), there are $N$ solutions (degenerate for $N>1$) with eigenvalue equal to zero that, in keeping with orthonormality, can be taken to be
(7)$${{}^{N} \!} \textrm{V}_{\textrm{d}, i}(x) = \sqrt{\frac{i + 1/2}{L/2}} P_i \! \left( \frac{x}{L/2} \right), \quad \textrm{for} \quad i = 0, \dots, N-1,$$
where the Legendre polynomial $P_i(x)$ is taken to be defined over $[-1,1]$. Note that these degenerate solutions are also alternately even and odd and automatically satisfy the boundary conditions of Eq. (4). In the limit $N \rightarrow \infty$, only this degenerate polynomial subspace remains, with eigenvalue $\alpha = 0$. This observation reveals an explicit connection of this family of bases between polynomials ($N\rightarrow \infty$) and the Fourier basis (similar to the case of $N=1$). For intermediate values of $N$, the non-degenerate subspace of solutions given in Eqs. (5) has properties that transition between orthogonal polynomials and Fourier series.Since Eq. (4) is a Sturm-Liouville differential equation, it is well-known that the values of $\alpha _{\textrm {e}, i}$ and $\alpha _{\textrm {o}, i}$ alternate in ascending order. Hence, it is straightforward to assign a global ordering to RDF elements simply through the order of their associated eigenvalues. The degenerate orthogonal polynomials (which all have zero eigenvalue) can be ordered as the first $N$ elements through their polynomial degree. Once ordered, the basis elements can simply be referred to through ${{}^{N} \!} \textrm {V}_n(x)$ with eigenvalue $\alpha _n$, where $n$ corresponds to the number of zeros. Figure 2 shows some of these functions for $N=1,2,3$, and $4$, along with the standard Fourier basis and the Legendre polynomials ($N\rightarrow \infty$) for comparison. It is evident that larger $N$ values give rise to non-degenerate basis elements that are steeper near the boundaries of the region.
RDF functions are orthogonal under two different inner products:
(8a)$$ \int_a {{}^{N} \!} \textrm{V}_n (x) \, {{}^{N} \!} \textrm{V}_{n'}(x) \, {\textrm{d}}x = \delta_{nn'},$$
(8b)$$ \int_a {{}^{N} \!}\textrm{V}^{(N)}_n (x) \, {{}^{N} \!}\textrm{V}^{(N)}_{n'}(x)\, {\textrm{d}}x = \alpha_n^{2N} \, \delta_{nn'},$$
While the first of these relations is a natural result of the fact that these are the eigenfunctions of a Hermitian operator, namely $({\textrm {d}}/{\textrm {d}}x)^{2N}$ along with Eq. (4), the proof of the second requires $N$ iterations of integration by parts as well as Eq. (4). It follows from the first of these orthogonality relations that the RDF expansion of any function $\mu (x)$ (e.g. the height of a surface with MSF structure) can be written as (9)$$\mu(x) = \sum_{n = 0}^{\infty} {{}^{N} \!}c_n \, {{}^{N} \!} \textrm{V}_n (x), \quad \textrm{where} \quad {{}^{N} \!}c_n = \int_{{-}L/2}^{L/2} \mu(x) \,\,{{}^{N} \!} \textrm{V}_n(x)\, {\textrm{d}}x.$$
Additionally, the squared norm of $\mu (x)$, referred to here as $M$, can be expressed as the sum of the square of the coefficients following Parseval’s relation as (10)$$M \triangleq \langle \mu^{2} (x) \rangle = \sum_n \left({{}^{N} \!}c_n \right)^{2}.$$
By using Eq. (4) and integration by parts it is possible to rewrite the integral in Eq. (9) as (11)$${{}^{N} \!}c_{n} = \frac{1}{\alpha^{2N}_{n}} \left[ \mu^{(N)}\, {{}^{N} \!}\textrm{V}_n^{(N-1)}{\bigg|}_a - \int_a \mu^{(N+1)}(x) \, {{}^{N} \!} \textrm{V}_n^{(N-1)}(x) \, {\textrm{d}}x \right].$$
As can be seen from Eqs. (5), the first term inside the square brackets of Eq. (11) has a leading behavior that is proportional to $\alpha _n^{N-1}$. Therefore, in the limit of large $n$, the spectral coefficients $|{{}^{N} \!}c_n |^{2}$ generally decay as $\alpha _n^{-(2N+2)}$, so the lines in the log-log plot in Fig. 1(a) approach a slope of $-(2N+2)$ for large $n$. Note that continuing to integrate Eq. (11) by parts leads to a closed-form expression for the RDF expansion whenever $\mu (x)$ is a polynomial.The second orthogonality relation, Eq. (8b), gives access to a simple calculation for the normalized mean square of the $N$-th derivative of $\mu (x)$ akin to the form appearing in Eq. (1):
(12)$$\begin{aligned} \frac{1}{M}\left\langle \left[ \frac{\partial^{N}}{\partial x^{N}} \mu(x) \right]^{2} \right\rangle &= \frac{1}{M}\int_a \left[ \sum_{n} {{}^{N} \!}c_{n} {{}^{N} \!}\textrm{V}_{n}^{(N)} (x) \right]^{2} \, {\textrm{d}}x\\ &= \frac{1}{M} \sum_{n} \sum_{n'}{{}^{N} \!}c_{n}\, {{}^{N} \!}c_{n'}\int_a {{}^{N} \!}\textrm{V}_{n}^{(N)}(x) \, {{}^{N} \!}\textrm{V}_{n'}^{(N)}(x) \, {\textrm{d}}x,\\ &= \frac{1}{M}\sum_{n} \left({{}^{N} \!}c_{n}\right)^{2} \alpha_{n}^{2N}, \end{aligned}$$
where Eq. (8b) was used in the last step and M is defined in Eq. (10). This equation reveals that the mean square $N$-th derivative corresponds simply to a normalized spectral moment of order $2N$. Furthermore, the $n$-th RDF element has $n$ zeros across this region; this provides an intuitive understanding as to why the RDF spectra shown in Fig. 1(a) peak around the seventeenth basis element (the sinusoid in the inset of Fig. 1 has 17 zeros). In Fig. 1(b), which shows the same coefficients plotted against their eigenvalues, the peak is approximately at $2\pi \,8.8\approx 55.3$.To give a clearer appreciation, the RDF basis for each of the cases of $N = 1$ and $N = 2$ is derived in the following subsections.
2.1 $N = 1$
For $N=1$, Eq. (4) can be solved to give even and odd orthonormal solutions of the form
(13a)$$ {}^{1}\! \textrm{V}_{\textrm{e}, i} (x) = \sqrt{\frac{2}{L}} \cos ( \alpha_{\textrm{e}, i} x) ,$$
(13b)$$ {}^{1}\!\textrm{V}_{\textrm{o}, i}(x) = \sqrt{\frac{2}{L}} \sin (\alpha_{\textrm{o}, i} x) ,$$
where the allowed values of $\alpha _{\textrm {e}, i}$ and $\alpha _{\textrm {o}, i}$ are given by (14)$$\alpha_{\textrm{e}, i} = \frac{2\pi }{L} i \quad \textrm{and} \quad \alpha_{\textrm{o}, i} = \frac{2\pi}{L} \left(i+\frac{1}{2} \right),$$
for $i \ge 1$. The zero eigenvalue has the solution ${}^{1} V_{\textrm {d},0}=1/\sqrt {L}$. After ordering the even, odd, and degenerate solutions in the manner described earlier, the eigenvalues can be written as $\alpha _n = \pi n /L$, with $n \ge 0$. Note that the solutions given by Eqs. (13) satisfy Neumann boundary conditions, whereas the standard Fourier basis functions satisfy periodic boundary conditions (see Fig. 2).2.2 $N=2$
For $N = 2$, Eq. (4) can be solved to give orthonormal, non-degenerate solutions of the form
(15a)$$ {}^{2} \textrm{V}_{\textrm{e}, i}(x) = \sqrt{\frac{1}{L}} \left[ \frac{\cos \left( \alpha_{\textrm{e}, i} x \right) }{\cos( \alpha_{\textrm{e}, i} L/2 ) } + \frac{ \cosh \left( \alpha_{\textrm{e}, i} x \right) }{\cosh(\alpha_{\textrm{e}, i} L/2)} \right],$$
(15b)$$ {}^{2}\textrm{V}_{\textrm{o}, i}(x) = \sqrt{\frac{1}{L}} \left[ \frac{\sin \left( \alpha_{\textrm{o}, i} x \right) }{\sin( \alpha_{\textrm{o}, i} L/2 ) } + \frac{ \sinh \left( \alpha_{\textrm{o}, i} x \right) }{\sinh(\alpha_{\textrm{o}, i} L/2)} \right] ,$$
where the allowed values of $\alpha _{\textrm {e}, i}$ and $\alpha _{\textrm {o}, i}$ correspond to the roots of the following transcendental equations: (16)$$\tan(\alpha_{\textrm{e}, i} L/2) ={-} \tanh(\alpha_{\textrm{e}, i} L/2) \quad \textrm{and} \quad \tan(\alpha_{\textrm{o}, i} L/2) = \tanh(\alpha_{\textrm{o}, i} L/2) .$$
The roots to Eq. (16) are well approximated by the simple expressions (17)$$\alpha_{\textrm{e}, i} \approx \frac{2\pi }{L} \left( i - \frac{1}{4}\right) \quad \textrm{and} \quad \alpha_{\textrm{o}, i} \approx \frac{2\pi}{L} \left( i + \frac{1}{4}\right),$$
for $i \ge 1$. The two degenerate solutions are given by ${}^{2} \textrm {V}_{\textrm {d},0}(x) = \sqrt {1/L}$ and ${}^{2} \textrm {V}_{\textrm {d},1}(x) = \sqrt {12/L^{3}} x$.It should be noted that the RDF elements for $N=2$ correspond directly to the vibrational modes of an unclamped one-dimensional bar [11–13]. As alluded to in Sec. 1, this basis is also a natural choice for the assessment of the validity of the perturbation model for the propagation of one-dimensional MSF structures in imaging systems and – as given in Eq. (12) – it has a valuable fourth spectral moment in that context [9].
3. RDF elements in two dimensions: flat space
The analysis in Secs. 2.1 and 2.2 is readily extended to find two-dimensional RDF bases for the case of a circular region $a$ with radius $R$, where plane polar coordinates $(r,\phi )$ are a natural choice. For this case, there are multiple ways to generalize higher derivatives from one to two dimensions. In keeping with the motivation of minimizing RMS-derived quantities of the surface structure, consider the following functional that is analogous to Eq. (2):
(18)$$\mathcal{L}(\textrm{V},\alpha) = \iint_a \left| \nabla^{N} \textrm{V}(r,\phi) \right|^{2} r \, {\textrm{d}}r \, {\textrm{d}} \phi - \alpha^{2N} \iint_a \textrm{V}^{2} (r,\phi) r \, {\textrm{d}}r \, {\textrm{d}} \phi,$$
where we adopt the notation (19)$$\nabla^{N} \triangleq \begin{cases} (\nabla^{2})^{N/2}, &N \textrm{ even}, \\ \nabla (\nabla^{2})^{(N-1)/2}, &N \textrm{ odd}. \end{cases}$$
This means that the vertical bars in the first integral of Eq. (18) represent either a scalar ($N$ even) or a vector ($N$ odd) absolute value.Similar variational calculations, described in Appendix A, can now be performed to ultimately arrive at the following $2N$-th order Sturm-Liouville differential (eigenvalue) system:
(20)$$\nabla^{2N} \textrm{V}(r,\phi) = (-\alpha^{2})^{N} \textrm{V}(r,\phi), \quad \hat{\mathfrak{D}}_j \textrm{V}(r,\phi) {\bigg|}_{r= R} = 0, \quad \textrm{for } j = N, \dots, 2N-1,$$
where the differential operator $\hat {\mathfrak {D}}_j$ is given by (21)$$\hat{\mathfrak{D}}_j \triangleq \begin{cases} (\nabla^{2})^{j/2}, &j \textrm{ even},\\ \partial_r (\nabla^{2})^{(j-1)/2}, &j \textrm{ odd}. \end{cases}$$
The differential equation in Eq. (20) can be solved by first performing a separation of variables where the solution is expressed as ${{}^{N} \!} \textrm {V}_{n}^{m}(r,\phi ) = {}^{N} \! \mathcal {R}_{n}^{m}(r)\Phi _m(\phi )$. [Once again, the Sturm-Liouville nature of Eq. (20) guarantees that the solutions are orthogonal and form a basis.] For all values of $N$, the azimuthal part of the solution can be written according to (22)$$\Phi_m (\phi) = \frac{1}{\sqrt{\pi}} \begin{cases} 1/\sqrt{2}, &m = 0, \\ \cos(m \phi), &m > 0, \\ \sin(|m|\phi), &m < 0, \end{cases}$$
where $m \in \mathbb {Z}$. It is always possible to label and order $n$ such that $n \in \mathbb {Z}^{+}$. The remaining ordinary differential equation for ${{}^{N} \!}\mathcal {R}_{n}^{m}(r)$, see Appendix B, has the solution (23)$${}^{N} \! \mathcal{R}_{n}^{m}(r) = {}^{N} \! \Gamma_{n}^{m} \sum_{\ell=0}^{N-1} c_\ell J_m \left[ ({-}1)^{\ell/N} \alpha_n^{m} r \right],$$
where $c_0 = 1$ and ${{}^{N} \!} \Gamma _{nm}$ is a normalization constant determined by (24)$$\iint_a \left| {{}^{N} \!} \textrm{V}_{n}^{m}(r,\phi) \right|^{2} r \, {\textrm{d}}r \, {\textrm{d}} \phi = 1.$$
The remaining values of $c_\ell$ and the allowed values of $\alpha _n^{m}$ are determined by the boundary conditions in Eq. (20).The differential (eigenvalue) system in Eq. (20), for a given value of $m$ and $N > 1$, has $\lfloor N/2 \rfloor$ solutions (degenerate for $N > 3$) with eigenvalue $\alpha = 0$. These are given by
(25)$${{}^{N} \!}\textrm{V}_{\textrm{d}, j}^{m}(r,\phi) = \sqrt{\frac{j+1}{\pi R^{2}}} \frac{r^{|m|}}{ R^{|m|}} Z^{m}_j \left( \frac{r^{2}}{R^{2}}\right) \Phi_m(\phi), \quad \textrm{for} \quad j = 0, \dots, \lfloor N/2 \rfloor - 1,$$
where the radial Zernike polynomial $Z^{m}_j(x)$ is defined as in [16]. In addition, for odd $N$ and $m=0$, there is one other solution for $\alpha =0$ of the form (26)$$\sqrt{\frac{(N+1)/2}{\pi R^{2}}} Z^{0}_{\frac{N-1}{2}} \left(\frac{r^{2}}{R^{2}} \right).$$
Once again, this degenerate subspace is all that remains as the solutions for Eq. (20) in the limit of $N \rightarrow \infty$. All these solutions satisfy the following orthogonality relations: (27a)$$ \iint_a {{}^{N} \!} \textrm{V}_n^{m} (r,\phi) \, {{}^{N} \!} \textrm{V}_{n'}^{m'}(r,\phi) \, r\, {\textrm{d}}r\,{\textrm{d}}\phi = \delta_{nn'}\,\delta_{mm'},$$
(27b)$$ \iint_a \nabla^{N}\left[{{}^{N} \!}\textrm{V}_n^{m} (r,\phi)\right]\cdot\nabla^{N}\left[ {{}^{N} \!}\textrm{V}_{n'}^{m'}(r,\phi)\right]\, r\, {\textrm{d}}r\,{\textrm{d}}\phi = \left(\alpha_n^{m}\right)^{2N} \, \delta_{nn'}\,\delta_{mm'}.$$
We again consider the cases of $N=1$ and $N=2$ explicitly.
3.1 N = 1
For $N=1$, the two-dimensional version of Eq. (2) becomes
(28)$$\mathcal{L}(\textrm{V},\alpha) = \iint_{a}|\nabla \textrm{V}(r,\phi)|^{2} r \, {\textrm{d}}r \, {\textrm{d}} \phi - \alpha^{2} \iint_a \textrm{V}^{2}(r,\phi)r \, {\textrm{d}}r \, {\textrm{d}} \phi.$$
The analogous variational calculation (see Appendix A) in two dimensions leads to the following second-order Sturm-Liouville partial differential system: (29)$$\nabla^{2} \,\textrm{V}(r, \phi) ={-} \alpha^{2} \textrm{V}(r,\phi), \quad \frac{\partial}{\partial r} \textrm{V}(r,\phi) {\bigg|}_{r=R} = 0.$$
Equation (29) can be solved to yield orthonormal solutions of the form ${}^{1} \textrm {V}_{n}^{m}(r,\phi ) = {}^{1}\mathcal {R}_{n}^{m}(r) \Phi _m(\phi )$, where (30)$${}^{1} \mathcal{R}_{n}^{m}(r) = \frac{ \sqrt{2} \thinspace J_m(\alpha_n^{m} r)}{ R \sqrt{ J_m^{2}(\alpha_n^{m} R) - J_{m-1}^{2}(\alpha_n^{m} R)}} ,$$
where $J_m$ is the $m$-th order regular Bessel function of the first kind. The allowed values of $\alpha _n^{m}$ are given by the roots of the following transcendental equation: (31)$$J_{m-1}(\alpha_n^{m}R) = J_{m+1} (\alpha_n^{m} R).$$
For a given value of $|m|$, these eigenvalues are approximately given by (32)$$\alpha_n^{m} \approx \frac{\pi}{R} \left( n + \frac{|m|}{2} - \frac{3}{4} \right),$$
for large $n$. The $\alpha = 0$ solution, as indicated in Eq. (26), is the constant function (33)$${}^{1} \textrm{V}_{\textrm{d},0}^{0} (r,\phi) = \frac{1}{\sqrt{\pi R^{2}}} Z_0^{0} \left( \frac{r^{2}}{R^{2}} \right) = \frac{1}{\sqrt{\pi R^{2}}}.$$
3.2 N = 2
For $N=2$, the two-dimensional version of Eq. (2) becomes
(34)$$\mathcal{L}(\textrm{V},\alpha) = \iint_a |\nabla^{2} \textrm{V}(r,\phi)|^{2} r \, {\textrm{d}}r \, {\textrm{d}} \phi - \alpha^{4} \iint_a \textrm{V}^{2}(r,\phi) r \, {\textrm{d}}r \, {\textrm{d}} \phi.$$
Once again, variational calculations can be performed to yield a fourth-order Sturm-Liouville partial differential equation: (35)$$\nabla^{4} \, \textrm{V}(r,\phi) = \alpha^{4} \textrm{V}(r,\phi), \quad \nabla^{2} \textrm{V}(r,\phi) {\bigg|}_{r = R} = 0, \quad \frac{\partial}{\partial r} [\nabla^{2} \textrm{V}(r,\phi)]{\bigg|}_{r=R} = 0.$$
Again with Eq. (22), Eq. (35) can be solved to yield orthonormal solutions of the form ${}^{2} \textrm {V}_{n}^{m}(r,\phi ) = {}^{2} \mathcal {R}_{n}^{m} (r)\,\Phi _m (\phi )$, where (36)$${}^{2} \mathcal{R}_{n}^{m}(r) = \frac{1}{R} \left[\frac{J_m (\alpha_n^{m} r)}{J_m (\alpha_n^{m} R)} + \frac{I_m (\alpha_n^{m} r)}{I_m (\alpha_n^{m} R)} \right],$$
for $n \ge 1$, where $I_m$ is the $m$-th order modified Bessel function of the first kind. The allowed (discrete) values of $\alpha _n^{m}$ are the roots of the following transcendental equation: (37)$$\frac{I_{m-1}(\alpha_n^{m} R)}{I_{m}(\alpha_n^{m} R)} = \frac{J_{m-1}(\alpha_n^{m} R)}{J_{m}(\alpha_n^{m} R)}.$$
For a given value of $|m|$, these eigenvalues satisfy (38)$$\alpha_n^{m} \approx \begin{cases} (n+|m|)\pi/R, &m \textrm{ even},\\ (n+ |m| - 1/2)\pi/R, &m \textrm{ odd}, \end{cases}$$
for large $n$. The solutions with $\alpha = 0$, as indicated by Eq. (25), are (39)$${}^{2} \textrm{V}_{\textrm{d},0}^{m}(r,\phi) = \frac{1}{\sqrt{\pi R^{2}}} \frac{r^{|m|}}{ R^{|m|}} Z^{m}_0 \left( \frac{r^{2}}{R^{2}}\right) \Phi_m(\phi) = \frac{1}{\sqrt{\pi R^{2}}} \frac{r^{|m|}}{R^{|m|}} \Phi_m(\phi).$$
Figures 3 and 4 shows some of these basis elements. In analogy with the one-dimensional solutions in Eqs. (15), the solutions given by ${}^{2} \textrm {V}_n^{m}(r,\phi )$ happen to describe the elastic modes of unclamped thin disks [11–13], and have been used to study deformable mirrors [14]. Furthermore, in two-dimensional paraxial imaging systems with MSF structures, these solutions are a natural basis for the assessment of the perturbation model [15]; the basis required for analogous studies of non-paraxial imaging systems is obtained by solving Eq. (35) on a sphere rather than on flat space, and this is considered in Sec. 4.4. RDF bases in two dimensions: spherical space
Equation (35) can also be solved with the spherical polar coordinate system $(\vartheta ,\phi )$, where $\vartheta$ is the colatitude (polar) angle and $\phi$ is once again the azimuthal angle. We define the region $a$ to be a spherical cap centered at the north pole, with colatitude boundary $\vartheta _0$. The derivation of the Sturm-Liouville differential system, analogous to that in Eq. (20), is shown in Appendix A. Furthermore, orthogonality relations like those in Eqs. (27) of the corresponding solutions to this differential system (the RDF basis elements for a spherical cap) are also stated in Appendix A.
For brevity, we analyze only the case $N = 2$, whose elements happen to correspond to the vibrational modes of an unclamped solid spherical cap. As mentioned at the end of Sec. 3.2, this basis arises naturally in the analysis of the perturbation model for MSF structures in non-paraxial imaging systems. For the coordinate system $(\vartheta ,\phi )$, the RDF elements are solutions to the following differential system:
(40)$$\nabla^{4} \, \textrm{V}(\vartheta,\phi) = \alpha^{4} \textrm{V}(\vartheta,\phi), \quad \nabla^{2} \textrm{V}(\vartheta,\phi) {\bigg|}_{\vartheta = \vartheta_0} = 0, \quad \frac{\partial}{\partial \vartheta} [\nabla^{2} \textrm{V}(\vartheta,\phi)] {\bigg|}_{\vartheta = \vartheta_0} =0,$$
where here $\nabla ^{2}$ denotes the Laplacian in the spherical space. In this case, we find that the RDF basis elements can be written in terms of the regular and modified spherical harmonics: (41)$${}^{2} \textrm{V}_{n}^{m}(\vartheta,\phi) = {}^{2} \Gamma_n^{m} \left\{ \frac{P^{m}_{f_+(\alpha_n^{m})} [\cos(\vartheta)]}{P^{m}_{f_+(\alpha_n^{m})} [\cos(\vartheta_0)]} + \frac{P^{m}_{f_-(\alpha_n^{m})}[\cos(\vartheta)]}{P^{m}_{f_-(\alpha_n^{m})} [\cos(\vartheta_0)]} \right\} \Phi_m(\phi),$$
where $m \in \mathbb {Z}$, $n \ge 1$, $P^{m}_{f_\pm }$ is an associated Legendre function with (42)$$f_\pm(\alpha_{nm}) \triangleq \frac{1}{2} \left[{-}1 + \sqrt{1 \pm 4 \left(\alpha_n^{m} \right)^{2}} \right].$$
The allowed values of the eigenvalue $\alpha _n^{m}$ correspond to the roots of the following transcendental equation: (43)$$\frac{(f_+{-}m + 1) P^{m}_{f_+{+}1}[\cos (\vartheta_0)]}{P^{m}_{f_+}[\cos (\vartheta_0)]} = (f_+{-} f_-)\cos (\vartheta_0) + \frac{(f_-{-}m + 1) P^{m}_{f_-{+} 1}[\cos (\vartheta_0)]}{P^{m}_{f_-}[\cos (\vartheta_0)]},$$
where the explicit dependence of $f_\pm$ on $\alpha _n^{m}$ in Eq. (43) has been left off for convenience. For a given value of $|m|$, these eigenvalues satisfy (44)$$\alpha_n^{m} \approx \begin{cases} (n+|m|)\pi/\vartheta_0, &m \textrm{ even},\\ (n+ |m| - 1/2)\pi/\vartheta_0, &m \textrm{ odd}, \end{cases}$$
for large $n$. The $\alpha = 0$ solutions are given by (45)$${}^{2} \textrm{V}^{m}_\textrm{d,0}(\vartheta,\phi) = \frac{\tan^{m}(\vartheta/2) \Phi_m(\phi)}{\sqrt{ 1 + \cos(\vartheta_0) + \frac{2(1-2m)}{1+2m}{}_2 F_1\left[1,\frac{1}{2} +m, \frac{3}{2} + m, - \tan^{2} (\vartheta_0/2) \right]\tan^{1+2m}(\vartheta_0/2) }},$$
where the denominator is a normalization constant with ${}_2F_1(a, b, c, z)$ being a hypergeometric function. Figure 5 shows some of these basis elements.5. Concluding remarks
We were motivated to pursue the ideas reported above by the fact that the fourth moment of an RDF-spectrum gives intuitive access to a physically significant characteristic length in investigations of MSF structure on optical surfaces [9,15]. As indicated in Eqs. (12) and more generally from (27b), these RDF bases allow the determination of the mean square $N$-th derivative of a function in terms of a spectral moment of order $2N$. For example, mean square slopes can be determined from the second moment of the spectra associated with the RDF bases presented in Secs. 2.1 and 3.1. Notice that this is quite different from a slope-orthogonal polynomial basis (akin to the Q polynomials [17], for example) where the mean-square slope is found by summing the squares of the spectral coefficients themselves. As pointed out in the Sec. 1, however, those polynomials are well matched to other aspects of the investigation of the characterization and modelling of MSF structure. By contrast, it follows from Eqs. (8a) and (27a), that the sum of squares of the RDF spectral coefficients is the mean square value of the decomposed function itself, see Eq. (10). In some contexts, it may also be of value to exploit the curious fact that the RDF bases are also bi-orthogonal in the more general sense of Eq. (55). Such properties are indicators of strong similarities to Fourier bases and give RDF bases potential value for a variety of applications that involve functions defined over finite domains. Indeed, certain subsets of RDF bases are already used in the study of elastic deformations of finite domains in different geometries [11–15].
Appendix A. Details of RDF derivation in two dimensions
In this appendix, we show the details of several steps in the derivation of the RDF bases in two dimensions. This will be done in general curvilinear coordinates $(u, v)$, and can be specialized to treat both the plane polar and spherical polar coordinates used in Secs. 3 and 4, respectively. We begin with
(46)$$\mathcal{L}(\textrm{V},\alpha) = \iint_a \left| \nabla^{N} \textrm{V}(u, v) \right|^{2} \, {\textrm{d}}A - \alpha^{2N} \iint_a \textrm{V}^{2} (u, v) \, {\textrm{d}}A,$$
where
${\textrm {d}}A$ is the infinitesimal area element of the corresponding metric associated with
$(u, v)$.
As in the one-dimensional case, one perturbs $V(u, v)$ by $\varepsilon (u, v)$ and requires $\mathcal {L}$ to be stationary (unchanged to first order), hence
(47)$$\mathcal{L}(V+\varepsilon,\alpha) = \mathcal{L}(V,\alpha) + 2 \left[ \iint_a \left(\nabla^{N} V\right) \cdot \left(\nabla^{N} \varepsilon \right) \! {\textrm{d}}A - \alpha^{2N} \!\!\iint_a V \epsilon \, {\textrm{d}}A \right] + O(\varepsilon^{2}),$$
where the
$\cdot$ operation is a scalar product between two scalar quantities or two vector quantities, depending on whether
$N$ is even or odd, respectively. Stationarity then requires the quantity within the square brackets of Eq. (
47) to vanish for any
$\varepsilon$. The first step in this process is to transfer the derivatives from
$\varepsilon$ over to
$V$ in the first integral within the square brackets. That integral can be written, using the product rule, as
(48)$$\iint_a \left(\nabla^{N} V\right) \cdot \left(\nabla^{N} \varepsilon \right) \, {\textrm{d}}A = \iint_a \nabla \cdot \left[ \left(\nabla^{N} V \right) \left( \nabla^{N-1} \varepsilon\right) \right] \, {\textrm{d}}A - \iint_a \nabla^{N+1}V \cdot \nabla^{N-1} \varepsilon \, {\textrm{d}}A.$$
(Note that, depending on the parity of
$N$, either the product within the integrand of the left-hand side or that in the second term of the right-hand side are a scalar product of two vectors while the other is a multiplication of two scalars.) We now use Gauss’s theorem to express the first integral on the right side of Eq. (
48) as a boundary term, hence
(49)$$\iint_a \left(\nabla^{N} V\right) \cdot \left(\nabla^{N} \varepsilon \right) \, {\textrm{d}}A = \int_{\partial a} \left[ \left(\nabla^{N} V \right) \left( \nabla^{N-1} \varepsilon\right) \right] \cdot {\textrm{d}} \vec{S}_\perp{-} \iint_a \nabla^{N+1}V \cdot \nabla^{N-1} \varepsilon \, {\textrm{d}}A.$$
At this point, one derivative has been transferred from
$\varepsilon$ to
$V$, but we are left with a boundary term that must vanish for any
$\varepsilon$. This condition is given by
(50)$$\hat{\mathfrak{D}}_N V {\bigg|}_{\partial a} =0$$
with
(51)$$\hat{\mathfrak{D}}_j \triangleq \begin{cases} (\nabla^{2})^{j/2}, &j \textrm{ even},\\ \hat{n} \cdot \nabla (\nabla^{2})^{(j-1)/2}, &j \textrm{ odd}, \end{cases}$$
where
$\hat {n}$ a unit vector that is parallel to
${\textrm {d}} \vec {S}_\perp$. Note that Eq. (
51) is consistent with Eq. (
21) when
$a$ is chosen to be a circular region of radius
$R$ and
$(u, v) = (r,\phi )$, the plane polar coordinates.
Repeating the process from Eqs. (48) to (49) $N$ times and forcing each boundary term to vanish, the quantity in the square brackets in Eq. (47) ultimately becomes
(52)$$\iint_a \left[({-}1)^{N} \nabla^{2N} V - \alpha^{2N} V \right] \varepsilon \, {\textrm{d}}A.$$
By requiring Eq. (
52) and the boundary terms to vanish for any
$\varepsilon$, we end up with the following differential system:
(53)$$\nabla^{2N} V = (-\alpha^{2})^{N} V, \quad \hat{\mathfrak{D}}_j V {\bigg|}_{\partial a} = 0, \quad \textrm{for } j=N,\dots,2N-1.$$
Equation (
53) is a coordinate-free form of Eq. (
20) in Sec.
3, which uses plane polar coordinates. For Sec.
4, the coordinates
$(\vartheta ,\phi )$ are used and
$\hat {n} \cdot \nabla = \partial _\vartheta$ in the definition of
$\hat {\mathfrak {D}}_j$.
Since Eq. (53) is a Sturm-Liouville problem, it has solutions ${{}^{N} \!} \textrm {V}_n^{m}(u, v)$ with eigenvalue $\alpha _n^{m}$, where $n$ and $m$ are integers. These solutions are orthogonal with respect to two inner products:
(54a)$$ \iint_a {{}^{N} \!}\textrm{V}_n^{m}(u, v) \, {{}^{N} \!}\textrm{V}_{n'}^{m'}(u, v)\, {\textrm{d}}A = \delta_{nn'}\delta_{mm'} ,$$
(54b)$$\iint_a \nabla^{N}\left[{{}^{N} \!}\textrm{V}_n^{m}(u, v)\right] \cdot \nabla^{N} \left[{{}^{N} \!}\textrm{V}_{n'}^{m'}(u, v)\right]\, {\textrm{d}}A = \left(\alpha_n^{m}\right)^{2N}\delta_{nn'}\delta_{mm'}.$$
We note for completeness that the solutions
${{}^{N} \!} \textrm {V}_n^{m}(u, v)$ also satisfy the following more general orthogonality relation:
(55)$$\iint_a \nabla^{N_1}\left[{{}^{N} \!}\textrm{V}_n^{m}(u, v)\right] \cdot \nabla^{N_2} \left[{{}^{N} \!}\textrm{V}_{n'}^{m'}(u, v)\right]\, {\textrm{d}}A = \left(\alpha_n^{m}\right)^{(N_1+N_2)}\delta_{nn'}\delta_{mm'},$$
valid for any non-negative integer
$N_1$ and
$N_2$ for which
$N_1 + N_2$ is an integer multiple of
$2N$. This relation includes Eqs. (54) as special cases. The one-dimensional analogue for these relations arises by replacing
$\nabla$ with a simple derivative.
For completeness, we state the gradient and divergence for plane polar coordinates $(r,\phi )$ and the spherical polar coordinates $(\vartheta ,\phi )$. The former is given by
(56)$$\nabla g(r,\phi) = \left(\partial_r g, \frac{1}{r} \partial_\phi g \right),$$
(57)$$\nabla \cdot [f_r(r,\phi), f_\phi(r,\phi)] = \frac{1}{r} \partial_r \left( r f_r \right) + \frac{1}{r} \partial_\phi f_\phi,$$
and the latter by
(58)$$\nabla g(\vartheta,\phi) = \left(\partial_\vartheta g, \frac{1}{\sin \vartheta} \partial_\phi g \right),$$
(59)$$\nabla \cdot [f_\vartheta(\vartheta,\phi), f_\phi(\vartheta,\phi)] = \frac{1}{\sin (\vartheta)} \partial_\vartheta \left[ \sin (\vartheta) f_\vartheta \right] + \frac{1}{\sin ( \vartheta)} \partial_\phi f_\phi.$$
The corresponding area elements are given by
${\textrm {d}}A = r {\textrm {d}}r {\textrm {d}} \phi$ and
${\textrm {d}}A = \sin (\vartheta ) {\textrm {d}} \vartheta {\textrm {d}} \phi$, respectively.
Appendix B. Separation of variables
Here, we show that the method of separation of variables may be used to solve the differential equation in Eq. (53) in the specific case of plane polar coordinates $(r,\phi )$; the procedure in spherical polar coordinates is similar since both coordinate systems have the same azimuthal coordinate $\phi$.
It follows from Eqs. (56) and (57) that the differential equation in Eq. (53) is
(60)$$\left[\frac{1}{r} \partial_r \left( r \partial_r \right) + \frac{1}{r^{2}} \partial^{2}_\phi \right]^{N} V(r,\phi) = (-\alpha^{2})^{N} V(r,\phi).$$
By assuming a solution of the form
$V(r,\phi ) = \mathcal {R}(r)\,\Phi _m(\phi )$, Eq. (
60) can be written as
(61)$$\begin{aligned} \left[\frac{1}{r} \partial_r \left( r \partial_r \right) + \frac{1}{r^{2}} \partial^{2}_\phi \right]^{N} \mathcal{R}(r)\,\Phi_m(\phi) &= (-\alpha^{2})^{N} \mathcal{R}(r)\,\Phi_m(\phi)\\ \Rightarrow \left[\frac{1}{r} \partial_r \left( r \partial_r \right) - \frac{m^{2}}{r^{2}} \right]^{N} \mathcal{R}(r) &= (-\alpha^{2})^{N} \mathcal{R}(r). \end{aligned}$$
Equation (
61) is an ordinary differential equation for
$\mathcal {R}(r)$ and its solutions are given by Eq. (
23).
Funding
National Science Foundation (NSF) I/UCRC Center for Freeform Optics (IIP-1338877); Aix-Marseille Université.
Acknowledgments
MAA received funding from the Excellence Initiative of Aix-Marseille University – A*MIDEX, a French “Investissements d’Avenir” programme.
Disclosures
The authors declare that there are no conflicts of interest related to this article.
References
1. L. Cohen, Time-Frequency Analysis (Prentice Hall PTR, 1995).
2. C. K. Chui, An Introduction to Wavelets (Academic, 1995).
3. J. J. Healy, M. A. Kutay, H. M. Ozaktas, and J. T. Sheridan, Linear Canonical Transforms, Theory and Applications (Springer Science+Business Media, 2016), Chaps. 1–3
4. A. Zettl, Sturm-Liouville Theory (American Mathematical Society, 2005).
5. F. Z. Fang, X. D. Zhang, A. Weckenmann, G. X. Zhang, and C. Evans, “Manufacturing and measurement of freeform optics,” CIRP Ann. 62(2), 823–846 (2013). [CrossRef]
6. K. Garrard, T. Bruegge, J. Hoffman, T. Dow, and A. Sohn, “Design tools for freeform optics,” Proc. SPIE 5874, 58740A (2005). [CrossRef]
7. R. E. Fischer, B. Tadic-Galeb, and P. R. Yoder, Optical System Design (McGraw Hill, 2000), Chap. 7
8. G. W. Forbes, “Fitting freeform shapes with orthogonal bases,” Opt. Express 21(16), 19061 (2013). [CrossRef]
9. K. Liang, G. W. Forbes, and M. A. Alonso, “Validity of the perturbation model for the propagation of MSF structure in 2D,” Opt. Express 27(3), 3390 (2019). [CrossRef]
10. H. Wang and S. Xiang, “On the Convergence Rates of Legendre Approximation,” Math. Comp. 81(278), 861–877 (2011). [CrossRef]
11. T. L. Schmitz and K. S. Smith, Machining Dynamics: Frequency Response to Improved Productivity (Springer Science & Business Media, 2008), Chap. 7.4
12. F. A. Khasawneh and D. Segalman, “Exact and Numerically Stable Expressions for Euler-Bernoulli and Timoshenko Beam Modes,” arXiv:1811.03222v3 [physics.class-ph] 28, (2019).
13. M. Avcar, “Free Vibration Analysis of Beams Considering Different Geometric Characteristics and Boundary Conditions,” Internat. J. Mechan. Appl. 4(3), 94–100 (2014). [CrossRef]
14. H. Wang, M. Zhang, Y. Zuo, and X. Zheng, “Research on elastic modes of circular deformable mirror for adaptive optics and active optics corrections,” Opt. Express 27(2), 404–415 (2019). [CrossRef]
15. K. Liang, G. W. Forbes, and M. A. Alonso, “Validity of the perturbation model for the propagation of MSF structure in 3D,” Opt. Express 27(3), 3390 (2019). [CrossRef]
16. G. W. Forbes, “Robust and fast computation for the polynomials of optics,” Opt. Express 18(13), 13851 (2010). [CrossRef]
17. G. W. Forbes, “Characterizing the shape of freeform optics,” Opt. Express 20(3), 2483 (2012). [CrossRef]
Open Access
Add to BibSonomy
