Accurate Direction-of-Arrival Estimation Method based on Space-Time Modulated Metasurface

A metasurface-based Direction of Arrival (DoA) estimation method is presented. The method exploits the properties of space-time modulated reflective metasurfaces to estimate in real-time the impinging angle of an illuminating monochromatic plane wave. The approach makes use of the amplitude unbalance of the received fields at broadside at the frequencies of the two first-order harmonics generated by the interaction between the incident plane wave and the modulated metasurface. Here, we first describe analytically how to generate the desired higher-order harmonics in the reflected spectrum and how to realize the breaking of the spatial symmetry of each order harmonic scattering pattern. Then, the one dimensional (1D) omnidirectional incident angle can be analytically computed using +1st and -1st order harmonics. The approach is also extended to 2D DoA estimation by using two orthogonally arranged 1D DoA modulation arrays. The accuracy of 1D DoA estimation is verified through full-wave numerical simulations. Compared to conventional DoA estimation methods, the proposed approach simplifies the computation and hardware complexity, ensuring at the same time estimation accuracy. The proposed method may have potential applications in wireless communications, target recognition, and identification.

is moving with respect to the other. Moreover, it is also needed for the next generations of mobile communication, such as 5G and beyond 5G (B5G), and the future sixth-generation (6G) communication, enabling location services for the mobile Internet [7] [8]. In this framework, array signal processing methods are the most common techniques to estimate DoA. Among them, it is worth mentioning techniques based on minimum variance distortion-less response (MVDR) [9], estimation of signal parameters via a rotational invariant (ESPRIT) [10] [11], multiple signal classification (MUSIC) [12] [13], and Bayesian compressive sensing (BCS) [14]. The high performances reached by these techniques are, however, counterbalanced by the need of massive arrays and the use of multiple sensors or channels, intense calculation, and expensive hardware. In this framework, a low in transmission and reflection, such as amplitude [15], polarization [16], propagation direction [17] [18], and, more recently, the frequency content [19]- [22] -cost technology and low-complexity DoA estimation system would represent an important advance in those applications where the angular localization of the transmitter, interfering signals, and users is mandatory.
In the last decades, metasurfaces have demonstrated to be a breakthrough technology in a number of applications, spanning from microwave to optical frequency ranges [23]- [25]. They typically consist of a bi-dimensional (2D) periodic array of subwavelength meta-atoms able to control the properties of the incident electromagnetic field. In particular, the control of the spatial and/or temporal characteristics of the interacting fields was made possible thanks to the development of dynamically reconfigurable metasurfaces, using different reconfiguring technologies based on pin diodes [26] [27], varactors [29]- [31], MEMS [32], graphene [33], and liquid crystal [34]. The control can be performed through analogical control signals [20]- [22] or via digital control [35]- [37]. In both cases, the resulting metasurfaces belong to the wider family of space-time metamaterials [38]- [45], which exhibit unprecedented capabilities and find applications in different operative scenarios, such as non-reciprocal structures [19], [46]- [49], arbitrary amplitude and phase programmable systems [50], wireless communications [52]- [54], and spread-spectrum radar camouflaging [55]. In this paper, we present and discuss a metasurface-based DoA method that exploits the properties of space-time modulated reflective metasurfaces to estimate in real-time the impinging angle of an illuminating monochromatic plane wave at frequency 0  . The system is composed by a space-time modulated metasurface and a detecting antenna placed in the far-field region just above it. The DoA estimation is here performed by measuring the amplitude unbalance of the reflected fields at the frequencies of the two first-order harmonics, i.e., 0 p    , generated by the interaction between the incident plane wave and the metasurface modulated at frequency However, we are unable to achieve the DoA estimation by just using the two first-order harmonics without any spatial difference. Due to the typical frequency used for the modulation, i.e. ωp<<ω0, in fact, the two harmonics exhibit an almost identical amplitude, being governed by the Manley-Rowe relations [56] [57]. Here, we exploit the spatial modulation as an additional degree-of-freedom for enhancing the required unbalancing, inducing a pattern asymmetry for the fields scattered at the two first-order harmonics. The proposed approach exhibits significant advantages with respect to the mentioned common array signal processing methods (MVDR, ESPRIT, MUSIC, BCS). Indeed, these systems require massive array antennas for achieving good performances. More recently, time modulated arrays with a limited number of antenna elements have been also proposed for DoA estimation [58]- [61] with the goal to reduce the number of antennas and overall computational complexity. However, the hardware of the backend network is still complex, involving properly connected phase shifters, switches, and power combiners. On the contrary, the metasurface-based DoA system shifts most part of the computational efforts at the electromagnetic level, reducing, thereby, the latency in detection and removing the need of a complex circuity [62]. Both the frequency modulation and the spatial unbalancing needed for estimating the DoA is here performed mainly at the metasurface level, whereas the remaining computation is extremely low-cost and consists only in the evaluation of the relative difference in amplitude of the two first-order harmonics received by the fixed detecting antenna. In the next sections, we describe the analytical model used for the DoA estimation and demonstrate through a proper set of numerical simulations the high accuracy reached by this method.
The paper is organized as follows. After having illustrated the basic idea, in Sec. II, we derive analytically the temporal modulation scheme to apply to a space-time modulated metasurface for generating the two desired first-order harmonics and the spatial modulation scheme for achieving the breaking of the spatial symmetry of each order harmonic scattering pattern. In Sec. III, the analytical DoA estimation method based on the double-sideband (DSB) space-time-modulated metasurface is derived for both 1D and 2D case. In Sec. IV, a realistic metasurface is used for demonstrating numerically the accuracy of the proposed method for DoA estimation. Different configurations have been tested, showing the operative bounds of the proposed system. Finally, in Sec. V, a short conclusion of our work is given.
II. SPACE-TIME MODULATED METASURFACE FOR DOA ESTIMATION In this Section, we present the proposed DoA estimation method based on space-time coding modulated metasurfaces. After having introduced the basic idea, we analytically describe the realization of the DSB time modulation by using a 1-bit coded metasurface, which allows exciting the required 1 st  harmonics in reflection. The breaking of the spatial symmetry of each order harmonic scattering pattern is achieved by applying a delayed copy of the modulating signal to adjacent portions of the space-time modulated metasurface.

A. Basic idea of DoA estimation
Let us consider the operative scenario reported in Fig. 1. A space-time modulated metasurface is illuminated by a monochromatic plane wave at frequency 0  propagating with an angle  with respect to the positive z-direction.
We assume that the metasurface-based DoA system operates only on the xz-plane. This constrain will be relaxed later. Regardless the technology used for imparting the space-time modulation on the metasurface, when its surface properties are temporally modulated with frequency ωp<<ω0, a frequency mixing between the incident and modulating signals takes place, spreading the energy of the incident field over a number of frequency harmonics located at 0 . The detecting antennas placed in the far-field region just above the metasurface can now receive a different amplitude contribution at the two frequencies ω0±ωp and estimate the DoA of the incident plane wave based on the relative unbalancing of such values. This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2022.3184556

B. DSB modulation in time coding metasurfaces
The main purpose of this part is to show analytically how a time-varying reflective metasurfaces can generate the required harmonics in the reflected spectrum for applying the proposed method for DoA estimation. Since we are here interested in the temporal modulation of the incident wave, we can relax the constrain on the oblique illumination condition and assume that the metasurface is normally illuminated by the plane wave. In the next sections, we will consider also the effect of the oblique incidence. Figure 2. Double-side band modulating metasurface illuminated by a normally incident plane wave able to reflect a plane wave consisting of a superposition of the two first-harmonics.
In Fig. 2, we report an ideal non-penetrable (i.e. characterized by a zero-transmission) reflective metasurface, whose macroscopic response can be modeled through its . The metasurface is  . This implies that half of the overall impinging energy is distributed between the two first order harmonics at 0 whereas the rest is necessary dissipated by the intrinsic reflective response of the metasurface with cosine temporal profile due to the zeros of the cosine function. Despite the cosine temporal profile of the reflection coefficient ensures theoretically the required response from the metasurface, such a temporal profile needs a continuous modulation of the surface properties of the metasurface, as demonstrated in [20] [21], which is not easy to achieve.
(a) (b) Figure 3. (a) Temporal signal with cosine profile to be used for achieving a perfect double-side band excitation at frequencies approximate three state signal used to match the curves in (a).
To simplify the metasurface implementation, the cosine temporal profile of the reflection coefficient can be approximately replaced with acceptable tolerance by a train of rectangular pulses [63], as shown in Fig. 3. The time-period Tp is the modulation period corresponding to the angular frequency 2 . The rectangular pulses of amplitudes "+1" and "-1" in Fig. 3(b) correspond to the reflection phases "0" or "π", respectively. Now, the start time t1 of the "+1" state, the start time t2 of the "-1" state, and the pulse duration τ characterizing the shape of the rectangular pulses must be analytically derived to match the cosine time profile as better as possible [63], and ensure the excitation of only the two first-order harmonics. For the first time-period p T starting at 0 t  , the rectangular time sequence can be expressed as: The periodic time sequence is decomposed by using the Fourier series to highlight the contribution of the different frequency harmonics as: where h is the order of the specific harmonic, and h a is the corresponding Fourier coefficients, whose expression is evaluated as follows: Considering the periodicity of the modulation, we can normalize the temporal axis with respect to the time-period p T and obtain the dimensionless quantities: Eq. (7) can now be easily evaluated as follows: This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication.
Considering the term 2 1 sin , and so on. This configuration ensures the highest isolation between the first order harmonic and the other order harmonics achievable with a binary time modulated metasurface. However, it is worth noticing that, to exactly match the temporal profile shown in in Fig. 3(b), an absorption state is required between two consecutive opposite pulses.
To further simplify the implementation, we remove the possibility to have an absorption state, setting 0.5    . This leads to a time-varying reflection coefficient that can exhibit only "+1" and "-1" reflection states, which corresponds to the reflection phases "0" and "π", respectively, as shown in Fig.  4(a), over a time-period p T . Consequently, the presence of the , the complete set of the h a coefficients is: where sgn(h) is the sign function. In Fig. 4(b), we report the amplitude of the harmonic coefficients as a function of the harmonic order. It can be found that the amplitudes of the first-order harmonics are still three times bigger than the other order harmonics, letting our approach to still benefit of the presence of strong 1 st  -order harmonics. In the next Section, the required spatial modulation is imposed on the metasurface in order to achieve the splitting of the radiation patterns of the reflected 1 st -order harmonics, which allows estimating the DoA of the incident wave.

C. Pattern unbalancing in DSB modulated metasurface
Once the proper temporal modulation of the 1-bit space-time metasurface is properly defined (see Section II-B for further details), the spatial asymmetry must be introduced for splitting the radiation patterns of the 1 st  -order harmonics and enable the DoA estimation capability of the system. The metasurface is partitioned in areas of width D in x-direction, called sub-macrocells, each of which covers the entire extension along y-direction, including Nx  Ny lines of space-time modulated inclusions, as shown in Fig. 5(a). The spatial asymmetry is here introduced by applying different modulation signals to adjacent sub-macrocells, i.e., "Sub-macrocell 1" modulated by the signal 1 ( ) u t , and "Sub-macrocell 2", modulated by the signal 2 ( ) u t , reported in Fig. 5(b) and 5(c), respectively. The two signals are identical, except for a minor time delay t  . They compose the DoA-MTS macro unit-cell responsible for the space-time asymmetry in the reflected patterns. It is worth mentioning that, due to the uniformity along y-direction, the scheme in Fig. 5(a) can perform DoA estimation only for plane waves propagating on the xz-plane.
Since the proposed method is valid regardless the specific technology and subwavelength inclusion composing the metasurface, in the following we describe its scattering pattern in terms of the Array Factor. The complete scattering response can be easily obtained by multiplying the array factor we derive in the following and the radiation pattern of the single subwavelength inclusion composing the metasurface [64]. The array factor of the single DoA-MTS macro unit-cell for the h-th harmonic is indicated as   (1) 1 This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2022.3184556 In (12), Substituting Eq. (12) into the first double-summation term of Eq. (11) we obtain: As shown in Fig. 5(a), the entire DoA estimating metasurface consists of many DoA-MTS macrocells, say M, and, thus, the full array factor describing the scattering pattern of the entire metasurface is: where M is the number of DoA-MTS macro-cells. We can notice now an interesting term in Eq. (13), and thus also in Eq. (14), that allows us enabling the spatial asymmetry in the scattered patterns for the two first-order harmonics of the reflected field: depending on the value of

A. 1D DoA estimation method based on DSB metasurface
Let us start deriving the DoA estimation equations in the case of illuminating plane waves, whose wavevectors lie on the xz-plane, i.e., 0 inc   . The goal is to estimate the performances of the proposed system in identifying the angle inc  of the incident plane wave, i.e., 1D DoA estimation.
In this scenario, the uniformity along the y-direction relaxes the need to consider the other scattering directions besides the ones laying on the xz-plane. Setting, thus, 0   , and observing the metasurface in broadside ( 0)   , where the detecting antenna is located, the amplitude ratio between 1 Using trigonometric angle sum and difference identities to expend the right hand side of Eq. (15), and it can be simplified to obtain the DoA estimation curve as: from which, we can estimate the DoA angle inc  as: This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication.
To avoid angle ambiguity when Eq. (17) is used for DoA estimation, the curve should be monotonic among the estimation range, thus the conditions for no angle ambiguity are expressed as: In the range of direction estimation , the sine function sin inc  is real and within   (18), to avoid angle ambiguity, the conditions can be written as: Eq. (19) defines a bound for the time delay between the modulating signals to be applied for a given electrical dimension 0 D  of the Sub-macrocell. However, it is worth remembering that, due to the periodicity of the modulating signals 2 ( ) u t , the time delay is always within a period, i.e., This sets also an upper-bound for the sub-microcell dimension D, that cannot exceed the half-a-wavelength width.
As to chromatic plane waves, the proposed approach is still effective as long as the bandwidth B of the incident signal is smaller than the modulation frequency fp = 1/Tp (fp < B). If fp ≥ B, the harmonic spectra will not overlap, and, thus, we can use the harmonic bands to estimate the DoA like in the case of the monochromatic incident wave. Though theoretically fp can be large enough to obtain wideband DoA estimation, due to hardware limitation, in the performed experiments [45] fp can reach only the MHz. For practical linear frequency modulation (LFM), incident waves (fp < B), in [59], the LFM signal with long duration can be separated into multiple short durations, and the bandwidth of each short duration is smaller than the modulation frequency. Thus avoids overlapping and relaxes the limitation of LFM bandwidth. In [60], the pulse compression technology is used to calculate the modulation harmonics, the harmonic coefficients of the modulated signals can also be obtained by the output of the matched filter.

B. 2D DoA estimation extension
The above method can be also extended to 2D DoA estimation. The 2D detection scheme requires that the scattering pattern splitting of the 1 st -order harmonics of the reflected field is introduced for both x-and y-directions of the metasurface to estimate both polar coordinates   This can be obtained by modulating the metasurfaces separately along the two main axes, x-and y-directions. For example, the 1D modulation scheme reported in Sec. II-C can be applied in sequence along the x-and y-directions by properly controlling in real-time the modulation profile on the metasurface or, alternatively, it is possible to use two identical DoA metasurfaces arranged orthogonally, as shown in Fig. 6.
(a) (b) Figure 6. (The 2D DoA estimation arrays (a) the first array along x-direction, (b) the second array along y-direction.
Following the same procedure illustrated in Sec. III-A for the 1D DoA estimation, we can derive the estimating function for the spatial modulation along the x-direction (Fig. 6(a)) considering an incident plane wave impinging from the where the detecting antenna is located, we can write: that returns: Moreover, the estimating function for the spatial modulation along the y-direction ( Fig. 6(b)) is: that, in turn, returns: From Eq. (21) and (23), inc  and inc  can be derived in closed forms as: The 2D DoA estimation results can be expressed as Eq. (26), when θ is positive, θ = θ, ϕ = ϕ, but if θ is negative, θ = -θ, ϕ = ϕ+π.
This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. The signal-to-noise ratio (SNR) is set as 10dB with respect to the amplitude of the total reflective field from the metasurface. For each incident angle, a 1000-times Monte Carlo simulation is used to calculate the MSE of the estimation. The simulation results are shown in Fig.  7(a), where it can be found that the error is lower than 0.4°over the entire field-of-view. Much lower error is found near ±60˚, whereas relative larger errors are found near 0˚, ±90˚, due to the changes of the slope of the estimation curve as explained above in Fig. 7(b).
The accuracy of the DoA estimation is here discussed by comparing the proposed method and the conventional MUSIC algorithm. To keep the same conditions of the Monte Carlo simulation, we use the same array size, center frequency, sampling frequency, SNR, etc. In the MUSIC method [12] with two omnidirectional antennas (dimensions 2×0.44 λ), 5000 continuous data points of time domain far-field are sampled from each antenna. In the MUSIC method, thus, 10000 (5000×2) data points are sampled in total. In the proposed method, the modulation frequency is set as  [66]. in our simulation, when the PIN diode is at the ON state, it can be equivalent as a series circuit of parasitic inductance (1.61 nH) and resistance (0.7Ohm); if the PIN diode is at the OFF state, a series circuit of parasitic inductance (1.39nH), capacitance (0.41pF), and resistance (1.41Ohm) is used to describe the physical model. In addition, a Floquet port is used to produce y-polarized waves incident onto the meta-atom and receive the reflected waves. Periodic boundary condition is set to its four sides to model the infinite array. It exhibits the binary 0- phase response at 3.3 GHz, with a reflection amplitude larger than 0.82 over the angular range [0˚-70˚] of the incident angles (Fig. 8(c)-(d)).

A. 1D DoA estimation using different widths of the DoA macrocell
To demonstrate the DoA estimation using the proposed DSB space-time-modulated metasurface, a simulation is implemented considering a metasurface with one DoA-MTS microcell, composed by two sub-macrocells, each of which including 4 x N  and 10 y N  unit-cells shown in Fig. 8(b) along the x-and y-directions, respectively. The operative central frequency is 0  The normalized scattering patterns of +1 st -and -1 st -order harmonics are shown in Fig. 9(a), and the DoA estimation curve can be calculated from Eq. (16) by the scattering patterns of +1 st -and -1 st -order harmonics. In Fig. 9(b), we show that the analytical and numerical curves agree very well except for the boundary of the incident angle range.
The DoA estimation results are summarized in Table Ⅰ  It is interesting to also consider a narrower field-of-view and derive the corresponding modifications to the modulation scheme at the metasurface level. For example, limiting the field-of-view to the range  

,30
   , using Eq. (19) we obtain that DoA estimation can be done using a larger width of the sub-macrocell. In particular, using again 0.05 p t T   , the width D can be doubled, allowing to include 8 x N  unit cells along the x-direction and 10 y N  unit cells along the y-direction within a single sub-microcell. The metasurface, therefore, is composed by 10×16 unit-cells.
The normalized scattering patterns of +1 st -and -1 st -order harmonics and the DoA estimation curve are reported in Fig. 10.
In the angular range of interest, the curve is monotonic and single-valued, but it is steeper than the one in Fig. 9(b). This implies an enhancement of the estimation accuracy, as shown in Table Ⅱ.  In practical implementation, there are two ideas to reduce the blockage effect of the receiving antenna: i. miniaturize the detecting antenna to reduce the aperture; ii. place the detecting antenna with a certain offset angle [67], to reduce the blockage to the incident wave from the angles we are interested in. Moreover, 1-bit refracting time-modulated metasurface can be introduced to implement the proposed method to avoid the occlusion effect [68].

B. 1D DoA estimation using different time delays
According to Eq. (16), the DoA estimation curve is not only related to the dimensions of the DoA element, but also to time delay t  of the modulating signals 1 ( ) u t and 2 ( ) u t applied to the two sub-macrocells of the DoA-MTS microcell. To demonstrate the impact of the time delay in the DoA estimation accuracy, a simulation is implemented considering a metasurface composed by just one DoA-MTS microcell, whose sub-microcells include

unit-cells along the
x-and y-directions, respectively.
In Fig. 11, we report the DoA estimation curves for different time delays t  . The light blue curve represents all the results for any value of 0.05 p t T   that satisfies Eq. (19). In this case, the curve is monotonic, and there is no ambiguity on the DoA. On the contrary for larger time delays (dashed curves in Fig.  11), the operational field of view shrinks according to Eq. (19), where the angle  must be smaller for satisfying the constrain. The DoA estimation accuracy is obviously affected by the error in time delays, being this quantity vitally important to the accuracy of the presented method. To analyze the effect of potential errors on time delays, we consider here the delay of periodic temporal modulations as Δt=0.05Tp, whose corresponding DoA estimation curve is shown in Fig.12 (a) (blue dotted line). If there are errors in time delays within the range of ±20%, i.e., between Δt = 0.04Tp (Δt error is -20%) and Δt = 0.06Tp (Δt error is 20%), the corresponding estimation curves shown in Fig. 12(a) are modified, with an incorrect range and slope of the DoA estimation curve. This will cause an error in the DoA estimation, correspondingly. For example, when θinc =30°, the DoA estimation error due to the Δt error is shown in Fig. 12(b). As expected, the DoA estimation error increases with the Δt error.
C. Analysis 1D DoA estimation error due to phase error For eq. (17), which is derived from eqs (13) and (15) under the condition of two ±1 harmonics exhibiting an identical amplitude, we just use ±1 harmonics from the spectrum to achieve the DoA estimation. Eq. (17) is still suitable for modulation methods with ±1 harmonics of identical amplitude. The different temporal modulations just affect the purity of ±1 harmonics that we substitute into eq. (17), but have no influence on the accuracy of DoA estimation theoretically. As shown in Fig. 13, when 1-bit temporal modulation (see in Fig.  4(a)) with different phase errors, the ±1 harmonics are with identical amplitude while having different efficiencies. The DoA estimation curves keep unchanged, namely the phases errors of 1-bit have theoretically no influence on the accuracy of the DoA estimation. When increasing the phase errors, the efficiencies of the ±1 harmonics decrease, and the anti-noise ability decreases correspondingly. Please, note that the DoA estimation method still works when the modulation generates two ±1 harmonics with different amplitude. However, in this case, eqs. (15) and (17) must be rederived starting from eq. (13) that must take into account the new modulation scheme. To better clarify this point, let us consider a modulation scheme that introduces a phase error generated during the time delay between u1(t) and u2(t), because of the different boundary conditions of the sub-macrocells (as shown in Fig. 14(a)). This error will lead to a different amplitude of ±1 harmonics, as shown in Fig. 14(b). The corresponding DoA estimation curves are shown in Fig.  4(c), clearly demonstrating that they do not overlap anymore, leading to an error in the estimation. In particular, when θinc =30°, the DoA estimation error due to phase error during the time delay is shown in Fig. 14(d). Obviously, the DoA estimation error increases with the phase error.

V. CONCLUSION
In this paper, we have proposed a DoA estimation method based on space-time modulated metasurfaces that exploits their capabilities to generate several harmonics in the scattered field, each of which is radiated towards different direction. By detecting the field in the broadside direction by using a single receiving antenna the metasurfaces-based DoA system can estimate in real-time the impinging angle of an illuminating monochromatic plane wave. The proposed approach exploits the amplitude unbalance of the received fields at broadside at only the two first-order harmonic frequencies generated by the interaction between the incident plane wave and the modulated metasurface. The analytical model for estimating the incident angle has been presented and successfully used in combination with a realistic 1-bit reflective metasurface. The accuracy of 1D DoA estimation is verified through full-wave numerical simulations. Compared to conventional DoA estimation methods, the proposed one simplifies the computation and hardware complexity, ensuring a very good estimation accuracy. The proposed approach may have potential applications in wireless communications, target recognition and identification. In the last ten years, his main research interests have been focused on the analysis and design of cloaking metasurfaces for antenna systems, on the modeling and applications of (space and) time-varying metasurfaces, on the topological-based design of antennas supporting structured field, on the modeling, design, implementation, and application of reconfigurable metasurfaces, on the concept of meta-gratings and related applications in optics and at microwaves, on the modeling and applications of optical metasurfaces. The research activities developed in the last 20 years has resulted in more than 500 papers in international journals, conference proceedings, book chapters, and three patents.
Prof This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2022.3184556