1. INTRODUCTION

Soon after the invention of the laser, its applicability to unexcelled measurements of lengths and displacements by interferometry was recognized [1]. Plane or spherical waves were assumed to relate the phase of the interference fringes to the measurand, but, as early as 1976, Dorenwendt and Bönsch pointed out that this is not correct and that diffraction brings about systematic errors [2]. Since then, metrologists have routinely applied corrections based on scalar and paraxial approximations of the interferometer operations [3 –13]. In recent years, cutting-edge interferometric measurements carried out to determine the Avogadro constant are looking at $1\times {10}^{-9}$ relative accuracies over propagation distances of the order of 1 m and propagation differences of many centimeters. Consequently, metrologists are drawing attention to nonparaxial vector models, with a goal to verifying the validity of the scalar and paraxial approximations or to improving the correction calculation. Since these measurements use continuous-wave lasers, which are stabilized and have an extremely narrow bandwidth, monochromatic light propagation can be assumed.

Optical interferometers may combine many optical components in a potentially complex geometry. Wolf and Richards carried out early work in 1959 on diffraction in optical systems, but that work provided only a description of the image or focal region [14,15]. Among the tools based on the paraxial approximation, a common one is the Collins integral, which allows the electromagnetic field at the output plane to be calculated from the field in the input plane and the ray matrix of the optical system [16]. Ray matrices have been also applied to propagate Gaussian beams, as well as higher-order Laguerre–Gaussian or Hermite–Gaussian beams, which are solutions of the paraxial scalar wave equation [17,18]. Finally, Fourier optics is used together with the thin lens approximation [19].

In regards to nonparaxial models, in 1980, Byckling and Simola [20] proposed a sequential application of spherical harmonic Green’s functions and the evaluation of the relevant diffraction integral by the stationary phase method, but their method is limited to the scalar approximation and to spherical interfaces.

Geometric optics and ray tracing allow nonparaxial propagation to be considered, but usually they do not include diffraction [21]. Some aspects of diffraction can be caught by associating phase information, in terms of the optical path length, to each ray. Douglas et al. simulated the operation of an interferometer by ad hoc assumptions that allowed computability to be simplified, but this waived the exact description of diffraction [22]. Improvements were made by Riechert et al., but at the cost of sampling issues and relatively large computation times [23].

There exist several commercial software packages, where the segments of an optical system are tackled by different propagation techniques, e.g., ZEMAX, GLAD, and LightTrans VirtualLab. The latter is described in detail by Wyrowski and Kuhn [24]. While free-space propagation is handled by Fourier optics methods, segments with interfaces are modeled either by the thin lens approximation or the geometrical optics field tracing. At the interfaces, the input field is locally decomposed, assuming local plane waves propagating in a direction perpendicular to the wavefront. These plane waves define the launch directions of ray tubes, which are traced through the subsequent subdomain, taking the amplitude changes into account by Fresnel equations and actual geometrical convergence or divergence. However, because they do not interact at the output plane, diffraction is not modeled. For this reason, LightTrans VirtualLab does not continuously model diffraction throughout the optical system; therefore, it is not suited for our purposes.

A different approach, which is common in optical design and is the basis of at least two software packages, was proposed by Greynolds in 1985 [25]. This approach is based on the input-field decomposition into Gaussian beams [26], which, in turn, are represented by base rays along with a set of characteristic parabasal rays. After tracing the base and parabasal rays, the beams are recombined to form the output field [27,28]. However, Gaussian beams do not form a complete orthogonal set of basis functions; the decomposition is not unique and produces artifacts in the output field. Moreover, this method is reliable only if the parabasal rays are paraxial everywhere, and it is not clear how, for instance, coma can be modeled by basis beams that, at most, describe astigmatism.

The method called stable aggregate of flexible elements (SAFE) relies on a similar but improved concept, which allows the estimation of its own accuracy. Alas, so far it has not yet been extended to 3D problems [29,30]. Recently, a quantum-mechanical method has been proposed, by using Feynman path integrals and stationary-phase approximation to link geometric and scalar wave optics [31,32].

When the measurements aim at subnanometer accuracy levels, it is possible that the description of interferometer operation by paraxial and scalar approximations is not sufficient. The software package $\text{Code}V$ exploits a method called beam synthesis propagation that might meet the case. However, to the authors’ knowledge, no detailed description of this method is available to carry out an uncertainty analysis when applied to compensate for phase delays due to diffraction in state-of-the-art dimensional metrology. Extensive investigations of the propagation of electromagnetic waves and diffraction can be found in the literature of the radio wave and microwave communities [33], where significant effort was put into the efficient calculation of large-scale antenna fields as well as of radar cross sections and where computationally efficient ray tracing techniques were developed [34 –43].

In order to provide an efficient, vectorial, and nonparaxial model of the operation of laser interferometers, as applied in length and displacement metrology, and to calculate the relevant diffraction corrections, we developed a ray-based method to integrate diffraction integrals (VRBDI, vectorial ray-based diffraction integral). An early scalar version of this method is described in [44]. Although it was independently developed, the main concept is the same as proposed by Ling et al., who, in 1989, published a comprehensive ray method, later known as SBR-PO for shooting and bouncing rays and physical optics [39,43].

Ling et al. approximated the ray-based field on the detector grid locally by plane waves [39], although in 1988, they described the far-field contribution of a ray tube by an approximate solution of an integral taken on the ray-tube wavefront [40]. In this paper, we calculate the field on the detector grid by using a ray aiming approach together with a local representation of the wavefront, which is based on matrix optics [16 –18,45] and differential ray tracing [34,46]. In [39], the input field was assumed to be a single plane wave. Therefore, a single set of parallel rays entered the investigated optical system. Here, we decompose the input field into numerous components, which are either spherical or plane waves and are represented by sets of divergent or parallel ray tubes, respectively. The output field is then obtained by an integral superposition of the traced components. Ling et al. applied the phase matching technique to describe the local curvature and, hence, the divergence assigned to a single ray [34]. In [39], it is noted that an alternative description by differential ray tubes, which, in turn, consist of base and parabasal rays, does not deliver a correct description of the Gouy phase as, potentially, too many ray path crossings that are difficult to detect could not be considered. In Appendix B of this paper, it is shown how the correct Gouy phase can indeed be obtained from differential ray tubes without the tracking of ray path crossings. In addition to the perfectly reflecting surfaces considered in [39], we extended the material equations to include transparent dielectrics. Finally, much effort was put into a stepwise, i.e., surface to surface, integral method, for which numerical tests of energy conservation were carried out in order to check physical consistency.

The VRBDI method is described in Section 3. In Section 2, the exact treatment of vector diffraction theory in linear, homogeneous, and isotropic dielectrics is reviewed to set the basis for the VRBDI and a touchstone of its performance. Section 4 compares this more rigorous method to the VRBDI methods by applying them to a simple toy model. Finally, to illustrate the applicability of the VRBDI, it is used in Section 5 to simulate an interferometer setup. Furthermore, these results are compared to an earlier paper where paraxial methods were utilized to describe the same setup.

2. VECTORIAL DIFFRACTION THEORY

A. Derivation of Vectorial Diffraction Integrals

For the following treatment, the same boundary conditions as, e.g., in Chen et al. [47] are assumed. The derivation of the diffraction integrals follows the one found in a book by Smith [48]. The harmonic time dependence $\exp (\mathrm{i}\omega t)$ is consequently omitted.

A monochromatic continuous electromagnetic wave in a linear, homogeneous, isotropic, and transparent dielectric is fully characterized by knowing the components of either the electric field $\mathbf{E}$ or the magnetic field $\mathbf{H}$ tangential to a plane. For the sake of simplicity, let this plane be the $x$, $y$ plane of a Cartesian coordinate system and let $E_x(x_0,y_0)$ as well as $E_y(x_0,y_0)$ be initially known. Then the angular plane-wave spectrum can be obtained by Fourier transform [19]:

The following derivation is done solely for the electric field, as it is fully analogous to the one for the magnetic field. Setting Eq. (2) back into Eq. (4) and using the versors $\widehat{\mathbf{x}}$, $\widehat{\mathbf{y}}$, and $\widehat{\mathbf{z}}$ of the coordinate system, one obtains

With the short notation $\mathbf{E}({\mathbf{P}}_j)={\mathbf{E}}_j$, the versor $\widehat{\mathbf{r}}=\mathbf{r}/r$, and, after pulling the curl operator into the integral, one obtains

Because we are interested in the simulation of physical light beams, we also request that the field on $S_0$ is a section of a field satisfying the Maxwell equations. Practically, this condition can only be met asymptotically, because a finite and discrete calculation window where the initial field is calculated by Eqs. (1 )–(4) from two given components involves a cut at the border of the calculation window and finite sampling resolution. However, with sufficient resolution and a large enough computation window, which encloses a square integrable field, the error can be made negligible.

A necessary but not sufficient check of physical correctness can be done by applying the law of energy conservation. For a time-averaged continuous wave, this implies conservation of power. The respective numerical tests will be treated after the next subsection. For these tests, we also require that the initial field is a restriction of a solution of the Maxwell equations to $S_0$ in order to calculate the initial power.

B. Propagation through an Interface

Although the boundary conditions for the initial field applied to Eqs. (12) and (13) require the definition on a plane, arbitrarily curved interfaces are assumed here. The impact of this approximation is checked later by numerical tests.

An alternative method, based on the Stratton–Chu diffraction integrals [54], was introduced by Guha and Gillen [55], but it lacked the possibility to test power conservation. We compared numerically the Stratton–Chu diffraction integrals with Eqs. (12) and (13) by checking power conservation. We found that Eqs. (12) and (13) ensure power conservation several orders of magnitude more accurately than the Stratton–Chu diffraction integrals. A thorough report of this issue will be given in a separate paper; for now, we do not see advantages in using the more complex Stratton–Chu diffraction integrals in our kind of propagation problems.

A sketch of the situation is shown in Fig. 1. For each versor $\widehat{\mathbf{r}}$, the respective reflected versor ${\widehat{\mathbf{r}}}_{\mathrm{r}}$ and refracted one ${\widehat{\mathbf{r}}}_{\mathrm{t}}$ are calculated by assuming for an infinitesimally small patch around the intersection point local planarity for the wavefront and interface. Furthermore, in transparent isotropic media, the Poynting vector, which describes the transport of energy by electromagnetic radiation [49,50],

 

Fig. 1. Propagation through interface. The versors for reflection ${\widehat{\mathbf{r}}}_r$ and refraction ${\widehat{\mathbf{r}}}_t$ are calculated for each versor $\widehat{\mathbf{r}}$.

Download Full Size | PDF

It is useful to define the reference frames $(\widehat{\mathit{\xi }},\widehat{\mathit{\eta }},\widehat{\mathbf{r}})$, $({\widehat{\mathit{\xi }}}_{\mathrm{r}},\widehat{\mathit{\eta }},{\widehat{\mathbf{r}}}_{\mathrm{r}})$, and $({\widehat{\mathit{\xi }}}_{\mathrm{t}},\widehat{\mathit{\eta }},{\widehat{\mathbf{r}}}_{\mathrm{t}})$, with

The dynamic properties of plane waves lead to the Fresnel equations [49,50]:

It is convenient to express the field calculations at the interface by differentials and perform the integration afterward. With the quantities defined above, one can write

The next step is the propagation from $S_1$ to $S_2$, where Eqs. (12) and (13) can be used by replacing ${\mathbf{E}}_0$, ${\mathbf{H}}_0$, $k$ and $\lambda$ with ${\mathbf{E}}_{1,\mathrm{t}}$, ${\mathbf{H}}_{1,\mathrm{t}}$, $k_0{n}_2$ and ${\lambda }_0/n_2$, respectively. If an integration is performed on a curved interface, the surface element $\mathrm{d}A$ must be changed. For the $(x,y)$ plane, one has $\mathrm{d}A=\mathrm{d}x\mathrm{d}y$. To define a surface, the $z$ coordinate can be given by the 2D function $z=f(x,y)$. Then, the surface element can be calculated according to $\mathrm{d}{A}_f(x,y)=\sqrt{1+{(\partial f/\partial x)}^2+{(\partial f/\partial y)}^2}\mathrm{d}x\mathrm{d}y$ [56]. For example, a hemisphere is uniquely described by

C. Numerical Tests of Power Conservation

The transport of energy by electromagnetic radiation is described by Eq. (14). Taking the scalar product with a surface normal ${\widehat{\mathbf{N}}}_i$ yields the irradiance,

The simulated toy models are depicted in Fig. 2. In accordance with [47], the input field on surface $S_0$ is defined by

 

Fig. 2. Toy models 1 (top) and 2 (bottom) for the numerical test of power conservation for reflection, refraction, and propagation of electromagnetic fields by vectorial diffraction integrals. An interface separates two transparent, isotropic, homogeneous, and nonmagnetic media with refractive indices $n_1$ and $n_2$. The fixed parameters are listed in Table 1; the variable parameters are given in Table 2.

Download Full Size | PDF

Table 1. Fixed Simulation Parameters for the Toy Models Shown in Fig. 2

View Table | View all tables in this article

Table 2. Variable Simulation Parameters for Toy Models Shown in Fig. 2 ^a

View Table | View all tables in this article

The field values on $S_2$ from the direct propagation are presented in Fig. 3. For sake of brevity, in the main article only the electric field components are shown. The corresponding magnetic components can be found in Appendix C. It can be seen that the $y$ component of the electric field, $E_{2^{\prime },y}$, vanishes everywhere.

 

Fig. 3. Directly from $S_0$ to $S_2$ [Fig. 2 (top), Test 1 in Tables 1 and 2] propagated electric field components of a linearly $x$-polarized Gaussian beam by use of vectorial diffraction integrals. The viewing angle is perpendicular to $S_2$. The corresponding magnetic components can be found in Appendix C.

Download Full Size | PDF

The field that is propagated over $S_1$ does not vanish exactly, as can be deduced from the relative deviations plotted in Fig. 4. However, the residual values are small and likely numerical artifacts. Because the directly propagated $E_{2^{\prime },y}$ is exactly zero, the relative deviation for $E_{2,y}$, as it is defined in Fig. 4, is formally up to 1. However, since, in absolute terms $|E_{2,y}|<4·{10}^{-15}{\mathrm{Vm}}^{-1}$, it is of no physical relevance. The irradiance from the directly propagated beam and its relative deviation to the indirectly propagated one is presented in Fig. 5. One can see that the size of the calculation window is indeed amply chosen. The peak-to-valley (PV) deviation of the irradiances is found to be $\mathrm{\Delta }{I}_{\mathrm{PV}}=[\max (I_{2^{\prime }}-I_2)-\min (I_{2^{\prime }}-I_2)]/\max (I_2)=4.8·{10}^{-13}$. One can state that this first test reveals the magnitude of the residual numerical errors because no systematic errors of the theory can be identified.

 

Fig. 4. Relative deviations of directly [$S_0$ to $S_2$, Fig. 2 (top), Test 1 in Tables 1 and 2] propagated electric field components of a linearly $x$-polarized Gaussian beam by use of vectorial diffraction integrals from indirectly propagated ones ($S_0$ via $S_1$ to $S_2$). The viewing angle is perpendicular to $S_2$. The large relative deviation for the $y$-component results from the fact that the $y$-component of the directly propagated field, i.e., $E_{2^{\prime },y}$, is exactly zero. The corresponding magnetic components can be found in Appendix C.

Download Full Size | PDF

 

Fig. 5. Irradiance (left) and its relative deviation of the directly [$S_0$ to $S_2$, Fig. 2 (top), Test 1 in Tables 1 and 2] propagated linearly $x$-polarized Gaussian beam by use of vectorial diffraction integrals from the indirectly propagated one ($S_0$ via $S_1$ to $S_2$). The viewing angle is perpendicular to $S_2$.

Download Full Size | PDF

For Test 2, the tilted plane is replaced by a spherical surface while again $n_1=n_2$. As mentioned earlier, in this case the found integrals in Eqs. (12) and (13) can only be approximations whose viability is tested now. The relative deviation of the field components is given in Fig. 6. Again, the $y$ component of the electric field is nonzero. Its absolute value as well as the relative deviation of the irradiance from the directly propagated one is shown in Fig. 7.

 

Fig. 6. Relative deviations of directly [$S_0$ to $S_2$, Fig. 2 (bottom), Test 2 in Tables 1 and 2] propagated electric field components of a linearly $x$-polarized Gaussian beam by use of vectorial diffraction integrals from indirectly propagated ones ($S_0$ via $S_1$ to $S_2$). The viewing angle is perpendicular to $S_2$. The large relative deviation for the $y$-component results from the fact that the $y$ component of the directly propagated field, i.e., $E_{2^{\prime },y}$, is exactly zero. The corresponding magnetic components can be found in Appendix C.

Download Full Size | PDF

 

Fig. 7. Left: from $S_0$ via $S_1$ to $S_2$ [Fig. 2 (bottom), Test 2 in Tables 1 and 2] propagated $y$ component of the electric field of a linearly $x$-polarized Gaussian beam by use of vectorial diffraction integrals. Right: relative deviation of the irradiance between the directly from $S_0$ to $S_2$ and the indirectly ($S_0$ via $S_1$ to $S_2$) propagated beam. The viewing angle is perpendicular to $S_2$.

Download Full Size | PDF

Now, with $S_1$ curved, the deviations can no longer be classified as numerical artifacts. These regular patterns of larger magnitude must be caused by systematic effects. By inspection of Eq. (12) and comparison to Test 1, it turns out that two reasons can be recognized as necessary and sufficient to explain at least the nonzero $y$ component shown in Fig. 7 (left):

The first reason is also true for the arbitrarily tilted plane. However, the second, which enables the (almost) utter cancellation, is not. For the sphere case, cancellation is not perfect due to the different normal orientations; thus, the residual regular pattern is allowed (Fig. 7, left). The regular deviation patterns for the other components are likely caused in an analogue way, i.e., by parasitic contributions, which did not cancel out because of the different normal directions and the vectorial nature of the integral kernels. However, the impact on the irradiance is still quite small: $\mathrm{\Delta }{I}_{\mathrm{PV}}=9.1·{10}^{-13}$ (Fig. 7, right). Since the residual irradiance error is quite symmetrically modulated around zero, upon integration the positive and negative parts almost cancel each other out. Therefore, the respective relative error of power conservation ${\delta }_{2^{\prime },2}$ is even smaller:

Thus, in rotationally symmetric setups, relying on power conservation alone can potentially lull us into a false sense of security. This symmetry can be broken when $S_2$ is tilted, as is the case for the remaining tests.

Table 3 lists the power conservation errors for the Tests 3 (plane interface) and 4 (curved interface). The relevant test parameters are given in Table 2. Now, the refractive indices $n_1$ and $n_2$ are different and $S_1$ becomes a real interface. The largest error appears in the calculation of reflected and refracted power and is more pronounced with a curved interface. Strangely, the error for the propagation step to or from the curved interface is smaller than for the tilted plane case by 1 order of magnitude. This may be explained by the symmetry effect mentioned above. On the other hand, these values are already quite small and in the range of the numerical limit. Therefore, it is questionable if they are interpretable this way.

Table 3. Numerical Test of Power Conservation for Propagation through a Real Interface between Two Transparent, Isotropic, Homogeneous, and Nonmagnetic Media^a

View Table | View all tables in this article

To check the influence of the sample intervals on the error, Test 4 is repeated for different sampling resolutions. Table 4 shows the results together with the needed calculation times on a personal computer with 64 bit MATLAB, 24 GB RAM and two Intel Xeon X5680 @3.33 GHz, i.e., altogether 12 cores, each working in parallel at 100% workload by use of MATLAB’s Parallel Processing Toolbox. Also, all computation times given later refer to the above hardware and software conditions. While ${\delta }_{\mathrm{1,0}}$ and ${\delta }_{\mathrm{2,1}}$ always are in the range of the numerical limit, $|{\delta }_1|$ decreases with increasing resolution to $1.3·{10}^{-9}$, which is likely a lucky value beyond the reachable numerical limit for this type of calculation in this setup, because the values for ${444}^2$ and ${555}^2$ pixels are larger than for ${333}^2$ pixels. Therefore, this numerical error appears to be in the range of $\approx 5·{10}^{-9}$, and it is reached at a certain resolution threshold, somewhere between ${255}^2$ and ${333}^2$ pixels.

Table 4. Numerical Test of Power Conservation for Different Sampling Resolutions for Propagation through an Interface between Two Transparent, Isotropic, Homogeneous, and Nonmagnetic Media^a

View Table | View all tables in this article

The dependence on the sphere radius $R$ is also checked. In order to keep the calculation time low, the sampling resolution was chosen to be $199\times 199$. The results are listed in Table 5.

Table 5. Numerical Test of Power Conservation for Different Radii of Curvature of an Interface between Two Transparent, Isotropic, Homogeneous, and Nonmagnetic Media^a

View Table | View all tables in this article

Now, ${\delta }_1$ is always $\approx -4·{10}^{-8}$, while ${\delta }_{1,0}$ and ${\delta }_{\mathrm{2,1}}$ are again governed by numerical noise.

The purpose of the last numerical test in this section is twofold. First, it checks whether the field in a focal plane can be calculated. Second, it delivers error values for a stronger refractive index contrast. The refractive index values are changed to $n_1=1.05$ and $n_2=3.17$. Then, the paraxial focus is obtained at $z_2$ for $R=z_1{z}_2(n_2-n_1)/(z_1{n}_2+z_2{n}_1)=20.113852\mathrm{mm}$. The results are shown in Table 6 for various sampling resolutions.

Table 6. Numerical Test of Power Conservation for Different Sampling Resolutions for Propagation in a Focusing Toy Model with Stronger Refractive Index Contrast $n_2/n_1\approx 3$ ^a

View Table | View all tables in this article

The values for ${\delta }_{\mathrm{1,0}}$ and ${\delta }_{\mathrm{2,1}}$ are again in the range of the numerical limit. Therefore, it would be daring to obtain significant trends from three data points. Comparison of the column at ${199}^2$ pixels to the $R=20\mathrm{mm}$ column of Table 5 reveals for $|{\delta }_1|$ a strong proportionality to the refractive index contrast. Due to this scaling effect, the decreasing of $|{\delta }_1|$ with increasing resolution is now more resolvable. However, again there appears to be a limit in the range of ${10}^{-9}$. As mentioned earlier, the Fresnel equations are only exact for infinite plane waves at infinite plane boundaries. Thus, the violation of this condition may reveal itself here.

It could be shown that the above stepwise method describes the transport of radiation power through an interface between two transparent, isotropic, homogeneous, and nonmagnetic media with reasonable accuracy for quite large wavelengths, small curvature radii, and strong refractive index contrasts. However, it should be noted here that this method is not immune to undersampling, which can always prevent proper results. Due to the sampling requirements of the $\exp (-\mathrm{i}kr)$ term in Eqs. (12) and (13), which have already been extensively discussed in the literature, e.g., see [57], the simulation of a macroscopic optical system with many (close-by) interfaces can easily inflate the computation time to astronomic scale. This consideration has led to a different approach that is not limited by the number and distance of interfaces in an optical system. This method is presented in the next section. Fortunately, it can be compared to the stepwise method in simple toy models, such as the ones shown above. This is done after the next section.

3. VECTORIAL RAY-BASED DIFFRACTION INTEGRAL

In this section, the vectorial ray-based diffraction integral (VRBDI) is presented. In order to increase readability, it is structured into several subsections. First, a brief introduction into the used geometrical ray tracing is given. Then, a method to get close to specific sampling points behind an optical system is briefly presented while the details can be found in Appendix A.

Afterward, the extension of ray tracing to gather all relevant information for the VRBDI is introduced while the details of the method, i.e., how matrix optics [16 –18,45] and differential ray tracing [34,46] can be used to calculate the field at the sampling grid points, have been moved to Appendix B. The VRBDI itself concludes this section.

A. Ray Tracing

Geometrical ray tracing is a common technique used in optical design software [21,50]. A ray is defined by position ${\mathbf{P}}_0$ and versor $\widehat{\mathbf{r}}$. Then, the question as to how to scale the versor in order to intersect the next surface, i.e., ${\mathbf{P}}_1={\mathbf{P}}_0+l\widehat{\mathbf{r}}$, leads to an equation system that needs to be solved. One has to distinguish between sequential ray tracing, where the next surface is given by the next entry in a list, and nonsequential ray tracing, where the next surface has to be determined for each ray individually. The latter one prevents effective parallel implementation, i.e., the tracing of many rays at once rather than one after another. Therefore, we chose sequential ray tracing in order to implement a parallel computation.

A plane is fully characterized by pivot vector $\mathbf{T}$ and normal ${\widehat{\mathbf{N}}}_z$. Then, the condition $({\mathbf{P}}_0+l\widehat{\mathbf{r}}-\mathbf{T})•{\widehat{\mathbf{N}}}_z=0$ directly yields

When the intersection point is found, the reflected or refracted versor is calculated by Eq. (16) or (17), respectively. Then, the next ray is defined by the intersection point and the respective versor. If a ray intersection lies outside of allowed positions, i.e., outside an aperture, it is marked as terminated. The first and last surfaces are usually planes, the input plane and the detector plane, respectively.

B. Ray Aiming

A ray starting from the input plane of an optical system is defined by position ${\mathbf{P}}_0$ and versor $\widehat{\mathbf{r}}$. There are two possible cases, depending on the chosen decomposition of the input field: which parameter is fixed and which one is free for aiming. For a plane wave decomposition, $\widehat{\mathbf{r}}$ is fixed and the starting point on a plane, which is normal to $\widehat{\mathbf{r}}$ and shares its origin with the input plane, is the free parameter. In contrast, for spherical wave decomposition, ${\mathbf{P}}_0$ is a fixed sampling grid point and $\widehat{\mathbf{r}}$ is free. Both cases are sketched in Fig. 8. Initially, a grid on the aperture plane defines the ray spacing.

 

Fig. 8. Types of used input field decompositions in the ray picture. Only one component is shown, respectively. Initially, a grid on the aperture plane defines the ray spacing. Top: decomposition into plane waves: each component is represented by parallel rays starting from the wave front, which intersects the origin of the input plane, i.e., the plane where the decomposition is done by Eqs. (1) and (2). Bottom: decomposition into spherical waves: each component is represented by divergent rays starting from one of the sampling points where the input field is given.

Download Full Size | PDF

By tracing test rays through the optical system, various starting conditions can be related to positions on the detector plane. A fit by Zernike polynomials [50] is then utilized to obtain a functional relation between the starting condition and the detector position. Then, it is possible to aim for specific points on the detector sampling grid. A detailed explanation of the used ray-aiming procedure and the description of a test tool for the reachability of detector positions and quality of the ray aiming is given in Appendix A. Generally, with this kind of ray aiming, the targeted positions on the detector cannot be reached exactly. There is always a smaller or larger displacement between the targeted grid point and the ray intersection. This displacement will also be important for the next subsection.

There are also situations where certain regions in the area of the detector cannot be reached by rays, while in other regions the relation between starting condition and position on the detector plane is not unique, i.e., where the different test rays intersect each other. This is the case at caustics, focal, and image regions.

C. Tracing the Electromagnetic Field

From each aimed ray, i.e., a ray that gets close to an assigned detector pixel and is found by ray aiming, a differential ray tube is constructed at an input position [34,46]. To each ray tube, an electric field contribution is assigned by

D. VRBDI

The VRBDI itself is then given by

 

Fig. 9. Flow chart of the VRBDI algorithm.

Download Full Size | PDF

4. COMPARISON BY A TOY MODEL

In this section, the VRBDI is compared to the stepwise propagation method of the second section. The second toy model from Fig. 2 is chosen, and the parameters are changed according to Table 7. Again, the input field is given by Eq. (30) with $w_0=0.5\mathrm{mm}$, $E_{0,y}=0$, and using Eqs. (1 )–(4) for the other components. It is worth noting that the optimal conditions for the plane wave decomposition are reciprocal to the optimal conditions for the spherical wave decomposition. For instance, in Table 7 the input field for the plane wave decomposition is chosen on a wider spatial mesh (the number of pixels is maintained while the overall dimensions of the mesh are increased), which yields the needed higher spatial frequency resolution. In contrast, for the spherical decomposition, a sufficient spatial resolution is most important. Therefore, a denser spatial mesh is better.

Table 7. Simulation Parameters for the Comparison Between the Vectorial Ray-Based Diffraction Integral and Stepwise Method Based on Vectorial Diffraction Integrals^a

View Table | View all tables in this article

The stepwise method is again checked for power conservation:

From the plane wave decomposition, 1069 components, i.e., plane waves, and from the spherical decomposition, 59,381 points were chosen for the VRBDI. In both cases, a relative threshold, ${10}^{-9}$ for the spherical wave decomposition and ${10}^{-10}$ for the plane wave decomposition, in relation to the strongest component, defines a cut-off level. The following calculation times are obtained for the three methods:

Obviously, choosing lower cut-off levels increases calculation time but, of course, also accuracy. Therefore, a trade-off between accuracy and computation time has to be found. However, most important is the proper choice of the decomposition type. In the presented case, VRBDI-PW is roughly 57 times faster than VRBDI-SW, which corresponds to the number of the respective input components: $59,381/1069\approx 56$. The results are quite comparable, as can be seen in Fig. 10. However, some distinctions are observable: in the VRBDI-PW deviations from the stepwise method, one can see some numerical noise in the low-field regions.

 

Fig. 10. Comparison between different propagation methods by toy model 2 shown in Fig. 2 with the parameters from Table 7. Left: electric field components propagated by a vectorial ray-based diffraction integral where a linearly $x$-polarized Gaussian beam is decomposed by spherical waves (VRBDI-SW). Middle: relative deviation from a stepwise propagation by vectorial diffraction integrals. Right: relative deviation of the electric field components propagated by a VRBDI with plane wave decomposition (VRBDI-PW) from the ones obtained by the stepwise method. The viewing angle is perpendicular to $S_2$. The corresponding magnetic components can be found in Appendix C.

Download Full Size | PDF

The irradiance calculated by VRBDI-SW and its relative deviation from the stepwise method as well as the relative deviation of VRBDI-PW from VRBDI-SW are shown in Fig. 11. The latter is quite small, $\mathrm{\Delta }{I}_{\mathrm{PV}}\approx \pm 6·{10}^{-8}$, and seems to be caused by sampling issues. Therefore, the uncertainty for both is represented by the relative deviation of VRBDI-SW from the stepwise method, which is $\mathrm{\Delta }{I}_{\mathrm{PV}}\approx -6·{10}^{-6}$. From the shape of the deviation, one can deduce that the beam spot generated by the stepwise method is slightly broader. By consideration of power conservation,

 

Fig. 11. Comparison between different propagation methods by toy model 2 shown in Fig. 2 with the parameters from Table 7. Left: irradiance obtained by a vectorial ray-based diffraction integral where a linearly $x$-polarized Gaussian beam is decomposed by spherical waves (VRBDI-SW). Middle: relative deviation from a stepwise propagation by vectorial diffraction integrals. Right: relative deviation of the irradiance obtained by a vectorial diffraction integral with plane wave decomposition (VRBDI-PW) from the one obtained by VRBDI-SW. The viewing angle is perpendicular to $S_2$.

Download Full Size | PDF

 

Fig. 12. Comparison between different propagation methods by toy model 2 shown in Fig. 2 with the parameters from Table 7: relative deviation of the electric field components propagated by a vectorial diffraction integral with spherical wave decomposition (VRBDI-SW) of a linearly $x$-polarized Gaussian beam from the ones propagated by a VRBDI with plane wave decomposition (VRBDI-PW). The viewing angle is perpendicular to $S_2$. The corresponding magnetic components can be found in Appendix C.

Download Full Size | PDF

5. INTERFEROMETER SIMULATION

In order to illustrate the versatility of the VRBDI method, it is applied to the calculation of the diffraction correction for an interferometer, as it was used to measure the diameter of silicon spheres [11 –13]. The parameters of the simulated setup are identical to the ones in [13], where a Gaussian beam-tracing algorithm, known as vectorial general astigmatic Gaussian beam (VGAGB) tracing, was utilized.

As a rule of thumb, when deciding which decomposition type for the input field is preferable, i.e., the one that needs the fewest components, it is suggested to consider two extreme situations:

Following this rule, VRBDI-SW is chosen because the source, a Gaussian beam with waist radius $w_0=0.28562\mathrm{mm}$, can be sufficiently sampled by $31\times 31$ pixels at its waist location, which is 683.6 mm away from the first interferometer component. Furthermore, a relative threshold of ${10}^{-6}$ yields 793 components for each VRBDI run.

The resolution and size of the detectors were chosen to be $255\times 255$ pixels over $8\times 8{\mathrm{mm}}^2$. The interferometer is operated by wavelength shifting. As in [13], the full procedure is simulated here as well but only for the ideal configuration of the interferometer (no Monte Carlo simulations). For each measurement, i.e., for an empty and loaded cavity, indicated as E and L, respectively, seven interferograms $I_{\mathrm{E},\mathrm{L}}(\lambda )$, one for each wavelength, are captured and fed into the phase retrieval formula of [13]. For each wavelength, three beam paths are considered: corresponding to zero, one, and two round trips in the empty or loaded gap. Therefore, $2\times 7\times 3=42$ VRBDI runs have to be performed. With Eqs. (14) and (28), the interferograms $I_{\mathrm{E},\mathrm{L}}(\lambda )$ are obtained from the respective sum of the (complex-valued) fields from each path. In Fig. 13, the phase biases $\mathrm{\Delta }{\varphi }_{\mathrm{E},\mathrm{L}}$ are plotted. We refer to bias as the difference between measured value and correct value of a quantity. Therefore, a correction to a measured value is the respective negative bias. The difference to the VGAGB calculation, $\delta \mathrm{\Delta }{\varphi }_{\mathrm{E},\mathrm{L}}=\mathrm{\Delta }{\varphi }_{\mathrm{E},\mathrm{L},\mathrm{VRBDI}}-\mathrm{\Delta }{\varphi }_{\mathrm{E},\mathrm{L},\mathrm{VGAGB}}$, is shown in Fig. 14. Table 8 compares the resulting measurement biases for empty gap, loaded gap, and resulting sphere diameter from VRBDI-SW and VGAGB. These values are calculated directly from $\delta \mathrm{\Delta }{\varphi }_{\mathrm{E},\mathrm{L}}(\mathrm{0,0})$. Since an odd sampling was chosen, the central value already exists and needs not to be interpolated. The agreement is quite remarkable. By inspection of Fig. 14, one can see that this is only true in a relatively small central area, i.e., the paraxial regime, where the VGAGB method can be considered to be accurate. It should be noted here that the meaning of paraxial and nonparaxial depends on the desired accuracy. On that note, nonparaxial means that the paraxial approximation can no longer be applied because the desired accuracy is no longer met. For a metrologist working in the field of ultraprecision interferometry, the paraxial approximation can potentially be much more restrictive than commonly assumed. However, it is shown here that it is well justified for the desired accuracy at the extremal position of the interferogram phases in the described interferometer, i.e., the location where the interferograms obtained from this setup are usually evaluated. Dividing the calculation time of 16.6 h by 42 yields a mean calculation time per run of ca. 24 min; further division by 793 results in ca. 1.8 s per input field component.

 

Fig. 13. Phase bias for empty and loaded cavity in the ideal configuration of the interferometer described in [11 –13] but calculated by the vectorial ray-based diffraction integral with spherical wave decomposition (VRBDI-SW).

Download Full Size | PDF

 

Fig. 14. Difference of the phase bias for empty (top) and loaded cavity in the ideal configuration of the interferometer described in [11 –13] between the calculation by a vectorial ray-based diffraction integral with spherical wave decomposition (VRBDI-SW) and a vectorial general astigmatic Gaussian beam-tracing algorithm (VGAGB) [13].

Download Full Size | PDF

Table 8. Measurement Bias for Empty Gap $\mathrm{\Delta }{D}_{\mathrm{E}}$, Loaded Gap $\mathrm{\Delta }{D}_{\mathrm{L}}$, and Resulting Sphere Diameter $\mathrm{\Delta }D$ ^a

View Table | View all tables in this article

6. DISCUSSION

The stepwise method, as it is introduced in the second section, represents to our knowledge the most rigorous and still practically computable propagation method applicable to the optical interface problem, where it can be successfully used to test other propagation methods such as the VRBDI.

To provide a more physically motivated and intuitive interpretation, we refer to the Huygens–Feynman–Fresnel principle [31,32]. Provided that the Huygens–Fresnel principle and the Maxwell equations [48 –50] are equivalent to the quantum-mechanical description of a single photon, the stepwise method solves the path integral problem for each segment of the optical system individually. It constitutes the full brute-force approach: for every point on the target surface, the coherent sum of all wavelets emerging from all possible input positions is calculated numerically, although spatial discretization is unavoidable.

Following this picture, the VRBDI would be an approximation whose main assumption is that the light interaction at the interfaces of an optical system can be approximately described deterministically. In this sense, the traced rays represent the different possibilities that remain when the influence of an aperture on the photon’s momentum uncertainty is neglected [59]:

More generally, any confinement of photon trajectories should increase the momentum uncertainty. This could explain why focal regions cannot be described by deterministic and straight photon trajectories. The question as to whether it is possible to properly describe the situation by curved trajectories goes beyond the scope of this paper [60].

Beside these fundamental considerations, there are also technical issues that prevent the proper usage of the VRBDI in focal regions. However, by using Eqs. (12) and (13) or, alternatively, Eqs. (1 )–(4), one can easily close this reachability gap. Alternatively, it seems natural to assume that proper usage of respective angle-to-angle and point-to-angle eikonals [58] could also solve this issue. However, this has not been attempted by the authors so far.

For proper handling of small apertures, an intermediate field has to be sampled on the apertured surface. Approximately, by assuming Kirchhoff boundary conditions [50], the field in the aperture can simply be obtained by multiplication to the respective aperture function. Then, further propagation can be started from there.

7. CONCLUSION

In order to simulate the operation of laser interferometers, a vectorial ray-based diffraction integral method has been developed. The tangential electric components of a continuous monochromatic wave field on a plane, which is immersed in a transparent nonmagnetic homogeneous medium, are necessary and sufficient to obtain the other electromagnetic components by use of the angular plane wave spectrum. From the latter, vectorial diffraction integrals are deduced, which satisfy Maxwell and Helmholtz equations. By exploitation of the kinematic and dynamic boundary conditions of the Maxwell equations, a stepwise surface-to-surface propagation method is established, which is still a good approximation when the interface between two linear, transparent, nonmagnetic, and homogeneous media is curved.

For application of the VRBDI method, an input field for a general optical system is decomposed either into spherical waves, i.e., point sources, or plane waves. By means of ray aiming, an input position of a point source or, alternatively, an input angle of a plane wave can be related to a certain sampling position on a detector plane at the end of this optical system. Differential ray tubes, each consisting of an aimed base ray and four close, i.e., paraxial, parabasal rays, are launched and traced in parallel through the optical system to their respective sampling point on the detector plane. The optical path length and a ray matrix are obtained for each ray tube. By use of the respective paraxial eikonals and further adaptation of the diffraction integrals, the electromagnetic field at all sampling points is calculated and summed to all the other traced contributions from the decomposed input field. For the location of the detector plane, caustics as well as focal and image regions have to be avoided. However, propagation to these regions can be carried out by common free-space propagation methods.

Direct comparison of the VRBDI with the stepwise method in a simple curved interface toy model using power conservation as a criterion reveals that the stepwise method is most exact, but the VRBDI is a good approximation when the apertures are much larger than the wavelength. In the chosen toy model, the remaining error of the simulated irradiance is at the parts per million level, which should be almost unresolvable by experiment. By proper choice of the input field decomposition, the VRBDI is already faster than the stepwise method. Improper choice can lead to longer computation times, though, albeit this drawback applies only to simple systems, which can only just be modeled by the stepwise method. Macroscopic optical systems with many (closely spaced) surfaces cannot be modeled with the stepwise method in reasonable computation times.

The simulation of an interferometer, which was previously described by a paraxial Gaussian beam tracing method, showed very good agreement in the paraxial regime, whereas deviations between the VRBDI and the paraxial method at more distant locations from the optical axis already indicate the expectable failure of the paraxial method in the nonparaxial regime.

APPENDIX A: DETAILED DESCRIPTION OF RAY AIMING AND TEST TOOL

For the sake of brevity, only the ray aiming for the spherical wave decomposition case is explained here. Then, ${\mathbf{P}}_0$ is a fixed sampling grid point and $\widehat{\mathbf{r}}$ is free. More exactly, the continuous independent variables are ${\widehat{r}}_x$ and ${\widehat{r}}_y$ when ${\widehat{r}}_z=\sqrt{1-{\widehat{r}}_x^2-{\widehat{r}}_y^2}$ is chosen to be dependent. When the intersection point ${\mathbf{P}}_i$ on the detector plane is expressed in the local detector frame defined by the orthonormal set $({\widehat{\mathbf{N}}}_x,{\widehat{\mathbf{N}}}_y,{\widehat{\mathbf{N}}}_z)$ and pivot vector $\mathbf{T}$,

(A1)$$\begin{gathered} x^{\prime }=({\mathbf{P}}_i-\mathbf{T})•{\widehat{\mathbf{N}}}_x, \\ {y}^{\prime }=({\mathbf{P}}_i-\mathbf{T})•{\widehat{\mathbf{N}}}_y, \end{gathered}$$

${\widehat{r}}_x$ and ${\widehat{r}}_y$ can formally be expressed as functions of these coordinates:

(A2)$$\begin{gathered} {\widehat{r}}_x={\widehat{r}}_x(x^{\prime },y^{\prime }), \\ {\widehat{r}}_y={\widehat{r}}_y(x^{\prime },y^{\prime }). \end{gathered}$$

Now, (A2) can be approximated by proper interpolation functions. For this purpose, we chose Zernike polynomials [50]. The definition space of the Zernike polynomials requires

(A3)$$\begin{gathered} \rho =\frac{\sqrt{{(x^{\prime })}^2+{(y^{\prime })}^2}}{{\rho }_{\max }}, \\ \phi =\mathrm{atan}2(y^{\prime },x^{\prime }). \end{gathered}$$

The largest reachable radius ${\rho }_{\max }$ is used to confine $\rho$ to [0,1]. Our goal is to write $\widehat{\mathbf{r}}$ as $\mathrm{\Sigma }{Z}_j(\rho ,\varphi )a_j$. The following equation system must then be solved:

(A4)$$[\mathbf{Z}]\mathbf{a}=\left[\begin{array}{ccc}{Z}_0^0({\rho }_1,{\phi }_1)& Z_1^{-1}({\rho }_1,{\phi }_1)& \dots \\ \dots & \dots & \dots \\ Z_0^0({\rho }_{\nu },{\phi }_{\nu })& Z_1^{-1}({\rho }_{\nu },{\phi }_{\nu })& \dots \end{array}\begin{array}{c}{Z}_u^u({\rho }_1,{\phi }_1)\\ \dots \\ Z_u^u({\rho }_{\nu },{\phi }_{\nu })\end{array}\right]\left[\begin{array}{c}{a}_1\\ \dots \\ a_{\zeta }\end{array}\right]=\left[\begin{array}{c}F({\rho }_1,{\phi }_1)\\ \dots \\ F({\rho }_{\nu },{\phi }_{\nu })\end{array}\right]=\mathbf{F},$$

where $F={\widehat{r}}_x$ or ${\widehat{r}}_y$. Thus, $\mathbf{F}$ must contain uncorrelated values, e.g., randomly chosen or discretely sampled, for ${\widehat{r}}_x$ or ${\widehat{r}}_y$, respectively. In Eq. (A4), the number of Zernike polynomials $\zeta$, i.e., the number of columns in $[\mathbf{Z}]$, should match the number of detected ray intersections $\nu$, i.e., the number of rows in $[\mathbf{Z}]$. However, a least-squares solution for $\nu \gg \zeta$ is readily obtained by [61]

(A5)$$\mathbf{a}={({[\mathbf{Z}]}^{\mathrm{T}}[\mathbf{Z}])}^{-1}{[\mathbf{Z}]}^{\mathrm{T}}\mathbf{F}.$$

When $\mathbf{a}$ is determined, the sought versor components to reach the sampling grid points are obtained directly from

(A6)$${\mathbf{F}}_{\mathrm{sgp}}={[\mathbf{Z}]}_{\mathrm{sgp}}\mathbf{a},$$

where ${[\mathbf{Z}]}_{\mathrm{sgp}}$ is constructed by use of the local coordinates $(x_{\mathrm{sgp}},y_{\mathrm{sgp}})$ of the chosen detector sampling grid points and (A3). Then, with ${\widehat{r}}_{x,\mathrm{sgp}}$ and ${\widehat{r}}_{y,\mathrm{sgp}}$, new rays are traced that should “hit” the sampling grid points.

In order to evaluate the quality of the ray aiming and the reachability of the detector sampling grid points, a test tool has been developed. Figure 15 illustrates the test method. To the aperture plane a sampling grid is assigned, holding either zero for blocking or one for passing a ray. From each source point, whose intensity corresponds to the numbering in Fig. 15 (left), a ray is directed onto each of the nonblocking aperture grid points and traced to the detector. Then, a simple incoherent additive sampling is done on the detector grid.

 

Fig. 15. Left: distribution of chosen source points on the input plane for the test of reachability of the detector sampling grid positions of an optical system. The numbering corresponds to the respective source intensity. Right: example of an optical system under test (only marginal rays from a single source point are shown).

Download Full Size | PDF

For the central source point, i.e., point 9 in Fig. 15, ray aiming is carried out and the positioning error is evaluated according to

(A7)$${\delta }_{\mathrm{ra}}=\frac{\sqrt{{(x^{\prime }-x_{\mathrm{sgp}})}^2+{(y^{\prime }-y_{\mathrm{sgp}})}^2}}{\sqrt{\mathrm{\Delta }{x}_{\mathrm{sgp}}^2+\mathrm{\Delta }{y}_{\mathrm{sgp}}^2}},$$

where $(x^{\prime },y^{\prime })$ are already the ray intersection coordinates and $\mathrm{\Delta }{x}_{\mathrm{sgp}}$ and $\mathrm{\Delta }{y}_{\mathrm{sgp}}$ are the sampling intervals of the detector-plane grid. Therefore, successful ray aiming should at least result in ${\delta }_{\mathrm{ra}}<0.5$ but preferably into a very small number thereof. Thus, when ${\delta }_{\mathrm{ra}}\ge 0.5$ the ray aiming is poor.

In Fig. 16, the result for a simple imaging system [distance $g=200\mathrm{mm}$ from input ($b_{0,x}=b_{0,y}=5\mathrm{mm}$) to lens middle plane], consisting of a biconvex spherical lens (focal length $f_1=50\mathrm{mm}$, diameter $d_{\mathrm{L}1}=40\mathrm{mm}$, thickness $t_1=10\mathrm{mm}$) and a circular aperture (diameter $d_{\mathrm{A}1}=4\mathrm{mm}$) directly in front of the lens, is shown at $\mathrm{\Delta }z\approx 51.7\mathrm{mm}$ behind the rear lens surface. One can see that the nine source points (Fig. 15, left) yield separated spots. Thus, there is almost no overlap. The values for $b_{x,y}$ and ${\overline{b}}_{x,y}$ given in Fig. 16 are quite different. For a good overlap, these values should become almost equal. The ray aiming error for the central source point is quite small. Figure 17 shows the same for different image distances and a larger aperture $d_{\mathrm{A}2}=30\mathrm{mm}$. Directly behind the lens ($\mathrm{\Delta }z=1\mathrm{mm}$) and far away ($\mathrm{\Delta }z=1\mathrm{m}$) overlap and ray aiming are good, although the ray aiming for the larger aperture is not as good as for the small aperture. Close to the image plane at $\mathrm{\Delta }z\approx 61.7\mathrm{mm}$, there is insufficient to no overlap, and the ray aiming error becomes very large.

 

Fig. 16. Test of overlap and ray aiming for a simple imaging system with a 4 mm aperture. Left: additive sampling of the rays from the nine input source points yields nine separated spots on the detector grid. The plotted quantity is arbitrary and corresponds to the ray density times the respective source intensity. The depicted values are $b_{x,y}$, total point spread width including all sources; ${\overline{b}}_{x,y}$, mean of the point spread widths taken from the individual source points; $b$, largest width of the point spread of the central source point, i.e., source point 9 in Fig. 15, which is also the calculation window size of the right picture; $\mathrm{\Delta }z$, distance from the last lens surface apex to the detector plane. Right: the ray aiming error according to Eq. (A7) for the central source point. The viewing angle is perpendicular to the detector plane.

Download Full Size | PDF

 

Fig. 17. Test of overlap and ray aiming for a simple imaging system with a 30 mm aperture. Left: additive sampling of the rays from the nine input source points on the detector grid. The plotted quantity is arbitrary and corresponds to the ray density times the respective source intensity. The depicted values are $b_{x,y}$, total point spread width including all sources; ${\overline{b}}_{x,y}$, mean of the point spread widths taken from the individual source points; $b$, largest width of the point spread of the central source point, i.e., source point 9 in Fig. 15, which is also the calculation window size of the right picture; $\mathrm{\Delta }z$, distance from the last lens surface apex to the detector plane. Right: the ray aiming error according to Eq. (A7) for the central source point. The viewing angle is perpendicular to the detector plane.

Download Full Size | PDF

In order to illustrate the application to a different optical system and to show that the optimal overlap and ray aiming are not always close to or far away from the last lens surface, a second lens ($f_2=100\mathrm{mm}$, $d_{\mathrm{L}2}=40\mathrm{mm}$, $t_2=10\mathrm{mm}$) is placed confocally to the image position of the first lens. The results for various distances are shown in Fig. 18. At $\mathrm{\Delta }z=50\mathrm{mm}$, the ray aiming and overlap are better than at $\mathrm{\Delta }z=5\mathrm{mm}$. For $\mathrm{\Delta }z=5\mathrm{m}$, both are poor.

 

Fig. 18. Test of overlap and ray aiming for a telescope with a 30 mm aperture. Left: additive sampling of the rays from the nine input source points on the detector grid. The plotted quantity is arbitrary and corresponds to the ray density times the respective source intensity. The depicted values are $b_{x,y}$, total point spread width including all sources; ${\overline{b}}_{x,y}$, mean of the point spread widths taken from the individual source points; $b$, largest width of the point spread of the central source point, i.e., source point 9 in Fig. 15, which is also the calculation window size of the right picture; $\mathrm{\Delta }z$, distance from the last lens surface apex to the detector plane. Right: the ray aiming error according to Eq. (A7) for the central source point. The viewing angle is perpendicular to the detector plane.

Download Full Size | PDF

 

Fig. 19. Situation after ray tracing: the points ${\mathbf{Q}}_{\pm x,y}$ on the plane normal to $\widehat{\mathbf{r}}$ have to be found from known intersections $\mathbf{P}$ and ${\mathbf{P}}_{\pm x,y}$ and known versors $\widehat{\mathbf{r}}$ and ${\widehat{\mathbf{r}}}_{\pm x,y}$.

Download Full Size | PDF

At caustics, image, or focal regions, certain locations in the area of the detector cannot be reached by rays, while in other locations rays accumulate in vicinity to certain sampling coordinates, i.e., the overlap is poor. Therefore, the relation between starting condition and position on the detector plane is not unique, i.e., Eq. (A2) is no longer bijective. Use of Eq. (A6) yields distances much larger than $\mathrm{\Delta }{x}_{\mathrm{sgp}}$ or $\mathrm{\Delta }{y}_{\mathrm{sgp}}$ between $(x^{\prime },y^{\prime })$ and $(x_{\mathrm{sgp}},y_{\mathrm{sgp}})$ and ${\delta }_{\mathrm{ra}}$ gets poor.

One can conclude that the test tool enables one to find positions where the features, which are necessary for the calculation of a ray-based diffraction integral, i.e., overlap and hitting of sampling points, are established. Generally, it should now be clear that caustics as well as image and focal regions have to be avoided.

APPENDIX B: DETAILED DESCRIPTION OF FIELD TRACING

From each aimed ray, i.e., a ray that gets close to an assigned detector pixel and is found by ray aiming, a ray tube is constructed at an input position. Such a ray tube consists of the base ray, which is identical to the aimed ray, and four parabasal rays. The ray tubes are generated by either translating the base ray along orthogonal directions or by tilting the base ray around two orthogonal axes. The finite differences, i.e., the translation of magnitude $D$ or the tilt angle $\beta$, should be chosen to be as small as numerically reasonable [46]. The orthogonal directions can be obtained from the arbitrary orthonormal set $({\widehat{\mathbf{e}}}_x,{\widehat{\mathbf{e}}}_y,\widehat{\mathbf{r}})$, which is defined by

(B1)$${\widehat{\mathbf{e}}}_x=\left[\begin{array}{c}{\widehat{r}}_z\\ 0\\ -{\widehat{r}}_x\end{array}\right]\frac{1}{\sqrt{{\widehat{r}}_x^2+{\widehat{r}}_z^2}},{\widehat{\mathbf{e}}}_y=\widehat{\mathbf{r}}\times {\widehat{\mathbf{e}}}_x.$$

Then, the parallel rays start from

(B2)$$\begin{gathered} {\mathbf{P}}_{0,\pm x}={\mathbf{P}}_0\pm D·{\widehat{\mathbf{e}}}_x, \\ {\mathbf{P}}_{0,\pm y}={\mathbf{P}}_0\pm D·{\widehat{\mathbf{e}}}_y, \end{gathered}$$

while the divergent ray directions are given by

(B3)$$\begin{gathered} {\widehat{\mathbf{r}}}_{\pm x}=\widehat{\mathbf{r}}\cos \beta \pm {\widehat{\mathbf{e}}}_x\sin \beta , \\ {\widehat{\mathbf{r}}}_{\pm y}=\widehat{\mathbf{r}}\cos \beta \pm {\widehat{\mathbf{e}}}_y\sin \beta . \end{gathered}$$

Parallel ray initialization is used when a plane wave decomposition of the input field is chosen, while a spherical wave decomposition implies the use of the divergent one.

When the rays of a ray tube are traced through an optical system, they stay sufficiently paraxial to the base ray if $D$ or $\beta$ are chosen sufficiently small, so that the positions and directions of the output rays can be described by matrix optics [45,46]:

(B4)$$\begin{gathered} {\left[\begin{array}{c}\breve{x}\\ \breve{y}\\ {\breve{x}}^{\prime }\\ {\breve{y}}^{\prime }\end{array}\right]}_2=\left[\begin{array}{cc}\mathbf{A}& \mathbf{B}\\ \mathbf{C}& \mathbf{D}\end{array}\right]{\left[\begin{array}{c}\breve{x}\\ \breve{y}\\ {\breve{x}}^{\prime }\\ {\breve{y}}^{\prime }\end{array}\right]}_1 \\ =\left[\begin{array}{cccc}{A}_{11}& A_{12}& B_{11}& B_{12}\\ A_{21}& A_{22}& B_{21}& B_{22}\\ C_{11}& C_{12}& D_{11}& D_{12}\\ C_{21}& C_{22}& D_{21}& D_{22}\end{array}\right]{\left[\begin{array}{c}\breve{x}\\ \breve{y}\\ {\breve{x}}^{\prime }\\ {\breve{y}}^{\prime }\end{array}\right]}_1. \end{gathered}$$

It is worth noting that, in this paper, the reduced slope, e.g., ${\breve{x}}^{\prime }=n(\breve{z})\mathrm{d}\breve{x}(\breve{z})/\mathrm{d}\breve{z}$, is used [18]. The coordinates $\breve{x}$ and $\breve{y}$ live in a plane orthogonal to the base ray $\widehat{\mathbf{r}}$, i.e., the base plane, while their slopes ${\breve{x}}^{\prime }$ and ${\breve{y}}^{\prime }$ are taken along $\widehat{\mathbf{r}}$. Obviously, ${\widehat{\mathbf{e}}}_x$ and ${\widehat{\mathbf{e}}}_y$ can be readily used to initially define the base plane. Thus, the orthonormal set $({\widehat{\mathbf{e}}}_x,{\widehat{\mathbf{e}}}_y,\widehat{\mathbf{r}})$ defines the base ray frame, i.e., the coordinate system that is fixed to the respective base ray. However, when the base ray undergoes refractions or reflections its direction changes. Therefore, at every interface of the optical system, when ${\widehat{\mathbf{r}}}_{\mathrm{r},\mathrm{t}}=\mathbf{M}\widehat{\mathbf{r}}$ the same rotation $O_M$ has to be applied to ${\widehat{\mathbf{e}}}_x$ or ${\widehat{\mathbf{e}}}_y$. For example, ${\widehat{\mathbf{e}}}_{x,\mathrm{r},\mathrm{t}}=\mathbf{M}{\widehat{\mathbf{e}}}_x$, while ${\widehat{\mathbf{e}}}_{y,\mathrm{r},\mathrm{t}}={\widehat{\mathbf{r}}}_{\mathrm{r},\mathrm{t}}\times {\widehat{\mathbf{e}}}_{x,\mathrm{r},\mathrm{t}}$ ensures the right-hand chirality of the orthonormal set. In order to find $O_M$, Rodrigues’ rotation formula [62] is used, which describes how a vector is rotated around an arbitrary rotation axis that is given by a versor $\mathbf{v}$. For the present problem, one can obtain $\mathbf{v}$ by

(B5)$$\mathbf{v}=\frac{\widehat{\mathbf{r}}\times {\widehat{\mathbf{r}}}_{\mathrm{r},\mathrm{t}}}{|\widehat{\mathbf{r}}\times {\widehat{\mathbf{r}}}_{\mathrm{r},\mathrm{t}}|},$$

while the rotation angle is

(B6)$$\theta =\mathrm{arc}\cos (\widehat{\mathbf{r}}•{\widehat{\mathbf{r}}}_{\mathrm{r},\mathrm{t}}).$$

The special case ${\widehat{\mathbf{r}}}_{\mathrm{r},\mathrm{t}}=\widehat{\mathbf{r}}$ would yield division by zero in Eq. (B5), which has to be prevented by explicitly setting $\mathbf{v}=\mathbf{0}$ in this case. There are two possible rotation directions, which have to be considered by probing with $\pm \mathbf{v}$ and taking the solution that minimizes $|\mathbf{M}(\mathbf{v})\widehat{\mathbf{r}}-{\widehat{\mathbf{r}}}_{\mathrm{r},\mathrm{t}}|$, where $\mathbf{M}(\mathbf{v})$ is obtained by Rodrigues’ rotation formula [62]:

(B7)$$\mathbf{M}(\mathbf{v})=\mathbf{I}+\sin \theta {\mathbf{V}}_{\times }+(1-\cos \theta ){\mathbf{V}}_{\times }^2,$$

where $\mathbf{I}$ is the identity matrix and

(B8)$${\mathbf{V}}_{\times }=\left[\begin{array}{ccc}0& -v_z& v_y\\ v_z& 0& -v_x\\ -v_y& v_x& 0\end{array}\right].$$

Extensive matrix multiplication yields

(B9)$$\mathbf{M}(\mathbf{v})=\left[\begin{array}{l}1+(1-\cos \theta )(-v_y^2-v_z^2)\\ v_z\sin \theta +(1-\cos \theta )v_x{v}_y\\ -v_y\sin \theta +(1-\cos \theta )v_x{v}_z\end{array}\begin{array}{l}-v_z\sin \theta +(1-\cos \theta )v_x{v}_y\\ 1+(1-\cos \theta )(-v_x^2-v_z^2)\\ v_x\sin \theta +(1-\cos \theta )v_y{v}_z\end{array}\begin{array}{l}{v}_y\sin \theta +(1-\cos \theta )v_x{v}_z\\ -v_x\sin \theta +(1-\cos \theta )v_y{v}_z\\ 1+(1-\cos \theta )(-v_x^2-v_y^2)\end{array}\right].$$

It is worth noting that, in the special case where $\mathbf{v}=\mathbf{0}$, it follows that $\mathbf{M}(\mathbf{v})=\mathbf{I}$.

While the initial ray vectors are simply given by [46]

(B10)$$\begin{gathered} {\mathbf{p}}_{1,\pm x}=\left[\begin{array}{c}\pm D\\ 0\\ 0\\ 0\end{array}\right],{\mathbf{p}}_{1,\pm y}=\left[\begin{array}{c}0\\ \pm D\\ 0\\ 0\end{array}\right], \\ {\mathbf{q}}_{1,\pm x}=\left[\begin{array}{c}0\\ 0\\ \pm n_0\tan \beta \\ 0\end{array}\right],{\mathbf{q}}_{1,\pm y}=\left[\begin{array}{c}0\\ 0\\ 0\\ \pm n_0\tan \beta \end{array}\right], \end{gathered}$$

where the ${\mathbf{p}}_{1,\pm x,y}$ are needed for the plane wave decomposition and the ${\mathbf{q}}_{1,\pm x,y}$ for the spherical wave decomposition of the input field, the final ray vectors must be obtained from the intersections of the parabasal rays with the $\breve{x}$, $\breve{y}$ plane of the base ray frame where, in turn, the base ray intersects the detector plane. Figure 19 illustrates the problem of finding these intersection points. In analogy to Eq. (31), one can obtain the sought intersection points from any known ray positions and versors by

(B11)$${\mathbf{Q}}_{\pm x,y}=\frac{(\mathbf{P}-{\mathbf{P}}_{\pm x,y})•\widehat{\mathbf{r}}}{{\widehat{\mathbf{r}}}_{\pm x,y}•\widehat{\mathbf{r}}}{\widehat{\mathbf{r}}}_{\pm x,y}+{\mathbf{P}}_{\pm x,y}.$$

The points ${\mathbf{Q}}_{\pm x,y}$ are still 3D vectors in the global frame. The local coordinates on the $\breve{x}$, $\breve{y}$-plane can be obtained from

(B12)$$\begin{gathered} {\breve{x}}_{\pm x,y}=({\mathbf{Q}}_{\pm x,y}-\mathbf{P})•{\widehat{\mathbf{e}}}_x, \\ {\breve{y}}_{\pm x,y}=({\mathbf{Q}}_{\pm x,y}-\mathbf{P})•{\widehat{\mathbf{e}}}_y. \end{gathered}$$

The versors ${\widehat{\mathbf{r}}}_{\pm x,y}$ also have to be transformed into the base ray frame:

(B13)$${\widehat{\mathbf{s}}}_{\pm x,y}=\left[\begin{array}{c}{\widehat{\mathbf{r}}}_{\pm x,y}•{\widehat{\mathbf{e}}}_x\\ {\widehat{\mathbf{r}}}_{\pm x,y}•{\widehat{\mathbf{e}}}_y\\ {\widehat{\mathbf{r}}}_{\pm x,y}•\widehat{\mathbf{r}}\end{array}\right].$$

Now the final ray vectors can be obtained:

(B14)$$\begin{gathered} {\mathbf{p}}_{2,\pm x},{\mathbf{q}}_{2,\pm x}=\left[\begin{array}{c}{\breve{x}}_{\pm x}\\ {\breve{y}}_{\pm x}\\ n_{m-1}{\widehat{s}}_{\breve{x},\pm x}/{\widehat{s}}_{\breve{z},\pm x}\\ n_{m-1}{\widehat{s}}_{\breve{y},\pm x}/{\widehat{s}}_{\breve{z},\pm x}\end{array}\right], \\ {\mathbf{p}}_{2,\pm y},{\mathbf{q}}_{2,\pm y}=\left[\begin{array}{c}{\breve{x}}_{\pm y}\\ {\breve{y}}_{\pm y}\\ n_{m-1}{\widehat{s}}_{\breve{x},\pm y}/{\widehat{s}}_{\breve{z},\pm y}\\ n_{m-1}{\widehat{s}}_{\breve{y},\pm y}/{\widehat{s}}_{\breve{z},\pm y}\end{array}\right]. \end{gathered}$$

The $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, and $\mathbf{D}$ matrices are directly related to the input and output ray vectors, whose components are denoted by the pre-superscripts [46]:

(B15)$$\left[\begin{array}{cc}\mathbf{A}& \mathbf{B}\\ \mathbf{C}& \mathbf{D}\end{array}\right]=\left[\begin{array}{cc}{}^1{p}_{2,\pm x}{/}^1{p}_{1,\pm x}& {}^1{p}_{2,\pm y}{/}^2{p}_{1,\pm y}\\ {}^2{p}_{2,\pm x}{/}^1{p}_{1,\pm x}& {}^2{p}_{2,\pm y}{/}^2{p}_{1,\pm y}\\ {}^3{p}_{2,\pm x}{/}^1{p}_{1,\pm x}& {}^3{p}_{2,\pm y}{/}^2{p}_{1,\pm y}\\ {}^4{p}_{2,\pm x}{/}^1{p}_{1,\pm x}& {}^4{p}_{2,\pm y}{/}^2{p}_{1,\pm y}\end{array}\begin{array}{cc}{}^1{q}_{2,\pm x}{/}^3{q}_{1,\pm x}& {}^1{q}_{2,\pm y}{/}^4{q}_{1,\pm y}\\ {}^2{q}_{2,\pm x}{/}^3{q}_{1,\pm x}& {}^2{q}_{2,\pm y}{/}^4{q}_{1,\pm y}\\ {}^3{q}_{2,\pm x}{/}^3{q}_{1,\pm x}& {}^3{q}_{2,\pm y}{/}^4{q}_{1,\pm y}\\ {}^4{q}_{2,\pm x}{/}^3{q}_{1,\pm x}& {}^4{q}_{2,\pm y}{/}^4{q}_{1,\pm y}\end{array}\right].$$

Therefore, the parallel initial ray vectors are chosen if $\mathbf{A}$ and $\mathbf{C}$ are sought, while the divergent initial ray vectors lead to $\mathbf{B}$ and $\mathbf{D}$. Fortunately, as is shown below, the VRBDI also needs no more than either $\mathbf{A}$ and $\mathbf{C}$ in the plane wave or $\mathbf{B}$ and $\mathbf{D}$ in the spherical wave decomposition case.

With Eqs. (33) and (34) repeatedly using Eqs. (21) or (22) after $m-1$ surfaces and with Eq. (35), a field contribution on the base plane can now be given either as [16,58]

(B16)$$\begin{gathered} {(\mathrm{d}\mathbf{E},\mathrm{d}\mathbf{H})}_{\mathrm{bp}}=\frac{\exp (-\mathrm{i}{k}_0\mathrm{OPL})}{\sqrt{\det (\mathbf{A})}}\exp \{-\frac{\mathrm{i}{k}_0}{2\det (\mathbf{A})} \\ ·[{\breve{x}}^2(C_{11}{A}_{22}-C_{12}{A}_{21})+{\breve{y}}^2(C_{22}{A}_{11}-C_{21}{A}_{12}) \\ +\breve{x}\breve{y}(C_{12}{A}_{11}-C_{11}{A}_{12}+C_{21}{A}_{22}-C_{22}{A}_{21})]\} \\ ·{(\mathrm{d}\mathbf{E},\mathrm{d}\mathbf{H})}_{m-1} \end{gathered}$$

for the plane or

(B17)$$\begin{gathered} {(\mathrm{d}\mathbf{E},\mathrm{d}\mathbf{H})}_{\mathrm{bp}}=\frac{\mathrm{i}}{{\lambda }_0}\frac{\exp (-\mathrm{i}{k}_0\mathrm{OPL})}{\sqrt{\det (\mathbf{B})}}(1-\frac{\mathrm{i}}{k_0\mathrm{OPL}})\exp \{-\frac{\mathrm{i}{k}_0}{2\det (\mathbf{B})} \\ ·[{\breve{x}}^2(D_{11}{B}_{22}-D_{12}{B}_{21})+{\breve{y}}^2(D_{22}{B}_{11}-D_{21}{B}_{12}) \\ +\breve{x}\breve{y}(D_{12}{B}_{11}-D_{11}{B}_{12}+D_{21}{B}_{22}-D_{22}{B}_{21})]\} \\ ·{(\mathrm{d}\mathbf{E},\mathrm{d}\mathbf{H})}_{m-1} \end{gathered}$$

for the spherical wave decomposition and the optical path OPL obtained by Eq. (36).

However, the field in the base plane ${(\mathrm{d}\mathbf{E},\mathrm{d}\mathbf{H})}_{\mathrm{bp}}$ is not the sought quantity. Additionally, the opportunity to choose only one of four possible $\pm$-combinations for the ray vectors implies an ambiguity. Figure 20 illustrates these problems and offers solutions. First, the nearby sample-grid point ${\mathbf{P}}_{\mathrm{sgp}}$ is expressed in the base ray frame:

(B18)$$\left[\begin{array}{c}{\breve{x}}_{\mathrm{L}}\\ {\breve{y}}_{\mathrm{L}}\\ {\breve{z}}_{\mathrm{L}}\end{array}\right]=\left[\begin{array}{c}({\mathbf{P}}_{\mathrm{sgp}}-\mathbf{P})•{\widehat{\mathbf{e}}}_x\\ ({\mathbf{P}}_{\mathrm{sgp}}-\mathbf{P})•{\widehat{\mathbf{e}}}_y\\ ({\mathbf{P}}_{\mathrm{sgp}}-\mathbf{P})•\widehat{\mathbf{r}}\end{array}\right].$$

 

Fig. 20. Situation after ray tracing. Left: the intersection $\mathbf{P}$ of the base ray with the detector plane defines the orthogonal base plane, which does, in general, not contain the nearby sample-grid point ${\mathbf{P}}_{\mathrm{sgp}}$. Further propagation by $\mathrm{\Delta }z$ along the base ray is necessary. Right: the four possible ±- combinations of ray vectors are decided by the position of ${\mathbf{P}}_{\mathrm{sgp}}$.

Download Full Size | PDF

Depending on the signs of ${\breve{x}}_{\mathrm{L}}$ and ${\breve{y}}_{\mathrm{L}}$, the initial and final ray vectors are chosen as indicated by Fig. 20. Then, by use of Eq. (B15), the corresponding submatrices $\mathbf{A}$ and $\mathbf{C}$ or $\mathbf{B}$ and $\mathbf{D}$ are calculated. Further propagation into base ray direction can simply be described by

(B19)$$\left[\begin{array}{cc}\tilde{\mathbf{A}}& \tilde{\mathbf{B}}\\ \mathbf{C}& \mathbf{D}\end{array}\right]=\left[\begin{array}{cc}\mathbf{I}& \mathrm{\Delta }z/n_{m-1}\mathbf{I}\\ \mathbf{0}& \mathbf{I}\end{array}\right]\left[\begin{array}{cc}\mathbf{A}& \mathbf{B}\\ \mathbf{C}& \mathbf{D}\end{array}\right].$$

Therefore, only

(B20)$$\tilde{\mathbf{A}}=\mathbf{A}+\mathrm{\Delta }z/n_{m-1}\mathbf{C}\text{or}\tilde{\mathbf{B}}=\mathbf{B}+\mathrm{\Delta }z/n_{m-1}\mathbf{D}$$

has to be calculated. It is worth noting that Eq. (B20) holds also if $\mathrm{\Delta }z<0$, which corresponds to a backward propagation. For example, this is the case when ${\mathbf{P}}_{\mathrm{sgp}}$ lies on the other side of $\mathbf{P}$ in Fig. 20 (left).

To properly include the well-known Gouy phase anomaly, which causes a $\pi /2$ shift at each astigmatic focus ($\pi$ at a symmetric double focus) [39,63] and which was also confirmed experimentally [64], it is still necessary to modify $\det {(\tilde{\mathbf{A}})}^{-1/2}$ or $\det {(\tilde{\mathbf{B}})}^{-1/2}$, respectively. It can be shown that

(B21)$$\begin{gathered} \frac{1}{\sqrt{\det (\tilde{\mathbf{A}})}}\to \frac{\sigma ({À}_{11})}{\sqrt{{À}_{11}}}\frac{\sigma ({À}_{22})}{\sqrt{{À}_{22}}}, \\ \frac{1}{\sqrt{\det (\tilde{\mathbf{B}})}}\to \frac{\sigma ({\mathrm{B̀}}_{11})}{\sqrt{{\mathrm{B̀}}_{11}}}\frac{\sigma ({\mathrm{B̀}}_{22})}{\sqrt{{\mathrm{B̀}}_{22}}}, \end{gathered}$$

where $\mathbf{À}$ and $\mathbf{B̀}$ are the diagonalized matrices of $\tilde{\mathbf{A}}$ and $\tilde{\mathbf{B}}$, respectively, describes the Gouy phase correctly. It should be noted that, e.g., in the case of a higher-order Hermite–Gaussian beam, additional order-dependent contributions to the Gouy phase, as they typically appear in respective analytical formulas, are implicitly included by the integration of all input field components at the output plane and need not to be considered explicitly here.

Eventually, a field contribution on a detector sampling grid point can be written as

(B22)$$\begin{gathered} {(\mathrm{d}\mathbf{E},\mathrm{d}\mathbf{H})}_{\mathrm{sgp}}=\frac{\sigma ({À}_{11})}{\sqrt{{À}_{11}}}\frac{\sigma ({À}_{22})}{\sqrt{{À}_{22}}}\exp (-\mathrm{i}{k}_0\mathrm{OPL}) \\ ·\exp \{-\frac{\mathrm{i}{k}_0}{2\det (\tilde{\mathbf{A}})}[{\breve{x}}_{\mathrm{L}}^2(C_{11}{\tilde{A}}_{22}-C_{12}{\tilde{A}}_{21}) \\ +{\breve{y}}_{\mathrm{L}}^2(C_{22}{\tilde{A}}_{11}-C_{21}{\tilde{A}}_{12}) \\ +{\breve{x}}_{\mathrm{L}}{\breve{y}}_{\mathrm{L}}(C_{12}{\tilde{A}}_{11}-C_{11}{\tilde{A}}_{12}+C_{21}{\tilde{A}}_{22}-C_{22}{\tilde{A}}_{21})]\} \\ ·{(\mathrm{d}\mathbf{E},\mathrm{d}\mathbf{H})}_{m-1} \end{gathered}$$

for the plane or

(B23)$$\begin{gathered} {(\mathrm{d}\mathbf{E},\mathrm{d}\mathbf{H})}_{\mathrm{sgp}}=\frac{\mathrm{i}}{{\lambda }_0}\frac{\sigma ({\mathrm{B̀}}_{11})}{\sqrt{{\mathrm{B̀}}_{11}}}\frac{\sigma ({\mathrm{B̀}}_{22})}{\sqrt{{\mathrm{B̀}}_{22}}}\exp (-\mathrm{i}{k}_0\mathrm{OPL}) \\ ·(1-\frac{\mathrm{i}}{k_0\mathrm{OPL}})\exp \{-\frac{\mathrm{i}{k}_0}{2\det (\tilde{\mathbf{B}})} \\ ·[{\mathrm{x̆}}_{\mathrm{L}}^2(D_{11}{\tilde{B}}_{22}-D_{12}{\tilde{B}}_{21})+{\mathrm{y̆}}_{\mathrm{L}}^2(D_{22}{\tilde{B}}_{11}-D_{21}{\tilde{B}}_{12}) \\ +{\mathrm{x̆}}_{\mathrm{L}}{\mathrm{y̆}}_{\mathrm{L}}(D_{12}{\tilde{B}}_{11}-D_{11}{\tilde{B}}_{12}+D_{21}{\tilde{B}}_{22}-D_{22}{\tilde{B}}_{21})]\} \\ ·{(\mathrm{d}\mathbf{E},\mathrm{d}\mathbf{H})}_{m-1} \end{gathered}$$

for the spherical wave decomposition. Therefore, $A{F}_{\mathrm{PW},\mathrm{SW}}$ and $S_{\mathrm{PW},\mathrm{SW}}$ in Eqs. (37) and (38) are

(B24)$$\begin{gathered} A{F}_{\mathrm{PW}}=\frac{\sigma ({À}_{11})}{\sqrt{{À}_{11}}}\frac{\sigma ({À}_{22})}{\sqrt{{À}_{22}}}, \\ {S}_{\mathrm{PW}}=\frac{1}{2\det (\tilde{\mathbf{A}})}[{\breve{x}}_{\mathrm{L}}^2(C_{11}{\tilde{A}}_{22}-C_{12}{\tilde{A}}_{21}) \\ +{\breve{y}}_{\mathrm{L}}^2(C_{22}{\tilde{A}}_{11}-C_{21}{\tilde{A}}_{12}) \\ +{\breve{x}}_{\mathrm{L}}{\breve{y}}_{\mathrm{L}}(C_{12}{\tilde{A}}_{11}-C_{11}{\tilde{A}}_{12}+C_{21}{\tilde{A}}_{22}-C_{22}{\tilde{A}}_{21})] \end{gathered}$$

for the plane or

(B25)$$\begin{gathered} A{F}_{\mathrm{SW}}=\frac{\sigma ({\mathrm{B̀}}_{11})}{\sqrt{{\mathrm{B̀}}_{11}}}\frac{\sigma ({\mathrm{B̀}}_{22})}{\sqrt{{\mathrm{B̀}}_{22}}}, \\ {S}_{\mathrm{SW}}=\frac{1}{2\det (\tilde{\mathbf{B}})}[{\mathrm{x̆}}_{\mathrm{L}}^2(D_{11}{\tilde{B}}_{22}-D_{12}{\tilde{B}}_{21}) \\ +{\mathrm{y̆}}_{\mathrm{L}}^2(D_{22}{\tilde{B}}_{11}-D_{21}{\tilde{B}}_{12}) \\ +{\mathrm{x̆}}_{\mathrm{L}}{\mathrm{y̆}}_{\mathrm{L}}(D_{12}{\tilde{B}}_{11}-D_{11}{\tilde{B}}_{12}+D_{21}{\tilde{B}}_{22}-D_{22}{\tilde{B}}_{21})] \end{gathered}$$

for the spherical wave decomposition. The term $S_{\mathrm{SW}}$ is the paraxial point-to-point and $S_{\mathrm{PW}}$ the paraxial angle-to-point eikonal, respectively [58].

APPENDIX C: MAGNETIC FIELD COMPONENTS

In this section, the magnetic field components corresponding to the main text figures showing the electric ones are collected in Figs. 21 –25 the sake of completeness.

 

Fig. 21. Magnetic field components corresponding to Fig. 3.

Download Full Size | PDF

 

Fig. 22. Magnetic field components corresponding to Fig. 4.

Download Full Size | PDF

 

Fig. 23. Magnetic field components corresponding to Fig. 6.

Download Full Size | PDF

 

Fig. 24. Magnetic field components corresponding to Fig. 10.

Download Full Size | PDF

 

Fig. 25. Magnetic field components corresponding to Fig. 12.

Download Full Size | PDF

FUNDING INFORMATION

EMRP Joint Research Project (SIB08); European Union.

This work is a study for the EMRP Joint Research Project SIB08: Traceability of sub-nm length measurements (subnano). The EMRP is jointly funded by the European Union. The authors gratefully acknowledge funding by the EMRP.

REFERENCES

1. J. N. Dukes and G. B. Gordon, “A two-hundred-foot yardstick with graduations every microinch,” Hewlett-Packard J. 21, 2–8 (1970).

2. K. Dorenwendt and G. Bönsch, “Über den Einfluß der Beugung auf die interferentielle Längenmessung,” Metrologia 12, 57–60 (1976). [CrossRef]  

3. J.-P. Monchalin, M. J. Kelly, J. E. Thomas, N. A. Kurnit, A. Szöke, F. Zernike, P. H. Lee, and A. Javan, “Accurate laser wavelength measurement with a precision two-beam scanning Michelson interferometer,” Appl. Opt. 20, 736–757 (1981). [CrossRef]  

4. G. Bönsch, “Wavelength ratio of stabilized laser radiation at 3.39 μm and 0.633 μm,” Appl. Opt. 22, 3414–3420 (1983). [CrossRef]  

5. G. Mana, “Diffraction effects in optical interferometers illuminated by laser sources,” Metrologia 26, 87–93 (1989). [CrossRef]  

6. A. Bergamin, G. Cavagnero, and G. Mana, “Observation of Fresnel diffraction in a two-beam laser interferometer,” Phys. Rev. A 49, 2167–2173 (1994). [CrossRef]  

7. A. Bergamin, G. Cavagnero, L. Cordiali, and G. Mana, “Beam-astigmatism in laser interferometry,” IEEE Trans. Instr. Meas. 46, 196–200 (1997). [CrossRef]  

8. A. Bergamin, G. Cavagnero, L. Cordiali, and G. Mana, “A Fourier optics model of two-beam scanning laser interferometers,” Eur. Phys. J. D 5, 433–440 (1999). [CrossRef]  

9. R. Masui, “Effect of diffraction in a Saunders-type optical interferometer,” Metrologia 34, 125–131 (1997). [CrossRef]  

10. G. Cavagnero, G. Mana, and E. Massa, “Aberration effects in two-beam laser interferometers,” J. Opt. Soc. Am. A 23, 1951–1959 (2006). [CrossRef]  

11. N. Kuramoto, K. Fujii, and K. Yamazawa, “Volume measurement of ²⁸Si spheres using an interferometer with a flat etalon to determine the Avogadro constant,” Metrologia 48, S83–S95 (2011). [CrossRef]  

12. B. Andreas, L. Ferroglio, K. Fujii, N. Kuramoto, and G. Mana, “Phase corrections in the optical interferometer for Si sphere volume measurements at NMIJ,” Metrologia 48, S104–S111 (2011). [CrossRef]  

13. B. Andreas, K. Fujii, N. Kuramoto, and G. Mana, “The uncertainty of the phase-correction in sphere-diameter measurements,” Metrologia 49, 479–486 (2012). [CrossRef]  

14. E. Wolf, “Electromagnetic diffraction in optical systems I. An integral representation of the image field,” Proc. R. Soc. A 253, 349–357 (1959). [CrossRef]  

15. B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems II. Structure of the image field in an aplanatic system,” Proc. R. Soc. A 253, 358–379 (1959). [CrossRef]  

16. S. A. Collins, “Lens-system diffraction integral written in terms of matrix optics,” J. Opt. Lett. 60, 1168–1177 (1970).

17. H. Kogelnik and T. Li, “Laser beams and resonators,” Appl. Opt. 5, 1550–1567 (1966). [CrossRef]  

18. A. E. Siegman, Lasers (University Science Books, 1986).

19. J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1968).

20. E. Byckling and J. Simola, “Wave optical calculation of diffraction effects in optical systems,” Opt. Acta 27, 337–344 (2010). [CrossRef]  

21. A. S. Glassner, An Introduction to Ray Tracing (Morgan Kaufmann, 1989).

22. N. G. Douglas, A. R. Jones, and F. J. van Hoesel, “Ray-based simulation of an optical interferometer,” J. Opt. Soc. Am. A 12, 124–131 (1995). [CrossRef]  

23. F. Riechert, F. Dürr, U. Rohlfing, and U. Lemmer, “Ray-based simulation of the propagation of light with different degrees of coherence through complex optical systems,” Appl. Opt. 48, 1527–1534 (2009). [CrossRef]  

24. F. Wyrowski and M. Kuhn, “Introduction to field tracing,” J. Mod. Opt. 58, 449–466 (2011). [CrossRef]  

25. A. W. Greynolds, “Propagation of generally astigmatic Gaussian beams along skew ray paths,” Proc. SPIE 560, 33–50 (1985).

26. M. J. Bastiaans, “The expansion of an optical signal into a discrete set of Gaussian beams,” Optik 57, 95–102 (1980).

27. R. Herloski, S. Marshall, and R. Antos, “Gaussian beam ray-equivalent modeling and optical design,” Appl. Opt. 22, 1168–1174 (1983). [CrossRef]  

28. J. Arnaud, “Representation of Gaussian beams by complex rays,” Appl. Opt. 24, 538–543 (1985). [CrossRef]  

29. G. W. Forbes and M. A. Alonso, “The Holy Grail of ray-based optical modeling,” Proc. SPIE 4832, 186–197 (2002).

30. M. A. Alonso, “Ray-based diffraction calculations using stable aggregate of flexible elements,” J. Opt. Soc. Am. A 30, 1223–1235 (2013). [CrossRef]  

31. A. V. Gitin, “Huygens-Feynman-Fresnel principle as the basis of applied optics,” Appl. Opt. 52, 7419–7434 (2013). [CrossRef]  

32. R. P. Feynman, “Space-time approach to non-relativistic quantum mechanics,” Rev. Mod. Phys. 20, 367–387 (1948). [CrossRef]  

33. R. Mittra, ed., Computer Techniques for Electromagnetics (Pergamon, 1973).

34. G. A. Deschamps, “Ray techniques in electromagnetics,” Proc. IEEE 60, 1022–1035 (1972). [CrossRef]  

35. S. W. Lee, P. Cramer Jr., K. Woo, and Y. Rahmat-Samii, “Diffraction by an arbitrary subreflector: GTD solution,” IEEE Trans. Antennas Propagat. AP-27, 305–316 (1979).

36. R. Mittra and A. Rushdi, “An efficient approach for computing the geometrical optics field reflected from a numerically specified surface,” IEEE Trans. Antennas Propagat. 27, 871–877 (1979).

37. S. W. Lee, M. S. Sheshadri, V. Jamnejad, and R. Mittra, “Refraction at a curved dielectric interface: geometrical optics solution,” IEEE Trans. Microwave Theory Trans. MTT-30, 12–19 (1982).

38. S. W. Lee and P. T. Lam, “Correction to diffraction by an arbitrary subreflector: GTD solution,” IEEE Trans. Antennas Propagat. AP-34, 272 (1986).

39. H. Ling, R. C. Chou, and S. W. Lee, “Shooting and bouncing rays: calculating the RCS of an arbitrarily shaped cavity,” IEEE Trans. Antennas Propagat. 37, 194–205 (1989).

40. S. W. Lee, H. Ling, and R. C. Chou, “Ray-tube integration in shooting and bouncing ray method,” Microwave Opt. Technol. Lett. 1, 286–289 (1988). [CrossRef]  

41. J. Baldauf, S. W. Lee, L. Lin, S. K. Jeng, S. M. Scarborough, and C. L. Yu, “High frequency scattering from trihedral corner reflectors and other benchmark targets: SBR versus experiment,” IEEE Trans. Antennas Propagat. 39, 1345–1351 (1991).

42. Y. Tao, H. Lin, and H. Bao, “GPU-based shooting and bouncing ray method for fast RCS prediction,” IEEE Trans. Antennas Propagat. 58, 494–502 (2010).

43. C. Y. Kee and C. F. Wang, “Efficient GPU implementation of the high-frequency SBR-PO method,” IEEE Antennas Wireless Propagat. Lett. 12, 941–944 (2013). [CrossRef]  

44. B. Andreas, G. Mana, E. Massa, and C. Palmisano, “Modeling laser interferometers for the measurement of the Avogadro constant,” Proc. SPIE 8789, 87890W (2013).

45. A. Gerrard and J. M. Burch, Introduction to Matrix Methods in Optics (Wiley, 1975).

46. M. Harrigan, “General beam propagation through non-orthogonal optical systems,” Proc. SPIE 4832, 379–389 (2002).

47. C. G. Chen, P. T. Konkola, J. Ferrera, R. K. Heilmann, and M. L. Schattenburg, “Analyses of vector Gaussian beam propagation and the validity of paraxial and spherical approximations,” J. Opt. Soc. Am. A 19, 404–412 (2002). [CrossRef]  

48. G. S. Smith, An Introduction to Classical Electromagnetic Radiation (Cambridge University, 1997).

49. J. D. Jackson, Classical Electrodynamics (Wiley, 1999).

50. M. Born and E. Wolf, Principles of Optics (Pergamon, 1964).

51. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University, 2008).

52. W. R. Smythe, “The double current sheet in diffraction,” Phys. Rev. 72, 1066–1070 (1947). [CrossRef]  

53. A. Drezet, J. C. Woehl, and S. Huant, “Diffraction of light by a planar aperture in a metallic screen,” J. Math. Phys. 47, 072901 (2006). [CrossRef]  

54. J. A. Stratton and L. J. Chu, “Diffraction theory of electromagnetic waves,” Phys. Rev. 56, 99–107 (1939). [CrossRef]  

55. S. Guha and G. D. Gillen, “Vector diffraction theory of refraction of light by a spherical surface,” J. Opt. Soc. Am. B 24, 1–8 (2007). [CrossRef]  

56. S. R. Ghorpade and B. V. Limaye, A Course in Multivariable Calculus and Analysis (Springer, 2010).

57. M. Hillenbrand, D. P. Kelly, and S. Sinzinger, “Numerical solution of nonparaxial scalar diffraction integrals for focused fields,” J. Opt. Soc. Am. A 31, 1832–1841 (2014). [CrossRef]  

58. A. Walther, The Ray and Wave Theory of Lenses (Cambridge University, 1995).

59. R. P. Feynman, R. B. Leighton, and M. Sands, The Feynman Lectures on Physics (California Institute of Technology, 2010), http://www.feynmanlectures.caltech.edu.

60. K. Y. Bliokh, A. Y. Bekshaev, A. G. Kofman, and F. Nori, “Photon trajectories, anomalous velocities and weak measurements: a classical interpretation,” New J. Phys. 15, 1–17 (2013).

61. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, 3rd ed. (Cambridge University, 2007).

62. http://en.wikipedia.org/wiki/Rodrigues'_rotation_formula.

63. T. D. Visser and E. Wolf, “The origin of the Gouy phase anomaly and its generalization to astigmatic wavefields,” Opt. Commun. 283, 3371–3375 (2010). [CrossRef]  

64. J. P. Rolland, K. P. Thompson, K. S. Lee, and J. Tamkin Jr., T. Schmid and E. Wolf, “Observation of the Gouy phase anomaly in astigmatic beams,” Appl. Opt. 51, 2902–2908 (2012). [CrossRef]