Near real-time full-wave inverse design of electromagnetic devices

Domain: neuroscience
Published: 2026-02-14
Editorial outcome: major_revision

Abstract
Inverse design enables automating the discovery and optimization of devices achieving performance significantly exceeding that of traditional human-engineered designs. However, existing methodologies to inverse-design electromagnetic devices require computationally expensive and time-consuming full-wave electromagnetics simulation at each iteration or generation of large datasets for training neural-network surrogate models. This work introduces the Precomputed Numerical Green Function method, an approach for ultrafast electromagnetic inverse design. The static components of the design are incorporated into a numerical Green function obtained from a single fully-parallelized precomputation step, reducing the cost of evaluating candidate designs during optimization to only being proportional to the size of the region under modification. A low-rank matrix update technique is introduced that further decreases the cost of the method to milliseconds per iteration without any approximations or compromises in accuracy. This method is shown to have linear time complexity, reducing the total runtime for an inverse design by several orders of magnitude compared to using conventional electromagnetics solvers. The design examples considered demonstrate speedups of up to 16,000x, shortening the design process from multiple days to weeks down to minutes. The approach enables practical and ultrafast design of complex structures that are prohibitively time-consuming for prior inverse design methods.

Introduction
Electromagnetic devices are an indispensable part of daily life, playing key roles in telecommunications, radar, sensors, biomedical devices, and more. The conventional process of electromagnetic design is heavily reliant on human intuition and experience, and the iterative nature of design is time-consuming and resource-intensive. As such, inverse design techniques—algorithmic approaches for the discovery and optimization of devices or structures yielding desired functional properties—have attracted significant focus across many disciplines, including radio frequency (RF) or millimeter-wave (mmWave)1–14 nanophotonics and optics15–25 and materials and structural engineering26–30. The properties of interest are encoded as objective functions that are extremized via optimization methods. The appeal of the inverse design paradigm owes to its capability for broad exploration of design spaces with many degrees of freedom, enabling the synthesis of novel devices achieving performance superior to that of conventional designs.

In gradient-based inverse design approaches, optimization is performed by iteratively following the gradient of the objective function computed over the space of input parameters. Such methods are liable to converge to local extrema, and many inverse design runs at random starting configurations may be required before satisfactory results are attained. Moreover, a gradient may not be available due to discrete-valued input parameters, such as metal conductivities and substrate dielectric constants; allowing such parameters to vary continuously may result in physically infeasible designs. As a result, gradient-free optimization approaches, such as genetic algorithms1,2,15 and particle swarm optimization3–5,31, have been introduced, enabling wider design space coverage. However, a prominent limitation of both gradient-based and gradient-free techniques for electromagnetic design is that full-wave field simulations are required to evaluate the objective function at each optimization iteration. Even with the fastest commercially-available solvers, such as Ansys HFSS or CST Studio Suite, single simulations often take tens of minutes to hours to run, even for structures of only moderate complexity.

As such, objective function evaluation is typically the rate-limiting factor for design throughput, and mitigating this has been the subject of much work. For instance, adjoint methods6,16,17 for gradient-based approaches allow the gradient to be computed with only two field simulations per iteration. Alternatively, to dispense with simulation entirely during optimization, machine-learning techniques, which construct surrogate models that allow performance to be predicted from the input parameters, have garnered widespread attention within and beyond electromagnetic design7–9,11–14,18,20,23,24,27–29. While neural-network surrogates can greatly reduce optimization time, the overall process is constrained by computationally expensive training as well as the large number of simulations needed for adequate design space coverage when generating training datasets. Although approaches such as transfer learning have been introduced to enhance training efficiency, the training phase, inclusive of the generation of the large dataset (on the order of 10,000 to >1 million simulations7,9,11,18,29), may nonetheless require multiple days to weeks8,11,13,23. Furthermore, over the full design space, there is no guarantee of accuracy given valid input parameters23,32; that is, a design predicted as optimal by the surrogate may yield completely different performance when verified with full-wave simulation or measurement of a fabricated device.

This paper introduces the precomputed numerical Green function (PNGF) method, an approach for inverse design of electromagnetic devices with rapid, direct, approximation-free objective function evaluation at each optimization iteration. Structure(s) to be optimized in the design are replaced with equivalent electric current densities, allowing the interaction between the static, unchanging portions of the device and the dynamic optimization region to be represented by a numerical Green function matrix obtained from a single fully parallelized precomputation step. Once found, the numerical Green function may be used to evaluate the performance of any candidate device within the design environment.

Reusing computations to accelerate solutions to electromagnetic problems is a widely used paradigm. For example, domain decomposition techniques—speeding up the analysis of a computational domain by solving subdomains and coupling adjacent subdomains together with boundary conditions at the interfaces—are often accelerated by reusing previously computed quantities. Applications have included perturbation approaches with surface currents in a method-of-moments (MoM) formulation for analyzing large aperiodic antenna arrays33, optimization of photonic crystal structures with a Dirichlet-to-Neumann discretization34, and speeding up the simulation of large metallic platforms with local variations35, among others. While the PNGF method shares broadly similar motivations, it is distinct from domain decomposition approaches. The numerical Green function reduces the problem size to only the number of unknowns in the optimization region, and the remainder of the domain, including the static portions of the device, is incorporated into the numerical Green function itself.

While this work demonstrates PNGF on pixelated metallic structures using the direct binary search (DBS) global optimization algorithm, the method is readily applicable to dielectric optimization problems and can also be used with other optimization approaches. Many optimizers, including DBS, make a small modification to the designate each iteration. As such, low-rank update methods may be utilized to accelerate objective function evaluation significantly without trade-offs in accuracy. The cost of evaluation is linear with the size of the optimization region, which reduces the total runtime of the full inverse design by multiple orders of magnitude, from several days or weeks (with commercial solvers) to minutes. Unlike neural-network-based surrogate models, the PNGF method yields a highly accurate solution, which matches those obtained from the full-wave electromagnetics solver used for precomputation with multiple digits of precision and is correct for every design in the feasible set, without training.

The PNGF method is applied to design three example devices: an ultrawideband 30GHz substrate antenna with 40% fractional bandwidth, a 6GHz planar switched-beam antenna (SBA) whose beam is switchable over a 70° angle, and a broadband short-length transition between a microstrip feedline and a substrate-integrated waveguide (SIW). When PNGF is utilized as the solver for DBS, speedups of up to four orders of magnitude in the optimization time versus DBS using the fastest commercial software (e.g., HFSS and CST) are obtained. The SBA and microstrip-SIW transition are fabricated, and the measured scattering parameters of the devices and the radiation pattern of the SBA agree closely with predicted simulation results.

Results
Current equivalence
An electromagnetic device to be optimized typically encompasses a predefined optimization region, where the goal of design is to find a configuration of structure(s) in the optimization region that yields desired electromagnetic properties. Additional components of the device, such as dielectric structures, ground planes, feedlines, and air gaps, that lie outside the optimization region are constrained to be static and are excluded from being modified. To model the device, 3D space is discretized into a set of voxels (typically by using a finite-difference or finite-element algorithm) in a simulation environment that encompasses the device and its immediate surroundings, such as an air box, as well as models of excitation source(s), such as lumped ports and incident fields.

Each optimization iteration would be significantly accelerated if candidate designs could be assessed as linear combinations of solutions to simpler designs. However, such solutions do not obey superposition. Although a structure can be viewed as a combination of simpler parts, the sum of the fields scattered by each part alone (under the same excitation) would not match the fields of the full structure due to scattering interactions between substructures. As such, traditional optimization approaches have required full-wave simulations to reevaluate the entire environment, including the static components, from scratch at each iteration.

To overcome this limitation, we apply the current equivalence theorem to represent any given configuration of the optimization region with polarization densities, which create fields identical to those that would be generated by the original device in response to an excitation. Arbitrary current densities and the fields that they produce are linear and obey superposition. This allow candidate designs to be evaluated by solving a linear system whose number of unknowns is equal to the size of the optimization region only, as opposed to the full simulation environment. The interactions between the optimization region and the static components of the design are encoded into a numerical Green function, a matrix that maps current densities to electric fields in the optimization region.

We introduce the current equivalence concept using pixelated metallic devices, in which a grid of tiles, each of which may be filled with metal or left open, forms the optimization region. Such devices are readily modeled via finite-difference methods, where each tile comprises a rectangular array of Yee cells, as shown in Fig. 1a, and the polarization densities and electric fields found with current equivalence are co-located on the edges of Yee cells. The concept of solving for the unknown polarization densities with a candidate design is depicted in Fig. 1b, and utilizing the numerical Green function is shown in Fig. 1c. Although we only consider metal in our optimization region for practical demonstration, this technique can be used with any arbitrary complex permittivity distribution to represent lossy metals and/or dielectrics, as discussed in Supplementary Note 4.

Numerical Green functions
A list of symbols used throughout this work is provided in Supplementary Table 1. We seek to replace any given configuration of the optimization region with an equivalent effective polarization density Jp(r) = ε0(εr(r)−1)E(r), where ε0 is the free-space permittivity and εr(r) is the material permittivity, such that the fields E produced by Jp in an empty optimization region are identical to those with the original metallic tiles in response to an excitation Einc. The conductivity σ(r) at points r throughout the optimization region is either zero (free space) or infinity (metal represented by perfect electrical conductor (PEC) material). Defining an auxiliary quantity p(r) such that σ(r) = p(r) / (1−p(r)), it can be shown (see Supplementary Note 1) using the electric field volume integral equation36 that
p(r)Einc = (1−p(r))Jp(r) + p(r) G0(r,r0)Jp(r0)dV0, (1)
where G0 is the dyadic free-space Green’s function, and regions r in the domain V with p(r)=1 correspond to metal and p(r)=0 correspond to free space. The solution Jp is unique and results in zero tangential electric field wherever there is metal, satisfying the PEC boundary conditions.

Equation (1) is strictly valid when the design comprises only metallic tiles and vacuum, as it uses the dyadic free-space Green’s function. However, G0 may be replaced with the Green’s function for any particular simulation environment, where additional materials representing arbitrary dielectric or metallic structures (e.g., substrates or feedlines) outside the optimization region are generally present. While closed-form analytical Green’s functions are usually not available, it is known that a Green’s function exists for every linear system, and a discrete numerical Green function matrix G can be obtained using a full-wave EM solver.

To solve for the current density numerically for a given design under consideration, by choosing a suitable basis to represent Jp (e.g., rooftop functions) and testing functions37, Eq. (1) can be discretized into
[(I−P) + PG]j = Cj = Peinc, (2)
where j and einc are the discretized polarization density and incident electric field vectors over the optimization region, and the design is encoded in P, a diagonal matrix whose entries indicate metal (1) or an empty location (0). The vector indices correspond to the degrees of freedom in the optimization region as discretized with the solution method of choice (e.g., spatial points on the Yee grid when a finite-difference method is utilized), with a total of Nopt such degrees of freedom. The matrix G, a discretized form of the Green’s function integral operator, needs only to be precomputed once for a given simulation environment, and then any candidate design may be evaluated by solving the linear system of Eq. (2) over only the optimization region, as illustrated in Fig. 2a. The number of unknowns Nopt is considerably smaller than the number of unknowns Nsim comprising the full simulation environment. This involves no approximations and incurs no loss of accuracy compared to a conventional full-wave simulation of the full system.

Precomputation
In general, an electromagnetic field solver finds the inverse of a matrix A that satisfies Ae = j, where the electric fields e and currents j encompass the entire simulation environment. The matrix A corresponds to the discretized Maxwell operator (e.g., the finite-difference frequency-domain (FDFD) matrix) of the simulation environment, comprising the static region and an empty optimization region (Supplementary Note 2). However, to obviate the need to compute the full A−1, a tall logical 0−1 projection matrix B may be defined to map vectors e and j in the optimization region to the corresponding vectors in the full simulation environment; that is, eopt = BTe and j = Bjopt, where BTB = I. The matrix BT extracts the fields from the appropriate indices in the full simulation result corresponding to the optimization region, and B maps the currents in the optimization region jopt into their corresponding locations in the larger j vector, filling the remaining entries with zeroes. From the definition of B, it follows that
eopt = BTA−1j = BTA−1Bjopt: (3)
The matrix BTA−1B corresponds to G in Eq. (2), as it maps currents jopt to fields eopt in the optimization region, while considering the static parts of the simulation outside of the optimization region via the full system matrix A.

An iterative solver may be used to obtain G column-by-column, where each simulation yields the fields due to a current density at a single discretized spatial location in the optimization region (i.e., each column of B). The Nopt simulations are completely independent of each other and, as such, may be run in parallel across many nodes. Additionally, time-domain methods such as finite-difference time-domain (FDTD) may be used, where the frequency-domain information in G is obtained from discrete Fourier transforms after each simulation. For multi-frequency optimization, this allows columns of multiple G matrices (one at each desired frequency) to be obtained from a single FDTD simulation, provided that the excitation has a sufficiently wide bandwidth to excite all frequencies of interest.

Alternatively, a sparse direct solver may be used with a frequency-domain formulation to obtain G efficiently in a single shot with the recently-introduced augmented partial factorization (APF) technique38. The full system matrix A is constructed using the FDFD formulation (see Supplementary Note 2). Then, an augmented sparse matrix K is set up such that A comprises the upper left block. K can be partially factorized as shown in Supplementary Fig. 1.
K = [A B; BT 0] = [KL11 0; KL21 I] [KU11 KL12; 0 KU22], (4)
where KL11 and KU11 are the LU factors and KL21 and KL12 are additional blocks not used for precomputation. The matrix KU22, known as the Schur complement39, is given by KU22 = −BTA−1B. The Schur complement has been previously applied to electromagnetic problems to reduce the size of the simulation domain, such as in domain decomposition methods, as with refs. 34,35. The G matrix is obtained as −KU22, avoiding the need to apply the LU factors for A to find BTA−1B.

Field quantities outside the optimization region are often required to evaluate the objective function. A vector xobj of quantities needed for evaluation may be defined, and a matrix Gobj that maps current densities in the optimization region to xobj may be precomputed and utilized in optimization together with G. Each of the Nobj elements of xobj is a field value or a linear combination of field values. In practical scenarios, Nobj is rarely significantly larger than 1; for example, to compute a scalar mode amplitude with a discrete mode overlap integral, xobj would have Nobj = 1 element that is a linear combination of the fields at the evaluation points40, and for an near-field to far-field transformation, two linear combinations are needed41, giving Nobj = 2. A wide logical 0−1 projection matrix WT may be defined to obtain xobj from the full simulation environment solution esim = A−1Bjopt:
xobj = WT A−1Bjopt + einc = Gobj jopt + xinc, (5)
where Gobj = WTA−1B and xinc = WTeinc. Using the above methods, Gobj may be precomputed together with G directly without additional computational cost. With an iterative solver, the number of excitations to be solved is still Nopt, since Gobj has Nopt columns. If APF is employed, the augmented system becomes
Kobj = [A B 0; BT 0 I; 0 WT 0],
and the Schur complement yields the vertical concatenation of G and the rows of WT.

It should be noted that although we have discussed obtaining the numerical Green function using FDTD and FDFD methods, any solution method of choice can be used, including finite-element and integral equation methods (e.g., Supplementary Note 7). Once precomputation has been performed for a given simulation environment and discretization scheme, the G and Gobj matrices may be used for any number of optimization runs, with no requisite prior training or no pre-existing libraries of simulated designs.

Direct binary search optimization
The DBS optimization algorithm starts with an initial design configuration P0, which may be randomly-generated or based on a priori design insight. The inverse of the initial system matrix, C−1 = [(I−P) + PG]−1, is found and stored, and the objective function is evaluated. At each iteration, a randomly chosen tile in the optimization region is flipped from empty to filled or vice versa. The objective function is evaluated, and if an improvement is obtained, the tile flip is retained for the next iteration. Otherwise, another random tile is flipped. Should all possible flips be tested without improvement, the optimization has converged. The DBS process utilizing PNGF is illustrated in Fig. 2b, and a flowchart is
xobj,n = Gobj jopt,n + xinc = Gobj C−1n-1 Pn einc + xinc (7)
must be found to evaluate the objective function. Although C−1n-1 = [(I−Pn) + PnG]−1 could be obtained by solving Eq. (2) from scratch, since DBS flips only a single tile per iteration, low-rank update methods can be used instead to avoid recomputing the inverse directly, as illustrated in Fig. 2c and discussed in the following section.

It is important to note that the update methods introduced in the following section are not restricted to DBS. Indeed, many optimization algorithms (e.g., level set methods, as discussed in Supplementary Note 6) modify a small (≪Nopt) number of unknowns in the design matrix P at each iteration of inverse design, and the low-rank update techniques are readily applied to such approaches.

Low-rank update evaluation
After a tile flip, the number M of modified elements in the diagonal design matrix P is, in the case where a finite-difference discretization is used, the number of field components comprising a tile on the Yee grid. Let dPn = Pn − Pn−1 represent the change to P. A wide logical 0−1 projection matrix Q may be constructed such that dPn = QTHP,n Q, where HP,n is a diagonal M-by-M matrix whose entries are the nonzero elements of dPn. Let U = QTHP,n and V = Q(G−I). Then, the update dCn = Cn − Cn−1 may be expressed as
dCn = dPn (G−I) = QTHP,n [Q(G−I)] = UV. (8)
Since C−1n-1 is the inverse of a low-rank (rank-M) update to Cn−1, the Woodbury matrix identity42 may be employed to obtain C−1n-1 through a rank-M update to the inverse of Cn−1. In computational electromagnetics, the Woodbury identity has seen use for efficient low-rank updates to a system matrix, particularly with MoM-based formulations. Prior work has employed the identity for applications including efficient computation of the gradient of an objective function via perturbations to the impedance matrix for optimization of metasurfaces43, photonic crystal design through low-rank updates to a Wannier basis analysis of the device44, and inverse design with a method of moments perturbation approach based on removing basis functions45.

Using the Woodbury identity,
C−1n-1 = (Cn−1 + UV)−1 = C−1n-1 − C−1n-1 U [I + VC−1n-1 U]−1 VC−1n-1. (9)
However, for many tile flips, the objective function will be worse than that of the previous iteration, and C−1n-1 would be discarded once found. Further performance may be obtained by instead finding xobj,n directly using C−1n-1. Substituting Eq. (9) into xobj,n = Gobj C−1n-1 Pn einc + xinc yields
xobj,n = Gobj C−1n-1 [I − C−1n-1 U [I + VC−1n-1 U]−1 V] Pn einc + xinc. (10)
Let
Rn−1 = Gobj C−1n-1, (11)
Sn−1 = C−1n-1 Pn−1 einc, (12)
and xobj,n-1 = Gobj C−1n-1 Pn−1 einc + xinc. (13)
Equation (10) becomes
xobj,n = Rn−1 dPn einc + xobj,n-1 − Rn−1 U [I + VC−1n-1 U]−1 V C−1n-1 dPn einc + Sn−1. (14)
Since Rn−1, Sn−1, and xobj,n-1 do not depend on the current candidate design Pn, they may be computed once at the start of a new iteration (after each successful tile flip) and used to rapidly evaluate xobj,n for new flips until the objective function is improved. Once this occurs, the tile flip is retained, C−1n-1 of the following iteration is obtained with Eq. (9), and Rn−1 and Sn−1 are updated. It can be shown (Supplementary Note 3) that computing xobj,n for a tile flip and updating C−1n-1 for a successful tile flip cost O(Nopt) and O(N2opt) operations, respectively, owing to the sparsity of dP and U.

Computational efficiency
The cost of objective function evaluation due to flipping a tile is completely independent of the size (Nsim) or complexity of the overall simulation domain and grows only linearly with respect to the number of unknowns (Nopt) inside the optimization region. Thus, for an optimization region of fixed size, the evaluation cost remains unchanged regardless of the simulation environment size. If a tile flip improves the objective function, the update to the system matrix using Eq. (9) requires O(N2opt) operations, which is significantly fewer than the O(N3opt) needed to invert C directly and also substantially smaller than the O(N3sim) operations required to invert the original full electromagnetic system.

Two examples to illustrate this performance are presented in Fig. 3. First, a planar optimization region of fixed size is considered on the surface of a dielectric substrate, and the evaluation time versus the simulation domain size (Nsim) is plotted for PNGF and compared with full-wave electromagnetic solvers as the substrate and simulation dimensions are increased. The second case keeps the simulation (Nsim) and substrate size fixed and plots the evaluation time versus the optimization region size (Nopt). The simulation environment is a 3D region with a finite dielectric substrate (λ0 = 10mm, εr = 3.5, thickness: 1.389mm) defined such that its sides are spaced a fixed distance (λ0/2) from Perfectly Matched Layer (PML) absorbing boundary layers41. The optimization domain is a square region centered on the substrate, and a ground plane covers the substrate’s bottom. Tiles (0.5 × 0.5mm, 3×3 Yee cells) populate the optimization region in a checkerboard pattern, a 2D lumped port (3×2 Yee cells) is defined in the center, and the objective function is the reflection coefficient. In comparison to FDFD (using APF as a solver), a custom FDTD solver, HFSS46, and CST47, PNGF achieves ultrafast (<100ms for all problem sizes tested) performance, faster than all other approaches by multiple orders of magnitude (10,900×–1,680,000× for the simulation sizes in Fig. 3b). Furthermore, PNGF provides an approximation-free solution, in common with the full-wave solvers considered, that is accurate for every possible optimization region configuration. For any given design, the PNGF results match with multiple digits of precision with those obtained by the solver used for precomputation (see Supplementary Note 8).

If multiple frequencies are of interest in optimization, a G matrix, with corresponding Gobj and C−10, may be precomputed for each frequency. The low-rank update evaluation procedure may then be applied to each system. Since each system is independent, finding each xobj after attempted tile flips and updating each C matrix after successful flips may be performed in parallel without any communication overhead.

It should be noted that, as long as G represents the simulation environment for the discretization scheme utilized for the optimization region and Gobj obtains the optimization region quantities of interest for evaluating the objective function, there are no limitations on the use of Eq. (14) for efficient updates to xobj,n. In particular, there is no upper limit on M. However, when M ~ Nopt, the evaluation of Eq. (14) may be slower than directly inverting C and using that to obtain xobj,n with Eq. (7).

Design studies
For each design study, PNGF is performed with two precomputation approaches: iterative, employing a custom GPU-accelerated FDTD solver, and direct, utilizing APF. The frequencies at which optimization is performed, and the objective functions for each case, are detailed in Supplementary Note 9. Simulations to verify the final designs are performed with HFSS and a custom FDTD solver. A comparison of the runtime performance of DBS optimization for each design study utilizing PNGF compared to full-wave solvers is shown in Table 1. The total times for PNGF represent the duration needed to design each device from scratch, with no requisite prior training or pre-existing libraries of simulated designs.

Progress in wireless systems, such as ultrawideband technologies, sub-terahertz/terahertz communications, and Internet of Things, has placed ever-increasing demands on antenna capabilities, performance, and size48–60. Planar antennas have been the subject of particular interest48–52,55,58 owing to their ease of fabrication and integration where space is limited. However, traditional topologies, such as patch antennas, are often narrowband or have strongly frequency-dependent radiation patterns. We design a broadband 30GHz center-fed substrate antenna with a wide fractional bandwidth and highly uniform pattern. The design is shown in Fig. 4 with simulated reflection coefficient and radiation patterns. The performance obtained at the final iteration of PNGF is also plotted, showing agreement with finite-difference simulation (see Supplementary Note 8). The design exhibits a 10dB return-loss bandwidth of approximately 13GHz, corresponding to a 40% fractional bandwidth. The radiation pattern remains largely unchanged over the entire frequency range, where the gain in the broadside direction is greater than 8.7dBi with a peak of 12.1dBi at 35GHz.

Advancements in cellular networks have placed ever-growing requirements on antennas for transmitting and receiving multi-directional, ultrawideband signals50,51,53–57. The multiple-input multiple-output functionality of current cellular technology is typically realized using phased arrays or multiple antenna elements51,54–56, whose large electrical size restricts miniaturization. As such, reconfigurable antennas, whose properties may be altered dynamically with inputs (e.g., switches), have garnered substantial attention57–60. We design a 30GHz SBA for 5G applications, with a switch for selecting between two target beam directions (θ = 45°, ϕ = 90° with switch open; θ = 45°, ϕ = 270° with switch closed). Due to practical equipment limitations, we scaled the inverse-designed antenna up in all dimensions by a factor of 5× before fabrication to shift the center frequency to 6GHz and facilitate measurement. The fabricated SBA is shown in Fig. 5 with the simulated and measured reflection coefficient and radiation patterns versus θ for ϕ = 90°. The simulations are performed with the fabricated upscaled design, and the measurements agree closely, with a 10dB return-loss bandwidth of 0.3GHz (switch closed) and simulated peak gains of 8.2dBi (switch open) and 10.6dBi (closed). The measured angle of beam switching when viewed in the yz plane is approximately 70°.

SIWs comprise a planar substrate enclosed by metal cladding and sidewalls formed by vias. Owing to their compatibility with printed circuit board fabrication processes, SIWs have attracted considerable interest61–65. The fundamental mode is typically excited with a tapered transition from a microstrip feed to the SIW. For compactness, it is desirable to decrease the length of the transition, but this often yields decreased performance. We design a taperless transition from a 50Ω-impedance transmission line to a broadband SIW. The length of the optimization region is more than 4× shorter than the length of a linear taper required to achieve comparable bandwidth, as in ref. 65. Optimization is performed to minimize the insertion loss over the bandwidth of interest. The fabricated waveguide section with transitions and the simulated and measured S11 and S21 are shown in Fig. 6, demonstrating a wide 10dB return-loss bandwidth of approximately 7.7GHz. As with the preceding SBA design, minor differences between the simulated (with feed connector) and the measured results are visible. The discrepancy likely arises from uncertainties and variations in fabrication and material properties, such as permittivity, substrate losses, substrate thickness, and surface roughness of the metallization, as well as imperfections associated with the measurement setup and instrumentation.

Methods
Finite-difference discretization
For each test case in Fig. 3 and each design study, the same Yee lattice is used for both FDTD and APF precomputations as well as the optimization with PNGF. The optimization regions are rectangular areas made up of faces of adjoining Yee cells. Equivalent polarization density components, as discussed in the “Results” section (Precomputation), in x and y are defined on these faces. For example, with tiles comprising 3×3 Yee cell faces as illustrated in Fig. 1a, each tile comprises 12 jx and 12 jy components, and the number of modified elements M in a tile flip is 24. The current components are co-located with the electric field x and y components on the edges of the Yee cells. Such flat optimization regions are appropriate for modeling the copper clad utilized in these examples, whose thickness is much smaller than the wavelength. For other applications, however, the optimization region may encompass multiple layers of Yee cells, within which z-components of the current density would also be present. The simulation environment for each design is truncated using Perfectly Matched Layer absorbing boundaries.

Computational resources
A custom solver that constructs the FDFD system and performs APF is used to generate the G and Gobj matrices. The MUMPS package66 is used to carry out partial factorization and compute the Schur complement. For large simulation environments, it may be infeasible to use a direct solver owing to the required amount of random access memory, and iterative solvers must be employed, incurring the cost of running Nopt full-wave simulations. However, this is not the case for the three design studies considered in this work, and APF achieves up to 12× decreases in execution time compared to the GPU-accelerated FDTD for precomputation.

For each design study, APF precomputations are run on three 128-core AMD EPYC 7763 nodes, where each precomputation uses 64 cores on one node. These are run in parallel, one for each frequency of optimization in each design. For GPU-accelerated FDTD precomputations, 24 FDTD simulations are run in parallel using the EPYC 7763 nodes with one Nvidia A100-SXM4-80GB GPU per simulation. For timing comparison, HFSS and CST are run with 128 cores on one node, and the PNGF implementation utilizes the BLAS and LAPACK linear algebra packages with the same number of CPU cores. As optimization is performed at multiple frequencies for each case, evaluation with PNGF is performed sequentially for each frequency during each iteration. For HFSS, the HFSS solution type (finite-element method-based) is used with a multi-frequency solution at each optimization frequency with the direct solver with first-order basis functions and adaptive meshing. For CST, the time-domain solver (finite integration technique-based) with a hexahedral mesh, which is similar to FDTD, is used without adaptive meshing.

For the objective function evaluation benchmarking of Fig. 3, HFSS and CST are utilized as solvers using the same resources as above, respectively. The FDTD solver is a custom multithreaded implementation using OpenMP and 128 cores on a single node. FDFD is performed with APF as a solver, using 128 cores on a single node. PNGF is used for objective function evaluation for each case with a single core and also with 128 cores on a single node.

Optimization parameters
The substrate antenna design utilizes a 1.39mm-thick substrate (εr = 3.5 representing Rogers RO3035) cladded with 13.9μm-thick copper. The bottom copper layer is fully filled as a metal ground reflector, and the optimization region is defined on the top layer. The optimization region comprises a 21×21 grid of tiles, where each tile is 3×3 Yee cells of 0.5×0.5mm each. This results in 4032 ex/jx and 4032 ey/jy components. In view of maximizing the gain in the broadside direction, x and y symmetry are enforced; as such, only 2048 simulations are required when precomputations are performed with FDTD, and 4 tiles are flipped at a time during optimization. A 50Ω x-directed lumped port in the center is used as the excitation source for field simulations.

The 30GHz SBA design utilizes a 0.508mm-thick Rogers TC350 (εr = 3.5) cladded with 18μm-thick copper. The bottom layer comprises the ground plane, whereas the top design domain is approximately one wavelength square with 21×20 tiles each comprising 3×3 Yee cells of 0.5×0.5mm each, giving 3843 ex/jx and 3840 ey/jy components. With axial symmetry, the number of simulations necessary for precomputations with FDTD is 3872, and 2 tiles are flipped per optimization iteration. The switch is modeled as an ideal metallic via connecting the top and the bottom. The antenna is edge-fed with a microstrip feed; for field simulations, a z-directed 50Ω lumped port is attached between the feed and the bottom ground. The objective function is set to minimize the reflection coefficient and increase the directivity in the directions of interest for each configuration of the switch.

The microstrip-SIW transition design utilizes a 0.508mm-thick Rogers RT/duroid 5880 substrate (εr = 2.2) with 35μm-thick copper clad. The bottom layer is completely filled with copper as a ground plane. Each of the two optimization regions, which are constrained to be identical, comprises 52×13 tiles. Since the fundamental mode and the structure are longitudinally symmetric, the optimization region can be reduced in half, to 26×13. Each tile comprises 3×3 Yee cells of 0.125×0.125mm each, corresponding to a total of 6240 ex/jx and 6123 ey/jy components in the optimization region, and 6201 simulations are required for FDTD precomputations. In field simulations, a z-directed 50Ω lumped port is attached to each microstrip feed end, connecting the ground plane (bottom layer) to the microstrip (top). To design a broadband device, the objective function is set to minimize the insertion loss at five different frequencies.

Device fabrication
The scaled 6GHz SBA is fabricated from a 2.5mm-thick Rogers TC350 substrate with ε = 3.5 cladded with 1oz copper. To implement the reconfiguration switch of the fabricated SBA, two antennas are fabricated, which differ only in whether the switch via is present (switch closed) or absent (open). A 2.92mm end-launch RF connector (Withwave SM03FS017) is soldered to pads at the end of the microstrip feed of each antenna to provide excitation.

The SIW section is fabricated from 0.508mm thick RT/duroid 5880 laminate cladded with 1oz copper. Two 2.92mm end-launch RF connectors (Withwave SM03FS007) are soldered at both microstrip feed ends for the S11 and S21 measurements.

Measurement system
To measure the reflection coefficient and radiation pattern of the SBA, a measurement setup is established in an anechoic chamber with a vector network analyzer (VNA) (Keysight N5247). For pattern measurements, the SBA is affixed to a two-axis rotary positioner (Diamond Engineering DCP252) driven by stepper motors, and an excitation horn antenna (Com-Power AH-118) is positioned facing the SBA. The pattern is obtained by recording transmission coefficient data with the VNA connected to both antennas, while the elevation and azimuth are swept.

The S11 and S21 measurements of the microstrip-to-SIW transition are obtained with a VNA (Rohde & Schwarz ZVA50), with each microstrip connection attached to a VNA port.

Data availability
Source data are provided with this paper.

Code availability
The source code for the PNGF method is publicly available at: https://github.com/ACME-Lab-Stanford/PNGF67.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this paper, the authors present a methodology for the rapid inverse design of electromagnetic structures. The proposed technique exploits the use of Green's functions and utilizes certain matrix handling schemes to achieve its objectives. The analysis is applied to three examples involving passive pixelated structures, two of which are fabricated and measured; something positive for the credibility and and verification of the results. The paper is adequately organized, while the supplementary material is proven to be helpful particularly, the provided flow chart. Finally, the language and the grammatical structure are very satisfactory without any subtle points or vague parts.\nHaving carefully read both the manuscript and the supplementary material my most important comments can be summarized as follows:\n1. The first concern is the fact that the contributions of this paper would be best suited to a more \"engineering\"-based journal. In fact, the novelties of the specific work are basically in the area of computational electromagnetics and not in the physics of the problem or the description of a new physical phenomenon, which could be simulated by the proposed technique. This can be also derived by the majority of the actually related references of the paper [coming basically from IEEE Transactions]. Thus, it is important for the authors to re-express their contribution and clearly attest on its suitability with the specific journal in terms of its contribution. Perhaps, structure of the Introduction and the presentation of the methodology could be described in a different way which would be more persuasive on its suitability.\n2. The second point refers to the applicability of the proposed technique. Although the authors mention in several places [and mostly in the title] that the their method can be applied to an assortment of problems, they mainly focus on the study of pixelated structures. I do understand that such examples can be more easily fabricated and tested, yet some additional examples with more complicated applications could be included. Specifically, what would be the Green's functions in the case of unit whose axes were not aligned to the axis of the coordinate systems or contained curved parts? What would be the efficiency of the proposed technique and the most important how fast would be the inverse design? Perhaps these structures could be analyzed via FEM.\n3. Following the above comment: The paper deals with PEC pixelated unit cells; It would be nice to include an example with some dielectric structures and provide the corresponding Green's functions along with the achieved acceleration.\n4. Since, nowadays, the analysis of smart metasurfaces or even RIS devices is rather significant, it would be interesting if the authors provided an example of a pixelated structure, yet with active elements [such as those used for a reconfigurable beam-forming procedure or the control of anomalous reflection from the metasuface].\n5. The title of Section 2 [i.e. Results] is somewhat misleading, since this the section where the authors analyze and describe the main methodology of their work. Furthermore, the limitations [if any] in the use of eq. (14) should be provided. Are there any scenarios or configurations that could not be treated via this expression or the related scheme?\n6. As a final recommendation: Please try to provide the solution of dome of your problems via the FEM or the integral equation method, since there are some important differences in the handling of a problem via a time-domain or a frequency-domain method that could affect the solution speed or even lead to loss of convergence.\nA more detailed README file would be appreciated from the interested readers.", "reviewer_2": "Summary: The paper presents a promising framework for accelerating the inverse design of pixelated metallic electromagnetic structures. Such structures are widely used in antennas, radio-frequency integrated circuits, etc. The authors introduce the precomputed numerical Green functon(PNGF) method, which precomputes the interaction between a static subdomain and a modifiable design subdomain. The work then performs iterative design updates using direct binary search (DBS) combined with a low-rank matrix update strategy based on the Woodbury identity. The work is well-motivated and includes measured validations. However, the reviewer recommends a major revision before the paper can be considered for publication, for the following key reasons:\n1. Clarification on Superposition and Equivalent Currents\nIn both Figure 1 and the technical writing in Sec. 2.1, the authors state that equivalent polarization currents obey superposition, in contrast to the original metallic tiles. It may be misleading without clarification.\nIf the equivalent currents j_1 and j_2 are obtained from two configurations, simply summing equivalent currents does not yield the correct solution for the new configuration, as shown in Figure 1(c). This is because the total field must satisfy PEC boundary conditions at both metallic tile locations, and the original equivalent current is not designed for the new setting.\nSuperposition holds only when the equivalent currents are solved from a system that simultaneously enforces all boundary conditions. This distinction is subtle but important and should be clearly stated to avoid confusion.\n2. Simplification of the Writing and Notations\nThe current presentation introduces a number of matrix symbols, G, A, B, C, P, R, U, V, etc, which are not easy for the readers to follow. The authors may consider grouping or renaming matrices to highlight their physical roles.\nThe authors could also consider introducing a simplified derivation based on the free-space Green’s function method, and then discuss the numerical Green’s function case.\n3. Suggestions for the Numerical Comparisons\nTable 1 compares the proposed PNGF method (with low-rank updates and APF) against commercial solvers (HFSS, CST) in terms of runtime. However, the benchmark is incomplete and potentially unfair.\nThe design problems all involve pixelated PEC structures, which are naturally suitable for the FDTD solvers. There is no baseline provided for a conventional FDTD + DBS implementation, which is a fairer comparison to highlight the benefit of APF and low-rank updates in the PNGF approach. The reviewer strongly encourages the authors to include an FDTD-based DBS implementation as a runtime and accuracy baseline. This would strengthen the numerical studies and comparisons.\n4. Suggested Title Refinement\nGiven the current scope of the paper, which focuses exclusively on binary optimization of pixelated PEC structures, a more accurate and informative title would be: Ultrafast Inverse Design of Pixelated Metallic Electromagnetic Structures\nThe authors have made the code available. The reviewer was able to download and review it. The package includes the core implementation of the PNGF framework, the scripts, and data to reproduce the numerical experiments presented in the manuscript. A readme file is included and provides basic instructions for running the code.\nSuggestions: The package could be a usable resource for the community with more detailed usage instructions, particularly on how to set up a new problem or new objective functions.", "reviewer_3": "Sun et al reports in “Ultrafast Inverse Design of Electromagnetic Devices”, a way to utilize precomputed simulated values to avoid having to simulate the whole electromagnetic structures, leading to shorter time in simulation of devices, as the effective simulations region made smaller to regions that change. The authors claim 16000-fold reduction in speedups in their demonstrations and claim in the report that the proposed method is relevant for general inverse design method where electromagnetic devices are optimized. I find the benefit of such speed reduction in computational simulation to be clear, and the authors find right problem space to demonstrate such simulation speed up to verify its functionalities. However, there are some places where I find this work to be somewhat misleading due to its broad claim and presentation in the manuscript.\nFirst of all, authors lack explicit accuracy validation, which is critical aspect of this work. The authors repeatedly claim that the PNGF method incurs “no approximations or compromises in accuracy.”. However the paper does not provide any explicit benchmarking of the solution accuracy. There is no clear quantitative comparison between results obtained via the PNGF approach and those from a standard full-wave solver on the same design, to verify they match exactly. Such a validation is critical for the community to see that such speedups do not come at the cost of errors. I strongly recommend including an accuracy benchmark.\nAuthors should consider discussing more about the generalization beyond PEC structure. All the examples and discussions in the paper focus on cases with perfect electric conductor (PEC) boundaries (metal antennas, waveguide walls) and a fixed background medium. It remains unclear how well the proposed method generalizes to more complex material scenarios, such as arbitrary permittivity distributions or dispersive materials. The PNGF approach relies on precomputing a Green’s function for the “static” part of the environment. If the design involves altering dielectric regions or materials with dispersion as a function of the frequency, one would need to recompute the Green’s function or apply additional adjustments, which could reduce the simulation speed-up efficiency gains. The manuscript would benefit from a discussion of this point. Are there fundamental limitations on using PNGF for, say, optimizing dielectric meta-optics or multi-material devices? If the method is principally demonstrated for PEC inclusions in a homogeneous or static dielectric host, the authors should clarify this scope limitation beyond what they have provided in the supplementary information. Or authors could consider more explicit acknowledgment that extension to those cases is non-trivial and left for future work. This is important so readers do not over-generalize the method’s applicability beyond what has been shown/proven.\nThe two main design examples (the switched-beam antenna in Fig. 5 and the broadband transition in Fig. 6) show experimental results, including simulated vs. measured performance. First, I would suggest authors to plot those measured /simulated data in linear scale instead of log only. This would provide readers to compare the accuracy of the simulation and how measurements are turning out.\nSecond, it is unclear if the authors are employing their PNGF method in those simulation plots, i.e. it is difficult to know how simulation of these structures are performing with respect to the PNGF simulation method. Did the authors use the PNGF solver to perform the entire optimization for Fig. 5 and Fig. 6, and if so, how long did it take and what speedup was achieved versus a normal full-wave approach? Similar to above accuracy bench marking point, authors should compare simulation result of PNGF to the results of other simulation method of this device and compare against the experimental result.\nLastly, the idea of reusing computation (Green’s functions, system matrices, etc.) to speed up EM simulations is not entirely new. The authors should acknowledge relevant prior approaches in the literature. For instance, citing Ludick et al. (IEEE Trans. Antennas Propag., 2014), who introduced a Domain Green’s Function Method to efficiently analyze large antenna arrays by reusing precomputed interactions, and Budhu & Grbic (IEEE Trans. Antennas Propag., 2023), who applied the Woodbury matrix identity to rapidly optimize metasurfaces in a Method-of-Moments framework would be a good start. These works are closely related in spirit – Ludick et al. demonstrate the reuse of domain-specific Green’s functions to reduce solve times, and Budhu & Grbic achieve speedups in iterative design by updating an inverted system matrix rather than recomputing it from scratch. Citing these will also clarify to readers that the authors are aware of existing relevant acceleration techniques and how their approach differs. For instance, the authors could consider explaining how PNGF is conceptually different from domain decomposition or how the specific application of the Woodbury identity here is novel in comparison to the prior art.", "reviewer_4": "In this paper, the electric field volume integral equation and the definition of an optimization region much smaller than the full simulation region (many less variables overall) are used to reduce the computational cost of direct binary search global optimization for pixelated metallic structures.\nI provide some comments for the Authors in the text below:\n1) Reading the title of the article leads to think of a very general methodology for fast inverse design for electromagnetic structures. However, the proposed algorithm uses multiple assumptions and is validated for specific structures, namely pixelated metallic structure. I think the title should be more specific and tailored to the actual main contributions of this work.\n2) The acronym PNGF is not defined in the Introduction.\n3) Regarding the sentence \"In particular, the fields generated by such equivalent polarization current densities do obey superposition, whereas the original metallic tiles do not.\". I think it can be important for the readers to further elaborate on this aspect, which is fundamental in this work.\n4) A comparison (Table 1) with commercial solvers is provided, namely Ansys HFSS and CST Microwave Studio. Could the Authors briefly describe the EM methods used in these solvers and clarify the main internal software settings used for the simulations?Among many possible settings, is adaptive meshing used?\n5) Could the Authors elaborate in detail on the expected performances of the proposed algorithm for other cases where the optimization region is not so much smaller than the simulation region?As one possible example, let us imagine a multilayer metasurface with potentially different dielectric materials per layer, where each voxel can be an optimization variable, an excitation wave comes from the substrate and another plane to record electromagnetic fields is located on top of the structure.\n6) Could the Authors elaborate in detail on how the complexity of the proposed algorithm and steps like the low-rank update change when other optimization algorithm as genetic optimization are used where multiple tiles can change per iteration without a specific defined pattern?\n7) Table 1, the Authors specify that “The speedup figure compares only the optimization time, since precomputation for the PNGF method needs only to be performed once for a given simulation environment, and the precomputed Green function matrices may be reused for any subsequent optimization runs. The smaller of the HFSS and CST times is used for each case.” However, if we think about an absolute speed-up measure I think it might be more appropriate to include all CPU time costs needed for the proposed method and re-evaluate the speed-up figures of merit.\n8) Are specific parallel computing internal software settings used for the HFSS and CST simulations on 128 cores?Are the results in Table 1 based on the use of 128 cores for all compared techniques?\n9) Looking at the differences between measured and simulated data as in Fig. 6, is there a systematic complexity for the computational electromagnetic methods to handle layout structures made of pixels with respect to more slowly varying classical parametric designs?\n10) Could the Authors improve the readability of figures as Figs. 4c, 5c, 5d, 6b?\n11) How does the CST Microwave studio solver perform considering the results in Figs. 4b, 5c, 5d, 6b?\n12) How do the electromagnetic results obtained by the PNGF method compare with HFSS and CST regarding the results in Figs. 4b, 5c, 5d, 6b?", "reviewer_5": "The manuscript aims at the large amount of electromagnetic simulation required in the inverse design of electromagnetic devices and proposes a fast reduced-order method to address the issue of long computation time in traditional numerical modeling methods. The proposed method first decomposes the computational domain into a small variable region and a large invariant region. Then, it leverages the Schur complement operation to construct a numerical Green’s function for eliminating the large invariant region. In this way, the simulation of the original large-scale domain is transformed into the simulation of the local small region that needs to be optimized. More importantly, the authors claim that local simulations no longer necessitate the solution of matrix equations. Experimental results demonstrate its superior computational efficiency compared to commercial software, highlighting its potential applicability. However, the manuscript requires further refinement in theoretical rigor, methodological clarity, and structural logic. Specific concerns are detailed below:\n1.The structure of the manuscript is not very reasonable. Chapter 2 is titled Results and Chapter 4 is titled Methods. But Chapter 2 also introduces a lot of methods that don't match the title. It is suggested to introduce the method completely first, and then describe the experiment design and results.\n2.While the core innovation lies in the electromagnetic modeling method, the manuscript presents only conclusion-oriented formulas (e.g., Eqs. 1 and 3) without detailed derivations. The readers need to refer to external references to judge the correctness of the formula. It is suggested to add a detailed derivation to illustrate formulas (1) and (3).\n3.The manuscript applies FDFD to generate matrix A in formula (3), and then inverts it as a numerical Green function to replace the Green function in formula (1). As far as I know, the matrix A generated by FDFD method contains the material information. Matrix A is different after material changes in the computing domain, so A should be recalculated and inverse after each optimization. However, from the introduction of this paper, A does not need to be recalculated, please explain the reasons in detail.\n4.The manuscript mentioned that FDTD method can be used instead of FDFD method to calculate numerical Green function. As far as I knows, the resulting matrix of FDTD is implicit, and the electromagnetic field in time and space is deduced by leapfrog algorithm. How to get the numerical Green function displayed?\n5.The current title, “Electromagnetic Devices,” is overly broad. As the method targets pixelated metallic structures, the title should be revised to reflect this specificity.\n6.I suggest that the author make major revision to the manuscript. Focusing only on the FDFD method, the derivation of numerical Green function is given in detail, with emphasis on explaining why the local matrix is no longer computed after optimization."}, "author_response": "We thank the editors and reviewers for their insightful remarks and valuable feedback, which have helped us significantly improve the quality of our paper. Our responses to the remarks and questions are below. In addition to the modifications to the manuscript that we have made according to the reviewer feedback, we have made several additional improvements to the paper, which are documented here: (1) Several symbols in the original submission were used to denote multiple different quantities, such as U (for both the LU-factor in Eq. 4 and the U matrix in the PNGF update equations). We have renamed these symbols to avoid as much re-use as possible. Please refer to the new Supplementary Table ST.1 for a list of symbols and their definitions to benefit readers’ clarity in following the derivations in the text. (2) In the original submission, we had accidentally mislabeled the S simulation results of 11 the switched-beam antenna from HFSS without the connector in Fig. 5(b) as FDTD data. This has been corrected in our revision. (3) We have found a few minor errors in Supplementary Note SN.3 (“Cost of system matrix update and objective function evaluation”), in the derivations of the complexity of computing the objective function and performing the low-rank update. The number of multiplication operations required to accomplish each of these can be decreased slightly by using the sparsity of the U matrix (for objective function evaluation) and by reusing the quantity C-1 U (for C-1 update). n-1 These changes do not affect the conclusion of the Supplementary Note; that is, the complexity of objective function evaluation is O(N ) and that of updating the system opt matrix inverse is O(N2 ) with respect to the size of the optimization region N . opt opt [Detailed point-by-point responses follow...]"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have made a serious effort and succeeded in elevating the quality of their work.", "reviewer_2": "The authors have successfully addressed most of my comments and concerns. As a minor correction, CST Studio Suite’s time-domain solver is based on the Finite Integration Technique (FIT), not the Yee-grid FDTD; please adjust the wording. In addition, please consider citing the following related prior work on Platform Green’s Functions for electrically large PEC platforms with local design variations: H.-W. Gao, S. Wang, X.-Q. Sheng and Z. Peng, \"Rapid Numerical Analysis of Electrically Large PEC Platforms With Local Variations via a Platform Green’s Function Method,\" IEEE Trans. Antennas and Propagation, vol. 70, no. 10, pp. 9544-9556, Oct. 2022, doi: 10.1109/TAP.2022.3181016. This work follows a similar offline-online methodology: a platform-aware integral operator is precomputed via Schur-complement reduction and skeletonization, so the online cost depends mainly on the local design region. A brief mention in the Introduction/Related Work would help situate the contribution.", "reviewer_4": "Thank you for the replies to the comments. I have some remaining comments that I list below:\n1) In Figs. 5 and 6, differences between simulations and measurements are visible. The reply to my comment 9) (in the rebuttal letter) should be included in the manuscript. It is valuable information to the readers.\n2) I think it is also valuable to include the figures used in the reply to my comment 11) (rebuttal letter) into the supplementary notes;\n3) I also encourage the Authors to include the reply to comment 5) of Reviewer #1 (about potential limitations of the proposed technique) into the main manuscript.", "reviewer_5": "The author has revised the paper according to the reviewers' suggestions, but I still have the following comments:\n1. The author still chooses a title with a relatively broad scope, emphasizing that the proposed PNGF method is applicable to medium optimization problems. In my opinion, it is necessary to present a numerical example related to medium optimization, rather than just textual discussion.\n2. The author compares the performance of the PNGF method with those of HFSS and CST, demonstrating the computational speedup of the proposed method. However, since the methods used in HFSS and CST differ from the FDFD algorithm adopted by the author, I believe objectivity is lacking. I suggest that the author compare the proposed fast algorithm with a traditional scheme that combines FDFD with an optimization algorithm to calculate the overall objective. When the author uses FDFD to compute PNGF multiple times, it is necessary to calculate the background environment of a large region, and the computational load is not much less than that of FDFD calculating the overall objective including optimization and background. If the number of optimization iterations is not far greater than the number of PNGF calculations, the time saved by the proposed method compared to the traditional FDFD method will not be so significant.\nAdditional comments regarding the responses to Reviewer #3:\nThe authors claim that they proposed PNGF method can be used with any suitable discretization approach, including finite element methods (FEM) and method of moments (MoM) integral equations. In SN.7, they present a proof-of-concept example of utilizing the PNGF method with a MoM solver to rapidly optimize a multilayer RF beamforming metasurface. However, even though I am proficient in MoM, I still do not understand the methodology of this PGNF based on surface integral equations. I have the following questions for the authors to address:\n1, From equation (S43), the numerical Green’s function matrix \\mathbf(G) is denoted as \\mathbf(A)^O-\n\\mathbf(A)^{OS}\\mathbf(A)^{S}^{-1} \\mathbf(A)^{SO}. To my knowledge, matrices \\mathbf(A)^O, \\mathbf(A)^{OS} and \\mathbf(A)^{SO} change along with the variations of the metal patterns on the upper and lower layers. Thus, the numerical Green’s function matrix \\mathbf(G) changes along with the optimization scheme of the metal pattern. Why can the matrix \\mathbf(G) be pre-calculated?\n2, To my knowledge, the calculation of the far-field characteristics of the metasurfce requires the electric and magnetic currents on all regions, including the optimized region and the background region. After solving Eq. (S43), only the electric current on the optimized region can be obtained and the electromagnetic current on the background region remains unknown. If you want to obtain the far-field characteristics of a certain optimization scheme, you need an additional process to obtain the electromagnetic current on the background region. How was this achieved, and was the required time accounted for the current numerical examples?\n3. I have doubts about such a high acceleration rate of the IE based PNGF method comparing to the conventional IE method. Does FEKO use an iterative method or a direct method to solve Eq. (S42)? Does you use an iterative method or a direct method to compute \\mathbf(A)^{S}^{-1} and solve Eq. (S43)? Employing a direct method, as opposed to FEKO's iterative method, will undoubtedly reduce the computation time. However, it will also require a substantially larger memory. The performance of a computational method must be assessed not only in terms of time but also memory usage. Please clarify the solver's applicability and suggest adding a comparison of memory consumption."}, "author_response": "We thank the editors and reviewers for their helpful feedback and comments on our revised manuscript. Our responses to the remarks and questions are below. [Detailed point-by-point responses follow...]"}, {"round": 3, "reviewer_comments": {"reviewer_5": "The authors have addressed most of my comments and concerns."}, "author_response": ""}]

Decision letter