Mutual Coupling Aware Time-Domain Characterization and Performance Analysis of Reconfigurable Intelligent Surfaces

Reconfigurable intelligent surfaces (RISs) are planar, almost 2-D, structures that can intelligently manipulate electromagnetic waves by low-cost near passive reflecting elements. RISs are considered a potential enabling technology for the sixth-generation (6G) wireless communication systems due to their capability to tune wireless signals, thus smartly controlling propagation environments. This work has a twofold objective: firstly, we present a systematic methodology for characterizing RISs in the time domain, which rigorously includes propagation delays and mutual coupling among RIS elements; secondly, we analyze the RISs through a convolution-based solver that is enriched by a specialized convolution scheme that enables the analysis of systems with complex linear, nonlinear as well terminations. The proposed approach is validated by extensive numerical evaluations against well-established methods based on both frequency- and time-domain analyses, with a varying number of RIS elements, interdistance, and frequency of operation.


I. INTRODUCTION
R ECONFIGURABLE intelligent surfaces (RISs) repre- sent one of the key enabling technologies for the nextgeneration wireless communication systems, making it possible to create an intelligent electromagnetic (EM) environment [1], [2], [3], [4].Specifically, RISs are almost-planar structures engineered to dynamically control reflections and refractions of impinging EM waves.Therefore, the use of RISs allows making the wireless propagation environment part of the network design parameters that can be optimized in order to gracefully integrate RIS with fundamental enabling technologies such as multiple-input and multiple-output (MIMO), and nonorthogonal multiple access (NOMA) [5], [6], [7], in both indoor and outdoor applications [8], and also in nonterrestrial networks [9].
Typically, a RIS consists of a large number of nearly-passive scattering elements and behaves as a reflectarray with subwavelength interdistances between the array elements.This feature makes RISs modeling and design extremely challenging, and thus, a systematic characterization of such systems becomes crucial [10], especially from an EM standpoint.From the modeling perspective, the system comprising the transmitting antennas, the RIS, and the receiving antennas can be regarded as a multiport system.Indeed, due to the closely spaced scattering elements, the mutual coupling among them cannot be ignored [11] and needs to be taken into account since the design stage to guarantee optimal performances, e.g., super-directive designs in which the gain scales with the square of the number of antenna elements, provided that the Ohmic losses are sufficiently low [12].
To consistently integrate the RIS paradigm into wireless networks, the availability of an end-to-end communication model that is EM-compliant is essential, and this requires a proper description of the mutual coupling that can be accomplished by relying on accurate modeling methods.For instance, in [13], rigorous EM modeling of a RIS has been presented, which permits explicitly to take the mutual coupling among the scattering elements into account.The transmitting antennas, the receiving antennas, and the passive scatterers are cylindrical thin wires of perfectly conducting material whose length is and whose radius a is finite but negligible compared to (thin wire regime).Despite its simplicity, this model is general enough for the characterization of passive scatterers and spatial distributions, provided that the radiating elements are minimum scattering antennas [14].In [13], the transmitting, receiving, and RIS dipoles are modeled in the frequency domain (FD) considering the assumption of minimum scattering radiating elements and thin wire regime, moreover, approximating the currents with sinusoidal functions as in [14]; the resulting model is then coupled to the circuit model of the terminations assumed as linear, and the analysis is completed in the FD as well.This approach is improved in [15], [16] where the end-to-end channel transfer function matrix is formulated in a closed-form expression in terms of the exponential integral function, greatly simplifying its computation.More recently, such a methodology has been extended to include also scatterers modeled as loaded dipoles and used to optimize the RIS in [17].
Recently, an FD macromodeling approach of RISs has been presented in [18], where the microwave network theory is used to compute the accurate reflection coefficient of RISs elements under the illumination of incident waves from arbitrary directions.Hence, in this work, only the RIS is modeled and transmitters and receivers are not considered in the overall model.
Different from the other works presented in the literature that are based on FD approaches to describe RISs, our goal is to present a time-domain (TD) characterization of RIS systems.Indeed, a direct TD description of linear time-invariant (LTI) EM systems may satisfy the causality principle in a more rigorous way than FD techniques [19], [20].Furthermore, when the EM systems have to be interfaced with nonlinear terminations, a TD model becomes mandatory, and this may be the case of more complex RIS structures based on a planar array of unit cells different from dipoles [21].A TD model is also desirable for time-varying (TV) scenarios such as urban ones with moving scatterers and antennas or when TV materials are used for the RIS.It is also to point out that, when using FD techniques, a further step is required to restore a TD model of the system.This can be done by relying on the inverse fast Fourier transform (iFFT), which is very sensitive to sampling and typically requires a large number of samples whose processing may result in cumbersome computation.Another option is represented by state-space or circuit synthesis techniques, which require the identification of poles and residues of the transfer functions [22].Nevertheless, this kind of approach still requires a fine sampling of the transfer functions resulting in a computational burden.Furthermore, this type of technique requires the identification of propagation delays, which are important in the case under consideration, and, in any case, it is not compatible with TV scenarios.
A preliminary effort of applying TD characterization has been conducted in our previous work [23], where the RIS has been analyzed through the partial element equivalent circuit (PEEC) method [24].Specifically, the TD scattering parameters have been computed by taking the mutual coupling among the RIS elements into account.Moreover, the propagation delay between the transmitter, the RIS elements, and the receiver has been rigorously modeled by considering a near-field scenario.Then, the port impulse response has been exploited within the framework of a convolution-based solver by assuming ohmic-inductive terminations, in accordance with the model provided in [13].

A. Article Contributions and Organization
By leveraging on the approach presented in [23], in this article, we provide the following contributions.
1) We propose a comprehensive TD model that extends the methodology presented in [23] by considering a more general description of the transmitters-RIS-receivers structure based on a state-space formulation.This allows for a more systematic computation of the TD impulse response of the RIS through the PEEC method; since transmitting and receiving antennas are modeled along with the RIS, propagation delays are rigorously considered.Hence, the application of the PEEC method results in an equivalent circuit with retarded controlled sources.2) We extend the convolution-based solver to consider more general linear terminations.When the transmitting antenna is excited and the RIS and receiving antennas are loaded by passive terminations, the previously computed S-parameters model is used to compute the port voltages.3) Since convolutions can be time-consuming to be computed, we propose a rationale to accelerate the solver, by resorting to the so-called segment fast convolution (SFC) method described in [25].4) With respect to the methodology presented in [25], we explore a further acceleration of the convolution-based solver that is obtained by exploiting the potential of matrix-vector computation and allows us to substantially reduce the computation time required for the RIS characterization.5) We validate our methodology via extensive numerical evaluations against well-established methods based on both frequency and TD analysis.It is to be remarked that, although the methodology is applied to RISs constituted by dipoles, it can be applied to any type of RISs whose elementary scatterers can be more complex.Furthermore, although in this work, we will consider LTI passive terminations, the proposed approach can be extended to consider nonlinear and TV terminations, which will be considered in forthcoming reports.
The rest of the article is organized as follows: the system model is outlined in Section II by using the impedance matrix multiport representation of the system comprising the transmitter, the RIS, and the receiver.Furthermore, the figures of merit of the end-to-end channel gain, H E2E , are the line-of-sight (LOS) link and the virtual LOS (VLOS).Section III presents the TD PEEC solver that is used to compute the TD S-parameters of the multiport system.The time-stepping methods that can be used to solve the delayed differential equations are also described in Section III.The convolution-based solver with arbitrary linear terminations is presented in Section IV.Section V briefly summarizes the SFC method and investigates the impact of the breakpoint selection on the accuracy of the method.The proposed methodology is validated in Section VI by several tests.Finally, Section VII concludes the article.

II. SYSTEM MODEL
We consider a wireless communication system as in Fig. 1, in which a transmitter (e.g., a base station) and a receiver (e.g., a mobile user) communicate directly and through a RIS that is made of N ris elements.
Following the model provided in [13], each element of the system can be thought of as a dipole antenna.So, the transmitter could be one, or more, dipoles fed by an independent voltage source V s (t).Likewise, the receiver could be one, or more, dipoles connected to load impedances Z R .Then, the RIS is assumed as an n × n dipole array with each element terminated on a tunable load Z .Indeed, it has been proved in [16], [26], and [27] that the loaded wire dipole model, besides offering tractability and EM consistency, is also well suited for modeling the EM scattering from arbitrarily shaped objects.
Port voltages and currents V and I of the entire system including the transmitters, the RIS, and the receiver, are related by the impedance matrix where Z is where subindexes T, S, and R refer to the transmitter, RIS, and receiver, respectively.Furthermore, we denote with Z T , Z ris , and Z R the matrix impedances of the port terminations.
We denote as end-to-end channel gain H E2E the channel gain that accounts for the transmit and receive antenna elements, the passive scatterers of the RIS, the load impedances and, more importantly, the tunable load impedances of the RIS.The endto-end channel gain can be thought of as the algebraic sum of the contribution of the LOS link and the VLOS due to the presence of the RIS.We denote LOS channel gain, H LOS , the contribution to the channel gain in the direct link and VLOS channel gain, H VLOS , the contribution to the channel gain in the secondary link that includes the RIS.The end-to-end channel gain can be formulated as where The impedance matrix Z can be easily recovered by the S-parameters representation through a standard transformation [28].It is common practice to use the S-parameters representation either because they can be measured or also because they are less sensitive to numerical issues when computed through a numerical simulation.From (3), it is clear the end-to-end channel gain is given by two main contributions: 1) the LOS channel gain, defined as H LOS = y 0 Z RT and 2) the VLOS channel gain, that is equal to A more detailed derivation and description of the discussed quantities can be found in [13].

III. PEEC-BASED S-PARAMETERS IMPULSE RESPONSE
In a typical wireless communication system, the involved elements (e.g., transmitter, RIS, and receiver) are usually configured to be very far apart.Consequently, propagation delays become significant and they should be properly taken into account.In [13], a discrete dipole approximation has been used to compute the figures of merit described in the previous section.Although the FD modeling is relatively simple to be implemented, it prevents considering TV scenarios involving moving scatterers or TV materials, which would require a TD approach.This suggests adopting a TD methodology, which allows for solving the same problem by surmounting the limitations of FD techniques.This can be achieved by using the PEEC method [24] for the EM analysis in the TD.In particular, the integral form of Maxwell's equations in the TD is solved by resorting to the methodology described in [24] and [29], where the enforcement of Kirchhoff voltage and current laws to the PEEC model yields the following set of neutral delayed differential equations [24]: The vector of the unknowns x(t) ∈ n u ×1 is given as where i(t) are the branch currents, φ sr (t) are the scalar potentials for surface nodes, φ i (t) are the scalar potentials for internal nodes, v d (t) are the excess capacitance voltages for dielectric branches, and q s (t) represent the surface charges.Furthermore, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the state space matrices C(t), G(t), and B are where * represents the convolution operator, and n b , n ns , n ni , n bd , and n p represent the cardinality of branches, surface nodes, internal nodes, dielectric cells, and surface cells, respectively.Moreover, C d is the excess capacitance matrix [30], R is the branches resistance matrix, A s is the incidence matrix for the surface nodes, A i is the incidence matrix for the internal nodes, Γ is the dielectric region selection matrix, M is the surface-tonode reduction matrix, and G le is the load conductance matrix (assuming by now for simplicity that only resistive lumped elements are connected to the PEEC model).
Also, the time-dependent partial inductance matrix L p (t) and the coefficient of potential matrix P(t) are considered as impulsive, that is where the superscripts DL and D stand for delay-less and delayed, respectively; N τ Lp is the number of significant delays between elementary volumes while N τ P is the number of significant delays between elementary surfaces; τ cc,i = R cc,i /c 0 , i = 1, . . ., N τ Lp and τ cc,q = R cc,q /c 0 , q = 1, . . ., N τ P denote the delays between the centers, identified by R cc,i and R cc,q respectively, of the spatial supports of the basis functions of currents and charges; c 0 is the speed of the light in the background medium.
Finally, the source vector u(t) is given as where v s (t) and i s (t) are the voltage and current sources, which are applied to branches and nodes, respectively.It is also worth observing that the proposed PEEC formulation, which is based on charge and current basis functions, has been found to be more robust towards the so-called low-frequency breakdown, as confirmed in [31] and [32], resulting to be a full spectrum methodology.More refined modeling approaches of time-dependent zero-thickness partial inductances and coefficients of potentials are presented in [33], [34], and [35]; nevertheless, they are not considered in this work, and their application to our problem is left for future investigations.
In the following, the PEEC method is used to compute the impulse response of the multiport system represented by the transmitter, the receiver, and the RIS.This task is performed by computing the TD scattering parameters s(t).This is done in the following four steps.
1) All the ports are terminated on a reference resistance R 0 , typically 50 Ω.
2) The PEEC model is excited at every single port with a step current source i s (t) with a small rise time and the voltages v(t) and currents i(t) are computed at all the ports.3) The incident and reflected waves at the ports, a(t) and b(t), respectively, are computed as 4) The TD scattering parameters s(t) are finally obtained by the time derivative of the reflected waves b(t) that can be regarded as the step response of the system in the wave domain.The steps for the computation of the TD S-parameters can be summarized, as shown in the flowchart in Fig. 2. It is important to remark that the choice of the rise time of the step current source is to be consistent with the bandwidth of the system application.The bandwidth of interest for RIS applications is 0-60 GHz, which implies that the rise time should not exceed a few picoseconds.
As already mentioned, TD scattering parameters s(t) are preferred to other possible multiport representations, e.g., the impedance or the admittance multiport ones, because of their milder boundary conditions, represented by resistances rather than open-circuit or short-circuit terminations.Indeed, this choice makes the computation of port voltages and currents less prone to numerical issues than the open-or short-circuit conditions that are to be enforced to compute the impedance or admittance models.This allows using safely the parameters s(t) within the framework of a convolution-based TD solver, as described in the next section.

A. Time-Stepping Analysis
The PEEC method, summarized in the previous section, leads to a set of delayed differential equations involving convolutions representing the magnetic and electric field couplings.
Let us consider (5), which, depending on the delays and the choice of the time step, can be recast separating the interactions that can be considered instantaneous from those that can be handled as known because took place in the past where the state space matrices C DL (t), C D (t), G DL (t), and 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Quantities C D (t) dx(t) dt and G D (t)x(t) are assumed as known because they depend on the state vector x(t) and its derivative dx(t) dt evaluated in the past due to delays τ cc,i , i = 1, . . ., N τ Lp and τ cc,q , q = 1, . . ., N τ P .Then, (13) can be expressed by moving the known quantities on the right-hand side which, can be rewritten in a more compact form as where One and two-steps integration methods can be represented in the following general form for the case where uniform time steps are used: where the subscript p identifies the time step and homogenous initial conditions are set for x p−1 , x p−2 , and . The coefficients for each of the methods are given in Table I.
Applying the scheme (20) to (19), we obtain Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where K p is the time discrete counterpart of K in (19) and is expressed as where m i and n q represent the number of time steps corresponding to delays τ cc,i and τ cc,q and the sampling property of the Dirac function in ( 15) and ( 17) has been exploited.It is important to remark that, assuming a time step smaller than the smallest time delay, K p in ( 22) is computed from the past values of x(t) and d dt x(t), and thus, it is moved to the right-hand side of (21).Equation ( 21) can be rewritten in a more compact form as where A reads and has to be updated at each time step.

IV. CONVOLUTION-BASED TD SOLVER
As stated in Section I, the transmitter (Tx), the receiver (Rx), and the elements of the RIS can be thought of as a multiport system with N p = N T x + N Rx + N RIS ports.As explained in Section II, the multiport system is characterized in the TD in terms of its scattering parameters s(t).
Through (12), the port voltages v(t) and currents i(t) can be expressed in terms of the incident and the reflected waves as where R 0 is the reference resistance.The reflected waves b(t) are computed by the convolution of the incident waves a(t) and the S-parameters s(t), as Furthermore, the Kirchhoff voltage law (KVL) at the N p ports is enforced as where v s (t) is the vector containing the port voltage sources, and for the problem considered, it will have nonzero entries only for the N T x ports.While ṽ(t) represents the voltages across the passive part of the terminations, which are assumed to be LTI systems.Without loss of generality, they admit a state-space representation where x(t) is the state vector of the passive part of the terminations.Equations ( 24), ( 25), (26), and ( 27) can be recast as which is a well-posed problem with unknowns a(t), b(t), ṽ(t), and x(t).The discrete-time domain counterpart of (28), at the kth time step, is shown in (29) shown at the bottom of next page, where h is the time step and μ is a general coefficient of the integration method.The coefficients of the most commonly used methods are reported in Table I.The time discrete port voltages and currents are obtained from (24) as is the contribution of the reflected waves at the ports up to the k − 1 time step.
In brief, once the TD characterization of the RIS is available as described in the previous section, the multiport system with different port terminations or different sources can be studied through (29) using a time-stepping approach.
The main advantages of the proposed approach are as follows.
1) The TD computation of the impulse response is in general much faster than the FD one that requires the solution of a dense system of equations over a wide and finely sampled frequency range.2) The impulse response is evaluated only once at the beginning and then reused with different terminations using the convolution-based solver (28), which returns the port voltages and currents.3) As already mentioned, the proposed approach validity is not limited to LTI terminations considered in this work but can be easily extended to consider nonlinear and TV terminations.On the other side, the step response from which the impulse response is finally recovered by the time derivative is computed using a steep step and this causes an approximation at a very high frequency.Hence, the rise time of the finite step and the time sampling need to be carefully chosen depending on the frequency range of interest.

V. ACCELERATION OF THE CONVOLUTION
The computation of the convolution in (25), can be accelerated by resorting to a piecewise linear (PWL) representation of the step response.As shown in [25] and [36], this entails a piecewise constant (PWC) approximation of the impulse response.In this section, the PWLFIT technique is briefly revised.Then, an extended version of this approach is illustrated, and finally, it is shown how this approximation method is used to accelerate the convolution.

A. PWL Approximation
The PWLFIT technique [25] allows computing the PWL model of each element of the step response of the system under study, from which the corresponding PWC model of the impulse response is obtained by its time derivative.
Let us consider a finite set of time samples t m with m = 1, . .., M defined in t ∈ [t min , . .., t max ] and the corresponding values of the step response S S,m = S S (t m ).Then, it is possible to identify a subset of N − 1 linear segments, where the extreme points of each segment consist of the pairs (t n , S S,n ), with n = 1, . .., N, for each S S parameter.We call breakpoints the endpoints of the linear segments and SS the corresponding PWL model.So, a PWL model consists of the pairs (t n , SS,n ), with n = 1, . .., N, for a given scattering step response parameter SS .To this aim, it is necessary first to identify the N breakpoints through an iterative approach, summarized in Algorithm 1, until sufficiently accurate results are obtained.
The PWLFIT algorithm starts from two breakpoints: (t min , SS,min ) and (t max , SS,max ).Then, a new breakpoint is chosen in order to minimize the maximum absolute error between the data and the PWL model computed based on the chosen breakpoints.This procedure is repeated until the desired rootmean-square (rms) error is reached.Additionally, it is possible to increase the accuracy of a PWL model without increasing the number of breakpoints as a postprocessing step, by solving a suitable linear system aimed at reducing the error of the PWL approximation in the least square sense [25].This feature of the PWLFIT algorithm is not adopted in this work, so additional information is not provided here.For the interested reader, a complete description of the PWLFIT modeling strategy is given in [25].
As an example, Fig. 3 shows the PWL approximation of a step response using a reduced number of breakpoints.Note that the breakpoints are concentrated on the area where the step response shows a dynamic behavior (< 3 ns) while no breakpoints are chosen where the step response assumes near-constant values, apart from the last breakpoint at 100 ns.
Starting from the PWL model of the step response, a PWC description of a generic scattering parameter s(t) in the TD can Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.be analytically derived as where N − 1 is the number of segments of the PWC model, A i the amplitude of s(t) in the ith interval, D i its center, and ΔT i the temporal length of each constant segment.Moreover, an improvement of the PWLFIT technique, called PWLFIT+, is presented in [36], where it is illustrated how to generate a PWL approximation using the frequency spectra input.Indeed, the Fourier transform of (31) reads A sample of the impulse response spectra of self and mutual S-parameters are illustrated in Fig. 4(a) and (b) for a different number of breakpoints.It can be noticed that increasing the number of breakpoints results in a more accurate spectrum reconstruction, approaching the Fourier transform of the original data.

B. Segment Fast Convolution
As introduced in [25], the computation of the convolution can be significantly accelerated by the SFC while preserving the accuracy of the results.Let us consider the convolution of the single elements in (25), which is repeated for clarity Then, the scattering parameters of the step response are piecewise linearly approximated as described in Section V-A.As it will be more clear in the following, it is convenient to ascribe to the first sample s 1 a single step time window.This implies increasing the number of breakpoints by one.The PWC model of the impulse response of the corresponding S-parameters, shown in (31), can be rewritten as [25] s(τ where It should be noted that (t i , SS,i ) is the ith breakpoint, and W denotes the time interval where time t falls.It is evident that , we obtain Splitting up the W th time window from the previous ones, it can be rewritten as Assuming a time step h for the convolution that, in general, can be different from the sampling rate used to estimate the TD scattering parameters, the time discrete counterpart of (37) is where k is the kth time step, e i is the number of elements within the ith time window and Note that e 0 = 0 and p 0 = k, due to the first single-step time window.The contribution of rect i is omitted for the sake of brevity.Moreover, it is worth noticing that the contribution of the first term of the impulse response impacts the kth time step input contribution, so (38) can be expressed as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where we can identify the direct contribution of the kth time step, that one related to the previous time window with φ i = hs i+1 p i −1 n=p i −e i a(n).Furthermore, the last term of (39) is linked to the W th time window and it is Therefore, in a more compact form, the discrete-time convolution of (33) reads: so, it is evident that the kth sample of the reflected wave b(k) is computed by adding the contribution of the a(k) to the past history.

C. Enhanced SFC
The PWLFIT approach is applied to each parameter s ij with i, j = 1, . ..N p [25], [36], where N p indicates the number of ports of the system under analysis.This guarantees a compact model representation, where the number of breakpoints for each element of the scattering matrix is minimized for the desired accuracy.It is suitable for an efficient implementation based on the parallel computation of the PWL models required.However, this approach will likely result in a different number (and values) of breakpoints t n for each entry of the S-parameters matrix.We indicate the time-value pairs by adding the numbering of the corresponding input/output ports as (t ij−n , s ij−n ).When different breakpoints are defined for each element of the scattering matrix, the convolutions between the S-parameters matrix s and the vector of the incident waves a described in Section V-B must be operated elementwise.
To exploit the potential of the matrix-vector computation of software designed to handle this type of calculation, such as Matlab, it is convenient to adopt a common set of breakpoints, shared by all the parameters s ij , i, j = 1, . ..N p .This can be obtained by considering a sorted vector of all the unique time samples t ij−n , with n = 1, . .., N, i, j = 1, . ..N p .Thus, the vector t b ∈ [t min , . .., t max ] with b = 1, . .., B is obtained, where B indicates the number of samples that constitute the vector of the shared breakpoints.Note that, for reciprocal systems, it is sufficient to consider only the upper (or lower) triangular scattering parameters when computing t b , due to the symmetry of the scattering matrix.
Once the vector t b is obtained, it is sufficient to sample all values of the TD scattering parameters at the values in t b to obtain the new breakpoints: (t ij−b , s ij−b ).This choice entails an increase in the number of breakpoints but allows the acceleration of the convolution.
Indeed, the SFC can now be computed by using the matrixvector product.Similarly to the equation presented in Section V-B, the time discrete counterpart of the convolution in where the W time windows are common for all the Sparameters.Hence, its compact form is

VI. NUMERICAL RESULTS
In this section, some examples are presented to highlight the benefits of the proposed techniques and to validate them by comparing the results with those of well-established techniques.First, the TD characterization of different configurations of the RIS system, performed applying the methodology illustrated in Section III-A, is presented in Section VI-A.Then, the analysis of a different configuration of the RIS with port terminations and sources different from those used to compute the S-parameters is presented in Section VI-B, exploiting the advantage of the preliminary characterization and highlighting the benefits and the validity of the methods described in Sections IV and V.

A. Characterization of the RIS
Similar to the work in [13], a system of dipoles is shown in Fig. 5, and it is constituted by a transmitter, a receiver, and a RIS with a different number of elements.
The main goal is to analyze the overall system using the time S-parameters impulse response described in Section III.

1) First Case of Study:
The distance between the transmitter and the receiver is set to 3 m; the transmitter is located in r T = (5, −1.5, 3) and the receiver is situated in r R = (5, 1.5, 1).In addition, a 4 × 4 dipole array constitutes the RIS, and it is centered in zero.Every dipole is L = λ/32 long, with a square cross-section of side w = λ/250 and a negligible gap between the two branches.The associated wavelength, λ, corresponds to Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.a frequency of 28 GHz.Each element of the RIS is terminated on an RL series impedance Z = R + jωL , with R = 1 Ω and L = 1 nH.The transmitter and receiver are loaded with a 50-Ω resistor.
The system is characterized considering an interdistance between the RIS elements of [ λ 16 , λ 8 , λ 4 , λ 2 ].The choice of subwavelength interdistance among the elements of the RIS offers the possibility of exploiting the mutual coupling between the elements of the RIS from the design, with the possibility to realize super-directive MIMO-antennas [8].The methodology presented in Section III is validated by comparing the spectra of the S-parameters obtained by the fast Fourier transform (FFT) of the results expressed in the TD, denoted in the following with TD, with the S-parameters evaluated directly with a well-established frequency PEEC solver [24], denoted with FD.The S-parameters relative to λ 2 interdistance are illustrated in Fig. 6 and the results show a very good agreement.
Moreover, the virtual LOS channel gain and end-to-end channel gain, introduced in Section II, are illustrated in Figs.7 and 8, respectively.As well, there is an excellent degree of accordance between the data obtained by the FFT of the results in the TD and the results evaluated directly in the FD as exhibited in Figs.7 and 8. Furthermore, similar to what is illustrated in [13], the trend of the virtual LOS channel gain, at the frequency of 28 GHz, for the different values of interdistances is shown in Fig. 9.The aim is to highlight how the interdistance of the elements of the RIS impacts the H VLOS and, therefore, on the H E2E .2) Second Case of Study: The system is analogous to that described in Section VI-A1, but the number of the elements constituting the RIS is equal to 64.For the sake of brevity, in this case, the characterization is carried out by considering only an interdistance between the elements of the RIS of λ 16 .The S-parameters are shown in Fig. 10.Note that the Tx self S-parameter is not reported because it is always the same.
Then, the parameter that best describes how the RIS affects the communication channel: virtual LOS channel gain, is illustrated in Fig. 11.As expected, the results are all in very good agreement.
3) Third Case of Study: In the third case, the configuration is the same as that described in Section VI-A1, but in this case, the distance between the transmitter and the receiver is doubled to 6 m.As for the previous cases, Fig. 12 shows the S-parameters for the case in which the interdistance of the RIS elements is λ 16 .Then, Figs. 13 and 14 illustrate the virtual LOS channel gain and end-to-end channel gain, respectively.
Despite the distance between the transmitter and the receiver is doubled, the correspondence of the FFT of TD results with the FD ones is always very good.
4) Fourth Case of Study: Finally, in the fourth and last case, the system is analogous to that described in Section VI-A1, but the number of the elements constituting the RIS is variable: 4,16,36,64.In this case, the characterization is carried out by considering only an interdistance of λ 16 between the elements of the RIS.The aim is to highlight how the number of the RIS elements affects the virtual LOS channel gain.As expected, Fig. 15 shows that the virtual LOS channel gain increases with the number of the RIS elements.

B. Analysis of the RIS
Let us consider the configuration described in Section VI-A3 that is already characterized.If the power supply of the transmitter is changed and the port terminations are also changed, then the whole system should be studied.However, with the method proposed in Sections IV or V, only (29) should be studied.The transmitter is fed by a modulated trapezoidal pulse with a rise and fall times of τ r = τ f = 1 ns, a width of w = 10 ns, and a modulation carrier frequency of f c = 28 GHz.Furthermore, the RIS elements are terminated on an impedance Z with R = 3 Ω and L = 5 nH.Fig. 16 shows the voltage at the transmitter, the receiver, and one element of the RIS ports.
The results extracted by applying the methods proposed in Sections IV (conv), V-B (SFC), and V-C (ESFC) are validated by the comparison with the nodal analysis (NA), evaluated in the FD and recovering the TD results by the iFFT, and with the time-step second-order integration method (BD2).From Fig. 16, it is clear that the BD2 method becomes unstable; however, the comparison highlights an excellent agreement among the other results.A zoom of the receiver dipole port voltage is shown in Fig. 16(f); in particular, the inset illustrates more clearly that the results of the iFFT-NA are noncausal because they are affected by inaccuracies due to the sampling and frequency spectrum truncation.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Moreover, to better appreciate the main advantages of the proposed convolutive methods, the computational time cost of a single run of the traditional convolution, the SFC e, the enhanced SFC, is reported in Table II which shows that the enhanced segment convolution (ESFC) scheme allows guaranteeing a significant speed-up compared to the standard convolution (conv) and the standard SFC.Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

VII. CONCLUSION
RISs represent one of the most important emerging technologies for next-generation communication systems for their to make the channel more flexible.The outcome of this work is twofold: 1) First, a methodology for accurate mutual coupling aware, end-to-end characterization in the TD of a system composed of a transmitter, a RIS, and a receiver based on the PEEC method was proposed; 2) second, a convolution-based solver that allows including arbitrary linear terminations has been presented.This solver is accelerated by using the SFC technique, resulting in a significant reduction in the computational cost compared to the state of the art.
The proposed method is accurate for the analysis of static linear system scenarios compared to other well-established techniques and, besides, is fully compatible with nonlinear as well as TV terminations.Moreover, the proposed modeling methodology will be the starting point to consider the case of TV scenarios such as those that are present in an urban environment involving moving unintentional scatterers and with the use of TV materials in the realization of the RIS.Future works will also include the analysis of additional acceleration techniques for the TD solver.

Fig. 2 .
Fig. 2. Flowchart for the evaluation of the TD s-parameters.

Fig. 16 .
Fig. 16.Port voltage of the dipoles.(a) Transmitter port voltage.(b) Zoom of the transmitter port voltage.(c) Element of the RIS port voltage.(d) Zoom of the RIS element port voltage.(e) Receiver port voltage.(f) Zoom of the receiver port voltage.

TABLE II COMPUTATIONAL
TIME COST FOR 50 000 TIME SAMPLES