In conventional scalar wave optics [1], optical elements are approximated using a thin-element model that is characterized by a simple transmittance function, a process which is referred to as the thin-element approximation (TEA) [2]. A major limitation of the TEA model is its directional invariance, meaning that the transmittance function of a TEA element is invariant to the direction of an incident optical wave denoted by the two-dimensional complex function T(x, y). In practice, this limitation is acceptable under paraxial conditions and for physically thin elements (within a few wavelengths). In general, however, the modeling of optically thick elements with inhomogeneous refractive index profiles, a high numerical aperture (NA), or a curved optical surface becomes less accurate under the TEA, because in those cases the transmittance function varies with the direction of the incident wave. Thick elements that are susceptible to this problem include high-NA lenses, high-NA freeform lenses, holographic optical elements (HOEs), diffractive optical elements (DOEs), and thick turbid media. If we model these elements using the TEA, significant errors are observed in the optical-field distribution. Although these elements are important, their accurate modeling is beyond the capability of current wave optics. In practice, vectorial electromagnetic analysis methods are well established, but their application to mesoscopic optical elements is practically impossible, due to limitations in computational resources. A scalar wave-equation solver is thus a necessity. Geometric optics is widespread in the field of optical design, but it does not consider diffraction. When the size of an optical element is mesoscopic or microscopic, and thus affected by the diffraction effect, geometric optics breaks down. For this reason, simulations based on scalar waves, which offer acceptable computational efficiency at these scales, is necessary, particularly for modern applications.

In previous research, the beam-propagation method (BPM) [3–5] and wave-propagation method (WPM) [6] have been proposed to deal with nonparaxial diffraction fields and thick elements, but these approaches struggle with the complex optical structures of high-contrast or gradient inhomogeneous index profiles [6]. To characterize complex optical systems consisting of thick, volumetric inhomogeneous elements, we propose a frequency-domain method capable of calculating a full transmittance function with directional T(x, y, k_inc,x, k_inc,y; λ), which reflects the essential properties of the optical functions of thick, high-NA elements. The computation of the full transmittance function is essential for the design, analysis, and optimization of these elements [7–9]. For instance, an imaging system consisting of a thick lens (such as the lens of the human eye) has a directionally varying point-spread function (PSF). To compensate for all of the optical aberrations of the human eye’s lens, the full transfer function of the ocular optical system should be calculated or measured [10, 11]. There is also no suitable method for large-scale modeling of HOEs with complicated grating profiles, which have been identified as essential elements in augmented-reality (AR) displays [12, 13].

In this paper, we investigate the essence of the Fourier modal method (FMM) [14–16], and propose the scalar Fourier modal method (SFMM) as a new numerical framework for the analysis of the scalar wave equation:

(1)∇2U(r)+k02ε(r)U(r)=0

where k₀ and ε(r) are the wavenumber and the three-dimensional permittivity function, respectively. In theory, Fourier modal analysis identifies the Bloch eigenmodes of the permittivity system ε(r), and constructs the total wave field distribution using the linear superposition of these eigenmodes. Because a Bloch eigenmode cannot be obtained from an explicit function, we take the pseudo-Fourier series representation [14], which is much more compatible with the angular spectrum representation of Fourier optics. The pseudo-Fourier series representation of an optical field in an inhomogeneous medium can be viewed as a generalized angular spectrum representation, in the context of Fourier optics or a scalar version of the FMM.

This paper is organized as follows. Section 2 establishes the numerical framework for the proposed SFMM, and numerical examples are presented in Section 3 with the introduction of absorbing boundaries. Concluding remarks are provided in Section 4.

The construction of the SFMM follows steps similar to those in the conventional vectorial FMM [16]. In the first step, the model function of the optical field based on the pseudo-Fourier series converts the scalar wave equation of Eq. (1) to an algebraic eigenvalue equation.

According to the conventional structure-modeling method, we take a staircase representation of the 3D permittivity profile ε_r(r) along the optical axis, in the form of a multilayered structure, as illustrated in Fig. 1 [16]. Details of the propositions related to the structural modeling are presented in Section 3. For a single thin layer with permittivity profile ε_r(x, y) in the range z₋ ≤ z ≤ z₊, the wave equation is

Figure 1.Staircase representation method of the 3D permittivity profile and related variable for each slab. (a) Conceptual diagram of staircase modeling and (b) the variables related to each slab, according to the direction of the incident wave. An arbitrarily shaped element is sliced into multiple layers. Each slab has a permittivity profile that varies only in the x and y directions. Using the boundary conditions between the slabs, the coupling coefficients C can be derived for each slab.

(2)∇2U+k02εr(x,y)U=0.

The models of the permittivity profile and the field representation take the form of a Fourier series and a pseudo-Fourier series respectively:

(3a)εr(x,y)=∑m=−MM∑n=−NN εm,nej(mGxx+nGyy),

(3b)U(x,y,z)=ej(kx,0x+ky,0y+kzz)∑ m=−MM∑n=−NNUm,nej(mGxx+nGyy),

where ε_m,n and U_m,n are the Fourier coefficients of the two-dimensional relative permittivity profile and the wave function respectively. G_x and G_y are the x- and y- directions primitive reciprocal vectors respectively, i.e. G_x = 2π / T_x and G_y = 2π / T_y. T_x and T_y are the x- and y- directions periods of the permittivity profile respectively, while M and N denote the x- and y- directions truncation numbers for the Fourier harmonics respectively. The propagation wavevector k = (k_x,0, k_y,0, k_z) is specified by the transverse wave-vector components (k_x,0, k_y,0).

Here the longitudinal wave-vector component k_z is an eigenvalue extracted from the scalar wave equation. The substitution of Eqs. (3a) and (3b) into Eq. (2) transforms the wave equation into an algebraic eigenvalue matrix equation for ε_m,n and U_m,n:

(4)  ∇2U+k02εr(x,y)U =∑ m=−MM∑n=−NNU m,n[− (k x,0+mGx)2− (k y,0+nGy)2− (k z,0)2]e j((k x,0+mGx)x+(k y,0+nG y)y+kzz) +k02∑ m=−MM∑n=−NN ∑s=−M M ∑t=−NN ε m−s,n−t U s,t e j((k x,0+mGx)x+(k y,0+nG y)y+kzz)=0.

Equation (4) is set so that it is satisfied for any x and y. The following eigenvalue equation is then derived as:

(5)[(kx,0+mGx)2+(ky,0+nGy)2]Um,n−k02∑s=−MM∑t=−NN εm−s,n−t Us,t=−(kz,0)2Um,n.

The above eigenvalue equation can be expressed in a compact matrix form:

(6)(K__+k02ε__)U_=−( kz,0)2U_.

The detailed expressions for K__, ε__, and U_ are provided in the Appendix. Comparing the SFMM to the vectorial FMM in terms of the size of the eigenvalue problem, the computational dimension of the SFMM is 1/4 that of the FMM, which means that the required memory quantity and calculation speed are reduced by a factor of 4. More importantly, given the same computational resources, the SFMM can deal with larger structures that the FMM cannot afford. By solving the eigenvalue equation, we can obtain the eigen-pair and form the Bloch eigenmode for Eq. (3b). The solution to this algebraic eigenvalue equation is interpreted as a pseudo-Fourier series representation of the optical Bloch eigenmodes of the scattering system ε(r).

The Bloch eigenmodes obtained for each region can be classified into two groups of positive (U^(g)+) and negative (U^(g)−) eigenmodes [16]:

(7)U(g)±=ejkz(g)±(z−z∓)∑ m=−MM∑n=−NNU m,n (g)±ej(kx,mx+ky,ny),

where k_x,m and k_y,n denote k_x,m = k_x,0+ mG_x and k_y,n = k_y,0+ nG_y, respectively. The classification of positive and negative eigenmodes is according to the rule that the amplitude of the eigenmode does not grow along the propagation direction. An eigenmode decreasing along the positive z direction is a positive mode, while the eigenmode decreasing along the negative z direction is a negative mode. The total optical scattering field is determined by the linear superposition of the Bloch eigenmodes with the coupling coefficients. The wave function in the thin layer is represented by a linear combination of these eigenmodes:

(8)Ua(z)=∑g=1M+C a,g+U (g)+(z)+∑g=1M−C a,g−U (g)−(z),

where Ca,g+ and Ca,g- denote the coupling coefficients for the positive and negative eigenmodes respectively. Here the subscript a is used for left-to-right characterization [16]. Right-to-left characterization is presented below with the subscript b. M⁺ and M⁻ represent the number of positive and negative eigenmodes respectively.

To determine the coupling coefficients, we employ the scattering-matrix (S-matrix) method, which is an essential element of the FMM. Figure 1(b) presents a schematic diagram of left-to-right characterization for a block of interest. The incident wave function is denoted by the incident wave U→i, the reflected wave U→r, and transmitted wave U→t. Those wave functions are represented by the following pseudo-Fourier series:

(9a)U→i=∑m=−MM∑n=−NN U→i,m,nej(kx,m,nx+ky,m,ny+kz,m,n(z−z−)),

(9b)U←r=∑m=−MM∑n=−NN U←r,m,nej(kx,m,nx+ky,m,ny−kz,m,n(z−z−)),

(9c)U→t=∑m=−MM∑n=−NN U→t,m,nej(kx,m,nx+ky,m,ny+kz,m,n(z−z+)).

Under the scalar approximation, the energy of an optical wave is completely conserved at a lossless dielectric interface, stating that the sum of the squares of the reflection and transmission coefficients is equal to the square of the input field coefficient, when the loss in both regions can be ignored. The coupling coefficients Ca,g+ and Ca,g- are determined by applying the boundary condition, which states that the wave field should be continuous and smooth at the interface. From a vector-field perspective, the transverse electric field under TE polarization has the same wave-field characteristics. While the electric field under TM polarization exhibits a discontinuity at the interface, the magnetic field profile is continuous and smooth. From a scalar-field perspective, the wave function U and its derivative with respect to the propagation direction (here the z axis) dU / dz is supposed to be continuous across the interface between the regions, which is the same boundary condition as in the Schrödinger equation. The continuity of U and dU / dz should be secured at the interface [17]. For a thin-slab structure with the left and right boundaries located at z = z₋ and z₊ respectively, at z = z₋ the boundary condition is given by

(10a) U→i+U←rz=z−= Uaz=z−,

(10b)dU→idz+dU←rdzz=z−=dUadzz=z−,

The substitution of Eqs. (9a), (9b), and (9c) into Eqs. (10a) and (10b) leads to

(11a) ∑m=−MM∑n=−NN (U→i,m,n+U←r,m,n)ej(kx,m,nx+ky,m,ny)= ∑m=−MM∑n=−NN ∑g=1M+ Ca,g+ ejkz(g)+(z−−z−) Um,n(g)++∑g=1M − Ca,g− ej kz(g)−(z−−z +) Um,n(g)−ej(k x,m,nx+k y,m,ny),

and

(11b) ∑m=−MM∑n=−NN (jkz,m,nU→i,m,n−jkz,m,nU←r,m,n)ej(kx,m,nx+ky,m,ny) =∑m=−MM∑n=−NN ∑g=1M+ jkz(g)+Ca,g+ej kz(g)+(z−−z−) Um,n(g)++∑g=1M− jkz(g)−Ca,g−ej kz(g)−(z−−z +) Um,n(g)−ej(kx,m,nx+ky,m,ny).

The term z₋ − z₋ is used above instead of 0, to make it easier to follow the process. Because the scalar wave function is represented in the linear matrix of Eq. (6), the boundary condition at z = z₋ can also be represented in an algebraic form. The boundary-matching condition can be read as

(12) II jk z,m,n −jk z,m,nU→ i,m,n U← r,m,n = Um,n(1)+ j kz(1)+ Um,n(1)+ ⋯⋯ Um,n(M+)+j kz(M+)+ Um,n(m+)+ Um,n(1)−ejkz(1)−(z−−z+) jkz(1)− Um,n(1)−ejkz(1)−(z−−z+) ⋯ ⋯ U m,n (M−)−e jkz(1)−(z−−z+) jkz (M−)− U m,n (M−)−e jkz(1)−(z−−z+) [Ca+][Ca−],

where [Ca+]=( Ca,1 +  ⋯  Ca,M + + )t,[Ca−]=( Ca,1 −  ⋯  Ca,M + − )t, and −M ≤ m ≤ M and −N ≤ n ≤ N. At z = z₊, the boundary condition is given by

(13a) Uaz=z+=U→tz=z+,

(13b)dUadzz=z+=d U→tdzz=z+

The corresponding matrix form of the boundary condition is given by

(14) II jkz,m,n −jkz,m,n U→ t,m,n0 = Um,n(1)+ej kz (1)+(z+−z−)j kz(1)+Um,n(1)+ej kz (1)+(z+−z−) ⋯ ⋯ Um,n(M+)+ej kz(M+)+(z+−z−) j kz(M+)+ Um,n(M+)+ej kz(M+)+(z+−z−) Um,n(1)− jkz(1)− Um,n(1)− ⋯ ⋯ U m,n (M−)− jkz (M−)− U m,n (m−)− [Ca+][Ca−].

By solving Eqs. (12) and (14), the reflection coefficient R←__, the transmission coefficient T→__, and the coupling coefficients C__a+ and C__a- can be obtained [16]. In a similar manner the right-to-left characterization is established, in which U←__ is the incidence coefficients, R→__ is the reflection coefficients, T←__ is the transmission coefficients, and C__b+ and C__b+ are the coupling coefficients. The boundary conditions at z = z₊ and z = z₋ can then be derived in a linear algebraic form:

(15a) II jkz,m,n −jkz,m,n0 U← t,m,n = Um,n(1)+j kz(1)+Um,n(1)+ ⋯ ⋯ Um,n(M+)+ j kz(M+)+ Um,n(m+)+ Um,n(1)−ejkz(1)−(z−−z+) jkz(1)− Um,n(1)−ejkz(1)−(z−−z+) ⋯ ⋯ U m,n (M−)−e jkz(1)−(z−−z+) jkz (M−)− U m,n (M−)−e jkz(1)−(z−−z+) [Cb+][Cb−],

(15b) II jkz,m,n −jkz,m,n U→ r,m,n U← i,m,n = Um,n(1)+ej kz (1)+(z+−z−)j kz(1)+ Um,n(1)+ej kz (1)+(z+−z−) ⋯ ⋯ Um,n(M+)+ej kz(M+)+(z+−z−) j kz(M+)+ Um,n(M+)+ej kz(M+)+(z+−z−) Um,n(1)− jkz(1)− Um,n(1)− ⋯ ⋯ U m,n (M−)− jkz (M−)− U m,n (m−)− [Cb+][Cb−],

where and [Cb+]=( Cb,1 +  ⋯  Cb,M + + )t and [Cb−]=( Cb,1 −  ⋯  Cb,M + − )t. A detailed theoretical explanation of the matrix formalism can be found in [16]. The scattering-matrix coefficients and coupling coefficients given by Eqs. (12), (14), (15a), and (15b) provide complete information for a single block. The S-matrix method states that the total response of a system composed of several blocks is obtained by combining the coefficients of neighboring blocks [16]. From the associative rule of S-matrices, the block S-matrix with the convolution of the multiblock M^(1,N) is obtained by

(16a)M(1,N)=S(1)*S(2)*⋯*S(N−1)*S(N).

The associative rules enable the parallel computation of the internal coupling coefficients:

(16b)(Ca,Cb)=( C__a(1), C__b(1))*( C__a(2), C__b(2))*⋯*( C__a(N−1), C__b(N−1))*( C__a(N), C__b(N)).

For instance, when the incident field is coming from the left side, the reflection field in the left half-infinite block and the transmitted field in the right half-infinite block can be calculated respectively by

(17a)U←r= R←_ _(0,N+1)U→i,

(17b)U→t= T→_ _(0,N+1)U→i,

The total field in the left half-infinite block, where z is less than 0, can be visualized by

(18a)U=U→iej(kxx+kyy+kz,(0)(z−z−))+U←rej(kxx+kyy+kz,(0)(z−z−)) for z≤0,

and the fields of the k^th layer in the multiblock region (0 ≤ z ≤ l_N) can be calculated by

(18b)U=U=∑g=1M+C a,(1),g+U (1) (g)++∑g=1M+ C a,(1),g +U (1) (g)− for 0≤z≤l 1,1 U=∑g=1M+C a,(2),g+U (2) (g)++∑g=1M+ C a,(2),g +U (2) (g)− for l 1,1≤z≤l 1,2 ⋮U=∑g=1M+C a,(N),g+U (N) (g)++∑g=1M+ C a,(N),g +U (N) (g)− for l 1,N−1≤z≤l 1,N

where l_n−1 and l_n denote the z axis coordinates of the left and right ends of the nth layer respectively. Finally, the field of the right half-infinite block, where z is greater than l_N, is given by

(18c)U=U→iej(kxx+kyy+kz,(N+1)(z−z+)) for z≥lN.

Because we already have the same coefficients for the opposite direction, it is straightforward to visualize the right-to-left case. This outcome is the key result of this paper, and we hereby refer to this overall approach as the SFMM.

In the modeling of thick optical elements with smooth surfaces, the conventional staircase approximation of the structure can generate physical errors, and may be limited by computational inefficiency. Figure 2(a) compares the conventional and proposed staircase approximation methods for a smooth optical element. Using the conventional method, unexpected physical phenomena can be observed, such as diffraction at the edges of the slabs, which will lead to errors if the number of slabs is not high enough to model the curved surface.

Figure 2.Proposed modeling method and sampling diagram. (a) Schematic diagram of the proposed staircase modeling method using a gradient index. (b) Diagram of the sampling scheme used to calculate the effective-refractive-index profile of the lens. f(x) is the curve of the air-glass boundary. Along this boundary function, equidistant sampling along the x axis is used to determine the thickness of each slab. The effective refractive index n(x) is calculated as the ratio of the area of each index region.

Figure 3 presents a numerical experiment for the calculation of the field, in which the number of staircase layers in a convex lens profile is varied. The diameter of the lens aperture is 1 mm, the focal length of the lens is 12 mm, and the curvature of the lens is 24 mm. The refractive index of the lens is 1.5, and that of the surrounding material is 1.0. The number of Fourier harmonics is set to 301 in each direction, which means that the truncation numbers are set to 150. A plane wave with a wavelength of 633 nm is normally incident from the left.

Figure 3.The amplitude of the field, after passing through the lens. The lens is modeled with (a) 10 layers, (b) 140 layers, and (c) 200 layers, using the conventional method.

The staircase structure with only 10 layers exhibits a significant edge-diffraction effect, losing the essential optical properties of the convex lens [Fig. 3(a)]. In the 140-layer model of the same convex lens, some defects caused by the edge diffraction are still observable [Fig. 3(b)]. More than 200 layers are required to model the lens with convergent results [Fig. 3(c)]. A staircase approximation with sharp edges requires significant computational resources for smooth surface modeling because the z-directional layer thickness approaches the wavelength scale. However, this heavy computational burden can be effectively reduced by approximating the sharp edges.

The proposed approximation is based on the assumption that a smoothly varying permittivity profile in the transverse direction at the edges of the slabs, i.e. a thin gradient-index (GRIN) profile, reduces the edge effect and leads to a physically more accurate model with a considerably lower number of staircase layers. This thin-GRIN-layer-model concept is schematically illustrated in Fig. 2(a). We take sampling points on the curved surface of the target structure at equidistant intervals, Δx along the x axis, and the nonuniform thickness of each layer Δz_n is determined by the surface profile [Fig. 2(b)]. The permittivity profile around the edges is approximated by the effective-refractive-index function n_eff (x). Bilinear interpolation is used to calculate the effective-refractive-index profile n_eff (x) in this region. For instance, the effective index of the first layer in Fig. 2(b) is calculated by

(19)neff(x)={n1(f(x)−f(x1))+n2(f(x2)−f(x))}Δz1,

where n₁ and n₂ are the refractive indices of the structure and free space, f(x) is the function of the curved surface of the target structure, and Δz₁ is the thickness of the first layer. Note that each slab can have an arbitrary permittivity profile in the transverse direction, in theory. The spatial resolution of the effective-index profile depends on the number of Fourier harmonics, which should be selected carefully to accurately represent the spatial variation of the index profile. Figure 4 presents the field-calculation results for the same target structure and the same sampling conditions as used in Fig. 3. Compared to Fig. 3(b), which is modeled with 140 layers, the noise due to diffraction at the slab’s boundary is not observed even in the 10-layer model [Fig. 4(a)]. This illustrates the superiority of the proposed model, in terms of both model accuracy and computational efficiency. For comparison, the field distribution obtained using the conventional TEA model of the structure is presented in Fig. 4(c). In this example, the NA of the lens is low enough to be approximated by the TEA. The field calculated using the SFMM is close to that of the TEA model.

Figure 4.The amplitude of the field, after passing through the lens. The lens is modeled with (a) 10 layers and (b) 20 layers, using the proposed modeling method. (c) The numerical results calculated for the same field using the angular spectrum method with the thin-element approximation (TEA).

The TEA can only approximate light transmittance, while the SFMM can model light absorption, diffraction, and multiple reflections, for an arbitrarily complex refractive index. In Fig. 5, a cylindrical lens in free space is modeled using the SFMM with various materials. In this example, the diameter of the lens is 5 µm, and 152 layers and 201 Fourier harmonics are used for the simulation. The angle of the incident plane wave is set to 30°, and the wavelength is 633 nm. In Figs. 5(b) and 5(c), the absorption due to the complex refractive index is observable by the field inside the lens, and the diffraction is also observable by the field behind the lens. On the front face of the lens, the interference patterns between the incident light and the reflected light on the glass-air boundary is shown in Figs. 5(b) and 5(c). The TEA is often used to model optics based on the scalar diffraction theory, such as in the angular spectrum method. The TEA does not consider multiple reflections, which is essential for some applications. The SFMM is thus an efficient tool for calculating a scalar field with the rigorous consideration of multiple internal reflections.

Figure 5.The real value of the field, with an incidence angle of 30°, after passing through cylinders of various materials: (a) glass (n = 1.5), (b) silicon (n = 3.877 + 0.001i), and (c) gold (n = 0.196 + i3.258). Diffraction behind the shadow of the object is clearly seen in (b) and (c), and interference patterns with the reflected wave are observable.

Another advantage of the SFMM is that it can obtain the directionally variant transfer function, unlike the TEA method. As mentioned in the introduction, practical optics often have a high NA and nonnegligible thickness, leading to directionally variant characteristics that need to be modeled accurately. To confirm these characteristics of the SFMM, a bi-convex lens of sufficient thickness and a high NA is modeled using the SFMM and TEA. The curvature radius of the bi-convex lens is 66.67 µm, and the focal length of the approximated thin lens is also 66.67 µm. Because the diameter of the lens is 100 µm, the NA of the lens is 0.60. To realize the simulation of the aperiodic optical structure, the additional numerical technique of an absorbing boundary condition (ABC) [16] is used on both side boundaries. Figure 6(a)–6(c) show field distribution by SFMM calculation, and (d)–(f) show field distribution by TEA. Comparing Figs. 6(c) and 6(f), the difference between the SFMM and TEA is noticeable. In Fig. 6(c), a different focal length for the various angles of the incident wave is confirmed. The simulated convex lens has a field-curvature aberration whose focal length varies depending on the direction of the incident plane wave, which the SFMM can model. Because the transfer function of TEA is directionally invariant, the focal length cannot vary with the angle of incidence.

Figure 6.The amplitude of the field passing through a bulky lens, simulated by (a)–(c) SFMM and (d)–(f) TEA. Field-curvature aberration, in which the focal length varies depending on the direction of the incident wave, can be observed in the SFMM results (c), but not in the TEA results (f). In addition, reflected waves on the curved surface of the lens can only be observed in the SFMM.

In conclusion, we have proposed the SFMM for wave-optic optical-element modeling. The SFMM is an appropriate method for the calculation of a full transfer matrix, for practical applications and for fundamental research. In addition to thick optical elements, HOEs (which have been spotlighted as an optical element for AR displays [12, 13]) and light scattering and transportation through disordered media [8, 9] may benefit from the application of the SFMM.

For matrix operations, the Toeplitz matrix is means of

(A1a) ε__=ε¯ 0ε¯ −1⋯ε¯ −2Mε¯ 1ε¯ 0ε¯ −2M+1⋮⋮ε¯ 2Mε¯ 2M−1⋯ε¯ 0,

where

(A1b)ε¯M=ε¯ m,0ε¯ m,−1⋯ε¯ m,−2Nε¯ m,1ε¯ m,0ε¯ m,−2N+1⋮⋮ε¯ m,2Nε¯ m,2N−1⋯ε¯ m,0.

(A2a) K__=k¯ 0k¯ −1⋯k¯ −2Mk¯ 1k¯ 0k¯ −2M+1⋮⋮k¯ 2Mk¯ 2M−1⋯k¯ 0,

where

(A2b)k¯m=( k x,m)2+( k y,0)2( k x,m)2+( k y,−1)2⋯( k x,m)2+( k y,−2N)2( k x,m)2+( k y,1)2( k x,m)2+( k y,0)2( k x,m)2+( k y,−2N+1)2⋮⋮( k x,m)2+( k y,2N)2( k x,m)2+( k y,2N−1)2⋯( k x,m)2+( k y,0)2,

and, k_x,m and k_y,n denote k_x,m = k_x,0 + mG_x and k_y,n = k_y,0 + nG_y in Eq. (5) respectively. The column vectors for the Fourier coefficients are expressed respectively by

(A3)U_={( U−M,−N,⋯, U−M,N)( U−M+1,−N,⋯, U−M+1,N) ⋯ ( UM,−N,⋯, UM,N)}t

This study was supported by the National Research Foundation of Korea (NRF) (NRF-2019R1A2C1010243).