Surrogate-Assisted Design of Checkerboard Metasurface for Broadband Radar Cross-Section Reduction

Metasurfaces have been extensively exploited in stealth applications to reduce radar cross section (RCS). They rely on the manipulation of backward scattering of electromagnetic (EM) waves into various oblique angles. However, arbitrary control of the scattering properties poses a significant challenge as a design task. Yet it is a principal requirement for making RCS reduction possible. This article introduces a surrogate-based approach for rapid design optimization of checkerboard metasurfaces. Our methodology involves fast metamodels, and a combination of surrogate-assisted global optimization with local, gradient-based tuning. It permits an efficient control of the EM wave reflection characteristics, and ensures arriving at that the globally optimum solution within the assumed parameter space. The design procedure is fully automated. The framework is employed to develop a novel broadband checkerboard metasurface, where the RCS reduction is fundamentally based on the backward scattering manipulation carefully controlled by simultaneous adjustment of the unit cell dimensions. The properties of the structure are demonstrated using simulated monostatic and bistatic RCSs. The proposed metasurface exhibits 6 dB RCS reduction within the frequency range from 16 to 37 GHz. The numerical results are validated using physical measurements of the fabricated prototype. Experimental data indicates that the relative RCS reduction bandwidth is 83 percent, which makes the proposed structure outperforming the designs reported in the literature.


I. INTRODUCTION
The advancements in the field of metamaterial technology have opened the new paths to numerous applications, such as invisibility cloaks, gradient index lenses, polarization converters, holograms, unique antenna designs, and many others [1]- [4]. Metasurfaces, two-dimensional equivalents of metamaterials, are planar patterned surfaces composed of subwavelength periodic arrays of unit cells [5]. Owing to their extraordinary capability of manipulating the scattering behavior of the electromagnetic (EM) waves, the popularity of metasurfaces has been steadily increasing in the field of stealth technology [6]. Therein, the primary concern is to reduce the radar cross-section (RCS) to evade from the enemy's radar, which can be achieved by diminishing The associate editor coordinating the review of this manuscript and approving it for publication was Jenny Mahoney. back-scattered EM waves from the metallic objects [7]. The four leading practical approaches extensively used in the literature to achieve RCS reduction include [8], [9]: (i) utilization of radar absorbing materials (RAM), which transforms the incident EM wave into heat; (ii) reshaping the geometry of a target to redirect the incident EM energy away from the source; (iii) redirecting (or deflecting) the incident EM wave around the object (invisibility cloaking); (iv) phase cancellation, both active and passive. However, all of the aforementioned approaches predominantly exhibit narrow RCS reduction bandwidth, suffer from design complexity, and extreme losses.
Quite recently, considerable interest emerged in utilizing metamaterials for wideband RCS reduction. On a generic level, there are two strategies for reducing RCS by means of metamaterials. The first one is the usage of a perfect metamaterial absorber [10]- [14]. Such materials can absorb EM waves and convert the energy into heat. Nevertheless, the RCS reduction band remains limited. The second strategy is to exploit the reflection phase controlling property of metasurfaces. Two types of surfaces have been presented that capitalize on this concept, i.e., electromagnetic gradient surface (EGS) [15], and checkerboard metasurface [6]. In EGS, the metal part of the surface is replaced by the unit cells of artificial magnetic conductors (AMC), and perfect electric conductors (PEC). The primary requirement in EGS is to maintain equal phase difference between the unit cells [16]. When the plane wave is incident from the normal direction, the EGS reflects back the tilted beam pattern, hence reducing the RCS. Due to non-linear relationship between the reflection phase curves and frequency, it is difficult to meet the equal phase difference condition over a wide frequency range. In a checkerboard metasurface, AMCs and PECs are arranged in an alternate fashion. The idea is to keep 180 • phase difference between the AMC and PEC unit cells. Such a combination successfully diffuses the scattering energy at four lobes in the diagonal plane [17]. The EGS and checkerboard metasurfaces are low profile, robust and simple to manufacture [18]. Their major drawback is the narrowband performance of the AMC structure. Outside the working bandwidth, the AMC properties are similar to those of PEC, and the condition for 180 • phase difference no longer holds. To overcome this drawback, PEC unit cell is substituted by another AMC unit cell operating at a different resonant frequency. Consequently, a dual-band design can be obtained [19], [20]. The idea of employing two AMCs in a checkerboard configuration was originally presented and developed by de Cos et al. [21], [22]. To achieve RCS reduction over a broad frequency band using this configuration, the phase difference between the two AMC unit cells should be 180 • when their reflection amplitudes are the same and equal to one [18], [23]. In terms of electrical characteristics, the phase reflection curves of the two unit cells should remain parallel (i.e., equidistant) over the frequency band of interest. Notwithstanding, the reflection amplitudes of the combined unit cells are not always the same due to losses. On the other hand, it has been shown that −10 dB RCS reduction can be maintained over a frequency band if the phase difference between the two unit cells remains within 180 • ± 37 • range [25]. In a related vein, the concepts of coding metasurfaces [26], [27], diffusion metasurfaces [28], [29], programmable metasurfaces [30], Huygens' metasurfaces [31], as well as cloaking structures [32], have been proposed, which offer a control over the wavefront in a more sophisticated manner. The primary advantage of coding and diffusion metasurfaces over the checkerboard type surfaces is that it scatters the incident EM waves into all directions. In addition to that, coding metasurfaces are also exploited as an absorptive surface to realize essential RCS reduction [33].
Until now, numerous novel designs have been proposed for attaining wideband RCS reduction using metasurfaces [18], [21]- [29]. In the absence of reliable analytical methods, the design process in the above-mentioned works typically relies on iterative full-wave EM simulations. Although such methods ensure accurate evaluation of the system response, they are time consuming and laborious due to a considerable amount of designer's interaction involved in the process. Furthermore, the design procedure relies mostly on empirical reasoning, physical intuition, or trial-and-error, which raises questions about the reliability and efficacy of such methods, as well as their capability of identifying truly optimum designs. From the perspective of hands-on design procedures, the problem is additionally aggravated by highly nonlinear input-output relationships. New and more sophisticated methods should be conceived to make the design process of metasurfaces computationally efficient, robust, and automated. In the recent years, data-driven techniques have emerged as promising tools, applicable to solving problems in many areas of science and engineering. Their advantages include the ability to yield acceptable solutions under time constraints and limited computational resources [34]- [41]. Some of the recent alternative approaches include phylogram analysis-based optimization method [42], island-based cuckoo search with polynomial mutation [43], hybrid swarm algorithm (a combination of the strengths in self-assembly and the particle swarm optimization) [44], and grey wolf optimizer-based method to tune pi-fuzzy controllers [45]. However, this work adopts some specific methods such as surrogate modeling frameworks and global optimization routines as the components of the developed metasurface design procedure.
The main objective of this paper is to enhance the RCS reduction bandwidth along with addressing the key challenges at the design level of a metasurface. The considered metasurface architectures are periodic arrays of two different AMC unit cells on the same ground plane in a checkerboard configuration. A surrogate-based framework proposed in this work involves fast kriging metamodels as well as a surrogateassisted global search algorithm. The metamodels are trained using sampled EM simulation data, and used as the unit cell phase characteristic predictors at the optimization stage. Our procedure allows for identifying the optimum geometries of the individual unit cells (concurrently for the cell pairs) in a given parameter space. Optimality is understood in the sense of ensuring the maximum possible RCS bandwidth. The cell optimization is implemented as a grid-confined exhaustive search followed by local tuning. This approach is computationally feasible due to low dimensionality of the unit cell parameter space. It guarantees global optimality, and eliminates the need for the employment of stochastic search routines. At the same time, excellent accuracy of the metamodel ensures good agreement with EM simulation data over broad frequency range.
The presented approach allows for fully automated and globally optimum metasurface design within the assumed unit cell topology and the parameter space. It has been used to develop a novel checkerboard metasurface featuring 6 dB RCS reduction in a frequency range from 16 to 37 GHz. VOLUME 9, 2021 The design is validated both numerically and experimentally, and shown to outperform the state-of-the-art benchmark structures with respect to the RCS reduction bandwidth. The technical novelty and the major contributions of this paper can be summarized as follows: (i) the development of a surrogate-assisted framework for reliable and efficient design optimization of checkerboard metasurfaces; (ii) the numerical verification of the framework as well as demonstration of its utility in the context of metasurface design, and (iii) the development of a novel high-performance checkerboard metasurface for broadband RCS reduction. It should be emphasized that the presented framework is-to the authors knowledgethe first comprehensive approach proposed in the literature for globally-optimum design of the unit cell geometries by means of fast metamodels.
The remaining part of the article is organized as follows. In Section II, the motivation for the proposed design framework is discussed, followed by the design and modeling of the unit cell, later used to illustrate the operation of the procedure, and the development of the broadband RCS reduction metasurface. In Section III, the description of the proposed surrogate-based approach and surrogate-assisted global optimization algorithm is provided. In Section V, a novel checkerboard metasurface is implemented and its scattering performance is investigated using full-wave EM simulations and physical measurements of the fabricated prototype. Section VI concludes the paper.

II. PROPOSED DESIGN APPROACH
This section briefly discusses the challenges of EM-driven metasuface design, and provides a motivation for the development of novel techniques that are not only more efficient than the traditional methods in computational terms, but also more reliable. Furthermore, a specific example of a unit cell (metasurface building block) is introduced to be used for the purpose of explaining the proposed machine-learning-based design methodology, and to develop a new high-performance metasurface featuring broadband RCS reduction.

A. MOTIVATION
Metasurface development necessarily involves full-wave EM analysis as the only tool capable of accurate evaluation of scattering properties of geometrically complex structures. Needless to say, the critical stage of the process, i.e., tuning of the unit cell geometry parameters to obtain desired phase characteristics has to be carried out at the level of EM simulation models. The fundamental challenges associated with parameter adjustment include: • High simulation cost of the building blocks and the entire metasurface; • Potential multi-modality of the optimization task resulting from the necessity of considering broadband responses, as well as mutual relationship between the unit cells of different geometries (zero/one cells); • The lack of reasonable initial designs.
The last two factors generally lead to a situation where yielding satisfactory design requires the employment of global search routines, which are extremely expensive when executed directly the level of EM simulation models. Clearly, optimum design of metasurfaces requires the development of novel procedures, capable of addressing the aforementioned difficulties. This work proposes utilization of data-driven modeling techniques to expedite the design process and to improve the optimization reliability. Towards this end, we utilize fast surrogate models (here, kriging interpolation [38]), as well as a combination of global and local optimization algorithms. The details of the framework will be presented in Section III, whereas its performance will be demonstrated in Section IV through the design of a chessboard metasurface featuring broadband RCS reduction. B. UNIT CELL GEOMETRY Figure 1(a) illustrates the geometry of the unit cell design utilized in this work. As shown, the topology resembles the crusader cross. The function f (t) parameterizing the cross arm has the following analytical form where p, b, and d, are the adjustable parameters of the cell that determine its overall shape and size. The specific data concerning the parameter space (lower/upper bounds) will be provided in Section III. This particular geometry has been chosen in order to ensure sufficient flexibility of the unit cell (cf. Fig. 1(b)) while limiting the number of adjustable parameters (here, three).
The latter facilitates the metamodeling-based optimization process, especially the construction of fast surrogate model.
A ground-backed Arlon AD250 lossy substrate (ε r = 2.5, h = 1.5 mm, tanδ = 0.0018) is used in the unit cell design. During the simulations, metallization is represented as perfect electrical conductor (PEC). The overall size of the unit cell is W s × L s = 6 × 6 mm 2 .
It should be noted that the geometries in Fig. 1(b) are for illustration purposes only, and they do not correspond to the final design. Notwithstanding, they are selected to illustrate the unit cell topologies in the assumed parameter space, and, thereby, to demonstrate the topological flexibility of the cell design.
It should be emphasized that the conventional design approaches are not reliable when optimizing such a topology where a small change in the design parameters drastically changes the cell geometry, and, consequently, the reflection phase. This applies to both interactive methods relying on parameter sweeping, but also direct EM-driven optimization techniques, the application of which is hindered by the entailed computational expenses.

III. OPTIMUM UNIT CELL DESIGN BY SURROGATE MODELING
In this section provides a description of the proposed datadriven approach to design optimization of the unit cell. We start by outlining the complete methodology, followed by a detailed explanation of the important components of the procedure. Utility of the proposed framework in the design process of unit cells is also considered. Demonstration of the novel metasurface based on the optimized cell geometries will be provided in Section IV.
The optimization procedure proposed in this paper accounts for geometrical flexibility of the unit cells, which makes global search necessary. At the same time, it capitalizes on the fact that the considered parameter spaces are of low dimensionality, which allows for a construction of fast metamodels, and realization of the global search process in a deterministic manner. As a result, it guarantees identification of a globally optimum design within a reasonable timeframe and it is fully deterministic. The latter alleviates the difficulties pertinent to poor repeatability of solutions, featured by nature-inspired algorithms (the latter currently being the methods of choice for solving this type of problems). At the same time, utilization of surrogates speeds up the search process when compared to direct EM-driven optimization using, e.g., population-based methods.

A. DESIGN METHODOLOGY
The goal of the proposed metamodeling-based design approach is to find a pair of unit cell geometries featuring the phase difference within the range of 180 • ± α max over a possibly broad frequency range F. Here, α max is set to 37 • , which is the value recommended in the literature (e.g., [24]). The operation of the optimization framework is outlined below, whereas the details concerning its major components are provided in Sections III.B through III.E. The vector of adjustable variables of the unit cell, and the response of its EM simulation model will be denoted as x = [x 1 . . . x n ] T ∈ X , and P(x), respectively. The latter represents the phase reflection characteristics. The parameter space X is determined by the user-defined lower and upper bounds l = [l 1 . . . l n ] T and The unit cell optimization is carried out over the Cartesian product X × X and aims at finding the vector that represents a pair of cell geometries corresponding to the maximum (continuous) range of frequencies for which the condition mentioned at the beginning of the section, i.e., 180 In plain words, we strive to determine the dimensions of both unit cells so that the aforementioned phase condition is satisfied for as broad frequency range as possible. The cells have to be optimized concurrently, because the phase difference simultaneously depends on both parameter vectors. Consequently, all dimensions are aggregated into a single vector x p . Formally, the design problem can be stated as follows: The analytical form of the objective function U has been given in Section III.C. The algorithmic flow of the optimization process is as follows: 1. Uniformly allocate N samples x (k) , k = 1, . . . , N , within X and acquire the responses P(x (k) ) from the EM simulation model; 2. Construct a Kriging surrogate S(x) in X using {x (k) , P(x (k) )} k=1,...,N , as the training dataset (cf. Section III.B); 3. Find the initial approximation x (0) p of the global optimum of the surrogate S (in an exhaustive manner) on the structured grid (cf. Section III.D); 4. Find the refined design x * p by solving (1) using x (0) p as a starting point. The refinement process is realized using local search routines (cf. Section III.E). In Step 1, the algorithm starts by uniformly allocating samples within the parameter space and acquiring the training data through EM simulation of the unit cells. The purpose of the training data acquisition is to gather information about the properties of the unit cells in terms of their phase characteristics across the parameter space. This knowledge will be then encoded for further use in the form of a fast surrogate model, which will replace expensive EM simulation in the design optimization process.
In Step 2, a kriging metamodel is constructed to be used as a predictor of the cell phase characteristics over the space X . The metamodel makes predictions about the unit cell phase characteristics as functions of the geometry parameters of the cell. Because it is essentially an analytical model (kriging surrogates are combinations of low-order polynomialbased regression models and linear combinations of kernel functions, e.g. Gaussian), it is fast to evaluate. Furthermore, it is interpolative, i.e., it agrees perfectly with the EM simulation data at the training locations.
In Step 3 of the procedure, the metamodel is employed in the global search. This step, described in detail in Sections III.C and III.D, employs the objective function (3) and carried out exhaustive search over a dense rectangular grid defined over the parameter space. This way of implementing the search process is justified by low dimensionality of the problem, the availability of fast metamodel. It has significant advantages over, e.g., nature-inspired populationbased procedures for the considered case because it is fully deterministic and guarantees identification of the optimum design when coupled with the local refinement.
In Step 4, the resolution of the design found through grid-constrained search is refined through conventional local (gradient-based) optimization. The details are provided in Section III.E. At this stage, the objective function (3) is used as well. As mentioned before, the utilization of the surrogate allows for expediting the optimization procedure to a great extent as compared to direct EM-driven optimization. The flow diagram of the proposed surrogate-based design framework has been shown in Fig. 2.
An alternative approach in this venture could be the utilization of physics-based surrogate models [47], which have become popular in high-frequency design over the last years. Physics-based methods exploit the problem-specific knowledge, typically, in the form of low-fidelity EM or equivalent network models. Some of popular techniques of this class include space mapping [48], and response correction methods (e.g., shape preserving response correction [49], adaptive response scaling [50]). However, in the considered case of unit cell optimization, the employment of data-driven surrogates seems more appropriate having in mind low dimensionality of the parameter space as well as the fact that global exploration is needed. These, along with the lack of convenient candidates for fast low-fidelity representation makes physics-based surrogates impractical.

B. SURROGATE MODELING
The surrogate model S is constructed within x ∈ X using kriging interpolation [38]. The surrogate is identified using the training data samples {x (k) , P(x (k) )} k=1,...,N , where P(x) represents the response of the EM-simulation model, whereas N denotes the total number of samples. The design of experiments strategy is a rectangular grid 7×12×7 (thus, N = 588), which is a suitable arrangement due to low-dimensionality of the parameter space.
The number of grid nodes in each direction is determined based on the large-scale sensitivity analysis with a larger number of nodes set up for the second variable, which has been found to affect the unit cell phase characteristics in a more significant manner than the remaining variables. The kriging model is set up with the first-order polynomial regression model used as a trend function, and a Gaussian correlation function.

C. OBJECTIVE FUNCTION DEFINITION
The design task has been formulated in Section III.A (cf. (1)) as identification of a pair of unit cell geometries x * p = [(x (1) * ) T (x (3) * ) T ] T that maximize the frequency range for which the phase difference satisfies the condition 180 •α max ≤ P(x (1) * , x (3) * ) ≤ 180 • + α max . The analytical form of the objective function U is defined as where f L and f R are the frequencies determining the largest frequency interval for which the phase difference condition is satisfied for all frequencies f ∈ [f L , f R ]. The minus sign in (3) allows for turning the maximization task into the minimization one according to (1). It should be noted that both frequencies are extracted from the phase characteristics of the unit cells using a postprocessing routine implemented in Matlab.

D. GLOBAL OPTIMIZATION
Step 3 of the optimization procedure (cf. Section III. p of the global optimum of S is found as (1) ) T (x (3) ) T ])) (4) 46748 VOLUME 9, 2021 In other words, x (0) p is obtained by searching through all possible pairs of unit cell geometries x (1) ∈ M m1...mn and x (2) ∈ M m1...mn and determining the one that minimizes U . Note that this is an exhaustive search but its computational cost is negligible because the surrogate model S is fast, and the number of parameters is low. Additionally, the entire process is vectorized to further speed-up the operation. In this work, we use m k = 9 for k = 1, . . . , n.
The optimization procedure is governed by the following control parameters: • The number N of the training data points to construct the surrogate model. This number is adjusted to ensure that the surrogate model accuracy in terms of the relative RMS error is at the level of one percent (which gives almost perfect visual agreement between the EM simulated data and the metamodel outputs; • Density of the search grid m k , k = 1, . . . , n. This parameter is of secondary importance because the objective of global search is only to provide a starting point for design refinement (cf. Section III.E), i.e., to ensure that the grid-constrained optimum is sufficiently close to the global optimum. The value used in this work (m k = 9) by far exceeds this requirement. It can be observed that the optimization procedure has only two control parameters, both of which can be easily adjusted to ensure the reliability of the process. This is one of important advantages of the method. The local design refinement uses off-the-shelf algorithm (cf. Section III.E) with default setup, i.e., no control parameters have to be adjusted.

E. DESIGN REFINEMENT
The last stage of the optimization process (Step 4) is local refinement, using x (0) p found in Step 3, as a starting point. The refinement is executed using Matlab's fmincon procedure [46], which is a variation of the sequential quadratic programming (SQP) method [46]. Again, the computational cost of this stage is negligible because it is executed at the level of the kriging metamodel. Figure 1(b) provides a general idea about the type of structures under consideration. For the unit cell of Fig. 1(b), there are three parameters, p, b, and d, that determine the shape of the unit cell. Hence, the vector of designable variables is  deviation of 1.7 degrees. This means that the metamodel is very reliable, especially when considering the typical range of the unit cell phase response (>400 degrees). Figure 3 shows the surrogate and EM-simulated cell responses at the selected test locations. The visual agreement between the two data sets is excellent, which corroborates the design utility of the metamodel.  Figure 4 shows the cell geometries, for convenience, labeled as Cell 0 and Cell 1. Verification of these designs has been conducted by comparing their phase characteristics with the EM simulation data. As demonstrated in Fig. 5, the agreement between the surrogate and EM-simulated responses is excellent. This confirms the efficacy of the proposed machine-learning-based design framework.

F. OPTIMIZATION RESULTS
The reflection phase and amplitude of the unit cells along with the reflection phase difference between the two cells are shown in Fig. 6. It can be observed that the condition 180 • -37 • ≤ P(x (1) * , x (3) * ) ≤ 180 • + 37 • is satisfied for the frequencies from 16 GHz to 35 GHz. Hence, more than 19 GHz RCS reduction bandwidth can be anticipated [25]. It should be reiterated that the objective of the optimization procedure is to find a globally-optimum design of the unit cells that maximizes the RCS reduction bandwidth, i.e., a pair of designs featuring the phase difference of 180 ± 37 degrees over possibly a broad frequency range. The outcome of the optimization procedure (pair of unit cells) serve as a building VOLUME 9, 2021  block of a high-performance RCS reduction metasurface as described in Section IV.

IV. NOVEL METASURFACE CONFIGURATION. NUMERICAL AND EXPERIMENTAL VALIDATION
This section introduces the configuration of a novel metasurface design. The monostatic and bistatic RCS performance of the proposed structure is discussed in detail. The experimental setup is also illustrated, along with the comparison of simulation and measurement results of the checkerboard measurface. Finally, benchmarking against the state-of-theart designs is discussed.

A. CHECKERBOARD METASURFACE PERFORMANCE
The operating principle of a checkerboard metasurfaces is to interleave the two structures featuring 180 • phase difference so that the backscattered fields are cancelled out, and a distinct scattering patterns are produced. Theoretically, monostatic and bistatic RCS reduction can be approximated by the array theory [51]. The concept of the RCS reduction can be understood by recalling a planar array having a progressive phase shift of 180 • among elements within a particular frequency band. In other words, the checkerboard measurface exploits the anti-phase reflection property of periodic arrays to manipulate the scattering behavior. In order to enable the aforementioned property, in the first step, the optimum unit cell designs, i.e., Cell 0 and Cell 1, featuring a phase difference of 180 • ± 37 • are obtained (cf. Section III). Hereafter, the periodic arrays containing multiple copies of Cell 0 and Cell 1 as the building blocks are employed in an alternate manner to realize a checkerboard metasurface. Figure 7 illustrates the proposed checkerboard metasurface comprising thirty-six elements: eighteen 4 × 4 periodic arrays of Cell 0 and eighteen 4 × 4 periodic arrays of Cell 1. Subsequently, the resulting 6 × 6 checkerboard surface is characterized. Note that the size of the periodic arrays is decided by considering the fact that diffractions due to discontinuities among the neighboring arrays do not significantly contribute when the overall size of a single array is greater than half wavelength [52]. The total size of the surface is W s × L s = 144 × 144 mm 2 . The inter-element spacing of individual unit cells in an array is s = 6 mm.
The surface is implemented on a ground-backed Arlon AD250 lossy substrate (ε r = 2.5, h = 1.5 mm, tanδ = 0.0018). To test the RCS performance of a proposed metasurface, a PEC surface of a similar size is also implemented to be utilized as a reference surface. The time-domain solver of CST Microwave Studio is used for both the monostatic and bistatic RCS analysis.
In order to validate the anticipated broadband RCS reduction of the proposed metasurface, its monostatic RCS performance for normal incidence has been determined. Figure 8 shows the reflection characteristics of the PEC surface along 46750 VOLUME 9, 2021  with the proposed metasurface. It is apparent that the RCS reduction occurs in a broad frequency range, i.e., from 15.7 GHz to 38 GHz, which confirms the low observable property of the metasurface.
The 3-D bistatic RCS patterns of the proposed metasurface and the metallic surface of same size has been presented in Fig. 9. It can be observed that the reflected waves from the proposed surface, under normal incidence, scatter into four diagonal planes. It corroborates minimum reflections from the metasurface in the boresight direction, as the incident waves are reflected into different directions. On the contrary, the metallic surface features strong reflections in the boresight direction, in a single lobe, when the plane wave impinges on it from the normal direction.  The scattered field versus the elevation angle theta θ along the principal and the diagonal planes are demonstrated in Figs. 10 and 11, respectively. The bistatic RCS performance of the proposed metasurface is compared with the PEC surface. The results indicate that the maximum RCS in the principal planes is 16.0 dB lower than the maximum RCS for the PEC ground plane, at both considered frequencies. Subsequently, in the diagonal planes, the maximum RCS of the proposed surface is 15.2 dB lower than a PEC ground plane. Hence, a significant RCS reduction has been observed for the proposed metasurface in the principal as well as the diagonal planes. This reduction occurs because the reflected fields are redirected into four main lobes, instead of the single main lobe of the PEC surface, (cf. Fig. 9).

B. MEASUREMENT SETUP AND EXPERIMENTAL VALIDATION
Following the EM-simulation-based verification, the prototype metasurface has been fabricated and measured. Figure 12 show a photograph of the structure. The RCS has been measured in terms of reflectivity, owing to limited amenities. The same size PEC surface has been used as a reference to determine the RCS reduction of our metasurface. For the sake of measurements, two PE9850/2F-15 horn antennas, operating from 26.5 GHz to 40.0 GHz, have been utilized as a transmitter and a receiver. The monostatic RCS characteristics of a checkerboard and a PEC surface has been evaluated by measuring the antenna transmission coefficients. The block diagram of the measurement setup has been provided in Fig. 13. The measurements have been carried out using the anechoic chamber of Reykjavik University (cf. Fig. 14). The comparison between the simulated and measured RCS reduction is depicted in Fig. 15. The agreement between the datasets is very good. A certain discrepancy can be attributed to the fabrication tolerances, as well as the misalignment of the transmitter/receiver antenna with respect to metasurface during measurements. The latter is essential considering that the experimental setup is for capturing reflections. A slight misalignment could  lead to relatively high inaccuracies. Nevertheless, the measurements corroborate 6-dB RCS reduction within the frequency range of 26.5 GHz and 38 GHz. As mentioned before, the lower edge is limited by the available hardware. The measured RCS reduction bandwidth of the proposed checkerboard metasurface and the expected bandwidth anticipated from the phase difference curves (cf. Fig. 6) are similar.
The above findings allow us to conclude that the proposed checkerboard metasurface features low scattering property in a broadband frequency range, and, therefore, it has the potential to replace the metallic surfaces in the applications where high stealthiness is essential.

C. BENCHMARKING
For the sake of benchmarking, the performance of the proposed checkerboard metasurface has been compared with the recent metasurfaces from the literature, see Table 1. The comparison is carried out in terms of the RCS reduction bandwidth. It can be observed that the proposed metasurface outperforms other designs with respect to fractional/relative RCS reduction bandwidth. It should be emphasized that apart from proposing a novel metasurface, an efficient surrogateassisted design framework is also provided-for the first time-to facilitate the design procedure of such surfaces. As a matter of fact, it is rigorous optimization that provides a competitive edge over less formal design approaches, and manifests itself through better properties of the resulting metasurface. As mentioned before, the crucial components of the procedure are those that take into account the specifics of the problem: the parameter space dimensionality, expensive (EM-based) evaluation of the unit cell characteristics, and the need for global search.

V. CONCLUSION
This article proposed a surrogate-assisted framework for rapid design of high-performance metasurfaces featuring broadband RCS reduction. Low RCS of a surface translates to its low observable nature, which is highly desirable for the stealth technology. Our procedure involves a construction of a fast metamodel that replaces the CPU-intensive EM simulations in both stages of the design process, i.e., the global search, and local (gradient-based) refinement. The optimization is executed to identify the optimum until cell geometries within the user-defined bounds. By employing the proposed methodology, a computational burden of the design process can be significantly reduced. Finally, a novel checkerboard metasurface, enabling broadband RCS reduction, has been developed using our framework. The monostatic and bistatic performance of the proposed checkerboard metasurface has been validated both numerically and experimentally. The numerical results indicate that the metasurface features low scattering property in a broadband frequency range, i.e., from 15.7-38 GHz. The experimental data confirms these findings starting from 26.5 GHz, which is due to the limitations of the available hardware. The proposed metasurface has been benchmarked against state-of-the-art designs demonstrated to be superior in terms of the RCS reduction bandwidth. This also validates the design utility of the presented metamodeling-based procedure in the context of metasurface development. As a matter of fact, the design of the above structure provides a link between the theory (here, a simulation-based design optimization procedure), and application, which is the development of high-performance metasurface with the intended use in stealth technology.
The authors believe that this study is a step toward exploring the data-driven techniques in the design of high-performance metasurfaces for RCS reduction, where intuition-inspired methods are still widespread although generally lack the ability to yield truly optimum results. Application of surrogate-based methods, including fast metamodels, improves reliability, enables global optimization in reasonable timeframe, eventually leading to the improvement of metasurface performance figures, as demonstrated through the specific design proposed in this work.