Intelligent Reflecting Surface Assisted Beam Index-Modulation for Millimeter Wave Communication
Sarath Gopi Sheetal Kalyani Lajos Hanzo
eess.SP


 

Intelligent Reflecting Surface Assisted Beam Index-Modulation for Millimeter Wave Communication
    Sarath Gopi and Sheetal Kalyani and Lajos Hanzo 	2020 IEEE.  Personal use of this material is permitted.  Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.
    December 30, 2023
====================================================================================================================================================================================================================================================================================================================================================================================================================================================




Millimeter wave communication is eminently suitable for high-rate
wireless systems, which may be beneficially amalgamated with
intelligent reflecting surfaces (IRS), while relying on beam-index
modulation. Explicitly, we propose three different architectures based
on IRSs for beam-index modulation in millimeter wave communication.
Our schemes are capable of eliminating the detrimental line-of-sight
blockage of millimeter wave frequencies.The
  schemes are termed as single-symbol beam index modulation,
  multi-symbol beam-index modulation and maximum-SNR single-symbol
  beam index modulation. The principle behind these is to embed the
  information both in classic QAM/PSK symbols and in the transmitter
  beam-pattern. Explicitly, we proposed to use a twin-IRS structure to
  construct a low-cost beam-index modulation scheme. We
conceive both the optimal maximum likelihood detector and a
low-complexity compressed sensing detector for the proposed
schemes. Finally, the schemes designed are evaluated through extensive
simulations and the results are compared to our analytical bounds.




    

§ INTRODUCTION

   
 Next-generation systems are expected to satisfy substantially improved specifications. Furthermore, new solutions, such as the Internet of Things (IoT), massive machine type communications (MTC) also contribute to the escalating mobile data traffic, as predicted by the International Telecommunication Union (ITU) <cit.>. Hence researchers aim for increasing the degrees of design-freedom in support of these ambitious requirements. 
   The 30-300 GHz so-called millimeter wave (mmWave) frequency band has substantial hitherto unexploited bandwidth resources for supporting Gigabit per seconds (Gb/s) data rates <cit.>. For example, in indoor scenarios a data rate of upto 6.7 Gbps is achieved by the IEEE 802.11ad standard developed at 60 GHz frequency <cit.>. This result has ignited research interest in this frequency range also for outdoor scenarios. In an early experiment, it has been shown that mmWave communication is capable of achieving a peak data rate of 1 Gbps in an outdoor environment for a communication range of upto 1.7 km at moderate Bit Error Rates (BERs) <cit.>.  This system used only 500 MHz of bandwidth at 28 GHz. Naturally, there are a number of propagation challenges to be overcome, since typically only line-of-sight (LOS) communication is possible at these frequencies, which also suffer from fading, significant absorption losses in the atmosphere and building-penetration losses <cit.> <cit.>. 
  Furthermore,  researchers are also aiming for reducing both the power consumption and hardware cost. Intelligent Reflecting Surfaces (IRS) offer a viable solution for meeting these requirements <cit.>. Explicitly, IRSs constitute passive reflecting surfaces equipped with integrated electronic circuits, which are capable of imposing carefully controlled amplitude and/or phase shifts on the incident signals <cit.>. The concept has been earlier proposed in <cit.> and its employment as a phase-shifter has become popularized by <cit.>. IRSs are eminently suitable for  energy-efficient solutions in a wide variety of applications, such as signal-to-noise-ratio (SNR) maximization <cit.>, rate-maximization <cit.>, <cit.>, for improving the energy efficiency  <cit.>, for minimizing transmit power <cit.>, for providing secure communication <cit.>, multi-cell MIMO communication <cit.>, <cit.>, over the air computation <cit.>, low latency mobile edge computing <cit.>, index modulation<cit.> and so on.
  In <cit.>, analog beamforming based beam-index modulation has been proposed as an extension of spatial modulation <cit.>. Inspired by these results, we conceive IRS assisted beam-index modulation for mmWave communications. Beamforming techniques have been exploited in mmWave communication for mitigating their path loss <cit.>, for achieving directional transmission <cit.>, for avoiding inter-carrier-interference <cit.> and also for safeguarding against eavesdroppers <cit.>. However, there is a paucity of contributions on beamforming-aided index modulation in IRS-assisted mmWave communication.  Our main contributions are:
 
 	
  * We propose IRS assisted beam index modulation for mmWave communication. Beamforming solutions proposed for mmWave frequencies tend to rely on either analog beamforming <cit.> or on hybrid techniques <cit.>.  In <cit.> digital beamforming is proposed, which relies on complex hardware. As a remedy, IRS has been proposed for imposing phase shifts on the incident signal, which can be exploited for beamforming. As a further benefit, they are capable of  circumventing the predominantly LOS nature of mmWave propagation. Hence, our proposed scheme has at least three appealing features: it supports non-LOS communication at mmWave frequencies at a low cost, whilst conveying extra information via beam-index modulation. 
 	
  * We propose three different architectures for IRS assisted beam-index modulation. The first is termed as single-symbol beam index modulation, where the information is carried both by classic QAM/PSK symbols and by the transmitter beam-pattern. This idea has also been extended for further improving the data rate in our Scheme 2, which is a multi-symbol beam-index modulation arrangement. In the third scheme, we provide an architecture for improving the SNR of the proposed beam-index modulation. 
 	
  * The optimal maximum likelihood (ML) detector is derived for the schemes conceived. Additionally, a low complexity compressed sensing assisted detector is also developed. 
 	
  * An upper bound of the average BER is obtained for the optimal ML detector. Finally, the proposed scheme is evaluated through extensive simulations and its performance is compared to the theoretically obtained bound.
 

In Table <ref>, we provide a bold
  summary and contrast our new contributions to the seminal
  literature. The key contribution of our scheme is a unique twin-IRS
  architecture, in which one of the IRSs can be positioned farther
  away from the transmitter. IM on this IRS is activated wirelessly
  using the other IRS. This architecture benefits in terms of
  accomplishing non-LOS communication by two LOS paths, additionally
  achieving a substantial beamformer gain and hence an SNR gain. 
The rest of the paper is organized as follows. Section <ref>
details the proposed IRS assisted beam-index modulation schemes. The
implementation aspects and parameter design of the schemes are
detailed in Section <ref>, while our detectors are developed
in Section <ref>. In Section <ref>, the error analysis of
the proposed scheme is provided. Our simulation results are given in
Section <ref> and we conclude in Section <ref>.  
 Notations: Throughout the paper, unless otherwise specified, bold
lower case and bold upper case letters are used to represent vectors
and matrices, respectively.  ^H, Tr{} and
λ_min() represents the hermitian, trace and minimum eigen
value of , respectively. . stands for L_2- norm. |a| and
ℛe{a} is the absolute and real value of scalar a,
respectively.  is the identity matrix of appropriate
dimension. (μ, 𝐂) represents the complex Gaussian
distribution with mean vector μ and covariance matrix
𝐂. ⌊ b ⌋ is the largest integer not greater
than b. Γ(.) is the Γ -function, i.e., Γ(z) =
∫_0^∞x^z-1e^-xdx and for integer z, Γ(z) =
(z-1)! and Γ(a,b) is the Gamma distribution with a and b
are shape and rate parameters, respectively.
   

§ PROPOSED IRS ASSISTED BEAM-INDEX MODULATION SCHEMES

   
   We propose three different IRS assisted beam-index modulation
   schemes. All these schemes are single input multiple output (SIMO)
   arrangements, containing a transmitter antenna (TA), two sets of
   IRSs and N_R receiver antennas (RAs). Each IRS has one or more
   reflecting surfaces (RS) and each RS has many elements. The first
   IRS, namely IRS_1, can be directly accessed by the transmitter
   and it is used for selectively activating the elements in the
   second IRS, i.e. in IRS_2. The block
     diagram of the proposed scheme is sketched in Fig. <ref>,
     which is further elaborated on using Fig. <ref>. The steps
     from 1 to 5 in the
     Fig. <ref> is detailed below.
    
   
   	
  * The incoming bit sequence is split into two groups. The first group is used for selecting the classic PSK/QAM symbols, while the second set is used for beam-index modulation.
   	
  * The TA and IRS_1 are kept close to each other. They have both wired and wireless connections. Based on the first group of bits, an appropriate PSK/QAM symbol (s) is selected at the transmitter, which is transmitted wirelessly to each RS in IRS_1. 		
   	
  * The wired connection is used for mapping the second group of bits onto beam-index modulation. These bits are converted to the appropriate phase vector, which are then forwarded to the elements of the RSs in IRS_1.
   	
  * Based on the received phase vector, the elements in IRS_1 impose the required phase shift on the incident signal, which are then forwarded to IRS_2. The phase is specifically adjusted for ensuring that only the desired elements in IRS_2 receive the signal. This specific selection is determined based on the information bits reserved for beam-index modulation.
   	
  * Each element of IRS_2 induce a constant phase  to reflect the signal towards the RAs. This is captured by the RAs. The information detected at the RAs includes both the conventional PSK/QAM symbols and the specific element indices of IRS_2, which reflect the symbols. This is done jointly by N_R RAs.
   
    
For detailing the schemes, we will make the following assumptions. 
   
   	
  * The channel between the TA and
          IRS_1 may be deemed to be a low-noise AWGN channel. A
          conventional horn antenna having a few centimetre length can
          be used as the TA, where IRS_1 is positioned, say 4-5
          m away from the TA <cit.>.  Hence, the
          distance satisfies the far field condition, albeit this
          is not actually necessary, since the position of the TA is
          perfectly known at IRS_1, hence near-perfect delay
          compensation can be arranged for this location. The
          TA is designed in such a way that the signal is pointed
          exactly towards IRS_1. Finally, IRS_1, which is a
          passive device, introduces a phase shift and reflects the
          signal towards IRS_2. Hence, the only source of noise,
          that can affect the signal is the one, which is added at
          IRS_1, and this is negligible.
   	
  * There is only LOS communication between IRS_1 and
          IRS_2. Typically, a Rician
          channel model is used for modelling IRS-assisted
          communication systems <cit.>. However, for
          the proposed scheme the elements in IRS_1 adjust the phase
          in such a way that it forms a directional beam and only the
          specifically selected elements of IRS_2 receive the
          signal.  Moreover, IRS_2, which is close to IRS_1, is
          carefully positioned for ensuring that there is no blockage
          between IRS_1 and IRS_2.  Hence, the channel between
          IRS_1 and IRS_2 is assumed to be an AWGN channel.
   	
  * Between IRS_2 and the receiver, we have a Rician channel model.
   
 A beneficial application of the proposed
   architectures can be found in the Internet of Things(IoT), where
   the desired information has to be collected by sensors and
   delivered to either a distant server or to a user. The various applications
   include smart homes, industrial and environmental monitoring,
   building and home automation etc. Let us assume that the
   information collected from a home or an industrial cite should be
   communicated to a BS, from where the information can be
   communicated to the destination through the cellular network. In
   this case, the TA and IRS_1, which only belong to the specific
   user, can be placed in the terrace of the home or at the industrial
   cite. The IRS_2 can be situated at the top of a tall building in
   the vicinity, which can be shared among many such users, who have
   orthogonal resources. The details of the schemes are given
 below.

   

 §.§ Scheme 1: Single-Symbol Beam-Index Modulation

   This scheme is shown in Fig. <ref>. In this scheme, both the IRSs  have only a single RS. IRS_1 is directly connected to the transmitter, whereas IRS_2 is kept at a distance, say D, from the first IRS. Let IRS_1 be a (N_1H× N_1W) element array, while IRS_2 be an (N_2H× N_2W) array and let N_1 = N_1HN_1W and N_2 = N_2HN_2W be the total number of elements in IRS_1 and IRS_2, respectively. The TA sends the symbols to IRS_1, where each element applies a specific phase shift to the incident wave so that only one of the elements in IRS_2 receives the signal. Hence, in this scheme the total number of bits per channel use (bpcu) is log_2 M + ⌊log_2 N_2 ⌋, where the first term corresponds to the QAM/PSK symbols, while the second term corresponds to the selection of the element in IRS_2.   Finally, IRS_2 reflects the signal and it is received at the RAs. 
   Let s be the transmitted symbol.  The symbol received at IRS_1 is s+w_1, where w_1 ∼(0,σ_1^2). However, under Assumption 1, we have σ_1^2  ≈ 0 and the contribution w_1 can be discarded. Therefore, the vector received at IRS_2 is: 
   
    _2 = s̱ + _2,

     
  where _2 ∼(0,σ_2^2) under Assumption 2) and $̱ is anN_2 ×1vector.  Ideally,$̱ should have only a single non-zero entry corresponding to the index of the beam (or equivalently corresponding to the selected element in IRS_2). However, this will not happen in practice, since a finite power will be dispersed on other directions also and this power distribution depends on the beampattern. The vector at the receiver can be written as:
   
    = Θ̋_2 + _R 
       = Θ̋s̱ + ,

   where $̋ is theN_R ×N_2channel matrix as defined under
   Assumption 3,Θis an(N_2 ×N_2 )diagonal matrix of phase shifts given by the elements inIRS_2and_R ∼(0,σ_R^2). Note that= Θ̋_2 +
   _Ris the additive noise component having a distribution of(0, Σ), whereΣ= ̋̋^Hσ_2^2 +
   σ_R^2. Finally, the receiver has to detect both$̱ and
   s from  to decode the transmitted bits. The detection schemes
   will be discussed in Section <ref>.  This scheme predominantly uses only a fraction of elements of the second reflecting surface, instead of exploiting
     all of them to improve the attainable beamformer gain. However,
     it should be noted that most of the transmitted energy is
     focussed on the intended elements, while the power impinging on
     all other elements is negligibly small, since beamforming is
     used in the first stage. Hence, with the aid of the proposed
     scheme, we will get the dual advantages of both a beneficial
     beamforming gain and the additional advantage of an increased
     data rate.
   

 §.§ Scheme 2: Multi-Symbol Beam-index Modulation

   In the second scheme, the first scheme is extended to multi-symbol communication. The architecture is shown in Fig. <ref>. In this case, there are N_T RSs in IRS_1 contrast to a single RS in Scheme 1. The modulator identifies N_T different phase-vectors depending on the bit sequence corresponding to the beam-index modulation and each vector is fed to different RSs in IRS_1. Therefore, N_T RSs focus the conventional QAM/PSK symbol onto N_T different elements of IRS_2. Hence, in this case, the total number of bpcu is log_2 M + ⌊log_2 N_2N_T⌋. Therefore, this scheme provides a higher data rate than scheme 1. The choice of the elements to be activated can be organized using a look up table method or the combinatoric approach <cit.>, <cit.>.
    
Mathematically, this scheme can be represented using Equations (<ref>) and (<ref>). However, the difference is that in this case, ideally there will be N_T non-zero entries in $̱.

   

 §.§ Scheme 3: Maximum-SNR Single-Symbol Beam-Index Modulation

    Fig. <ref> shows the architecture of this scheme. This is similar to Scheme 1, except that in this case each element ofIRS_2is replaced by an RS havingN_3elements. Hence, there will be a total ofN_2N_3elements inIRS_2. Both the TA andIRS_1function in the same way as in the case of single-symbol beam-index modulation. Hence, the signal received atIRS_2can be written using (<ref>).  However, in contrast to the other two cases, here the elements inIRS_2apply a phase shift to the incident signal.  The phase shift inIRS_2is adjusted in such a way that the SNR at the receiver is maximized.  
     
  
  Note that in (<ref>)$̋ and Θ is an (N_R × N_2N_3) and (N_2N_3 × N_2N_3) matrix, respectively in this case. Furthermore, Θ is not a constant matrix, but depends on $̋. The overall SNR in this case is defined as:
    
    SNR = Θ̋s̱^2/Var(),
 
    where the denominator is the variance of the norm of the vector. The SNR can be maximized by maximizing the numerator of (<ref>), since the denominator is independent ofΘ, which is the maximization variable.  Letθ_1:N_2N_3represents the entries of the diagonal ofΘ. Hence, the SNR maximization can be written as: 
   

    θ_lmax Θ̋s̱^2 s.t. |θ_l| =1 , ∀ l = 1, 2, ..., N_2N_3.

  However, the above optimization problem has the following challenges.IRS_2is a passive device and it may not be practical to solve a complex optimization problem there. Hence, the optimization should ideally be carried out at transmitter or receiver and the resultant information has to be communicated toIRS_2. Therefore, if the optimization depends on the data to be transmitted (s̱),Θhas to be updated in every time slot, which is a substantial communication overhead. Hence, the optimization should preferably only depend on either an average value ofs̱or indeed ideally should be independent of  it. Accordingly, we will propose the following solutions for (<ref>).
    
    

  §.§.§ Solution 1

    This solution is based on the assumption that an ideal beam pattern exists, i.e., all elements in the selected RS ofIRS_2receives the same power, while all other elements receive no power. Without loss of generality, let this constant be1. Hence, the optimization function in (<ref>) can be written as:
    
    Θ̋s̱^2 = Θ̋1^2 = θ^2,

 where1is a vector of 1s andθis a vector formed from the diagonal elementsθ_1:N_2N_3ofΘ. Hence, the maximization problem (<ref>) becomes:
    
    θ_lmax θ^H^̋Hθ̋ s.t. |θ_l| =1 , ∀ l = (Î-1)N_3+1, ..., ÎN_3,

whereÎis the specifically selected RS inIRS_2. Letθ_l = e^jα_l, since|θ_l| =1. Bearing this in mind and noting that^̋H$̋ is a Hermitian matrix, (<ref>) is reformulated as the following unconstrained optimization problem.
     
    α_lmax ∑_i=(Î-1)N_3+1^ÎN_3∑_j=(Î-1)N_3+1^ÎN_3ℛe{e^j(α_i-α_j)(^̋H)_ij},

  where (^̋H)_ij is the (i,j)^th element of ^̋H$̋. Since (<ref>) is not a concave function, it can only be solved using some iterative technique for finding its local maximum. 
     

  §.§.§ Solution 2

    This solution relies on the assumption thatN_R ≥N_3, i.e. there are more number of RAs than the number of elements in the RS ofIRS_2. In order to develop the solution, let us  state and prove Lemma <ref>.
      
     Let _̋Q and Θ_Q be the (N_R × N_3) and (N_3 × N_3) sub-matrices of $̋ andΘcorresponding to the selected RS, respectively and let_̱Qbe the corresponding sub-vector of$̱. If N_R ≥ N_3,  with probability 1, the bound	
      	
    _̋QΘ_Q _̱Q s^2 ≥λ_min(Θ_Q^H_̋Q^H_̋QΘ_Q ) Tr{(_̱Q s) (_̱Q s)^H}
    		
      	is non-trivial, which equivalently leads to λ_min(Θ_Q^H_̋Q^H_̋QΘ) > 0.
          	
      	 	See Appendix <ref> for proof.
      	 
 Lemma <ref> can be used  for solving the optimization problem (<ref>). The idea is to maximize the non-trivial lower bound instead of the actual function. Hence, the optimization problem (<ref>) becomes:
   
    θ_lmax λ_min(Θ_Q^H_̋Q^H_̋QΘ)
        s.t. |θ_l| =1 , ∀ l = (Î-1)N_3+1, ..., ÎN_3.

We know thatλ_min() = =1min  <cit.>. Therefore (<ref>) can be rewritten as: 
 	 
    θ_lmax min Θ_Q^H_̋Q^H_̋QΘ 
           s.t.  = 1,   |θ_l| =1 , ∀ l = (Î-1)N_3+1, ..., ÎN_3,
  
	where (<ref>) is a constrained non-linear minimax optimization problem. This can be solved directly <cit.>. Alternatively, it can be converted into a non-linear maximization problem by introducing an additional variable and then solved using standard techniques. 
      
It should be noted that in both solutions of the SNR maximization problem,Î, i.e. the selected data dependent RS ofIRS_2that has to be optimized, depends on the information bits. In order to avoid the dependence of optimization on the information bits, each RS is optimized separately whenever there is considerable change in the channel.  The optimized phase information is passed toIRS_2, which applies phase shifts to all elements instead of the selected RS. This scheme can be extended to the case of multi-symbol beam-index modulation (Scheme 2), where there will beN_TRSs inIRS_1, which activateN_TRSs inIRS_2. Finally, all the activated RSs inIRS_2can apply a phase shifts for improving the SNR. Thus the scheme will have both an improved data rate and improved SNR. Practically, the optimal phase shifts have to be estimated at the receiver and then communicated toIRS_2whenever there is significant change in the channel characteristics. 
   


§ IMPLEMENTATION OF BEAM-INDEX MODULATION

The principle behind the proposed beam-index modulation is the data-dependent activation of the elements inIRS_2. This is achieved by appropriately choosing the phase shifts applied by the elements inIRS_1. In order to estimate the phase shifts, it is assumed that there is only LOS communication betweenIRS_1andIRS_2. Therefore, the phase shifts only depend on the geometry of the pair of IRSs. The estimation of phase shifts is detailed below. 
 
Let the centre ofIRS_1be the origin co-ordinate(0,0,0)andbe the position vector of elements ofIRS_1.  Let then^thelement ofIRS_2be activated byIRS_1according to the input bit sequence and let(ϕ_n^h, ϕ_n^v)represents the azimuth and the elevation angle pair for this element with respect to the origin. Then, the phase-vector to be given by the elements ofIRS_1to choose then^thbeam isψ= 2πfτ, whereτ= ^H_̆n/cwithcandfbeing the speed of the light and the carrier frequency, and_̆n = (sinϕ_n^hcosϕ_n^v   cosϕ_n^hcosϕ_n^v   sinϕ_n^v ). Finally, in the case of multi-bit beam-index modulation, these phase shifts have to be calculated for each of IRSs according to the input bit sequences. 
 
Additionally, if the IRS elements can modify the amplitude of the incident signal along with the phase, one can modify the relative weighting of each element. Since IRSs constitute passive devices, amplification may be difficult to achieve and will not be a cost effective solution. However, attenuation can be readily applied to the incident signal <cit.>. The attenuation can be adjusted in such a way that it acts as a window function for the beamforming and the beam pattern can be accordingly modified. This can be used to reduce signal received by unwanted elements inIRS_2. In Section <ref>, the design of two IRSs is detailed. 



 §.§ Parameter Design

The parameters to be designed are the number of elements and the corresponding inter-element spacing inIRS_1, as well as inIRS_2and the distance between two IRSs.  Letλbe the wavelength corresponding to the highest frequency of operation. We will fix the design parameters as follows <cit.>.

  * Inter-element spacing in IRS_1 (d_1): The elements in each RSs of IRS_1 should be spaced at distances ≤λ/2. For example, if the maximum operating frequency is f_max = 60  GHz, then the spacing between the elements in IRS_1 is d_1 ≤ 2.5 mm. Now, if IRS_1 is a 100 × 100 element system, then its dimension is going to be as compact as 0.25 × 0.25 m.  This spacing is important for avoiding grating lobes in the beams formed using IRS_1 <cit.>. 
	
  * Distance between two IRSs (D): The distance (D) between IRS_1 and IRS_2 should meet the far field condition of D > 2L^2/λ, where L is the length of the RS <cit.>. In the above example D > 25  m. If this condition is met, it can be assumed that the wave front travelling from IRS_1 to IRS_2 is planar, so that the phase shifts can be computed as detailed in Section <ref>.
	
  * Inter element spacing in IRS_2 (d_2): The width of each element in IRS_2 (d_w) should be less than Dθ_BW, where θ_BW is the beam-width of IRS_1 and the separation d_2 between elements in IRS_2 should be higher than this value. These conditions ensure that the intended element and only the intended element receives the signal reflected by IRS_1. For a rectangular window, the approximate beam-width is θ_BW≈50λ/L <cit.>. In the example we have considered θ_BW≈1^0 corresponding to the minimum frequency. Therefore d_w < 48 cm and d_2 > 48 cm for D = 25  m. Now, if  N_1H= N_1W= 8, IRS_2 has an approximate dimension of 4  m × 4  m. 
	
  * Number of elements in IRS_1 (N_1): The number
          of elements in IRS_1 determines the length L of the
          array and its beam-width θ_BW, where these
          parameter decide the spacing between two IRSs and the
          inter-element spacing in IRS_2. A larger value of N_1 with accurate phase
            shifting and a larger value of d_2 direct the beam to
            the desired direction, hence reducing the potential beam
            misalignment problem of our beam index modulation scheme.
	
  * Number of elements in IRS_2 (N_2): This determines the data rate of the system. A large value of N_2 gives a higher data rate. However, this will make the size of IRS_2 large. Hence, N_2 is restricted by the maximum affordable array dimension. 


§ DETECTOR

 The detector has to recover the bits embedded both into the QAM/PSK symbol and the TA activation pattern inIRS_2. Explicitly, it has to detectsand$̱ from  in (<ref>). Let = s̱ and = Θ̋. Now, (<ref>) can be written as:
 
    =  + .

 It is assumed that  is known at the receiver. 
 

 §.§  Optimal Detector

 
 We first derive the optimal ML detector for the single symbol cases, i.e. for Scheme 1 and Scheme 3. Then we extent it to the multi symbol case of Scheme 2. 


  §.§.§ Single-Symbol Schemes
 
Consider the vector  in (<ref>). Ideally in single-symbol schemes only one of the entries in  should be a non-zero value, since only one element receives the symbol. However, this will not be the case in practice, since the beamformer will introduce a non-zero power also in directions other than the required one. Hence, practically more than one element of IRS_2 receives the symbol. However, the power in the undesired beam-indices is much lower than that in the intended index and these powers depend on the beampattern of IRS_1.  The data-dependent beam-index changes can be represented approximately by a beampattern rotation. Hence, (<ref>) can be written as:








  
    = Π_ + ,

  where  represents the vector of powers in the various indices of  and Π_ represents a particular permutation of the power pattern. Therefore, in order to identify the beam-index, we have to identify the power pattern permutation Π_. Now, ∼(Π_,Σ). Hence, the ML detector of this problem is formulated as:
  
    Π_min    ( - Π_)^HΣ^-1( - Π_) 
    Π_max    ℛe{( -1/2Π_)^HΣ^-1Π_}.

   In general, the search problem (<ref>) is NP-hard. However, in our case, there are only N_2M different patterns corresponding to N_2 different beam indices and M QAM/PSK symbols. Hence, a moderate-complexity search will give the optimal solution to the ML problem (<ref>).  
   

  §.§.§ Multi-Symbol Scheme

   The ML detector (<ref>) is also suitable for  multi-symbol case. However, in this case, since there are N_T desired beam-indices at a time, which interact with each other and thereby produce a large number of possible combinations Π_. Explicitly, N_2N_TM different patterns hypothesis must be tested for N_T RSs in IRS_1. Hence the ML detector may no longer be a computationally  attractable solution. Hence, in Section <ref>, we will be proposing a suboptimal compressed sensing (CS) aided detector, which can be used for any of the proposed schemes at a lower computational complexity.
 

 §.§ Suboptimal Compressed Sensing Detector

 
 The transmitted vector  in (<ref>) is sparse, when the number of active elements (i.e., elements that receive the symbol) is much less than the total number of elements in IRS_2. Therefore, one can use an efficient sparse reconstruction algorithm <cit.> for identifying the non-zero components in , which can be used to estimate $̱. However, it should be noted that for the successful recovery of the sparse vector, there should be a sufficient number of measurements. This can be either achieved by having a sufficient number of RAs (N_Rshould be sufficiently large) or taking multiple measurements, which would naturally reduce the data rate. Finally,scan be obtained from the estimated$̱ as:
 
    ŝ = s ∈ℳmin  - s̱^2,

 where ℳ is the constellation used.


 §.§ Complexity

The optimal ML detector has to compute
  (<ref>) for all possible combinations, which requires
  approximately on the order of (N_R^3 + N_RN_2)
  multiplications. This has to be done for each possible symbol. For
  the multi-symbol case, there are N_2N_TM possible
  symbols. Hence, the total computational complexity is approximately
  on the order of N_2N_TM(N_R^3 + N_RN_2), which reduces
  to the order of N_2M(N_R^3 + N_RN_2) for single-symbol
  cases. On the other hand, if a greedy compressed sensing based
  suboptimal algorithm is used, the complexity will be reduced to the
  order of N_2N_RN_T, which is much lower than that of the optimal
  ML detector.


§ AVERAGE BIT ERROR RATE ANALYSIS


In this section, we will estimate an upper bound for the average bit error rate (BER) of the optimal ML detector of Section <ref>. Let {Π_^i →Π_^j } represent the probability that the pattern Π_^i is identified as Π_^j and ν_i,j represent the number of bits in error between the two permutations Π_^i and Π_^j. Then the average BER is formulated as:

    B̂ÊR̂ = ∑_i=1^Ω∑_j  ij =1^Ων_i,j/n_bΩ{Π_^i →Π_^j },

where n_b is the total number of bits per channel use and Ω is the total number of possible permutations. Equation (<ref>) assumes that all permutations are equally likely. The probability of symbol error {Π_^i →Π_^j } in (<ref>) can be found as follows. When Π_^i is transmitted, the detector identifies Π_^j as the transmitted symbol  based on:

    kmax ℛe{( - 1/2Π_^k)^HΣ^-1Π_^k } = j,

where  is given in (<ref>) in conjunction with Π_= Π_^i. Let us define r_k = ℛe{( - 1/2Π_^k)^HΣ^-1Π_^k }. Hence, we have

    {Π_^i →Π_^j } = {k  j⋂(r_j > r_k) }.

The computation of the probability of intersection of the event in (<ref>) is very difficult. Hence, it is bounded using Fretchet's inequality <cit.> as follows:

    {k  j⋂(r_j > r_k) }≤kmin { r_j > r_k }.

In order to estimate the bound, the probabilities of { r_j > r_k } have to be calculated for each k  j. Theorem <ref> stated below gives an expression of the probability { r_j > r_k }.

	
Let us assume the Rician channel
          model with {}̋= and that the row vectors of
          $̋ are independent and identically distributed with
          covariance matrixΣ̃_cand also letσ_2^2 ≪σ_R^2. WhenΠ_^iis
        the actual signal transmitted andΘis a constant
        matrix independent of$̋, the probability of the events r_j
        > r_k, i.e.  { r_j > r_k } can be
        approximated as:
	
		
  * For k = i:
		
    { r_j > r_i  }  =1/2[1 - √(β_1/1+β_1)∑_n=0^N_R-12nn( 1/4(1+β_1))^n],


		where β_1 = ( Π_^i -Π_^j )^H Θ^H(Σ̃_c + 1/N_R^H)Θ( Π_^i -Π_^j )/4σ_R^2.
		
  * For k  i and when q_R = ℛe{q} 0, where q is defined in (<ref>):

		
    { r_j > r_k  } =   1/2[1 - 1/2√(β_2/1+β_2)∑_n=0^N_R-12nn( 1/4(1+β_2))^n ],


		where β_2 = q_R^2/(2+σ_κ^2)σ_z_1^2  and the constants σ_z_1^2 and σ_κ^2 are defined in  (<ref>) and (<ref>), respectively. 
				
  * For k  i and when q_R =0, { r_j > r_k  } = 1/2.
			
	
		See Appendix <ref> for proof.
	
 

Finally, for each transmitted symbol Π_^i, the minimum value of { r_j > r_k }, ∀ k  j is computed using Theorem <ref> and it is substituted for {Π_^i →Π_^j } into (<ref>) for achieving the bound of the average BER. Note that (<ref>) and (<ref>) also hold for the Rayleigh channel, in which case the matrix (Σ̃_c + 1/N_R^H) is replaced by  to calculate β_1 and β_2.
  The bounds derived for the average BER can be used for both Scheme 1 and Scheme 2. In the case of Scheme 1 Ω = N_2M and n_b = log_2 M + ⌊log_2 N_2 ⌋, whereas for Scheme 2, the corresponding values are Ω = N_2N_TM and n_b = log_2 M + ⌊log_2 N_2N_T⌋. For Scheme 3, Θ is no longer independent of $̋ and therefore the bounds in Theorem <ref> do not hold. However, the conditional probabilities derived in  Appendix <ref> can be used for Scheme 3 also. Based on this the unconditional probabilities can be derived using sampling method for computing the bound.



§ SIMULATION RESULTS

Extensive simulations have been carried out to establish the performance of the proposed scheme on the system parameters. Explicitly, we studied the average BER of the proposed schemes vs. the SNR, the number of elements (N_2) inIRS_2, the number of receivers (N_R), the number of RSs inIRS_1and the Rician factor, denoted byK. We have considered both the optimal ML detector and the low complexity compressed sensing detector in our performance evaluation. The system parameters used are given below. 

  * IRS_1 : A 100 × 100 rectangular array with spacing 2.5 mm. This corresponds to half wavelength of the frequency. 
	
  *  Distance between IRSs (D): 30 m. This distance satisfies the far-field condition for IRS_1. 
	
  * IRS_2: In general, an 8 × 8 rectangular array is used with inter element spacing of d_2=60 cm. However, these parameters are changed for the various performance studies, which is mentioned in the corresponding discussions. 

We used16levelQAMat60 GHzin all simulations. The results are shown in Fig. <ref>-<ref> for10^5Monte Carlo runs. We used both Rician and Rayleigh fading channels. Throughout the simulations, it is assumed that the channel is perfectly known at the receiver. Channel estimation in IRS-aided systems is quite a challenge, since the IRS is passive and has no signal processing capability. However, the schemes adopted in <cit.> could be used as a solution to the channel estimation problem.
In Fig. <ref>-<ref>, the proposed schemes 1 and 2 are referred  asS1andS2, respectively.S3represents the solution of (<ref>), whereasS301andS302represent the solution of (<ref>) and (<ref>), respectively.S1-Err represents the results of the proposed scheme 1 in the presence of channel estimation errors. ML and CS represent the results of optimal ML and compressed sensing based detectors, respectively and UB, the theoretical upper bound given in (<ref>).


Fig. <ref> shows the performance of
  schemeS1at various SNRs for different values of the Rician factor
  (K). Observe from Fig. <ref> that the performance of
  the low-complexity CS detector is inferior to that of the ML
  detector. As for the ML detector, the average BER tends to zero above10 dBSNR, whereas it exhibits an error floor near0.1for
  the CS detector. The performance of the CS detector is improved,
  when the Rician factorKdecreases. This is because asKdecreases, the projection matrixof the signal recovery
  becomes more random and the restricted isometric property
  <cit.> is improved and it is best for the Rayleigh
  channel. However, the performance of the ML detector is improved
  withK, since the LOS component is increased withK. This is
  also reflected by the upper bound seen in Fig.  <ref>. 

Fig. <ref>  shows the performance of single-symbol beam-index modulation (S1) against the Rician factor (K) for different SNRs. The  average BER obtained through simulations is compared against the upper bound derived in Section <ref>. The BER improvement vs.Kis due to increasing the LOS component. However, the variation in BER is only moderate, because the SNR is kept constant upon increasingK. Furthermore, the optimal ML detector has perfect channel knowledge. However, asKincreases, the gap between the simulation results and the corresponding upper bound is reduced.  This is also observed at high SNRs. Therefore, it can be concluded that for both these cases, our bound becomes tighter.

        
The performance of the CS detector can be improved by increasing the number of RAs as shown in Fig. <ref>, where the average BER is plotted againstN_R. This is plotted for Rayleigh channel, which gives the best performance for CS detector. The curves are shown for different number of elements inIRS_2(i.e.N_2). Observe that for the same number of RAs, the performance degrades, asN_2increases.  However, asN_2increases, the data rate will increase. 
                     

        Fig.  <ref> compares the ML detector's performance for the three proposed schemes in terms of their average BER for Rician factorK=0 dB. ForS2, we usedN_T=2, i.e. the number of RSs inIRS_1is two. This is because, ifN_Tis large, the complexity of optimal ML decoding will escalate. For fair comparison, the number of elements in all three cases are kept the same. Therefore, forS3, where the optimization is to be carried out in an array, a2 ×2element array is considered to form a single RS. Hence, the effective dimension ofIRS_2inS3is4 ×4, while it is8 ×8in the case ofS1andS2.  Hence, the data rate will be lowest forS3, whereas it is the highest forS2, since there are more RSs inIRS_1. The data rate forS1,S2andS3are 10, 14 and 8bpcu, respectively. In Fig. <ref>, the legendsS3O1andS3O2represent the results of two optimization methods, i.e. the solution of (<ref>) and that of (<ref>), respectively. Both these schemes perform better thanS1andS2. This is because there is an increase in the received SNR due to optimization. In addition in S3, the modulating symbol is embedded in4elements, which gives an additional performance improvement. This makes the BER gap between the curves ofS3and the other schemes substantial. Observe thatS3O1performs marginally better thanS3O2. This is becauseS302maximizes the lower bound, whereasS301operates on the exact function. Note thatS301performs almost similar to the solution of the exact equation (<ref>). The performance ofS2is approximately2 dBworse than that ofS1. There are two differences between these two schemes. InS2, the same QAM/PSK symbol is carried by more than one elements inIRS_2. Hence, the probability of error in decoding the modulating symbol is reduced compared toS1. However, the information carried by the beam-index is higher in the case ofS2, whose probability of decoding error will be higher than that ofS1. The average BER reflects these two opposite effects. 
                             
     Fig. <ref> shows the effect of the number of RSs (N_T) inIRS_1on the BER performance inS3for different number of receivers (N_R) for CS detector. Similar to Fig. <ref>, this is also under Rayleigh channel condition. The average  BER increases asN_Tincreases, which can be reduced by increasing the number of receivers. However, asN_Tincreases, the data rate increases. In this case, for the single RS case (which is equivalent toS1), the data rate is10 bpcu, while it is  14, 23, 30 and 36bpcuforN_T = 2, 4, 6  &  8, respectively. Finally, in Fig. <ref>, the effect of channel estimation errors is demonstrated. The true channel coefficients are corrupted by adding noise having a variance ofσ_R^2, which affects both the  optimization as well as detection. The average BER is shown in the figure both with and without channel estimation error. It can be seen that bothS1andS2have approximately2-3 dBperformance degradation owing to the channel estimation error, whereas this gap is in excess of4 dBforS3. This is because, inS3, the contaminated channel information is used both for optimization and detection.
                                    
   
                               

§ CONCLUSIONS
   
We proposed beam index modulation for millimeter wave communication exploiting the benefits of IRSs. The proposed scheme has three main advantages: 1) It achieves low-cost beamforming by using IRS for applying phase shifts, 2) it is capable of achieving reliable communication with the help of multiple IRSs in non-LOS scenarios, and 3) it sends additional information using beam-index modulation without any additional cost. Furthermore, we developed the optimal ML detector and a low-complexity compressed sensing detector for the proposed schemes. An upper bound of the average BER of the optimal ML detector is also achieved. Finally, the performance of the proposed schemes was evaluated through extensive simulations. 


§ PROOF OF LEMMA <REF>


    _̋QΘ_Q _̱Q s^2    =  (_̱Q s)^HΘ_Q^H_̋Q^H_̋QΘ_Q (_̱Q s) 
       = Tr{ (_̱Q s)^HΘ_Q^H_̋Q^H_̋QΘ_Q (_̱Q s)} = Tr{Θ_Q^H_̋Q^H_̋QΘ_Q (_̱Q s) (_̱Q s)^H}.

  Note that Θ_Q^H_̋Q^H_̋QΘ_Qand(_̱Q s) (_̱Q s)^Hare positive definite matrices. Therefore applying <cit.> on  (<ref>) results in (<ref>). Now, sinceH_Qis a random matrix and ifR ≥Q,_̋Q^H_̋Qwill be a full-rank matrix with probability 1 and consequentlyλ_min(Θ_Q^H_̋Q^H_̋QΘ_Q )  > 0. 



§ PROOF OF THEOREM <REF>

Lemma <ref> of Appendix <ref> gives the conditional probability{ r_j > r_k |}. Explicitly, the conditional probabilities are in the form of complementary error function (Q-functions). The distributions of the arguments of theseQ-functions are derived in Lemma <ref> in Appendix <ref>. Therefore, the unconditional probabilities can be calculated by taking expectation of conditional probabilities (<ref>) with respect to the corresponding distributions of their arguments in (<ref>) and (<ref>).  
 
Fork =i, the argument of the conditional probability is aΓ-distributed random variable with parametersΓ(N_R, 2σ_R^2/( Π_^i -Π_^j )^HΘ^HΣ_cΘ( Π_^i -Π_^j )). The unconditional probability in this case is 

    { r_j > r_i  }  = 	∫_0^∞Q(√(γ_ij))  (2σ_R^2/( Π_^i -Π_^j )^HΘ^HΣ_cΘ( Π_^i -Π_^j ))^N_Rγ_ij^N_R-1e^-2σ_R^2/( Π_^i -Π_^j )^HΘ^HΣ_cΘ( Π_^i -Π_^j )γ_ij/Γ(N_R)dγ_ij.
 
 The closed-form expression for (<ref>) given in <cit.> can be applied to get (<ref>). This proves the first part of the theorem. 
 
 Fork i, the probability is computed as follows. First we will consider the case ofq_R 0. By exploiting the relationshipQ(x) = 1- Q(-x), the probability{ r_j > r_k  } can be written as:


    { r_j > r_k  }   = 	∫_-∞^∞Q(κ)f_κ(κ)dκ 	= 	∫_-∞^0(1-Q(-κ))f_κ(κ)dκ + 	∫_0^∞Q(κ)f_κ(κ)dκ
       = 	∫_0^∞f_κ(-κ)dκ + 	∫_0^∞Q(κ)(f_κ(κ)-f_κ(-κ))dκ = C̃(I_1+I_2),


whereC̃ = 2Γ(2N_R)/Γ(N_R)√(πσ_κ^2)(v-1/2v)^N_Ris a constant term and 

    I_1    = ∫_0^∞ e^- 2v-1/2vσ_κ^2κ^2D_-2N_R(√(2/vσ_κ^2)κ)dκ,

while 

    I_2 = ∫_0^∞ Q(κ)e^- 2v-1/2vσ_κ^2κ^2(D_-2N_R(-√(2/vσ_κ^2)κ) - D_-2N_R(√(2/vσ_κ^2)κ))dκ.

NoteD_(.)(.)is the parabolic cylinder function <cit.> and it can be written in terms Kummer's confluent hypergeometric function_1F_1(.)as <cit.>
    D_K(z) = 2^K/2√(π)e^-z^2/4(1/Γ(1-K/2)_1 F_1(-K/2;1/2;z^2/2)-z/√(2)Γ(-K/2)_1 F_1(1-K/2;3/2;z^2/2)).

Let us substituteκ= +√(t)into (<ref>) and expandD_-2N_R(.)using (<ref>). Note thatdκ= 1/2√(t)dt. ThereforeI_1becomes:

    I_1    = √(π)2^-(N_R+1)/Γ(N_R + 1/2)∫_0^∞t^-1/2e^-t/σ_κ^2 _1F_1(N_R;1/2;t/vσ_κ^2)dt 
       - 2^-(N_R+1)/Γ(N_R)√(π/vσ_κ^2)∫_0^∞ e^-t/σ_κ^2 _1F_1(N_R+1/2;3/2;t/vσ_κ^2) dt.

Now, the difference in (<ref>) is formulated as: 

    Δ   = D_-2N_R(-√(2t/vσ_κ^2)) - D_-2N_R(√(2t/vσ_κ^2)) 
       =  2e^-t/2vσ_κ^22^-N_R/Γ(N_R)√(π t/vσ_κ^2) _1 F_1(N_R+1/2;3/2;t/vσ_κ^2).

	In order to evaluateI_2, first we express theQ-function in terms of the complimentary error function asQ(x) = 1/2erfc(x/√(2))<cit.>, and then subsequently it is expressed in terms of the hypergeometric function as <cit.>:
	
    Q(x) = 1/2 - x/√(2π) _1 F_1(1/2;3/2;-x^2/2).

	Upon substituting (<ref>) and (<ref>) into (<ref>),I_2becomes:
	
    I_2    = 2^-(N_R+1)/Γ(N_R)√(π/vσ_κ^2)∫_0^∞ e^-t/σ_κ^2 _1F_1(N_R+1/2;3/2;t/vσ_κ^2)dt 
       - 2^-N_R/Γ(N_R)√(1/2vσ_κ^2)∫_0^∞ t^1/2e^-t/σ_κ^2 _1 F_1(1/2;3/2;-t/2) _1F_1(N_R+1/2;3/2;t/vσ_κ^2) dt.

	Note that the second term of the RHS in (<ref>) and the first term of RHS in (<ref>) will get cancelled. Hence{ r_j > r_k  } will become:
	
    { r_j > r_k  }   = C_1 ∫_0^∞t^-1/2e^-t/σ_κ^2 _1F_1(N_R;1/2;t/vσ_κ^2)dt 
       - C_2 ∫_0^∞ t^1/2e^-t/σ_κ^2 _1 F_1(1/2;3/2;-t/2) _1F_1(N_R+1/2;3/2;t/vσ_κ^2) dt 
       = C_1I_3 -C_2I_4,

	whereC_1 = 1/2√(πσ_κ^2)(v-1/v)^N_RandC_2 = 1/πσ_κ^2√(2v)Γ(N_R+1/2)/Γ(N_R)(v-1/v)^N_R. Note that we have exploited the relationshipΓ(x)Γ(x +1/2) = √(π)/2^2x-1Γ(2x)<cit.> for reducingC_1andC_2. Now,I_3is expressed using equation <cit.> and it is given below:
	
    I_3    = ∫_0^∞t^-1/2e^-t/σ_κ^2 _1F_1(N_R;1/2;t/vσ_κ^2)dt 
       = Γ(1/2)√(σ_κ^2) _2F_1(N_R, 1/2;1/2;1/v) = √(πσ_κ^2)(1 - 1/v)^-N_R,

	where_2 F_1(a,b;c;z)is Gauss' Hypergeometric function <cit.>. Observe that_2F_1(a, b;b;z ) = _1F_0(a; z ) = (1- z)^-a<cit.>. In order to evaluate the integralI_4, we make the substitutionz = t/vσ_κ^2and use <cit.> and the integral becomes:
	
    I_4    = ∫_0^∞ t^1/2e^-t/σ_κ^2 _1 F_1(1/2;3/2;-t/2) _1F_1(N_R+1/2;3/2;t/vσ_κ^2) dt 
       = (vσ_κ^2)^3/2∫_0^∞ z^1/2e^-vz _1 F_1(1/2;3/2;-vσ_κ^2/2z) _1F_1(N_R+1/2;3/2;z) dz 
       = √(π/2(2+σ_κ^2))(σ_κ^2)^3/2(v-1/v)^-(N_R+1/2) _2F_1(N_R+1/2,1/2;3/2;-σ_κ^2/(v-1)(2+σ_κ^2)).

	Now we apply the transformations _2 F_1(a,b;c;z) = (1-z)^-b  _2 F_1(c-a,b;c;z/z-1)<cit.> and_2 F_1(a,b;b+1;z) = bz^-bℬ_z(b,1-a)<cit.>, when (<ref>) becomes:
	
    I_4 = σ_κ^2/2√(π v/2)(v-1/v)^-N_Rℬ_σ_κ^2/2(v-1)+vσ_κ^2(1/2,N_R),

whereℬ_z(a,b) = ℬ(a,b) z^a  ∑_k=0^b-1(a)_(k)/k! (1 -z)^kis the incomplete Beta function <cit.> with ℬ(a,b) = Γ(a)Γ(b)/Γ(a+b)being the Beta function and(a)_(k)is the Pochhammer symbol. Finally upon substituting (<ref>) and  (<ref>) into (<ref>), we arrive at:

    { r_j > r_k  }  = 1/2 - 1/4√(σ_κ^2/2(v-1)+vσ_κ^2)∑_n=0^N_R-1(1/2)_(n)/n!(1 - σ_κ^2/2(v-1)+vσ_κ^2)^n.

Now substituting forvand the Pochhammer symbol(a)_(n) = Γ(n+a)/Γ(a)<cit.> will give the second term in the RHS of (<ref>). 
 
Finally, whenq_R =0, the distribution ofκis zero mean Gaussian. Hence, in this caseI_2in (<ref>) will be zero, whileI_1 = 1/2and the constantC̃ =1, which completes the proof.



§ CONDITIONAL PROBABILITY {R_J > R_K |}

		Let Π_^i be the transmitted signal. Then, we have 
			
    { r_j > r_k |}=
    			
    			Q(√(γ_ij)),    if k = i 
    
    			Q(κ),    otherwise,

			where Q(x) = 1/√(2π)∫_x^∞^-t^2/2dx is the complementary error function and
			
    γ_kj =1/2( Π_^k -Π_^j )^H ^HΣ^-1( Π_^k - Π_^j ),

			and 
			
    κ = γ_ik - γ_ij/√(γ_kj).

					First the conditional probability { r_j > r_k |} is estimated for k  i when Π_^i is transmitted as follows.The event r_j > r_k is 			
					
    ℛe{( - 1/2Π_^k)^HΣ^-1Π_^k } < ℛe{( -1/2Π_^j)^HΣ^-1Π_^j }.

					Using (<ref>), (<ref>) can be written as:
					
    ℛe{( + Π_^i )^H Σ^-1( Π_^k -Π_^j ) } < ℛe{-1/2( Π_^j -Π_^k )^H ^HΣ^-1( Π_^j + Π_^k )}.

					Hence, { r_j > r_k |} = {η < g_jk | },




					where η = ℛe{( + Π_^i )^H Σ^-1( Π_^k - Π_^j )}  and g_jk is the RHS of  (<ref>). Finally, note that η∼(μ_jk,γ_jk), where μ_jk = ℛe {(Π_^i)^H ^HΣ^-1( Π_^k - Π_^j )} and γ_jk is defined in (<ref>). Now, let us make a substitution η̃ = η -μ_jk/√(γ_jk). Clearly, η̃ is a standard normal random variable and hence the probability in (<ref>) can be written in the form of the Q-function as:
					
    {r_j > r_k | } = {η̃  <  -s_jk/√(γ_jk)} = Q(s_jk/√(γ_jk)),

					where we have
						
    s_jk   = -(g_jk-μ_jk)=ℛe {1/2( Π_^j -Π_^k )^H ^HΣ^-1( Π_^j + Π_^k -2Π_^i )}
       = ℛe {1/2( (Π_^j - Π_^i )-(Π_^k -Π_^i ) )^H ^HΣ^-1((Π_^j - Π_^i )+(Π_^k -Π_^i ) )}
       = γ_ij - γ_ik.
 
					When k=i, we have γ_ik=0 and therefore (<ref>) becomes {r_j > r_i | } = Q(√(γ_ij)).




					This completes the proof. 
			

§ DISTRIBUTION OF Γ_MN AND Κ

	Under the assumptions of Theorem <ref>,
	
  *   The random variable γ_kj defined in (<ref>) is a Γ-distributed random variable, i.e., 
		
    γ_kj∼Γ(N_R, 2σ_R^2/( Π_^k -Π_^j )^HΘ^HΣ_cΘ( Π_^k -Π_^j )).

  *  When q_R  0, the distribution of the random variable κ defined in (<ref>) is 
			
    f_κ(κ) =2Γ(2N_R)e^-2v-1/2vσ_κ^2κ^2/Γ(N_R)√(πσ_κ^2)(v-1/2v)^N_RD_-2N_R(-√(2/vσ_κ^2)κ),

		where  D_K(.) is the Parabolic cylinder function <cit.>  and the parameters are defined in Theorem <ref>.  
		
  *  For q_R =0, κ∼(0,σ_κ^2/2), where σ_κ^2 is defined in (<ref>)
		 
	
			When we have σ_2^2 ≪σ_R^2, γ_kj in (<ref>) can be approximated as:
			
    γ_kj≈1/2σ_R^2( Π_^k -Π_^j )^H ^H( Π_^k - Π_^j ).

		Recall that = Θ̋. Under the assumptions of Theorem <ref>, ^̋H$̋ obeys a non-central complex Wishart distribution, i.e.,^̋H∼̋𝒲(N_R, ,Σ̃_c). This can be approximated as a central complex Wishart distribution having the covariance matrix ofΣ_c = Σ̃_c + 1/N_R^H<cit.>[For Rayleigh channel, ^̋H$̋ exactly follows the central complex Wishart distribution associated withΣ_c = . If we use this approximation, it can be assumed that  each row of$̋ is distributed according to (0,Σ_c). Hence, the random variable = 1/√(2σ_R^2)( Π_^i -Π_^j ) is distributed according to (0,1/2σ_R^2( Π_^i -Π_^j )^HΘ^HΣ_cΘ( Π_^i -Π_^j )).  Hence, γ = ^H is a Gamma distributed variable having the distribution function of (<ref>). This proves the first part of the Lemma.
			 
			The  rest of the Lemma is proved as follows. Define _1 = 1/√(2σ_R^2)( (Π_^j - Π_^i )-(Π_^k -Π_^i ) ) and _2 = 1/√(2σ_R^2)( (Π_^j - Π_^i )+(Π_^k -Π_^i ) ). Similar to the case of first part of this Lemma, the distributions of _1 and _2 can be approximated as _1 ∼(0,σ_z_1^2 ) and _2 ∼(0,σ_z_2^2), where 
			
    σ_z_1^2 =  (Π_^j - Π_^k )^HΘ^HΣ_cΘ(Π_^j - Π_^k )/2σ_R^2,

			and 
			
    σ_z_2^2 =(Π_^j + Π_^k -2Π_^i)^HΘ^HΣ_cΘ(Π_^j + Π_^k -2Π_^i)/2σ_R^2.

			Also note that {_1_2^H }  = q and {_2_1^H }  = q^H, where
			
    q = 1/2σ_R^2Tr{(Π_^j - Π_^k )( Π_^j + Π_^k  -2Π_^i )^H Θ^HΣ_cΘ} = ( Π_^j + Π_^k  -2Π_^i ) ^HΘ^HΣ_cΘ(Π_^j - Π_^k )/2σ_R^2.

			Now we can rewrite γ_kj in (<ref>) as:
			
    γ_kj≈1/2σ_R^2( (Π_^j - Π_^i)-(Π_^k - Π_^i ) )^H ^H ( (Π_^j - Π_^i)-(Π_^k - Π_^i ) ) = _1^2.

		Similarly κ≈ℛe{_1^H_2/_1}.
		First we derive the distribution of κ given _1. Note that _1 and _2 are complex Gaussian distributed random vectors. Hence, the distribution f_κ|_1(κ|_1) is Gaussian. Let κ̃ = _1^H_2/_1. Therefore, using <cit.>, we have μ_κ̃ = {κ̃ |_1 = } =q^H/σ_z_1^2.
		Note that κ = ℛe{κ̃} and hence {κ|_1 = } =  ℛe{μ_κ̃} = q_R/σ_z_1^2, where q_R=ℛe{q}. In order to compute variance of the κ given _1, let us expand κ as:
		
    κ =  ℛe{_1^H_2/_1} = ℛe{_1}^T	ℛe{_2}/_1 +  ℐm{_1}^Tℐm{_2}/_1 = u_1 + u_2.

		Hence,
		
    var(κ |_1 = ) = var(u_1|_1 = )+var(u_2|_1 = ) + Cov(u_1,u_2|_1 = )+Cov(u_2,u_1|_1 = ).

		Using <cit.>, it can be shown that: 
		
    var(u_1|_1 = ) =ℛe{}^TCov{ℛe{_2} |_1= }ℛe{}/^2=  ℛe{}^2/^21/2(σ_z_2^2  - q^2/2σ_z_1^2).

		Similarly, var(u_2|_1 = ) =  ℐm{}^2/^21/2(σ_z_2^2  - q^2/2σ_z_1^2) and 
		
    Cov(u_1,u_2|_1 = ) =	-Cov(u_2,u_1|_1 = ) =-jq^2/4σ_z_1^2ℛe{_1}^Tℐm{_1}/_1^2.

		Therefore, var(κ |_1 = ) = 1/2( σ_z_2^2 - q^2/2σ_z_1^2). Hence, κ|_1 ∼(q_R/σ_z_1^2,σ_κ^2/2), where
		
    σ_κ^2 = ( σ_z_2^2 - q^2/2σ_z_1^2).

		If q_R =0, f_κ|_1 is independent of _1 and hence the unconditional distribution of κ is the same as the conditional distribution, i.e., κ∼(0,σ_κ^2/2).
			 
			
			When q_R  0,  the unconditional distribution can be obtained by eliminating the conditioning with respect to the distribution of z̃, which is Nakagami distributed Nakagami(N_R, N_Rσ_z_1^2), since  z̃^2 ∼Γ(N_R,1/σ_z_1^2). Explicitly, f_Z̅(z̅) = 2z̅^2N_R-1e^-z̅^2/σ_z_1^2/Γ(N_R)(σ_z_1^2)^N_R<cit.>, where z̅ = z̃. Hence, the unconditional distribution of κ is:
			
    f_κ(κ) = 2∫_0^∞ e^-g(κ,z̅) z̅^2N_R-1dz̅/(σ_z_1^2)^N_RΓ(N_R)√((πσ_κ^2)) ,

			where 
			
    g(κ,z̅) = (κ -q_R/σ_z_1^2z̅)^2/σ_κ^2 +z̅^2/σ_z_1^2 = κ^2/σ_κ^2+ q_R^2+σ_z_1^2σ_κ^2/(σ_z_1^2)^2σ_κ^2z̅^2-2q_Rκ/σ_z_1^2σ_κ^2z̅.

			Now using <cit.>, we arrive at: 
	
			
    I_κ   = ∫_0^∞e^-g(κ,z̅) z̅^2N_R-1dz̅ =e^-κ^2/σ_κ^2∫_0^∞e^-q_R^2+σ_z_1^2σ_κ^2/(σ_z_1^2)^2σ_κ^2z̅^2+2q_Rκ/σ_z_1^2σ_κ^2z̅z̅^2N_R-1dz̅
       = e^-κ^2/σ_κ^2(σ_z_1^2(v-1)/2v)^N_RΓ(2N_R)e^κ^2/2σ_κ^2vD_-2N_R(-√(2/σ_κ^2v)κ),

		
	where v = 1+ σ_κ^2σ_z_1^2/q_R^2. Finally, substituting (<ref>) into (<ref>) will give (<ref>). 
	ieeetr]