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

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.


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 *Correspondence: coetmel@coria.fr 1 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 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 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/λ 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 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.

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 z l 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 , 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 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 wavelenght 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 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 where r 2 = (x − ξ) 2 + (y − η) 2 + z 2 l with z l = z − δ > 0 due to the δ-shifting along the optical axis. If the third component is undetermined, it could be derived from (2) and (3) by means of the Maxwells'equation div E z l = 0 (see Appendix A). We get where k = 2π/λ is the wave number and is the vectorial field in the sensor plane. The wave E z l satisfies the Helmholtz equation in the half space z l > 0. The uniqueness of the solution E i z l 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 and in any solid sector π/2 + < θ < π/2 − of the domain z l > 0, there exists a constant B( ) such that for all points (x, y, z l ) of the sector we have the second condition It is important to note that the conditions in Eq. (5) apply on the function W i δ (ξ , η) too. Indeed, for z l = 0, we have the boundary limit (see Appendix B) 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) with and for Eq. (3)

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 L z l then 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 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 aŝ wheref (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 tô The evaluation of the Fourier transform of the Luneberg kernel is realized by means of [9], Eq. (A8), p. 112 and gives us In a similar way, the propagation ofÊ By introducing Eq. (14) into Eq. (16), we obtain the amplitude spectrum at the distance z l such aŝ Knowing that z l1 > 0 and z l2 > 0, in the spatial domain, we have and this is Eq. (12). Note that 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 where E i z l and W i δ are scalar components and 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 under the following condition, see Appendix C 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. 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 − circ(r/R), with circ(r/R) = 1 if r = (x 2 + y 2 ) 1/2 ≤ 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 on 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 z l = 18.98mm. 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 z l 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.

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 where the velocity of light in the vacuum is approximately equal to c = 3 · 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 The magnetic vector and electric scalar potentials are expressed versus the Hertz electric potential, denoted e (r, t), and the magnetic potential, denoted m (r, t) by writing that e = −div e (r, t). The expressions of the vectorial electromagnetic fields (E δ , B δ ) versus ( e , m ) in exp(−iωt)-Fourier space are the following where μ 0 , ε 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 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: By introducing the Eqs. (30) and (29) into (28) and (27), the expression of the matrices M e and M m , are and The lower index z = δ indicates that the derivations in the matrices have to be taken at the value z = δ with positive or negative value of δ. The authors of [10] have been developed these matrices in circular coordinates (ρ, θ) and taken a negative value of δ (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 M e and M m in Appendix D. Nevertheless, we will choose the particular case illustrated in Fig. 1 for our illustrations, i.e. δ > 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 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. E δ (ξ , η). From the point of view of the definition of the object, we have previously mentioned that the amplitude function P(ξ , η) = 1 − circ(ρ/R), R = D/2. Thus, the amplitude of the wave in the transverse plane just next the object is given by with i = x, y and ρ 2 = ξ 2 + η 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 The first term of Eq. (34) is the reference vectorial wave. It will be denoted E r 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 L z l to the matrices M e and 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) with and The function 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 |δ| + z l . It is clear that when δ > 0 then and E rz z l of the electric field where and 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 δ = 5λ and it propagates over z l = 3λ. 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 i δ (ξ , η) circ(ρ/R) as in which and where n and m in the summation in Eq. (39) and in Eq. (40) are integers such that n − |m| is even and nonnegative. The Zernike coefficients γ m ni are obtained by using the orthogonality of the Zernike circle polynomials Z m n , so that  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λ localized at δ = 5λ. The polarization vectors p e = (1, 0, 0) and 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 , and imaginary part, denoted , 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Ẽ x δ (ξ , η) 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Ẽ x δ (ξ , η) component has a very slight rotational asymmetry. On the other hand, the decomposition of theẼ with σ = R √ u 2 + v 2 and ϕ = arg(u + iv) allows us to obtain the spectrum, denotedÊ i z l (σ , ϕ), of the second term of Eq. (34) as For σ ≥ R/λ, 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 σ = 0. The inverse Fourier transform with a circular symmetry (Hankel transform) is By noting that 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  , ϕ), of E z z l (x, y) with the same processes of calculus as previously, we obtain By noting that 2π 0 γ m nx cos ϕ +γ m ny sin ϕ e imϕ exp (i2πσ s cos (θ −ϕ)) dϕ the expression of the longitudinal component in the spatial domain of the vectorial field is defined by 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, +∞) separately.

Evaluation of the integrals of the Bessel product
The problem of the study of the Bessel function J ν (x) in which the order ν and argument x are positive, and then the product of the Bessel functions under the integral, are usually linked to the ratio x/ν 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 ω and, on the other hand, on the radial argument s and of ω for the other. It is clear that when s and ω 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 I nm (β, ω, s), the integral in Eq. (47) so that The evaluation of this integral is given in (F) and yields with β 0 and 0 ≤ k < ∞, 0 ≤ t < ∞ and An another point of view from [24,36], the outcome of the integral in Eq. (51) for β = 0 equals to for 0 ≤ s < 1 and 0 for s > 1. For the second integral in Eq. (50), denoted as J nmk (β, ω, s), so that 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.
with l = 0, 1, · · · , t = 0, 1, · · · , k = 0, 1 and with A q (β) 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 ω = 3.5 and β = 40.84 that correspond to the optical parameter λ = 633nm, D = 7λ, z l = 6.5λ. 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 chose 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 J nmk (β, ω, s), I nm (β, ω, 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.

Simulation of the holograms
In this section, we will show the intensity distribution of the hologram in the quadratic sensor plane, denoted I(s, θ). Its expression is as follows The intensity distribution has been simulated for various polarization states of the incident vectorial fields in the object. The authors of [10] 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 θ from the horizontal in the plane before the object, at the output of the polariser, we will have where 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 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 objet 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 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. approximation, for small angles. To introduce the numerical parameter such as the number of sampling N which is the the number of pixel of the CCD camera, it is necessary to introduce the instantaneous frequency, denoted f i such as where the phase φ 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 = x, with x is the spatial support of the image. Then f imax = x λz . The sampling frequency, denoted f e , must exceed twice the maximal frequency. Consequently, by written that f e ≥ 2f imax and by considering that the support x = Nδ e with δ e = 1/f e , we obtain If we suppose that the support y along y-axis is equal to x, then r max = 1 2 ( x 2 + y 2 ) 1/2 and with the condition in Eq. (75) combined with the Eq. (77), we get Finally, It is unlikely that the first term is greater than the second. Consequently, we can consider only the second term.

Appendix E: Free luneberg propagation of the electromagnetic field
In this appendix, we will only illustrate the development of one term of the matrices M e and M m . Before, let us note that the Luneberg transforms of the propagation matrices M e and M m is the matrix of the Luneberg transforms of each term of M e and M m and To prove the results in Eqs. (36) and (37) with z l > 0. Next, we substitute Z = |z| + z l and ∂ ∂z = sign(z) ∂ ∂Z in Eq. (84). We get then Finally, we have 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), Firstly, the recurrence relation of the Bessel functions [3], Eqs. (9.1.27) on p. 361, is necessary. In this case, the integral I nm becomes The evaluation ofĨ 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) and 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 ≥ 2πωs, n + 2k ≥ 2πω. Next, the important second relation for the evaluation, is (see [32], Eq. (82)) where the function B r (β), is given by with j r and h (1) r the spherical Bessel and Hankel functions, see [23], 10.4.7, of order r. The B r (β) are bounded, which is obvious, and actually decay as β/(2r 2 ) when r exceeds 1 2 β, 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 This latter part follows from the part that R n is a polynomials with even powers ≤ n +|m|+2k + 2t of X while 1 0 X 2j R 0 2r (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 We shall give below a method for systematically writing where LC stands for linearization coefficient. As a consequence and by the orthogonality of the radial polynomials, the remaining integrals in Eq. (96) are equal to where δ rq is the the Kronecker delta, equal to 1 if r = q, and 0 otherwise. The LC |m|,n,0 tkq are evaluated recursively in Recursion step below, where we have writ-ten n = |m| + 2L, L = 0, 1, · · · . The final expression is thusĨ with k = 0, 1, · · · , t = 0, 1, · · · , and 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 LC
where A's coefficients are linked to the C's Clebsch-Gordan coefficients according to 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, we have from Eqs. (97)  The computation of j 1 = q j 2 = 1 2 n = 1 2 |m| + t j 3 = 1 2 h = 1 2 |m| + k m 1 = 0 m 2 = 1 2 |m| m 3 = − 1 2 |m| , q = 0, 1, · · · , can be done recursively according to [32] , The square root of the result of the recurrence allow to obtain the result of Eq. (109).

Cases of evanescante waves
Now, we must treat the integral correspond to the evanescent wave. From Eq (89), we note As β is large, greater than 30, exp(−β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 This approximation can be used to evaluate the integral over X ∈ [ 1, +∞) 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 I n,m in Eq. (90).