Performance Analysis of Wireless Mesh Backhauling Using Intelligent Reflecting Surfaces

This paper considers the deployment of intelligent reflecting surfaces (IRSs) technology for wireless multi-hop backhauling of multiple basestations (BSs) connected in a mesh topology. The performance of the proposed architecture is evaluated in terms of outage and symbol error probability in Rician fading channels, where closed-form expressions are derived and demonstrated to be accurate for several cases of interest. The analytical results corroborated by simulation, show that the IRS-mesh backhauling architecture has several desired features that can be exploited to overcome some of the backhauling challenges, particularly the severe attenuation at high frequencies. For example, using IRS with four elements, <inline-formula> <tex-math notation="LaTeX">$N=4$ </tex-math></inline-formula>, provides a symbol error rate of about 10<sup>−5</sup> at a signal-to-noise ratio of about 0 dB, even for a large number of hops. Moreover, the obtained analytical results corroborated by Monte Carlo simulation show that the gain obtained by increasing <inline-formula> <tex-math notation="LaTeX">$N$ </tex-math></inline-formula> decreases significantly for <inline-formula> <tex-math notation="LaTeX">$N>5$ </tex-math></inline-formula>. For example, increasing <inline-formula> <tex-math notation="LaTeX">$N$ </tex-math></inline-formula> from 1 to 2 provides about <inline-formula> <tex-math notation="LaTeX">$8\tilde{\text {d}}\text { B}$ </tex-math></inline-formula> of gain, while the increase from 3 to 4 provides about <inline-formula> <tex-math notation="LaTeX">$48\tilde{\text {d}}\text { B}$ </tex-math></inline-formula>. Moreover, the degradation caused by the relaying process becomes negligible when the number of IRS elements <inline-formula> <tex-math notation="LaTeX">$N=$ </tex-math></inline-formula> 3.

the demand for data services is continuously increasing as shown by the latest statistics of wireless networks, which revealed that data traffic for mobile users grew 68% in 2019, reaching 38 Exabyte (EB) per month as compared to 27 EB per month in 2018 [1], [2]. Ericsson predicts that the mobile data traffic will show an upsurge of 27% annually till 2025. Network densification (NeDe) is one of the prominent techniques that can be used to improve the capacity of wireless networks, and it is an integral part of the 5G architecture [4], [5]. The main concept of NeDe is to deploy a large number of small cells with low power to complement the macro cell functionality in areas with poor signal quality, such as indoor environments and urban areas with tall buildings [6], [7]. The term small-cell covers a variety of cell types such as minicells, micro-cells, pico-cells, femto-cells, and mobile-small cell. Wireless coverage at the street level is a key solution in urban areas where the access points are placed on poles and walls [8], [9].
While indoor small-cell backhauling is typically based on wired technologies, outdoor small cells backhauling is mostly performed through wireless technologies, which may increase the macro cell base station (BS) backhauling traffic rate to several gigabit/s (Gbps) [10]. Although existing fiber infrastructure provides reliable and large capacity communications with rates that may reach terabit/s, installing new fiber infrastructure has limitations such as the high installation and maintenance cost, considerable installation time, and trenching is prohibitively expensive and might be impossible in some areas. For example, according to [11], [12], about 85% of the fiber links cost is due to trenching and installation. Consequently, high dense small cells with wireless backhauling, which is the main focus in this work, are widely considered to combat the wired solutions limitations. The main advantages of wireless over wired backhauling are the scalability, flexibility, cost efficiency, and low-complexity maintenance process. To reduce the traffic load on the macro BS, wireless mesh backhauling is indispensable, and is considered an integral technology of the current 5G, future sixth generation (6G) and 6G + cellular communication networks [8], [9]. Moreover, wireless backhauling is the key enabler for the integrated access and backhauling (IAB) technology, which is a promising approach to reduce the deployment cost in ultra-dense networks by utilizing a portion of the available radio resources for wireless backhauling [13]- [16]. Recently, wireless backhauling has been investigated and standardized 1536-1276 © 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
One of the prominent efforts on wireless backhauling is the Telecom Infra Project (TIP), which consists of more than 500 member organizations that devote their efforts to accelerate the development of new infrastructure solutions for future networks [18]- [22]. Wireless backhauling can be realized with several network topologies such as ring, tree and mesh, with the later being the most attractive because it provides backhauling with low-cost, flexible configuration, maintainable, and long distance coverage [23]- [26]. However, wireless backhauling suffers from some limitations as compared to optical fiber, particularly the relatively low capacity, disruptive interference, and severe signal attenuation at high frequencies. Therefore, several researchers have proposed transmission schemes and technologies to overcome one or more of the wireless backhauling limitations. Examples for such solutions include millimeter-wave (mmWave) transmission, free-space optical (FSO), interference management protocols [27]- [31] and massive multiple-input-multiple-output (mMIMO) [32]- [36].
The lack of spectrum supporting wide channel bandwidths has been identified as a potential bottleneck for wireless backhaul. This has given spectrum administrators the opportunity to introduce wider channels in currently used frequencies.
An additional possibility is to open new frequency bands such as the 90 GHz band [10]. The European Telecommunications Standards Institute (ETSI) and TIP recommend mmWave communications as a core technology [18]- [20]. Realistic experiments have been conducted for one Gbps average peak user throughput for a maximum range of 250 m [19]. In addition, the TIP backhaul team has proposed different methods for assessing the performance of mmWave, and provided a guidance for the installation process [21], [22].
Intelligent reflecting surfaces (IRSs), also called metasurfaces, have been introduced recently with the aim of controlling the propagation medium to enhance the quality of service (QoS) by boosting the energy and spectral efficiencies of wireless networks. The IRS technology is expected to play a significant role in the future, where smartness, energy efficiency and spectral efficiency are the main requirements for the forthcoming wireless networks. IRS applies a large number of passive antenna elements which introduce a phase-shift to the received signals, and reflect them back to the destination. For efficient transmission, multiple reflectors are used to target a certain destination, and the introduced phase shifts are selected to ensure that the reflected signals add coherently at the receiver. As a result, the signal-to-noise ratio (SNR) increases considerably, and consequently, the spectral efficiency is boosted [37]- [52]. In the context of wireless backhauling, particularly when small cells are deployed in urban areas, using IRSs is crucial to improve the signal quality in the absence of line-of-sight (LoS) connectivity between certain BSs. In this case, a virtual LoS can be created by inserting IRS panels between such BSs, where the locations of IRS panels are selected to ensure a strong LoS for transmitter-IRS and IRS-receiver paths. Such scenarios are expected when the small cells are deployed at the street level in urban areas [8]. Moreover, even if LoS exists, IRS can significantly boost the signal quality by focusing the radiated energy towards the destination BS.

A. Related Work
IRS has recently attracted extensive research attention. For example, the authors of [37], [38] presented detailed IRS technology overview, and discussed state-of-art solutions and theoretical performance limits. Energy-efficient approaches for the transmit power allocation and the phase shifts of the IRS elements is introduced in [39], where an accurate model for IRS power consumption is presented. A realistic implementation in outdoor environments has shown that the methods proposed in [39] for IRS power allocation may provide up to 300% higher energy efficiency when compared to multi-antenna and amplify-and-forward systems. Joint active and passive beamforming is considered in [40], where some recommendations are provided for optimal deployment. Other research efforts are dedicated to study the performance of IRS with other existing signalling and communication technologies such as index modulation (IM) [41], space-shift-keying (SSK) [42], and non-orthogonal multiple access (NOMA) [45]. In [43], IRS assisted MIMO is investigated, where efficient algorithms for the phase shifts at the IRS and precoding at the transmitter are proposed to minimize the symbol error rate (SER). Since the optimum values for the phase shifts depend on the instantaneous channel side information (CSI), channel modeling and estimation are considered in [44]- [47]. In [49], the far-field pathloss model is derived for IRS based links using optical physics techniques, and it is shown that each reflecting element acts as a diffuse scatterer. A practical IRS implementation is introduced in [50], where a high-gain and low-cost IRS with 256 reflecting elements is designed, in which positive intrinsic negative (PIN) diodes are used to design 2-bit phase shifters, and it was shown that a gain of 19.1 dBi can be achieved using mmWave. In [51], the joint design of the beamforming matrix and phase shift matrix at the BS and IRS is investigated using deep reinforcement learning algorithms. An overview about holographic MIMO surfaces is provided in [53] including the hardware architectures for reconfiguring the surfaces.

B. Motivation and Contribution
Future mobile networks are expected to confine more of the radio communication needs of current users, and support many new users and industries, with applications that require ultra high data rates. Examples for such applications include holographic communication, online gaming, and streaming of super high definition immersive 3D videos, etc. To be able to support such applications for a large number of users, the network should be able to handle traffic volumes in the range of terrabit/s per square kilometer, coming from a massive simultaneous connections [10]. Consequently, the backhaul traffic between the base stations (BSs) is expected to experience enormous data rate and traffic volume increase, and thus the design of a backhaul that satisfies such scenarios is indispensable. This work focuses on the wireless portion of the backhaul network where data volumes are forwarded to a nearby BS which has a fiber connection with the corenetwork. Since the link to the fiber-connected BS may have multiple hops, the analysis is generalized to multiple hop scenario, where IRS panels with N reflectors are deployed between communicating BSs in each hop to enhance the connection reliability. The obtained analytical results corroborated by Monte Carlo simulation demonstrate that the gain achieved by increasing the number of IRS elements is inversely proportional to the number of IRS elements N , particularly for N > 5. Moreover, the degradation caused by the relaying process becomes negligible when the number of IRS elements is more than 3. The contribution of this paper can be summarized as follows: 1) Although IRSs have been considered widely in the recent literature [37], [38], to the best of the authors' knowledge, there is no work reported that considers the application of IRS in wireless backhauling with multiple hops, where BSs relay the traffic until it arrives to the core network. 2) Because wireless backhauling typically requires LoS connectivity to avoid severe signal fading, the channel gain is assumed to follow the Rician fading model. 3) An accurate approximation for the probability density function (PDF) of the received SNR is derived. 4) Analytical expressions are derived for the SER and outage probability (OP) for single and multi hop scenarios. 5) Because some of the analytical results are represented in terms of infinite series, numerical results for the truncation error are presented to demonstrate the convergence behavior of the derived solution. The obtained results showed that using about 30 terms in the summations provides a truncation error of less than 10 −7 for most cases of interest.

C. Paper Organization
The rest of the paper is organized as follows. Sec. II presents the system model with IRS based wireless backhauling. A derivation for an accurate approximation for the PDF of the received SNR per hop is introduced in Sec. III. The derivations for the SER and OP are provided in Secs. IV and V, respectively. Sec. VI shows and discusses the numerical results of the proposed system whereas the conclusion and future work are provided in Sec. VII.
II. SYSTEM MODEL This work considers a heterogenous wireless network where a macro cell is overlaid with multiple small cells. The macro BS (mBS) and small-cells BS (sBS) are configured in a mesh topology. In such networks [54], each sBS might be able to communicate with multiple other sBSs and choose a particular route to backhaul its data to the core network. Fig. 1 shows a simplified mesh network with IBA using four sBSs and one mBS. As shown in the figure, some routes may only be utilized through IRS. Once a particular route is selected, the sBSs cooperate, act as relays, to route the backhaul traffic to the mBS, and then to the core network. Therefore, the traffic of a particular sBSs may arrive to the mBS through multiple hops, where the signal in each hop is decoded and forwarded to the next sBS, or to the mBS. The route selection process is typically performed to optimize certain performance parameters such as delay, interference, power, energy or load balance [55]- [59]. Therefore, the number of hops is generally a random variable which has been the subject of some studies dedicated to model the hop count [60], [61]. To provide LoS links between adjacent sBSs, IRS panels, each of which consists of N reflecting elements, are placed between each two BSs. The backhaul links are realized using a single broadband antenna at each BS. Using the region three-dimensional maps, advanced ray tracing techniques [62] can be used to predict the CSI between adjacent BSs and the IRS, which allow selecting the optimum location for the IRS panels'. The channel fading is considered to be Rician to capture the LoS signal component and other small reflections from the surrounding environment. The LoS is considered only between the BSs and the IRS, while the direct path between BSs is typically blocked by large buildings or other obstacles.
In each hop, a sBS transmits its backhaul traffic to one of the adjacent sBSs, or to the mBS, through an IRS panel, where each reflecting element introduces a phase shift ϕ n and amplitude gain β n . The process is repeated L times until the data arrive at the corenetwork or mBS. It is worth noting that an IRS panel can be shared by multiple BSs simultaneously, where the available reflecting elements can be assigned to BSs that have LoS link with the IRS panel. However, an IRS element can be assigned only to one BS to avoid interference.
In this work, we initially consider the single hop scenario, and then the model is extended to multiple hops. Consequently, the hop index is dropped unless it is necessary to include it. Given that a particular sBS transmits a data symbol s, the signal reflected from the nth IRS element can be expressed as where H n = β n e jϕn , ϕ n is the phase shift and β n ≤ 1 is the reflection coefficient of the nth IRS element, P 0 is the sBS transmission power, and the channel coefficient n | n | e jθn ∼ CN m n , 2σ 2 n . In urban environments, sBSs do not typically have LoS link with each other due to the dense and large structures in such environments. Therefore, IRSs can be deployed in specific locations such that a LoS connection is established between the IRS and each of the concerned sBSs. Therefore, the PDF of the channel envelop α n | n | can be considered Rician [63], where I 0 (·) is the modified Bessel function of the first kind and zero order, Ω n E α 2 n = μ 2 n + 2σ 2 n , μ n = |m n |, the average channel fading m n ≥ 0, and K n = μ 2 n 2σ 2 n is the Rician factor that determines the link quality, K n ∈ (0, ∞).
For small values of K n , the channel fading becomes severe, which indicates that the LoS signal component is weak. However, according to some experimental measurements, it has been shown that the received signal amplitude in urban and suburban areas follows the Rician distribution with K n 9 dB due to strong LoS propagation [64]- [69]. The first moment of α n with parameters K n and Ω n is given by [70] where E [·] denotes the expected value and 1 F 1 (·, ·, ·) is the confluent hypergeometric function of the first kind.
In the reflecting phase, the signal will go through a second fading channel h n before arriving to the destination BS. Therefore, the received signal at the destination BS is given by where the channel coefficient h n |h n | e jζn ∼ CN m n , 2σ 2 n and w is the additive white Gaussian noise . By selecting ϕ n = − (θ n + ζ n ), all the reflected signal components are added coherently in the channel, and thus, the received signal BS can be rewritten as where | n | |h n | n . Therefore, n is formed by the product of two Rician random variables, and hence, its PDF is given by [71, eq. 6.67] whereK is the Rician factor for channel h n and K m (·) is the modified Bessel function of the second kind with order m.

III. SNR DISTRIBUTION
To derive the SER and OP, the SNR distribution should be evaluated. Towards this goal, the following two subsections present the PDF derivation of the instantaneous SNR at the destination BS for a single and multiple reflectors, respectively.

A. Single Reflector (N = 1)
Based on the received signal in (5), the instantaneous SNR of the received signal at the BS for a signal reflected from the nth reflector, can be written as where P P 0 /σ 2 w . Using f n ( n ) given in (6), the PDF of γ n can be obtained by applying random variable transformation. Consequently, f γn (γ n ) can be written as where

B. Multiple Reflectors (N > 1)
Based on the received signal in (5), the instantaneous SNR for N > 1 can be written as Unlike the single reflector case, deriving a closed-form expression for the PDF of γ is intractable. Consequently, an approximate PDF will be derived to enable analytical performance evaluation. To this end, the cascaded channels for the nth reflector can be written as In general, the product of two independent Gaussian random variables is not Gaussian [72]. However, if μX σX , μY σY 1, then PDF of the product , [74]. For Rician fading channels, such condition can be satisfied when the Rician factors K , K 1 dB, and thus the approximation is very accurate for strong LoS environments. Consequently, the PDFs of Z I n and Z Q n can be approxi- Given that σ I n = σ Q n = σ n andσ I n =σ Q n =σ n , the statistics of Z I n and Z Q n are given by In addition, as can be noted from the statistics of Z I n and Z Q n in (13)- (16), Z I n and Z Q n have the same variance but different mean values. Consequently, the distribution of n can be approximated by a Rician PDF with mean and variance that are given by Consequently, the Rician PDF parameters Ω n and K n can be written as Based on (20), a cascaded channel of two Rician channels having Rician factorsK and K can be approximated as a Rician channel with Rician factor K n < min K n , K n given thatK n and K n are large. By simple random variable transformation, the PDF of λ n = √ Pβ n n can be found as Rician with the following distribution.
where K λn = K n and Ω λn = Ω n Pβ 2 n . However, the distribution of the sum of multiple Rician random variables, Λ = n λ n , does not have a closed form expression. Therefore, an accurate closed-form approximation will be used as described in [75] to enable analytical performance evaluation. It is worth noting that there is a typo in [75,Eq. 5], where the first term should be 1/ 2πσ 2 λn instead of 1/ √ 2π. To simplify the discussion, it is assumed that β n = β ∀n, and a new variable is defined asP 0 β 2 P 0 andP P 0 /σ 2 w , to simplify the notation. Finally, by simple random variable transformation for the distribution of Λ given in [75], the distribution of γ=Λ 2 can be provided as where μ λn and σ 2 λn are the mean and variance of λ n , respectively, and f c (γ) is a correction where N a 1 , and the coefficients a 0 , a 1 and a 2 are used to control the amplitude, spread and shift, respectively. In [75], the values of these coefficients depend on the Rician factor K and the number of random variables N . The method for evaluating the coefficients is introduced in [75], which is based on the least squares fitting with the exact cumulative distribution function (CDF). A table is provided in [75] for certain cases, however, these coefficients can be evaluated using nonlinear curve-fitting in least-squares sense, where the exact CDF can be obtained by convolving N Rician PDFs and then performing numerical integration.

C. The SNR Gain
The SNR improvement gained by increasing N can be measured using the effective SNR of the received signal, which can be defined for the lth hop as Substituting (12) into (24) gives where the expected value of the channel envelope is given in (3). For independent and identically distributed (i.i.d.) channels and equal reflection coefficients, β n = β ∀n, thenγ N can be expressed as N , and thus, Therefore, the effective received SNR for large number of reflectors becomes linearly proportional to N 2 . It is also worth noting that E [α] is an increasing function versus the Rician factor K, which implies that increasing K provides additional improvement toγ N,l .

IV. SYMBOL ERROR RATE (SER) ANALYSIS
To evaluate the SER of signalling over fading channels, the general method is to model the channel as conditionally Gaussian, obtain the SER, and then to eliminate the conditioning by averaging the conditional SER over the instantaneous SNR γ. Therefore, the average SERP S can be computed asP For most widely used modulation schemes such as quadrature amplitude modulation (QAM) and phase shift keying (PSK), the conditional SER for a given γ in the presence of AWGN can be approximated as where Q (·) is the complementary cumulative distribution function of the Gaussian distribution, and the values of A and B depend on the modulation scheme [76,

A. Single Reflector (N = 1), Single Hop (L = 1)
The SER for the single reflector case can be obtained by substituting (8) and (29) into (28), and defining √ γ y, which yields As can be noted that the closed form for the integral is very difficult to obtain. Therefore, approximating either the Q-function or the modified Bessel function can result in a solvable integration. However, to make the integration feasible, tight and tractable approximation for Q (x), which has been proposed in [77], is applied. The approximation is based on a truncated series and given by where the number of terms in the summation n a can be selected depending on the desired tightness. Substituting (31) in (30) yields The integral can be solved as described in [78, eq. (2.16.8.4)], which yields the SER in (33), as shown at the bottom of the page, where W is the Whittaker hypergeometric function.

B. Multiple Reflectors (N ≥ 2), Single Hop (L = 1)
The average SERP S for this case can be derived by substituting (29) and (22) into (28), which after some manipulations, as shown in Appendix I, gives where the integrals T 1 , T 2 and T 3 are evaluated in Appendix I and the final solutions are given by 1 NBP , and D −n (·) is the parabolic cylindrical function.

C. Single and Multiple Reflectors (N > 1), Multiple Hops (L ≥ 2)
To extend the SER analysis for the general case where {N, L} ≥ 2, we define the vectorŝ = [ŝ 1 ,ŝ 2 , . . . ,ŝ L ], which contains the decoded symbols of each hop. Given that symbol s 0 is transmitted over a route of L total number of hops, the SER can be expressed as Pr (ŝ L = s 0 |s 0 ) Pr (s 0 ) where M is the modulation order, Pr (ŝ 1 |s 0 ) is the conditional probability for the first hop, and Pr (ŝ l |ŝ l−1 ) is the conditional probability for the remaining L − 1 hops. The vectorŝ ∈ S, where S is the set of all vectors that hasŝ L = s 0 . Evaluating (38) can be performed using the approach described in [79], [80], which for the special case of binary phase shift keying (BPSK) is given bȳ whereP S,l is the average symbol/bit error rate for the lth hop.

D. Performance Analysis With Random L
The routing protocol for a mesh network is typically designed to optimize the network parameters, i.e., maximize the total traffic, minimize the average error rate or the outage probability. In addition, most of the routing strategies permit a certain maximum limit for the number of hops based on QoS measures such as latency and time delay. Accordingly, the number of hops in a mesh network is a random variable, where the probability mass function (PMF) of the hop count depends on the system parameters and the adopted routing technique [60], [61]. For example, empirical and analytical PMFs have been compared in [60] for unit disk planar network (UDPN) and unit disk linear network (UDLN) routing protocols where the furthest neighbor (FN) and nearest neighbor (NN) strategies are considered. On the other hand, different path selection strategies for IAB are discussed in [61] including the highest quality first (HQF), wired first (WF), position aware (PA) and the maximum local rate (MLR). In this paper, the WF routing strategy is adopted [61]. According to the WF approach, if the sBS has a link with a wired mBS with SNR greater than a certain threshold, then this link is selected. Otherwise, the sBS selects the link with the highest SNR by employing HQF strategy. The same process is repeated by the receiving sBS to select the appropriate BS for the next hop. Therefore, the WF strategy tends to select paths with low number of hops. In general, given the PMF of the number of hops, the probability of error can be expressed as Pr (L)P S|L (40) where Pr (L) is the probability that the message travels over a route with L hops.
V. OUTAGE PROBABILITY The OP is defined as the probability that the SNR is below a certain threshold γ th , will be also denoted as χ for notational simplicity,P A. Single Reflector (N = 1), Single Hop (L = 1) By substituting (8) into (41), and applying the series expansion of the modified Bessel function, the OP can be found as By substituting y = γ/χ, and after some simplificationsP O,1 can be expressed as Based on [81, Eq. 6.592.2],P O can be solved in a closed-form in terms of the Meijer G function, The OP for this case is derived by substituting (22) in (41) and evaluating the integral. The resultant expression can be written asP where the expressions ofP where erf (·) is the error function,χ = Given that the signal will undergo L hops, the outage event occurs if one or more of the L hops go through an outage. Therefore, OP can be formulated as [82], [83], whereP O,l is the OP for the lth hop, which is given by (44) and (45) for N = 1 and N ≥ 2, respectively.
C. Outage Probability With Random L Similar to Sec. IV-D, the effect of a random number of hops on the outage probability can be computed as Pr (L)P O|L (50) where Pr (L) is the probability that the message travels over a route with L hops.

VI. NUMERICAL RESULTS
This section presents the analytical and simulation results for OP and SER of the considered system. The simulation results are obtained using Monte Carlo simulation, where each simulation run consists of 10 7 realizations. The links TxBS-IRS and IRS-RxBS are considered i.i.d. flat Rician fading channels where K n =K n = K, unless mentioned otherwise. The average transmission power of the TxBS P 0 and the reflection coefficient β are normalized to unity, i.e., β = P 0 = 1, and the SNR is defined as P = P 0 /σ 2 w . Figs. 2 and 3 show the SER for the single hop case, L = 1, using various values of N using BPSK and quadrature phase shift keying (QPSK), respectively, and using K = 10 dB. As can be noted from the two figures, the analytical and simulation results for the single reflector case, N = 1, match very well for the considered range of SNR because the cascaded channel PDF is exact, and the only approximation used is for the Q-function. For the multiple reflectors case, the PDF approximation of the cascaded channel resulted in  some discrepancies for the N = 2, 3 cases, at high SNRs. Such results are obtained because the approximation error is relatively more significant for N < 4, particularly at high SNRs where the effect of the AWGN is small when compared to the approximation error. As can be noted from the figures, increasing the number of reflectors N effectively improves the SNR, or equivalently, it enhances the SER considerably. Nevertheless, the improvement gained by increasing N decreases for large values of N . For example, the gain obtained using N = 2 as compared to the case of N = 1 is about 8 dB atP S = 10 −4 , increasing N from 2 to 3 provides only and additional 6 dB, and so forth. Such behavior is obtained because SNR increases nonlinearly versus N . As an example, for the special case where {β n , n } = 1 ∀n in (5), the SNR improvement becomes 2 log 10 (N ). Interestingly, increasing N to 30 provides about 36 dB of SNR improvement. The significant SNR improvement can be exploited to increase the modulation order, and thus enhance the spectral efficiency. Fig. 4 shows the analytical and simulated SER versus SNR using QPSK for various values of N and K. As can be noted from the figure, the derived approximation matches the simulation results very well for high values of K for the considered SNR range. For small K values, the difference between the approximated and simulation results is noticeable for N = 10, and it decreases by increasing N . It can be also noted that the approximation error becomes less significant at low SNRs because the performance in such regions is dominated by the AWGN. The impact of the channel quality on SER decreases by increasing N , which is demonstrated by the condensed SER curves for large values of N .  [61] with gNodeB BSs density of 30 gNB/km 2 . As can be seen from the figure, the SER for a random number of hops is bounded by the SER of L max = 1 and 6, but closer to L max = 1 because the considered routing protocol tends to select the path with minimum number of hops. Fig. 6 presents the results of the OP for various values of γ th , i.e., χ, and N . The Rician factor K = 10 dB, γ th = [5,10,15] dB, and N = [10,20]. As the figure shows, derived OP shows an excellent match with simulation results for all of the considered scenarios. Clearly, increasing γ th for a given value of N deteriorates OP, however, the degradation can be compensated by increasing N . It is interesting to note that increasing γ th by a certain value degradesP O approximately  by the same value. For example, when N = 10, increasing γ th from 5 to 10 dB increases the SNR required to achievē P O = 10 −5 from −5 to −10 dB. Fig. 7 shows the OP versus SNR for various values of K, where N = 10, L = 1, and γ th = 10 dB. As shown in the figure, the difference between the analytical and simulation results is about 0.3 dB atP O = 10 −4 , it decreases by decreasing SNR. The approximation becomes more accurate for larger values of K. The figure also shows the impact of K onP O , where increasing K can significantly improveP O , particularly when K has originally small values. For example, the SNR required to obtain P O = 10 −5 can be reduced by about 2.5 dB when K increases from 4 dB to 7 dB, while the gain is reduced to less than 1 dB when K increases from 16 dB to 20 dB.   As can be noticed from the figure, the analytical and simulation results match very well for N = 1 becauseP O is exact for this case, and for N ≥ 5 because the approximation accuracy improves versus N . Moreover, it can be noted that increasing N improvesP O significantly, and dilutes the degradation caused by increasing L. More specifically, the effect of N becomes negligible for N ≥ 10. Fig. 9 shows the normalized truncation error that results when using a small number of terms in the summations used to evaluate the BER in (33) and (34), which is defined as  where l max is the number of terms used in the summations, andP e | 50 is the probability of error for l max = 50. TheP e | 50 is taken as the reference because TE for l max > 50 ≈ 0.
The results are obtained for various values of the number of reflectors N and SNR values. As can be noted from the figure, TE decreases as l max increases, and it becomes less than 10 −7 for l max ≥ 30. Fig. 10 shows the theoretical and simulated effective SNR γ N versus N for K = 0, 10 and 30 dB, σ 2 w = −10 dB, Ω = 1, and all channels are considered i.i.d. As can be seen from the figure, the relation betweenγ N in dB and N is nonlinear due to the log function. Moreover, the figure shows that increasing K improvesγ N . Fig. 11.a and Fig. 11.b show the outage probability versus N and L, respectively, using SNR= −34.5 dB, K = 20 dB, Ω = 1, and γ th = 5 dB. The results in Fig. 11.a show that using L > 10 may provideP O > 10 −2 for N < 98. Therefore, using large number of hops should be avoided unless the number of elements in each IRS is large. Moreover, it can be noted that P O decreases sharply for N > 98. For the case of Fig. 11.b, it can be noted thatP O increases severely for L < 20, and then the increase becomes moderate. Moreover, it can be seen thatP O is very sensitive to the variations of N .

VII. CONCLUSION AND FUTURE WORK
This paper have investigated the performance of the promising IRS technology for wireless mesh backhauling, where data traffics for each BS reaches the corenetwork through multiple wireless hops. The proposed system model was analyzed for the single hop, and then generalized to multiple hops. The system performance was evaluated in terms of SNR, SER and OP were exact and approximated analytical expressions were derived for IRS systems over Rician fading channels. The analytical results were verified using Mote Carlo simulation, and the extensive comparisons confirmed the accuracy of the derived solutions, particularly when the number of reflectors N and the Rician K factors are high. In addition, the results showed that the proposed IRS based backhauling can be considered as an attractive solution for wireless backhauling because it can boost the effective SNR, and reduce OP and SER considerably. Consequently, IRS is an energy and spectrum efficiency enabler for wireless backhauling.
Our Future work will focus on investigating the IRS backhauling system in multihop scenarios where the channel estimation and compensation at each hop is imperfect, and the limitations for the timing alignment and its impact on the bandwidth will be evaluated. In addition, the design of routing protocols and path selection mechanisms based on the derived formulas will be performed to minimize the OP for a given energy constraint will be considered. The use of multi IRSs between adjacent BSs will be also considered.
The SER for N ≥ 2 is obtained by substituting (23) in (22), and thus the complete expression f γ (γ) is easily obtained. The obtained f γ (γ) and (29) are then substituted in (28), where the result is expressed in (52), as shown at the bottom of the page.
After some mathematical manipulations, the resulting integral can be expressed as where These three integrals are solved individually in the three subsections below.
A. Evaluating the Integral T 1 By substituting y = γ P , the integral T 1 in (54) is reduced to The approximation of the Q-function provided in (31) is applied, and then the resulting integral can be written as where ρ n = 1 σ 2 λn N + BP, Γ (·) is the Gamma function, and D −n (·) is the parabolic cylinder function.

B. Evaluating the Integral T 2
By substituting x = √ γ, the integral T 2 in (55) is reduced to The definition of the Q-function is then substituted in (60), i.e., Q (x) " ∞ x e − y 2 2 dy, and thus T 2 can be written as Subsequently, by changing the order of the integrals, thus T 2 can be rewritten as where T 2,1 is given by Therefore, T 2,1 should be solved first, then the result is substituted in (62) to solve T 2 . By using the change of variables Ψ (x 2 ) = (63), T 2,1 can be written as which can be solved as [81] T 2,1 = NPa 1 (2 +á) e − 1 2á 2 + Ψ 2 Consequently, (65) is substituted in (62), and thus T 2 is reduced to By using the definition of the Q-function, it can be easily shown that " ∞ 0 e − y 2 2 = π 2 . After some algebraic operations, T 2 can be simplified to (67), as shown at the bottom of the previous page. Thereafter, the integration can be evaluated in closed form as given in (68) [81, 3.462.1].
Finally, the expression in (68), as shown at the top of the page, can be written using the summation notation, where the result is expressed in (36) where δ and ν are defined below (37).

C. Evaluating the Integral T 3
The steps followed to evaluate T 3 is very quiet similar to the procedure used to evaluate T 2 , and therefore the derivation of T 3 is not included for the sake of brevity.