# Assessment about Luneberg integrals and application to digital in-line holography

Sébastien Coëtmellec, Denis Lebrun, Marc Brunel, A.J.E.M. Janssen

► To cite this version:

Sébastien Coëtmellec, Denis Lebrun, Marc Brunel, A.J.E.M. Janssen. Assessment about Luneberg integrals and application to digital in-line holography. Journal of the European Optical Society : Rapid publications, 2021, 10.1186/s41476-021-00154-x. http://hal-03275428

HAL Id: hal-03275428

https://hal.science/hal-03275428v1

Submitted on 12 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Assessment about Luneberg integrals and application to digital in-line holography

<figure><img src="image_2.png" /></figure>

Check for
updates

S. Coëtmellec¹* , D. Lebrun¹, M. Brunel¹ and A. J. E. M. Janssen²

**Abstract**

In this publication, the Luneberg integrals are revisited and the conditions of the using of such integrals have been recalled. Additivity law of Luneberg's integrals and the link with the Frenel kernel for the propagation are discussed. By means of the definition of the Luneberg's integrals, the propagation of a vectorial electromagnetic field (Hertz potentials) is developed and a new approach of the computation have been proposed based on Zernike polynomials. With this new approach simulations of holograms is illustrated in the case of the digital in-line holography with an opaque disk.

**Keywords:** Luneberg's integrals, Digital in-holography, Zernike polynomials

**Introduction**

Digital in-line holography often uses Fresnel and Fourier integrals, either to describe the light propagation and to calculate to the different expressions of amplitude distribution or to reconstruct the image of the object from the recorded holograms. In the case of in-line holography, these integral operators consider three main essential parameters for the holographists: the distance z between the object and the recording plane of the hologram, the diameter D of the object [12] or the wavelengths λ in the case of multi-wavelength holography [11, 15, 27]. However, these descriptions do not include finer parameters of optics such as light polarization. The models used in holography are indeed often based on scalar diffraction models. Some methods using machine learning [29] appeared recently in three-dimensional vector holography studies. Although all these studies yields to remarkable results, the vectorial nature of light is considered too complex to propose and use explicit descriptions. However, works have already been carried out and are yielding remarkable results [6, 8, 31]. Taking into account polarization in digital holography would enrich the domain

of applicability and the accuracy of this technique. The Rayleigh-Sommerfeld integrals, type I or type II depending on the circumstances (near field, far field) are the most naturally used integral forms in vector diffraction models [30].

In this article, we develop a vectorial model applied to the digital in-line holography. We use the Luneberg integral that tends to the second Rayleigh-Sommerfeld diffraction integral [20] in the far-field approximation [17]. Luneberg adressed rigorously and with the highest level of attention the boundary value problem of the equation $\Delta u + k^2 u = 0$ for a plane boundary. Two conditions were obtained and will be recalled and analysed in our article as they are very often neglected in the studies. For example, Luneberg conditions are not checked in the case of Gaussian beam.

For several years now, the objects under study by holography have become even smaller in size to reach nanometric scales [26, 35]. It is therefore clear that the ratio $D/\lambda$ becomes a central point of these holography studies. The Gaussian model cannot be used in the context of wavelength or sub-wavelength particles [22] and the scalar model is not sufficient [4].

This is the reason why, based on the work of the authors of [10], we develop here a vector electromagnetic field model allowing to describe as precisely as possible the

<figure><img src="image_3.png" /></figure>

<footnote>*Correspondence: coetmel@coria.fr
¹Département d'optique, UMR-6614, Av. de l'Université, 76801 Saint-Etienne du Rouvray, France
Full list of author information is available at the end of the article</footnote>
holographic patterns delivered by such small objects. In Ref. [10], a first order approximation of the Luneberg's kernel is proposed to obtain a solution to the propagation of the electromagnetic field vector. This is why in most publications authors use Luneberg, respectively Rayleigh-Sommerfeld, and then carry out a far field approximation. The integrals obtained lose their interest in the near field with objects whose dimensions do not exceed a few wavelengths. Here, we preserve the structure of the Luneberg integrals while using a vector model of the electromagnetic field that respects the conditions of Luneberg and thus allows to describe the interaction with objects of a few wavelengths. The Zernike polynomials are used as a basis of decomposition of the electromagnetic field.

This article is organized as follows. First, the Luneberg's integrals and their properties are recalled. The conditions of our study are recalled as well. The formalism is then illustrated with the description of the transition from Luneberg to Fresnel models. The third part is dedicated to recall the results of [10] by explaining all the terms of the electromagnetic field. In our study, this field interacts with an opaque particle in the sense of holography, which is here approximated by an opaque disc. The question of the Luneberg transform of the electromagnetic field on a finite and infinite opening is treated and we propose semi analytical and analytical solutions. These solutions are then used to simulated holograms in various polarization configurations.

## Luneberg's integrals

### Background and definition

It should be noted that an optical system for digital in-line holography (DIH) is composed of two parts. Each part is delimited by the source, the object and a CCD sensor. An illustration is given in Fig. 1. The center of the Cartesian coordinates system (x, y, z) is placed for example at the beam waist behind a focusing lens. Around this center of

the system, a dashed line circle, denoted Γ, is drawn and illustrates the constant phase of a spherical wave of radius r. The spherical wave is thus centered at zero. In this article, the spherical Hertz potentials are used to define the vectorial spherical illuminating the object, illustrated in Fig. 1, by a black disk of diameter D is located at the coordinates (ξ, η, δ). The CCD sensor records is placed at the distance zl from the object. This sensor records the intensity of the field (i.e. the hologram of the object) hologram of the object.

A stationary regime of monochromatic electromagnetic waves in a homogeneous, isotropic, linear, medium without charges and without currents is considered. Let P(ξ, η) be the complex amplitude of the object. This object is then supposed to be illuminated by a vectorial monochromatic electromagnetic wave, denoted by $E_δ = (E_δ^x, E_δ^y, E_δ^z)$, that propagates through linear and homogeneous medium in regions free of currents and charges. The notation in bold denotes the vectorial nature of the function and linked to the polarization state of the wave. From that, the amplitude of the vectorial wave, denoted $W_δ = (W_δ^x, W_δ^y, W_δ^z)$ just next the object is given by

$$W_{\delta}(\xi, \eta) = P(\xi, \eta) E_{\delta}(\xi, \eta). \quad (1)$$

Note that if the object P(ξ, η) does not modify the polarization state of the incident wave and if its size is in order to ten times of the wavelength or more, the taking into account of the polarization is of no interest from the Luneberg operator point of view. Nevertheless, two situations present a major interest in the taking into account of the polarization and the Hertz potentials method. The first is when the object is deformed for example following a stretching or when its orientation is modified with respect to the wave vector. In the first case, a birefringence may appear and in the second the polarization of the incident wave can be modified by reflexion even if the object is considered as a scalar. The polarization state

<figure><img src="image_4.png" /><figcaption>Fig. 1 Records hypothetical experiment here with δ > 0</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 2 Profiles of the intensity distribution of holograms simulated with Fresnel and Luneberg integrals with D = 20µm, zl = 18.98mm, δ = 0mm, N=1024 and λ=633nm</figcaption></figure>

of the incident wave will therefore be modified and justifies the using of the Luneberg operator. For the second, the Hertz potentials method is justified by the fact that if the object is very small in size, of the order of the wavelength dimension, scalar models are not appropriate. In the previous example, the object becomes smaller and smaller during stretching and its dimension may become equal to a few wavelengths. In all these cases, the function $W_δ(ξ, η)$ becomes vectorial. In the case where the polarization is taken into account and in the presence of discontinuous plane, the Luneberg integrals, derived from the integral theorem of Helmholtz and Kirchoff [5], Eq. 7 on p. 419, can be used to describe the propagation next the object and their definitions in [20], Eqs. (45.75) on p. 320, according to the system illustrated in Fig. 1 are the following

$$ \mathcal{E}_{z_l}^x(x, y) = -\frac{z_l}{2\pi} \iint W_\delta^x(\xi, \eta) \left(\frac{ikr-1}{r^3}\right) \exp(ikr) d\xi d\eta, (2) $$

$$ \mathcal{E}_{z_l}^y(x, y) = -\frac{z_l}{2\pi} \iint W_\delta^y(\xi, \eta) \left(\frac{ikr-1}{r^3}\right) \exp(ikr) d\xi d\eta, (3) $$

where $r^2 = (x - \xi)^2 + (y - \eta)^2 + z_l^2$ with $z_l = z - \delta > 0$ due to the $\delta$-shifting along the optical axis. If the third component is undetermined, it could be derived from (2) and (3) by means of the Maxwell's equation div $\mathcal{E}_{z_l} = 0$ (see Appendix A). We get

$$ \mathcal{E}_{z_l}^z(x, y) = \frac{1}{2\pi} \iint \left[ (x - \xi) W_\delta^x(\xi, \eta) + (y - \eta) W_\delta^y(\xi, \eta) \right] \left( \frac{ikr - 1}{r^3} \right) \exp(ikr) d\xi d\eta, \quad (4) $$

where $k = 2\pi/\lambda$ is the wave number and $\mathcal{E}_{z_l} = (\mathcal{E}_{z_l}^x, \mathcal{E}_{z_l}^y, \mathcal{E}_{z_l}^z)$ is the vectorial field in the sensor plane. The wave $\mathcal{E}_{z_l}$ satisfies the Helmholtz equation in the half space $z_l > 0$. The uniqueness of the solution $\mathcal{E}_{z_l}^i$ depends on conditions as mentioned in [20], Eqs. (45.14) and (45.141) on p. 312. In the domain of the half space $z_l > 0$ there exist a constant $A$ such that the first conditions are

$$ |\mathcal{E}_{z_l}^i| < \frac{A}{r}, \quad \left|\frac{\partial \mathcal{E}_{z_l}^i}{\partial r}\right| < \frac{A}{r}, \qquad (5) $$

and in any solid sector $\pi/2 + \Omega < \theta < \pi/2 - \Omega$ of the domain $z_l > 0$, there exists a constant $B(\Omega)$ such that for all points $(x, y, z_l)$ of the sector we have the second condition

$$ \left| \frac{\partial \mathcal{E}_{z_l}^i}{\partial r} - ik \mathcal{E}_{z_l}^i \right| < \frac{B(\Omega)}{r^2}, \quad r = \sqrt{x^2 + y^2 + z^2}, \quad z \ge \delta. \qquad (6) $$

It is important to note that the conditions in Eq. (5) apply on the function $W_\delta^i(\xi, \eta)$ too. Indeed, for $z_l = 0$, we have the boundary limit (see Appendix B)

$$ \lim_{z_l \to 0} \mathcal{E}_{z_l}^{x,y}(x, y) = W_{\delta}^{x,y}(\xi, \eta). \qquad (7) $$

Then, these boundary values are not arbitrary and they must satisfy the conditions established by Eqs. (5) and (6). A Gaussian beam can be used but that corresponds to the paraxial Maxwell equation [19] and do not satisfy the conditions. In the case of a very near field (of the order of a few wavelengths) or for very focused beams, the Gaussian
beams are not a solution [13, 14]. The use of Hertz potentials allow to obtain a vectorial spherical wave which satisfy the Luneberg conditions. From simulation point of view and by means of fast Fourier transform (FFT), it is convenient to write the Luneberg integrals as convolution products. So, by noting $L_{z_l}(x, y)$ the Luneberg kernel, we have for Eqs. (2) and (3)

$$ \mathcal{E}_{z_l}^i(x, y) = W_{\delta}^i(x, y) * L_{z_l}(x, y), \quad i = x, y, \quad (8) $$

with

$$ \begin{aligned} L_z(x, y) &= -\frac{z}{2\pi} \left( \frac{ikr - 1}{r^3} \right) \exp(ikr) \\ &= -\frac{1}{2\pi} \frac{\partial}{\partial z} \left( \frac{\exp(ikr)}{r} \right), \end{aligned} \quad (9) $$

and for Eq. (3)

$$ \begin{aligned} \mathcal{E}_{z_l}^z(x, y) &= -W_{\delta}^x(x, y) * \left[ \frac{x}{z_l} \cdot L_{z_l}(x, y) \right] \\ &\quad - W_{\delta}^y(x, y) * \left[ \frac{y}{z_l} \cdot L_{z_l}(x, y) \right], \end{aligned} \quad (10) $$

### Additive law

Let us recall here an important relation called additive law. This relation is interesting because the axial distance $z_l$ can be viewed as an index from point of view of the Luneberg integrals. The continuity of the propagation of the wave corresponds to the continuity of the index and thus the additive law of the indices. If the Luneberg operator on the transverse components is denoted by $\mathcal{L}_{z_l}$ then

$$ \mathcal{E}_{z_l}^i(x, y) = \mathcal{L}_{z_l} [W_{\delta}^i(\xi, \eta)] := W_{\delta}^i(x, y) * L_{z_l}(x, y). \quad (11) $$

The propagation of a spherical wave between two planes separated by a distance $z_l > 0$ could be decomposed into two propagations of $z_{l1} > 0$ and $z_{l2} > 0$ with $z_l = z_{l1} + z_{l2}$

$$ \mathcal{L}_{z_l} [W_{\delta}^i] = \mathcal{L}_{z_{l1}+z_{l2}} [W_{\delta}^i] = \mathcal{L}_{z_{l2}} [\mathcal{L}_{z_{l1}} [W_{\delta}^i]]. \quad (12) $$

We have the continuity of wave propagation. To prove this, we start from Eq. (8) in the Fourier domain. The used definition of the Fourier transform is expressed as

$$ \begin{aligned} \hat{f}(u, v) &= \mathcal{F}[f(x, y)](u, v) \\ &= \iint_{-\infty}^{+\infty} f(x, y) \exp[-i2\pi(ux + vy)]dx dy, \end{aligned} \quad (13) $$

where $\hat{f}(u, v)$ is the spectrum amplitude of $f(x, y)$ and $(u, v)$ are the spatial frequency coordinates. According to the property of the Fourier transformation of a convolution product, in the spectral domain, Eq. (11) over the distance $z_{l1}$ leads to

$$ \hat{\mathcal{E}}_{z_{l1}}^i(u, v) = \mathcal{F}[W_{\delta}^i(x, y)] \cdot \mathcal{F}[L_{z_{l1}}(x, y)]. \quad (14) $$

The evaluation of the Fourier transform of the Luneberg kernel is realized by means of [9], Eq. (A8), p. 112 and gives us

$$ \mathcal{F}[L_{z_{l1}}(x, y)] = \text{sign}(z_{l1}) \exp\left[i|z_{l1}|\sqrt{k^2 - 4\pi^2(u^2 + v^2)}\right]. \quad (15) $$

In a similar way, the propagation of $\hat{\mathcal{E}}_{z_{l1}}^i(u, v)$ over the distance $z_{l2}$, gives

$$ \hat{\mathcal{E}}_{z_{l2}}^z(u, v) = \text{sign}(z_{l2}) \hat{\mathcal{E}}_{z_{l1}}^i(u, v) \exp\left[i|z_{l2}|\sqrt{k^2 - 4\pi^2(u^2 + v^2)}\right]. \quad (16) $$

By introducing Eq. (14) into Eq. (16), we obtain the amplitude spectrum at the distance $z_l$ such as

$$ \hat{\mathcal{E}}_{z_l}^z(u, v) = \mathcal{F}[W_{\delta}^i(x, y)] \cdot \text{sign}(z_{l1}) \text{sign}(z_{l2}) \exp\left[i(|z_{l1}| + |z_{l2}|)\sqrt{k^2 - 4\pi^2(u^2 + v^2)}\right]. \quad (17) $$

Knowing that $z_{l1} > 0$ and $z_{l2} > 0$, in the spatial domain, we have

$$ \mathcal{E}_{z_l}^i(x, y) = W_{\delta}^i(x, y) * L_{z_{l1}+z_{l2}}(x, y) = \mathcal{L}_{z_{l1}+z_{l2}}[W_{\delta}^i]. \quad (18) $$

and this is Eq. (12). Note that $\mathcal{L}_0$ is the identity operator.

### From luneberg integrals to fresnel integral

In digital in-line holography (DIH), the wave propagation is commonly described by the Fresnel integral in the convolution form. The definition is given by

$$ \mathcal{E}_{z_l}^i(x, y) = W_{\delta}^i(x, y) * H_{z_l}(x, y), \quad (19) $$

where $\mathcal{E}_{z_l}^i$ and $W_{\delta}^i$ are scalar components and

$$ H_z(x, y) = \frac{\exp(ikz)}{i\lambda z} \exp\left(i\frac{\pi}{\lambda z}(x^2 + y^2)\right). \quad (20) $$

To study the Luneberg integrals in the far field approximation, it is sufficient to compare their kernel. We often see in the literature the case where only one axis of polarization is considered. So, here this situation is studied to illustrate this transition between the simulated hologram with Fresnel integral and Luneberg integrals. The first order Taylor series expansion of the Luneberg kernel gives

$$ L_z(x, y) \approx H_z(x, y), \quad (21) $$

under the following condition, see Appendix C

$$ z \gg \max\left(\frac{3\lambda}{4\pi}, \frac{1}{4}N\lambda\right), \quad (22) $$

where the first term is a physical parameter and the second is linked to digital consideration where N corresponds to the number of samples contained in the image.
<figure><img src="image_6.png" /><figcaption>Fig. 3 Comparison of the electric field profiles obtained from the convolution and the analytical calculus with δ = 5λ, zl = 3λ and λ = 633nm</figcaption></figure>

This condition can be used to delimit the two Fresnel and Luneberg domains. To illustrate the previous consideration about the kernels, holograms with Fresnel and Luneberg integrals have been simulated where the object is an opaque disk of diameter D illuminated by a wave of amplitude unity, i.e. $E_δ = (E_δ^x, 0, 0) = (1, 0, 0)$ with $δ = 0$. The amplitude function $P(ξ, η) = 1 - \text{circ}(r/R)$, with $\text{circ}(r/R) = 1$ if $r = (x^2 + y^2)^{1/2} \le R$ and 0 otherwise, $R = D/2$. Of course, the choice of an amplitude unity do not satisfies the conditions in Eqs. (5) and (6) but this allows us to demonstrate the behavior of the Luneberg integrals as regard to one of the two transverse components. So, only one polarization state of the wave has been considered to retrieve the scalar condition of the Fresnel integral. The hologram represented in Fig. 2 has been simulated with $zl = 18.98$mm. It can be seen from Fig. 2 that the difference between the holograms simulated with Fresnel and Luneberg is negligible in the far field domain. The offset difference is probably due to the amplitude differences at the origin as mentioned by the authors of [21] in Fig. (9) on p. 515. Nevertheless, if we consider a shorter distance $zl$ of several wavelengths, the Fresnel integral applied to scalar waves is not adapted compared to the Luneberg integral. The next section is then devoted to combine the vectorial, spherical waves and Luneberg integrals.

## Holography with vectorial spherical waves

### Definition of the vectorial spherical waves

The Luneberg conditions require that the object is illuminated by a spherical wave. This polarized spherical wave must also satisfy the Maxwell's equations as mentioned above. The authors of [10] have proposed a solution that allows both to satisfy Maxwell's equations but also to satisfy the conditions of Luneberg whatever the polarization state. This solution uses the Hertz potentials. So, here we recall the results obtained by adding some supplements and a few modifications. The definition of the vectorial wave in [10], Eqs. (11) and (12) on p. 734 are the following

$$E_δ(ξ, η) = [M_e(ξ, η, δ) \cdot p_e + M_m(ξ, η, δ) \cdot p_m], \quad (23)$$

$$B_δ(ξ, η) = \frac{1}{c} [-M_m(ξ, η, δ) \cdot p_e + M_e(ξ, η, δ) \cdot p_m], \quad (24)$$

where the velocity of light in the vacuum is approximately equal to $c = 3 \cdot 10^8 m/s$ and the electric vector and magnetic induction are denoted by ($E_δ$, $B_δ$). As shown by Fig. 1, the radial coordinate $r = (ρ^2 + δ^2)^{1/2}$ are considered with the circular coordinates $(ξ, η) = (ρ \cos θ, ρ \sin θ)$. The polarization states are defined by $p_e^T = (p_{ex} \ p_{ey} \ p_{ez})$ for electric vector and $p_m^T =$
<figure><img src="image_7.png" /><figcaption>Fig. 4 Components of the electric and magnetic fields ($E_δ$, $B_δ$) with $δ = 5λ$, $D = 7λ$ and the polarization $p_e = (1, 0, 0)$, $p_m = (0, 1, 0)$</figcaption></figure>

($p_{mx}$, $p_{my}$, $p_{mz}$) for magnetic induction. The upper symbol T denotes the transpose. The vectorial electromagnetic fields satisfy the Maxwell's equation in regions free of currents and charges. For this to happen many different steps are needed. Firstly, the Maxwell's equation must be rewriting in terms of the magnetic vector and electric scalar potentials $Φ_e$ and $A_m$ under the Lorentz gauge. Consequently, and next straightforward process, we have

$$E_{\delta}(\xi, \eta) = -\operatorname{grad}\Phi_{e} - \frac{\partial A_{m}}{\partial t}, \qquad (25)$$

$$B_{\delta}(\xi, \eta) = \operatorname{curl} A_{m}, \qquad (26)$$

with $A_m = \mu_0\varepsilon_0 \frac{\partial \Pi_e}{\partial t} + \operatorname{curl} \Pi_m$. The magnetic vector and electric scalar potentials are expressed versus the Hertz electric potential, denoted $\Pi_e(r, t)$, and the magnetic potential, denoted $\Pi_m(r, t)$ by writing that $Φ_e = -\operatorname{div} \Pi_e(r, t)$. The expressions of the vectorial electromagnetic fields ($E_δ$, $B_δ$) versus ($\Pi_e$, $\Pi_m$) in $\exp(-i\omega t)$-Fourier space are the following

$$E_{\delta}(\xi, \eta) = \operatorname{curl} \operatorname{curl} \Pi_{e}(r) + i\omega \operatorname{curl} \Pi_{m}(r), \qquad (27)$$

$$B_{\delta}(\xi, \eta) = \operatorname{curl} \operatorname{curl} \Pi_{m}(r) - i\omega\mu_{0}\varepsilon_{0}\operatorname{curl} \Pi_{m}(r), \qquad (28)$$

where $\mu_0$, $\varepsilon_0$ are the magnetic permeability and permittivity in the vacuum. The Hertz potentials have to satisfy the Helmholtz equation on the one hand and the other

satisfy the Luneberg conditions in Eqs. (5) and (6). Consequently, the possible solutions for the expressions of the Hertz potentials are

$$\Pi_e(r) = \frac{\exp(ikr)}{r} \cdot p_e, \qquad (29)$$

$$\Pi_m(r) = \frac{\exp(ikr)}{cr} \cdot p_m, \qquad (30)$$

which have a spherical amplitude weighted by a polarization vector. The parameter $ikr$ is taken positive to obtain a divergent spherical wave, for example next focusing. A negative value allows to have a convergent wave. The celerity $c$ has been introduced to respect the units of the fields: $[E_δ] = V \cdot m^{-1}$ and $[B_δ] = V \cdot s \cdot m^2$ with $[\Pi_e] = V \cdot m$ and $[\Pi_m] = V \cdot s$. By introducing the Eqs. (30) and (29) into (28) and (27), the expression of the matrices $M_e$ and $M_m$, are

$$M_e(x, y, z) = \begin{bmatrix} -\left(\frac{\partial^2 \Pi_e}{\partial y^2} + \frac{\partial^2 \Pi_e}{\partial z^2}\right) & \frac{\partial^2 \Pi_e}{\partial y \partial x} & \frac{\partial^2 \Pi_e}{\partial z \partial x} \\ \frac{\partial^2 \Pi_e}{\partial x \partial y} & -\left(\frac{\partial^2 \Pi_e}{\partial x^2} + \frac{\partial^2 \Pi_e}{\partial z^2}\right) & \frac{\partial^2 \Pi_e}{\partial z \partial y} \\ \frac{\partial^2 \Pi_e}{\partial x \partial z} & \frac{\partial^2 \Pi_e}{\partial y \partial z} & -\left(\frac{\partial^2 \Pi_e}{\partial x^2} + \frac{\partial^2 \Pi_e}{\partial y^2}\right) \end{bmatrix}_{z=\delta} \qquad (31)$$
and

$$ \mathcal{M}_m(x, y, z) = i\omega \cdot \begin{bmatrix} 0 & -\frac{\partial \Pi_m}{\partial z} & \frac{\partial \Pi_m}{\partial y} \\ \frac{\partial \Pi_m}{\partial z} & 0 & -\frac{\partial \Pi_m}{\partial x} \\ -\frac{\partial \Pi_m}{\partial y} & \frac{\partial \Pi_m}{\partial x} & 0 \end{bmatrix}_{z=\delta} \quad (32) $$

The lower index $z = \delta$ indicates that the derivations in the matrices have to be taken at the value $z = \delta$ with positive or negative value of $\delta$. The authors of [10] have been developed these matrices in circular coordinates $(\rho, \theta)$ and taken a negative value of $\delta$ (focal length). In our case, the general case is considered for our theoretical developments. This is the reason that absolute values of $z$ appear in the matrices $\mathcal{M}_e$ and $\mathcal{M}_m$ in Appendix D. Nevertheless, we will choose the particular case illustrated in Fig. 1 for our illustrations, i.e. $\delta > 0$. Thus Eqs. (80) and (79) together with Eqs. (24) and (23) give the expression of the components of the electromagnetic field illuminating the object then give the transverse components ($W_\delta^x, W_\delta^y$) of the field $W_\delta(\xi, \eta)$ used in the Luneberg integrals Eqs. (2), (3) and (4).

## Theoretical development of the luneberg integrals with the vectorial spherical waves

The Luneberg operator applies in the same way to the magnetic field as to the electric field. We will therefore only deal with the case of the electric field, i.e. $\mathcal{E}_\delta(\xi, \eta)$. From the point of view of the definition of the object, we have previously mentioned that the amplitude function $P(\xi, \eta) = 1 - \text{circ}(\rho/R)$, $R = D/2$. Thus, the amplitude of the wave in the transverse plane just next the object is given by

$$ W_{\delta}^{i}(\xi, \eta) = E_{\delta}^{i}(\xi, \eta) - E_{\delta}^{i}(\xi, \eta) \text{circ}(\rho/R) \quad (33) $$

with $i = x, y$ and $\rho^2 = \xi^2 + \eta^2$. The third component, i.e. Eq. (4) of the field depends on the first two, i.e. Eqs. (2) and (3). Next propagation, in the plane of the CCD sensor localized at the distance $z_l$ from the object, we have

$$ \mathcal{E}_{z_l}^i(x, y) = \mathcal{L}_{z_l}[E_{\delta}^i(\xi, \eta)] - \mathcal{L}_{z_l}[E_{\delta}^i(\xi, \eta) \text{circ}(\rho/R)], \quad (34) $$

The first term of Eq. (34) is the reference vectorial wave. It will be denoted $\mathcal{E}_{z_l}^r = (\mathcal{E}_{z_l}^{rx}, \mathcal{E}_{z_l}^{ry}, \mathcal{E}_{z_l}^{rz})$. Only the two first components $x$ and $y$-axis will be used in the propagation. The second term of Eq. (34) is the object vectorial wave. The following “Luneberg propagation of the electromagnetic fields: the reference vectorial waves” section deals with the first term of Eq. (34) and the second will be treated in “Luneberg propagation of the electromagnetic fields through a pupil” section.

## Luneberg propagation of the electromagnetic fields: the reference vectorial waves

We will be interested here in the propagation of the electric field through the Luneberg operator without the boundary plane of the object. Although the integrals of Luneberg apply on the $x$- and $y$-components and the third one is deduced only from the other two, we let's nevertheless apply the Luneberg operator $\mathcal{L}_{z_l}$ to the matrices $\mathcal{M}_e$ and $\mathcal{M}_m$ which define the vectorial spherical waves in order to obtain the interesting components. By means of Eqs. (8) and (23), we obtain the following result for the reference vectorial field (see Appendix E)

$$ \begin{aligned} \mathcal{E}_{z_l}^r(x, y) &= \mathcal{L}_{z_l}[\mathcal{E}_\delta(\xi, \eta)] \\ &= \mathcal{L}_{z_l}[\mathcal{M}_e(x, y, \delta)]\mathbf{p}_e + \mathcal{L}_{z_l}[\mathcal{M}_m(x, y, \delta)]\mathbf{p}_m, \end{aligned} \quad (35) $$

with

$$ \mathcal{L}_{z_l}[\mathcal{M}_e(x, y, \delta)] = \begin{bmatrix} -\left(\frac{\partial^2 \Pi_e}{\partial y^2} + \frac{\partial^2 \Pi_e}{\partial z^2}\right) & \frac{\partial^2 \Pi_e}{\partial y \partial x} & \text{sign}(\delta) \cdot \frac{\partial^2 \Pi_e}{\partial z \partial x} \\ \frac{\partial^2 \Pi_e}{\partial x \partial y} & -\left(\frac{\partial^2 \Pi_e}{\partial x^2} + \frac{\partial^2 \Pi_e}{\partial z^2}\right) & \text{sign}(\delta) \cdot \frac{\partial^2 \Pi_e}{\partial z \partial y} \\ \text{sign}(\delta) \cdot \frac{\partial^2 \Pi_e}{\partial x \partial z} & \text{sign}(\delta) \cdot \frac{\partial^2 \Pi_e}{\partial y \partial z} & -\left(\frac{\partial^2 \Pi_e}{\partial x^2} + \frac{\partial^2 \Pi_e}{\partial y^2}\right) \end{bmatrix}_{z=|\delta|+z_l}, \quad (36) $$

and

$$ \mathcal{L}_{z_l}[\mathcal{M}_m(x, y, \delta)] = i\omega \cdot \begin{bmatrix} 0 & -\text{sign}(\delta) \cdot \frac{\partial \Pi_m}{\partial z} & \frac{\partial \Pi_m}{\partial y} \\ \text{sign}(\delta) \cdot \frac{\partial \Pi_m}{\partial z} & 0 & -\frac{\partial \Pi_m}{\partial x} \\ -\frac{\partial \Pi_m}{\partial y} & \frac{\partial \Pi_m}{\partial x} & 0 \end{bmatrix}_{z=|\delta|+z_l}. \quad (37) $$

The function $\text{sign}(x) = x/|x|$ comes from the $z$-partial derivations in the matrices and physically, allows us to take into account the sign reversal on both sides of the focused point of the field, at the center of the sphere. The solutions obtained in Eqs. (36) and (37) are the same expressions as in Eqs. (31) and (32). The theoretical developments of [10] can be reinvested by adding the sign function and by shifting along the $z$ axis of the quantity $|\delta| + z_l$. It is clear that when $\delta > 0$ then

$$ \mathcal{E}_{z_l}^r(x, y) = \mathcal{M}_e(x, y, \delta + z_l) \cdot \mathbf{p}_e + \mathcal{M}_m(x, y, \delta + z_l) \mathbf{p}_m, \quad (38) $$

Then, the Luneberg transformation over $z_l$ of the matrices fields $\mathcal{M}_e$ and $\mathcal{M}_m$ at the location $\delta$ is the matrices field at $\delta + z_l$. The choice of the potential in the resolution of the Maxwell's equations is important from Luneberg operator point of view. Only the Hertz potentials allows us to obtain these properties about the propagation. This results are illustrated in Fig. 3 about the components $\mathcal{E}_{z_l}^{rx}$
<figure><img src="image_8.png" /><figcaption>Fig. 5 Double exponential numerical evaluation of the Zernike coefficients for the transverse components of the vectorial field $\tilde{E}_{\delta}^x$ and $\tilde{E}_{\delta}^y$ with $\delta = 5\lambda, D = 7\lambda$</figcaption></figure>

and $E_{z_l}^{rz}$ of the electric field where $\Re$ and $\Im$ represent the real and imaginary parts. These two components are simulated by means of the convolution product in Eqs. (8,10) and by means of the matrix expression in Eq. (38). The initial electric field is localized at $\delta = 5\lambda$ and it propagates over $z_l = 3\lambda$. The polarizations are $p_e = (1, 0, 0)$ and $p_m = (0, 1, 0)$.

The good agreement between these two results allows us to extract easily from Eq. (38), the transverse components next the propagation which will be used in Eq. (34) without additional theoretical calculus.

### Luneberg propagation of the electromagnetic fields through a pupil

**Transverse components** This subsection deals with of the Luneberg propagation of the transverse components of the vectorial field diffracted by the pupil, i.e. Eqs. (2) and (3). The product of the field by the pupil will be treated as a generalized pupil function. Thus, to treat the second term of Eq. (34), it is now commonplace to expand the generalized pupil $E_{\delta}^i(\xi, \eta)$ circ($\rho/R$) as

$$E_{\delta}^{i}(\xi, \eta) \text{circ}(\rho/R) = \sum_{n,m} \gamma_{ni}^{m} \cdot Z_{n}^{m}\left(\frac{\rho}{R}, \theta\right), \quad (39)$$

in which

$$Z_{n}^{m}(\varrho, \theta) = R_{n}^{\vert m \vert}(\varrho) \cdot e^{i m \theta}, \quad 0 \le \varrho \le 1, \quad 0 \le \theta \le 2\pi, \quad (40)$$

and where $n$ and $m$ in the summation in Eq. (39) and in Eq. (40) are integers such that $n - |m|$ is even and non-negative. The Zernike coefficients $\gamma_{ni}^m$ are obtained by using the orthogonality of the Zernike circle polynomials $Z_n^m$, so that

$$\begin{aligned} \gamma_{ni}^{m} &= \frac{n+1}{\pi R^2} \int_{0}^{R} \int_{0}^{2\pi} E_{\delta}^{i}(\xi, \eta) \\ & \qquad \circ \text{circ}(\rho/R) \cdot \frac{Z_{n}^{m}\left(\frac{\rho}{R}, \theta\right)}{\overline{Z_{n}^{m}(\varrho, \theta)}} \rho \, d\rho \, d\theta \\ &= \frac{(n+1)}{\pi} \int_{0}^{1} \int_{0}^{2\pi} E_{\delta}^{i}(R\varrho \cos\theta, R\varrho \sin\theta) \\ & \qquad \cdot \overline{Z_{n}^{m}(\varrho, \theta)} \, \varrho \, d\varrho \, d\theta, \end{aligned} \quad (41)$$

with $\varrho = \rho/R$ and the overhead bar indicates complex conjugation. The Zernike coefficient $\gamma_{nx}^m$ for $i = x$ is linked to the component field $E_{\delta}^x(\xi, \eta)$ and $\gamma_{ny}^m$ for $i = y$ to $E_{\delta}^y(\xi, \eta)$. Note that $E_{\delta}^i(\xi, \eta) = E_{\delta}^i(R\varrho \cos\theta, R\varrho \sin\theta) = \tilde{E}_{\delta}^i(\varrho, \theta)$. The Zernike coefficients do not have a closed form. Consequently, the Zernike coefficients $\gamma_{ni}^m$ in Eq. (39) should be numerically evaluated, for instance by
<figure><img src="image_9.png" /><figcaption>Fig. 6 Comparison of the Zernike coefficients between the integral evaluation and the double exponential method of $\tilde{E}_{\delta}^{x}(\xi, \theta)$ and $\tilde{E}_{\delta}^{y}(\xi, \theta)$ with $\delta = 5\lambda, D = 7\lambda$ and the polarization $\mathbf{p}_e = (1, 0, 0), \mathbf{p}_m = (0, 1, 0)$. For the x-component of the electric field $n = 0(2)20, m = 0$ and for the y-component $n = 2(2)20, m = 2$</figcaption></figure>

means of the double exponential formula (DE) of numerical integration [12, 23]. To illustrate the numerical evaluation of the Zernike coefficients, we have been considered the vectorial electromagnetic field of wavelength defined in Eq. (23) by extracting the all components limited by pupil of diameter $D = 7\lambda$ localized at $\delta = 5\lambda$. The polarization vectors $\mathbf{p}_e = (1, 0, 0)$ and $\mathbf{p}_m = (0, 1, 0)$ which are linear polarizations. All the electric field components are represented in Fig. 4 but only the first components of the electric field will be considered in the paper. To treat the magnetic field, the same development process shall be carried out.

The illustrations in Fig. 5 give us the results of the evaluation of the Zernike coefficients by the DE method of the vectorial electric field.

The real part, denoted $\Re$, and imaginary part, denoted $\Im$, of the Zernike coefficients are represented. The illustrations themselves in Fig. 6 represent the comparison between direct integral evaluation from Eq. (41) and the DE method.

As we can see, the decomposition of the $\tilde{E}_{\delta}^{x}(\xi, \eta)$ field gives us dominant Zernike coefficients for $m = 0$. A few azimuth coefficients are linked to the small astigmatism of the field. The interpretation is that the $\tilde{E}_{\delta}^{x}(\xi, \eta)$ component has a very slight rotational asymmetry. On the other hand, the decomposition of the $\tilde{E}_{\delta}^{y}(\xi, \eta)$ component gives us dominant Zernike coefficients for $m \neq 0$. This implies rotational asymmetry. From the spectral point of view, the combination of the Eqs. (8) and (68) and by noting that from [16], Eq. (5) on p. 2

$$ \mathcal{F}[Z_n^m(\frac{\rho}{R}, \theta)] = 2\pi R^2 (-i)^m (-1)^{\frac{n-m}{2}} \cdot \frac{J_{n+1}(2\pi\sigma)}{2\pi\sigma} \exp(im\varphi), \quad \sigma \ge 0, \ 0 \le \varphi \le 2\pi \quad (42) $$

with $\sigma = R\sqrt{u^2 + v^2}$ and $\varphi = \arg(u + iv)$ allows us to obtain the spectrum, denoted $\hat{\mathcal{E}}_{z_l}^i(\sigma, \varphi)$, of the second term of Eq. (34) as
$$ \hat{\mathcal{E}}_{z_l}^i = \frac{\pi D^2}{2} \mathrm{sign}(z_l) \sum_{n,m} \gamma_{ni}^m (-i)^m (-1)^{\frac{n-m}{2}} \frac{J_{n+1}(2\pi\sigma)}{2\pi\sigma} e^{im\varphi} \\ \times \exp\left(i \frac{4\pi |z_l|}{D} \sqrt{\frac{D^2}{4\lambda^2} - \sigma^2}\right) \quad (43) $$

where

$$ \left(\frac{D^2}{4\lambda^2} - \sigma^2\right)^{1/2} = \begin{cases} \sqrt{\frac{D^2}{4\lambda^2} - \sigma^2} & \text{if } 0 \le \sigma \le R/\lambda, \\ i\sqrt{\sigma^2 - \frac{D^2}{4\lambda^2}} & \text{if } R/\lambda \le \sigma < \infty. \end{cases} \quad (44) $$

For $\sigma \ge R/\lambda$, this corresponds to evanescent waves in Weyl's integral where the evanescent wave attenuate exponentially with increasing $z_l$. This is the basis for the usual argument that evanescent plane waves can be neglected sufficiently far away from the $z_l = 0$ plane. The numerical simulation of Eq. (43) could be realized by using the recurrence relations of the Bessel function as solution when

$\sigma = 0$. The inverse Fourier transform with a circular symmetry (Hankel transform) is

$$ \begin{aligned} \mathcal{E}_{z_l}^i(s, \theta) &= \mathcal{F}^{-1} \left[ \hat{\mathcal{E}}_{z_l}^i(\sigma, \varphi) \right] \\ &= \int_0^{+\infty} \int_0^{2\pi} \hat{\mathcal{E}}_{z_l}^i(\sigma, \varphi) \exp[i2\pi\sigma s \cos(\varphi - \theta)] \sigma d\sigma d\varphi. \end{aligned} \quad (45) $$

By noting that

$$ \begin{aligned} & \int_0^{2\pi} \exp(i m \varphi) \cdot \exp(i 2\pi \sigma s \cos(\varphi - \theta)) d\varphi \\ &= 2\pi i^m J_{-m}(2\pi s \sigma) \exp(-i m \theta), \end{aligned} \quad (46) $$

with $J_m$ the Bessel function of the first kind and of order $m$, the expression of the transverse components of the vectorial field are defined by

$$ \begin{aligned} \mathcal{E}_{z_l}^i(s, \theta) &= 2\pi\omega \mathrm{sign}(z_l) \sum_{n,m} \gamma_{ni}^m (-1)^{\frac{n-m}{2}} \exp(i m \theta) \\ & \quad \int_0^{+\infty} \exp(i\beta\sqrt{1-X^2}) J_{n+1}(2\pi\omega X) J_{|m|}(2\pi\omega sX) dX, \end{aligned} \quad (47) $$

<figure><img src="image_10.png" /><figcaption>**Fig. 7** Comparison of $\mathcal{I}_{nm}(\beta, \omega, s)$ between Eq. (51) and Eq. (52-53) and comparison of $\mathcal{J}_{nmk}(\beta, \omega, s)$ between Eq. (57) and Eq. (58-59) for $\omega = 3.5$, $\beta = 40.84$, $n = 4$, $m = 2$ and $k = 1$</figcaption></figure>
in which $\beta = 2\pi \frac{|z|}{\lambda}$ and $\omega = \frac{D}{2\lambda}$. The integral in (47) will be evaluated by considering the ranges $[0, 1]$ and $[1, +\infty)$ separately.

**Longitudinal component** This paragraph deals with the Luneberg propagation of the longitudinal component of the vectorial field diffracted by the pupil from Eq. (4). Firstly, by combining Eqs. (4), (39) and considering the Fourier transform in circular coordinates, the expression of the spectrum, denoted $\hat{\mathcal{E}}_{z_l}^z(\sigma, \varphi)$, of $\mathcal{E}_{z_l}^z(x, y)$ with the same processes of calculus as previously, we obtain

$$ \hat{\mathcal{E}}_{z_l}^z(\sigma, \varphi) = -\frac{D^2}{4} \sum_{n,m} (-i)^m (-1)^{\frac{n-m}{2}} \left( \gamma_{nx}^m \cos \varphi + \gamma_{ny}^m \sin \varphi \right) \exp(im\varphi) \frac{\exp\left(i \frac{4\pi |z_l|}{D} \sqrt{\frac{D^2}{4\lambda^2} - \sigma^2}\right)}{\sqrt{\frac{D^2}{4\lambda^2} - \sigma^2}} \quad (48) $$

By noting that

$$ \int_0^{2\pi} \left( \gamma_{nx}^m \cos \varphi + \gamma_{ny}^m \sin \varphi \right) e^{im\varphi} \exp(i2\pi\sigma s \cos(\theta - \varphi)) d\varphi \\ = \pi (i)^m \sum_{k=0,1} \left( \gamma_{ny}^m + i^{2k+1} \gamma_{nx}^m \right) J_{1+m-2k}(2\pi s \sigma) e^{i(1+m-2k)\theta}, \quad (49) $$

the expression of the longitudinal component in the spatial domain of the vectorial field is defined by

$$ \mathcal{E}_{z_l}^z(s, \theta) = \pi\omega \sum_{nmk} (-1)^{\frac{n-m}{2}+1} (i)^{1+m-2k-|1+m-2k|} \left( \gamma_{ny}^m + i^{2k+1} \gamma_{nx}^m \right) e^{i(1+m-2k)\theta} \int_0^{+\infty} \frac{\exp\left(i\beta\sqrt{1-X^2}\right)}{\sqrt{1-X^2}} J_{n+1}(2\pi\omega X) J_{|1+m-2k|(2\pi\omega sX)} dX. \quad (50) $$

As the same way of the integral in Eq. (47), the evaluation of the remaining integral in Eq. (50) will be done by considering the ranges $[0, 1]$ and $[1, +\infty)$ separately.

**Evaluation of the integrals of the Bessel product** The problem of the study of the Bessel function $J_\nu(x)$ in which the order $\nu$ and argument $x$ are positive, and then the product of the Bessel functions under the integral, are usually linked to the ratio $x/\nu$ that is less than, nearly equal to, or greater than unity. In our case, the Bessel functions depend, on the one hand, on the argument $\omega$ and, on the other hand, on the radial argument $s$ and of $\omega$ for the other. It is clear that when $s$ and $\omega$ are large, the power series expansion of Bessel function converges slowly. Results have been obtained by the author in [28], Appendix E, pp. 1278-1280 on the integration of a product of two Bessel

functions. However, the experimental context means that variables such as the radial coordinate of the Bessel functions are less than or equal to unity, limited by a pupil. This context does not correspond to ours because the variables are spatial and therefore not limited by the unit. It is therefore necessary to choose another approach to evaluate the integral involving these two Bessel functions. The chosen approach here is to apply a Zernike expansion of the two Bessel functions because the excellent convergence and on the exponential function. In this manner, the diameters of the object will be able to be increased compared to the wavelength. Denote by $\mathcal{I}_{nm}(\beta, \omega, s)$, the integral in Eq. (47) so that

$$ \mathcal{I}_{nm}(\beta, \omega, s) = \int_{0}^{+\infty} \exp\left(i\beta\sqrt{1-X^2}\right) J_{n+1}(2\pi\omega X) J_{|m|}(2\pi\omega sX) dX. \quad (51) $$

The evaluation of this integral is given in (F) and yields

$$ \mathcal{I}_{nm}(\beta, \omega, s) = \frac{\pi \omega}{(n+1)} \sum_{t,k} \left[ G_{tk}^{n,|m|}(\beta, \omega) + G_{tk}^{n+2,|m|}(\beta, \omega) \right] \frac{J_{|m|+2t+1}(2\pi \omega s)}{2\pi \omega s} + \frac{1}{\beta^2} J_{n+1}(2\pi \omega) J_{|m|}(2\pi \omega s), \quad (52) $$

with $\beta \gg 0$ and $0 \le k < \infty$, $0 \le t < \infty$ and

$$ \mathcal{G}_{tk}^{N,M}(\beta, \omega) = 4(-1)^{k+t}(N+2k+1)(M+2t+1) \frac{J_{N+2k+1}(2\pi\omega)}{2\pi\omega} \sum_{q=0}^{\frac{N+M}{2}+k+t} \frac{\mathcal{L}C_{tkq}^{M,N,0}}{2(2q+1)} \mathcal{B}_q(\beta). \quad (53) $$

The expression of $\mathcal{B}_q(\beta)$ is in Eq. (95), Appendix (F) and linked to the $R_{2r}^0$-expansion. The linearization coefficients $\mathcal{L}C_{tkq}^{M,N,0}$ are defined in Eq. (124) in Appendix (F). It can be computed by recursion, the first two terms of Eq. (52) are linked to the integration of $\mathcal{I}_{nm}(\beta, \omega, s)$ in the range $[0, 1]$ and the third is linked to the evanescent wave in the ranges $[1, +\infty)$. This third term is evaluated by means of the steepest descent method. The Zernike expansion of the Bessel function and a linearisation of the product of radial polynomials are used to obtain the $\mathcal{L}C_{tkq}^{M,N,0}$ coefficients. The combination of these steps with the orthogonality of the radial polynomials gives us $G_{tk}^{N,M}$. Notwithstanding the condition on $\beta$, if $\beta = 0$, the Zernike expansion of the exponential in Eq. (51) is not necessary. As (see Eq. (97))

$$ R_{M+2t}^{M}(X) R_{N+2k}^{N}(X) = \sum_{q=0}^{\frac{N+M}{2}+k+t} \mathcal{L}C_{tkq}^{M,N,0} \cdot R_{2q}^{0}(X), \quad (54) $$

thus by orthogonality of the Zernike polynomials

$$ \int_0^1 R_{N+2k}^N(X) R_{M+2t}^M(X) X dX = \frac{1}{2} \mathcal{L}C_{tk,q=0}^{M,N,0} \quad (55) $$
An another point of view from [24, 36], the outcome of the integral in Eq. (51) for $\beta = 0$ equals to

$$ \mathcal{I}_{nm}(0, \omega, s) = \frac{1}{2\pi\omega}(-1)^{\frac{n-|m|}{2}} R_{n}^{|m|}(s), \quad (56) $$

for $0 \le s < 1$ and $0 \le s > 1$. For the second integral in Eq. (50), denoted as $\mathcal{J}_{nmk}(\beta, \omega, s)$, so that

$$ \mathcal{J}_{nmk}(\beta, \omega, s) = \int_{0}^{+\infty} \frac{\exp(i\beta\sqrt{1-X^2})}{\sqrt{1-X^2}} J_{n+1}(2\pi\omega X) J_{|1+m-2k|}(2\pi\omega sX) X \, dX, \quad (57) $$

this type of integral appears in numerous papers which have treated the problems of sound radiation from vibrating circular plates in particular King [18] in which he has exhibited some useful integral expressions for the acoustic pressure. The case of integral like in Eq. (57) arise when one considers generalization of King's integral to the non-radially symmetric pistons. The authors of [1, 2] have been proposed a new unified method to resolve King's integral where the product of two Bessel functions are involved. In this context, the unified solution are based on the case of different orders for the Bessel functions but with of same argument. Here, with the previous results, it is now possible to propose a new solution for this difficult issue. Consequently, its evaluation in Appendix (F) gives us the result

$$ \mathcal{J}_{nmk}(\beta, \omega, s) = \sum_{t,l} \mathcal{H}_{tl}^{n+1, |1+m-2k|}(\beta, \omega) \frac{J_{|1+m-2k|+2t+1}(2\pi\omega s)}{2\pi\omega s} \\ - i \frac{1}{\beta} J_{n+1}(2\pi\omega s) J_{|1+m-2k|+2t+1}(2\pi\omega s) \quad (58) $$

with $l = 0, 1, \dots, t = 0, 1, \dots, k = 0, 1$ and

$$ \mathcal{H}_{tl}^{N,M}(\beta, \omega) = 4(-1)^{t+l}(N+2l+1)(M+2t+1) \frac{J_{N+2l+1}(2\pi\omega)}{2\pi\omega} \sum_{q=0}^{\frac{N+M}{2}+l+t} \frac{\mathcal{L}\mathcal{C}_{tlq}^{M,N,0}}{2(2q+1)} \mathcal{A}_q(\beta) \quad (59) $$

with $\mathcal{A}_q(\beta)$ defined in Eq. (126). We will now illustrate and compare the results. The numerical calculations, illustrated in Fig. 7, are based on the comparison of Eqs (51) and (52-53) as well as Eqs (57) and (58-59). The truncations of the series are based on the work in [33]. These simulations are performed in the case where the parameters $\omega = 3.5$ and $\beta = 40.84$ that correspond to the optical parameter $\lambda = 633$nm, $D = 7\lambda$, $z_l = 6.5\lambda$. We choose as degree $n = 4$, azimuthal $m = 2$ and $k = 1$ in Eqs. (57) and (51). As we can see in Fig. (7), the radial coordinate takes a maximal value $s = 4$. The maximum value of $s$ impacts the computation time and even the possibility of performing the computation. Indeed, the choice of the approach to compute the integrals in Eqs. (51) and (57) allows or not to realize the computation. Compared to the approach of many linearisation of products of radial polynomials by Clebsch-Gordan coefficients, the approach by means of Eq. (54) and with using the recursion allows to substantially decrease the computation time, typically 116ms for Eq. (52) and 141ms for Eq. (58), and to increase the maximal value of $s$. With $s = 2.5$, the approach by many linearisation is not optimal in term of computation time because the time is around 198s for (51) and 90s for (57).

The Fig. 7 demonstrate that a good agreement is obtained between the functions $\mathcal{J}_{nmk}(\beta, \omega, s)$, $\mathcal{I}_{nm}(\beta, \omega, s)$ and their integral forms evaluated under Mathematica software. Thus, this concordance will allow us to simulate a few holograms of opaque particles that we consider to be common in holography like opaque disks.

<figure><img src="image_11.png" /><figcaption>Fig. 8 Holograms for different polarizations: linear, transverse electric and transverse magnetic with $\lambda = 0.633\mu m$, $\delta = 5\lambda$, $D = 7\lambda$, $z_l = 6.5\lambda$, $\omega = 3.5$ and $\beta = 40.84$</figcaption></figure>
## Simulation of the holograms

In this section, we will show the intensity distribution of the hologram in the quadratic sensor plane, denoted $I(s, \theta)$. Its expression is as follows

$$I(s, \theta) = |\mathcal{E}_{\delta+z_l}^x|^2 + |\mathcal{E}_{\delta+z_l}^y|^2 + |\mathcal{E}_{\delta+z_l}^z|^2. \quad (60)$$

The intensity distribution has been simulated for various polarization states of the incident vectorial fields in the object. The authors of [10] studied the following different states of polarization: the linear polarisation ($p_e = (1, 0, 0)$ and $p_m = (0, 1, 0)$), the transverse electric ($p_e = 0, p_m = (0, 0, 1)$) and the transverse magnetic ($p_e = (0, 0, 1), p_m = 0$) mode of the incident light. For example, in the case of the transverse electric mode, from Eqs. (23), (31), (32) and (80), the incident electric field $E_\delta(\xi, \eta)$ is defined as

$$E_{\delta} = i\omega \cdot \begin{bmatrix} \frac{\partial \Pi_m}{\partial y} \\ -\frac{\partial \Pi_m}{\partial x} \\ 0 \end{bmatrix}_{z=\delta} = \frac{\exp(ikr)}{2r^5} \begin{bmatrix} -\rho S(\rho, \delta) \sin(\theta) \\ \rho S(\rho, \delta) \cos(\theta) \\ 0 \end{bmatrix} := \begin{bmatrix} E_{\delta}^x \\ E_{\delta}^y \\ 0 \end{bmatrix}. \quad (61)$$

The propagation of the field is along the z-axis as mentioned by the Luneberg integrals. The components of the field are perpendicular to the optical axis. These components illuminate the object and in the same time are the reference wave. By means of the Zernike decomposition, the Lunberg's integrals and the results obtained and presented previously, the intensity distribution of the hologram can be evaluated.

As illustrated in Fig. 8, the structure of the holograms are versus the polarisation states in near field compared to the far field. Once again, a comparison is made with the use of the fast Fourier Transformation (see Fig. 9) tool to check adequacy in the simple manner but now with the know the functions which compose the results of the simulations.

Of course, it is possible to use the Jones matrices too. The x and y-components of the vectorial electromagnetic field are determined versus $p_e$ and $p_m$ as previously seen. From this example, if we place a linear polariser with axis of transmission angle $\theta$ from the horizontal in the plane before the object, at the output of the polariser, we will have

$$\begin{bmatrix} E_{\delta}^{x,o} \\ E_{\delta}^{y,o} \end{bmatrix} = \begin{bmatrix} \cos^2 \theta & \frac{1}{2} \sin(2\theta) \\ \frac{1}{2} \sin(2\theta) & \sin^2 \theta \end{bmatrix} \begin{bmatrix} E_{\delta}^{x} \\ E_{\delta}^{y} \end{bmatrix}, \quad (62)$$

where $(E_{\delta}^{x,o}, E_{\delta}^{y,o})$ are the components illuminating the opaque disk and which will propagate just next the object. Then, the same process as seen in this paper is used to simulate the hologram represented in Fig. 10.

## Conclusion

We have recalled the definitions of the Luneberg integrals and the conditions of the uniqueness of the solution of the boundary value problem of the equation $\Delta u + k^2 u = 0$ for a plane boundary. These conditions must be taken into account each time we want to use a definition of a wave function. We clearly mentioned that the solution of the Hertz potentials allows us to obtain a definition of the vectorial waves that satisfy the Luneberg conditions. Moreover, the Luneberg integrals applied to the free propagation of the vectorial electromagnetic wave in matrix form give us the same matrix form to define the vectorial electromagnetic wave after the propagation. This central point of the free propagation throughout the Luneberg integrals allow us to simulate the holograms because, in the case of the DIH, the reference wave is necessary to be associated with the object wave. Like this, the reference wave conserves the vectorial nature. Concerning the object wave, the approach by means of the Zernike polynomials has been proposed to evaluate the propagation integrals where Bessel functions products appear. Linear and connection coefficients are defined and recursion relations are proposed to obtain the results for

<figure><img src="image_12.png" /></figure>
the Luneberg integrals of the propagation. The recursion steps are important to ensure stable, quick and accurate computations. The knowledge of the model will allow us in the future to interpret the holographic results. Furthermore, taking into account the vectorial property of the light allows to access to new analyses by means of the DIH such as the birefringent objects.

However, it will be necessary to adapt the reconstruction operators of the image of the object to access the usual informations such as the size and 3D location. Indeed, the reconstructions are currently carried out by means of Fresnel integral or by means of the fractional Fourier transform. These operators are adapted to the case of the plane wave while here we are spherical waves. A metrology on the indices as well as the interpretation of the reconstructed holograms will require a theoretical calculation reference standard.

Comparisons between theoretical and numerical results are illustrated all the paper. Finally, we illustrate digital in-line holograms with different polarization states and we use briefly the Jones matrix. Matrices that can be introduced into the mathematical definition of the object.

## Appendix A: Longitudinal component

About longitudinal component, the vector field $\mathcal{E}_{z_l} = (\mathcal{E}_{z_l}^x, \mathcal{E}_{z_l}^y, \mathcal{E}_{z_l}^z)$ must satisfy the Maxwell's equations in the vacuum, ie div $\mathcal{E}_{z_l} = 0$, then

$$ \frac{\partial \mathcal{E}_{z_l}^z}{\partial z} = - \left( \frac{\partial \mathcal{E}_{z_l}^x}{\partial x} + \frac{\partial \mathcal{E}_{z_l}^y}{\partial y} \right). \qquad (63) $$

By introducing Eqs. (2)and (3) into Eq. (63) and by inverting the partial derivations, we have by integration

$$ \mathcal{E}_{z_l}^z = \frac{1}{2\pi} \iint \left[ \mathbf{W}_{\delta}^x(\xi, \eta) \frac{\partial}{\partial x} \left( \frac{\exp(ikr)}{r} \right) + \mathbf{W}_{\delta}^y(\xi, \eta) \frac{\partial}{\partial y} \left( \frac{\exp(ikr)}{r} \right) \right] d\xi d\eta. \qquad (64) $$

With $r^2 = (x - \xi)^2 + (y - \eta)^2 + z_l^2$ and

$$ \frac{\partial}{\partial x} \left( \frac{\exp(ikr)}{r} \right) = (x - \xi) \left( \frac{ikr - 1}{r^3} \right) \exp(ikr), \qquad (65) $$

Eq. (64) can be re-written in the form of convolution product such as

$$ \mathcal{E}_{z_l}^z = - \mathbf{W}_{\delta}^x(x, y) * \left[ \frac{x}{z_l} \cdot L_{z_l}(x, y) \right] - \mathbf{W}_{\delta}^y(x, y) * \left[ \frac{y}{z_l} \cdot L_{z_l}(x, y) \right]. \qquad (66) $$

## Appendix B: Boundary values

we start by the two important relations: firstly, the Fourier transform of the spherical wave which is a derivation of

Weyl's plane wave expansion of the spherical wave [9, 34]

$$ \mathcal{F}\left[\frac{\exp(ikr)}{r}\right](u, v) = i2\pi \frac{\exp\left(\frac{i|z|\sqrt{k^2 - 4\pi^2(u^2 + v^2)}}{\sqrt{k^2 - 4\pi^2(u^2 + v^2)}}\right)}{\sqrt{k^2 - 4\pi^2(u^2 + v^2)}}, \qquad (67) $$

where $r^2 = x^2 + y^2 + z^2$ and secondly the Fourier transform of the Luneberg kernel

$$ \begin{aligned} \mathcal{F}[L_z(x, y)](u, v) &= \mathcal{F}\left[-\frac{1}{2\pi}\frac{\partial}{\partial z}\left(\frac{\exp(ikr)}{r}\right)\right](u, v) \\ &= -\frac{1}{2\pi}\frac{\partial}{\partial z}\mathcal{F}\left[\frac{\exp(ikr)}{r}\right] \\ &= \operatorname{sign}(z)\exp\left(|z|\sqrt{k^2 - 4\pi^2(u^2 + v^2)}\right) \end{aligned} \qquad (68) $$

where $\frac{\partial}{\partial z}|_{z|} = \operatorname{sign}(z)$. From Eq. (8) in the spectral domain, we have

$$ \begin{aligned} \hat{\mathcal{E}}_{z_l}^i(u, v) &= \mathcal{F}[W_{\delta}^i(x, y) * L_{z_l}(x, y)](u, v) \\ &= \hat{W}_{\delta}^i(u, v) \cdot \operatorname{sign}(z_l) \exp\left[|z_l|\sqrt{k^2 - 4\pi^2(u^2 + v^2)}\right], \end{aligned} \qquad (69) $$

where $\operatorname{sign}(x) = x/|x|$. The limit case give

$$ \lim_{z_l \to 0^+} \hat{\mathcal{E}}_{z_l}^i(u, v) = W_{\delta}^i(u, v). \qquad (70) $$

This is Eq. (7) by an inverse Fourier transformation.

## Appendix C: Conditions about the far field approximation of the luneberg's kernel

As previously written in Eq.(9), the Luneberg's kernel is the following

$$ L_z(x, y) = \frac{1}{2\pi z} \left( \frac{1}{z\left(1 + \frac{x^2+y^2}{z^2}\right)^{3/2}} - i \frac{2\pi}{\lambda} \frac{1}{1 + \frac{x^2+y^2}{z^2}} \right) \times \exp \left[ i \frac{2\pi}{\lambda} z \left( 1 + \frac{x^2+y^2}{z^2} \right)^{1/2} \right]. \qquad (71) $$

With the first order Taylor series expansion

$$ (1 + X)^\alpha = 1 + \alpha X + \mathcal{O}[X^2], \qquad (72) $$

where $X = (x^2 + y^2)/z^2$, Eq. (71) becomes

$$ L_z(x, y) = \frac{1}{2\pi z} \left( \frac{x^2 + y^2}{z^2} \left[ i \frac{2\pi}{\lambda} - \frac{3}{2z} \right] + \left( \frac{1}{z} - i \frac{2\pi}{\lambda} \right) \right) \times \exp[ikz] \exp \left[ i \frac{\pi}{\lambda z} (x^2 + y^2) \right]. \qquad (73) $$

Now, if

$$ \left| i \frac{2\pi}{\lambda} \right| \gg \frac{3}{2z}, \quad \text{and} \quad \left| -i \frac{2\pi}{\lambda} \right| \gg \frac{1}{z} \qquad (74) $$

then the Luneberg's kernel equal to Fresnel's kernel with the most restrictive conditions

$$ z \gg \frac{3}{2k}, \quad \text{and} \quad z \gg r_{\max}, \qquad (75) $$

with $r_{\max} = \max(\sqrt{x^2 + y^2})$ is the maximal transverse radius. This second condition comes from the series expansion where $X \ll 1$ and corresponds to the paraxial
<figure><img src="image_13.png" /><figcaption>Fig. 10 Holograms with linear polarizer $\theta = \frac{\pi}{4}$, $\lambda = 0.633\mu m$, $\delta = 5\lambda$, $D = 7\lambda$, $z_l = 6.5\lambda$, $\omega = 3.5$ and $\beta = 40.84$</figcaption></figure>

approximation, for small angles. To introduce the numerical parameter such as the number of sampling N which is the number of pixel of the CCD camera, it is necessary to introduce the instantaneous frequency, denoted $f_i$ such as

$$f_i(x) = \frac{1}{2\pi} \frac{\partial\phi(x)}{\partial x}, \qquad (76)$$

where the phase $\phi$ is the quadratic phase of the Fresnel's kernel. Note that, it is possible to refer to this quadratic phase because we consider that the most restrictive conditions in Eq. (75) are satisfied. Then, the maximal frequency, denoted $f_{imax}$ is reached when $x = \Delta x$, with $\Delta x$ is the spatial support of the image. Then $f_{imax} = \frac{\Delta x}{\lambda z}$. The sampling frequency, denoted $f_e$, must exceed twice the maximal frequency. Consequently, by written that $f_e \ge 2f_{imax}$ and by considering that the support $\Delta x = N\delta_e$ with $\delta_e = 1/f_e$, we obtain

$$\Delta x \le \left( \frac{1}{2} N \lambda z \right)^{1/2}. \qquad (77)$$

If we suppose that the support $\Delta y$ along y-axis is equal to $\Delta x$, then $r_{max} = \frac{1}{2}(\Delta x^2 + \Delta y^2)^{1/2}$ and with the condition in Eq. (75) combined with the Eq. (77), we get

$$z \gg \max \left( \frac{3\lambda}{4\pi}, \frac{1}{4} N \lambda \right) \qquad (78)$$

Finally, It is unlikely that the first term is greater than the second. Consequently, we can consider only the second term.

## Appendix D: Definition of the functions of the matrices $M_e$ and $M_m$

From [10], Eqs. (13) and (14) on p. 734-735 we have, the expression of the matrices $M_e$ and $M_m$ are

$$M_e(x, y, z) = \frac{\exp(ikr)}{2r^5} \times \begin{bmatrix} \mathcal{P}_1(\rho, z) + \rho^2 \mathcal{P}_2(\rho, z) \cos(2\theta) & \rho^2 \mathcal{P}_2(\rho, z) \sin(2\theta) & 2\rho \mathcal{P}_2(\rho, z) \cos(\theta) \\ \rho^2 \mathcal{P}_2(\rho, z) \sin(2\theta) & \mathcal{P}_1(\rho, z) - \rho^2 \mathcal{P}_2(\rho, z) \cos(2\theta) & 2\rho \mathcal{P}_2(\rho, z) \sin(\theta) \\ 2\rho \mathcal{P}_2(\rho, z) \cos(\theta) & 2\rho \mathcal{P}_2(\rho, z) \sin(\theta) & \mathcal{Q}(\rho, z) \end{bmatrix}, \qquad (79)$$

where $x + iy = \rho \exp(i\theta)$, $r^2 = \rho^2 + z^2$ and

$$M_m(x, y, z) = \frac{\exp(ikr)}{2r^5} \begin{bmatrix} 0 & zS(\rho, z) & -\rho S(\rho, z) \sin(\theta) \\ -zS(\rho, z) & 0 & \rho S(\rho, z) \cos(\theta) \\ \rho S(\rho, z) \sin(\theta) & -\rho S(\rho, z) \cos(\theta) & 0 \end{bmatrix}. \qquad (80)$$
with

$$
\begin{align}
\mathcal{P}_1(\rho, z) &= 2k^2 z^4 + \rho^2 \left(1 + k^2 \rho^2 - ik\sqrt{\rho^2 + z^2}\right) \nonumber \\
&\quad + z^2 \left(-2 + 3k^2 \rho^2 + i2k\sqrt{\rho^2 + z^2}\right), \\
\mathcal{P}_2(\rho, z) &= 3 - i3k\sqrt{\rho^2 + z^2} - k^2 (\rho^2 + z^2), \nonumber \\
\mathcal{Q}(\rho, z) &= 4 (\rho^2 + z^2) - i4k (\rho^2 + z^2)^{3/2} \nonumber \\
&\quad + \rho^2 \left(-6 + i6k\sqrt{\rho^2 + z^2} + 2k^2 (\rho^2 + z^2)\right), \tag{81} \\
\mathcal{S}(\rho, z) &= 2k (\rho^2 + z^2) (i + k\sqrt{\rho^2 + z^2}). \nonumber
\end{align}
$$

Appendix E: Free luneberg propagation of the electromagnetic field

In this appendix, we will only illustrate the development of
one term of the matrices $\mathcal{M}_e$ and $\mathcal{M}_m$. Before, let us note
that the Luneberg transforms of the propagation matrices
$\mathcal{M}_e$ and $\mathcal{M}_m$ is the matrix of the Luneberg transforms of
each term of $\mathcal{M}_e$ and $\mathcal{M}_m$

$$
\mathcal{L}_{z_l}[\mathcal{M}_e(x, y, \delta)] =
\begin{bmatrix}
\mathcal{L}_{z_l} \left[ -\left( \frac{\partial^2 \Pi_e}{\partial y^2} + \frac{\partial^2 \Pi_e}{\partial z^2} \right) \right] & \mathcal{L}_{z_l} \left[ \frac{\partial^2 \Pi_e}{\partial y \partial x} \right] & \mathcal{L}_{z_l} \left[ \frac{\partial^2 \Pi_e}{\partial z \partial x} \right] \\
\mathcal{L}_{z_l} \left[ \frac{\partial^2 \Pi_e}{\partial x \partial y} \right] & -\mathcal{L}_{z_l} \left[ \left( \frac{\partial^2 \Pi_e}{\partial x^2} + \frac{\partial^2 \Pi_e}{\partial z^2} \right) \right] & \mathcal{L}_{z_l} \left[ \frac{\partial^2 \Pi_e}{\partial z \partial y} \right] \\
\mathcal{L}_{z_l} \left[ \frac{\partial^2 \Pi_e}{\partial x \partial z} \right] & \mathcal{L}_{z_l} \left[ \frac{\partial^2 \Pi_e}{\partial y \partial z} \right] & \mathcal{L}_{z_l} \left[ -\left( \frac{\partial^2 \Pi_e}{\partial x^2} + \frac{\partial^2 \Pi_e}{\partial y^2} \right) \right]
\end{bmatrix}_{z=\delta},
\quad (82)
$$

and

$$
\mathcal{L}_{z_l}[\mathcal{M}_m(x, y, \delta)] = i\omega \cdot \begin{bmatrix} 0 & \mathcal{L}_{z_l}\left[-\frac{\partial \Pi_m}{\partial z}\right] & \mathcal{L}_{z_l}\left[\frac{\partial \Pi_m}{\partial y}\right] \\ \mathcal{L}_{z_l}\left[\frac{\partial \Pi_m}{\partial z}\right] & 0 & \mathcal{L}_{z_l}\left[-\frac{\partial \Pi_m}{\partial x}\right] \\ \mathcal{L}_{z_l}\left[-\frac{\partial \Pi_m}{\partial y}\right] & \mathcal{L}_{z_l}\left[\frac{\partial \Pi_m}{\partial x}\right] & 0 \end{bmatrix}_{z=\delta} \quad (83)
$$

To prove the results in Eqs. (36) and (37), resolve Eqs. (82)
and (83) is possible in Fourier space. Let us write the
Fourier transform of one the terms in Eq. (83)

$$
\begin{align}
\mathcal{F}\left[\mathcal{L}_{z_l}\left[-i\omega \frac{\partial \Pi_m}{\partial z}\bigg|_{z=\delta}\right]\right](u, v) &= -i\omega \mathcal{F}\left[\frac{\partial \Pi_m}{\partial z}\bigg|_{z=\delta} * L_{z_l}(x, y)\right], \nonumber \\
&= -i\omega \frac{\partial}{\partial z} \mathcal{F}[\Pi_m]_{z=\delta} \cdot \mathcal{F}[L_{z_l}(x, y)], \nonumber \\
&= -i\omega \frac{\partial}{\partial z} (i2\pi) \frac{\exp\left(i\left|z\right|\sqrt{k^2 - 4\pi^2(u^2 + v^2)}\right)}{\sqrt{k^2 - 4\pi^2(u^2 + v^2)}} \bigg|_{z=\delta} \cdot \operatorname{sign}(z_l) \nonumber \\
&\quad \exp\left(i\left|z_l\right|\sqrt{k^2 - 4\pi^2(u^2 + v^2)}\right), \nonumber \\
&= -i\omega \frac{\partial}{\partial z} (i2\pi) \frac{\exp\left(i\left|z_l+z\right|\sqrt{k^2 - 4\pi^2(u^2 + v^2)}\right)}{\sqrt{k^2 - 4\pi^2(u^2 + v^2)}} \bigg|_{z=\delta'} \tag{84}
\end{align}
$$

with $z_l > 0$. Next, we substitute $Z = |z| + z_l$ and $\frac{\partial}{\partial z} = \text{sign}(z) \frac{\partial}{\partial Z}$ in Eq. (84). We get then

$$
\begin{align}
\mathcal{F}\left[\mathcal{L}_{z_l}\left[-i\omega \frac{\partial \Pi_m}{\partial z}\bigg|_{z=\delta}\right]\right](u, v) & \nonumber \\
&= -i\omega \operatorname{sign}(\delta) \frac{\partial}{\partial Z} (i2\pi) \frac{\exp\left(iZ\sqrt{k^2 - 4\pi^2(u^2 + v^2)}\right)}{\sqrt{k^2 - 4\pi^2(u^2 + v^2)}} \bigg|_{Z=|\delta|+z_l}, \nonumber \\
&= -i\omega \operatorname{sign}(\delta) \frac{\partial}{\partial z} \mathcal{F}[\Pi_m]_{z=|\delta|+z_l}, \nonumber \\
&= \mathcal{F}\left[-i\omega \operatorname{sign}(\delta) \frac{\partial \Pi_m}{\partial z}\bigg|_{z=|\delta|+z_l}\right] \tag{85}
\end{align}
$$

Finally, we have

$$
\mathcal{L}_{z_l} \left[ -i\omega \frac{\partial \Pi_m}{\partial z} \bigg|_{z=\delta} \right] = -i\omega \operatorname{sign}(\delta) \frac{\partial \Pi_m}{\partial z} \bigg|_{z=|\delta|+z_l} \quad (86)
$$

The same process could be realized to all terms of (82) and
(83) and give the results in Eqs. (36), (37) and (38).

**Appendix F: Integrals of the bessel product**

We outline an approach for computing the integral
defined in Eq. (51),

$$
\mathcal{I}_{nm}(\beta, \omega, s) = \int_{0}^{+\infty} \exp \left( i\beta \sqrt{1 - X^2} \right) J_{n+1} (2\pi \omega X) J_{|m|} (2\pi \omega s X) \, dX. \quad (87)
$$

Firstly, the recurrence relation of the Bessel functions [3],
Eqs. (9.1.27) on p. 361,

$$
J_{n+1}(2\pi\omega X) = \frac{\pi\omega X}{(n+1)} (J_n(2\pi\omega X) + J_{n+2}(2\pi\omega X)), \quad (88)
$$

is necessary. In this case, the integral $\mathcal{I}_{nm}$ becomes

$$
\mathcal{I}_{nm}(\beta, \omega, s) = \frac{\pi \omega}{(n+1)} \left[ \tilde{\mathcal{I}}_{n,m}(\beta, \omega, s) + \tilde{\mathcal{I}}_{n+2,m}(\beta, \omega, s) \right] + \frac{\pi \omega}{(n+1)} \int_{1}^{+\infty} \exp \left( -\beta \sqrt{X^2 - 1} \right) J_{n+1}(2\pi \omega X) J_{|m|}(2\pi \omega s X) \, dX. \quad (89)
$$

with

$$
\tilde{\mathcal{I}}_{n,m}(\beta, \omega, s) = \int_0^1 \exp \left( i\beta \sqrt{1 - X^2} \right) J_n \left( 2\pi\omega X \right) J_{|m|} \left( 2\pi\omega s X \right) X dX. \quad (90)
$$

The evaluation of $\tilde{\mathcal{I}}_{n,m}$ where $n$ and $m$ are integers such that $n - |m|$ is even and non-negative, need some relations. The first ones are the Zernike expansion of Bessel functions (see [7], Eq. (54) in Appendix B)

$$
J_n(2\pi\omega X) = \sum_{k=0}^{\infty} 2(n+2k+1)(-1)^k \frac{J_{n+2k+1}(2\pi\omega)}{2\pi\omega} R_{n+2k}^{n}(X), \quad 0 \le X \le 1, \quad (91)
$$

and

$$
J_{|m|}(2\pi\omega s X) = \sum_{t=0}^{\infty} 2(|m|+2t+1)(-1)^t \frac{\Gamma(|m|+2t+1)(2\pi\omega s)}{2\pi\omega s} R^{|m|}_{|m|+2t}(X), \quad 0 \le X \le 1. \tag{92}
$$
The coefficients of the expansion, in terms of Bessel functions, are bounded and exhibit super-exponential decay in $t$ and $k$ as soon as $|m| + 2t \ge 2\pi\omega s$, $n + 2k \ge 2\pi\omega$. Next, the important second relation for the evaluation, is (see [32], Eq. (82))

$$ \exp[i\beta\sqrt{1-X^2}] = \sum_{r=0}^{\infty} B_r(\beta) \cdot R_{2r}^0(X), \qquad (93) $$

where the function $B_r(\beta)$, is given by

$$ \mathcal{B}_r(\beta) = (2r+1) \frac{1}{2}\beta \left[ j_{r-1}\left(\frac{1}{2}\beta\right) h_r^{(1)}\left(\frac{1}{2}\beta\right) - j_r\left(\frac{1}{2}\beta\right) h_{r+1}^{(1)}\left(\frac{1}{2}\beta\right) \right], \quad \beta > 0, \qquad (94) $$

with $j_r$ and $h_r^{(1)}$ the spherical Bessel and Hankel functions, see [23], 10.4.7, of order $r$. The $B_r(\beta)$ are bounded, which is obvious, and actually decay as $\beta/(2r^2)$ when $r$ exceeds $\frac{1}{2}\beta$, which is not obvious. The convergence of the triple series is rapid, due to the rapid decrease of the two jinc-functionss (in $k$ and $t$, respectively) and the fact that the remaining integrals in Eq. (96) are bounded and are non-vanishing only when

$$ r \le \frac{1}{2}(n + |m|) + k + t. \qquad (95) $$

This latter part follows from the part that $R_{n+2k}^n(X)R_{|m|+2t}^{|m|}(X)$ is a polynomials with even powers $\le n+|m|+2k+2t$ of $X$ while $\int_0^1 X^2 R_{2r}^0(X)XdX = 0$ when $r > j$. Hence, when $r$ is large, then at least one of $k$ and $t$ is large as well and so the corresponding jinc-function is exponentially small. Thus we get

$$ \begin{aligned} & \int_0^1 \exp(i\beta\sqrt{1-X^2}) J_n(2\pi\omega X) J_{|m|}(2\pi\omega s X) XdX \\ & = \sum_{t,k,r} 4(n+2k+1)(|m|+2t+1)(-1)^{k+t} \times \\ & \qquad \mathcal{B}_r(\beta) \frac{J_{n+2k+1}(2\pi\omega) J_{|m|+2t+1}(2\pi\omega s)}{2\pi\omega 2\pi\omega s} \\ & \int_0^1 R_{2r}^0(X) R_{n+2k}^n(X) R_{|m|+2t}^{|m|}(X) XdX. \end{aligned} \qquad (96) $$

We shall give below a method for systematically writing

$$ R^{|m|}_{|m|+2t}(X) R^n_{n+2k}(X) = \sum_{q=0}^{\frac{1}{2}(n+|m|)+k+t} \mathcal{L} C^{|m|,n,0}_{tkq} \cdot R^0_{2q}(X), \qquad (97) $$

where $\mathcal{L}C$ stands for linearization coefficient. As a consequence and by the orthogonality of the radial polynomials, the remaining integrals in Eq. (96) are equal to

$$ \int_0^1 R_{2r}^0(X) R_{n+2k}^n(X) R_{|m|+2t}^{|m|} (X) XdX = \sum_{q=0}^{\frac{n+|m|}{2}+k+t} \frac{\mathcal{L} C^{|m|,n,0}_{tkq}}{2(2q+1)} \delta_{rq}, \qquad (98) $$

where $\delta_{rq}$ is the Kronecker delta, equal to 1 if $r = q$, and 0 otherwise. The $\mathcal{L}C^{|m|,n,0}_{tkq}$ are evaluated recursively in **Recursion step** below, where we have written

$n = |m| + 2L$, $L = 0, 1, \dots$. The final expression is thus

$$ \tilde{\mathcal{I}}_{n,m}(\beta, \omega, s) = \sum_{t,k} \mathcal{G}_{tk}^{n,|m|}(\beta, \omega) \frac{J_{|m|+2t+1}(2\pi\omega s)}{2\pi\omega s} \qquad (99) $$

with $k = 0, 1, \dots, t = 0, 1, \dots$, and

$$ \mathcal{G}_{tk}^{N,M}(\beta, \omega) = 4(-1)^{t+k} (N + 2k + 1)(M + 2t + 1) \frac{J_{N+2k+1}(2\pi\omega)}{2\pi\omega} \sum_{q=0}^{\frac{N+M}{2}+k+t} \frac{\mathcal{L}C_{tkq}^{M,N,0}}{2(2q+1)} \mathcal{B}_q(\beta). \qquad (100) $$

Note that, the number of the summation over $t, k, r, q$ has been reduced to $t, k, q$ due to the presence of the Kronecker delta. The summation over $q$ is a finite sum. This evaluation in Eq. (99) stays valid for the second term of Eq. (89) for $n+2$.

### Initialization of $\mathcal{L}C^{|m|,n,0}_{tkq}$ for $L=0$

For $L=0$, we have then the equality $n = |m|$. Since $\mathcal{L}C^{|m|,|m|,0}_{tkq} = \mathcal{L}C^{|m|,|m|,0}_{tkq}$, we assume that $k \le t$. Now, (see [32] Eq. (140) in Appendix C)

$$ R^{|m|}_{|m|+2k}(X) R^{|m|}_{|m|+2t} = \sum_{q=0}^{|m|+t+k} A^{|m|,-|m|,0}_{|m|+2t,|m|+2k,2q} R^0_{2q}(X), \qquad (101) $$

where $A$'s coefficients are linked to the $C$'s Clebsch-Gordan coefficients according to

$$ \begin{aligned} A^{|m|,-|m|,0}_{|m|+2t,|m|+2k,2q} &= \left| C^{\frac{1}{2}|m|+t,\frac{1}{2}|m|+k,q}_{\frac{1}{2}|m|,-\frac{1}{2}|m|,0} \right|^2 \\ &= (2q+1) \left| \begin{pmatrix} \frac{1}{2}|m|+t & \frac{1}{2}|m|+k & q \\ \frac{1}{2}|m| & -\frac{1}{2}|m| & 0 \end{pmatrix} \right|^2. \end{aligned} \qquad (102) $$

The last term at the right hand side of Eq. (102) corresponds to the 3j-symbol in the notation of [25], Chap. 34. By the symmetry properties,

$$ \begin{pmatrix} j_1 & j_2 & j_3 \\ m_1 & m_2 & m_3 \end{pmatrix} = \begin{pmatrix} j_3 & j_1 & j_2 \\ m_3 & m_1 & m_2 \end{pmatrix}, \qquad (103) $$

we have from Eqs. (97) and (101)

$$ \begin{aligned} \mathcal{L}C^{|m|,|m|,0}_{tkq} &= A^{|m|,-|m|,0}_{|m|+2t,|m|+2k,2q} \\ &= (2q+1) \left| \begin{pmatrix} q & \frac{1}{2}|m|+t & \frac{1}{2}|m|+k \\ 0 & \frac{1}{2}|m| & -\frac{1}{2}|m| \end{pmatrix} \right|^2. \end{aligned} \qquad (104) $$

The computation of

$$ \begin{pmatrix} j_1 = q & j_2 = \frac{1}{2}n' = \frac{1}{2}|m|+t & j_3 = \frac{1}{2}h = \frac{1}{2}|m|+k \\ m_1 = 0 & m_2 = \frac{1}{2}|m| & m_3 = -\frac{1}{2}|m| \end{pmatrix}, \quad q = 0, 1, \dots, \qquad (105) $$

can be done recursively according to [32], Eqs. (63-68) and (150-161), where we may restrict to the case that
$m' = |m| \le h = |m| + 2k \le n' = |m| + 2t$ since we have assumed that $k \le t$. In this case, we have

$$ \begin{pmatrix} \frac{1}{2}|n' - h| - 1 & \frac{1}{2}n' & \frac{1}{2}h \\ 0 & \frac{1}{2}|m| & -\frac{1}{2}|m| \end{pmatrix} = 0, \qquad (106) $$

$$ \begin{pmatrix} \frac{1}{2}|n' - h| & \frac{1}{2}n' & \frac{1}{2}h \\ 0 & \frac{1}{2}|m| & -\frac{1}{2}|m| \end{pmatrix} = (-1)^{\frac{n'-|m|}{2}} \left( \frac{(n' - h)! h!}{(n' + 1)!} \right)^{1/2} \cdot \begin{pmatrix} \left( \frac{n' + |m|}{2} \right)! \left( \frac{n' - |m|}{2} \right)! \\ \left( \frac{h + |m|}{2} \right)! \left( \frac{h - |m|}{2} \right)! \end{pmatrix}^{1/2} \frac{1}{\left( \frac{n' - h}{2} \right)!} \qquad (107) $$

and

$$ \begin{aligned} |m|(2q + 1) \begin{pmatrix} q & \frac{1}{2}n' & \frac{1}{2}h \\ 0 & \frac{1}{2}|m| & -\frac{1}{2}|m| \end{pmatrix} = & \left( q^2 - \left( \frac{n' - h}{2} \right)^2 \right)^{1/2} \\ & \left( \left( \frac{n' + h}{2} + 1 \right)^2 - q^2 \right)^{1/2} \begin{pmatrix} q - 1 & \frac{1}{2}n' & \frac{1}{2}h \\ 0 & \frac{1}{2}|m| & -\frac{1}{2}|m| \end{pmatrix} + \\ & \left( (q+1)^2 - \left( \frac{n' - h}{2} \right)^2 \right)^{1/2} \left( \left( \frac{n' + h}{2} + 1 \right)^2 - (q+1)^2 \right)^{1/2} \\ & \begin{pmatrix} q + 1 & \frac{1}{2}n' & \frac{1}{2}h \\ 0 & \frac{1}{2}|m| & -\frac{1}{2}|m| \end{pmatrix} \end{aligned} \qquad (108) $$

for $q = \frac{1}{2}|n' - h|, \frac{1}{2}|n' - h| + 1, \dots, \frac{1}{2}(n' + h) - 1$. The quantity

$$ (-1)^{\frac{n'-|m|}{2}} \left( \frac{(n' - h)! h!}{(n' + 1)!} \right)^{1/2} \cdot \left( \frac{\left( \frac{n'+|m|}{2} \right)! \left( \frac{n' - |m|}{2} \right)!}{\left( \frac{h+|m|}{2} \right)! \left( \frac{h-|m|}{2} \right)!} \right)^{1/2} \frac{1}{\left( \frac{n'-h}{2} \right)!} \qquad (109) $$

is easily computed recursively by noting its simple form for $h = n'$ and the simple form of the ratio of two consecutive quantities (109) as a function of $h = n', n' - 2, \dots, m'$. This recursive computation of (109) is as follows. Let $n'$ and $m'$ be non-negative integers with $n' - m'$ even and non-negative. Consider

$$ T_h = \frac{(n' - h)! h!}{(n' + 1)!} \left( \frac{\left( \frac{n' + m'}{2} \right)! \left( \frac{n' - m'}{2} \right)!}{\left( \frac{h + m'}{2} \right)! \left( \frac{h - m'}{2} \right)!} \right) \left( \frac{1}{\left( \frac{n' - h}{2} \right)!} \right)^2, \qquad (110) $$

for $h = m', m' + 2, \dots, n'$. Then

$$ T_{n'} = \frac{1}{n' + 1}, \qquad (111) $$

and, for $h = m' + 2, m' + 4, \dots, n'$,

$$ \begin{aligned} T_{h-2} &= \frac{(n' - h + 2)!(h - 2)!}{(n' + 1)!} \frac{\left( \frac{n' + m'}{2} \right)! \left( \frac{n' - m'}{2} \right)!}{\left( \frac{h + m'}{2} - 1 \right)! \left( \frac{h - m'}{2} - 1 \right)!} \\ & \qquad \left( \frac{1}{\left( \frac{n' - h}{2} + 1 \right)!} \right)^2 \qquad (112) \\ &= T_h \cdot \frac{n' - h + 1}{n' - h + 2} \cdot \frac{(h + m')(h - m')}{h(h - 1)}, \qquad (113) \end{aligned} $$

The square root of the result of the recurrence allow to obtain the result of Eq. (109).

Expressing $R^{|m|+2L+2}|m|+2L+2B}$ as a linear combination of $R^{|m|+2L}|m|+2L+2A$, $A = 0, 1, \dots, B + 1$.

Let $L = 0, 1, \dots, B = 0, 1, \dots$. There is an explicit series representation

$$ R^{|m|+2L+2}|m|+2L+2B(X) = \sum_{A=0}^{B+1} CC_{BA}^{|m|+2L+2, |m|+2L} R^{|m|+2L}|m|+2L+2A(X), \qquad (114) $$

where CC stands for connection coefficient. This follows from [32], Eq. (200) with $m' = |m| + 2L + 2$, $h = |m| + 2L + 2B + 2$; so that $p = \frac{1}{2}(h - m') = B$, $q = \frac{1}{2}(h + m') = |m| + 2L + B + 2$ that thus takes the form

$$ \begin{aligned} R^{|m|+2L+2}|m|+2L+2B(X) = & \frac{B+1}{|m|+2L+B+2} R^{|m|+2L}|m|+2L+2B+2(X) \\ & + \frac{|m|+2L+1}{(|m|+2L+B+2)(B+1)} \sum_{D=1}^{B+1} \\ & \qquad \left( |m|+2L+2(B+1-D)+1 \right) \frac{\binom{B+1}{D}}{
\binom{|m|+2L+B+1}{D} } \\ & (-1)^{D+1} R^{|m|+2L}|m|+2L+2(B+1-D)(X). \end{aligned} \qquad (115) $$

The ratio of the two binomials occurring in the series over $D$ can be conveniently computed recursively in $D = 0, 1, \dots, B + 1$:

$$ \frac{\binom{B+1}{D}}{\binom{E+1}{D}} = 1, \quad D=0; \quad \frac{\binom{B+1}{D+1}}{\binom{E+1}{D+1}} = \frac{\binom{B+1}{D}}{\binom{E+1}{D}} \cdot \frac{B+1-D}{E+1-D}, \quad D=0, 1, \dots, B. \qquad (116) $$

Thus, letting $Q = |m| + 2L + B + 2$, we have

$$ CC_{BA=B+1}^{|m|+2L+2, |m|+2L} = \frac{B+1}{Q}, \qquad (117) $$

and, writing $A = B + 1 - D = 0, 1, \dots, B$ in the series (115)

$$ CC_{BA}^{|m|+2L+2, |m|+2L} = \frac{(|m|+2L+1)(|m|+2L+2A+1)}{Q(B+1)} (-1)^{B-A} \frac{\binom{B+1}{A}}{\binom{Q-1}{B+1-A}}, \quad A=0, 1, \dots, B. \qquad (118) $$

### Recursion step

Let $L = 0, 1, \dots$ and assume that we have for a given $t = 0, 1, \dots$ available all quantities

$$ \mathcal{L}C_{tq}^{|m|, |m|+2L, 0}, \quad t = 0, 1, \dots, \quad q = 0, 1, \dots, |m| + L + t + k. \qquad (119) $$

Then

$$ \begin{aligned} R^{|m|}_{|m|+2t}(X) R^{|m|+2L+2}_{|m|+2L+2+2k}(X) &= \sum_{A=0}^{k+1} CC_{kA}^{|m|+2L+2, |m|+2t} R^{|m|}_{|m|+2t}(X) R^{|m|+2L}_{|m|+2L+2A}(X), \\ &= \sum_{A=0}^{k+1} CC_{kA}^{|m|+2L+2, |m|+2L} \sum_{q=0}^{|m|+L+t+A} \mathcal{L}C_{tq}^{|m|, |m|+2L, 0} R^0_{2q}(X). \end{aligned} \qquad (120) $$
We rearrange the repeated series, writing $E = |m| + L + t$ according to

$$ \sum_{A=0}^{k+1} CC_{kA}^{|m|+2L+2,|m|+2L} \cdot \sum_{q=0}^{E+A} CC_{tAq}^{|m|,|m|+2L,0} R_{2q}^0(X) \quad (121) $$

$$ = \sum_{A=0}^{k+1} \sum_{q=0}^{E+A} CC_{kA}^{|m|+2L+2,|m|+2L} \cdot CC_{tAq}^{|m|,|m|+2L,0} R_{2q}^0(X) \quad (122) $$

$$ = \sum_{q=0}^{E+k+1} \left[ \sum_{A=\max(0,q-E)}^{k+1} CC_{kA}^{|m|+2L+2,|m|+2L} CC_{tAq}^{|m|,|m|+2L,0} \right] R_{2q}^0(X). \quad (123) $$

Hence

$$ CC_{tkq}^{|m|,|m|+2L+2,0} = \sum_{A=\max(0,q-E)}^{k+1} CC_{kA}^{|m|+2L+2,|m|+2L} \cdot CC_{tAq}^{|m|,|m|+2L,0} \quad (124) $$

for $k = 0, 1, \dots$ and $q = 0, 1, \dots, |m|+t+k+L+1$ where $E = |m|+L+t$. Note that for $L = 0$, $CC_{tAq}^{|m|,|m|,0}$ is defined by Eq. (104) to initialize the recursion.

### Evaluation of the integral $\mathcal{J}_{nmk}(\beta, \omega, s)$ in Eq. (57)

To evaluate this integral, it is not necessary to use the recurrence relation of the Bessel functions. This integral is identical to the integral in Eq. (90) except for the exponential function. In this case, there is the $R_{2r}^0$-expansion, see [31], Eq. (35),

$$ \frac{\exp\left[\frac{i\beta\sqrt{1-X^2}}{\sqrt{1-X^2}}\right]}{\sqrt{1-X^2}} = \sum_{r=0}^{\infty} A_r(\beta) \cdot R_{2r}^0(X), \quad (125) $$

with

$$ A_r(\beta) = i\beta(2r+1)j_r\left(\frac{1}{2}\beta\right)h_r^{(1)}\left(\frac{1}{2}\beta\right), \quad \beta > 0. \quad (126) $$

We replace $B_q(\beta)$ by $A_q(\beta)$ in the final result Eq. (100). Thus,

$$ \mathcal{J}_{nmk}(\beta, \omega, s) = \sum_{t,l} \mathcal{H}_{tl}^{n+1,1|1+m-2k|}(\beta, \omega) \frac{J_{1+m-2k|+2t+1}(2\pi\omega s)}{2\pi\omega s} \quad (127) $$

with $l = 0, 1, \dots, t = 0, 1, \dots, k = 0, 1$ and

$$ \mathcal{H}_{tl}^{N,M}(\beta, \omega) = 4(-1)^{t+l}(N+2l+1)(M+2t+1) \frac{J_{N+2l+1}(2\pi\omega)}{2\pi\omega} \sum_{q=0}^{\frac{N+M}{2}+l+t} \frac{\mathcal{L}_{tlq}^{M,N,0}}{2(2q+1)} A_q(\beta) \quad (128) $$

with $A_q(\beta)$ defined in Eq. (126).

### Cases of evanescente waves

Now, we must treat the integral correspond to the evanescent wave. From Eq (89), we note

$$ \int_{1}^{+\infty} \exp(-\beta\sqrt{X^2-1}) J_{n+1}(2\pi\omega X) J_{|m|}(2\pi\omega sX) dX = I, \quad (129) $$

which becomes with $Y = (X^2 - 1)^{1/2}$, $dX = Y(1 + Y^2)^{-1/2}dY$,

$$ I = \int_{0}^{+\infty} \exp(-\beta Y) J_{n+1}(2\pi\omega\sqrt{1+Y^2}) J_{|m|}(2\pi\omega s\sqrt{1+Y^2}) (1+Y^2)^{-1/2} Y dY. \quad (130) $$

As $\beta$ is large, greater than 30, $\exp(-\beta Y)$ becomes zero everywhere except close to the point $Y = 0$. Thus, with an slowly variable envelop approximation, and by applying the steepest descent method, we have

$$ \begin{aligned} I &\approx J_{n+1}(2\pi\omega)J_{|m|}(2\pi\omega s) \int_{0}^{+\infty} \exp(-\beta Y)Y dY \\ &= \frac{1}{\beta^2} J_{n+1}(2\pi\omega)J_{|m|}(2\pi\omega s) \end{aligned} \quad (131) $$

This approximation can be used to evaluate the integral over $X \in [1, +\infty)$ in Eq. (57). Alternatively, the integrals in Eq. (129) and Eq. (130) can be done by numerical integration in a much easier way than, for instance, the integral $\tilde{I}_{n,m}$ in Eq. (90).

### Acknowledgements

Not applicable.

### Authors' contributions

All authors contributed equally in compiling this paper. The authors read and approved the final manuscript.

### Funding

This project has been funded by the CORIA.

### Availability of data and materials

Upon request to the authors.

### Declarations

### Competing interests

The authors declare that they have no competing interests.

### Author details

¹Département d'optique, UMR-6614, Av. de l'Université, 76801 Saint-Etienne du Rouvray, France. ²Department of Mathematics and Computer Science, Eindhoven University of Technology, MB Eindhoven, The Netherlands.

Received: 14 November 2020 Accepted: 24 May 2021

Published online: 30 June 2021

### References

1. Aarts, R. M., Janssen, A. J. E. M.: Sound radiation quantities arising from a resilient circular radiator. J. Acoust. Soc. Am. **126**(4), 1776-1787 (2009)

2. Aarts, R. M., Janssen, A. J. E. M.: On-axis and far-field sound radiation from resilient flat and dome-shaped radiators. J. Acoust. Soc. Am. **125**(3), 1444-1455 (2009)

3. Abramowitz, M., Stegun, I.: Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (1972). https://doi.org/10.2307/1266136

4. Andreas, B., Mana, G., Palmisano, C.: Vectorial ray-based diffraction integral. J. Opt. Soc. Am. A **32**(8), 1403-1424 (2015)

5. Born, M., Wolf, E., Bhatia, A. B., Clemmow, P. C., Gabor, D., Stokes, A. R., Taylor, A. M., Wayman, P. A., Wilcock, W. L.: Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light. 7 edition. Cambridge University Press (1999). https://doi.org/10.1063/1.1325200

6. Braat, J. J. M., Dirksen, P., Janssen, A. J. E. M., van de Nes, A. S.: Extended nijboer-zernike representation of the vector field in the focal region of an aberrated high-aperture optical system. J. Opt. Soc. Am. A **20**(12), 2281-2292 (2003)
7. Braat, J. J. M., Haver, S., Janssen, A. J. E. M., Pereira, S. F.: Image formation in a multilayer using the extended nijboer-zernike theory. J. Eur. Opt. Soc. Rapid Publ. **4**(09048) (2009). https://doi.org/10.2971/jeos.2009.09048

8. Braat, J. J. M., Dirksen, P., Janssen, A. J. E. M., van Haver, S., van de Nes, A. S.: Extended nijboer-zernike approach to aberration and birefringence retrieval in a high-numerical-aperture optical system. J. Opt. Soc. Am. A. **22**(12), 2635-2650 (2005)

9. Braat, J. J. M., van Haver, S., Janssen, A. J. E. M., Dirksen, P.: Assessment of optical systems by means of point-spread functions. Prog. Opt. **51**, 349-468 (2008)

10. Carretero, L., Acebal, P., Blaya, S.: Diffraction of convergent spherical waves with all possible polarization states using the luneburg integral method. J. Opt. Soc. Am. A. **30**(4), 733-740 (2013)

11. Coëtmellec, S., Remacha, C., Brunel, M., Lebrun, D., Janssen, A: Digital in-line holography with a spatially partially coherent beam. J. Eur. Opt. Soc. Rapid Publ. **6**(11060) (2011). https://doi.org/10.2971/jeos.2011.11060

12. Coëtmellec, S., Wichitwong, W., Gréhan, G., Lebrun, D., Brunel, M., Janssen, A. J. E. M: Digital in-line holography assessment for general phase and opaque particle. J. Eur. Opt. Soc. Rapid Publ. **9**, 14021, 06 (2014)

13. Davis, L. W.: Theory of electromagnetic beams. Phys. Rev. A. **19**, 1177-1179 (1979)

14. Gouesbet, G.: and Gréhan, G., Generalized Lorenz-Mie Theories. Physics, Classical Continuum Physics. Springer Berlin, Heidelberg (2011)

15. Grare, S., Coëtmellec, S., Allano, D., Grehan, G., Brunel, M., Lebrun, D: Dual wavelength digital holography for 3d particle image velocimetry. J. Eur. Opt. Soc. Rapid Publ. **10**(15009) (2015). https://doi.org/10.1364/dh.2015.1dm2a3

16. Janssen, A. J. E. M: New analytic results for the zernike circle polynomials from a basic result in the nijboer-zernike diffraction theory. J. Eur. Opt. Soc. Rapid Publ. **6**(11028) (2011). https://doi.org/10.2971/jeos.2011.11028

17. Kim, J., Wang, Y., Zhang, X.: Calculation of vectorial diffraction in optical systems. J. Opt. Soc. Am. A. **35**(4), 526-535 (2018)

18. King, L. V: On the acoustic radiation field of the piezo-electric oscillator and the effect of the viscosity on transmission. Can. J. Res. **11**(2), 135-155 (1934)

19. Lewis, W. E., Vyas, R.: Maxwell-gaussian beams with cylindrical polarization. J. Opt. Soc. Am. A. **31**(7), 1595-1603 (2014)

20. Luneberg, R. K.: Mathematical theory of optics. University of California Press, Berkeley and Los Angeles (1964)

21. Marathay, A. S., McCalmont, J. F: On the usual approximation used in the rayleigh-sommerfeld diffraction theory. J. Opt. Soc. Am. A. **21**(4), 510-516 (2004)

22. Masson, J. B., Gallot, G.: Diffraction from a subwavelength elliptic aperture: analytic approximate aperture fields. J. Opt. Soc. Am. A. **29**(9), 2005-2014 (2012)

23. Mori, M.: Discovery of the double exponential transformation and its developments. Publ. RIMS, Kyoto Univ. **41**, 897-935 (2005)

24. Noll, R. J.: Zernike polynomials and atmospheric turbulence. J. Opt. Soc. Am. **66**(3), 207-211 (1976)

25. Olver, F., Lozier, D., Boisvert, R., Clark, C.: Nist handbook of mathematical functions. Cambridge University Press, New York (2010)

26. Pejchang, D., Coëtmellec, S., Gréhan, G., Brunel, M., Lebrun, D., Chaari, A., Grosges, T., Barchiesi, D.: Recovering the size of nanoparticles by digital in-line holography. Opt. Express. **23**(14), 18351-18360 (2015)

27. Picart, P., Montrésor, S., Fages, M., Xia, H., Guo, R., Li, J., Solieman, O., Durand, J. C.: Analysis of computerized aided designed and manufactured dental occlusal ceramics with multi-wavelength digital holography. In: SPECKLE 2018: VII International Conference on Speckle Metrology, p. 50. SPIE, SPIE, Janów Podlaski, Poland, (2018)

28. Rdzanek, W. P.: Sound scattering and transmission through a circular cylindrical aperture revisited using the radial polynomials. J. Acoust. Soc. Am. **143**(3), 1259-1282 (2018)

29. Ren, H., Shao, W., Li, Y., Salim, F., Gu, M: Three-dimensional vectorial holography based on machine learning inverse design. Sci. Adv. **6**(16) (2020). https://doi.org/10.1126/sciadv.aaz4261

30. Romero, J. A., Hernández, L.: Diffraction by a circular aperture: an application of the vectorial theory of huygens's principle in the near field. J. Opt. Soc. Am. A. **25**(8), 2040-2043 (2008)

31. van Haver, S., Braat, J. M., Janssen, A. J. E. M., Janssen, O. T. A., Pereira, S. F.: Vectorial aerial-image computations of three-dimensional objects based

on the extended nijboer-zernike theory. J. Opt. Soc. Am. A. **26**(5), 1221-1234 (2009)

32. van Haver, S., Janssen, A. J. E. M: Advanced analytic treatment and efficient computation of the diffraction integrals in the extended nijboer-zernike theory. J. Eur. Opt. Soc. Rapid Publ. **8**(13044) (2013). https://doi.org/10.2971/jeos.2013.13044

33. van Haver, S., Janssen, A. J. E. M: Truncation of the series expressions in the advanced enz-theory of diffraction integrals. J. Eur. Opt. Soc. Rapid Publ. **9**(14042) (2014). https://doi.org/10.2971/jeos.2014.14042

34. Weyl, H.: Ausbreitung elektromagnetischer wellen über einem ebenen leiter. Ann. Phys. **60**, 481 (1919)

35. Yanagawa, T., Abe, R., Hayasaki, Y.: Three-dimensional mapping of fluorescent nanoparticles using incoherent digital holography. Opt. Lett. **40**(14), 3312-3315 (2015)

36. Zernike, F.: Diffraction theory of the knife-edge test and its improved version, the phase-contrast method. Physica (Amsterdam). **1**, 689-704 (1934). https://doi.org/10.1093/mnras/94.5.377

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Submit your manuscript to a SpringerOpen journal and benefit from:

* Convenient online submission

* Rigorous peer review

* Open access: articles freely available online

* High visibility within the field

* Retaining the copyright to your article

Submit your next manuscript at ► springeropen.com