Modeling of Finite Beams in Periodic FDTD Methods

This article presents a finite-difference time-domain (FDTD) framework that accelerates the modeling of the interaction of infinite periodic structures with practical sources (such as Gaussian or Bessel beams) by exploiting periodic boundary conditions (PBCs). The technique involves the introduction of arbitrary beam-generating sources in periodic simulations via the array-scanning method (ASM). The proposed method is both efficient and practical, being naturally broadband and highly parallelizable. We present examples of Gaussian beams in free space and interacting with grating geometries, efficiently replicating the results of multi-period brute-force simulations.

The finite-difference time-domain (FDTD) algorithm has emerged as a popular simulation technique.The FDTD method is inherently broadband and does not require computationallyexpensive matrix operations, making it an attractive alternative to frequency-domain methods such as the finite-difference frequency-domain (FDFD) algorithm.
Periodic structures are efficiently modeled in FDTD by simulating a single unit cell terminated by periodic boundary conditions (PBCs) in the periodic directions [8].Fields simulated in models with PBCs propagate as though the unit cell were infinitely replicated in a periodic fashion.
Simulations with PBCs duplicate sources in each periodic unit cell, thereby preventing the introduction of physical sources of finite extent.Additionally, since sources in each unit cell are identical (up to a phase shift), arbitrary line sources spanning multiple periodic unit cells cannot directly be realized.In particular, these limitations hinder simulations of electromagnetic beams-which are practically realizable in a lab setting-interacting with infinitely-periodic geometries.
In the literature, beams have been produced in periodic FDTD domains by simulating a finite number of periods without using PBCs [9], [10], [11], [12], [13].In this case, sources are not infinitely replicated due to PBCs, and beams can be generated without difficulty.Although the number of periods simulated is finite, it is chosen to be high enough to produce fields approximately equal to those generated in an infinitely periodic version of the medium.The main disadvantage of this technique is the prodigious computational cost associated with the large size of the simulated domain.
Another technique to simulate beam interactions with periodic structures in FDTD is the plane-wave expansion method (PWM), which involves decomposing an analytical beam formulation into a sum of constituent plane waves using a spatial Fourier transform.Each plane wave, and its interaction with the periodic geometry, can be simulated in a single unit cell terminated by PBCs [11], [14].Alternatively, FDTD can be used to measure reflection and transmission coefficients to model the scattering behavior of the constituent plane waves [15].An inverse Fourier transform of the scattered fields across different spatial frequencies approximates the beam interaction with the periodic structure [16].
The inverse Fourier integral used to produce the beam interaction generally spans the entire wavenumber plane when evanescent waves play a role, and the integral may therefore be computationally expensive to calculate [17].Even within the range of propagating wave numbers, the quadrature convergence rate of the Fourier integral depends on the geometry.In addition, since the PWM uses fields directly from the frequency domain, frequency-domain simulation techniques are more naturally suited for use in PWM implementations than are time-domain methods, such as FDTD.When modeling monochromatic waves, FDTD simulations need to run for many time steps before transient broadband signals decay.
The array-scanning method (ASM) is a technique used to remove periodic image sources outside the central unit cell in simulations terminated by PBCs [18], [19], [20].It was recently noticed that the ASM could be generalized to add sources to or remove sources from arbitrary unit cells [20].Using this observation, this article presents a novel algorithm that employs the ASM to introduce continuous 0018-9480 © 2023 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
sources-potentially spanning multiple unit periods-into FDTD simulations, thereby providing a means to simulate arbitrary, finite-sized practical sources, such as beams, interacting with periodic media.Section II describes the total-field/scattered-field (TF/SF) method, its relationship to the source-equivalence theorem, and how it can be used to produce beams in the FDTD algorithm.Section III outlines the background and implementation of FDTD PBCs, and explains how the ASM is used to produce finite-sized sources in periodic simulations.A discussion about efficient integration of the ASM integral is presented.The algorithm for simulating beams in periodic structures is derived in Section IV.Finally, Sections V-A and V-B provide examples to demonstrate the efficacy and efficiency of the algorithm of Section IV.

II. TF/SF SOURCE CONDITION
The TF/SF source formulation [21] provides a technique to efficaciously introduce arbitrary electric and magnetic fields into a computational domain.
The TF/SF scheme bifurcates the computational domain into a "total field" (TF) region and a "scattered field" (SF) region.The TF/SF boundary generates fields that enter the TF region.In the absence of a structure within the TF region, fields are not present in the SF region.Fields which are perturbed by the presence of a structure in the TF region will create scattered fields in the SF region.
The TF/SF formulation is in fact a discretized formulation of the equivalence principle [1], [22].According to the source equivalence principle, sources along a closed boundary ∂ can be chosen to produce any desired field within the interior of the boundary , and any field exterior to the boundary ¯ , so long as both fields satisfy Maxwell's equations in those regions.Let us define the electric and magnetic fields on and ¯ as E , H and E ¯ , H ¯ , respectively.The desired fields can be produced by the electric and magnetic currents J s and M s as for all r ∈ ∂ (here n is the outward-facing unit vector normal to ∂ ), as in Fig. 1.
When E ¯ = H ¯ = 0, (1a) and (1b) become which corresponds precisely to the case of TF/SF: the domain is the total-field (TF) region, the domain ¯ is the scatteredfield (SF) region, and the surface ∂ is the TF/SF boundary.Any fields whose values are known along ∂ -either analytically, or through a previous simulation-can therefore be simulated in FDTD simulations through the TF/SF method.

A. Simulation of Gaussian Beams Using the TF/SF Method
Analytical formulations of beams can be used to calculate the currents used by the TF/SF technique using (2a) and (2b).shows how the fields produced by the point source can be generated throughout , while maintaining an absence of fields in ¯ .Electric and magnetic currents along ∂ are defined using tangential field components of the fields via (1a) and (1b).A FDTD discretization of the current densities J s and M s leads to the TF/SF method.
To exemplify this concept, let us consider generating a general Gaussian beam using the TF/SF method.
An electric or magnetic field phasor ⃗ U(x, y, z) with a phase proportional to z can be written as for some ⃗ U ⊥ (x, y).The tilde ( • ) is used to denote phasor quantities, where a time dependence of e + jωt has been suppressed.Since ⃗ U satisfies the Helmholtz equation the function ⃗ U ⊥ satisfies the paraxial Helmholtz equation where [23].Gaussian beams arise as approximate solutions to the paraxial Helmholtz equation.Gaussian beams have a narrower bandwidth than Bessel beams (which solve the Helmholtz equation exactly) but are not strictly divergence-free.
A general x-polarized Gaussian beam propagating in the z-direction is defined by the phasor formula [23] where A 0 equals the beam amplitude at the origin, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 2. Cross section of a z-directed beam generated by a source (shown as a dashed green line) defined in the plane z = z 0 .The beam is produced in the TF region, which occupies z > z 0 .The volume z < z 0 defines the SF region, which remains devoid of fields as long as no scattering objects are present in the TF region.
is the radius of curvature of the wavefronts and ζ (z) = arctan(z/z 0 ) is the Gouy phase.The constant z 0 denotes the Rayleigh range, which equals the distance from the beam waist along the propagation direction at which point the cross-sectional area of the beam has doubled.
The beam can be angled in an arbitrary direction using the formula where T is the rotation matrix and ϕ and θ are azimuthal and polar angles, respectively.The magnetic field H R may be deduced via Faraday's law.
The plane z = z 0 is the most natural placement of a TF/SF boundary to generate a beam that travels along θ ≈ 0, since it is approximately orthogonal to the beam axis.In this case, the TF region is defined as the domain z > z 0 and the SF region is defined as the domain z < z 0 , as shown in Fig. 2. Defining n = +ẑ, (2a) and (2b) provide the TF/SF current sources to generate the Gaussian beam In the above formulas, ω is the electromagnetic field angular frequency and µ is the magnetic permeability.
Although J s and M s are defined throughout the entire plane z = z 0 , this cannot be implemented in practice.Nonetheless, this restriction does not pose a problem, since the source magnitudes decay rapidly away from the center of the source in the z = z 0 plane.Consequently, sources of finite extent can be used to effectively generate a Gaussian beam, so long as the source magnitudes have dropped to nearly zero at the source boundaries.
A broadband Gaussian beam pulse can be generated by a Blackman-Harris window function in a straightforward manner.While it is turned on, a Blackman-Harris pulse is comprised of a linear combination of discrete monochromatic sources.The Blackman-Harris window is notable for its ability to be turned on instantaneously without needing a time delay, as is required for Gaussian pulses [24].
To produce the Blackman-Harris window pulse, define coefficients and the period T b = 9.74/ω bw , where ω bw = 2π f bw is the desired half width at half maximum (HWHM) of the excitation.Given an electric or magnetic current source phasor S(ω), the broadband pulse can be expressed in the time domain where ω n = n/T b + ω c and where ω c = 2π f c is the central angular frequency [24].

III. PERIODIC BOUNDARY CONDITIONS (PBCS)
Periodic boundary conditions (PBCs) can be implemented in the usual FDTD scheme by applying a constant phase shift to waves in the periodic directions derived using the Bloch-Floquet theorem.The Bloch-Floquet theorem states that for an electric or magnetic field phasor U(r, ω) at any point r = (x, y, z) T in a 3-D periodic geometry Here, d = n x d x x + n y d y ŷ + n z d z ẑ is the general lattice vector of the structure, where d x , d y and d z are the periods in the x, y and z directions, respectively, and n x , n y and n z are integers.The vector k is the reciprocal (Bloch) wave vector of the fields.Assuming a frequency-independent k, a time-domain analog of the Bloch-Floquet theorem is derived via an inverse Fourier transform of (12) with respect to frequency This formula relates fields one period apart, and can therefore be used as a PBC [8], [25], [26].FDTD simulations with PBCs are commonly used to produce band diagrams, extract attenuation constants, and determine reflectivity of periodic structures [8].
The fact that (13) produces complex-valued fields does not pose a problem, since the linearity of the FDTD algorithm ensures that real and imaginary fields are simulated separately (except at the PBCs) and satisfy Maxwell's equations.The imaginary fields are necessary for the PBCs to enforce the complex phase shift in the frequency domain (12).Further discussions on the nature of these complex fields are found in [8].
To illustrate how periodicity is integrated into FDTD, consider for simplicity a 2-D FDTD discretization of a single unit cell of a structure that is periodic in the z direction, as in Fig. 3.A column of auxiliary H y nodes is added outside of the unit cell at z = d z + z/2.Using (13), the boundary field nodes are related using the Bloch-Floquet theorem as follows: The field E x (z = 0) should be updated via (14a) whenever E x (z = d z ) changes.Likewise, the H update equation (14b) should be invoked whenever H y (z = z/2) changes.The real and imaginary parts of these fields describe physical electromagnetic fields and they can thus be simulated by FDTD.The simulations of the real and imaginary parts are disjoint except at the periodic boundaries.
Although only one period with PBCs needs to be simulated to analyze an infinite periodic structure, (13) can be used to determine the fields in any period.

A. Array-Scanning Method (ASM)
Periodic boundary conditions replicate both the geometry and current sources in the unit cell simulated.Fields produced from simulations with PBCs thus correspond to fields generated by sources in each unit cell.The "image" sources S(r, t) in each cell are related to each other by a Floquet-like phase shift where d = n x d x x + n y d y ŷ + n z d z ẑ and k is the wave vector set by the PBCs.Consequently, the fields generated in simulations with PBCs, U ∞ (r, t, k), are a superposition of the fields generated by the sources in the primitive (0, 0, 0) period, U(r, t, k), phase-shifted and translated to each unit cell Simulating finite sources interacting with an infinitelyperiodic geometry is often desired.Suppose we want fields produced by a single electric or magnetic source A(r, t) in period (n x , n y , n z ) alone, without image sources.Then, a complex-valued source S(r, t) = A(r, t)e j k•d can be placed in the simulated unit cell (corresponding to cell (0, 0, 0)).
The fields in cell (N x , N y , N z ) produced by sources in the (n x , n y , n z ) cell alone are determined by by orthogonality.In this equation, The integration technique of ( 17) is known as the arrayscanning method (ASM) [8], [27].Note that the ASM integral removes the complex fields introduced by the PBCs, returning only the physical real component of the fields.
The ASM has been successfully implemented within the FDTD framework [18], [19], [20].The fields U ∞ (r, t, k) can be derived from a single FDTD simulation when the PBC wave vector is set to k.As such, the integral ( 17) can be approximated by running multiple independent simulations with PBCs initialized with different wave vectors.
The periods d x , d y , and d z -describing the size of the simulated unit period-are usually chosen to correspond to the smallest periods of the geometry under study in each Cartesian direction, but this does not need to be the case.The simulated period can be defined as a "supercell," comprising multiple primitive periods of the geometry.Increasing the size of the simulated volume increases the cost of each individual simulation, but commensurately decreases the number of simulations needed to accurately calculate (17).Thus, simulating a supercell may be practical when the number of available processing units limits the degree of parallelization.
1) ASM Quadrature Techniques: The integral in ( 17) is often approximated by numerical techniques.The midpoint and trapezoidal rules, with equally spaced abscissas, have been shown to be especially effective at numerically integrating smooth periodic functions, typically converging exponentially or faster [28], [29].These quadrature techniques are appropriate choices when integrating the ASM field (17) [17].Although both the trapezoidal and midpoint rules exhibit similar convergence patterns, we show in this section that the trapezoidal rule operates more efficiently than the midpoint rule when the number of sufficient integration points has yet to be determined.
The midpoint rule for integrating a 1-D ASM integral is where N is the number of integration points, x is the integrand of ( 17) in Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
one dimension.The trapezoidal rule is When N points are used in the approximation of the integral, unwanted sources in N − 1 cells to either side of cell n x are removed in the periodic direction exactly, by using the midpoint or trapezoidal rules [17].Setting N sufficiently high removes the nearby image sources, which generate fields that enter the region of interest of the simulation.Simulations with fields that propagate many periods away from the source from which they were generated typically require more integration points when using the ASM, since fields generated by distant periodic image sources can interfere with the fields of interest.
Typically, N cannot be determined a priori for a given simulation since it is dependent on the geometry and simulation time.The quadrature order is usually determined by progressively increasing N until the numerical integral has converged.
The abscissa points used in the midpoint and trapezoidal rule algorithms are necessarily equally-spaced to attain high convergence rates [28].As a consequence, no evaluations of f (k x ) may be re-used when refining the number of quadrature points from N to N + 1 points using the midpoint rule.The trapezoid rule behaves only marginally better, allowing the i = 0 term of the trapezoid rule sum (19) at N points to be re-used by the N + 1-point simulation.Since an evaluation of f (k x ) requires an independent simulation result, increasing N to N + 1 in order to observe convergence may be computationally intensive.
The process of increasing the number of integration points to observe convergence is made more efficient by finding a larger value by which to increase N .This allows us to reuse the fields from the N simulations that were previously run.In the case of the midpoint rule (18) we look for the smallest N M > N such that for all values of i there exists an integer 0 ≤ l < N M such that i + 1 2 Rearranging this formula, we find The quotient term must resolve to be the same constant for each i, so we take N M = 3N , which occurs as l = 3i + 1.
Proceeding similarly, we determine the smallest N T > N that enables the trapezoidal rule to reuse the N previouslycomputed fields Fig. 4. ASM can be used to convert a source S(x, y, t) from a single unit cell into one spanning multiple cells, A(x, y, t).This figure illustrates a component of a source A(x, y, t) spanning 7 × 7 unit cells in the x-y plane.Dotted green lines define the unit cells.
We take N t = 2N which occurs as l = 2i.The expediency of increasing N to 2N over increasing N to N +1 is considerable when using the trapezoidal rule; both scenarios require 2N total simulations, but the former produces far more accurate results.These calculations reveal a distinct advantage to using the trapezoidal rule over the midpoint rule.Although both quadrature techniques exhibit similar convergence behavior, the trapezoidal rule can reuse all the formerly-computed fields whenever the number of abscissas is doubled, whereas the midpoint rule requires that the number of abscissas is tripled.
Computation of the ASM integral can be further optimized in scenarios with even or odd symmetry, in which case N can be reduced by a factor of two [8].

IV. BEAMS IN PERIODIC GEOMETRIES USING THE ASM ALGORITHM
As described in Section II-A, the TF/SF technique can be used to introduce arbitrary beams into a FDTD domain, as long as the incident electric and magnetic fields are known a priori along the TF/SF boundary.
Beams generated by sources spanning multiple periodic unit cells can be simulated in a single period using the ASM method.As a simple and practical illustration of the principle, this section derives the TF/SF current sources in the x − y plane used to generate a z-directed beam.
Suppose that the geometry of the structure under consideration is periodic in the x and y directions, with periods d x and d y , respectively.Defining ρ = (x, y) T , consider an electric or magnetic source A(ρ, t) describing a beam at z = z 0 [see (9a)].Suppose the source A(ρ, t) spans B x = 2N x + 1 and B y = 2N y + 1 cells in the x and y directions, respectively, as in Fig. 4, so that the source covers the plane With d = n x d x x + n x d y ŷ, we define the source Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Fig. 5. Source A(x, y, t) (as in Fig. 4) from the unit cell (n x , n y ) = (2, 1) is simulated in the primitive periodic unit cell (0, 0).When the source is phase-shifted with exp(− j k • d), the ASM effectively moves the source to the desired unit cell (n x , n y ).Gray circles represent three rings of nulls of the source in Fig. 4. The expected correspondence between the gray circles and the nulls of the shifted source in cell (2, 1) is present. where and where is the region occupied by the central cell.The phase shift in ( 24) is applied so that the ASM effectively translates the source A(ρ + d, t) C (ρ) to the period (n x , n y ), thereby producing a source equal to A(ρ, t) C (ρ − d) = A(ρ, t) C+d (ρ) (as shown in Fig. 5).Superimposing all the sources S n x ,n y (ρ, t) produces the effect of the source A(ρ, t) The source S(ρ, t) is only defined in the primitive period (0, 0), but with the ASM, it creates fields that would be produced by the beam-generating source A(ρ, t) over its entire range F.
To see that ( 26) indeed defines a source producing the desired beam when the ASM is applied, let us define V n x ,n y (ρ, t) to be the electric or magnetic fields generated by S n x ,n y (ρ, t) interacting with the infinitely-periodic geometry.The field may be explicitly defined using the dyadic Green's function of the periodic structure G(ρ, ρ ′ ) Following (26), the function describes the field produced by S(ρ, t) interacting with the periodic geometry.The PBCs infinitely replicates the source S as a phased array, to produce fields equal to where d ′ = n ′ x d x x + n ′ y d y ŷ.These are the fields realized by a single FDTD simulation of the source S with PBCs operating at the wave vector k.
We apply the integral from (17) in two dimensions (31) to determine the fields obtained by the ASM in period (N x , N y ) the integral in (32) removes all the terms of V ∞ , except those satisfying d − d ′ − D = 0.This results in the evaluation The summand may be simplified The last equality follows from the translational invariance of Green's function.Using (34), we conclude Equation (38) implies that the fields recovered by the ASM are those in period D generated by the source A defined over the plane F.
As described in Section III-A1, the ASM integral in (32) can be effectively approximated using the trapezoid rule.The effect of the trapezoid rule on calculating the ASM integral can be assessed by observing that Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
From this formula, we see that after numerical integration, an infinite array of image sources are present, each producing a translated copy of V ASM .As the number of abscissa points N and M are increased, the image sources are distanced further from the desired source (corresponding to l = m = 0).The numerical integral effectively converges when the image sources are sufficiently far that their effects on the fields of interest are negligible.

V. GAUSSIAN BEAM SIMULATIONS
This section illustrates the beam-generating method introduced in Section IV.The first example shows how the ASM can be used to simulate an angled, broadband Gaussian beam propagating through a region of free space spanning many periods.The second example shows how the ASM can be used to simulate an obliquely angled Gaussian beam impinging on an infinitely-periodic metallic grating.
In both sections, full-sized simulations of the geometries, without using PBCs, were run to confirm the correctness of the ASM results.The full-sized simulations contained several periods and were terminated along all boundaries with halfwavelength perfectly-matched layers (PMLs).The difference between the fields determined using the ASM E ASM x and the fields from the full-sized simulations E full x is calculated using the relative root-mean-square (rms) error formula In this formula, the set TF denotes the TF region of the simulation, excluding PMLs.The evaluation of the ASM integral (17) requires several simulations of the same period, each with different PBC wavenumber components.All these simulations can be run in parallel, in which case the total simulation time would approximately equal the execution time of the slowest of the single-period simulations.In order to quantify execution times of the ASM simulations in Section V-A and V-B, each ASM simulation was run three times.In each of the simulations, the execution time of the slowest of the single-period simulations was recorded.The average of the slowest single-period execution times of the three ASM simulations is given in each section below.Each section also reports the average simulation time between three identical runs of the full-sized simulations, for comparison.
All simulations were performed using MATLAB on a computer with a 1.90 GHz Intel Core i7-8665U CPU, and with 16 GB of RAM.Simulations were run once before recording execution times to accommodate just-in-time (JIT) compiler optimizations.

A. Gaussian Beam Simulation: Free Space
A simulation of a transverse electric (TE) Gaussian beam traveling in free space at an angle of θ = π/16, ϕ = 0 was carried out using the technique described in Section IV.The source was given a broadband profile using a Blackman-Harris window with a central frequency of f c = 1 GHz and with a HWHM of f bw = 0.07 f c .The Gaussian beam parameters were set to z 0 = 2.5 m and A 0 = 1.
The simulation volume was discretized into square Yee cells with side lengths = λ/20, where λ = c/ f c is the wavelength corresponding to the central frequency of the beam fields.The time step was fixed at t = 0.99 /( √ 3c), corresponding to a Courant reduction factor of 0.99 [21].At this time step, we use (11) to determine that the Gaussian source is turned off after T b / t = 123 time steps.
The dimensions of the simulated unit period were (d x , d y , d z ) = (0.30, 0.30, 1.50 m) in the Cartesian axes x, y and z, with an additional half-wavelength PML terminating the faces normal to the z-axis.Periodicity was taken to be in the x and y directions.
For reference, define the origin (0, 0, 0) to be the center of the computational volume.The Gaussian beam sources were defined in the plane three Yee cells from one of the PMLs at z = −70.4cm.
The 2-D ASM integral (17) was approximated by the trapezoidal rule of order 11 × 11 in the x and y directions and was run for t v = 150 time steps.
For comparison, an equivalent full-sized simulation of 11 × 11 periods terminated by PMLs was also produced  and run for t v time steps.The Gaussian beam in the full-sized simulation had the same parameters and placement as the beam generated in the ASM simulation.
A snapshot of the E x field across 11 periods in the y = 0 plane at t v = 150 time steps using the ASM (E ASM x ) is shown in Fig. 6(a).The E x field from the full-sized simulation of 11 × 11 periods of free space (E full x ) is shown in Fig. 6(b).The field components E ASM x and E full x at t v = 150 time steps are compared along the cut (x, 0, 0.2 m) in Fig. 7(a), and along the cut (0, 0, z) in Fig. 7(b).
The relative error between the ASM simulation and full-sized simulation of the free-space Gaussian beam is ε free space = 0.00129, indicating strong agreement.The average Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.execution time of the ASM simulation was 6.19 s when it was run at full parallelism three times, and the simulation of each period consumed 355 MB of memory.Three simulations of the full-sized geometry were completed in 461 s, on average, and used 3.05 GB of memory.

B. Gaussian Beam Simulation: Metallic Grating
In this section, we present a demonstration of the interaction between a Gaussian beam and an infinite periodic metallic grating by simulating a single periodic unit cell using the ASM, and compare it with a full-sized simulation of the same structure.We added a square metallic grating element centered at (0, 0, −0.33 m) to the unit cell described in Section V-A.The metallic square was one Yee cell thick in the z direction and was λ/2 wide in the x and y directions.The simulated period corresponded to the period of the grating.Fig. 8 shows a schematic of the unit period simulated using PBCs with the ASM method; Fig. 9 shows a rendering of part of the simulation geometry in an equivalent full-sized simulation.The same 1 GHz angled Gaussian beam described in Section V-A was generated to impinge on the grating.
An ASM simulation of the structure was carried out by integrating fields at 11 × 11 wave numbers using the trapezoidal rule.Eleven periods of the grating in the y = 0 plane are shown in Fig. 10(a) with E x fields (E ASM x ) at t g = 200 time steps.A full-sized simulation of 11 × 11 periods of the grating was also performed without PBCs, with E x fields (E full x ) shown in the y = 0 plane in Fig. 10(b).In both figures, the Gaussian beam can be seen to diffract at the grating interface, producing a second beam traveling at θ = −53.6 • .Fig. 11(a) and (b) shows 1-D comparisons of fields from the ASM simulation and from the full-sized simulation, sampled along (x, 0, 0.2 m) and along (0, 0, z), respectively.The strong agreement between the fields from ASM simulation and the full-sized simulation is reflected quantitatively by the relative error (42) ε grating = 0.00768.
The ASM and full-sized simulations were each run three times.At full parallelism, the ASM simulations ran for an average of 6.62 s, with a peak memory usage of 333 MB by each simulated period.The full-sized simulations took an average of 558 s to complete, and consumed 3.03 GB of memory.

VI. CONCLUSION
Simulating physical beams interacting with periodic structures has been problematic in the FDTD algorithm.Periodic boundary conditions are typically used to simulate infinite periodic structures, but they also duplicate the sources in each period.To overcome this problem, a method for efficiently simulating generic electromagnetic beams interacting with periodic structures using the ASM was presented.
Electromagnetic beams can be produced using electric and magnetic currents, via the source-equivalence principle.Consequently, the ASM-a technique to simulate finite sources in infinitely periodic media-can be applied.The ASM works by numerically integrating fields across numerous parallelizable simulations of a single periodic unit cell.
The ASM integral is effectively estimated using midpoint and trapezoidal quadrature.Trapezoidal quadrature has the characteristic that whenever the number of abscissa points is doubled, all the previous simulation results can be reused.In this way, it offers an efficient means to refine the ASM results if they are not sufficiently accurate.The midpoint rule, on the other hand, is able to save the original simulation results only when the number of abscissa points is tripled.
In general, electromagnetic beam sources span multiple periodic unit cells.The ASM is able to model sources which extend beyond a single period by using the fact that a source in the simulated unit cell may effectively be moved to any arbitrary period by phase-shifting the source appropriately.The ASM can therefore unfold a superposition of phase-shifted sources to produce the effect of a continuous source spanning numerous periods.
Two simulations were presented to illustrate beam propagation in periodic environments.The first simulation modeled an obliquely-angled broadband Gaussian beam traveling through free space generated by a source spanning 11 × 11 unit cells.A similar angled Gaussian beam was simulated impinging on an infinite, 2-D metallic grating.Both simulations were compared against full-sized simulations without PBCs.Strong agreement between the results was exhibited, confirming the efficacy of the algorithm described in this article.Run time savings of the ASM simulations, in comparison to the equivalent full-sized simulations, were additionally exhibited.

Fig. 1 .
Fig.1.Illustration of the source equivalence principle.A green point source is shown in (a), which generates the fields depicted in red.A fictitious boundary ∂ is drawn (dotted gray), separating an internal region ¯ from an external region .Figure (b) shows how the fields produced by the point source can be generated throughout , while maintaining an absence of fields in ¯ .Electric and magnetic currents along ∂ are defined using tangential field components of the fields via (1a) and (1b).A FDTD discretization of the current densities J s and M s leads to the TF/SF method.

Fig. 3 .
Fig. 3. Unit cell showing 2-D FDTD PBC updates (in the z direction).The E x field at z = d z updates the E x node at z = 0 (red arrow).The H x field at z = z/2 is used to update the corresponding magnetic field at z = d z + z/2 (blue arrow).
When (32) is evaluated by a 2-D N × M-point trapezoid rule, all the terms of the sum are removed except those satisfying d−d ′ − D = D Trap l,m where D Trap l,m = Nld x x+Mmd y ŷ.Consequently, the fields obtained by the ASM using the trapezoid rule are equal to

Fig. 6 .
Fig. 6. y = 0 plane illustrating the E x field intensity of the angled Gaussian beam (as described in Section V-A) across 11 periods at t v = 150 time steps, generated (a) using the ASM and (b) by direct simulation without PBCs.The dotted green vertical lines on the left denote the current sources.The horizontal dashed green lines in (a) mark the boundaries of the periodic unit cells.The vertical and horizontal dashed purple lines indicate the paths along which fields are sampled in Fig. 7(a) and (b), respectively.PBCs are omitted from figures.

Fig. 7 .
Fig. 7. E x field produced by the Gaussian beam traveling through free space at t v = 150 time steps, sampled along (a) (x, 0, 0.2 m), and along (b) (0, 0, z) by the ASM technique and by a full-sized 11 × 11-period simulation.

Fig. 8 .
Fig. 8. Simulated periodic unit cell in the x-z plane.The blue regions on either side represent PMLs.The dotted vertical red line represents the TF/SF source used to produce the Gaussian beam; the TF region lies to the right of the line, and the SF region lies to the left.A single grating element exists in the center of the simulation.

Fig. 9 .
Fig. 9. Three-dimensional schematic illustration of 3 × 3 periods of a full-sized, finite simulation of a 2-D grating.The blue regions represent PMLs.The red surface represents the TF/SF source used to produce the Gaussian beam.

Fig. 10 .
Fig. 10.y = 0 plane illustrating the E x field intensity of an angled Gaussian beam impinging on a periodic grating (as described in Section V-B) across 11 periods at t g = 200 time steps generated (a) using the ASM and (b) by direct simulation without PBCs.The dotted green vertical lines on the left denote the current sources.The horizontal dashed green lines in (a) mark the boundaries of the periodic unit cells.The vertical and horizontal dashed purple lines indicate the paths along which fields are sampled in Fig. 11(a) and (b), respectively.PBCs are omitted from figures.In both figures, a diffracted wave can be seen traveling at θ = −53.6 • due to the effect of the grating.

Fig. 11 .
Fig. 11.E x field produced by the Gaussian beam impinging on a grating at t g = 200 time steps, sampled along (a) (x, 0, 0.2 m), and along (b) (0, 0, z) by the ASM technique and by a full-sized 11 × 11-period simulation.