Space–Time Adaptive Modeling and Shape Optimization of Microwave Structures With Applications to Metasurface Design

This article presents a time-domain modeling and shape optimization framework for microwave structures, including metasurfaces, based on a nodal discontinuous Galerkin time-domain (DGTD) method. In particular, we employ an unstructured mesh and the inherent mesh refinement ability of DGTD to model multiscale geometries via adaptive mesh refinement. More importantly, we integrate a multitier local time-stepping technique into the time integration of DGTD, which significantly alleviates the cost introduced by mesh refinement. We further present a flexible parameterization approach by defining the contours of metasurface unit cells by B-spline curves, which allows us to explore a wide range of smooth shapes by only a few design variables. Besides, we adopt a polynomial chaos-Kriging (PCK) surrogate method to approximate the full-wave DGTD model, reducing the computational cost for optimization problems that require time-consuming simulations by three orders of magnitude. The B-spline parameterization and the PCK surrogate model are combined with a pattern search-based Pareto front algorithm to optimize metasurfaces for multiple design objectives. The proposed optimization framework is also applicable to other microwave structures. We demonstrate the effectiveness of our approach through the optimization of an omega-bianisotropic Huygens’ metasurface unit cell and an $E$ -plane microwave T-junction.


I. INTRODUCTION
E LECTROMAGNETIC metasurfaces are thin structures consisting of a large number of metallic and dielectric objects embedded in a host medium [1]. The shape and topology of these objects are chosen to tune the ability of the surface to manipulate an impinging electromagnetic (EM) wave in unconventional ways that go beyond those available in natural media [2], [3], [4]. The design of such structures is a multiscale problem: both the shape and the arrangement of electrically small inclusions need to be determined based Manuscript  on the performance of multiple cells on an electrically large surface. One may take advantage of the quasi-periodicity of some realizations of metasurfaces to accelerate their analysis and optimization [5]. Yet, as more complex shapes are being explored to meet more ambitious specifications, the currently popular approximate techniques will have to be replaced by comprehensive full-wave analysis tools [6]. This class of problems poses several challenges to computational EM solvers. Modeling these structures requires flexibility for mesh generation and local mesh refinement in a computationally efficient way. Moreover, current developments in the area include the use of time-varying and nonlinear elements [6], which are more suitable for the time domain rather than frequency-domain techniques. Addressing these challenges is a task of high importance for a growing number of researchers within the microwave community.
In the past 15 years, the discontinuous Galerkin timedomain (DGTD) method [7] has emerged as a viable timedomain technique and is ideally suited to these challenges. Its unstructured mesh achieves better conformity to curved boundaries than the structured mesh used in the standard finitedifference time-domain (FDTD) technique. DGTD offers a stable framework for both h (mesh)-refinement and p (order of approximation)-refinement, allowing for the use of variable order and resolution for local field approximations. Being a time-domain technique, DGTD is well suited to handle timevarying and nonlinear elements used in some metasurface geometries [8], [9], [10]. More importantly, it offers a flexible framework for the adaptation of the time step in mesh regions of different resolutions by local time stepping (LTS) [11], [12], [13], [14]. This allows us to reduce the computational cost of h-p refinement in space. With these advantages, DGTD strikes an appealing balance between accuracy and efficiency for optimizing complex EM structures, components, and systems. In this article, we present a nodal DGTD method, which adopts a second-order leap-frog time-stepping scheme and a multitier LTS technique for the accurate modeling of metasurface geometries.
In addition to the modeling of metasurfaces, we propose to use B-spline curves to parameterize the geometry of metasurface unit cells. B-spline parameterization has been widely utilized and proved efficient in the design optimization of aerodynamic wings [15], [16], [17]. With B-splines, we can generate smooth, deformable shapes, varying a small number 0018-9480 © 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
of design variables. In this work, we parameterize the contours of metallic patterns of unit cells by B-spline curves; complex metallic patterns are generated and can be deformed by a small set of design variables. The proposed DGTD solver and the B-spline parameterization are connected to pattern search [18] toward a systematic shape optimization framework. However, such a framework is computationally expensive for optimization tasks, requiring thousands of hours of simulations. To alleviate this computational cost, different surrogate-based methods have been adopted in the literature. These include the polynomial chaos expansion (PCE) method for yield-driven optimization of microwave filters [19], [20], the Kriging method for optimization of metasurfaces and microwave circuits [21], [22], [23], and a neural-network-based active learning algorithm for metasurface design [24]. In this work, we develop surrogate models based on the polynomial chaos-Kriging (PCK) method [25], [26]. PCK is a surrogate modeling technique that combines the advantages of PCE and Kriging. The global behavior of the output function of interest with respect to design variables is captured by orthonormal polynomials of PCE. The local variations of the output function cannot be captured by a low polynomial order PCE and are interpolated by Kriging. The adopted PCK surrogate model approximates the DGTD results for output functions (such as scattering parameters) so that optimization of complex geometries can be completed within minutes.
Compared with the traditional design of metasurface unit cells, which only considered rectangular geometries [3], [4], [27] or other pixel-based topologies popularly adopted in the literature [28], [29], [30], this article presents an efficient parameterization technique for the design of metasurface unit cells by defining the metallization patterns by B-spline curves. This requires only a small number of design variables to generate a wide range of curvilinear and smooth candidate patterns for metasurface unit cells. Simulations including such curvilinear geometries are computationally expensive with standard FDTD, as fine meshes are needed to minimize the pronounced staircasing errors [31], [32]. On the other hand, DGTD is well suited to such problems, as it easily handles unstructured meshes, allowing for local mesh refinement, and high polynomial order approximation within elements. Nevertheless, the computational cost of DGTD increases dramatically when the global time step is restricted by the size of the smallest element in nonuniform grids. Therefore, we employ the multirate time-stepping scheme of [33] that alleviates this computational cost by adapting different time steps to different element sizes. Then, we use a pattern search-based Pareto front algorithm to optimize the curvilinear parametric geometries to meet design objectives with respect to the transmission phase and magnitude of the unit cell. To further accelerate the optimization process, we develop PCK surrogate models that approximate the DGTD results for output functions with respect to the design variables. The small number of design variables determines a low dimensionality of the input design space, which ensures a low computational cost for the construction of an accurate surrogate model. Hence, this article builds on the analysis approach of [33] to assemble a comprehensive design optimization technique for metasurface unit cells.
Overall, this article presents a comprehensive shape optimization framework for the design of metasurfaces. We present a flexible B-spline-based parsimonious parametric model for metallization patterns of metasurface unit cells. The generated geometries are solved by an accurate and efficient DGTD simulation engine. Then, the solved candidate geometries are optimized by pattern search, and the whole optimization process is accelerated by PCK surrogate models. Thus, this article significantly extends the approach outlined in [33] that only focused on exploring DGTD with multirate stepping for the modeling of a rectangular metasurface unit cell. We, finally, validate our proposed optimization framework by optimizing an omega-bianisotropic Huygens' metasurface unit cell. We also demonstrate the potential of this approach for other microwave circuits by optimizing a microwave T-junction.
The rest of this article is structured as follows. In Section II, we review the formulation of the nodal DGTD method and LTS. We also discuss the implementation of the perfectly matched layer (PML) in DGTD. In Section III, we present the principles of B-spline parameterization. We elucidate the process of parameterizing geometries of metasurfaces' unit cells through examples. In Section IV, we show how a surrogate model for metasurface unit cells can be combined with a pattern search-based Pareto front algorithm to accelerate shape optimization, and we illustrate our optimization framework through a flowchart. In Section V, we apply our optimization process to optimize unit cells of a Huygens' metasurface. In Section VI, we use our approach to optimize the geometry of a waveguide T-junction. Section VII summarizes our contributions.

II. NODAL DISCONTINUOUS GALERKIN TIME-DOMAIN METHOD
For completeness, we review the basic steps of the formulation of the nodal DGTD technique. Consider Maxwell's equations for an isotropic medium We can convert Maxwell's equations into the form of a conservation law [7] α ∂u ∂t where I 3×3 represents the 3 × 3 identity matrix. The flux term F = (F x , F y , F z ) is an array of three vector components F x , F y , F z . The three components are expressed as Now, we premultiply (3) by α −1 and move the flux term to the right-hand side to obtain ∂u ∂t = −α −1 ∇ · F(u). (6) This equation serves as the starting point for the formulation of DGTD. We assume that the simulation region is tiled using K nonoverlapping tetrahedral elements as The unknown local field solution inside each element k is expanded in terms of a set of N p Lagrange polynomials of order N, where E x,i represents the value of E x at the interpolation node x i , with the other field components defined similarly. Note that N p is the number of nodal points of each element and N is the order of each polynomial. In 3-D cases, these two are related as N p = (N + 1)(N + 2)(N + 3)/6. We test (6) with N p Lagrange polynomials and integrate it over a tetrahedral element k enclosed by a surface ∂ k to obtain Then, by applying the divergence theorem twice on the right-hand side of (9), we obtain wheren represents the normal unit vector of the surface ∂ k , pointing outward from volume k , and u * refers to the intermediate fields defined on the interface, which are determined by the fields on both sides of the interface. Also, F(u) ·n and F(u * ) ·n represent the flux on the surface ∂ k . F(u)·n is evaluated by the variables inside k and, thus, is easy to be determined, and F(u * ) ·n represents the numerical flux depending on the intermediate fields u * . These two can be rewritten in terms of the field vectors as follows: In this article, we adopt an upwind flux scheme based on the solution of a local Riemann problem [34], [35] to determine the flux terms in (11) where Z = √ μ/ε and Y = √ ε/μ represent the intrinsic impedance and admittance of the medium. The addition and jump operators are defined as where the superscripts − and + refer to the variables in element k and neighboring elements, respectively.

A. Field Update Equations
For time integration, we employ a second-order accurate, leap-frog time-stepping scheme. This results in a scheme that is similar to FDTD. For example, the update equations for E x and H x are given as In (14) and (15), f is the face index of the tetrahedral element, and N f equals four. M −1 is the inverse of the mass matrix that is defined as The stiffness matrix is defined as The evaluation of numerical flux is based on the face mass matrix The stability condition for this scheme using first-order basis functions is given by [13] where V i is the volume of the i th element and S i,k is the area of the kth face of the i th element, v i p = 1/ ε i μ i , and j refers to the element adjacent to element i .

B. Local Time Stepping
For many EM designs that comprise complex or subwavelength geometrical details, the computational cost will be significantly increased if taking the standard update scheme (14) and (15), as the time step is restricted by the size of the smallest element in nonuniform meshes. To alleviate this computational cost, we integrate a multitier LTS algorithm [12] into our solver. It allows for the adaptation of the time step in mesh regions of different resolutions. This flexibility is depicted in Fig. 1, where different time steps are assigned to regions of different mesh densities, for a 9:1 ratio of maximum-to-minimum time steps.
Here, we summarize the basic concepts of the LTS method. The mesh is partitioned into N tiers based on their local stability condition of (19). The elements are sorted from the smallest t i,max to the largest. The time step for k − th tier is chosen as (2m + 1) k t min , k = 0, 1, . . . , N − 1, where t min is the minimum time step determined by the smallest element. We choose m = 1 so that each tier has a time step three times larger than the previous tier. Note that it is unnecessary to define clear interfaces between different tiers; the partition of each element is determined only by its local stability condition. For example, element i is included in To explain the concept of the LTS update scheme, we assume only two tiers (N = 2), tier 0 and tier 1, with time step of t min and 3 t min . The update sequence is shown in Fig. 2.
For elements not located at the interface of the two tiers, update equations, such as (14) and (15), can be applied with the local time step of their tiers without any problem. However, if the elements are located at the interface of the two tiers, field values at unavailable time steps will be needed for the updates. Instead of interpolating, these field values are approximated by their closest in time values available. Notably, the information exchange between tiers is related to the field components E x , E y , E z and H x , H y , H z . Therefore, these six field components are stored as global matrices; dozens of other matrices, e.g., matrices for materials, face Update sequence for multitier time stepping with two tiers; tiers 0 and 1 are updated by t min and 3 t min , respectively. normals, and mass matrices, are all stored by tiers to mitigate the global computational cost.

C. Perfectly Matched Layer
The simplest way to truncate the computational domain is the first-order Silver-Müller absorbing boundary condition (ABC). However, the method only works well for normal incident waves. For open problems, we use the PML [36] to terminate the computational domain. Instead of using the Runge-Kutta method proposed in [36] for the time integration, we use the leap-frog update scheme, which is compatible with the presented multitier LTS.
In the FDTD method, the standard implementation of PML takes polynomially scaled parameters to reduce numerical reflections near the PML interface. The conductivity σ (k) is typically chosen as a graded profile in the form of where k = 0 is the coordinate of the PML interface, d represents the layer's width, σ max is the maximum value of σ at k = d, and m is a power typically chosen between 2 and 4. However, if the same profile is employed in DGTD, a late stability problem was observed in [37]. The stability issue can be solved by taking a constant profile (m = 0) for the conductivity. Besides, we found that a scaled profile introduces larger reflections from the PML layer than a constant profile. This conclusion has also been verified by other authors [13]. Thus, we take a constant profile (m = 0) for conductivity in the PML layer in this work.

III. B-SPLINE PARAMETERIZATION OF METASURFACE GEOMETRIES
In this section, we discuss the geometry parameterization of metasurface unit cells by B-spline curves. A brief introduction to B-spline curves is included in the Appendix.
In the following, we demonstrate the B-spline parameterization for a unit cell of an omega-bianisotropic Huygens' metasurface [4]. The metasurface unit cell is a strip and split ring resonator (SRR) structure, as shown in Fig. 3.  The strip is printed on the front face, and the SRR is printed on the back face of a dielectric slab. The unit cell utilizes the strip and SRR to control its electric and magnetic response, respectively. These are realized by modifying the width and position of the strip and the radius of the SRR. However, Chen and Eleftheriades [4] only considered rectangular geometries for the strip and the SRR, limiting the design space. In the following, we will demonstrate how to parameterize this standard design via B-spline curves.
We first consider the design of the metallic strip of Fig. 3. The initial rectangular topology is shown in Fig. 4(a), and the lines of symmetry are shown in dotted lines. In Fig. 4(b), we propose to represent the two sides of the rectangle as B-spline curves, which are plotted in dashed lines. In Fig. 4(c), we added five control points on one side to control the B-spline curve, and the geometry can be easily perturbed by moving the control points. In Fig. 4(d), we show a smooth hourglass design obtained by simply moving the control point in the middle. Due to symmetry, we use only three control points to parameterize this structure. To further reduce the number of design variables, we only allow the control points to move in the horizontal direction. As a result, the horizontal coordinates of the three control points are taken as the design variables.
Then, we introduce our strategy of parameterizing the SRR. Instead of taking the rectangular SRR in Fig. 3, we start with a curvilinear ring by replacing the contours with B-spline curves,  as shown in Fig. 5(a). Two B-spline curves are utilized to control the internal boundary and the external boundary of the SRR, respectively. The distance between the two B-spline curves is fixed for a constant width of the SRR. We use eight control points for each B-spline curve to modify the shape of the SRR. Due to symmetry, only five design variables are used to parameterize this structure. To further reduce the number of design variables, we only allow the control points to move along the designed trajectories, shown by green solid lines at 45 • from each other. Thus, the design variables of the SRR are the distances of the five control points from the origin. Three smooth design examples, as shown in Fig. 5(b)-(d), are generated by simply moving the control points.
Once we obtain the designed geometry, we need to generate the corresponding mesh needed for DGTD simulations. With the adopted meshing tool, an automated built-in meshing process for B-spline curves is implemented. In Fig. 6, we show the corresponding mesh of Figs. 4(d) and 5(c). We refine the mesh near the B-spline curves to capture the subtle modifications introduced by B-splines.

IV. SURROGATE MODEL AND OPTIMIZATION FRAMEWORK
In this section, we present our optimization framework by combining the DGTD method, B-spline parameterization, and a pattern search algorithm. Moreover, we derive a surrogate model that can accelerate the optimization process. To approach the desired objectives, we utilize the Hooke-Jeeves algorithm [18] also known as "pattern search." However, to optimize metasurface unit cells, thousands of hours of simulations are required. To solve this problem, we derive a surrogate model that approximates the full-wave DGTD simulation needed to compute the objective functions.

A. Construction of the PCK Surrogate Model
A surrogate model is an approximation function that mimics the original system but can be evaluated much faster [38]. A robust surrogate model can partly or completely substitute the full simulation to save computation time. In this work, we use the PCK surrogate model [25] to approximate the original system. The method combines Kriging [26] and PCE [39]. The pure PCE surrogate model requires highorder polynomials to capture fast, local variations of the output function of interest with respect to the design variables. However, the required number of training samples of PCE increases exponentially with the polynomial order. Instead, we enhance the resolution of a low-order PCE model via Kriging interpolation.
To compute the surrogate of the original system using PCK, we first need to construct a probabilistic input model that specifies the distribution of the input variables. Considering the parameterization of the metasurface unit cell that we have shown in Section III-B, we use three and five design variables to control the shape of the strip and SRR, respectively. These eight design variables are the input variables needed for the surrogate model. In this work, we specify uniform distributions for all input variables. We generate a number of random input vectors in the form of x = [x 1 , x 2 , . . . , x 8 ], which consists of the eight uniformly distributed input variables between the following bounds: 0.6 mil ≤ x i ≤ 15 mil, i = 1, 2, 3 12 mil ≤ x i ≤ 25.5 mil, i = 4 13.6 mil ≤ x i ≤ 25.5 mil, i = 5, 6, 7, 8 (22) in which x 1 -x 3 are the design variables of the strip, and x 4 -x 8 where S PCK (ξ Val i ) is either the phase or magnitude of S 21 at 20 GHz computed by the derived PCK surrogate model of the i th sample, and Y Val i is the corresponding output computed by full-wave DGTD simulation. In order to obtain an accurate surrogate model, we need to specify a suitable number of samples of the training dataset. To that end, the convergence of the error via (23) to the number of training samples is studied. After we obtain a convergent surrogate model, the full-wave DGTD analysis can be replaced by the surrogate model to evaluate the desired output of S 21 .

B. Optimization Process
The flowchart of the optimization framework is depicted in Fig. 7. We first parameterize the geometry by B-spline curves; this results in an array of design variables. Then, some initial design variables are chosen, and the corresponding geometry and mesh are generated. We compute the S-parameters of these structures by the surrogate model and determine the optimal geometry of the current iteration. If the objectives (typically the phase and magnitude of S 21 ) converge, the final geometry is computed by full-wave DGTD simulations for validation. Otherwise, this geometry is exported to pattern search as an initial shape for the next iteration until the objectives converge.

V. APPLICATIONS TO METASURFACE UNIT CELLS
In the following, we consider the analysis and optimization of an omega-bianisotropic Huygens' metasurface unit cell, as proposed in [4]. The unit cells are designed to achieve reflectionless refraction from normal incidence to 71.8 • at 20 GHz. For the design of metasurface unit cells, the key is to introduce a resonance near the target frequency where a transmission coefficient close to unity is achieved. More importantly, a sharp transmission phase change near the resonance can be introduced, and different desired phases can be obtained by perturbing the geometry of the unit cell.
In [4], several different geometries were designed to achieve the objective phase and magnitude of the transmission coefficient. The authors analyzed a large library of possible unit cells by modifying the length of the copper strips, the position and number of the air gaps, and the radius of the SRRs. Then, unit cells with specified transmission/reflection coefficients were used to assemble the metasurface. However, this process was conducted by trial and error, and the unit cells were designed only using rectangular shapes. Hence, the design space was constrained, limiting the possibility of achieving designs with broader bandwidth. To enlarge the design space, we introduce B-spline curves to define the boundaries of Fig. 7. Shape optimization process for metasurface unit cells: the optimization starts with the B-spline parameterization, followed by a process of geometry and mesh generation. A PCK surrogate model is used to compute the S-parameters of the generated geometries. The optimal shape of the current iteration is obtained to test the convergence of objectives. New update design variables will be generated if the objectives do not converge. the metallic strips. Besides, we derive surrogate models to accelerate the optimization process.
In this case study, we start with an initial topology similar to Fig. 3 where the front face includes a metallic strip and the back face includes an SRR. However, instead of using rectangular geometry, we utilize the B-spline parameterization process shown in Figs. 4 and 5. The 2-D layout of the structure is shown in Fig. 8, where the B-spline curves and the control points are marked. The thickness of the dielectric slab is 0.635 mm. For the strip, the control points are allowed only to move in the horizontal direction. For the SRR, the points are allowed to move along the green rays starting from the origin (see Fig. 8). Due to symmetry, we have eight design variables in total to control the shape of the metasurface unit cell.

A. LTS-DGTD Modeling of Metasurface Unit Cells
Before presenting the optimization of the metasurface unit cell, we would elucidate the mesh partition of LTS and the importance of this feature of the method for metasurface modeling. We consider the scattering parameter calculation of the initial unit cell of Fig. 8 inside a waveguide geometry, as shown in Fig. 9. We use a relative permittivity of 11.2 for the dielectric layer, and we represent the metallic strips by perfect electric conductors (PECs). The scattering parameters are computed by probing the vertical component of the electric field in front of and after the unit cell. The waveguide is terminated by perfect magnetic conductors (PMCs) in the lateral direction, PECs in the vertical direction, and PML absorbers in the longitudinal direction. A 0-30-GHz Gaussian pulse is used to excite the TEM mode.  For this problem, a mesh of 9601 tetrahedral elements is generated and shown in Fig. 10. Progressively refined elements are formed around the strip/SRR structure, especially between the air gap where fields change rapidly, while coarser elements occupy the space around the structure and in the PML area. These elements are divided into three tiers based on (20). The resulting three tiers are shown in Table I, where t min is the smallest time step determined by the smallest element. For this mesh, small elements near the strip/SRR are partitioned into Tier I, elements occupying the waveguide and the PML area are mostly partitioned into Tier III, and the rest of the elements are partitioned into Tier II. In the following simulation, we use the second polynomial order, i.e., p = 2 within each element.
In Table II, we compare the execution time with the standard DGTD and the DGTD with LTS. The DGTD with LTS takes 10 380 s for a single run of the simulation, whereas the standard DGTD requires almost two times more execution time.

B. PCK Surrogate Models for Metasurface Unit Cells
In the following, we derive PCK surrogate models for the design of metasurface unit cells. As discussed in Section IV-A, we generate a set of samples, where the input of each sample is a design vector that contains eight uniformly distributed random numbers in the form of x = [x 1 , x 2 , . . . , x 8 ], with a set of predefined limits, as in (22). For the output, the optimization of the unit cell requires both the phase and magnitude of the transmission coefficient. We construct two independent surrogate models for phase and magnitude, respectively. We take the phases and magnitudes of the transmission coefficient at 20 GHz as outputs for both surrogate models. We set the polynomial degree of PCE to one, which results in nine orthogonal polynomial basis functions.
The accuracy of the derived surrogate model is validated via (23). To that end, we generate another 40 test samples, independent of the training samples, to evaluate the difference between S 21 computed by the surrogate model and the full DGTD simulation. We derive several different surrogate models with an increasing number of training samples to study its convergence. The errors of magnitude and phase versus the number of training samples are plotted in Fig. 11. We require the errors to be less than 10% for both phase and magnitude. The error of the phase decreases below 10% after about 450 training samples, and the error of the magnitude is reduced to 10% after 300 training samples. In the following work, we take surrogate models trained by 800 samples to compute S 21 at 20 GHz for both phase and magnitude. Note that, with B-spline curves, we only need eight input variables to control a complex curvilinear shape. This small number of input variables makes the PCK surrogate model more efficient for this problem, as the number of training sets required for PCK increases with the number of input variables.

C. Shape Optimization of Metasurface Unit Cells
Once we obtain the desired surrogate models, we can optimize metasurface unit cells. The design goal of this unit cell is to achieve the desired phase of transmission coefficient and maximize the transmission coefficient at 20 GHz. For demonstration, we perform an optimization where the target phase is set to −40 • and the target magnitude to unity. The optimization can be conducted by minimizing the following two objectives: The optimization requires the two objectives to be minimized simultaneously. However, the standard pattern search  can only conduct optimizations with a single objective function. To solve this problem, we adopt a pattern search-based Pareto search algorithm to facilitate multiobjective optimization. Note that, in the Pareto search, we can apply weights to the objectives. In the following optimization, we set equal weights to the percentage error of the two objectives; the solution that minimizes the following equation is considered the optimal solution of the current iteration: Fig. 13 shows the variation of the computed phase and magnitude with iterations. The results converge after the third iteration, the convergent objective magnitude equals 0.998, and the objective phase equals −40.17 • . Note that the Pareto search multiobjective algorithm does not require an initial shape; the tradeoff is that each iteration takes hundreds of simulations. However, the surrogate model can compute the output of these simulations within seconds. For this optimization, the numbers of computations needed for the first, second, third, and fourth iterations are 60, 261, 295, and 362 respectively. In total, it takes at least 978 computation runs. If this optimization was conducted with full DGTD simulations, the total execution time needed would be around 2000 h. On the other hand, the surrogate model-based optimization only takes 2 min. The obtained surrogate models can be used to optimize unit cells with different target phases.
We show the shape modifications of the unit cell during the optimization process for the strip and the SRR in Fig. 12.
Note that geometries of the first and third iterations are different; the design vectors for the first and third iterations are .48] mil, respectively. Small changes in the fifth and sixth design variables lead to improvements in both phase and magnitude, as shown in Fig. 13.
Finally, we apply DGTD to validate the design obtained by the surrogate model of the third iteration, and the corresponding S-parameters are plotted in Fig. 14. The magnitude of S 21 at 20 GHz is 0.969, and the phase is −42.64 • . A good agreement has been achieved between the surrogate model and the full-wave DGTD simulation.
For comparison, we also compute the S-parameters of the rectangular unit cell in Fig. 3. The rectangular unit cell is also designed for a target phase of −40 • . As shown in Fig. 14, both designs reach a transmission angle close to −40 • . Moreover, our design achieves a much broader bandwidth than the rectangular design. In terms of S 11 , a −10-dB bandwidth close to 20% is achieved with our design, whereas S 11 of the rectangular design failed to stay below −10 dB through the simulated frequency range. On the other hand, our design does so for 17-21 GHz. This indicates that B-spline parameterization enlarged the initial design space, thus leading to a unit cell of superior performance.
In [4], a set of unit cells with transmission phases varying from −π to π were designed. These unit cells were connected to assemble a metasurface that could manipulate the wavefront from normal incidence to 71.8 • ; differences in transmission phases are essential to tune the ability of the metasurface to bend the wavefront in unconventional ways. Based on the proposed optimization approach and the initial B-spline topology of Fig. 8, we design eight unit cells with target phases ranging from −20 • to −90 • and target magnitudes of unity. Designs for other target phases, not discussed in this work, can be explored by varying the initial B-spline arrangement of Fig. 8. Then, we analyze the eight obtained designs via fullwave DGTD simulations to confirm their performance. The obtained unit cells and comparisons in phase and magnitude are shown in Fig. 15, where the horizontal axis refers to the index of the eight unit cells. For example, "case 3" corresponds to the unit cell with the target transmission phase of −40 • , as discussed above.
For the phase of S 21 , the solutions obtained from the surrogate model-based optimization match well with the target phase. A good agreement between the DGTD results and  the surrogate model is observed for all the designs. For a magnitude of S 21 , cases 1 and 2 achieve a magnitude close to unity, and the magnitude decreases as the target phase decreases. This is expected, as these target phases are achieved with the cost of going away from the resonance, where the highest transmission magnitude is achieved. A good agreement is observed between the DGTD results and the surrogate model for cases 1-5. The differences for cases 6-8 are within 2 dB.

D. Summary
In summary, metasurface unit cells with different transmission phases and broad bandwidth are obtained by the proposed approach. The LTS-DGTD method serves as an efficient and accurate simulation engine to model the multiscale and curvilinear unit cells. The B-spline parameterization provides a wide range of candidate shapes that are explored by a Pareto search optimizer to obtain designs with desired transmission phases and magnitudes. The whole optimization process is dramatically accelerated by a PCK surrogate model. The combination of these advanced features leads to a flexible and computationally efficient optimization process for the design of metasurface unit cells.

VI. APPLICATIONS TO OTHER MICROWAVE STRUCTURES: WAVEGUIDE T-JUNCTION
In this section, we explore the potential of the optimization process for other microwave circuits. We consider the optimization of a microwave T-junction (see Fig. 16), an example also treated in [40]. We show the 3-D geometry and its top view in Fig. 16(a) and (b), respectively. The T-junction, which is terminated at PMCs at the top and bottom, is bounded by PECs in all other directions. To truncate the computation domain of the T-junction, we terminate the three ports with the Silver-Müller ABC. The design goal is to minimize the return loss at the input port over the operating bandwidth. At the input port, we use a Gaussian electric current source in the y-direction to generate a TEM wave. A similar objective function as in [40] is adopted to estimate the return loss at the input port where E i y and E r y are the incident and reflected electric fields probed near the input port, m is the probing plane, and T is  Mesh of 29 741 tetrahedral elements generated for the initial T-junction. the total simulation time. The numerator represents the total reflected energy, and the denominator represents the incident energy. To apply the proposed B-spline parameterization technique, we choose a B-spline curve that consists of several control points on the PEC wall, as shown in Fig. 16(b). Due to symmetry, we take a B-spline on the left half of the structure and mirror it to the right. In this case, a B-spline curve consisting of 11 control points is taken, and five of them are chosen as the design variables due to symmetry. The endpoints are kept unchanged to connect the optimized curve to the fixed boundaries. The control points are only allowed to move along the x-direction. Once the B-spline curve is determined, we can project it to the z-direction with the thickness of the T-junction to form the full 3-D structure.

A. LTS-DGTD Modeling of T-Junction
For the initial T-junction, a mesh of 29 741 tetrahedral elements is generated and shown in Fig. 17. Progressively   TABLE III   DIVISION OF THE 29 741 ELEMENTS OF THE MESH OF FIG. 17   TABLE IV  EXECUTION TIME COMPARISON WITH THE 19. Contour comparison between the initial shape and the final shape. The shape obtained in [40] is also shown for comparison. refined elements are generated to reconstruct the geometrical details around the B-spline boundaries. Larger cell sizes, close to (λ/15), are adopted in the region away from the B-spline boundaries. In the following simulation, we use the first-order polynomial, i.e., p = 1 within each element. Applying the mesh partition of (20), we have three tiers shown in Table III, where t min is the time step determined by the smallest element. In Table IV, we compare the execution time of a single simulation for this problem with the standard

B. Shape Optimization of T-Junction
The optimization of this geometry requires less than 100 simulation runs, and each simulation takes about 7 min. As a result, the computational cost is not as expensive as the optimization of metasurface unit cells. Thus, we do not need to generate a surrogate model for this problem. We use the fullwave DGTD solver to compute the objective function of (26). Fig. 18 shows the variation of the objective function value with the iteration number. The objective function rapidly decreases and converges after 12 iterations.
In Fig. 19, we show the final shape and compare it to the optimal shape from [40], which was determined by a design sensitivity analysis (DSA) method. The two convex contours have similar shapes, and the endpoints at the contour minima are almost the same. However, Chung et al. [40] took 30 design points to achieve such a contour, whereas we only have five. Besides, we found two local bumps in our final design. The two local bumps are close and tend to point in opposite directions. Fig. 20 shows the corresponding mesh that we used to simulate the final design. The curvilinear boundaries are resolved by a relatively fine mesh to reconstruct the geometrical details.
In Fig. 21, we compare the scattering parameters for the initial shape, the final shape of this work, and the DSA [40]. The scattering parameters are computed by probing the vertical electric field at the input port and one of the output ports. With our final design, the return loss greatly decreases compared with that of the initial shape. Moreover, we found that the return loss of our final design is smaller than that of the design in [40]. This is expected because the B-spline parameterization helps enlarge the design space, which contributes to the discovery of better designs.

VII. CONCLUSION
This article presents a comprehensive approach to the shape design and optimization of microwave structures with an emphasis on metasurface unit cells. We chose DGTD as the simulation engine employed, leveraging its h-p adaptation and mesh flexibility along with local time stepping. These features of DGTD are well suited to the challenges that arise in the modeling of current metasurface geometries. We also used B-splines to parameterize the geometry of metasurface unit cells. This parameterization allowed us to produce a wide range of smooth candidate geometries that a Pareto optimization algorithm explored to meet design objectives with respect to the transmission phase and magnitude of the unit cell. The optimization was dramatically accelerated by a surrogate model that reduced the computation time of the optimizer to few minutes. Finally, we demonstrated the wide applicability of this approach to microwave circuit optimization via an example of a microwave T-junction. In all cases, we showed that the ability to explore a wide range of shapes with a small number of design variables, via B-spline parameterization, led to designs that met their functional requirements over a broad bandwidth in a computationally efficient manner.

APPENDIX BRIEF OVERVIEW OF B-SPLINE CURVES
The B-splines can be utilized to parameterize smooth curves, surfaces, and volumes. The parametric equation of a B-spline curve of the pth order is defined as a weighted sum of basis functions where {x i } N i=1 are the coordinates of the control points and N p i are the B-spline basis functions of order p (they are p − 1 degree spline polynomials) [16]. The range of the basis functions is determined by the knot location t i that is defined over the range [0, 1] = [t 1 , t p+N ]. The basis functions, N p i , are defined recursively using the parameter t and the knots . (28) In this work, we adopt fourth-order basis functions (cubic splines), which is a common choice to ensure smoothness. Besides, we want our B-spline curves to pass through the two end control points, which is useful when several B-spline segments are connected to generate a single curve. To that end, we consider an open uniform knot vector (0, 0, 0, 0, t 5 , t 6 , . . . , t N , 1, 1, 1, 1) with The multiplicity of 0 and 1 at the two ends of the knot vector ensures that the fourth-order B-spline curve passes through x 1 and x N . In between, the curve passes near the control points but does not go through them. Note that, in special cases, when all the control points line up as a straight line, the B-spline curve turns to a straight line that passes through all of them. Once we obtain the knot vector via (29), we can determine the corresponding basis functions by (28).
To explain how to generate a B-spline curve, we parameterize a simple B-spline curve consisting of five control points step by step. We first determine the open uniform knot vector by substituting N with 5 in (29), and we obtain the knot vector (0, 0, 0, 0, (1/2), 1, 1, 1, 1). Then, we determine the five corresponding basis functions by (28) with the obtained knot vector. Once we obtain the basis functions, a parametric equation for a B-spline curve can be generated via (27). In Fig. 22, we plot a B-spline curve consisting of five control points with the basis functions that we obtained. The curve can be easily rearranged by simply moving the control points without modifying the basis functions.