The Design and Analysis of Electronically Reconfigurable Liquid Crystal-Based Reflectarray Metasurface for 6G Beamforming, Beamsteering, and Beamsplitting

Reconﬁgurable intelligent surfaces (RISs) are becoming increasingly popular in the ﬁeld of wireless communications, given their potential to ameliorate the challenges faced in millimeter-wave wireless communications. Speciﬁcally, liquid crystal (LC)-based RISs have demonstrated numerous advantages over other types of RISs in terms of cost, complexity, and radiation efﬁciency in the high-frequency regime. This paper presents the design, algorithms, implementation routines, and simulation of a novel reconﬁgurable LC-based reﬂectarray metasurface operating at 108 GHz. The scanning range of the proposed device is ± 40 ◦ (azimuthally and horizontally), with an average scanning beamwidth of 8 . 6 ◦ . We present semi-analytical ﬁndings on the scalability and phase continuity of our design, showing how key performance indicators (KPIs) are affected by dimension and phase degree-of-freedom changes. We demonstrate agreement between our semi-analytical models and full-wave analysis, focusing on genetic algorithm (GA)-optimized beam manipulations. Our results present a feasible workﬂow that enables dynamic beamforming, beamsteer-ing, and multibeams at 108 GHz, and is easily scalable for applications in other 6G and beyond frequency spectra.


I. INTRODUCTION
In the world of wireless communications, the frequency spectrum of interest is constantly shifting higher as a result of spectrum occupation; just as cellular companies began rolling out 5G in 2019, beyond-5G investments and explorations have already begun.The Office of Communications (Ofcom) announced a proposal to increase access for the above 100 GHz spectrum for the development of new and innovative services [1]; the Federal Communications Commission (FCC) announced the opening of the 95 GHz to 3 THz spectrum for experimental 6G purposes [2]; the International Telecommunication Union (ITU) is already planning the spectrum management aspects of The associate editor coordinating the review of this manuscript and approving it for publication was Cunhua Pan .
6G wireless communications [3], for emerging applications such as holographic-type communications and multi-sensing networks.Not only are we progressing towards the higher generations, but also at a faster pace; as 5G infrastructure is being laid out, there are already considerable interests and actions taken in preparing for 6G, especially in fields such as high-fidelity holographic communications and connectivity for all things [4].
With an increasing network traffic demand that is expected to approach Tbps rates [5], there is no option but to begin utilizing the terahertz (THz) band (0.1-10 THz), which offers enormous potential in terms of the bandwidth and highdata-rate transmission [5].Despite the pouring interest in the THz band, two notable issues still pose significant challenges to its practical commercialization: non-line-of-sight (NLOS) and free-space path loss.With millimeter and submillimeter wavelengths, many previously ''transparent'' objects, such as tree leaves and human body parts, will become barriers to EM propagation.In addition, surfaces that are ''smooth'' in the lower frequency regimes become rough and dispersive, adding greater complexities to traditional reflective channel models.Furthermore, the free-space path loss is proportional to the frequency squared, implying that a much greater power source and beamforming capabilities will be needed if modern transceivers are to work at the same data rate over the same distance.The significant increase in power consumption and losses associated with the feeding network make current architectures impractical [6], [7].
RISs are being increasingly studied as potential mediators in addressing the problems of high-frequency wireless communications owing to their relatively low manufacturing cost, structural simplicity, and practicality.RISs usually function as programmable reflectors or transmitters to further strengthen and guide the EM wave toward the desired direction.As an intermediate bridge, RISs can significantly lower the extreme requirements of 6G transmitters while avoiding the high costs associated with employing greater densities of micro-cells.Recent studies [8]- [12] have already stepped into the signal propagation and channel estimation domains and have shown that the inclusion of RISs can assist the data transmission to cell-edge users and create more power-efficient propagation architectures.
In this work, we present the study of an LC-based reconfigurable reflectarray metasurface.Depending on the practicality and application scenarios, the proposed device can function both as a ''smart relay'', reflecting impinging signals from nearby base stations or as a standalone access point with an integrated feed source.
Expanding from our previous studies with binary phase LC-based reflectarrays [13], [14], in this work, we updated the unit cell structure and geometries to include a biasing circuit and LC spacers for realistic full-wave results.These updates are due to our plans to prototyping the device in future works.Therefore, we would like to include as many details as possible to obtain more accurate predictions of what would occur in reality.These unit cell updates resulted in negligible changes in the optical responses of the device (<1% deviations in the unit cell reflection phase and amplitude at the operating frequency).We improved the accuracy of the theoretical far-field approximation with a semi-analytical approach by incorporating a full-wave simulation-derived unit cell radiation pattern, which should include previously ignored effects of material losses, coupling, surface waves, and edge effects.The GA optimization algorithm has been matured and updated with the inclusion of probabilistic adaptations to achieve higher efficiency and quality pattern synthesis, improving the side-lobe level (SLL) and resolutions of multibeams that were previously either divergent or deficient.In addition, we present new semi-analytical and full-wave results from the scalability and phase continuity studies, where we analyze the relation between KPIs (directivity, SLL and beamwidth) and structural parameters (such as unit cell and supercell sizes), as well as between KPIs and phase continuity (the degree of control over reflected phases), demonstrating the possibilities of up/down-scaling the device to tailored application demand.Lastly, we unified the far-field plots in 2D surface plots and presented a full-wave simulation-derived table of KPIs for different scan angles.
The rest of the paper is organized as follows: in Section I, we introduce the unit cell structure, theoretical constructs on the far-field calculations, and the GA used in optimization.In Section II, the unit cell design principles, formulation for the far-field, implementation of the GA, and the full device are described in detail.In Section III, we present the scalability and continuity results from full-wave simulations, GA-optimized configurations and corresponding beam profiles, and the KPI table for different scanning angles.In Section IV, we summarize our findings, comment on future work, and suggest areas for improvements.

II. RELATED WORKS
Lumped element-based reconfigurable reflectarrays have been widely studied [15]- [17]; these devices make use of control elements such as PIN diodes and varactor diodes to form different circuits under electrical inputs, resulting in different electrical properties (such as electrical length variations) of the unit cells.Although practical, lumped elements often cannot perform as well in the higher frequency regimes (>60 GHz, for instance), owing to the intrinsic parasitic losses associated with diodes.Mechanical devices [18], [19] utilize micro-electrical-mechanical systems (MEMEs) to physically vary the dimensions and hence the electrical properties of the devices.Mechanical micromotors can be tailored to achieve high performance and accuracy, but these devices often suffer from high cost (owing to the manufacturing of MEMEs) and relatively low reliability when compared to electrically controlled metasurfaces.
LC-based reconfigurable reflectarray metasurfaces have recently gained much popularity [20]- [22], [30], as they have been demonstrated to be relatively efficient, cheap, and easy to design or manufacture in high-frequency regimes (>100 GHz), especially when compared to phase shifters at those frequencies.The fabrication of LC is straightforward, as the LC unit cell is often as simple as a resonator cavity enclosing the LC material and does not involve much soldering of micro-electrical components.The efficiencies of reflectarrays are commonly higher than other types of tunable materials, such as photo-induced Fresnel Zone Plate [23], [24], or transmit-type lumped element reconfigurable metasurfaces [25].The reflection coefficient of the device studied in this study is approximately 0.88.
The principle behind tuning LC metasurfaces lies in the nematic nature of LC molecules; these molecules behave like electrical dipoles, which will align with an applied electric field.By tuning the electric field, the molecular alignment can be controlled, which in turn determines the permittivity of the LC.Because LC is used as a substrate for the resonating VOLUME 9, 2021 patch antenna, the change in permittivity will directly affect the resonant frequency and thus the phase of the reflected EM wave.
Recent studies [21], [22] have established very strong agreement between the theoretical and measured relation of the applied voltage and the corresponding S-parameter changes (most notably, between the voltage and the phase of the reflected wave).This was achieved through the accurate modeling of the LC substrate, which ephasizes the anisotropic and inhomogeneous nature of LCs.With a strong theoretical model, researchers have been able to design multi-resonant antennas for LC reflectarrays operating in the F-band [26].These kinds of antenna ''unit cells'' consists of the LC substrate, which is sandwiched between the resonant antenna and the ground plane, silicate superstrate, spacers to enclose the LC, and the biasing circuit.
Despite the growing developments, most LC reflectarray metasurfaces [21], [22], [26], [27] are designed with a lineby-line addressing technique, which often limits the scanning of the device to two dimensions.This is due to the fact that in traditional liquid crystal displays (LCD) addressing techniques, the control of the individual elements is achieved on the basis of a threshold voltage, which is not valid in our scenario, as any non-zero voltage level is associated with a unique phase difference in the reflected EM wave.In our work, we propose a 20 × 20 individually biased structure, so that each unit cell antenna can be actively controlled to be in either the ON or OFF state, which gives the reflected EM wave the freedom of scanning both in azimuth and elevation, and the capability to generate multibeams.

III. METHODOLOGY
This section introduces the structural designs of the unit cell and the theoretical model for the far-field electric field pattern, on which we base our beam-manipulation optimization algorithm.We used the GA for the optimization problem, which is further explained in this section.
The general approach can be outlined as follows: We started with the design of a patch antenna given the operational frequency.Then, using full-wave simulations, we performed a parameter sweep on some of the antenna's physical parameters to obtain the desired optical responses.Subsequently, we confirmed the accuracy of the unit cell by performing full-wave full-device simulations with known ON/OFF configurations.Finally, the GA was used to solve the multibeam configurations.These processes are described in detail in the following sections.

A. UNIT CELL
In this study, we used GT3-23001 as the model for the LC substrate.The choice of LC is based on its permittivity frequency response at >100 GHz, which for GT3-23001 has been reported to be stable and suitable for achieving the required reflection phase change [22], [28]- [31].GT3-23001 also exhibited stable responses at temperatures ranging from −20 •  ∼ 100 • , making it suitable for outdoor applications.In Fig. 1, we present the unit cell antenna element's schematics.The dimensions of the unit cell are optimized using cross-platform routines to achieve the desired reflection phases and amplitudes for the ON and OFF states.The selection of the 108 GHz operational frequency is not of particularity, but to demonstrate a >100 GHz application scenario in beyond-5G systems.
The unit antenna element is essentially a voltage-biased patch antenna, forming a capacitative relation with the ground plane.The simulations of the ON and OFF states were achieved by adjusting the corresponding LC substrate permittivities, which directly affected the resonant frequency of the patch antenna.When using the same frequency source, the shift in the resonant frequency translates to drive or lag in the reflection phase.
The updated unit cell now includes the biasing circuits and spacer walls that enclose the LC substrate and the ground plane from neighboring unit cells.The biasing circuit is what will be used in the prototype to create the capacitative effect between the top patch antenna and the ground plane, which is how the LC's relative permittivity will be controlled in reality.The top (patch antenna) and the bottom (ground plane) plates will function as the two capacitative plates, and the two biasing lines will feed each of the plates separately with the required voltage difference.In simulations, we have simply adopted the outcomes of applying the voltage difference from the biasing circuit: a change in LC's relative permittivity, as shown in Table-1.The design of the biasing circuit also considers the surface current that forms when the patches are excited by the source, and the placement of the biasing line also aims to avoid the interference of DC and RF signals.
The periodicity of the unit cell antenna was 1 mm, and the full device had a 20 × 20 structure, which was 2 cm × 2 cm in dimensions.Depending on realistic application scenarios, the device can increase its aperture for more practical applications, either by 1) incorporating more antenna elements or 2) utilizing multiple devices as a collective and collaborative network.

B. THEORETICAL FAR-FIELD
The antenna far-field radiation pattern [32], E(θ, φ), is essentially a summation of EM radiation from the feed source to individual antenna elements and finally to the far-field region (where there are only dependencies on θ and φ).The feed source radiation pattern is modeled as cosine to the power of q.The exponential term is the spherical wave radiation from each individual antenna element, which can have unique amplitude ( mn ) and phase terms (e iφ mn ).In (1), we have the far-field radiation pattern, E(θ, φ): where φ mn is the configuration matrix of the ON/OFF state, which is multiplied by the phase difference between the ON and OFF states, φ : This antenna radiation pattern approximation is used for basic far-field calculations; however, for multibeam optimization with the GA algorithm, we adopted a more accurate semi-analytical model, which could include more intricate EM details such as coupling, edge and surface effects: mn • e iφ mn .
(3) Instead of approximating the unit cell antenna radiation pattern as cosine to the power of q, we performed a full-wave simulation of the unit cell radiation pattern under periodic conditions (with infinitely many unit cells as the repeating periodic boundary condition), as shown in Fig. 3.We then replaced the cosine unit cell radiation pattern approximation with the full-wave derived far-field radiation pattern, A(θ, φ).In this way, we have a more reliable far-field prediction, which is especially needed when synthesizing a higher number of multibeams, where noises/SLL's are higher.

C. DIRECTIVITY
The directivity of the metasurface is calculated as follows: where U (θ, φ) = |E(θ, φ)| 2 is the radiation intensity, which is divided by the total radiation power.The integration in θ is limited to π/2 because of the nature of reflectarray operations, where only the upper semi-sphere where EM waves are reflected is of interest.

D. PHASE DISTRIBUTION
The configuration of the ON/OFF states of the unit cells results in a specific distribution of reflection phases (specifically, a 180 • reflection phase difference between the ON and OFF states).The reflection phase distribution across the device's surface ultimately determines the behavior of the reflected waves according to the Huygens principle.A fundamental theoretical concept in synthesizing simple phase distributions is the Generalized Snell's Law (GSL).GSL is derived from Fermat's principle, whereupon momentum conservation, we arrive at the continuity of the phase condition at the reflection boundary, and when an extra phase profile is introduced, this momentum conservation leads to generalized Snell's Law for reflection: where n i is the refractive index, and the sub-indices represent either incident or reflection.From the generalized Snell's law, we can see how an arbitrary reflection phase (also known as anomalous reflection, when θ i = θ r ) can be achieved when a phase gradient d dx can be artificially engineered.It is relatively straightforward to see how beamsteering, namely anomalous reflection, can be achieved through the introduction of a constant phase gradient; in 2D array theory, it is well known that the continuous phase change required for beamsteering applications is as follows: where the indices b represent the beam-steered, i for incidence, and x i and y i for the coordinates of each of the unit cell antennas.For more complex beam profiles, such as multi-directional multibeams, the solution can be retrieved through numerical optimization methods, which in our case is the GA.

E. GENETIC ALGORITHM
In order to solve the inverse problem described in (1), where we have an intended far-field radiation profile E(θ, φ) and are looking for the input 20 × 20 configuration matrix φ mn .This problem can be addressed using GA, which is an optimization algorithm motivated by biological evolutionary processes.As shown in Fig. 4, we start the GA optimization algorithm with an initial population (where the population is the matrix of ON/OFF configurations); this can either be randomly generated or, in our case, loaded from a library of pre-stored known configurations, to fasten the optimization.The cost function, otherwise known as the fitness function, determines the quality of our optimized configuration, which is simply the difference between the intended far-field radiation profile and the optimized profile: where E target is the desired electric field and E the current outcome.

IV. IMPLEMENTATION ROUTINES A. UNIT CELL
With the fundamental patch antenna design formulas, known operational frequency of 108 GHz, and periodicity of approximately half a wavelength, we arrive at the preliminary dimensions of the unit cell structure.The half-wavelength element spacing was chosen for an optimal compromise between the beamwidth and steerability (and avoidance of grating lobes).With the elementary dimensions established, we can also check the reflected radiation pattern with known configurations of the ON and OFF states.The next step involves unit cell optimization in frequency domain solvers, such as Studio Suite, or the more traditional method of moments (MoM).With these solvers, we can enforce periodic boundary conditions and include coupling effects, which are often significant in subwavelength-spacing resonators.The S-parameters in reconfigurable reflectarrays are defined differently when compared to traditional antenna arrays, as reflectarrays are not fed in the traditional way, but rather through a feed horn that is usually situated on top of the reflecting surface; thus, S11 here refers to the proportion of EM reflection/radiation reflecting away from the plane of the patch antenna and towards the feed source in free space.
With the desired solver and appropriate boundary conditions, we perform parameter sweep over structural parameters (height h and width W ) and store the S11 amplitude and phase results for analysis.We can then post-process the data and retrieve optimal structural parameters, given the desired conditions, which are as follows: a reflection phase difference between the ON and OFF states of 180 • ; second, a maximal reflection amplitude of both states; and third, a minimal reflection amplitude difference between the two states.

B. THEORETICAL FAR-FIELD AND GA
As described in the previous section, the theoretical far-field is calculated with a given configuration matrix that represents the ON and OFF states of the individual unit cell antennas.The ON and OFF configurations govern the radiation pattern of the reconfigurable reflectarray.Thus, we seek to solve the inverse problem in (1), where we have a desired and known radiation pattern (a multibeam steering, for instance) but are looking for a non-trivial input matrix.
Together with the GA toolbox offered by MATLAB, we also developed our own custom probabilistic crossover and mutation functions to better tailor to the optimization problem that we have.The general GA routine is shown in Fig. 4.
As described in the previous section, the cost function determines the optimization subject.For simplicity, we defined it as the difference between the desired outcome and the current outcome, as shown in (7).In the case of multibeam optimization, the fundamental cost function can be expressed as (8), where the E pk s are the locations of the electric field peaks, and Es are the intended locations of the electric field peaks.So if the peaks are far away, then the square of the location difference will be large, whereas if the peaks are located as desired, then the square of the difference will be zero.With this general form of the cost function, we can also perform rudimentary SLL reductions.
It is important to note that the cost function (also known as the fitness function) has the strongest and most direct impact on the effectiveness of the GA optimization (or any optimization), and we always prioritized the adjustment of the cost function over that of the GA operators.This method of optimization is computationally costly, as can be seen from Fig. 5, where convergence is reached after more than a hundred generations.Therefore, this approach is more suitable when used offline to generate a codebook, which can then be stored in the proposed device.We are working on a Gerchberg-Saxton iterative method that can potentially be efficient enough to be deployed for online optimizations.
For the full cost function, one must add the individual weights for each cost for the GA to work optimally.The weights affect the prioritization in optimization, which is especially important for multi-objective optimizations.A generalized form for the complete cost function is expressed as follows: where C is the cost function with input φ mn , the binary configuration matrix of the metasurface, E dk is the desired electric field pattern for the k th electric field profile (such as a peak), E k is the current electric field given the input matrix, H is the Heaviside step function, indicating whether the desired beam is larger or smaller than the current beam, and W (n) is the weight operator that operates on the magnitude difference between the desired and current electric field strengths, providing both product and exponential terms based on the importance of the operand.The determination of the individual weights is the most time-consuming and challenging process, as often only trial and error or educated-guesses methods are available.Through experience, we established a matured optimization cost function weighting method with an improved convergence rate from adaptive crossover and mutation operations.

C. FULL DEVICE
With a known ON/OFF configuration, the 3D model can be generated in a full-wave simulation software that solves Maxwell's equations in the time domain.For this process, we wrote a routine in Visual Basics to automate the 3D modeling process, which reads the configuration matrix and updates the 3D model accordingly.In our case, if the matrix elements are read as numeric ''1'', then the LC permittivity will be adjusted to that of the ON state; likewise, if the matrix element is ''0'', the LC permittivity for that unit cell will be adjusted to that of the OFF state.

V. RESULTS AND EVALUATIONS A. SCALABILITY ANALYSIS
Ideally, the reflectarray would have optimum directivity if it can achieve continuous control of the reflected phases (0 • − 360 • ) from each antenna element, similar to the operation of a traditional phased array.However, achieving a full 360 • and continuous control often comes at the cost of increased power consumption, cost, and design complexity.
Thus, it is of interest to understand the relation between directivity and phase quantization error, which can result from the continuousness of the reflected phase, or the degree of spacing from the subwavelength elements.
In Fig. 6, in a 20 × 20 elements reflectarray, the continuous control over phase produces a directivity of 31.09dB, while the binary phase control produces 26.17 dB, which is approximately 5 dB difference.This significant gain advantage of the continuous phase control device could be of interest to certain applications.However, practically, it will come at the cost of higher power consumption (from the variable voltage controls), lower energy efficiency, and greater complexity in the design, which will ultimately lead to a more expensive device for manufacturing.
In Fig. 7, given the same aperture of the reflectarray, the directivity of λ/2 unit cells is 26.17 dB, while the directivity for λ/8 unit cells is 28.5 dB, approximately 2 dB difference.This shows that by increasing the number of elements while retaining the overall structure size, having more and smaller radiating elements enhances the intensity of the reflected field.Nonetheless, there are certain effects as a result of increasing element density.First, the bandwidth of the full device tends to decrease.Second, when the radiating elements are electrically very small, surface waves or even localized surface plasmons may result at such a scale, thus reducing the radiation efficiency.
Figure 8 shows the directivity of binary phase reflectarrays with three different aperture sizes.20 × 20 elements has a directivity of 31.09dB, 40 by 40 elements gives 36.91 dB and 80 by 80 elements gives 42.65 dB.The aperture of the device is directly related to the gain, so the only underlying consideration here is the application intention and budget available for the appropriate aperture.TABLE 2. Summary of the KPIs for different dimensions of the metasurface when using binary phase unit cells.
In TABLE-2, we notice that the 3 dB beamwidth is halved every time the device area is doubled, whereas the main-lobe gain increases by approximately 5 dB following each device area doubling.

B. UNIT CELL AND FEED SOURCE
The feed horn has a dipole embedded inside, which radiates omnidirectionally, while the horn cover serves as a guide to focus the EM energy toward the intended surface.We designed the horn to maximize the gain achieved and radiation coverage to the surface.The design and optimization of the feed horn are important; otherwise, the assumptions made in theoretical models would mismatch the full-wave simulation conditions and result in inaccurate far-fields.
In Fig. 9 part (a), we can see that the 3 dB beam-width is 40 • , this is designed according to the horn aperture and distance to the reflectarray surface, in order to maximize coverage and at the same time provide relatively uniform strength of radiation across, an assumption made in the theoretical modeling of the excitation source.
In Fig. 10, we demonstrate that the phase difference between the ON and OFF states of the new unit cell, simulated under periodic boundary conditions, is 189 • at the operational frequency of 108 GHz.The magnitude of the reflection is 0.855 for the ON state and 0.857 for the OFF state.Note that the unit cell characteristics are practically identical to our earlier works [13], [14], as the integration of the biasing line and FR-4 spacers were designed to minimize any interference with high-frequency resonances.Nonetheless, the reflection amplitudes of the ON and OFF states at 108 GHz are slightly different; thus, the difference is included in the theoretical far-field model.
In Fig. 11, we demonstrate the reflection phase and amplitude of the structure for oblique incident angles for both the TE and TM modes.We can see that for the reflection phase, both the TE and TM polarizations remain relatively unchanged from 0 • to 40 • incidence angle, after which the reflection phase becomes significantly different (more than 20 • different from the normal incidence phase), affecting the performance of the structure.Both the TE and TM modes have similar phase values until the angle of incidence reaches approximately 40 • .This implies that the device can possibly be programmed to operate for both polarizations up to an incidence of 40 • , without compromising the pattern synthesis capabilities.
The magnitudes of the reflections for both polarizations remained identical until approximately 55 • incidence.Because we adopted a phase-only pattern synthesis process, the amplitude variations did not significantly affect the accuracy of the peak locations, but the overall gain values.However, after 55 • incidence, there begins to exhibit a significant difference in the amplitude value of the two polarizations, which will also affect the accuracy of the pattern synthesis.
Thus, we can suggest that the device is capable of operating from −40 • to 40 • (azimuth and horizontal) for both TE and TM polarizations, as shown in Fig. 13.

C. RADIATION PATTERNS
In this section, we demonstrate the far-field radiation patterns for both the beam scanning and multibeam scenarios.For these results, we used the feedhorn antenna described in the section above at a distance of 10λ above the device.This specific distance is chosen as we intend to simulate an integrated feed horn, where excitation is near the reflectarray surface.
In Fig. 12, we show the general agreement between the theoretical normalized far-field electric field and the full-wave normalized far-field electric field.The ON/OFF states are shown as green/red unit cells in the figure.These far-field radiation patterns were also checked with those in other studies [15]- [17], where the same configurations of ON/OFF were used.
In Fig. 13, we show a collection of the full-wave beam-scanning simulation results, for φ = 135 • and θ scanning from 0 • to 57 • .At angles greater than 40 • , such as at θ = 57 • , not only is the beamwidth greatly increased (and the beam profile became asymmetrical to the peak), but also the directivity significantly decreased.However, from θ = 0 • to θ = 40 • , the directivities are almost identical, with a minor increase in the beamwidth as the scanning angle increases.
In Fig. 14, we show the theoretical and full-wave far-field radiation pattern when the configuration is set to steer the beam to θ = 15 • , φ = 135 • .By using (6), we can determine the phase distribution required to achieve a certain beamsteering angle, which can then be used to construct the  configurations of the ON and OFF states.Similarly, in Fig. 15, the for θ = 30 • , φ = 315 • is shown.
In Fig. 16, we demonstrate the theoretical and full-wave radiation patterns of a GA-optimized three-beam configuration.In Fig. 16  In Fig. 16 (b), we present the full-wave simulation plot for the same configuration that gives the 3-beam split.In this plot, we see that each of the three steered main lobes has a directivity of at least 14 dBi, which shows promising performance in terms of the beam-forming capabilities of our device, given that the horn antenna has a directivity of 13 dBi.
Similarly, in Fig. 17, we demonstrate the GA-optimized four-beam split.Although the main peaks are correct in terms of position, it is noticeable that the noise and SLL are more significant here, likely because of the limitations of the device's dimensions: the lower the number of unit cells, the lower the resolution of the multibeams.
In Table-3, we listed the KPIs that are derived from semi-analytical calculations and full-wave simulations.Although the beamwidths of the semi-analytical model are consistent with those of the full-wave simulations, the semi-analytical main lobes (and consequently sidelobes) are approximately 4 dBi stronger than those of the full-wave simulations.We believe that the inconsistencies in weaker radiation regions, and the overestimation in mainlobe/sidelobe strengths can be attributed to the absence of considerations in the theoretical model of the: 1) coupling effects between the individual antennas, which become significant in the subwavelength regime, and often manifest as resonant features, as observed in the full-wave plots, in terms of the oscillatory blue lines versus the smooth blue lines in theoretical plots; 2) material losses and surface effects, such as surface waves, which can lead to secondary radiation manifesting in the weaker radiation regions and decreasing the main-lobe intensity.However, it is important to note that the decibel-scale plot disproportionately magnifies any numerical effects and errors.
We would also like to point out that, likely due to the step-size and accuracy settings in the full-wave simulations, the lowest far-field radiation intensities recorded in full-wave simulations were much higher in comparison to the lowest radiation intensities in the theoretical far-field plot (−21 dBi in full-wave versus −46 dBi in theoretical, for beamsteering).To clarify the theoretical far-field contour plots in the areas of interest (mainlobe and sidelobes), we extended the range of lower intensity values for the lower color range (lower intensity values are inside the blue color range) while maintaining the higher intensities in a color range that is consistent with color mapping of full-wave plots.

VI. CONCLUSION
In this study, we performed computational studies on a proposed 20 × 20 element LC-based binary phase reconfigurable reflectarray metasurface, which operates at a central frequency of 108 GHz.The phase difference between the ON and OFF states achieved at the central frequency is 177 • , and the reflection amplitudes are 0.88 for both states.
We presented a scalability analysis, semi-analytical, and full-wave simulation of far-field radiation patterns, where we demonstrated the performance of the device in terms of beam scanning and GA-optimized multibeam capabilities through scanning plots and KPI tables, as well as detailed implementation techniques employed to achieve the cross-platform  simulation and optimization.In terms of pattern synthesis, we showed a generalized adaptive GA implementation and principle ideas for designing the cost function for specific beam profiles.
In future work, we would like to include amplitude synthesis, in addition to the phase synthesis.This will allow us to design a suitable feed horn with higher energy efficiency through a lower spillover.To do so, we will look into higher degree-of-freedom configuration possibilities (more than just binary control), as the requirement to achieve a specific reflection amplitude is beyond the capability of the binary phase control unit.We would also like to produce a prototype to test the performance of the device in practice.
XIAOMIN MENG received the B.Sc. degree in physics from the University of Rochester and the M.Sc.degree in theory and simulation of materials from Imperial College London.He is currently pursuing the Ph.D. degree with the Professor Nekovee's Group, University of Sussex.For his M.Sc.graduation thesis, he worked with Sir John Pendry on bragg-stacked graphene plasmonic metagratings, which has motivated him to pursuit further researches in the field of metasuraces and computational EM applied to wireless communications.

FIGURE 1 .
FIGURE 1.The schematics of unit cell design.v = 0.18mm, W = 0.714mm, h = 0.087mm.The patch antenna is positioned at the top, with a length of W .The LC is enclosed between the patch antenna and the ground plane.The thickness of the LC is h.The light green walls around the edges of the unit cell are the spacers designed to enclose the LC substrate.The biasing lines (depicted as tubes) function to create the capacitative effect between the patch antenna and the ground plane.

FIGURE 2 .
FIGURE 2. Conventions used for the far-field radiation pattern calculations.

FIGURE 3 .
FIGURE 3. Unit cell radiation pattern (A(θ, φ)) retrieved from infinite periodic conditions (only showing the nearest eight unit cells).The purple box is the simulation box, where the radiation pattern is recorded for the unit cell.

FIGURE 5 .
FIGURE 5.The 4-beam optimization convergence rates for the adaptive and non-adaptive GA algorithm.The adaptive GA converges (to an acceptable cost) roughly twenty generations faster, and the improvement is especially noticeable when an effective initial population is given.

FIGURE 6 .
FIGURE 6. Phase profile and directivity for θ = 15 • beamsteering from a f/d = 0.5 (5λ above) central feed source.The corresponding phase distribution for a 20 × 20 reflectarray with: (a) a continuous degree of change in the reflected phase (from 0 • − 360 • ) and (b) binary degree of change in the reflected phase (0 • or 180 • ), (c) directivity comparison.

FIGURE 9 .
FIGURE 9.The feed source horn antenna.(a) Far-field directivity in the polar plot at φ = 0 • with the red line, dark blue line for the main lobe direction, light blue for the 3 dB beamwidth (40.1 • ), and green circle for the SLL, (b) 3D radiation pattern and corresponding directivity in dBi.

FIGURE 10 .
FIGURE 10.Unit cell characteristics, reflection (a) phase and (b) amplitude of the ON and OFF states.

FIGURE 11 .
FIGURE 11. a) Reflection magnitude and b) reflection phase for oblique incidence.

FIGURE 12 .
FIGURE 12. (a) Theoretical and (b) full wave normalized far-field electric field.
MAZIAR NEKOVEE (Member, IEEE) received the bachelor's degree in electrical and electronic engineering from the Delft University of Technology, The Netherlands, and the Ph.D. degree in physics from the University of Nijmegen, The Netherlands.He is currently a Professor of telecoms and mobile technology and the Head of the Centre for Advanced Communication, Mobile Technology, and IoT, School of Engineering and Informatics, University of Sussex, U.K. Prior to joining Sussex in 2017, he was the Chief Engineer and the Head of 5G research at Samsung Research and Development, U.K. Prior to that, he was a Senior Scientist with British Telecom (BT) Research and Technology and subsequently a Team Leader.His research interests and expertise include 5G and beyond-5G/6G mobile communication RAN and core design, machine learning and AI applied to communication networks, spectrum sharing, cognitive radio, millimeter-wave, and THz communications.He has authored over 120 highly cited peer-reviewed papers, one best-selling book Cognitive Radio Communications and Networks: Principle and Practice and has 13 patents in telecommunication and mobile technologies.He is a sought-after speaker at high-level industry events, C-level meetings with European operators, policy makers, and international conferences.He was a recipient of a number of industry and academic awards, including the Royal Society (U.K. Academy of Science) Industry Fellowship, the BT Innovation Award, and the Samsung Research and Development Best Practice of Research Award.DEHAO WU received the bachelor's degree in optical and information engineering, the M.Sc.degree in microelectronics and telecommunication engineering, and the Ph.D. degree in optical wireless communications from Northumbria University, Newcastle, U.K., in 2007, 2009, and 2013, respectively.He is a Senior Lecturer at Bournemouth University (BU), U.K. Before joining BU, he held the role of a Sussex Research Fellowship and has been a Senior Research Fellow with the Department of Engineering and Design, University of Sussex, since April 2018.Before joining Sussex, he was a Postdoctoral Research Fellow at the Nanyang Technological University of Singapore and the University of Manchester.He has published more than 45 peer-reviewed papers in international conferences and journals, two book chapters, and two patents in indoor high-accuracy positioning technology.His current research interests include machine-learning algorithms, intelligent management of 5G and beyond-5G systems, artificial intelligence for 5G verticals, the IoT, optical wireless communications, free-space optics, visible light communications, RF communications, indoor high-accuracy positioning, optical sensing and detection, underwater RF and acoustic communications, information-centric networking, and information resilience.He is an Executive Committee Member of the IET RF and Microwave Technical and Professional Networks (TPN).He serves as a technical committee member for a number of IEEE conferences.

TABLE 1 .
LC permittivity and loss tangent values used for ON and OFF states.

TABLE 3 .
KPIs for different scan angles.