Open Access
Add to BibSonomy
Abstract
Light scattering from residual manufacturing errors of optical surfaces has a large impact on the image quality of optical systems. Classical ray-based methods to simulate surface scattering in optical systems depend on statistical models of surface errors and neglect the wave properties of light, which prohibit the integration of statistical surface error models with beam propagation methods. Additionally, the impact of multiple scattering from different frequency components of surface errors cannot be easily modelled by existing methods. Here we analyze the impact of different frequency components of surface errors induced by diamond-turned surface grinding on image quality, and we propose a Wigner function-based approach in which light is modelled as partially coherent. In this unified model, by selecting the proper definition of light coherence, we can combine the statistical and deterministic models of surface errors, enabling efficient, simultaneous simulation of multiple scattering from high- and mid-spatial frequency (HSF and MSF, respectively) surface errors, as well as the interference and edge diffraction of light.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
In real imaging optical systems, the residual manufacturing errors of optical surfaces can never be completely eliminated, and the estimation of their impact on the image quality is necessary for design and tolerancing of the optical system. Figure 1 shows three examples of typical surface errors. Here we see that the microstructures of different types of surface errors vary. While the HSF errors are highly irregular and unpredictable, the MSF errors show regular structures resulting from the manufacturing process. Therefore, rather than using exact profiles of all surface errors, the residual surface errors are usually characterized by a power spectral density (PSD), which is the Fourier transform of the autocorrelation function of surface variations [1]. Figure 2 shows an example of the PSD of an optical surface manufactured by surface grinding [2]. From this figure, we observe that the surface errors can be divided into three categories according to their spatial frequencies. Each of the frequency components of the PSD has a different impact on optical performance and thus various methods have been developed to model them effectively. For example, Harvey et al. proposed an analytical method to estimate the impact of HSF component of the surface errors (microroughness) on the transfer function of the optical system [1], Peterson proposed methods based on geometrical ray tracing [3], and Liang et al. proposed field propagation methods based on perturbation theory [4]. However, none of these methods can simultaneously handle statistical angular scattering and rigorous deterministic wave diffraction, nor can they deal with the propagation of partially coherent light in real optical systems with scattering surfaces, which limit their applications for optical systems in which the partial coherence of light and multiple scattering are of concern, such as the UV lithography systems in which partially coherent illumination is applied to increase resolution [5]. Furthermore, there lacks a unified model that is capable of modelling multiple scatterings from all frequency components of the PSD simultaneously. To solve these problems, we propose a Wigner function-based approach to propagate partially coherent light in real optical systems with scattering surfaces. The Wigner function has been widely applied to model light scattering in turbid medium [6,7], and a primary advantage of using the Wigner function is its capability to characterize and propagate partially coherent light field, which is usually encountered in reality. Another advantage of utilizing the coherence theory of light in the modelling of surface scattering is that it allows us to apply statistical and deterministic models of surface errors simultaneously. As we will see later, statistical and deterministic surface models are necessary to model the MSF and HSF components of the PSD, which have different influence on the spatial coherence of light, and the Wigner function is able to describe both of them in a single model. Furthermore, the phase space representation of light by the Wigner function allows direct coupling with the bidirectional scattering distribution function (BSDF), which is commonly used to characterize scattering properties of HSF errors.
Fig. 1. Illustration of surface height variations of three types of surface error. The surface variations of HSF errors are modelled as Gaussian processes with different RMS surface heights. The MSF errors is a typical structure resulting from diamond-turning. Only part of the full aperture is shown for clarity. σh and lc are the standard deviation and correlation length of surface height variation, T is the period of the MSF structure.
Fig. 2. An illustration of the log-log scaled PSD of optical surfaces manufactured by surface grinding [2].
2. Surface scattering and the coherence of light
2.1 Impact of different types of surface errors on image quality
Restricting our attention to monochromatic light, let us consider the modulation transfer function (MTF) of the optical system with real manufactured surfaces. Light from a monochromatic point object is fully coherent, and therefore, if the surface errors do not vary with time, the light field at the exit pupil is also fully coherent, even though it contains extremely complex phase variations due to the residual surface errors. Therefore, the amplitude distribution at the image plane can be considered as the coherent superposition of all the Huygens wavelets from the exit pupil.
As shown in Fig. 3, the complex field at the exit pupil is given by the intrinsic pupil function Pi(ξ,η) and the space-variant random phase deviation Ps(ξ,η) induced by the residual errors. The intrinsic part of the pupil function is inherent to the optical design and can be obtained by raytracing or more rigorous beam propagation methods [8], while the impact of the residual surface error on the complex pupil function is more complicated due to the presence of MSF and HSF surface errors.
Fig. 3. Schematic sketch of an imaging optical system with scattering surfaces, Pi(ξ,η) is the intrinsic pupil function that includes all the geometrical aberrations, apodizations and edge diffractions, while Ps(ξ,η) is the phase variation at the exit pupil induced by residual surface errors.
According to the theories of Fourier optics [9], the exact MTF of an optical system can be calculated by the autocorrelation of the complex pupil function, as indicated by the solid lines in Fig. 4. However, the fine structures that we observe in the MTFs of HSF errors correspond to the speckle structures and yield no useful information on the image quality. Alternatively, Harvey et al. have presented a more useful representation of the MTF as the ensemble average over different systems with the same design and statistical surface properties [1]. According to the Harvey model, if the surface error can be modelled by a zero-mean Gaussian process with wide-sense stationarity, the impact of surface error on the MTF is determined by the autocorrelation of surface errors. The dashed lines in Fig. 4 corresponds to the averaged MTFs calculated by the Harvey model. The differences between the exact and average transfer function indicate that the exact MTF induced by microroughness is of little use and the average MTF predicted by the Harvey model is adequate to characterize the image quality degradation. This is modeled as a constant ratio of contrast reduction across all spatial frequencies. However, the same conclusion cannot be drawn for MSF errors such as surface 3, whose exact MTF strongly oscillates around the predicted MTF with large deviations. This leads to the exact MTF falling to 0 at spatial frequencies smaller than the cut-off frequency, indicating resolution loss.
Fig. 4. MTF of the optical systems with residual surface errors. The solid lines correspond to the actual MTFs of the surface errors 1-3 shown in Fig. 1. The dashed lines indicate the MTFs predicted by the Harvey model.
On typical optical surfaces, the MSF and HSF errors coexist. Therefore, it is necessary to consider the impact of both types of surface error simultaneously when the image quality of as-built optical systems are to be assessed. This requires simultaneous representation of the surface errors by the statistical and deterministic models. Therefore, the impact of surface scattering on the coherence of light should be investigated, for which we need to utilize the time and ensemble-average-based definitions of light coherence simultaneously.
2.2 Ensemble and time-average coherence
The time and ensemble-average-based definitions of light coherence are usually applied separately. However, in order to model surface scattering as the change of light coherence, both definitions must be applied simultaneously. Therefore, it is necessary to clarify this combined definition of spatial coherence before we proceed. Here we restrict our discussion to the one-dimensional case for the sake of simplicity and without loss of generality.
The classical definition of spatial coherence is based on the time average [10], with the mutual intensity defined as
The above equation gives the average of the mutual intensity function of an ensemble of optical systems with the same design and surface properties. In this case, the statistical properties of the surface errors can be reflected in the mutual intensity function by the ensemble average, enabling us to interpret surface scattering as the change of coherence of the incident light.
2.3 Coherence of scattered light
According to Appendix A, for a partially coherent incident light whose mutual intensity function is ${\overline{\textbf{J}} _i}({P_1},{P_2})$ that passes through a non-absorbing scattering surface, the mutual intensity function of the beam after the scattering surface is calculated as
Substituting Eq. (4) into Eq. (3), we obtain the mutual intensity of the light after the scattering surface
The phase variations $\varphi ({P_1})$ and $\varphi ({P_2})$ can be considered as two random variables which are proportional to the surface height of the scattering surface at P1 and P2, and the mathematical expectation $E\{{\exp [ - j(\varphi ({P_1}) - \varphi ({P_2}))]} \}$ can be calculated by the joint characteristic function Φ(ω1,ω2) of the random variables [11]
If $\varphi ({P_1})$ and $\varphi ({P_2})$ are zero-mean Gaussian processes and jointly normal, the joint characteristic function can be analytically calculated. Substituting Eq. (6) into Eq. (5), we get
Equations (8)–(9) provides useful insights into the impact of the scattering surface on the spatial coherence of the incident light. We can see that the mutual intensity function after the scattering surface has two mutually incoherent components. One is the specularly transmitted part which has the same degree of spatial coherence as the incident light with attenuated intensity, while the other part is the diffusively transmitted part, whose spatial coherence is determined by the degree of coherence of the incident light together with the statistical property of the scattering surface. The percentage of the energy contained in each component is determined by the TIS, which depends on the roughness of the scattering surface.
From Eqs. (8)–(9) we also learn that for a coherent incident light field, the scattered part of the field after a surface shares the same autocorrelation coefficient as the profile of that surface. For high frequency microroughness with small correlation length, the scattered component can be considered to have very low spatial coherence, meaning that the scattered light is only partially coherent within very small areas of coherence defined by the correlation length of the scattering surface. In this case, the scattered irradiance distribution on a plane after the scattering surface can be obtained by incoherently superposing the light scattered from each area of correlation on the scattering surface. Therefore, Eq. (9) serves as the theoretical foundation for the ray-based method to simulate light scattering in optical systems, which are commonly used despite the fact that their theoretical background is rarely discussed. The ray-based methods assume that the correlation area of the scattered field is much smaller than the scattering surface, such that light scattered by different areas of optical surfaces can be incoherently superposed.
Although the assumption of small correlation area is usually fulfilled for high frequency surface roughness, it fails for MSF errors due to the large correlation length and because the wide-sense stationarity is no longer guaranteed. Therefore, the statistical characterization of the scattering property of surfaces using the BSDF can only be applied to surface roughness. For MSF errors, similar to what we have concluded in the last section based on the MTF, the finite degree of spatial coherence of the scattered field requires us to take the interference and diffraction effect into account. Consequently, rigorous beam propagation methods based on the exact metrology data of the MSF structures need to be applied to evaluate their scattering effect.
Therefore, under the ensemble and time average-based definition of light coherence, surface scattering can be interpreted as the splitting of the incident light field into specular and scattered parts which are mutually incoherent with each other, resulting in a reduction of coherence. This enables us to model the scattered light as partially coherent and to use both statistical and analytical models simultaneously to simulate light scattering due to residual surface errors.
3. Propagation of partially coherent light in optical systems with the Wigner function
Although the modelling of surface scattering by the coherence theory is convenient and provides more physical insights, propagation of the mutual intensity function requires the calculation of the propagation integral, which is rather computationally intensive [10]. Furthermore, light scattered by HSF errors possesses extremely small correlation lengths, which requires an exceptionally fine spatial sampling of the mutual intensity function and is extremely memory-consuming. In contrary, the Fourier transform of the mutual intensity function, which is the Wigner function, are much easier to be propagated by shearing or Radon transform, and most importantly, it can be directly coupled with the BSDF of scattering surfaces in the angular domain. Therefore, compared to the propagation of the mutual intensity function, it is much more convenient to propagate partially coherent light through scattering surfaces in the Fourier domain by the Wigner function. In this section, by applying the time and ensemble average-based definition of light coherence, we demonstrate the method to simulate surface scattering in optical systems by the propagation of partially coherent light based on the Wigner function. For simplicity, we restrict our discussion to 2D Wigner functions.
3.1 Wigner function
The Wigner function of a light field is defined as the Fourier transform of the cross-spectral density function with respect to the displacement Δx [12]
where Γ(x,Δx) is the cross-spectral density function3.2 Propagation of the Wigner function through a scattering surface
Consider an incident light field gi with cross-spectral density function Γi propagating through a thin scattering surface. Applying the TEA, we can write the cross-spectral density function of the scattered light field as
Substituting Eq. (12) into Eq. (10) and applying the convolution theorem of the Fourier transform, the Wigner function of the scattered field can be calculated as
where Wi(x,u) and Wt(x,u) are the Wigner functions of the incident field and the scattering surface, and ${\otimes}$ denotes convolution in the angular domain. Therefore, the Wigner function of the light field after a scattering surface can be calculated as the convolution of the incident Wigner function with the Wigner function of the surface transmission function. Next, we will apply the above derivation and demonstrate the modelling of surface scattering from different types of surface errors.3.3 Scattering from high spatial frequency surface roughness
According to the discussion in Sec. 2, the surface height variations of HSF errors can be modelled by Gaussian processes with wide sense stationarity. Based on this assumption and using the paraxial approximation, we substitute Eq. (13) into Eqs. (10)–(12), then follow the derivations from Eqs. (5)–(8) to obtain the Wigner function of the scattering surface
Here δ(u) is the delta function, and S(u) is the angle spread function which is linearly dependent on the PSD of the scattering surface according to the Harvey-Shack scattering model [1]
where n is the refractive index of the thin element, and λ is the wavelength of the incident light.Subsequently, by substituting Eq. (15) into Eq. (14), we can calculate the Wigner function of the light field after the scattering surface
Therefore, similar to the mutual intensity function, the Wigner function after the scattering surface is composed of two parts. The first is the specular transmission which can be calculated by scaling the incident Wigner function by (1-TIS), while the second is the scattered portion which is calculated by convolving the incident Wigner function with the angle spread function S(u). It should be noted that here S(u) is independent of the incident angles because the derivation of Eq. (17) is based on the paraxial approximation and the TEA, while in reality, if the paraxial approximation is violated, the transmission function t(x) is also dependent on the incident angle, so is the angle spread function. Therefore, the angle spread function S(u) in Eq. (17) is only a slice of the BSDF for normal incidence, while for general cases the Wigner function of the scattered light should be calculated by an integration of the BSDF over the incident angles
In case where the paraxial approximation is fulfilled or if the BSDF is shift-invariant for different incident angles, Eq. (18) reduces to Eq. (17). Besides surface scattering, the free space propagation of the Wigner function is calculated by paraxial ABCD transitions. According to the linear matrix theory, the free space transfer from plane 1 to plane 2 is described as follows:
The above equation indicates that the free space propagation, as well as any first order optical transformation of the Wigner function, are simply shears that can be efficiently calculated with minor computational effort in the Fourier domain [13,14]. This simple method of free space propagation contributes to the high efficiency of the Wigner function method. During the propagation, the intensity profile at an arbitrary surface can be recovered by the marginal integration of the Wigner function:
As an example, Fig. 5 shows the propagation of a top-hat Gaussian Schell beam through a non-absorbing scattering surface with HSF errors.
Fig. 5. Propagation of a Gaussian Schell beam through a scattering surface at z = 3 mm with HSF surface errors. The first three figures show the Wigner function of the beam before the scattering surface (a), after the scattering surface (b), and after 30 mm of free space propagation (c). (d) Shows the irradiance distribution of the beam.
From Fig. 5, we can see that the surface scatters part of the incident beam to large angles, and consequently the spatial extent of the beam increases with propagation in z. During the simulation, the integration in Eq. (18) are evaluated at every x-location, which can be done efficiently by employing vectorization of the calculations. The BSDF of the scattering surface is modelled by the generalized Harvey-Shack (GHS) theory [1] and the PSD of the surface is modelled by the K-correlation method [15]. The wavelength of the monochromatic beam is 1 µm, the correlation length of the surface error is 50 µm, the refractive index of the scattering surface is 1.55, and the standard deviation of the surface error is 0.2 µm.
3.4 Scattering from mid-spatial frequency surface errors
MSF errors are usually related to the specific manufacturing method and possess intrinsic patterns. Therefore, statistical methods are no longer ideal for modelling such errors, and the scattering effect should be modelled using the surface transmission function t(x) based on exact metrology data. In this case, the cross-spectral density function of MSF errors is given as
As an example, we consider a surface with concentric grooves, which is a typical structure resulting from diamond turning as described by surface 3 in Fig. 1. In the 1D case, along the radial direction, the concentric grooves can be modelled by a piecewise-parabolic function, but in reality, the MSF surface structure is not perfectly periodic due to the vibration of the tool tip [16]. The vibrations in the thrust and feed directions induce displacement of the parabolic function in the vertical and horizontal directions. Here we assume that these random displacements follow Gaussian distributions, based on which we have generated an example of the MSF structures as shown in Fig. 6. The period and the PV value of the MSF structure in Fig. 6 are T = 40 µm and hpv = 0.5 µm, and the standard deviations of the tool tip displacement in the thrust and feed directions are σt = 0.05 µm and σf = 4 µm.
Fig. 6. An example of surface structure of MSF surface error. Displacement of the tool tip due to vibrations in thrust and feed directions follows Gaussian distributions, and phase delay Δφ is calculated for the wavelength of 1 µm.
Based on the MSF structure stated above, we have simulated the paraxial propagation of a Gaussian Schell beam though the MSF structure and a single lens based on Eqs. (12)–(13), and the results are shown in Fig. 7 and Fig. 8. From Fig. 7(a) we can see that the MSF structures transform the top-hat beam into many beamlets. In the focal plane of the single lens, we observe strong side lobes on both sides of the central peak. Figure 7(b) shows the PSFs of the system for four MSF structures with different periods and PV values. By comparing the four PSFs, we can see that the PV values of the MSF structures govern the intensity of the side lobes, with larger PV values resulting in weaker central peak and stronger side lobes. On the other hand, the scattering angles are dependent on the periods of the MSF structures, with larger periods corresponding to smaller scattering angles.
Fig. 7. Propagation of a Gaussian Schell beam through a scattering surface with MSF structures and a single lens. (a) shows the irradiance distribution on different z-planes for the scattering surface described in Fig. 6, the single lens having a focal length of f = 40 mm. (b) Shows the irradiance distribution at the beam focus for four MSF structures with different periods and PV values.
Fig. 8. Wigner functions of the Gaussian Schell beam propagated through the optical system described in Fig. 7(a).
The blue curve in Fig. 7(b) corresponds to a critical frequency at which the valleys between the side lobes and central peak disappear and the central peak is broadened, resulting in resolution reduction of the optical system. Therefore, the critical frequency can be used to estimate the impact of MSF structures on image quality. For MSF structures below this critical frequency, the major impact is resolution reduction, while for MSF structures above the critical frequency, the resolution is not significantly reduced, but straylight structures such as ghost images or halos around bright objects will be observed due to the strong side lobes.
Figure 8 shows the Wigner functions of the beam through propagation. In Fig. 8(b)–(d) we observe negative values of the Wigner function, which is a consequence of phase space destructive interference [17]. Finally, the intensity distribution at the focal plane of the single lens results from the interference of all the beamlets formed by the piecewise-parabolic structures.
In real optical surfaces, the HSF and MSF errors coexist. Therefore, in order to calculate the Wigner function of a beam which is scattered by a real optical surface, the Wigner function of the incident beam must be convolved with the Wigner functions of both the HSF and MSF errors in the angular domain:
where WHSF(x,u) and WMSF(x,u) are the Wigner functions of the HSF and MSF errors.4. Application in real optical systems
In the last two sections we have shown the scattering of a partially coherent beam by a single scattering surface. However, real optical systems usually contain more than one surface, and therefore, multiple scattering of partially coherent beams must be considered. In classical straylight analysis tools, multiple scattering is usually modelled by ray-splitting, which is extremely time and memory-consuming because the number of possible ray paths grows exponentially with the number of scattering surfaces. In contrast, the Wigner function approach is more efficient because it models the surface scattering as angular spreading of the incident phase space distribution. Since the propagation of the scattered field is done by a shearing of the Wigner function, no ray-splitting is involved, and the calculation time is linearly dependent on the order of scattering. Furthermore, although the field scattered by surface roughness is incoherent, it becomes partially coherent after it has propagated any distance from the scattering surface according to the van Cittert-Zernike theorem [10]. Therefore, the diffraction of the scattered field from the subsequent edges and MSF structures needs to be considered. Since the Wigner function is able to model partially coherent light, the edge diffraction of the scattered and specular light can be simultaneously modelled with surface scattering. In this section, we demonstrate the Wigner function-based simulation of multiple scattering from HSF and MSF surface errors in real optical systems.
4.1 EUV Schwarzschild objective
As shown in Fig. 9, we consider a EUV Schwarzschild objective [18] composed of two mirrors on which HSF and MSF errors coexist. Considering the short wavelength of the EUV light at 13.5 nm, a large fraction of the incident light is scattered by the two mirrors. Therefore, multiple scattering must be considered to calculate the PSF of the system.
Fig. 9. A EUV Schwarzschild objective with an object space numerical aperture (NA) of 0.19 and an image space NA of 0.004. Only part of the ray bundle in the image space is shown for better visualization. The diameters of M1 and M2 are 50 mm and 10.6 mm, respectively.
In this example, we assume that both mirror surfaces of the Schwarzschild objective possess the same MSF structure with a period of T = 200 µm and a PV value of hpv = 2.5 nm, and the standard deviations of the tool tip displacement in the thrust and feed directions are σt = 0.3 nm and σf = 5 µm. Additionally, the microroughness properties of the two lens surfaces are identical, with a correlation length of lc = 27 nm and a RMS surface height of σs = 0.5 nm. The partially coherent EUV incident beam is modelled by a top-hat Gaussian Schell beam in addition to a parabolic wavefront.
As we can see in Fig. 9, the two mirrors are much larger than the scattering surfaces discussed in Secs. 3.3 and 3.4, and the wavelength of the EUV light is much shorter than the NIR light discussed above. Consequently, the maximum angle that can be described by the Wigner function is strongly constrained by the bandwidth limit of the discrete Fourier transform (DFT). For example, the maximum angle that can be described by a Wigner function with 14,000 × 14,000 sampling points is only 0.4 degree for this EUV Schwarzschild objective. The small angular range of the Wigner function is not a problem for the description of surface scattering since the MSF surface errors produce small-angle scattering within the angular range of the Wigner function, and the light scattered by microroughness to large scattering angles does not reach the detector. However, the small angular range of the Wigner function prohibits us from describing high-NA beams in the optical system since the marginal ray angles are much larger than the angular limit of the Wigner function. In order to solve this problem, we apply a special propagation method to remove the parabolic wavefront from the high-NA beams during propagation, which allows us to convert the convergent and divergent beams into quasi-collimated beams [13,19]. In this case, the angular range of the Wigner function only needs to cover the angular spreading induced by surface scattering and the partial coherence of the beam. In this simulation, the aberrations of the system are modelled by a phase plate at the exit pupil, which contains the residual wavefront aberration of the system at the exit pupil. The irradiance distribution of the beam during propagation is shown in Fig. 10, in which the mirrors have been unfolded to separate the beams for clear demonstration.
Fig. 10. (a) Irradiance distribution of the beam from the entrance window to M1 and M2. (b) Irradiance distribution from M2 to the image plane. (c) Irradiance distribution on the image plane.
From Fig. 10(a) we observe the beamlets formed by the MSF structures, as well as a bright background after M1 induced by microroughness scattering. The irradiance in the white area of Fig. 10(a) has no influence on the scattered field on the image plane and is not calculated in order to reduce memory consumption. The scattered field after M1 is propagated to M2, by which it is reflected and scattered again. The multiply scattered field is then propagated back and truncated by the central hole of M1. As a final step, the truncated field is propagated to the image plane, as shown in Fig. 10(b). The final PSF is shown in Fig. 10(c), from which we observe several side lobes around the central peak due to the MSF errors and a small non-zero offset due to the microroughness. By simulating the single scattering from M1 and M2 individually, we can determine that the side lobes near the central peak are generated by the MSF error on M2, while the side lobes further from the central peak originate from the MSF structure of M1. Since the periods of the MSF structures on M1 and M2 are the same, the angular extension of the scattered light from the two mirrors should also be identical. The different locations of the side lobes originate from the demagnification between M1 and M2, which leads to smaller periods of the phase modulations from M1 at the exit pupil, and thereby larger lobe separation.
The Wigner distributions of the quasi-collimated beam (with parabolic wavefronts removed) after the entrance window, M1, and M2 are shown in Fig. 11. From Fig. 11(a), we can observe the ripples generated by the edge diffraction at the entrance window, while in Fig. 11(b) we observe the complicated structure induced by surface scattering and edge diffraction at M1. Additionally, we observe in Fig. 11(c) that the structure after M2 becomes even more complicated due to multiple scattering. This demonstrates the ability of the Wigner function-based simulation to simultaneously model multiple scattering from MSF errors and surface roughness, as well as diffractions of the specular and scattered beams from the edges of the mirrors and stop.
Fig. 11. Wigner distributions of the beam in the Schwarzschild objective after the entrance window (a), M1 (b), and M2 (c). The parabolic wavefronts of the beams have been removed.
In order to accelerate the simulation and reduce memory consumption, we have applied various methods including parallelization and parabolic wavefront removal. Additionally, we have applied the Radon transform instead of shearing to model the free space propagation after M2, which further accelerates the calculation and eliminates the aliasing effect due to the spreading of the scattered field outside of the x-limit [13].
4.2 Runtime of the Wigner function-based method
One of the advantages of the Wigner function in straylight analysis is its high efficiency due to the absence of Monte Carlo raytracing and ray splitting. Table 1 summarizes the runtime of the simulations presented above. From the table we see that the runtime of the simulation is positively correlated with the number of elements and the resolution of the Wigner function.
Table 1. Runtime of the simulations shown in Fig. 5, Fig. 7 and Fig. 10. The runtime is defined as the time to calculate the intensity distribution on the image plane at the right side of the z-axis without unnecessary z-steps. The simulations are based on a computer platform with an 8-core Intel Xeon E5-2690 CPU @ 2.90 GHz.
4.3 Limitations of the Wigner function-based method
Although the Wigner function approach shows great advantage for modelling multiple scattering in optical systems, it also has clear limitations resulting from the paraxial approximation and TEA, which prevent us from simulating high-NA optical systems with the Wigner function.
According to the TEA, Eq. (13) is only valid for one incidence and scattering angle, and the phase modulation of a scattering surface on the incident field is dependent on the incident and scattering angles. Therefore, when the angular spreading of the beam is greatly increased by a scattering surface, applying the TEA to model the MSF errors of subsequent scattering surfaces under the assumption of normal incidence will induce errors. One solution to this problem is to decompose the scattered field into mutually incoherent modes and to propagate each mode through the scattering surface. However, mode expansion greatly increases the computational effort as the number of scattering surface increases. On the other hand, this problem is less critical for aperture-dominant systems in which only the scattered light inside a small angular range around the specular ray can reach the detector, allowing large scattering angles to be ignored.
Additionally, it should be noted that the Wigner function-based beam propagation method is completely paraxial and based on the ABCD matrix, including the parabolic wavefront removal method which we have used to propagate the high-NA beams in the last section. Therefore, errors are unavoidable during the propagation of high-NA beams. However, there also exist non-paraxial formulations of the Wigner function [20], as well as the corresponding propagation methods [21], which could be utilized to overcome the paraxial limitation.
Furthermore, although the aberrations induced by the optical surfaces can be modelled by a series of differential operators acting on the Wigner function or by a phase plate that contains all the wavefront aberrations [14,22,23], doing so requires the neglection of the surface sag and approximating the optical surfaces as thin sheets, which induces errors for optical surfaces whose surface sag cannot be neglected. One possible solution to overcome the limitation on the surface sag is to generalize the TEA by the local spherical interface approximation (LSIA) [24].
So far, we have limited our simulations in 2D, while simulation in 3D is necessary for optical systems without circular symmetry. In case of 3D beam propagation, the Wigner function becomes a 4D function and thus greatly increases the memory consumption and computational complexity. Therefore, the resolution of the Wigner function is limited in the 3D simulations, making it difficult to model large surfaces with many periods in the MSF structures.
5. Conclusion
We have proposed a Wigner function-based method to simulate the multiple scattering of partially coherent light in optical systems. The application of a time and ensemble average-based definition of light coherence enables us to simultaneously apply statistical and deterministic models of surface errors to model scattering from high and mid-spatial frequency surface errors. Consequently, the Wigner function-based method is able to simulate the edge diffraction and multiple scattering of partially coherent light from MSF and HSF surface errors simultaneously in a unified model. This method can be found useful for beam propagation and stray light analysis of optical systems in which partial coherence and multiple surface scattering is critical for system performance, such as x-ray imaging systems, synchrotron EUV mirror systems, and space telescopes. However, the paraxial approximation and TEA strictly limits the application of the method to paraxial systems, while the large memory consumption of the 4D Wigner function makes it difficult to extend the method to 3D systems with large apertures. Therefore, methods such as sparse 4D matrices should be applied to reduce the memory consumption and enable the application of the Wigner function method in 3D systems with realistic sizes.
Appendix A: derivation of mutual intensity function after a scattering surface
Consider a fully coherent beam incident on a scattering surface, its field distribution before the scattering surface is
while the mutual intensity function according to Eq. (2) isAfter the scattering surface, the field distribution of the beam is given by
and its mutual intensity function isEquation (27) is derived for a fully coherent beam, but considering the fact that any partially coherent light can be decomposed into a set of fully coherent but mutually incoherent modes, Eq. (27) is also valid for partially coherent light [25].
Disclosures
The authors declare no conflicts of interest.
References
1. J. E Harvey, Understanding Surface Scatter Phenomena: a linear systems formulation, (SPIE Press, 2019).
2. B. Dörband, H. Müller, and H. Gross, Handbook of Optical Systems, Volume 5: Metrology of Optical Components and Systems, (John Wiley & Sons, 2012).
3. G. L. Peterson, “Analytic expression for in-field scattered light distribution,” Proc. SPIE 5178 (2004).
4. K. Liang, G. W. Forbes, and M. A. Alonso, “Validity of the perturbation model for the propagation of MSF structure in 2D,” Opt. Express 27(3), 3390–3408 (2019). [CrossRef]
5. A. K. Wong, Resolution enhancement techniques in optical lithography, Vol. 47, (SPIE Press, 2001).
6. A. Wax and J. Thomas, “Measurement of smoothed Wigner phase-space distributions for small-angle scattering in a turbid medium,” JOSA A 15(7), 1896–19081998. [CrossRef]
7. H. Yura, L. Thrane, and P. Andersen, “Closed-form solution for the Wigner phase-space distribution function for diffuse reflection and small-angle scattering in a random medium,” J. Opt. Soc. Am. A 17(12), 2464–2474 (2000). [CrossRef]
8. S. Schmidt, T. Tiess, S. Schröter, R. Hambach, M. Jäger, H. Bartelt, A. Tünnermann, and H. Gross, “Wave-optical modeling beyond the thin-element-approximation,” Opt. Express 24(26), 30188–30200 (2016). [CrossRef]
9. J. W. Goodman, Introduction to Fourier Optics, (Roberts and Company Publishers, 2005).
10. J. W. Goodman, Statistical optics, (John Wiley & Sons, 2015).
11. A. Papoulis and S. U. Pillai, Probability, random variables, and stochastic processes, (Tata McGraw-Hill Education, 2002).
12. M. Testorf, B. Hennelly, and J. Ojeda-Castañeda, Phase-Space Optics: Fundamentals and Applications: Fundamentals and Applications, (McGraw Hill Professional, 2009).
13. M. Zhong and H. Gross, “Propagation algorithms for Wigner functions,” J. Eur. Opt. Soc.-Rapid Publ. 12(1), 7 (2016). [CrossRef]
14. H. Gross, “Numerical propagation of partially coherent laser beams through optical systems,” Opt. Laser Technol. 29(5), 257–260 (1997). [CrossRef]
15. M. G. Dittman, “K-correlation power spectral density and surface scatter model,” In Optical Systems Degradation, Contamination, and Stray Light: Effects, Measurements, and Control II, Vol. 6291, p. 62910R, International Society for Optics and Photonics, 2006.
16. K. Liang and M. A. Alonso, “Effects on the OTF of MSF structures with random variations,” Opt. Express 27(24), 34665–34680 (2019). [CrossRef]
17. D. Dragoman, “Phase-space interferences as the source of negative values of the Wigner distribution function,” J. Opt. Soc. Am. A 17(12), 2481–2485 (2000). [CrossRef]
18. U. Zastrau, C. Rödel, M. Nakatsutsumi, T. Feigl, K. Appel, B. Chen, T. Döppner, T. Fennel, T. Fiedler, L. B. Fletcher, E. Förster, E. Gamboa, D. O. Gericke, S. Göde, C. Grote-Fortmann, V. Hilbert, L. Kazak, T. Laarmann, H. J. Lee, P. Mabey, F. Martinez, K. H. Meiwes-Broer, H. Pauer, M. Perske, A. Przystawik, S. Roling, S. Skruszewicz, M. Shihab, J. Tiggesbäumker, S. Toleikis, M. Wünsche, H. Zacharias, S. H. Glenzer, and G. Gregori, “A sensitive EUV Schwarzschild microscope for plasma studies with sub-micrometer resolution,” Rev. Sci. Instrum. 89(2), 023703 (2018). [CrossRef]
19. A. Siegman, Lasers, (University Science Books, 1986).
20. S. Cho, J. Petruccelli, and M. Alonso, “Wigner functions for paraxial and nonparaxial fields,” J. Mod. Opt. 56(17), 1843–1852 (2009). [CrossRef]
21. Y. Zhang and B. Lü, “Propagation of the Wigner distribution function for partially coherent nonparaxial beams,” Opt. Lett. 29(23), 2710–2712 (2004). [CrossRef]
22. A. W. Lohnmann, J. Ojeda-Castañeda, and N. Streibl, “The influence of wave aberrations on the Wigner distribution,” Opt. Appl. 13(4), (1983).
23. D. Dragoman, “Wigner distribution function applied to third-order aberrations,” Appl. Opt. 35(1), 161–168 (1996). [CrossRef]
24. H. Lajunen, J. Tervo, J. Turunen, T. Vallius, and F. Wyrowski, “Simulation of light propagation by local spherical interface approximation,” Appl. Opt. 42(34), 6804–6810 (2003). [CrossRef]
25. E. Wolf, “New theory of partial coherence in the space–frequency domain. Part I: spectra and cross spectra of steady-state sources,” J. Opt. Soc. Am. 72(3), 343–351 (1982). [CrossRef]
