Giorgiani

H Bufferand
email: hugo.bufferand@cea.fr

J Bucalossi

G Ciraolo

G Falchetto

Ph A Gallo

N Ghendrih

P Rivals

H Tamain

Yang

G Giorgiani

F Schwander

M Scotto D'abusco

E Serre

Y Marandet

M Raghunathan

Progress in edge plasma turbulence modelling hierarchy of models from 2D transport application to 3D fluid simulations in realistic tokamak geometry

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Progress in edge plasma turbulence modelling hierarchy of models from 2D transport application to 3D fluid simulations in realistic tokamak geometry being still hard to simulate, especially for ITER size tokamaks. Being able to predict the SOL width or the power load imbalance between inner and outer divertor legs even for todays JET size tokamaks is still an open issue. First principle modelling of edge plasma turbulence is thus today a very active topic in the fusion community driving many dedicated research projects such as the TSVV European boundary code project in Europe.

Despite this effort towards first principle modelling of the turbulent transport, 2D transport codes where the cross-field turbulent transport is emulated by ad-hoc diffusion remain today the most popular tool for experiments analysis or engineering applications. In these codes, the prescription of the diffusion coefficients rely on empirical basis and are usually adjusted to match simulation and experimental at a given location (usually midplane profiles). In this empirical procedure, the nature of the turbulence behind these transport coefficients is most of the time not taken into account. The main drawback of this approach is to lack predictability by missing the basic physics of turbulence mechanisms.

In order to address both the first principle modelling of edge turbulence but also to feed information about turbulence into transport codes, a dedicated effort has been made at IRFM in the last two years to develop the new code SOLEDGE3X presented in Section 2 which aims at encompassing a hierarchy of models from standard 2D transport code to 3D first principle turbulence modelling. The development of this new code comes from the merging of the transport code SOLEDGE2D [START_REF] Bufferand | Numerical modelling for divertor design of west device with a focus on plasmawall interactions[END_REF] and the turbulence code TOKAM3X [START_REF] Tamain | The tokam3x code for the edge turbulence fluid simularions of tokamak plasmas in versatile magnetic geometries[END_REF] previously developed within the French fusion community. Integrating features of the two latter codes, SOLEDGE3X is able to simulate tokamak edge plasma either in 2D or 3D, including: the realistic wall geometry, neutrals since it is coupled to EIRENE and impurities since the code is fully multi-species, based on the state-of-theart Zhdanov collisional closure for multi-component plasmas. First results of turbulent plasma in WEST geometry including the sputtering of the wall have been produced in 2019 and represent a major step for first principle modelling of tokamak plasma in realistic conditions.

However, since full 3D simulations remain expensive in CPU time, it is still profitable to use SOLEDGE3X as a transport code in 2D to run fast interpretative simulations. Section 3 presents an original idea implemented in SOLEDGE3X to improve turbulent transport predictability and go beyond the standard empirical approach is to use a reduced model for turbulence in the same fashion as the "kepsilon" widely used in neutral fluid community. One or two equations are added to the standard mass, momentum and energy balance to describe the evolution of the turbulence intensity k and optionally the evolution of the turbulence dissipation ε. This model is used as a platform to include some ingredients of the turbulence physicssuch as the interchange instability in the framework of 2D transport codes. Of course, this reduced model is not first principle and require an empirical closure. We use the multi-machine scaling law for the SOL width λ q to do so, as a first approach of data assimilation. When SOLEDGE3X is used as a transport code with this reduced "k-epsilon" model for turbulence activated, the number of free parameters is drastically reduced since there is no need to prescribe transport coefficients, the "k-epsilon" model predicting a 2D map of cross-field diffusivities. Since the model is based on interchange instability, one recovers for instance the ballooning of radial transport in the low field side mid-plane. Also, since the model is closed using the multi-machine scaling law, the overall level of transport is automatically adjusted to get a SOL width compatible with this scaling law. This closure with the scaling law is thus powerful to obtain reasonable profiles, however it is also a weakness of the model since in principle, one could not reproduce an experiment where the SOL width does not follow the scaling law. This is the drawback of this kind of semi-empirical models which are not first principle and thus limited by the main assumptions behind their closure. In order to test the applicability of this reduced model for turbulence, a series of L-mode TCV shots have been simulated and the simulation results have been compared with experimental data. Even if the simulation does not recover exactly the SOL width measured in the experiment, the overall agreement in term of peak heat flux, density and temperature on the target is quite remarkable. This "k-epsilon" model has shown to be a promising first step toward integration of turbulence physics inside transport simulations.

To go further, the advantage of the SOLEDGE3X code is to be able to run also first principle turbulent simulations where the turbulence intensity can be directly measured, see Section 4. The comparison between turbulent simulations and the reduced model for turbulence should in the future provide a clear path to improve these reduced models for turbulence. In that perspective, first results of full turbulent simulation in TCV like geometry should help interpreting the results obtained with the reduced "k-epsilon" model and identify missing ingredients.

The SOLEDGE3X code

Developed from the merging of the codes SOLEDGE2D and TOKAM3X developed at IRFM, the new code SOLEDGE3X implements a fluid description of the edge plasma based on the drift-fluid approximation. Braginskii like equations can be solved either in 2D or in 3D. The 2D simulations can be used for two different purposes: Firstly, the code can be used as a transport code (similar to SOLPS-ITER [START_REF] Wiesen | The new solps-iter code package[END_REF] for instance); in this case, a symmetry is assumed (e.g toroidal symmetry for a tokamak or symmetry around the magnetic axis for a linear machine) and the turbulent cross-field transport is emulated using anomalous diffusivities. Secondly, the code can be used as a 2D turbulent code (similar to HESEL [START_REF] Nielsen | Simulation of transition dynamics to high confinement in fusion plasmas[END_REF] for instance) assuming flute approximation [START_REF] Bufferand | Three-dimensional modelling of edge multi-component plasma taking into account realistic wall geometry[END_REF]. The code can also be used in 3D. It can then be used for two purposed: as a 3D transport code, for instance including non-symmetric wall elements like antenna limiters or as a 3D turbulence code or as a 3D turbulent code (similar to GBS [START_REF] Halpern | The gbs code for tokamak scrape-off layer simulations[END_REF], STORM [START_REF] Riva | Three-dimensional plasma edge turbulence simulations of the mega ampere spherical tokamak and comparison with experimental measurements[END_REF] or GRILLIX [START_REF] Stegmeir | Grillix: a 3d turbulence code based on the flux-coordinate independent approach[END_REF]) modelling both the turbulent structure micro-scales and the large-scale plasma flows.

Progress in edge plasma turbulence modelling hierarchy of models from 2D transport application to 3D fluid s

The drift-fluid model implemented in SOLEDGE3X

Despite the different nature of the above-mentioned set-up, the code always solves the same system of equations: a mass balance, a parallel momentum balance and an energy balance for every charged species. In addition, the current balance is solved to compute the electric potential. In its present version, the code is electrostatic and the magnetic field is fixed and taken as an input of the code. The set of equations is summarized as follow: for all charged species:

∂ t n + ∇ • (n v) = S n (1)
∂ t (mnv ) + ∇ • mnv v = -∇ p + ZnE (2)
+ ∇ • ν ∇ v b + mnν ⊥ ∇ ⊥ v + S v ∂ t E + ∇ • (E v) = -∇ • pv b + Znv E (3)
+ ∇ • ν v ∇ v b + mnν ⊥ v ∇ ⊥ v + ∇ • κ ∇ T b + h b + χ ⊥ n ∇ ⊥ T + v R + Q + S E
Where n is the species density; v thespecies velocity (1sst order in the drift ordering) that is decomposed as a parallel velocity v , "E cross B" drift v E = E × B/B 2 where E = -∇φ denotes the electric field, diamagnetic drift v = -∇p × B/ZnB 2 and anomalous diffusive transport v D = -D ⊥ ∇ ⊥ n/n. m denotes the species mass and Z the species charge number. The species temperature T is used to compute species pressure p = nT and species total energy E = 3 2 nT + 1 2 mnv 2 . Parallel viscosity ν , parallel heat conductivity κ , parallel friction force in between species R , parallel heat flux due to inter-species collisions h and energy equipartition term Q take values computed by Zhdanov collisional closure for multi-component plasmas [START_REF] Bufferand | Three-dimensional modelling of edge multi-component plasma taking into account realistic wall geometry[END_REF][START_REF] Zhdanov | Transport Processes in Multicomponent Plasma[END_REF][START_REF] Raghunathan | Generalized collisional fluid theory for multi-component, multi-temperature plasma using the linearized Boltzmann collision operator for scrape-off layer/edge applications[END_REF][START_REF] Bufferand | Implementation of multi-component zhdanov closure in soledge3x (pre-print)[END_REF]. For a simple hydrogenic plasma, Braginskiis closure can also be used. Source terms S n , S ( v ) and S E gather ionisation-recombination-radiation sources as well as external sources such as additional power deposit by heating systems. Finally, cross-field diffusivities D ⊥ , ν ⊥ and χ ⊥ can either take classical value (∼ 10 -2 m 2 s -1 ) for turbulent simulations or "anomalous" value (∼ 1m 2 s -1 ) for transport simulations, the cross-field diffusion then emulating turbulent transport. For electrons, only the energy equation is solved, the electron density being computed assuming quasi-neutrality and electron velocity being computed assuming ambipolarity. The current balance is solved in addition to the above set of equations:

∇ • j b + j + j pola + j dif f = 0 ( 4
)
Where j is the parallel current given by the generalised Ohms law:

j = -σ ∇ φ + ∇ p e n e - R ,Te n e (5)
The diamagnetic current j is computed as j = e,i Zn v whereas the polarisation current is computed similarly (neglecting electron polarisation velocity) j pola = Progress in edge plasma turbulence modelling hierarchy of models from 2D transport application to 3D fluid s i Zn v pola . The polarisation velocity is a second order velocity in the drift ordering that takes a quite complex expression:

v pola = - 1 n ∂ t ω + ∇ • ( v ⊗ ω) (6)
where

ω = m ZB 2 n ∇φ + 1 Z ∇p (7)
One can define the vorticity Ω = ∇ • i Z ω, homogeneous to a charge. The current equation hence takes the form of a vorticity equation

∂ t Ω + ∇ • i Z ∇ • ( v ⊗ ω) = ∇ • j b + j + ζ ∇ ⊥ Ω (8)
The diffusion current being used to regularise the expression j dif f = ζ ∇Ω.

Boundary conditions

SOLEDGE3X uses finite volumes numerical scheme to solve the above cited mass, momentum, energy and current balances. The grid is structured and aligned with magnetic flux surfaces and cannot thus be also aligned with the wall. To deal with complex wall geometries, the grid is extended inside the wall and a mask function determines which points are in the plasma and which points are in the wall. Figure 1 shows an example of boundary condition mask for a JET case. Wall boundary conditions are applied on faces where the transition between points in the plasma and points in the wall. The following boundary conditions inspired by Bohm-Chodura boundary conditions apply there:

• Outgoing velocity normal to the wall greater that parallel sound speed normal to the wall:

| v • n wall | ≥ |c s b • n wall | (9)
This property guarantees that the total plasma velocity is oriented outward

• Sheath transmission factor between energy and particle fluxes:

φ E,se = γT + 1 2 mv 2 φ n,se (10)
Where for each species, φ E,se is the total energy flux at sheath entrance, φ n,se is the particle flux at sheath entrance and γ is the sheath transmission factor about 2.5 for ions and 4.5 for electrons.

• Total plasma current on the wall is given by

j wall = 1 -exp Λ - φ T e i Zφ n,se (11)
Where the ions saturation current is computed from ions particle fluxes i Zφ n,se and where Λ denotes the normalised potential drop in the sheath Λ ∼ 3.

Numerical implementation

SOLEDGE3X relies on a finite volumes numerical scheme. Divergences in every grid point are computed integrating fluxes on the faces of grid elements. The grid is structured and aligned with magnetic surfaces. A domain decomposition is used to treat X-point topology, see Figure 1 for a JET case where 6 domains are used to deal with the X-point. The following numerical methods are used to compute fluxes:

• Advection: 3rd order WENO interpolation and Donat, Marquina fluxes for modified Riemann solver [START_REF] Donat | Capturing shock reflections: An improved flux formula[END_REF] • Parallel diffusion: since the grid is aligned with magnetic flux surfaces, the parallel diffusion can be solved independently for each flux surface. For 2D cases, the parallel diffusion can be treated as N F S (number of flux surfaces) 1D problems aligned with the grid: simple 2nd order finite differences are used. For 3D cases, the parallel diffusion becomes N F S 2D problems. The grid is no longer aligned with parallel diffusion which occur across both the poloidal and toroidal direction.

Günter scheme [START_REF] Günter | Modelling of heat transport in magnetised plasmas using non-aligned coordinates[END_REF] is used to avoid numerical diffusion in this case.

• Perpendicular diffusion: perpendicular diffusion couples all the points of the domain. We use preserving monoticity scheme for anisotropic diffusion [START_REF] Sharma | Preserving monotonicity in anisotropic diffusion[END_REF].

The implementation of the different operators has been verified using the method of manufactured solution [START_REF] Bufferand | Implementation of multi-component zhdanov closure in soledge3x (pre-print)[END_REF].

Most of operators are treated explicitly in time. Only parallel diffusion which has the fastest dynamic is treated implicitly for mass, momentum and energy balances. The current equation takes a different form as other equation. The electric potential being the unknown, it can be rewritten schematically ∂ t ∆ ⊥ φ + ∆ φ = RHS where the time derivative applies to the perpendicular Laplacian of the unknown. This perpendicular Laplacian must then also be solved implicitly. The matrix used to write the electric potential equation hence couples all the points of the domain. Besides, the anisotropy between a fast dynamic in the parallel direction and a slower one in the perpendicular direction makes the matrix ill-conditioned. If direct solvers such as PASTIX [START_REF]Gitlab of the pastix library[END_REF] can be used for small 2D cases, they become too time consuming for 3D cases. The best approach found so far relies on iterative solvers, more precisely algebraic multigrid preconditioner followed by Krylov solvers (GMRES or BiCGstab). Different libraries have been coupled to the code such as AGMG [START_REF]Agmg webpage[END_REF], PETSC-GAMG [START_REF]Portable, extensible toolkit for scientific computation website[END_REF], or HYPRE-BoomerAMG [START_REF]Scalable linear solvers and multigrid methods website[END_REF].

Neutrals

The code SOLEDGE3X implements two neutral models:

• A crude fluid model for neutrals where neutrals transport is considered diffusive.

Despite its simplicity, this model provides a good approximation for the source of plasma generated by recycling as well as power losses by radiation in the divertor. It is used mainly for turbulent simulations where the detail of neutral transport

is not yet necessary and where a gross estimation of the sources is enough for the scoping studies

• A full kinetic description of neutrals including complex neutrals chemistry. This is done by a coupling of SOLEDGE3X with the neutral Monte-Carlo code EIRENE [START_REF]Eirene website[END_REF] widely used in the community. A grid based on triangles is used for neutrals, see Fig 1. Interpolations between the SOLEDGE3X quadrangles and the EIRENE triangles are done within the SOLEDGE3X-EIRENE interface, preserving mass and energy conservation.

k-epsilon model for turbulence prediction in transport mode

For experiments interpretation and for the preparation of future tokamaks operation, two-dimensional transport simulations remain very popular for integrating a broad range of physical aspects of the plasma wall interaction, from plasma recycling, wall sputtering to impurity transport. They are also quite fast to compute and make possible parameter scans for a reasonable computing time. However, the main drawback is the poor description of the cross-field transport since the plasma turbulent structure cannot be simulated assuming toroidal symmetry. In the mean field approach, the turbulent transport is emulated by an effective diffusion process of mass, momentum and energy, the ad-hoc diffusivities being to be determined most of the time empirically. In the diffusive assumption, the transport becomes local since the flux is determined only by the local gradient. Turbulence measurements and simulations show that turbulent structures can however propagate ballistically over a long range which may contradict the local approximation made in the diffusive assumption. Several contributions have proposed non-local formulations [START_REF] Dif-Pradalier | On the validity of the local diffusive paradigm in turbulent plasma transport[END_REF] to go beyond diffusion and the subject is still open. Another key missing ingredient in the mean field approach is by definition the fluctuations. In the averaging procedure of transport equations, several non-linearities are neglected, see the example below for ionization cross-section for instance:

S n = n n n e σv iz (T e ) ( 12
)
Averaging gives:

S n = n n n e σv iz (T e ) + n n n e σv iz (T e ) (13)
Most of the time, in mean field codes, the second term is neglected and the non-linearity of the ionisation cross-section in neglected giving:

S n = n n n e σv iz (T e ) (14)
The impact of the excursion of plasma parameters around their average due to fluctuation can be taken into account to compute "fluctuation dressed" cross-section for instance, see [START_REF] Guzman | Ionization balance of impurities in turbulent scrape-off layer plasmas i: local ionization-recombination equilibrium[END_REF][START_REF] Marandet | Assessment of the effects of scrape-off layer fluctuations on first wall sputtering with the tokam-2d turbulence code[END_REF]. However, to use these expressions in practice in a transport code, one needs at least to estimate the level of fluctuations in the plasma.

In order to feed some information about turbulence in mean field transport simulations, a heuristic turbulence model inspired from the "k-epsilon" model used in neutral fluid community has been implemented in SOLEDGE3X. It relies either on a single equation describing the production, saturation and transport of the turbulence intensity k ∼ ṽE 2 [m 2 s -2 ] or a coupled system of two equations describing the evolution of the turbulence intensity k and of its dissipation rate ε. The second model enables in principle a richer description of turbulence properties, see [START_REF] Baschetti | A k-epsilon model for plasma anomalous transport in tokamaks: closure via the scaling of the global confinement[END_REF][START_REF] Ph | Self-consistent cross-field transport model for core and edge plasma transport[END_REF]. Other models have been proposed in the community with one or two equations, see [START_REF] Wolf | Bayesian approach to parameter estimation and model validation for nuclear fusion reactor mean-field edge turbulence modelling[END_REF][START_REF] Coosemans | Bayesian analysis of turbulent transport coefficients in 2d interchange dominated exb turbulence involving flow shear[END_REF]. The field is still very open to find an efficient model for turbulence enabling a better predictability for cross-field transport in transport codes, going beyond the state-of-the-art empirical setting of cross field diffusivities. We present in this contribution the simplest model, the single k equation implemented in SOLEDGE3X. Despite its simplicity, this model computes self-consistently a map of diffusion coefficients. The latter is computed from k as D = τ k where the time τ is taken as τ = R/c s . In addition to the diffusion coefficient, the model also provides an estimation of the fluctuation level in the plasma enabling the use of "fluctuation dressed" cross sections (though not yet tested in SOLEDGE3X so far). The equation for k is the following:

∂ t k + ∇ • (k v) = S k -P k + ∇ • D ∇ ⊥ k (15)
Where S k denotes the source of turbulence and P k a saturation mechanism for turbulence intensity. Interchange instability has been considered as the main source of turbulence Progress in edge plasma turbulence modelling hierarchy of models from 2D transport application to 3D fluid s and thus the linear growth rate of interchange instability is used to compute S k :

S k = γ I k = c s max ∇p p • ∇B B , 0 k (16)
The saturation of turbulence remains a free parameter in the model. We chose it to be quadratic hence P k = αk The model has been applied with success to TCV L-mode discharges [START_REF] Baschetti | A k-epsilon model for plasma anomalous transport in tokamaks: closure via the scaling of the global confinement[END_REF]. Below is presented an example of application to a JET detached plasma from the M18-27 "Isotope effects on detachment in L-mode" experiment. The experimental set-up is summarized in Figure 2. The magnetic equilibrium used corresponds to JET shot #95235. This shot is a density ramp in L-mode with a fixed 1MW NBI heating during the ramp and an Ohmic power in the range of 2MW. Figure 5 compares radiation location between simulation and experimental tomographic reconstruction. The radiation seems more concentrated around the Xpoint in the simulation compared to the experimental radiation map. This may be due to the absence of radiating impurities in the simulation. Also, the density peaking in the closed flux surface region is not caught by the simulation. A radial pinch might be needed to recover this feature.

Full 3D turbulent simulations in X-point geometry

In addition to its ability to emulate turbulence in a 2D transport mode, either with adhoc diffusivities or with reduced models for turbulence, the SOLEDGE3X code is able to perform full 3D turbulent simulations including the main ingredients of the plasma-wall interactions such as recycling or wall sputtering. Here is an illustration of 3D simulation of the WEST tokamak, including the complex wall geometry of WEST. The diffusive neutral model is used to simulate plasma recycling. The size of the machine and the intensity of the magnetic field have been reduced for numerical reasons, the resolution needed for a real scale WEST simulation being out of reach of the code today. Figure 6 shows simulation results, in particular the nice plasma filaments propagating on the tokamak low-field side.

Wall sputtering by the plasma filaments has also been taken into account once again using the simple diffusive model for the neutral sputtered species (in this case, despite the wall is in Tungsten, we chose to sputter Carbon as a proxy for the light impurities trapped at the surface of the tungsten such as oxides, Carbon co-deposits...). One notices that despite its simplicity, the fluid neutral model generates plasma sources at the right location and with a reasonable intensity. It is thus a good candidate to start implementing neutral physics in edge plasma turbulent simulations and it has been implemented in other codes such as GRILLIX [START_REF] Zholobenko | Simulations of turbulence, its suppression and profile evolution across the edge and scrape-off layer of asdex upgrade tokamak[END_REF].

Conclusions

The new SOLEDGE3X code has been developed in the past few years for a versatile use both as a transport code and as a turbulence code including complex wall geometry, neutrals and impurities. Reduced model for an enhanced prediction of turbulent transport have been implemented, based on "k-epsilon" like models from the neutral fluid community. First results applied to JET and to WEST have been shown in this contribution. A dedicated effort will be made in the near future to improve numerical efficiency of the code, in particular to speed-up the inversion of the electric potential equation which remains the bottleneck in term of computation workload.

Figure 1 .

 1 Figure 1. Left: Example of domain decomposition for a JET case (single X-point configuration). Right: EIRENE grid used for neutrals. In black, triangles where plasma is computed -in red, triangles without plasma background (subdivertor).

2 . 2 SOL = 2πqR 2 c 2 sc 2 s 1 λ 2 scaling.

 2222 We decided to close it from experiments, ensuring that the model recovers experimental scaling laws for the SOL width. More precisely, neglecting transport terms in the k equation, at steady state source and saturation terms must compensate each other giving S k = P k hence k = γ I /α. Given the expression D = Rk/c s , and the theoretical SOL width for a diffusive model with Bohm boundary conditions at the target:λ 2 SOL = DL /c 2 s , one gets λ × γ I α .If one wants the system to recover λ SOL = λ scaling , the parameter α takes the value α = γ I 2πqR 2 For L-mode plasma, we usually consider the simple scaling law λ scaling = 4qρ L , the parameter α can then be rewritten as α = γ I

Figure 2 .

 2 Figure 2. Left: magnetic equilibrium for JET #95235 at t = 54s. Gas injection locations are shown. Right: Time traces of core density, power and gas injection during the discharge

Figure 3 .

 3 Figure 3. Comparison between midplane profiles computed by SOLEDGE3X-EIRENE and measured on JET

Figure 4 .

 4 Figure 4. Poloidal map of turbulence intensity k (left) and subsequent cross-field diffusivity D (right) computed by the reduced model

Figure 3

 3 Figure3shows a comparison between simulation results and experimental measurments. The overall gradients are well recovered by the k-epsilon model. Figure4shows quantities computed by the reduced turbulence model. One recovers a ballooning of the turbulence intensity on the low field side. Turbulence propagates on the high field side by parallel transport.Figure5compares radiation location between simulation and experimental tomographic reconstruction. The radiation seems more concentrated around the Xpoint in the simulation compared to the experimental radiation map. This may be due

  Figure3shows a comparison between simulation results and experimental measurments. The overall gradients are well recovered by the k-epsilon model. Figure4shows quantities computed by the reduced turbulence model. One recovers a ballooning of the turbulence intensity on the low field side. Turbulence propagates on the high field side by parallel transport.Figure5compares radiation location between simulation and experimental tomographic reconstruction. The radiation seems more concentrated around the Xpoint in the simulation compared to the experimental radiation map. This may be due

Figure 5 .

 5 Figure 5. Comparison between radiation map computed by SOLEDGE3X-EIRENE and tomographic reconstruction from bolometry. The strong density gradients near the separatrix are explained by the presence of a strong particle source around the X-point, confirmed by X-point radiation

Figure 6 .

 6 Figure 6. SOLEDGE3X simulation results for a 3D WEST case including recycling and Carbon sputtering. Top left: electron temperature. Top right: electron density. Bottom: Carbon neutral density.

Acknowledgements

This work was granted access to the HPC resources of CINES, under the allocation 2020-A0080510482 made by GENCI and to the HPC resources of Aix-Marseille University financed by the project Equip@Meso (ANR-10-EQPX-29-01). This work Progress in edge plasma turbulence modelling hierarchy of models from 2D transport application to 3D fluid s was supported by the EoCoE-II project, grant agreement 824158, funded within the EUs H2020 program. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission

Bibliography


Reconfigurable Intelligent Surface-Based Quadrature Reflection Modulation Shaoe Lin * , Miaowen Wen * , Marco Di Renzo † , and Fangjiong Chen * * School of Electronic and Information Engineering, South China University of Technology, Guangzhou 510641, China † Université Paris-Saclay, CNRS and CentraleSupélec, Laboratoire des Signaux et Systèmes, Gif-sur-Yvette, France Email: eemwwen@scut.edu.cn

Abstract-Reflection modulation based on reconfigurable intelligent surface (RIS) is considered to be a promising information transfer mechanism that conveys information by reconfiguring incident electromagnetic waves without using radio frequency chains. Existing reflection modulation schemes propose to convey additional information by manipulating the ON/OFF states of reflecting elements, which underuse the potential of the RIS. In this paper, we first propose a new reflection modulation scheme, called RIS-based quadrature reflection modulation (RIS-QRM), to achieve joint passive beamforming and information transfer for RIS-aided single-input single-output wireless communications by assuming that the direct link between the transceivers is blocked. We then extend the RIS-QRM scheme to the multiple-input single-output setting with the direct link between transceivers being taken into account. In this case, an optimization problem is formulated to jointly optimize the active beamforming and the phase shifts of the RIS in order to improve the received signal power at the user. Computer simulation results corroborate the effectiveness of the RIS-QRM scheme. It is shown that in contrast to ON/OFF-based schemes, the RIS-QRM scheme is able to improve the error performance of the additional bits delivered by the RIS without deteriorating that of the bits carried by the constellation symbols.

I. INTRODUCTION

Reconfigurable intelligent surfaces (RISs) have been envisioned as a promising technology for beyond fifth-generation networks [START_REF] Di Renzo | Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and road ahead[END_REF]. Specifically, an RIS is a planar surface composed of numerous passive elements, each being a reconfigurable scatterer that reflects the incoming signal with a adjustable phase shift. By properly adjusting the phase shifts, the reflected signals can be added constructively with the signal from the direct path at the receiver to improve the signal-to-noise ratio (SNR) or destructively to mitigate interference/eavesdropping. In the literature, this is referred to as passive beamforming. The advantages of RISs, such as enhancing the communication performance and reducing the energy consumption, have motivated a variety of RIS-aided wireless communication systems, e.g., [START_REF] Zheng | Intelligent reflecting surface-assisted multiple access with user pairing: NOMA or OMA?[END_REF]- [START_REF]Double-IRS assisted multi-user MIMO: Cooperative passive beamforming design[END_REF].

Spatial modulation (SM), which activates only one transmit antenna per time slot to convey one constellation symbol and embeds additional information in the index of the active antenna, has been receiving significant attention as it is capable of striking an attractive tradeoff between the energy efficiency and the spectral efficiency [START_REF] Di Renzo | Spatial modulation for generalized MIMO: Challenges, opportunities, and implementation[END_REF], [START_REF] Wen | A survey on spatial modulation in emerging wireless systems: Research progresses and applications[END_REF]. A recent generalization is to implement SM by an RIS [START_REF] Li | Single-RF MIMO: From spatial modulation to metasurface-based modulation[END_REF]. In [START_REF] Basar | Reconfigurable intelligent surface-based index modulation: A new beyond MIMO paradigm for 6G[END_REF], the author proposed to carry out SM on receive antenna indices by an RIS that steers a beam towards a particular receive antenna according to the input information bits. In [START_REF] Yan | Passive beamforming and information transfer design for reconfigurable intelligent surfaces aided multiuser MIMO systems[END_REF], the authors proposed a passive beamforming and information transfer (PBIT) scheme in which an RIS is used to improve communication performance in the uplink direction and simultaneously deliver additional information via the ON/OFF states of its elements. In order to minimize the communication outage probability, the authors of [START_REF] Lin | Reconfigurable intelligent surfaces with reflection pattern modulation: Beamforming design and performance analysis[END_REF] proposed an RIS-based reflection pattern modulation (RIS-RPM) scheme to achieve PBIT, where the number of RIS elements that are switched ON is fixed at all time slots.

Existing reflection modulation schemes [START_REF] Yan | Passive beamforming and information transfer design for reconfigurable intelligent surfaces aided multiuser MIMO systems[END_REF]- [START_REF] Guo | Reflecting modulation[END_REF] use the ON/OFF states of the RIS elements to carry additional information and therefore have two main drawbacks. Firstly, since only a portion of the available RIS elements are switched ON, the potential performance improvement offered by an RIS is underused, resulting in relatively low received signal power levels. Secondly, due to the inherent randomness of the RIS configuration that depends on the incoming random bits, the joint optimization of the beamforming at the transceivers and the RIS configuration is a stochastic optimization problem that is difficult to solve in general. To address the above issues, we consider in this paper an RIS-aided downlink multipleinput single-output (MISO) wireless system, and propose a new reflection modulation scheme, referred to as RIS-aided quadrature reflection modulation (RIS-QRM), to achieve PBIT by considering the full reflection capabilities of the RIS. We first consider the single-input single-output (SISO) setting and present the principle of the RIS-QRM scheme. Specifically, all the reflecting elements at the RIS are randomly partitioned into two subsets according to the RIS local data: (i) the phase shifts of the elements in the first subset are tuned to -θ and (ii) the phase shifts of the elements in the second subset are tuned to -θ + π/2, where θ is the phase of the cascaded channel associated with the corresponding reflecting element. In this manner, the RIS element partition carries the local data of the RIS. We then extend the RIS-QRM principle to the general MISO setting by taking the direct link between the transceivers into account. In this case, we formulate an optimization problem to jointly optimize the active beamforming at the access point (AP) and the phase shifts of the RIS to improve the received signal power at the user. Extensive computer simulation results are presented to corroborate the effectiveness of the proposed RIS-QRM scheme. In particular, compared to the RIS-RPM scheme, the RIS-QRM scheme is able to improve the error performance of the additional bits delivered by the RIS without deteriorating that of the bits carried on the constellation symbols.

Notation: Superscripts † , T , and H represent complex conjugate, transpose, and Hermitian transpose operations, respectively. 1 n×m denotes an all-one matrix of size n × m. | • | denotes the absolute value. For a complex-valued vector x, ∠x denotes a column vector with the elements being the phases of the elements in x, and diag(x) denotes a diagonal matrix with the elements in x on its main diagonal. • returns the largest integer less than or equal to the argument. • • denotes the binomial coefficient. A\B denotes the relative complement of set B in set A. E[•] denotes the expectation.

II. SYSTEM MODEL AND OVERVIEW OF EXISTING REFLECTION SCHEMES

A. System Model

As illustrated in Fig. 1, we consider a downlink RISaided MISO wireless system, where an RIS consisting of L 0 reflecting elements is deployed to assist the transmission from an AP equipped with N transmit antennas to a singleantenna user located in its proximity. To reduce the channel estimation overhead and implementation complexity, we recall the RIS-elements grouping method introduced in [START_REF] Zheng | Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization[END_REF], where the L 0 RIS elements are divided into L groups, denoted by L {1, 2, . . . , L}, each consisting of L 1 = L 0 /L adjacent elements that have the same reflection coefficient. Moreover, the RIS is connected to a smart controller, which is responsible for exchanging information between the AP and RIS via a feedback link and for the dynamic adjustment on the reflection of each RIS element. We consider that the RIS is inherently equipped with sensors for monitoring/collecting environmental data (typically low-rate bursty data). Therefore, the RIS needs to report sensed data to the user.

We consider a quasi-static block fading channel model for all individual links, where the channel coefficients remain constant within each fading block and are independent and identically distributed across different fading blocks. Without loss of generality, we focus on the downlink communication in one particular fading block. Let

h H d ∈ C 1×N denote the baseband channel of the AP→user direct link. Let G ∈ C L0×N G H 1 , G H 2 , . . . , G H L H and h H r ∈ C 1×L0 h H r,1 , h H r,2 , . . . , h H r,
= √ -1, β l ∈ [0, 1] and θ l ∈ [0, 2π
) denote the common reflection amplitude and phase shift of the RIS elements in group l, respectively. For maximizing the reflection power of the RIS and simplifying its hardware design, the reflection amplitude of each RIS element is set to be the maximum value, i.e., β l = 1, ∀l = 1, 2, . . . , L. Accordingly, the baseband received signal at the user is given by

y = L l=1 h H r,l e jθ l G l + h H d wx + n (1)
where w ∈ C N ×1 is the digital beamforming vector at the AP, x is the transmitted signal with power P , and n ∼ CN (0, N 0 ) is the additive white Gaussian noise (AWGN) with noise power N 0 . The transmitted signal x is drawn from an M -ary signal constellation, such as the pulse amplitude modulation (PAM), denoted by the set M of size M . Let h H l = h H r,l G l ∈ C 1×N with l = 1, 2, . . . , L denote the effective (aggregated) channel for the AP→group l→user cascaded link in the absence of the corresponding reflection coefficient. Therefore, we have

y = L l=1 e jθ l h H l + h H d wx + n = θ H H + h H d wx + n ( 2
)
where

H = [h 1 , h 2 , . . . , h L ]
H ∈ C L×N denotes the AP-RIS-user cascaded channel matrix in the absence of reflection coefficients. Accordingly, the received SNR is given by

SNR = P N 0 θ H H + h H d w 2 . ( 3
)
For simplicity, we assume in this paper that the perfect channel state information (CSI) of H and h H d is available at the AP.

B. Overview of Existing Reflection Schemes

A proper RIS configuration based on CSI has a great impact on the communication performance between the AP and the user. To be specific, when the phase shifts of the L groups are properly adjusted, the reflected signals by the groups are added constructively (with the signal coming directly from the AP) at the user, thus achieving the maximum SNR. This is referred to as passive beamforming. An RIS, however, despite its functionality of enabling the maximum received SNR, is unable to convey any information of its own. To deliver additional information, the PBIT scheme adjusts the ON/OFF state of each group according to the incoming random bits, while the RIS-RPM scheme intentionally switches off a fixed number of groups. Both schemes perform passive beamforming based on a portion of the available RIS elements, thus inevitably suffering from a received power/SNR loss.

III. PRINCIPLE OF RIS-QRM

To give an intuitive idea of the operating principle of RIS-QRM, in this section, we assume that the AP is equipped with only one transmit antenna, i.e., N = 1 with G ≡

g = g H 1 , g H 2 , . . . , g H L H
and

H ≡ h = [h 1 , h 2 , . . . , h L ] H .
Moreover, we assume that the AP→user direct link is blocked due to obstructions, i.e., h H d ≡ h † d = 0. We extend the RIS-QRM design to the MISO setting with the AP→user direct link being taken into account in the next section. The baseband received signal in ( 2) is now expressed as

y = θ H hx + n (4)
and thus the preferred phase shifts 1 for achieving coherent signal combination at the user are given by

θ * [θ * 1 , θ * 2 , . . . , θ * L ] = ∠h. (5)
For implicitly conveying the RIS information, the L groups are randomly partitioned into two subsets: the phase shifts of the first subset with p groups are tuned to their corresponding preferred phase shifts such that the reflected signals via these p groups are in phase, while the phase shifts of the second subset with the remaining L -p groups are tuned to their corresponding preferred phase shifts plus π/2 so that the reflected signals via the remaining L -p groups are in quadrature with those reflected via the first subset. By design, given L and p, there are L p group-partition realizations. This implies that at most b = log 2

φ l = θ * l , l∈ I θ * l + π/2, l ∈ L \ I . ( 6
)
1 The AP feedbacks the preferred phase shifts to the RIS which adjusts the reflections of its elements accordingly.

Accordingly, the RIS configuration that carries its own information bits can be expressed as φ e jφ1 , e jφ2 , . . . , e jφ L H = diag(θ * )s

where s [s 1 , s 2 , . . . , s L ] H carries the RIS information with each entry given by

s l = 1, l ∈ I j, l ∈ L \ I , l = 1, 2, . . . , L. (8)
Let S {s 1 , s 2 , . . . , s B } denote the set of all possible realizations of s. With the information embedded in the RIS configuration, the baseband received signal at the user is given by

y = φ H hx + n = s H diag(θ * ) H hx + n. ( 9
) Define v [v 1 , v 2 , . . . , v L ] H = diag(θ * ) H h,
where v l = |h l | denotes the effective cascaded channel from the AP to the user via the l-th group with the preferred phase shift, l = 1, 2, . . . , L. Then, we have

y = s H vx + n = ⎛ ⎝ l∈I v l + j l∈L\I v l ⎞ ⎠ x + n. ( 10
)
Suppose that the CSI of v is perfectly known at the user, which can be obtained by downlink pilot training. Then, the maximum-likelihood (ML) estimate of < x, s > is given by

< x, ŝ >= arg min x∈M, s∈S y -s H vx 2 . ( 11
)
It is worth pointing out that although the proposed RIS-QRM scheme has the same spectral efficiency as the RIS-RPM scheme (i.e., b + log 2 M ), the former improves the received power by reactivating the inactive elements of the latter, thus increasing the achievable rate in the low SNR regime.

IV. EXTENSION TO MISO SETTING

We extend the RIS-QRM design for the SISO setting to the MISO case with the AP→user direct link being taken into account. In this more general case, the active beamforming at the AP needs to be appropriately designed to improve the communication performance together with the passive beamforming at the RIS. Specifically, we jointly optimize the active beamforming at the AP and the phase shifts of the RIS based on the CSI of H and h H d to improve the received signal power at the user. Note that this is effective in improving the error performance of the bits carried on x as well as s, as it will be shown in the simulation results.

1) Joint Active and Passive Beamforming: Our objective is to maximize the received signal power at the user by jointly optimizing the active beamforming at the AP and phase shifts of the RIS, subject to the power budget of the AP and the unit-modulus reflection constraint of the RIS. The optimization problem is formulated as follows.

(P1): max

w,θ θ H H + h H d w 2 (12)
s.t.

w 2 ≤ 1 (13) 0 θ l < 2π, l = 1, 2, . . . , L. ( 14
)
It is easily verified that this problem is always feasible, since there exists at least one solution that satisfies the constraints in ( 13) and ( 14), e.g., θ = 1 L×1 and w = (11×LH+h H d )

H 1 1×L H+h H d .
However, it is non-convex since the objective function in [START_REF] Yan | Passive beamforming and information transfer design for reconfigurable intelligent surfaces aided multiuser MIMO systems[END_REF] is not jointly concave with respect to w and θ, and the unitmodulus constraint in ( 14) is non-convex. This problem can be suboptimally solved by leveraging the alternating optimization (AO) technique. Specifically, based on the triangle inequality, we have

θ H H + h H d w ≤ θ H Hw + h H d w (15)
which holds with equality if and only if ∠θ H Hw = ∠h H d w. Note that since the received signal power increases monotonically with the transmit power (i.e., P w

2 ), the inequality constraint in [START_REF] Lin | Reconfigurable intelligent surfaces with reflection pattern modulation: Beamforming design and performance analysis[END_REF] needs to be met with equality at the optimal solution, i.e., w * 2 = 1. Therefore, for a fixed active beamforming direction w, problem (P1) is equivalent to the following optimization problem.

(P2):

max θ θ H Hw 2 (16) s.t. 0 θ l < 2π, l = 1, 2, . . . , L (17)
∠θ H Hw = ∠h H d w. (18)
The optimal solution to problem (P2) is given by θ = e j(∠Hw-∠h H d w) , where the preferred phase shift of the RIS elements in the l-th group is given by

θ l = ∠h H d w -∠h H l w, l = 1, 2, . . . , L. ( 19
)
Then, for a given θ in (19), the optimal w is the wellknown maximum-ratio transmission (MRT) solution, i.e., w = (θ

H H+h H d ) H θ H H+h H d
. θ and w are alternately optimized until the fractional increase of ( 12) is less than a sufficiently small value > 0 or the maximum number of iterations, denoted by ζ, is reached. The iterative process is summarized in Algorithm 1. We note that this iterative process is guaranteed to converge since the objective value in ( 12) is non-decreasing over the iterations and the optimal objective value is upper-bounded.

2) Reflection Modulation Design:

We assume to use the optimal active beamforming and phase shifts obtained from problem (P1), denoted by w * and θ * . Note that the preferred phase shifts in (19) make the signals from the AP to the user via the RIS in phase with the signal over the AP→user direct link so as to achieve coherent signal combination at the user. For implicitly conveying the RIS information, the L groups are randomly partitioned into two subsets: the phase shifts of the first subset with p groups are tuned to their corresponding preferred phase shifts, while the phase shifts of the second subset with the remaining L -p groups are tuned to their corresponding preferred phase shifts plus π/2 so that the signals from the AP to the user via the remaining L -p groups are orthogonal to the signal over the AP→user Algorithm 1 AO Algorithm for Solving Problem (P1)

Input: , ζ 1: Initialize the active beamforming w (1) := h d h d and set the iteration number q := 1 2: repeat 3:

For given w (q-1) , obtain the phase shifts θ (q) according to (

For given θ (q-1) , perform MRT at the AP based on the combined AP-user channel to obtain the active beamforming w (q) 5:

Update q := q + 1 6: until The fractional increase of ( 12) is less than or q > ζ Output: w * and θ * direct link. The baseband received signal at the user is given by

y = φ H H + h H d w * x + n = s H diag(θ * ) H Hw * + h H d w * x + n. (20)

V. SIMULATION RESULTS AND DISCUSSION

We present computer simulation results to numerically evaluate the performance of the RIS-QRM scheme in terms of bit error rate (BER). The total number of the RIS elements is L 0 = 160. The channel path loss is modeled as γ = γ 0 D -α , where γ 0 is the path loss at the reference distance of 1 meter (m), D and α are the individual link distance and path loss exponent, respectively. Let D 1 , D 2 , and D 3 denote the distances of the AP→RIS, RIS→user, and AP→user links, respectively. For small-scale fading, we consider the Rician fading channel model for all individual links. The path loss exponents of these three links are denoted by α 1 , α 2 , and α 3 , respectively, and the Rician factors of these links are denoted by κ 1 , κ 2 , and κ 3 , respectively. We set D 1 = D 3 = 50 m and D 2 = 5 m. The reference path loss is γ 0 = 30 dB for all individual links. The path loss exponents of the AP→RIS and RIS→user links are α 1 = 2.2 and α 2 = 2.6, respectively. The Rician factors of the AP→RIS and RIS→user links are κ 1 = 1 and κ 2 = 1, respectively. Let γ 3 denote the average channel power of the AP→user direct link, which is given by γ 3 = γ 0 D -α3 3 unless otherwise specified. We set α 3 = 3.3 and κ 3 = 0. Other parameters are as follows: M = 256, N 0 = -80 dBm, ζ = 10 and = 10 -4 . The simulation results are averaged over 5, 000 independent fading channel realizations.

In Fig. 2, we show the error performance of the RIS-QRM scheme versus the SNR assuming N = 1, L = 2, p = 1, γ 3 = 0 (i.e., the AP→user direct link is blocked), and 256-PAM. Two existing ON/OFF-based reflection modulation schemes are considered for comparison: i) the RIS-RPM scheme in which the elements in the vector s indexed by L \ I take the value of 0, i.e., s l = 0, ∀l ∈ L\I; ii) the PBIT scheme in which the elements in the vector s are independently drawn from the set {0, 1} with equiprobability and s remains fixed over L/b consecutive time slots. The lower bounds (LBs) on the error performance of x and s are obtained by assuming that the knowledge of s and x are perfectly known when detecting x and s. We observe that the proposed scheme not only significantly improves the error performance of x, but also improves that of s, as compared to the RIS-RPM scheme. This is attributed to the use of the quadrature component of the combined AP-user channel, which enhances the received SNR and increases the minimum Euclidean distance between all possible realizations of s by adding additional degrees of freedom.

We disclose the impact of p on the error performance of the RIS-QRM and RIS-RPM schemes in Fig. 3, where p varies from 1 to 3 and other simulation parameters are N = 4, L = 4, p = {1, 2, 3}, and 256-PAM. For clarity, the lower bounds are omitted. We observe that the size of the first subset p has a minor impact on the performance of the proposed scheme, which is due to the reciprocity between the in-phase and quadrature components of the combined AP-user channel. In contrast, the size of the first subset has a significant impact on the error performance of the RIS-RPM scheme, which can be understood by the fact that the reflection power in this scheme is determined by the number of groups in the first subset. In particular, the proposed scheme is always better than the RIS-RPM scheme, regardless of the value of p. Compared with the PBIT scheme, the proposed scheme significantly improves the BER performance of x, which is attributed to the full reflection of the RIS. It is worth noting that the PBIT scheme gives rise to communication outage, which causes an irreducible error floor in the BER performance of x, as shown in Fig. 2(a). In the BER performance of s, on the other hand, we observe from Figs. 2(b) and 3(b) that the proposed scheme outperforms the PBIT scheme in the low-to-medium SNR region, while the PBIT scheme shows a slightly superior performance to the proposed scheme at high SNR. This can be understood by the fact that, in order to fairly compare those three schemes with the same information rate, the PBIT scheme adopts repetition coding.

In Fig. 4, we compare the error performance of the proposed RIS-QRM scheme by using different active beamforming designs versus the SNR, where N = 4, L = 4, p = 3 and 256-PAM. The following active beamforming schemes are considered for comparison: 1) the AP-user direct path MRT scheme where the active beamforming at the AP is steered toward the AP-user direct path, i.e., w = h d h d ; 2) the strongest path MRT scheme where the active beamforming at the AP is steered toward the path of the largest channel power gain among the cascaded channels of all groups and the AP-user direct channel, i.e., w = h max h max with h max = arg max{ h 1 2 , . . . , h L 2 , h d 2 }; 3) the weakest path MRT scheme where the active beamforming at the AP is steered toward the path of the smallest channel power gain among the cascaded channels of all groups and the APuser direct channel, i.e., w = h min h min with h min = arg min{ h 1 2 , . . . , h L 2 , h d 2 }. We observe that the AO algorithm outperforms the other schemes, which can be well understood by the fact that the objective of problem (P1) is to maximize the received SNR. Moreover, the performance of the strongest path MRT scheme is better than that of the AP-user direct path MRT scheme and the weakest path MRT scheme.

It is worth noting that the complexity of the strongest path MRT scheme is much lower than that of the AO algorithm.

VI. CONCLUSIONS

In this paper, we proposed the RIS-QRM scheme to achieve joint PBIT for RIS-aided MISO wireless communication assuming the direct link between the AP and user is not blocked. An optimization problem was formulated to jointly optimize the active beamforming at the AP and the phase shifts of the RIS for improving the received signal power at the user. Computer simulation results corroborated the effectiveness of the RIS-QRM scheme. In particular, the RIS-QRM scheme was shown to be able to improve the error performance of the additional bits delivered by the RIS without deteriorating that of the bits carried on the constellation symbol, as compared to existing ON/OFF-based schemes.

  The work was supported in part by the National Natural Science Foundation of China under Grant 61871190, in part by the National Natural Science Foundation of China under Grant U1701265, in part by the Natural Science Foundation of Guangdong Province under Grant 2018B030306005, in part by the Pearl River Nova Program of Guangzhou under Grant 201806010171, and in part by the Fundamental Research Funds for the Central Universities under Grant 2019SJ02.

Fig. 1 .

 1 Fig. 1. An illustration of the proposed RIS-QRM system, where the RIS focuses reflected signals towards the user while conveying additional information by the RIS element partition. By adjusting the phase shifts at the RIS, the signals reflected by the first subset differ from those by the second subset by a relative phase shift of π/2.

  L denote the elementwise baseband channels of the AP→RIS and RIS→user links, respectively, where G l ∈ C L1×N and h H r,l ∈ C 1×L1 with l = 1, 2, . . . , L denote the corresponding elementwise baseband channels associated with the l-th group. Let θ β 1 e jθ1 , β 2 e jθ2 , . . . , β L e jθ L H denote the reflection coefficients of the L groups, where j

  the RIS can be conveyed by group partitioning. Let B {B 1 , B 2 , . . . , B B } denote arbitrary B = 2 b group-partition realizations, where B k {I k , L \ I k } denotes the k-th group-partition realization with I k representing the indices of the p (in-phase) groups and L \ I k representing the indices of the remaining L -p groups, k = 1, 2, . . . , B. According to the incoming random bits, one group-partition realization is selected, e.g., {I, L \ I} ∈ B, such that the phase shift of the RIS elements in the l-th group (l = 1, 2, . . . , L) is given by

Fig. 2 .

 2 Fig. 2. Error performance of different reflection modulation schemes versus the SNR with N = 1, L = 2, p = 1, γ 3 = 0, and 256-PAM.

Fig. 3 .

 3 Fig. 3. Error performance of different reflection modulation schemes versus the SNR with N = 4, L = 4, p = {1, 2, 3}, and 256-PAM.

Fig. 4 .

 4 Fig. 4. BER of the RIS-QRM scheme using different active beamforming designs versus the SNR assuming N = 4, L = 4, p = 3 and 256-PAM.


Yan Zhang

Senior Member, IEEE Jiayi Zhang
email: jiayizhang@bjtu.edu.cn

Fellow, IEEE Marco Di Renzo
email: marco.di-renzo@universite-paris-saclay.fr

Huahua Xiao
email: xiao.huahua@zte.com.cn

Bo Ai
email: boai@bjtu.edu.cn

Performance Analysis of RIS-Aided Systems With Practical Phase Shift and Amplitude Response

Keywords: Performance analysis, practical phase shift and amplitude response, reconfigurable intelligent surface. I

Reconfigurable intelligent surfaces (RISs) have drawn significant attention due to their capability of controlling the radio environment and improving the system performance. In this paper, we study the performance of an RIS-assisted single-input single-output system over Rayleigh fading channels. Differently from previous works that assume a constant reflection amplitude, we consider a model that accounts for the intertwinement between the amplitude and phase response, and derive closed-form expressions for the outage probability and ergodic capacity. Moreover, we obtain simplified expressions under the assumption of a large number of reflecting elements and provide tight upper and lower bounds for the ergodic capacity. Finally, the analytical results are verified by using Monte Carlo simulations.

and low-cost scattering elements that can effectively control the wireless propagation environment [START_REF] Zhang | Prospective multiple antenna technologies for beyond 5G[END_REF]. By intelligently adapting the phase shifts and the amplitude response of the scattering elements of an RIS, the signals reflected from it can be added constructively or destructively with other signals so as to enhance the signal strength or to suppress the co-channel interference at the receiver [START_REF] Wu | Beamforming optimization for wireless network aided by intelligent reflecting surface with discrete phase shifts[END_REF]- [START_REF] Yang | On the performance of RIS-assisted dual-hop UAV communication systems[END_REF]. Thanks to these properties, RISs are considered to be a promising candidate technology for future wireless communication systems.

Several works have investigated the performance of RISassisted wireless systems [START_REF] Yang | Coverage, probability of SNR gain, and DOR analysis of RIS-aided communication systems[END_REF]- [START_REF] Trigui | Performance evaluation and diversity analysis of RIS-assisted communications over generalized fading channels in the presence of phase noise[END_REF]. In [START_REF] Yang | Coverage, probability of SNR gain, and DOR analysis of RIS-aided communication systems[END_REF], the authors studied the coverage, the delay outage rate, and the probability of the signal-to-noise-ratio (SNR) gain of an RIS-assisted communication system over a Rayleigh fading channel by using the central limit theorem (CLT). In [START_REF] Ferreira | Bit error probability for large intelligent surfaces under double-nakagami fading channels[END_REF], exact and accurate approximated expressions for the bit error rate (BER) were derived over a Nakagami-m fading channel. In [START_REF] Li | Ergodic capacity of intelligent reflecting surface-assisted communication systems with phase errors[END_REF] and [START_REF] Boulogeorgos | Ergodic capacity analysis of reconfigurable intelligent surface assisted wireless systems[END_REF], the authors studied the ergodic capacity (EC) of an RIS-assisted communication system. In [START_REF] Li | Ergodic capacity of intelligent reflecting surface-assisted communication systems with phase errors[END_REF], in particular, the impact of phase errors was analyzed. In [START_REF] Guo | Outage probability analysis and minimization in intelligent reflecting surface-assisted MISO systems[END_REF], the outage probability (OP) was computed and minimized by optimizing the phase shifts of an RIS over Rician fading channels. In [START_REF] Atapattu | Reconfigurable intelligent surface assisted two-way communications: Performance analysis and optimization[END_REF], an RIS was considered for assisting the communication between two users, and the OP and spectral efficiency were studied by using a Gamma approximation over Rayleigh fading channels. In [START_REF] Trigui | A comprehensive study of reconfigurable intelligent surfaces in generalized fading[END_REF], exact expressions of the OP and EC for an RIS-assisted system over Fox's H fading channels were provided. In [START_REF] Badiu | Communication through a large reflecting surface with phase errors[END_REF], the authors analyzed the impact of phase noise on the BER over Rayleigh fading channels. In [START_REF] Zhang | Reconfigurable intelligent surfaces assisted communications with limited phase shifts: How many phase shifts are enough?[END_REF], the authors quantified the impact of discrete phase shifts on the achievable rate over Rician fading channels. In [START_REF] Zhou | Spectral and energy efficiency of IRS-assisted MISO communication with hardware impairments[END_REF], the authors analyzed the impact of phase noise and hardware impairments for transmission over line-of-sight (LoS) channels. In [START_REF] Xu | Reconfigurable intelligent surfaces-assisted communications with discrete phase shifts: How many quantization levels are required to achieve full diversity?[END_REF], the authors studied the impact of discrete phase shifts for achieving the full diversity order over Rayleigh fading channels. In [START_REF] Xu | Ergodic secrecy capacity of RIS-assisted communication systems in the presence of discrete phase shifts and multiple eavesdroppers[END_REF], the authors investigated the ergodic secrecy capacity in the presence of discrete phase shifts and phase noise. In [START_REF] Yang | Secrecy performance analysis of RIS-aided wireless communication systems[END_REF], the authors studied the secrecy outage probability of an RIS-assisted communication system over Rayleigh fading channels. In [START_REF] Yang | Accurate closed-form approximations to channel distributions of RIS-aided wireless systems[END_REF], the authors introduced a tight approximation for the distribution of the SNR in RIS-assisted communications that is formulated in terms of the squared K G distribution. In [START_REF] Qian | Beamforming through reconfigurable intelligent surfaces in single-user MIMO systems: SNR distribution and scaling laws in the presence of channel fading and phase noise[END_REF], the authors proposed a general framework to calculate the distribution of the SNR in multiple-antenna RIS-assisted systems in the presence of phase noise. In [START_REF] Huang | Reconfigurable intelligent surfaces for energy efficiency in wireless communication[END_REF], the authors maximized the energy efficiency of an RIS-assisted downlink multi-user system by jointly designing the power allocation and phase shifts. In [START_REF] Huang | Multi-hop RIS-empowered terahertz communications: A DRL-based hybrid beamforming design[END_REF], the authors proposed a hybrid beamforming scheme for multi-hop RIS-assisted communication systems in order to improve the coverage range in the terahertz frequency band. In [START_REF] Trigui | Performance evaluation and diversity analysis of RIS-assisted communications over generalized fading channels in the presence of phase noise[END_REF], the authors analyzed the impact of phase noise on the OP of RIS-assisted communication systems over generalized fading channels.

Although there exist several contributions that analyze the performance of RIS-assisted systems, most of them are based on approximations. For instance, the CLT approximation is frequently used but it is only accurate when the number of reflecting elements of the RIS is sufficiently large. Additionally, to the best of the authors' knowledge, no analytical studies have been conducted to investigate the system performance of RIS-assisted transmission by taking into account the interplay between the amplitude and phase response of each reflecting element of the RIS [START_REF] Abeywickrama | Intelligent reflecting surface: Practical phase shift model and beamforming optimization[END_REF]. The above mentioned works, in particular, assume that the amplitude response of each reflecting element is independent of the applied phase shift, which is not always possible as discussed in [START_REF] Abeywickrama | Intelligent reflecting surface: Practical phase shift model and beamforming optimization[END_REF] and [START_REF] Gradoni | End-to-end mutual-coupling-aware communication model for reconfigurable intelligent surfaces: An electromagnetic-compliant approach based on mutual impedances[END_REF]. In [START_REF] Gradoni | End-to-end mutual-coupling-aware communication model for reconfigurable intelligent surfaces: An electromagnetic-compliant approach based on mutual impedances[END_REF], in particular, it is shown that the amplitude and phase responses are, in general, intertwined and depend on the circuital model of the tuning circuit that controls each reflecting element of the RIS.

Motivated by these considerations, in this paper, we present a detailed performance analysis of RIS-assisted single-input single-output (SISO) systems over Rayleigh fading channels, by taking into account the intertwinement between the amplitude and phase response of the reflecting elements of the RIS. In particular, exact closed-form expressions for the OP and EC are derived by considering the amplitude and phase model empirically derived in [START_REF] Abeywickrama | Intelligent reflecting surface: Practical phase shift model and beamforming optimization[END_REF]. In addition, in order to gain design insights, simplified expressions are obtained in some asymptotic regimes. Furthermore, we characterize the impact of key parameters on the system performance. The obtained findings show that an increase of the transmit SNR, the number of elements of the RIS as well as the minimum reflecting amplitude help improve the performance. It is shown, on the other hand, that the performance degrades with an increase of the steepness of the amplitude response. The main contributions of this paper can be summarized as follows.

r We provide a new analytical framework for the perfor- mance analysis of RIS-assisted systems. Considering a practical model for the phase shift and amplitude response, we derive exact closed-form expressions for the OP and the EC.

r In order to get additional insights into the impact of system parameters, we present asymptotic expressions for the OP and EC under the assumption of a large number of reflecting elements. In addition, we derive upper and lower bounds for the EC. The accuracy of the bounds becomes tighter when the number of reflecting elements of the RIS grows large.

r Capitalizing on the achieved analytical results, we analyze the impact of the transmit SNR, the number of elements of the RIS, the minimum reflecting amplitude as well as the steepness of the amplitude response on the RIS-assisted system performance. The remainder of this paper is organized as follows. In Section II, we introduce the system and channel models. In Section III, exact closed-form expressions for the OP and the EC are obtained. Also, asymptotic expressions when the number of reflecting elements of the RIS is sufficiently large are obtained. Moreover, tight upper and lower bounds for the EC are derived. In Section IV, some numerical and simulation results are presented to confirm the accuracy of the derived expressions. Finally, Section V concludes the paper and reports some future research directions.

II. SYSTEM MODEL

As illustrated in Fig. 1, we consider a SISO1 system in which a single-antenna base station (BS) serves a single-antenna user. The communication is assisted by an RIS which is installed on the facade of a building. The RIS comprises N nearly-passive reflecting elements and one micro-controller. The N reflecting elements are assumed to be spaced half of the wavelength apart. The signals transmitted from the BS are appropriately reflected by the RIS by dynamically adjusting the phase shifts of the reflecting elements.We assume that the direct link from the BS to the user is blocked due to the presence of obstacles, such as buildings. In addition, we make the assumption that full channel state information (CSI) is available at the BS in order to optimize the phase shifts. Also, we consider only the signal reflected by the RIS the first time and ignore the signals reflected by the RIS two or more times.

We denote by h i Δ = α i e -jϕ i (i = 1, 2, . . . , N) the baseband equivalent fading channel between the BS and the i-th element of the RIS, and by g i Δ = β i e -jψ i the fading channel between the i-th element of the RIS and the user. Furthermore, h i and g i are assumed to be independent and identically distributed (i.i.d.) complex Gaussian random variables (RVs) with zero mean and σ 2 0 = 1/2 variance. Therefore, the magnitudes of h i and g i (i.e., α i and β i ) follow the Rayleigh distribution and the phases ϕ i and ψ i are uniformly distributed in [-π, π). Hence, the received signal y at the user can be expressed as

y = P s N i=1 h i v i g i x + n, ( 1
)
where P s is the transmit power at the BS, x is the transmit signal with unit energy, and n is the zero-mean additive white Gaussian noise (AWGN) whose variance is N 0 . Additionally,

v i Δ = ρ i (φ i )e jφ i
is the reflection coefficient applied by the i-th reconfigurable element of the RIS, where φ i ∈ [-π, π) denotes the induced phase shift and ρ i (φ i ) ∈ [0, 1] denotes the induced amplitude which is dependent on the phase shift. In general, the relation between the phase and amplitude applied by the reconfigurable elements of the RIS depends on the load or surface impedances as elaborated in [START_REF] Renzo | Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead[END_REF] and [START_REF] Gradoni | End-to-end mutual-coupling-aware communication model for reconfigurable intelligent surfaces: An electromagnetic-compliant approach based on mutual impedances[END_REF]. In this paper, we consider the closed-form empirical model that was developed in [START_REF] Abeywickrama | Intelligent reflecting surface: Practical phase shift model and beamforming optimization[END_REF]. According to [START_REF] Abeywickrama | Intelligent reflecting surface: Practical phase shift model and beamforming optimization[END_REF], in particular, the amplitude response ρ i (φ i ) can be explicitly expressed, as a function of the phase shift, as

ρ i (φ i ) Δ = (1 -κ min ) sin (φ i -ϑ) + 1 2 ξ + κ min , (2)
where κ min ≥ 0 is the minimum amplitude, ϑ ≥ 0 is the horizontal distance between -π/2 and κ min , which corresponds to the difference between -π/2 and the phase shift φ i that minimizes sin(φ iϑ) (i.e., when φ i = -π/2 + ϑ), and ξ is the steepness of the function curve. It is worth noting that (2) yields an ideal phase shift model (i.e., ρ i (φ i ) = 1) if κ min = 1 or ξ = 0. From ( 1) and ( 2), the instantaneous end-to-end SNR at the user can be expressed as

γ = N i=1 α i β i ρ i (φ i ) e j(φ i -ϕ i -ψ i ) 2 γ t , (3)
where | • | denotes the absolute value, γ t Δ = P s /N 0 represents the SNR. By setting φ i = ϕ i + ψ i for i = 1, . . . , N to fully compensate the channel phases 2 , (3) simplifies to

γ = A 2 γ t , (4)
where

A Δ = N i=1 α i β i ρ i (ϕ i + ψ i ).
It is worth mentioning that the SNR in (4) is not necessarily the maximum achievable SNR, since ρ(•) depends on the phase shift [START_REF] Abeywickrama | Intelligent reflecting surface: Practical phase shift model and beamforming optimization[END_REF]. The SNR in (4) provides information on the actual SNR that is obtained by ignoring the interplay between the phase and amplitude response of the reflecting elements of the RIS at the optimization stage. The objective of the present paper is, in 2 As for an RIS-assisted communication system in the presence of direct link, the phase shifts at the RIS can be designed in a way that the reflected signals from the RIS and the direct signal are co-phased [START_REF] Yang | Secrecy performance analysis of RIS-aided wireless communication systems[END_REF]. More specifically, the approach proposed in this paper can be applied to analyze the performance of RIS-assisted communication systems in the presence of the direct link. This is left to a future research work. fact, to assess the performance of RIS-assisted systems under this considered mismatched design and to quantify the impact of ignoring the interplay between the amplitude and the phase response.

III. PERFORMANCE ANALYSIS

A. Outage Probability

Define Z i Δ = α i β i , the probability density function (PDF) of Z i is given by [33, Eq. (3)] f Z i (z) = 4zK 0 (2z), (5)
where

K ν (.)
f Z i (z) = 2G 2,0 0,2 z 2 - 1 2 , 1 2 . ( 6
) Define R i Δ = ρ i (ϕ i + ψ i ) = (1 -κ min )( sin(ϕ i +ψ i -ϑ)+1
2

) ξ + κ min . It is not difficult to show that ϕ i + ψ i is uniformly dis- tributed in [-π, π). Since R i is a function of ϕ i + ψ i ,
the PDF of R i for κ min = 1 and ξ = 0 can be formulated by using the transformation method between two RVs [36, Eq. (2.1.49)] as

f R i (r) = r-κ min 1-κ min 1/ξ-1 πξ (1 -κ min ) r-κ min 1-κ min 1/ξ -r-κ min 1-κ min 2/ξ , ( 7
)
where r ∈ (κ min , 1). If κ min = 1 and/or ξ = 0, then R i = and the PDF of R i is independent of ϑ.

Corollary 1: The cumulative distribution function (CDF) of the end-to-end SNR for an RIS-assisted system is given in [START_REF] Liu | Reconfigurable intelligent surfaces: Principles and opportunities[END_REF] at the bottom of the page, where H[•, . . . , •] is the multivariable Fox's H-function [37, Eq. (A.1)]. In [START_REF] Liu | Reconfigurable intelligent surfaces: Principles and opportunities[END_REF], in particular, we define τ

Δ = k mod 3 and have ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ τ = 1 : m k = 0, n k = 2, p k = 2, q k = 0; and c (k) 1 = 0, γ (k) 1 = 1, c (k) 2 = 0, γ (k) 2 = 1; z k = √ γ t κ min √ γ τ = 2 : m k = 0, n k = 1, p k = 1, q k = 1; and c (k) 1 = 0, γ (k) 1 = 1 2 ; d (1) 1 = 1 2 , δ (1)
1 = 1 2 , z k = 1 τ = 0 : m k = 0, n k = 1, p k = 1, q k = 0; and c (k) 1 = 0, γ (k) 1 = 1 2ξ ; z k = a -1 2ξ F γ (γ) = 1 √ πξ N H 0,2N :m 1 ,n 1 :...:m k ,n k :•••:m 3N ,n 3N
2N,N +1:p 1 ,q 1 :...:p k ,q k :...:p 3N ,q 3N × a j ; α

(1) j , . . . , α (k) j , . . . , α (3N ) j 1,2N : c (1) j , γ (1) j 1,p 1 ; . . . ; c (3N ) j , γ (3N ) j 1,p 3N (0; -2, 0, 0, . . . , -2, 0, 0) b j ; β (1) j , . . . , β (k) j , . . . , β (3N ) j 1,N : d (1) j , δ (1) j 1,q 1 ; . . . ; d (3N ) j , γ (3N ) j 1,q 3N +1 z 1 , . . . , z k , . . . , z 3N (8)
In particular, if j is even, i.e., j can be expressed as j = 2ω(ω = 1, 2, . . . , N), we have a j = 0, α

(3ω-1) j = -1, α (3ω) j = -1, and α (k) j = 0 for k = 3ω -1, 3ω. Otherwise, if j is odd, i.e., j = 2ω -1, we have a j = 1, α (3ω-2) j = -2, α (3ω) j = -1 2ξ , for k = 3ω -2, 3ω, α (k) j = 0. Furthermore, if k mod 3 = 1, we have b = 1 and β (k) = -2, for other k, β (k) = 0. Also, b j = -1, β (3ω-1) j = -1, β (3ω) j = -1 and β (3ω-2) j = 0. Proof: See Appendix A.
The multivariate Fox's H-function in (8) can be calculated by using widely used mathematical software tools, such as Mathematica [START_REF] Lei | Secrecy capacity analysis over α -μ fading channels[END_REF], [START_REF] Zhang | Performance analysis of 5G mobile relay systems for high-speed trains[END_REF] and Matlab [START_REF] Chergui | Rician K-factor-based analysis of XLOS service probability in 5G outdoor ultra-dense networks[END_REF]Appendix C]. The Python implementation of this function is provided in [START_REF] Alhennawi | Closed-form exact and asymptotic expressions for the symbol error rate and capacity of the H-function fading channel[END_REF].

The OP is defined as the probability that the instantaneous end-to-end SNR γ falls below a given threshold γ th . From ( 8), shown at the bottom of the next page, the OP can be directly obtained by setting γ = γ th , as follows

P out = Pr N i=1 Z i R i ≤ γ th γ t = F γ (γ th ) . ( 9
)
1) Large N Approximation: For a sufficiently large number of reflecting elements, i.e., N 1, according to the CLT, A converges to a Gaussian distributed RV. In this case, a simplified expression of the OP can be obtained by computing the mean and variance of A. To this end, we note that the mean and variance of Z i are E(Z i ) = π/4 and Var(Z i ) = (16π 2 )/16 [START_REF] Basar | Wireless communications through reconfigurable intelligent surfaces[END_REF], respectively. Accordingly, we provide the following corollary.

Corollary 2: The mean value ε and the variance δ of R i are given by

ε Δ = E (R i ) = 1 -κ min π B ξ + 1 2 , 1 2 + κ min , ( 10
)
δ Δ = Var (R i ) = (1 -κ min ) 2 π B 2ξ + 1 2 , 1 2
- 1 -κ min π 2 B 2 ξ + 1 2 , 1 2 , ( 11
)
where B(

E(A) = Nπε/4, Var(A) = N δ + ε 2 1 -π 2 /16 .
Therefore, A 2 can be approximated with a non-central chisquare RV with one degree of freedom whose PDF is [13, Eq. ( 4)]

f A 2 (x) = 1 2η 2 x λ -1 4 exp - x + λ 2η 2 I -1 2 √ xλ η 2 , ( 12
)
where I ν (.) is the modified Bessel function of the first kind [34, Eq. (8.406)], λ Δ = ( Nπε 4 ) 2 , and

η 2 Δ = N (δ + ε 2 (1 -π 2 16
)). From (4), therefore, the PDF of γ can be formulated as

f γ (γ) = 1 2γ t η 2 γ γ t λ -1 4 exp - γ + λγ t 2γ t η 2 I -1 2 √ γλ √ γ t η 2 , ( 13
)
Finally, the OP can be obtained by using [42, Eq. ( 3)], F γ (γ) = γ 0 f γ (γ)dγ and ( 9), as follows

P out = 1 -Q 1 2 √ λ η , √ γ th √ γt η , ( 14
)
where

Q m (•, •) is the Marcum Q-function [42,
Eq. ( 3)]. Equation ( 14) provides a more efficient approach than (9) to calculate the OP when N → ∞. Although ( 14) is obtained under the assumption of large N , it is sufficiently tight even for moderate values of N in the low-SNR regime, as it is shown in the numerical results.

B. Ergodic Capacity

In this subsection, we analyze the EC defined as

C = E [log 2 (1 + γ)] = ∞ 0 log 2 (1 + γ) f γ (γ) dγ. ( 15
)
Corollary 3: The ergodic capacity in [START_REF] Li | Ergodic capacity of intelligent reflecting surface-assisted communication systems with phase errors[END_REF] is given in ( 16), shown at the bottom of the next page, where for j = 1, . . . , 2N , the parameters a j , α (k) j are the same as in [START_REF] Liu | Reconfigurable intelligent surfaces: Principles and opportunities[END_REF], while for other values of j, we have

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ a 2N +1 = 1, α (3ω-2) 2N +1 = -2, α (3N +1) 2N +1 = 1; for k = 3ω -2, α (k) 2N +1 = 0 a 2N +2 = 0, α (3ω-2) 2N +2 = 1, α (3N +1) 2N +2 = -1; for k = 3ω -2, α (k) 2N +2 = 0 a 2N +3 = 0, α (3ω-2) 2N +3 = 1, α (3N +1) 2N +3 = -1; for k = 3ω -2, α (k) 2N +3 = 0 Furthermore, for j = 1, . . . , N, b j and β (k) j are the same as in (8), for j = N + 1, we have b N +1 = -1, β (3ω-2) N +1 = 1, β (3N +1) N +1
= -1, and for k = 3ω -2, we have β

(k) N +1 = 0. In addition, c (k) j , γ (k) j , d (k) j and δ (k) j
are the same as in (8) but z k = ( √ γ t κ min ) 2 when k mod 3 = 1. Also,

for k = 3N + 1, we have m 3N +1 = 0, n 3N +1 = 1, p 3N +1 = 1, q 3N +1 = 0; c (3N +1) 1 = 0, γ (3N +1) 1 = -1, z k = s. Proof: See Appendix C.

1) Bounds for the Ergodic Capacity:

The exact expression of the ergodic capacity in Corollary 3 is not simple enough to gain engineering insights. Therefore, we introduce upper and lower bounds according to Jensen's inequality as follows

C lb Δ = log 2 1 + (E (1/γ)) -1 ≤ C ≤ C ub Δ = log 2 (1 + E (γ)) . ( 17
)
Authorized licensed use limited to: IEEE Editors-in-Chief. Downloaded on October 11,2021 at 16:15:24 UTC from IEEE Xplore. Restrictions apply.

Corollary 4:

The ergodic capacity of the considered RISassisted communication system is upper bounded by

E (γ) = γ t N δ + ε 2 + γ t N (N -1) π 2 ε 2 16 . ( 18
)
The lower bound for the ergodic capacity can be approximated as

C lb ≈ log 2 1 + (1/E(γ) + V ar(γ)/[E(γ)] 3 ) (-1) , ( 19
)
where

Var (γ) = Nγ 2 t 4 -ε 2 + δ 2 + 2N (N -1) γ 2 t × δ + ε 2 2 - π 4 ε 4 256 + N (N -1) γ 2 t 9π 2 νε 16 - π 2 ε 2 ε 2 + δ 4 + N (N -1) (N -2) γ 2 t π 2 ε 2 ε 2 + δ 4 - π 4 ε 4 64 , Δ = (1 -κ min ) 4 π B 4ξ + 1 2 , 1 2
+ 4κ min (1 -κ min ) 3 π B 3ξ + 1 2 , 1 2
+ 6κ 2 min (1 -κ min ) 2 π B 2ξ + 1 2 , 1 2
+ 4κ 3 min (1 -κ min ) π B ξ + 1 2 , 1 2 + κ 4 min , ( 20
)
ν Δ = 3κ 2 min (1 -κ min ) π B ξ + 1 2 , 1 2 + κ 3 min + (1 -κ min ) 3 π B 3ξ + 1 2 , 1 2 + 3(1 -κ min ) 2 κ min π B 2ξ + 1 2 , 1 2 . Proof: See Appendix D.
When ξ is fixed, we rewrite [START_REF] Du | Millimeter wave communications with reconfigurable intelligent surfaces: Performance analysis and optimization[END_REF]

as ε = f (κ min ) Δ = (1 - Γ(ξ+ 1 2 ) √ πΓ(ξ+1) )κ min + Γ(ξ+ 1 2 )
√ πΓ(ξ+1) using [34, Eq. (8.384.1)]. Differentiating f (κ min ) with respect to κ min , we obtain f (κ min ) = 1 -Γ(ξ+ 1 2 )

√ πΓ(ξ+1) . Then, let us define

f (ξ) Δ = Γ(ξ+ 1 2 )
√ πΓ(ξ+1) . Its derivative can be expressed as

f (ξ) = Γ(ξ+ 1 2 )
√ πΓ(ξ+1) (ψ(ξ + 2 )ψ(ξ + 1) can be written as

∞ k=0 ( 1 ξ+1+k -1 ξ+ 1 2 +k
), from which we conclude that f (ξ) < 0 always holds. As such, we obtain f (ξ) < f(0) = 1 and f (κ min ) > 0, i.e., ε increases monotonically for κ min ∈ (0, 1) when ξ is a constant. Similarly, it can be proved that δ increases when κ min increases and ξ is fixed. In addition, when κ min remains unchanged, both ε and δ decreases as ξ increases. As a result, it can be observed form [START_REF] Atapattu | Reconfigurable intelligent surface assisted two-way communications: Performance analysis and optimization[END_REF] that the upper bound of the EC increases as N increases. Furthermore, increasing κ min or decreasing ξ with the other parameters being fixed, which corresponds to an increase of the reflection amplitude, further increases the upper bound of the EC. Moreover, it can be proved that the lower bound of the EC is improved when κ min increases or ξ decreases.

2) Large N Approximation: Corollary 5: For a large number of the reflecting elements, the EC can be approximated as

C = Bπ 2 √ η ln 2 γ 3 8λγ t 3 -1 4 exp - λ 2η 2 H 0,1:1,0:1,2 1,0:1,3:2,2 × 1 4 ; 1, 1 - 1 4 , 1 -1 4 , 1 , 1 4 , 1 , 1 4 , 1 (1, 1) , (1, 1) (1, 1) , (0, 1) λ 2η 2 , 2γ t η 2 . ( 21
)
Proof: See Appendix E. Although [START_REF] Zhang | Reconfigurable intelligent surfaces assisted communications with limited phase shifts: How many phase shifts are enough?[END_REF] is obtained under the assumption of large N , the asymptotic result provides a good tightness even for moderate values of N , as substantiated by the numerical results.

IV. NUMERICAL RESULTS

In this section, Monte Carlo simulations are presented to verify the analytical expressions obtained in the previous section. The outage threshold γ th in Fig. 2 and Fig. 3 is set equal to 10 dB.

In Fig. 2, the OP of the considered RIS-assisted system is shown for different values of N and κ min . It can be observed from this figure that the analytical results closely overlap with the simulated results, which proves the accuracy of the analysis. Moreover, as expected, the system performance can be improved

C = 1 2s ln 2 1 √ πξ N H 0,2N +3:m 1 ,n 1 :...:m k ,n k :...:m 3N +1 ,n 3N +1 2N +3,N +2:p 1 ,q 1 :...:p k ,q k :...:p 3N +1 ,q 3N +1 × a j ; α (1) j , . . . , α (k) j , . . . , α (3N +1) j 1,2N +3 : c (1) j , γ (1) j 1,p 1 ; . . . ; c (3N +1) j , γ (3N +1) j 1,p 3N
(1; -2, 0, 0, . . . , -2, 0, 0) b j ; β by increasing N . Furthermore, the larger κ min is, the smaller the OP. We can also observe that the system performance improves as the SNR increases. In addition, we observe that the CLT approximation is accurate even for moderate values of N in the low SNR regime. In Fig. 3, we analyze the OP as a function of N , κ min and ξ. From Fig. 3, we observe that increasing N and κ min improves the OP. On the other hand, the OP degrades by increasing ξ. In the considered setup, based on Fig. 3, we can compare the impact of N , κ min and ξ on the OP. Take the transmit SNR of 30 dBm for instance, we consider the OP when N = 8, κ min = 0.2 and ξ = 1.5 as a reference. When we set N = 16 with fixed values of κ min and ξ, the OP decreases by 83%. Similarly, the OP decreases by 49% when we double the value of κ min with fixed values of N and ξ. Finally, the OP increases by 39% when the value of ξ is doubled with fixed values of N and κ min . Therefore, we observe that N has the greatest impact on the system performance, followed by κ min , and finally by ξ.

(1) j , . . . , β (k) j , . . . , β (3N +1) j 1,N +1 : d (1) j , δ (1) j 1,q 1 ; . . . ; d (3N +1) j , γ (3N +1) j 1,q 3N +1 z 1 , . . . , z k , . . . , z 3N +1 (16)
Fig. 4 illustrates the tightness of the upper and lower bounds of the EC. We observe that the performance gap between the upper and lower bounds and Monte Carlo simulations decreases as the number of reflecting elements N increases, which confirms the tightness of ( 18) and [START_REF] Trigui | A comprehensive study of reconfigurable intelligent surfaces in generalized fading[END_REF]. Additionally, we observe that the EC increases as the SNR γ t increases when N is fixed. For example, if N = 32, setting the SNR to γ t = 20 dB results in 69.5% of improvement of the upper bound compared to γ t = 15 dB. Also, if SNR = 15 dB, the upper bound of the EC increases by about 18.4% as N increases from 16 to 32. Fig. 5 shows the impact of different system parameters on the EC performance. It is observed that the EC increases by increasing the minimum amplitude κ min . In addition, it is observed that the larger the value of ξ, the smaller the EC is. When ξ increases, in fact, the amplitude ρ i decreases for a given phase shift φ i [START_REF] Abeywickrama | Intelligent reflecting surface: Practical phase shift model and beamforming optimization[END_REF]. Furthermore, we observe that the effect of κ min on the system performance is less significant when ξ increases. We also observe that the CLT approximation is accurate even for a moderate number of elements, such as N = 32.

V. CONCLUSION AND FUTURE WORK

In this paper, we have studied the performance of an RISassisted SISO system under Rayleigh fading channels. Considering a practical phase shift model, we derived exact closed-form expressions for the outage probability and the ergodic capacity. Furthermore, upper and lower bounds for the ergodic capacity were derived, whose accuracy increases with the number of elements of the RIS. Moreover, simplified expressions were obtained under the assumption of a large number of reflecting elements. Our analysis reveals that the system performance improves with the SNR, the number of reflecting elements, the minimum values of the reflection amplitude and a decrease of the steepness of the amplitude function.

This paper can be generalized along several potential research directions. The proposed approach can be generalized to the analysis of multiple-antenna transceivers, to the impact of imperfect CSI, and to discrete phase shifts at the RIS. Furthermore, other performance metrics and their tradeoffs could be analyzed.

APPENDIX A PROOF OF COROLLARY 1

We first consider the PDF of

M i Δ = Z i R i by using the integral f M i (m) = ∞ 0 f (z, m z ) 1
z dz which is given by with u and applying [34, Eq. (9.31.5)] with some algebraic manipulations, we obtain

f M i (m) = 2 πξ (1 -κ min ) m κ min m 1 z G 2,0 0,2 z 2 -
f M i (m) = 2 πm 1 0 1 √ u √ 1 -u G 0,2 2,0 (1 -κ min ) u ξ + κ min 2 m 2 0, 0 - du. (A.2)
Using the definition of the Meijer's G-function [34, Eq. (9.301)], we write (A.2) as

f M i (m) = 1 2πi L 1 Γ (1 + s 1 ) Γ (1 + s 1 ) κ min m 2s 1 ds × 1 0 1 √ u √ 1 -u 1 + au ξ 2s 1 du I 1 , (A.3)
where a Δ = 1-κ min κ min and L i denotes the path of the integration. By employing the integration by parts method, I 1 can be expressed as

I 1 = 2 1 0 1 √ 1 -u 1 + au ξ 2s 1 d √ u. (A.4)
Substituting √ u with t, we arrive at the following result

I 1 = 2 1 0 1 -t 2 -1 2 1 + at 2ξ 2s 1 dt. (A.5)
Utilizing the Mellin transform of (1t 2 ) -1

2 [START_REF] Bateman | Tables of Integral Transforms[END_REF], we obtain

I 1 = 1 2πi L 2 Γ 1 2 Γ s 2 2 Γ 1+s 2 2 ds 2 1 0 t -s 2 1 + at 2ξ 2s 1 dt. (A.6)
With the help of [START_REF] Bateman | Tables of Integral Transforms[END_REF], (1 + at 2ξ ) 2s 1 can be expressed as

1 + at 2ξ 2s 1 = 1 2πi L 3 a -s 3 2ξ 2ξ Γ s 3 2ξ Γ -2s 1 -s 3 2ξ Γ (-2s 1 ) t -s 3 ds 3 , (A.7)
where 0 < Re(s 3 ) < -4ξRe(s 1 ). Plugging (A.7) and (A.6) into (A.3), we obtain

f M i (m) = 1 √ πξm H 0,2:0,2:0,1:0,1 2,1:3,0:1,1:1,0 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1; -2, 0, -1 2ξ (0; 0, -1, -1) (-1; 0, -1, -1) (0, 1)(0, 1)(0, 2) - 0, 1 2 1 2 , 1 2 0, 1 2ξ - κ min m 2 , 1, a -1 2ξ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (A.8)
Based on [START_REF] Rahama | On the sum of independent fox's H-function variates with applications[END_REF], we compute the moment generating function (MGF) of M i defined as E(e -sm ) with the help of [START_REF] Gradshteyn | Table of Integrals, Series and Products[END_REF]Eq. (3.326.2)]. Since the RVs M i are independent, we can obtain the MGF of the sum A Δ = N i=1 M i as fA (s) = N i=1 fM i (s). Then the CDF F A (x) can be obtained using the inverse Laplace transform of fA (s)/s. From (4) with some algebraic manipulations, we arrive at

F γ (γ) = 1 √ πξ N 1 2πi 3N L 1 . . . L 3N N n=1 Υ (s 3n-2 , s 3n-1 , s 3n ) × 1 Γ 1 -2 N m=1 s 3m-2 √ γt κ min √ γ 2 N m=1 s 3m-2 a -1 2ξ N m=1 s 3m ds 1 . . . ds 3N ,
where we define Υ(s 3n-2 , s 3n-1 , s 3n )

Δ = Γ 2 (1+s 3n-2 )Γ( s 3n-1 2 )Γ( s 3n 2ξ )Γ(-2s 3n-2 - s 3n 2ξ )Γ(1-s 3n-1 -s 3n ) Γ( 1+s 3n-1 2 )Γ(2-s 3n-1 -s 3n )
. Finally, the proof is completed using the definition of multivariable Fox's H-function.

APPENDIX B PROOF OF COROLLARY 2

We can derive the mean of the RV R i by using E(R i ) = ∞ 0 rf R i (r)dr and ( 7), which yields

E (R i ) = 1 πξ (1 -κ min ) 1 κ min r r-κ min 1-κ min 1/ξ-1 r-κ min 1-κ min 1/ξ -r-κ min 1-κ min 2/ξ dr. (B.1)
Substituting ( s-κ min 1-κ min ) 1/ξ with t, we can rewrite (B.1) as

E (R i ) = 1 π 1 0 (1 -κ min ) t ξ + κ min 1 √ t -t 2 dt = 1 -κ min π 1 0 t ξ-1 2 (1 -t) -1 2 dt+ κ min π 1 0 t -1 2 (1 -t) -1 2 dt. (B.2)
Then we can obtain [START_REF] Du | Millimeter wave communications with reconfigurable intelligent surfaces: Performance analysis and optimization[END_REF] with the aid of [START_REF] Gradshteyn | Table of Integrals, Series and Products[END_REF]Eq. (3.191.3)]. In order to obtain the variance of R i , we first derive E(R 2 i ) following the same approach as for the derivation of (10)

E R 2 i = (1 -κ min ) 2 π B 2ξ + 1 2 , 1 2
+ κ 2 min + 2κ min (1 -κ min ) π B ξ + 1 2 , 1 2 . (B.3) Using Var(R i ) = E(R 2 i ) -(E(R i )) 2
, we obtain [START_REF] Chen | Resource allocation for intelligent reflecting surface aided vehicular communications[END_REF] to complete the proof.

APPENDIX C PROOF OF COROLLARY 3

The PDF of γ can be obtained by taking the derivative of ( 8) and then inserting it into [START_REF] Li | Ergodic capacity of intelligent reflecting surface-assisted communication systems with phase errors[END_REF]. Thus, we obtain the ergodic capacity as

C = 1 2 1 √ πξ N 1 2πi 3N × L 1 . . . L 3N N n=1 Υ (s 3n-2 , s 3n-1 , s 3n ) √ γt κ min 2 N m=1 s 3m-2 × a -1 2ξ N m=1 s 3m Γ -2 N m=1 s 3m-2 ds 1 . . . ds 3N ∞ 0 log 2 (1 + γ) (γ) - N m=1 s 3m-2 -1 dγ I 2 .
Employing [35, Eq. (07.34.03.0456.01)], we can express I 2 as Using the definition of the multivariable Fox's H-function, we obtain [START_REF] Boulogeorgos | Ergodic capacity analysis of reconfigurable intelligent surface assisted wireless systems[END_REF] to conclude the proof.

I 2 = 1 ln 2 ∞ 0 G 1,2 2,2 ⎡ ⎢ ⎢ ⎣ γ - N m=1 s 3m-2 , - N m=1 s 3m-2 - N m=1 s 3m-2 , - N m=1 s 3m-2 -1 ⎤ ⎥ ⎥ ⎦ dγ.

APPENDIX D PROOF OF COROLLARY 4

By using ( 4) and ( 17), we can rewrite [START_REF] Li | Ergodic capacity of intelligent reflecting surface-assisted communication systems with phase errors[END_REF] as The integral I 3 can be decomposed as

C ≤ C ub = log 2 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 + E ⎛ ⎝ N i=1 α i β i R i 2 γ t ⎞ ⎠ I 3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (D.1) L ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ G 1,2 2,2 ⎡ ⎢ ⎢ ⎣ γ - N m=1 s 3m-2 , - N m=1 s 3m-2 - N m=1 s 3m-2 , - N m=1 s 3m-2 -1 ⎤ ⎥ ⎥ ⎦ ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ = 1 s G 1,3 3,2 ⎡ ⎢ ⎢ ⎣ 1 s 0, - N m=1 s 3m-2 , - N m=1 s 3m-2 - N m=1 s 3m-2 , - N m=1 s 3m-2 -1 ⎤ ⎥ ⎥ ⎦ = 1 s L 3N +1' Γ - N m=1
I 3 = γ t E N i=1 (α i β i R i ) 2 + 2 N -1 i=1 N k=i+1 α i β i R i α k β k R k . (D.2)
Applying ( 10) and ( 11) into (D.2), we obtain

E N i=1 (α i β i R i ) 2 = N δ + ε 2 , (D.3) E N -1 i=1 N k=i+1 α i β i R i α k β k R k = N (N -1) π 2 ε 2 32 . (D.4)
Thus, by substituting (D.2), (D.3) and (D.4) into (D.1), we obtain E(γ) and the upper bound for the ergodic capacity is derived as in [START_REF] Atapattu | Reconfigurable intelligent surface assisted two-way communications: Performance analysis and optimization[END_REF].

Next, according to [START_REF] Zhang | Power scaling of uplink massive MIMO systems with arbitrary-rank channel means[END_REF], we apply the Taylor series expansion of 1/γ around E(γ) and obtain E (1/γ) ≈ 1/E (γ) + Var (γ)/[E (γ)] 3 .

(D.5)

The variance Var(γ) can be obtained as follows. From (4), we can write

Var (γ) = Var N i=1 (α i β i R i ) 2 +2 N -1 i=1 N k=i+1 α i β i R i α k β k R k . Using Var( N i=1 X i ) = N i=1 Var(X i ) + 2 N -1 i=1 N k=i+1 Cov(X i , X k ), Var(X i ) = E(X 2 i ) -E(X i )E(X i ), Cov(X i , X k ) = E(X i X k ) -E(X i )E(X k ) and E(x n ) = ∞ 0
x n f X (x)dx, the variance of γ can be expressed as in [START_REF] Badiu | Communication through a large reflecting surface with phase errors[END_REF], which completes the proof.

APPENDIX E PROOF OF COROLLARY 5

The ergodic capacity can be written as

C = B ln 2 ∞ 0 ln (1 + γ) f γ (γ) dγ. (E.1)
Substituting ( 13) into (E.1), we obtain

C = B 2η 2 γ t ln 2 1 γ t λ -1 4 exp - λ 2η 2 × ∞ 0 ln (1 + γ) exp - γ 2γ t η 2 I -1 2 √ λ √ γ t η 2 √ γ γ -

Fig. 1 .

 1 Fig. 1. RIS-assisted communication system.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Outage probability versus γ t for different N and κ min (ξ = 1.5).

Fig. 4 .

 4 Fig. 4. Upper and lower bounds for the ergodic capacity versus γ t (κ min = 0.2, ξ = 1.5).

Fig. 5 .

 5 Fig. 5. Comparison of ergodic capacity performances under different κ min and ξ (N = 32).

(C. 4 )

 4 With the help of[START_REF] Gradshteyn | Table of Integrals, Series and Products[END_REF] Eq. (9.301)] and the Laplace transform of the Meijer's G-function[START_REF]The Wolfram functions site[END_REF] Eq. (07.34.22.0003.01)], we obtain (C.5) as shown at the bottom of this page. By defining L(p(z)) = P (s) and utilizing the final value theorem lim t→∞ ( t 0 p(z)dz) = lim s→0 + sL( t 0 p(z)dz) = lim s→0 + s P (s) s = lim s→0 + P (s), we can rewrite the ergodic capacity as (C.6), which is given at the top of the next page, where s is a number close to zero.

  is the modified ν-order Bessel function of the second kind [34, Eq. (8.432)]. By utilizing [35, Eq. (07.34.03.0605.01)], we can rewrite (5) in terms of the Meijer's G-function [34, Eq. (9.301)] as

  •, •) is the beta function defined in [34, Eq. (8.384)], E(•) and Var(•) represent the statistical expectation and the variance, respectively. Proof: See Appendix B. Since Z i and R i are independent of each other, the mean and variance of A are, respectively, equal to

  + s 3N +1 Γ (1s 3N +1 ) Γ 2 1 +

	C =	1 2s ln 2	1 √ πξ	N	1 2πi	3N +1	L 1	. . .	L 3N +1	N n=1	Υ (s 3n-2 , s 3n-1 , s 3n )
			Γ -	N	s 3m-2 N	s 3m-2 -s 3N +1
		×		m=1				N			N	m=1
						Γ -2		s 3m-2 Γ 2 +	s 3m-2 -s 3N +1
							m=1			m=1
		×	√ γκ min	2	N m=1	s 3m-2	a	-1 2ξ	N m=1	s 3m	s s 3N +1 ds 1 . . . ds 3N +1	(C.6)

s 3m-2 + s 3N +1 Γ (1s 3N +1 ) Γ 2 1 + N m=1 s 3m-2s 3N +1 Γ 2 + N m=1 s 3m-2s 3N +1 s s 3N +1 ds 3N +1 (C.5)

Authorized licensed use limited to: IEEE Editors-in-Chief. Downloaded on October 11,2021 at 16:15:24 UTC from IEEE Xplore. Restrictions apply.

  With the aid of[START_REF]The Wolfram functions site[END_REF] Eq. (07.34.26.0008.01)] and [34, Eq. (9.31.5)], we can rewrite I -1 Furthermore, by employing [35, Eq. (07.34.03.0228.01)], [35, Eq. (07.34.03.0456.01)] and [35, Eq. (07.34.26.0008.01)], the logarithm function and the exponential function can be written in terms of the Fox's H-function as Substituting (E.3), (E.4) and (E.5) into (E.2) and utilizing the definition of the multivariable Fox's H-function, the proof is completed and (21) is proved.

							1 4 dγ.
							(E.2)
				2	(	√ γλ √ γ t η 2 ) in terms of the Fox's H-
	function as
	I -1 2	√ γλ √ γ t η 2 = πH 1,0 1,3	λγ 4γ t η 4	1 4 , 1 4 , 1 , 1 -1 4 , 1 , 1 4 , 1	.
							(E.3)
		ln (1 + γ) = H 1,2 2,2 γ	(1, 1), (1, 1) (1, 1), (0, 1)	,	(E.4)
		exp -	γ 2γ t η 2 = H 1,0 0,1	γ 2γ t η 2	-(0, 1)	.	(E.5)

The generalization of this work to multi-antenna transceivers is a promising direction for future research.

This work was supported in part by the National Key R&D Program of China under Grant 2020YFB1807201, in part by the National Natural Science Foundation of China under Grants 61971027, U1834210, and 61961130391, in part by Beijing Natural Science Foundation under Grant L202013, in part by Frontiers Science Center for Smart High-speed Railway System, in part by Royal Society Newton Advanced Fellowship under Grant NA191006, in part by the State Key Lab of Rail Traffic Control and Safety under Grants RCS2018ZZ007 and RCS2019ZZ007, in part by Open Research Fund of the State Key Laboratory of Integrated Services Networks under Grant ISN20-04, and in part by ZTE Corporation, and the State Key Laboratory of Mobile Network and Mobile Multimedia Technology. Marco Di Renzo's work was supported in part by European Commission through the H2020 ARIADNE project under Grant 871464 and through the H2020 RISE-6 G project under Grant 101017011. The review of this article was coordinated by Dr. C. Yuen.