An Expedient Approach to FDTD-based Modeling of Finite Periodic Structures
aa r X i v : . [ c s . C E ] A ug An Expedient Approach to FDTD-based Modelingof Finite Periodic Structures
Aaron J. Kogon, and Costas D. Sarris,
Senior Member, IEEE
Abstract —This paper proposes an efficient FDTD techniquefor determining electromagnetic fields interacting with a finite-sized 2D and 3D periodic structures. The technique combinesperiodic boundary conditions—modelling fields away from theedges of the structure—with independent simulations of fieldsnear the edges of the structure. It is shown that this algorithmefficiently determines the size of a periodic structure necessaryfor fields to converge to the infinitely-periodic case. Numericalvalidations of the technique illustrate the savings concomitantwith the algorithm.
I. I
NTRODUCTION
Periodic structures are common geometries in electromag-netics, appearing in such forms as frequency selective surfaces,electromagnetic band gap media, photonic crystals, antennaarrays and metasurfaces, among others.A common way to simulate infinitely periodic structuresin the Finite-Difference Time-Domain (FDTD) method [1] isby employing periodic boundary conditions (PBCs) [2]. PBCscollapse an infinitely-long periodic domain into a single unitcell. PBC simulations have been successfully used to extracttransmission and reflection spectra [3], attenuation constants[4], [5], Brillouin diagrams [4], [6], among others.Although PBCs are effective at simulating infinitely-periodic structures, the assumption of infinite extent is non-physical. Fields obtained by simulations with PBCs are thusnecessarily approximations of real finite-sized periodic struc-tures. The approximations become increasingly worse closerto the edges of the structure, where the periodicity ends. Manyanalytical, semi-analytical and method-of-moments modelsfor analyzing truncated structures have been discussed inthe literature [7]–[12]. In general, full-scale simulations offinite periodic structures—without PBCs—are typically run toprecisely determine edge effects [13]. Full-sized simulationscan be computationally intensive.This paper advances a technique for accurately simulatingfields interacting with finite-sized periodic structures. This hy-brid method takes advantage of the periodicity of the structurewhile also effectively capturing the edge effects.First, the mathematical and theoretical background on pe-riodic boundary conditions in FDTD is outlined. Next, thearray scanning method (ASM), an algorithm used to removeunwanted periodic sources from PBC simulations, is pre-sented. The methodologies for efficiently simulating finitetwo- and three-dimensional periodic structures are describedand numerical results are provided. Convergence quantificationmetrics are defined to numerically measure accuracy of resultsobtained through the algorithm. Because of their computational efficiency, PBCs are heav-ily employed in the design of periodic structures. However,PBCs simulate infinite, practically unrealizable, structures. Apertinent question is how many unit cells of a finite periodicstructure are needed to approximate the behavior of an infiniteone. This question is addressed in the final section of thispaper.II. P
ERIODIC B OUNDARY C ONDITIONS (PBC S )Before discussing the proposed algorithm, a mathematicalbackground for FDTD PBCs is provided, with a discussion onthe array scanning method (ASM) [2].Let U ( r ) represent an electromagnetic field in a three-dimensional periodic structure (periods d x , d y and d z in the x, y and z directions respectively). Floquet’s theorem statesthat for phasor e U : e U ( r + d ) = e U ( r ) e − j k · d (1)where d = n x d x ˆ x + n y d y ˆ y + n z d z ˆ z is the lattice vector of thestructure ( n x , n y , n z ∈ Z ) and where k = k x ˆ x + k y ˆ y + k z ˆ z isthe reciprocal (Bloch) wave vector. Define the inverse Fouriertransform of e U as U s ( r , t ) = 12 π Z ∞−∞ e U ( r , ω ) e jωt dω (2)so that the the real and imaginary parts of U s correspondto physical time-domain fields in quadrature. The field U s additionally satisfies the Floquet condition (equation (1)).Both the real and imaginary components of U s may besimultaneously simulated using the FDTD algorithm. Then,the complex-valued Floquet condition, which relates fields oneperiod apart, may be used as a periodic boundary conditionover a unit cell. The wave vector k of the real and imaginaryparts of U s is set by the PBCs.To illustrate the implementation of FDTD PBCs, consider astructure that is periodic in the z direction with a unit cell d z ,as shown in Fig. 1. Introducing auxiliary H y nodes one half-Yee cell outside the periodic domain at z = d z + ∆ z/ , thePBC update equations may be written as E x ( x,
0) = E x ( x, d z ) e jk z d z (3) H y ( x, d z + ∆ z/
2) = H y ( x, ∆ z/ e − jk z d z (4)FDTD PBCs can impose periodicity in more dimensionssimilarly [2].As mentioned, PBCs can be used to efficiently determineuseful parameters and properties of periodic structures. Theycan be used to produce reflection and transmission spectra k z k θz = 0 z = d z E x ⊙ ⊙ H y ∆ zE x ⊙ H y E z Fig. 1. A period showing two-dimensional FDTD PBC updates (in the z direction). The E x field at z = d z is used to update the correspondingelectric field at z = 0 (red arrow). The H y field at z = ∆ z/ updates the H y node at z = d z + ∆ z/ (blue arrow). S ou r ce I m a g e P BC P BC Translate & phase shift
Fig. 2. An illustration of a unit cell (right square) with PBCs (blue lines)on the left and right faces. A source in the unit cell produces fields (red)which travel rightwards within the dashed lines emanating from the source.The fields are then phase-shifted and translated to the left face, where theycontinue to propagate, as shown. These fields then appear as though theywere generated by a phase-shifted image source in an adjacent unit cell (leftsquare). diagrams, Brillouin diagrams and determine attenuation con-stants. Of particular importance to this paper, however, is thatPBCs can be used to determine fields due to finite sourcesinteracting with periodic structures using the array scanningmethod.
A. Array Scanning Method (ASM)
PBCs operate by translating fields from one edge of theperiodic domain to the other side, with a complex phase shift.When fields produced by a source arrive at a PBC, they appearat the opposite side of the simulation with a phase shift (seeFig. 2). This produces the effect of having infinite periodswith sources in each period. The image sources have a linearlyprogressive phase shift from one unit cell to the next, set bythe wave vector k used in the PBC formulation. Explicitly, asimulation of a structure which is periodic in the z directiongenerates the field U ∞ such that: U ∞ ( x, y, t, k y ) = ∞ X n = −∞ U n ( x, y, t ) e − jk y nd y (5)where U n ( x, y, t ) = U ( x, y − nd y , t ) represents the fieldsgenerated by a source in unit cell n . · · · S ou r ce I m a g e ad Md − a Fig. 3. An illustration of the ASM method removing unwanted image sourcesin a periodic structure (period d ). The ASM integration order M determinesthe location of the nearest parasitic image source. By selecting M to besufficiently large, fields (shown in red) generated by the image do not enter aregion of interest a units to either side of the source within runtime. By theend of the simulation, the wavefronts from unwanted sources travel a distanceof at most Md − a . The orthogonality of the complex exponential may be usedto extract fields due to sources in unit cell n : U n ( x, y, t ) = d y π Z π/d y − π/d y U ∞ ( x, y, t, k y ) e jk y nd y dk y . (6)In three-dimensional periodic structure (size d x × d y × d z ), thefield U l,m,n due to a source in the ( l, m, n ) th unit cell is: U l,m,n ( r , t ) = d x d y d z (2 π ) Z π/d x − π/d x Z π/d y − π/d y Z π/d z − π/d z U ∞ ( r , t, k ) e j k · d d k (7)where d = ( ld x , md y , nd z ) .This integration procedure is known as the array scanningmethod [14]–[16], which provides a technique for removingunwanted source images.Equation (6) may be approximated by numerical integrationof fields determined by FDTD. U ∞ can be evaluated from asingle simulation, where k y is set by the PBCs. The varioussimulations used to approximate equation (6) are disjoint andmay be run simultaneously.One of the most robust quadrature techniques for estimatingequation (6) is the midpoint rectangular rule: U Mn = 1 M M − X m =0 I (cid:18) k y = ∆ k y (cid:18) m + 12 (cid:19) − M (cid:19) (8)where I = U ∞ ( x, y, t, k y ) e jk y nd y is the integrand of equa-tion (6) and ∆ k y = 2 π/ ( M d y ) . This formula removes theeffect of image sources in M unit cells to either side ofunit cell n exactly, up to a numerical error. The midpointrectangular rule has been additionally found to outperformhigher order integration methods [16].Using the midpoint rectangular rule, it is straightforwardto select quadrature order M to be large enough that fieldsfrom parasitic images do not enter a region of interest withina time t . Suppose fields in a region a units to either side of asource are desired in a periodic structure of period d (shownin Fig. 3). If M is sufficiently large, the wavefront from aparasitic image source will reach the boundary of the regionof interest at the end of the simulation, meaning ct = M d − a where c is the speed of light. Selecting the integration orderto be at least M = ⌈ ( t c + a ) /d ⌉ ensures that fields fromunwanted image sources will not enter the region of interestby t . Symmetry with respect to the periodic directions (as inFig. 6) may be exploited to reduce the ASM integration order.When the field generated by sources in the center period U is even: U n ( x, y, t ) = d y π Z π/d y (cid:0) U ∞ ( x, y, t, k y ) e jk y nd y + U ∞ ( x, − y, t, k y ) e − jk y nd y (cid:1) dk y . (9)When U is odd: U n ( x, y, t ) = d y π Z π/d y (cid:0) U ∞ ( x, y, t, k y ) e jk y nd y − U ∞ ( x, − y, t, k y ) e − jk y nd y (cid:1) dk y . (10)The domain of integration of these integrals is half of thatin equation (6). Therefore, the computation of these integralsrequires half as many simulations to be approximated at agiven precision.III. G ENERAL M ETHODOLOGY
Simulations with PBCs produce fields which behave asthough they interact with an infinitely periodic structure.Therefore, fields obtained through simulations with PBCs maybe reasonably used as approximations of fields interacting withfinite periodic structures far from the edges of the structurewhere the periodicity ends.Fields determined by PBC simulations may then be injectedinto a full-size FDTD simulation of the edge unit cells usingthe Total-Field Scattered-Field (TF/SF) source condition [1].With the TF/SF boundary “replaying” the fields from thePBC simulation, only the edges of the structure need to besimulated to sufficiently capture the edge effects. When thereis a sufficient number of periods included in the simulationof the edges of the structure, there is a smooth transitionbetween fields obtained by the ASM and fields obtained fromthe simulation of the edges. At this point, fields obtained by theASM and fields obtained from the simulation of the edges maybe patched together to obtain a continuous field. Fig. 4 showshow the infinite structure provides an estimate for the fields inthe inner region (far from the edges), and a simulation of theedges alone can estimate the edge fields. The fields of the innerand edge regions are pasted together at the interfaces betweenthe regions, as shown. The procedure is outlined schematicallyin Fig. 5.The TF/SF boundary reproduces the fields propagatingalong an infinitely periodic structure. If there were an infinitenumber of edge unit cells, there would be no fields in thescattered-field (SF) region. Since the simulation of the edgesis finite, however, some fields will not behave precisely asthe TF/SF boundary anticipates, so some fields will appearbehind the boundary. The fields in the SF region become lesspronounced as more periods along the edge are simulated.In any case, an absorber should be placed behind the TF/SFboundary in the SF region to absorb these fields and preventthem from being reflected back.Some time t s is required for wavefronts to reach the planesalong which sampling takes place in the ASM simulations.This means that time t s is required to elapse before the Infinite structure ...... EdgesMerged fieldsEdgeregion Innerregion Edgeregion
Fig. 4. The fields in the inner region (between the vertical blue lines) of afinite structure can be estimated by fields from an infinite structure (top). Thefields along the edges of the structure (outside the blue lines) can be calculatedby simulating the edges of the finite structure (middle). The fields can be fusedtogether to create a continuous field approximating those interacting with thefinite structure (bottom).
Set number of internal unit cells N I and number of edge unit cells N E Use ASM to simulate fields of interest in in theinternal N I unit cells from the first unit cellSimulate edges of the structure with TF/SFsources, using the tangential fields from theASM simulation; extract fields of interestMerge fields of interest from the ASMsimulations and the simulation of the edges Fig. 5. A flowchart describing the procedure of efficiently determining fieldsinteracting with finite periodic structures and estimating the accuracy of theresult. Note that the ASM step may be used to determine tangential fieldsalong the boundaries of many differently-sized regions simultaneously. Sub-sequently, variously-sized edge-unit cell simulations can be run simultaneouslyas well.
TF/SF boundaries start emitting any fields at all. For instance,consider a source positioned in a periodic structure as Fig. 6.If l is the shortest distance between the source and a one ofthe sides of the inner region (where sampling takes place),waves arrive arrive at the nearest inner region boundary attime t s = l/c or later. Thus, there is no need to record (orreplay) fields before time t s .It should be noted that only one set of ASM simulationsis necessary to extract tangential fields over internal regionsof several sizes. This is possible because the ASM may beused to extract fields anywhere along the infinitely-periodicstructure, given that the ASM integration order is sufficientlyhigh. Additionally, several simultaneous simulations of theedges may be carried out by parallelization.In a two-dimensional simulation, a periodic surface has twodisconnected edges (see Fig. 6), in which case each edgemay be simulated independently. A three-dimensional periodicstructure has edges continuously running along the edges ofa rectangle, shown in Fig. 11, which inhibits the separationof the edges from each other. The two- and three-dimensionalcases are discussed individually below.In sum, the proposed algorithm may be broken into a two-step process. First, the ASM is used to record fields tangentialto the periodic faces (several periods away from the source)as a function of time. Second, a simulation of the edgesis performed, where the fields previously determined by theASM are injected into the simulation using the TF/SF sourceformulation.All simulations below were performed using MATLAB ona computer with a quad-core .
00 GHz
Intel Core i5-7400processor and of RAM.IV. 2D F
INITE -S IZED S TRUCTURE S IMULATIONS
Two-dimensional simulations may be used to model struc-tures and sources with infinite length in one axis. In manysituations, two-dimensional simulations can be used to effi-ciently approximate physical structures.Consider a structure made up of N E edge unit cells at eitherside of the structure and N I + 1 internal unit cells as shownin Fig. 6. First, a single unit cell of the structure is simulatedwith PBCs (as in Fig. 7a). The ASM is used to determineelectric and magnetic fields perpendicular to the direction ofperiodicity, N I periodic unit cells to either side of the first unitcell as functions of time. The position on the computationalunit cell is shown as a red dotted line in Fig. 7a. That is, theASM may be used to record the fields E z ( x, ( N I + 1 / d y , t ) , H x ( x, ( N I + 1 / d y , t ) (11)and E z ( x, − ( N I + 1 / d y , t ) , H x ( x, − ( N I + 1 / d y , t ) . (12)The recorded fields can then be “replayed” in a simulationof a finite number of unit cells to model the edge of periodicstructure (as in Fig. 7b). Some space in the SF region may beused as buffer between the source and the PML. Source N I + 1 Periods N E Periods N E Periods
Fig. 6. A periodic grating in two dimensions. The inner region of N I + 1 unit cells may be approximated by the ASM, and the outer edge regions of N E unit cells may be simulated separately. PBCPBC P M L P M L G r a ti ng Sample S ou r ce H y ⊙ E z H x (a) PMLPML P M L P M L TF/SF Boundary G r a ti ng N E + . TF RegionSF Region H y ⊙ E z H x (b)Fig. 7. Schematic of (a) a single unit cell simulation and (b) of the edge effectssimulation of the subwavelength grating. The red dotted line in (a) displayswhere electric and magnetic fields are sampled. The red dotted line in (b)represents the TF/SF boundary at which the recorded fields are injected intothe simulation. The “total-field” (TF) and “scattered-field” (SF) regions aredenoted. Ten periods of the grating appear in the SF region, while N E + 0 . periods of the grating are in the TF region. A. Finite-Sized 2D Subwavelength Grating
In order to illustrate the procedure, a 2D transverse-electric(TE) simulation of a 200-unit cell subwavelength metallicgrating was performed. Subwavelength gratings have beenstudied for their ability to convert evanescent waves intopropagating waves, making them useful in lensing applications[17], [18]. They have also been used as reflectors [19] andoptical switches [20].First, a simulation of a single unit cell of the grating wascarried out, with PBCs (see Fig. 7a). The operating frequencyof the structure was set to f op = 12 GHz . The simulation wasrun for steps.The dimensions were set to d x = 4 .
37 cm and d y =1 .
25 mm in the x and y directions respectively. The perfectelectric conductor (PEC) making up the metal of the grating spanned half of the unit cell and was one Yee cell thick.The domain was discretized into 361 and 889 Yee cells inthe x and y directions respectively ( ∆ x = ∆ y = λ/ ) inorder to finely model the subwavelength features of the fieldsand grating. A Courant number of . was chosen to calculatethe time step.The boundaries parallel to the x axis were lined by 40-cell-thick perfectly-matched layers (PMLs) [1]. The boundariesparallel to the y axis were bordered by PBCs.The source excitation was defined as a Gaussian electriccurrent J z containing the operating frequency. The source waspositioned λ/ units to the side of the grating. The ASM wasused to record fields along one boundary parallel to x (denotedby a red dashed line on the simulated unit cell in Fig. 7a) asa function of time.Symmetry was exploited to improve ASM efficiency (seesection II-A). Fields from this simulation were used to ap-proximate the internal region of the grating, away from theedges.Another simulation was then carried out to approximate thefields by the edges (see Fig. 7b). The discretization of spaceand time of this simulation matches that of the PBC simulationdescribed above, as does the dimension parallel to the x axis,the positioning of the grating and the PML parameters.The simulation of the edge had PMLs along each boundary.Ten unit cells of the grating were placed between the TF/SFboundary and the PML, in the SF region. The number ofperiodic unit cells in the SF region does not appear to greatlychange the field patterns in the TF (total-field) region. Avariable number of unit cells of the grating were located inthe TF region. Empty space spanning was included inthe y direction, shown below the grating in Fig. 7b.Sampling was performed along lines positioned at λ and λ/ units away from the grating (on the opposite side of thesource) and Fourier transformed at the operating frequency.Fields obtained from a finite simulation, from the ASMsimulation (i.e. a simulation of an infinitely-periodic grating)and from the grating edge simulation are shown in Fig. 8. Thesimulations of the edges were carried out with 25 edge unitcells and without any edge unit cells. Fields to either side ofthe TF/SF boundary are included in these plots.Edge effects and the beginnings of side lobes are visible inboth plots. There are strong ripples around the edges Fig. 8a inparticular. Edge fields can be accurately estimated with veryfew, or even no edge unit cells. However, if one wants tosmoothly merge the fields obtained by the ASM and the fieldsobtained via the edge unit cells simulation to produce a singlecontinuous field, several unit cells may need to be included inthe edge unit cells simulation as the ASM may be insufficientto describe fields near the edges. In Fig. 8a, for example, theASM fields may be merged more readily with the 25-unit-cellsimulation of the edge than with the 0-unit-cell simulation. B. Computational Cost Analysis
The times taken to complete the full-sized simulations, theASM simulations and the edge unit cells simulations wererecorded. Due to symmetry, only half of the grating was . . . . Distance y from source (cm) E z A m p lit ud e Full ASM0 25 (a) . . . . Distance y from source (cm) E z A m p lit ud e Full ASM0 25 . . . . Distance y from source (cm) E z A m p lit ud e Full ASM (b)Fig. 8. Normalized magnitudes of frequency-domain E z field samples at (a) λ units from the grating and (b) λ/ units from the grating. Each plot showsfields from full-sized simulations of a 200-unit-cell grating (labelled in thelegend as “Full”), fields obtained from the ASM, and fields captured in thesimulation of edge unit cells. In the latter case, the legend entry is the numberof edge unit cells simulated ( N E ). A vertical dashed line marks the end ofthe grating. An inset shows a magnification of the region near the edge ofthe grating. The transitions of the fields of the edge-cell simulations from theSF region, where fields are low, to the TF regions, where fields track thoseof the full simulation, are distinguishable. simulated. A perfect magnetic conductor (PMC) was placedalong the symmetry axis to produce the desired field symmetryprofile.Table I describes the timing of the ASM and edge unit cellsimulations when N E = 0 and N E = 25 edge unit cellsare considered. The time taken to simulate the entire structureis included as well. The full simulation time assumes ASMparallelization. TABLE IT
IME TAKEN TO SIMULATE A
UNIT - CELL SUBWAVELENGTH GRATING USING THE PROPOSED ALGORITHM ( VARYING EDGE UNIT CELLS N E ) AND AFULL SIMULATION OF THE FINITE STRUCTURE . T HE ASM
INTEGRATION ORDER USED FOR ESTIMATING THE INTERNAL FIELDS AWAY FROM THE EDGESIS PROVIDED ALONG WITH THE RUNTIME OF EACH SIMULATION .Sim. Type N E ASM Integration Order ASM Sim. Time (s) Edge Sim. Time (s) Total Sim. Time (s)Edge .
52 26 . .
925 201 1 .
52 46 . . Full — — — —
PML a s a s Source y ⊙ z x Fig. 9. A schematic of the full-sized simulation of the three-dimensionalphotonic crystal from [21] ( x - y plane, z = 0 ). V. 3D F
INITE -S IZED S TRUCTURE S IMULATIONS
Consider a structure which is periodic in the x and z directions in a three-dimensional simulation, as in Fig. 9.The algorithm outlined below is used to estimate the fieldsthroughout the structure following the two-step process de-scribed previously (section III).First, a single unit cell is simulated with PBCs along the fourperiodic faces, and absorbers along the other two (aperiodic)faces (see Fig. 10). Let us define the inner region of thestructure, where edge effects are negligible, to be N Ix unit cellsto either side of the center unit cell in the x direction, and N Iz unit cells on either side of the center unit cell in the z direction(see Fig. 11; the center unit cell, with the source, appears inFig. 9). The ASM can then be used to record tangential fieldsalong the boundary of the inner region. That is, the ASM canbe used to record: U x ( x, y, ( N Iz + 1 / d z , t ) , U y ( x, y, ( N Iz + 1 / d z , t ) (13) U x ( x, y, − ( N Iz + 1 / d z , t ) , U y ( x, y, − ( N Iz + 1 / d z , t ) (14)where | x | ≤ ( N Ix + 1 / d x and U y (( N Ix + 1 / d x , y, z, t ) , U z (( N Ix + 1 / d x , y, z, t ) (15) U y ( − ( N Ix + 1 / d x , y, z, t ) , U z ( − ( N Ix + 1 / d x , y, z, t ) (16) PBCPBC P M L P M L Samples S ou r ce x ⊗ z y (a) PBCPBC P M L P M L Samples S ou r ce z ⊙ x y (b)Fig. 10. (a) The z = 0 plane and (b) the x = 0 plane of one unit cell of thethree-dimensional photonic crystal. The ASM is used to estimate the fields inthe inner region of the finite structure and to record the fields which are usedto estimate the fields along the edges of the structure. The faces along whichfields are sampled in the simulated unit cell are marked by dashed red lines. where | z | ≤ ( N Iz + 1 / d z . In these expressions, U i denotesboth electric and magnetic field components.The recorded fields are then used to estimate the fieldsnear the edges. As shown in Fig. 11, the recorded fields areinjected several unit cells away from the edges in the periodicdirections. An absorber should be placed behind the TF/SFboundary in the SF region to prevent reflections of fields in theSF region back into the TF region. As with the 2D simulation(section IV-A), some space should be left between the TF/SFboundary and the absorber as a buffer zone for the absorberand to facilitate accuracy quantification of the algorithm (seesection VI). The absorber in the SF region should not be astandard PML, since convex PMLs are not stable [22].In the region enclosed by the absorber, fields do not needto be simulated since the ASM simulation estimates the fieldswithin this region. No updates are required within this region. A. Numerical Results: Finite-Sized 3D Photonic Crystal
A 3D model of a periodic photonic crystal was simulatedto illustrate the procedure of edge-effects determination. Thephotonic crystal simulated is notable due to its exhibition ofnegative refraction at all angles of incidence, making it usefulin lensing applications [21].The photonic crystal was composed of a dielectric slab witha relative permittivity ε r = 12 . centered at the origin, with PMLTF/SF N Ix + 1 Unit Cells N I z + U n it C e ll s ABCNo FDTD updatesSF RegionTF Region z ⊗ y x Fig. 11. Schematic of the simulation of the edges of the three-dimensionalphotonic crystal ( x - z plane, y = − a s − r ). The photonic crystal is shownin yellow, and each unit cell is outlined in grey dotted lines. The simulationdomain comprises of (from inside to outside): (1) a block (grey) which is notupdated by FDTD, (2) an ABC (red dashed line) encasing the non-updatedblock, (3) a TF/SF boundary (blue dotted line) several unit cells from theABC, (4) the photonic crystal boundary and (5) the PML along the simulationboundary (blue). The “total-field” (TF) and “scattered-field” (SF) regions aredenoted. staggered rows of vacuum-cylinders (radius r = 33 . )in the z -direction (see Fig. 9). The period in the x and y directions was a s = 135 . . Since the infinite structurewas uniform in the z direction in the ASM simulation (seeFig. 10b), the period was chosen to be a s . The photonic crystalwas four periods thick in the y direction, while the size ofthe crystal may be varied in the other axes. Throughout thissection, the photonic crystal studied had unit cells in boththe x and z directions.The slab was excited by a Gaussian magnetic current pointsource M z located at a distance r away from the slab face(i.e. at y = − a s − r ) in the x = z = 0 plane. Theexcitation spanned frequencies between f min = 0 PHz and f max = 1 . in order to enclose the operating frequency f op = 600 THz (cyan visible light).The domain was discretized as ∆ x = ∆ y = 8 .
58 nm in the x and y directions. The z axis, along which thephotonic crystal is uniform, was discretized in intervals of ∆ z = 22 . .A single unit cell was simulated using the ASM to calculatethe fields in the internal region (see Fig. 10). In the simulationof the edge effects, the inner region (in which no FDTDupdates occur) was surrounded by first-order Mur’s absorbingboundary conditions (ABC) [1]. The border of the ABC blockwas located one period a s from the TF/SF sources in the x and z directions in all simulations, as illustrated in Fig. 11.The number of periodic unit cells in the SF region betweenthe ABC and the TF/SF boundary does not appear to greatly . . . . . . . . Distance x from source ( µ m) H z A m p lit ud e Full ASM0 12 3
Fig. 12. Normalized magnitudes of Fourier-transformed H z samples ata distance r from the photonic crystal. The plot shows fields from full-sized simulations of the finite structures (labelled in the legend as “Full”),fields obtained from the ASM, and fields captured in the simulation of edgeunit cells. In the latter case, the legend entry is the number of edge unitcells simulated ( N Ex = N Ez ). The transitions of the fields of the edge-cellsimulations from the SF region, where fields are low, to the TF regions, wherefields are higher, are clearly distinguishable in the N I = 2 and N I = 3 cases. change the field patterns in the TF region.The simulation region was terminated with PMLs (thickness . ), located . from the edges of the photoniccrystal.The Courant number was set to . and the simulationwas run for time steps. A sample of H z at ( x, a s +36 . , —lying along the image plane of the source [21]—was measured and Fourier transformed at the operating fre-quency f op . The sample magnitudes at various numbers ofedge unit cells is shown in Fig. 12, and are compared againstthe full-sized simulation of the finite photonic crystal. Fieldsto either side of the TF/SF boundary are shown in these plots. H z fields from full finite simulations and from the algorithmdescribed here are presented throughout the entire z = 0 planein Fig. 13.Good correspondence exists between the fields from thefull-sized simulations versus fields combined from the ASMsimulations and from the simulation of the edge unit cells,when two or more edge unit cells are simulated. In particular,the simulations of the edges are able to approximate fieldsoutside the slab. Simulations with PBCs alone cannot do thatmeaningfully. B. Computational Cost Analysis
The time taken to complete the full-sized simulations, theASM simulations and the edge unit cells simulations of an × photonic crystal are recorded in table II. The ASMintegration order in the periodic directions ( x and z ) and thetime taken to simulate each unit cell are noted in the table. TheASM integration orders are not the same in each dimension, − . − . . . − − . . y ( µ m) x ( µ m ) − − − − − H z (dB) (a) − . − . . . − − . . y ( µ m) x ( µ m ) − − − − − H z (dB) (b)Fig. 13. The Fourier transform of the normalized H z field amplitudes at f op on the z = 0 plane determined from (a) a full-sized simulation and (b) frommerging fields from the ASM simulation with fields from the simulation ofthe edges ( N Ex = N Ez = 3 edge unit cells). White dashed lines in (b) outlinethe boundary between the internal region and edge regions. due to discretization discrepancies. As expected, significanttime savings are possible when the ASM simulations areparallelized.VI. Q UANTIFYING C ONVERGENCE OF F IELDS G ENERATED BY THE P ROPOSED A LGORITHM VERSUS F IELDS FROM F ULL F INITE S IMULATIONS
As seen in the previous sections, the merger of the ASMfields and fields from simulations of the edge unit cells canapproximate fields from the full simulations of the finitestructures. The accuracy of the approximation may vary asthe number of unit edge unit cells N E changes. A metric fornumerically quantifying this accuracy is therefore desired.A criterion for quantifying convergence may be identifiedfrom the previous results. A smooth transition between thefields of the ASM solution and the fields from the edgesimulation is required for accurately merging fields. Fig. 8a may be used to illustrate this criterion. There is nosmooth transition between the ASM solution and the N E = 0 case. When N E = 25 , there is a smoother transition betweenthe edge simulation fields and the ASM fields.To quantify the convergence criterion, the Euclidean ( L )norm is used to describe the difference between the fieldsdetermined by the ASM simulations and the fields deducedfrom the edge simulations over several unit cells.Consider a 2D simulation. Define E ASM ( y ) = E ASM z ( λ, y, f op ) and E N E ( y ) = E N E z ( λ, y, f op ) where E ASM z ( x, y, f ) and E N E z ( x, y, f ) are the z -directed electricfields (in frequency-domain) determined from the ASMsimulation and from simulations of N E edges, respectively.Then, the difference between E ASM and E N E over P TF unitcells in the total-field (TF) region past the TF/SF boundary(located at y = (0 . N I ) d y ) may be expressed using theformula: || E ASM − E N E || TF = sZ ( P TF +0 . N I ) d y (0 . N I ) d y | E ASM − E N E | dy. (17)The 3D formula may be derived in the same fashion.The variable P TF should be chosen to span the lengthover which agreement between the edge-cell simulation andASM is desired. If P TF is too small, || E ASM − E N E || TF maysuggest convergence in certain cases even if it has not yetbeen reached. If E ASM and E N E happen to be similar arounda small domain past the TF/SF boundary in the TF region,the metric will suggest convergence, even if it has not beenreached.The convergence metric was applied to the fields of thethe subwavelength grating sampled λ units away from thegrating (section IV-A) and of the fields of the photonic crystal(section V-A) described above. The results are presented intables III and IV. In the results, cases with the largest numberof edge unit cells have the lowest scores, implying strongestconvergence. As anticipated, || E ASM − E N E || TF becomes moreconservative as P TF is increased. More precisely, the metric || E ASM − E N E || TF (as a function of N E ) tends to zero moreslowly when when P TF is increased.VII. C ONVERGENCE OF F IELDS ON F INITE P ERIODIC S TRUCTURES
Finite periodic structures approximate infinitely periodicstructures increasingly well as the number of unit cells in-creases. Typically, the number of unit cells needed for conver-gence is determined by monitoring simulations of increasinglylarge finite structures [23]–[25]. The present algorithm maybe used to determine the number of unit cells required forconvergence more efficiently. By selecting the inner regionof unit cells N I + 1 to be the region where convergence isdesired, the number of edge unit cells N E can be increasedfrom zero until convergence is reached.Suppose convergence is desired within a range of N I + 1 unit cells in the periodic direction (as Fig. 6). In other words, itis desired that the fields within the region of N I + 1 internal TABLE IIT
IME TAKEN TO SIMULATE AN × PHOTONIC CRYSTAL USING THE PROPOSED ALGORITHM ( VARYING EDGE UNIT CELLS N E ) AND A FULLSIMULATION OF THE FINITE STRUCTURE . T HE ASM
INTEGRATION ORDER USED FOR ESTIMATING THE INTERNAL FIELDS AWAY FROM THE EDGES ISPROVIDED ALONG WITH THE RUNTIME OF EACH SIMULATION .Sim. Type N E ASM Integration Order ( x × z ) ASM Sim. Time (s) Edge Sim. Time (s) Total Sim. Time (s)Edge ×
23 28 . ×
22 27 . ×
21 27 . ×
20 27 . Full — — — —
StartSet number of internal unitcells N I and number ofedge unit cells N E = 0 Use ASM to simulate fields ofinterest in in the internal N I unit cells from the first unit cellSimulate edges with TF/SFboundaries using the tangentialfields from the ASM simulation;extract fields of interestMerge fields of interest fromthe ASM simulations andthe simulation of the edgesAssess accuracy ofresults (section VI)Converged?Stop Increase N E yes no Fig. 14. A flowchart describing the procedure of efficiently determining howlarge a finite periodic structure must be for fields to behave as though thestructure were infinite over a region of interest. TABLE IIIC
ONVERGENCE QUANTIFICATION OF A
UNIT CELL SUBWAVELENGTHGRATING ( SEE SECTION
IV-A)
AS THE NUMBER EDGE UNIT CELLS N E VARIES . T
HE DATA IN EACH COLUMN IS NORMALIZED TO THE LARGESTVALUE . N E || E ASM − E N E || TF P TF = 1 || E ASM − E N E || TF P TF = 100 1 125 0 .
169 0 . TABLE IVC
ONVERGENCE QUANTIFICATION OF AN × PHOTONIC CRYSTAL ( SEE SECTION
IV-A)
AS THE NUMBER EDGE UNIT CELLS N E VARIES .T HE DATA IN EACH COLUMN IS NORMALIZED TO THE LARGEST VALUE . N E || E ASM − E N E || TF P TF = 1 || E ASM − E N E || TF P TF = 20 0 .
608 11 1 0 . .
386 0 . . . unit cells agree with the fields generated by the infinitely-periodic structure (at a given distance from the grating).The ASM—which assumes infinite periodicity—may beused to simulate the fields on the N I + 1 internal unitcells, since convergence has been assumed in this region.Then, as described above, TF/SF boundaries can be used toinject these fields onto a simulation of a finite number of unitcells. This simulation, which is much smaller than a full-sizedsimulation of the finite structure, may be run several times withincreasingly more unit cells until convergence is reached. Notethat this algorithm requires one set of ASM simulations. Thisiterative procedure is outlined schematically in Fig. 14. A. Numerical Results: Subwavelength Grating
The algorithm above was applied to analyze the convergenceof the fields of a finite subwavelength grating to its infiniteperiodic counterpart. The characteristics and simulation detailsof the grating are described in section IV-A. The ASM wasused to simulate N I = 29 periodic unit cells (i.e. . λ op )to either side of the center unit cell. The fields generated bythe ASM were used to simulate edge unit cells as describedin section IV. E z samples were taken one wavelength λ awayfrom the grating and were Fourier transformed at the operatingfrequency. TABLE VC
ONVERGENCE QUANTIFICATION OF A SUBWAVELENGTH GRATING N I = 29 AS THE NUMBER OF EDGE UNIT CELLS N E VARIES . T
HE DATAIN EACH COLUMN IS NORMALIZED TO THE LARGEST VALUE . N E || E ASM − E N E || TF P TF = 100 110 0 . . . . . Fig. 15a shows the field samples from the ASM simulationand from the simulation of the edge unit cells. Values ofthe convergence metrics (section VI) are tabulated in table Vfor various values of N E . Based on these results, reasonableconvergence appears at around N E = 30 , corresponding toa grating of unit cells in total. The plots in Fig. 15avisually suggest as well that convergence is attained at around30 edge unit cells. Convergence at this point was confirmedby a full simulation of finite-sized grating shown in Fig. 15b.The N E = 30 case is presented in Fig. 16, which comparesthe fields over the entire computational domain from the fullsimulation of the finite grating with the fields determined fromthe algorithm given here.Note that since sampling was performed one wavelengthaway from the grating (orthogonal to the direction of period-icity), convergence of fields within the region of interest wasassured at distance of at most one wavelength from the grating. B. Computational Cost Analysis
The time taken to simulate the edges in each configurationis provided in table VI, along with the corresponding timerequired to simulate the entire grating. The table also includesthe time taken to run full-sized simulations at the various sizes.As in section IV-B, symmetry was exploited both in the ASMsimulations and in the full-size simulations. As expected, thereare significant runtime savings at high ASM parallelizationrates. VIII. S
UMMARY AND C ONCLUSIONS
This paper presents a new FDTD technique to simulatefinite periodic structures. The technique involves a two-stepprocedure: First, a simulation of a single unit cell terminatedby PBCs is performed. The ASM method is used to removeunwanted image sources generated by the PBCs. Electricand magnetic fields tangential to the periodic faces are thensampled several unit cells away from the source as functionsof time. Second, a simulation of the edges of the structureis carried out. A TF/SF boundary is used to inject the fieldsrecorded from the ASM simulation onto the edge unit cells.When the number of edge unit cells in the simulation islarge enough, there is a smooth transition between fields fromthe ASM simulation and fields from the simulation of theedges. A metric for quantifying the success of the algorithmis suggested. . . Distance y from source (cm) E z A m p lit ud e . . Distance y from source (cm) E z A m p lit ud e (a) . . Distance y from source (cm) E z A m p lit ud e . . Distance y from source (cm) E z A m p lit ud e (b)Fig. 15. Normalized magnitudes of frequency-domain E z field samples at λ units from the grating ( N I = 29 ) as the number of edge unit cells N E is varied. Figure (a) illustrates how convergence ( ± N I unit cells from thecenter unit cell) may be obtained by selecting N E high enough such that theedge simulation fields agree with the ASM fields. An inset (bottom) magnifiesweak fields in the SF region. Full simulations of finite gratings of various sizesare shown in (b). Numbers in the legends refer to the number of edge unitcells N E . An inset (bottom) magnifies fields at the boundary between theinner unit cells and the boundary unit cells. The vertical dashed lines in (a)and (b) denote the boundary between the inner unit cells and the boundaryunit cells. TABLE VIA
COMPARISON OF THE TIME TAKEN TO SIMULATE THE EDGE UNIT CELLS OF A SUBWAVELENGTH GRATING ( N I = 29 INTERNAL UNIT CELLS ) VERSUSTHE TIME TAKEN TO SIMULATE THE CORRESPONDING FINITE STRUCTURE . N E ASM Integration Order ASM Sim. Time (s) Edge Sim. Time (s) Total Sim. Time (s) Full Structure Sim. Time (s)
178 1.52 . . .
010 36 . . .
030 52 . . .
950 75 . . . . − − − x (cm) y ( c m ) − − − − − H z (dB) (a) − − − x (cm) y ( c m ) − − − − − H z (dB) (b)Fig. 16. The Fourier transform of the normalized E z field amplitudes at f op determined from (a) a full-sized simulation and (b) from merging fieldsfrom the ASM simulation ( N I = 29 internal unit cells) with fields from thesimulation of the edges ( N E = 30 edge unit cells). White dashed lines in (b)outline the boundary between the internal region and edge regions. Observehow the edges are distant enough from the internal region that the edge effectsdo not enter the internal region. Surfaces in two-dimensional simulations have two disjointends, while three-dimensional surfaces have borders along fourcontiguous edges. They thus are handled in different ways.This paper’s algorithm was used to determine the fields ofa subwavelength grating in two-dimensions and to determinethe fields of a photonic crystal in three-dimensions. When theASM simulation was parallelized, the algorithm gave rise tocomputational savings.The algorithm described here may also be used to efficientlydetermine the number of unit cells required to attain conver-gence within a certain region of the periodic structure. TheASM can be used to estimate the fields within the convergentinner region. Then, successive simulations of more and moreedge unit cells may be used to analyze convergence.Further optimizations of the algorithm described here arepossible. For instance, ASM integration orders used in theexamples above are conservative and may be relaxed in manycases.The structures considered here have identical unit unit cellswithin their inner region and edge regions. However, there isfreedom to change the structure of the unit cells arbitrarily inthe simulations of the edges. For instance, one may use thepresented technique to simulate a periodic structure terminatedby an edge taper [26]–[28].R
EFERENCES[1] A. Taflove and S. C. Hagness,
Computational Electrodynamics: TheFinite-Difference Time-Domain Method . Artech House, 2005.[2] A. J. Kogon and C. D. Sarris, “FDTD Modeling of Periodic Structures:A Review,” arXiv preprint arXiv:2007.05091 , 2020.[3] F. Yang, J. Chen, R. Qiang, and A. Elsherbeni, “A Simple and EfficientFDTD/PBC Algorithm for Scattering Analysis of Periodic Structures,”
Radio Science , vol. 42, no. 04, pp. 1–9, 2007.[4] T. Kokkinos, C. D. Sarris, and G. V. Eleftheriades, “Periodic FDTDAnalysis of Leaky-Wave Structures and Applications to the Analysis ofNegative-Refractive-Index Leaky-Wave Antennas,”
IEEE Transactionson Microwave Theory and Techniques , vol. 54, no. 4, pp. 1619–1630,2006.[5] F. Xu, K. Wu, and W. Hong, “Finite-Difference Time-Domain Modelingof Periodic Guided-Wave Structures and its Application to the Analy-sis of Substrate Integrated Nonradiative Dielectric Waveguide,”
IEEETransactions on Microwave Theory and Techniques , vol. 55, no. 12, pp.2502–2511, 2007.[6] C. Chan, Q. Yu, and K. Ho, “Order-N Spectral Method for Electromag-netic Waves,”
Physical Review B , vol. 51, no. 23, p. 16635, 1995.[7] J. M. Usoff and B. A. Munk, “Edge Effects of Truncated PeriodicSurfaces of Thin Wire Elements,”
IEEE Transactions on Antennas andPropagation , vol. 42, no. 7, pp. 946–953, 1994.[8] A. Ishimaru, R. Coe, G. Miller, and W. Geren, “Finite Periodic StructureApproach to Large Scanning Array Problems,”
IEEE Transactions onAntennas and Propagation , vol. 33, no. 11, pp. 1213–1220, 1985. [9] W. L. Ko and R. Mittra, “Scattering by a Truncated Periodic Array,” IEEE Transactions on Antennas and Propagation , vol. 36, no. 4, pp.496–503, 1988.[10] D. R. Denison and R. W. Scharstein, “Decomposition of the Scatteringby a Finite Linear Array into Periodic and Edge Components,”
Mi-crowave and Optical Technology Letters , vol. 9, no. 6, pp. 338–343,1995.[11] C. Wu, “Analysis of Finite Parallel-Plate Waveguide Arrays,”
IEEETransactions on Antennas and Propagation , vol. 18, no. 3, pp. 328–334, 1970.[12] F. Capolino, M. Albani, S. Maci, and L. B. Felsen, “Frequency-DomainGreen’s Function for a Planar Periodic Semi-Infinite Phased Array. I.Truncated Floquet Wave Formulation,”
IEEE Transactions on Antennasand Propagation , vol. 48, no. 1, pp. 67–74, 2000.[13] H. Holter and H. Steyskal, “Some Experiences from FDTD Analysis ofInfinite and Finite Multi-Octave Phased Arrays,”
IEEE Transactions onAntennas and Propagation , vol. 50, no. 12, pp. 1725–1731, 2002.[14] B. Munk and G. Burrell, “Plane-Wave Expansion for Arrays of Ar-bitrarily Oriented Piecewise Linear Elements and its Application inDetermining the Impedance of a Single Linear Antenna in a Lossy Half-Space,”
IEEE Transactions on Antennas and Propagation , vol. 27, no. 3,pp. 331–343, 1979.[15] D. Li and C. D. Sarris, “Efficient Finite-Difference Time-DomainModeling of Driven Periodic Structures and Related Microwave CircuitApplications,”
IEEE Transactions on Microwave Theory and Techniques ,vol. 56, no. 8, pp. 1928–1937, 2008.[16] F. Capolino, D. R. Jackson, D. R. Wilton, and L. B. Felsen, “Comparisonof Methods for Calculating the Field Excited by a Dipole Near a 2-DPeriodic Material,”
IEEE Transactions on Antennas and Propagation ,vol. 55, no. 6, pp. 1644–1655, 2007.[17] M. Memarian and G. V. Eleftheriades, “Evanescent-to-Propagating WaveConversion in Sub-Wavelength Metal-Strip Gratings,”
IEEE Transac-tions on Microwave Theory and Techniques , vol. 60, no. 12, pp. 3893–3907, 2012.[18] Z. Liu, S. Durant, H. Lee, Y. Pikus, N. Fang, Y. Xiong, C. Sun, andX. Zhang, “Far-Field Optical Superlens,”
Nano Letters , vol. 7, no. 2, pp.403–408, 2007.[19] C. F. Mateus, M. C. Huang, Y. Deng, A. R. Neureuther, and C. J. Chang-Hasnain, “Ultrabroadband Mirror using Low-Index Cladded Subwave-length Grating,”
IEEE Photonics Technology Letters , vol. 16, no. 2, pp.518–520, 2004.[20] C. Min, P. Wang, C. Chen, Y. Deng, Y. Lu, H. Ming, T. Ning, Y. Zhou,and G. Yang, “All-Optical Switching in Subwavelength Metallic Grat-ing Structure Containing Nonlinear Optical Materials,”
Optics Letters ,vol. 33, no. 8, pp. 869–871, 2008.[21] C. Luo, S. G. Johnson, J. Joannopoulos, and J. Pendry, “All-AngleNegative Refraction without Negative Effective Index,”
Physical ReviewB , vol. 65, no. 20, p. 201104, 2002.[22] F. Teixeira and W. C. Chew, “On Causality and Dynamic Stability ofPerfectly Matched Layers for FDTD Simulations,”
IEEE Transactionson Microwave Theory and Techniques , vol. 47, no. 6, pp. 775–785, 1999.[23] H. Holter and H. Steyskal, “On the Size Requirement for Finite Phased-Array Models,”
IEEE Transactions on Antennas and Propagation ,vol. 50, no. 6, pp. 836–840, 2002.[24] J. Zhou, T. Koschny, M. Kafesaki, and C. M. Soukoulis, “Size Depen-dence and Convergence of the Retrieval Parameters of Metamaterials,”
Photonics and Nanostructures-Fundamentals and Applications , vol. 6,no. 1, pp. 96–101, 2008.[25] J. Parsons and A. Polman, “A Copper Negative Index Metamaterial inthe Visible/Near-Infrared,”
Applied Physics Letters , vol. 99, no. 16, p.161108, 2011.[26] B. A. Munk,
Finite Antenna Arrays and FSS . John Wiley & Sons,2003.[27] M. Gustafsson, “RCS Reduction of Integrated Antenna Arrays withResistive Sheets,”
Journal of Electromagnetic Waves and Applications ,vol. 20, no. 1, pp. 27–40, 2006.[28] S. G. Johnson, P. Bienstman, M. Skorobogatiy, M. Ibanescu, E. Li-dorikis, and J. Joannopoulos, “Adiabatic Theorem and ContinuousCoupled-Mode Theory for Efficient Taper Transitions in Photonic Crys-tals,”