Numerical visualization of aperiodic scalar optical wave fields

In numerical wave optic simulations using the conventional angular spectrum method, numerical inter-periodic interference occurs as a result of the inherent periodicity of the discrete Fourier transform. In this paper, we propose an aperiodic angular spectrum representation of a scalar optical wave field without numerical inter-periodic interference in an aperiodic area. The visualization of holographic three-dimensional image light field is presented for the validity of the proposed method. It is believed that the proposed method can be broadly applied to visualization of scalar optical wave field represented by the angular spectrum plane wave bases. © 2017 Optical Society of America OCIS codes: (050.1755) Computational electromagnetic methods; (070.0070) Fourier optics and signal processing; (050.1960) Diffraction theory. References and links 1. J. Goodman, Introduction to Fourier Optics, 3rd ed. (Reberts and Company Publishers, 2004). 2. K. Matsushima and S. Nakahara, “Extremely high-definition full-parallax computer-generated hologram created by the polygon-based method,” Appl. Opt. 48(34), H54–H63 (2009). 3. H. Kim, J. Hahn, and B. Lee, “Mathematical modeling of triangle-mesh-modeled three-dimensional surface objects for digital holography,” Appl. Opt. 47(19), D117–D127 (2008). 4. J. Hahn, S. Lim, K. Choi, R. Horisaki, and D. J. Brady, “Video-rate compressive holographic microscopic tomography,” Opt. Express 19(8), 7289–7298 (2011). 5. Y. Lim, J. Hahn, S. Kim, J. Park, H. Kim, and B. Lee, “Plasmonic light beaming manipulation and its detection using holographic microscopy,” IEEE J. Quantum Electron. 46(3), 300–305 (2010). 6. L. Yu and M. K. Kim, “Wavelength-scanning digital interference holography for tomographic three-dimensional imaging by use of the angular spectrum method,” Opt. Lett. 30(16), 2092–2094 (2005). 7. J. Hahn, H. Kim, S.-W. Cho, and B. Lee, “Phase-shifting interferometry with genetic algorithm-based twin image noise elimination,” Appl. Opt. 47(22), 4068–4076 (2008). 8. E. Ye, A. H. Atabaki, N. Han, and R. J. Ram, “Miniature, sub-nanometer resolution Talbot spectrometer,” Opt. Lett. 41(11), 2434–2437 (2016). 9. G. Gibson, J. Courtial, M. Padgett, M. Vasnetsov, V. Pas’ko, S. Barnett, and S. Franke-Arnold, “Free-space information transfer using light beams carrying orbital angular momentum,” Opt. Express 12(22), 5448–5456 (2004). 10. T. Shimobaba, K. Matsushima, T. Kakue, N. Masuda, and T. Ito, “Scaled angular spectrum method,” Opt. Lett. 37(19), 4128–4130 (2012). 11. C.-Y. Hwang, S. Oh, I.-K. Jeong, and H. Kim, “Stepwise angular spectrum method for curved surface diffraction,” Opt. Express 22(10), 12659–12667 (2014). 12. A. Sommerfeld, “Mathematische theorie der diffraction,” Math. Ann. 47(2-3), 317–374 (1896). 13. J. P. Berenger, “Perfect matched layer for the FDTD solution of wave-structure interaction problems,” IEEE Trans. Antenn. Propag. 44(1), 110–117 (1996). 14. J.-P. Hugonin and P. Lalanne, “Perfectly matched layers as nonlinear coordinate transforms: a generalized formalization,” J. Opt. Soc. Am. A 22(9), 1844–1849 (2005). 15. Y. Umul, “Modified theory of physical optics,” Opt. Express 12(20), 4959–4972 (2004). 16. Y. Z. Umil, “Three dimensional modified theory of physical optics,” Optik-Int. J. Light Elect. Opt. 127(2), 819– 824 (2016). 17. Y. Umul, “Diffraction by a black half plane: Modified theory of physical optics approach,” Opt. Express 13(19), 7276–7287 (2005). Vol. 25, No. 13 | 26 Jun 2017 | OPTICS EXPRESS 14715 #292774 https://doi.org/10.1364/OE.25.014715 Journal © 2017 Received 24 Apr 2017; revised 3 Jun 2017; accepted 14 Jun 2017; published 19 Jun 2017


Introduction
The numerical simulation of scalar optical waves is a fundamental problem in optical engineering [1].The numerical model of scalar wave field in free space is crucial as the main theoretical tool for various optical engineering topics such as holographic three-dimensional (3D) display [2,3], holographic microscopy [4][5][6][7], wave-optic based dispersive spectrometer [8], and free space communication [9].Although the fundamental analytic theory for scalar optical field in free space was already established, regarding numerical modeling 3D visualization of scalar optical wave fields, we have much room to research in numerical domain further because of practical limit in computational resource.
In this paper, we address a theoretical issues associated to the numerical modeling of scalar wave field in free space.The angular spectrum representation is the most popular formula for optical wave propagation in free space and known as a perfect analytic solution of arbitrary free space optical scalar wave field.In practice, the numerical model derived from the angular spectrum representation is calculated and visualized in the discrete spatial sampling grids [10,11].Theoretically, the angular spectrum representation is the linear and coherent superposition of infinitely extended plane wave basis with complex amplitude.In continuum, the analytic angular spectrum representation is a perfect solution of scalar wave equation in free space, however, when the spatial domain is discretized and finite-sized, the basic assumption of using infinitely extended plane wave basis and its validity becomes questionable.In practice, the inter-periodic interference that appears in the direct numerical simulation of the conventional angular spectrum representation of aperiodic scalar optical wave fields is problematic.In this paper, its cause is analyzed and a candidate solution for this issue is proposed.The proposed method is based on the famous Sommerfeld analysis of halfinfinite diffraction [12].An alternative form of the angular spectrum representation using the sommerfeld half-infinite diffraction wave function as the wave function basis is developed, which seems to be quite different from the conventional angular spectrum representation using infinite plane wave basis.For the validation of the proposed formula, we present a comparative study of numerical reconstructions of complex digital hologram using the proposed method and the conventional methods.
The paper is organized as follows.In section 2, the inter-periodic interference in conventional angular spectrum representation is addressed and its cause is analyzed.In section 3, a numerical method for numerical visualization of aperiodic scalar optical wave field is proposed and the numerical simulation results of digital hologram reconstruction are presented.The concluding remarks are followed in section 5.

Inter-periodic interference in conventional angular spectrum representation
The scalar wave equation in free space is where 0 k is the wavenumber 0 2 / k π λ

=
. Numerical modeling methods for this theory have been the subject of a great deal of research.In Fourier optics [1] and scalar optical wave propagation has been theoretically established.Various mathematical formulations of ( ) , , U x y z such as the Rayleigh-Sommerfeld formulation and the angular spectrum representation, are well-known.In particular, the angular spectrum representation is a rigorous and versatile formula for the numerical calculation of ( ) , , U x y z and various applications based on scalar wave optic theory.The optical field ( ) , , U x y z is represented by the following analytic formula in 3D space, A k k is the angular spectrum of ( ) , U x y .The form of the total field is interpreted to be the coherent superposition of many plane waves with their own direction in 3D space and complex amplitude ( ) A k k .The total field ( ) , , U x y z can be decomposed into propagating and evanescent components, which are partial integrals in the range of x y z is tightly bound near 0 z = and does not extend to the free space.In numerical computation, the angular spectrum representation is modeled by ( ) , , , , 4 where x Δ and y Δ are the sampling intervals in the x-y spatial domain.x-y spatial frequency domain, which is usually set to an odd number to preserve the origin point of this domain.In practice, Eq. ( 5) can be efficiently calculated using the fast Fourier transform (FFT).For this, , where x T and y T are the x-and y-direction widths of the computation region, respectively.The discrete angular spectrum Eq. ( 5) leads to a periodic pattern in the spatial domain due to the periodicity of the Fourier series.This periodicity produces numerical inter-periodic interference patterns in the resulting field visualization.For example, let us consider the optical field for a single slit aperture under normal incidence, ( ) U x , which is given by ( ) ( ) x W = (6) Figure 1 presents the numerical simulation result for the aperture diffraction using the angular spectrum representation of Eq. ( 5).The uniform plane wave basis in the computation region and the diffraction field contaminated by inter-periodic interference are shown, in Figs.1(a) and 1(b), respectively.The plane wave basis completely fills the rectangular computation region and this is assumed to induce the periodicity shown in Fig. 1(c).This means that the coherent superposition of the infinitely extended plane wave bases in Eq. ( 5) produces the total optical field distribution shown in Fig. 1(d).Infinitely extended plane wave bases influence the entire computational space, preventing the numerical isolation of computational sections.This strong interference pattern is a problem when we want to simulate real single slit diffraction.

Aperiodic angular spectrum representation without inter-periodic interference
In order to address the inter-periodic interference issue, many numerical electromagnetic simulation techniques such as the finite different time domain (FDTD) method and the Fourier modal method (FMM) employ a perfectly matched layer (PML).The PML is a necessary tool for numerical isolation in aperiodic structure analysis [13,14].To our knowledge, an effective PML for the scalar wave equation has yet to be proposed.This means that a scalar angular spectrum method for aperiodic wave field representation has not yet been developed for Fourier optics.The half-infinite barrier can play the same role of the PML for aperiodic visualization of scalar optical wave fields.To present aperiodic plane wave basis instead of the infinite plane wave basis, we employ the approach of the modified theory of physical optics (MTPO) [15,16].The MTPO wave function is the modified representation of the Sommerfeld half-infinite diffraction field.The function of the barrier can be understood in terms of a zero-thickness analytic PML effectively rejecting in-coming wave fields through the vertical side-wall.For a vertical half-infinite barrier ( 0 z ≤ ≤ ∞ , / 2 ), one of the angular spectrum wave bases for the longitudinal angle 0 θ , which corresponds to  With this wave basis generation algorithm, we can complete the set construction of the 3D angular spectrum wave basis and can span the general arbitrary scalar wave field in an isolated cylindrical volume of free space.An alternative to the wave basis is the Sommerfeld solution for the black half-plane problem with a perfect absorber [12].Interestingly, the Sommerfeld solution with a black half plane is an exact equivalent to the incident scattered field in the MTPO solution for a perfect conducting half plane [15][16][17].The only difference is the weighting of the phase ' 0 sinθ .For ' 0 / 2 θ π = , these formulas are the same.From Ref [16], the analytic formula for the incident scattered field of an MTPO wave function with a half-infinite barrier ( 0 z ≤ ≤ ∞ , / 2 ) for ( ) where, P is the observation point ( ) , , x y z , and the function [ ] In Eq. ( 9), [ ] ⋅ S is the internal sign function and can be represented by The basis for the angular spectrum is generated by rotating the reference wave function Eq. ( 7) by the azimuthal angle 0 φ .Therefore, the proposed angular spectrum method uses a rotated MTPO wave function basis in contrast to the conventional angular spectrum method, which uses an infinite plane wave basis.By taking the rotational variable change, 0 0 0 0 cos sin 0 sin cos 0 , 0 0 1 we can obtain the rotated MTPO wave function ( ) and ξ − is written as Here, one important step remains.The MTPO wave function Eq. ( 12) has a non-zero phase when compared with the corresponding uniform plane wave ( ) at the origin point ( ) 0, 0,0 P .The phase compensation term should be multiplied by Eq. ( 12) so that ( ) . The rotated MTPO wave function in Eq. ( 12) is a complex function consisting of amplitude function ( ) g P and phase function ( ) The phase at the origin point can be extracted using the delta function
In addition, for the spectrum basis for ( ) , we take a normal uniform plane wave as the spectrum basis because there is no adapted solution for the MPTO wave function.Consequently, we can complete the construction of the aperiodic angular spectrum basis by solving for x k and y k as , 0 The MTPO wave function is and ξ − can be expressed by ( Finally, the aperiodic angular spectrum integral of the 3D wave field takes the form In Fig. 3, the diffraction field calculated with the proposed aperiodic angular spectrum method in Eq. ( 19) is compared with that obtained by the conventional angular spectrum method in Eq. ( 2).It can be seen that, apart from the bottom aperture, the diffraction field in Fig. 3(b) does not reveal the inter-periodic interference pattern, which is contrasted to the appearance of the strong interference patterns in Fig. 3(b).In the simulation, x T and y T are set at 10μm and a rectangular aperture of 1μm by 1μm is placed at the origin.The 3D crosssectional field visualizations of the conventional angular spectrum method and the proposed aperiodic angular spectrum method are shown in Figs.3(a) and 3(b).The transversal crosssections of the 3D field distributions for the two methods are compared in Figs.3(c) and 3(d).As clearly proven in this comparison, the proposed aperiodic angular spectrum method successfully produces an aperiodic 3D scalar wave field without inter-periodic interference.Also, it is worth noting that, in the proposed aperiodic angular spectrum method, the angular spectrum basis is a solution to the scalar wave equation; therefore, the proposed aperiodic angular representation is also an exact solution to the scalar wave equation.For the evanescent field, the MTPO wave function is not defined.From a physical point of view, the evanescent field is confined to just around the aperture, thus evanescent fields in adjacent computational grids do not considerably interfere with each other, even though the infinitely extended plane wave basis in the evanescent spectrum is used in the numerical modeling.Therefore, the MTPO scattered wave function was only applied to the propagating wave component in the present study.

Numerical visualization of aperiodic holographic optical fields
In this section, the visualization of holographic optical field is presented for the validity and usefulness of the proposed method.In the fields of digital holography and holographic microscopy, the computational synthesis and reconstruction of digital holograms are the essential wave optic topic.The aforementioned inter-periodic interference causes considerable numerical problems in both of synthesis and reconstruction of digital holograms.Normally, to avoid this problem, computation domain is enlarged so that the inter-periodic interference dose not invade the region of interest.It could be a temporary expedient, but cannot be a sound solution in many practical numerical wave optic simulations.
In Figs.4(a) and 4(b), the visualizations of scalar diffraction field generated from a periodic binary amplitude grating on a 11μm width finite aperture is compared by the conventional angular spectrum method and the proposed aperiodic angular spectrum method.density of the wave bases at a specific position is improved.As the measurement point is farther from the x-y plane ( 0 um z = ), the superposition density becomes thinner, meaning that the accuracy of the field representation is degraded.Therefore, the field calculation accuracy is determined by the superposition density of the wave bases and the superposition density is controlled by the radius R of the computational domain 2 x y R + ≤ and varied with the z-directional position.As for the calculation time, the several steps of assignment, rotation, and multiplication of matrices necessary for the calculation of the MTPO wave function is u can affect the calculation time of the proposed method.When the number of sampling points are same, the calculation time is estimated to be about three times longer than that of the conventional method.Moreover, the further investigation for the fast calculation algorithm of the MTPO wave function is a necessity.
However, assuming a specific computation grid volume is set up, we can save the MTPO wave function is u as the basis set by the form of matrices in advance, then the computation load is posed on the linear combination of the pre-calculated bases with specified complex angular spectrum coefficients and the computation efficiency can be enhanced in practice.

Conclusion
In conclusion, we have proposed the aperiodic angular spectrum representation incorporating an analytic thin PML, which might be the first use of a PML for scalar optic waves in Fourier optics.In principle, the proposed method can be widely applied to general optical field models using plane wave bases and Fourier analysis.It is hoped that this method can be extend the method to solve the general inhomogeneous scalar wave equation and visualize aperiodic scalar optical wave field in inhomogeneous media.

Funding Information
This work was supported by Samsung Display Co. Ltd.

Fig. 1 .
Fig. 1.Optical field visualization using the conventional angular spectrum method: (a) a uniform directional plane wave basis in a rectangular computation region, (b) diffraction field ( ), , U x y z for a single slit aperture contaminated by numerical inter-periodic interference, (c) a plane wave basis and (d) diffraction field distribution extended in three-period domain.In (b) and (d), each silt has a width of 5 μm , and the wavelength is set to 633 nm .In (c) and (d), the x-axis period for the unit computational domain is set to 20 μm .
x-y spatial frequency domain, M and N indicate the number of samples in the in the manner illustrated in Fig.2(a).The rotation of the volumetric field of the reference basis with the azimuthal angle 0 φ generates the necessary complete angular spectrum basis in 3D space corresponding to ( ) ( ) in Fig.2(b).The rotational trajectory of the half-infinite barrier, therefore, specifies the computation volume in the cylindrical volume space as shown in Fig.2(c).

Fig. 2 .
Fig. 2. Numerical structure of the proposed aperiodic angular spectrum representation using an MTPO wave function for a half-infinite plane: (a) diffraction wave departing from a finite bottom region with a vertical half-infinite barrier ( 0 z ≤ ≤ ∞ , / 2 x x T = − ), (b) diffraction

Fig. 3 .
Fig. 3. Cross-sectional plots of ( ) , , U x y z 3D field visualization of optical wave fields obtained through (a) the conventional angular spectrum method and (b) the proposed aperiodic angular spectrum method.The x-z cross-section field distributions of each method are plotted in (c) and (d) respectively.The dashed red circle at the bottom indicates the trajectory of the tangent point of the barrier wall.

Fig. 5 .
Fig. 5. Numerical reconstruction of a digital hologram using (a) the conventional angular spectrum method and (b) the proposed aperiodic angular spectrum method.The complex hologram pattern is placed at 0um z = .The target images 'A' and 'B' are focused,