Matthieu Roy

Stéphane Paquelet

Luc Le Magoarou

Matthieu Crussière

MIMO Channel Hardening: A Physical Model based Analysis

Keywords: channel hardening, physical model, MIMO

In a multiple-input-multiple-output (MIMO) communication system, the multipath fading is averaged over radio links. This well-known channel hardening phenomenon plays a central role in the design of massive MIMO systems. The aim of this paper is to study channel hardening using a physical channel model in which the influences of propagation rays and antenna array topologies are highlighted. A measure of channel hardening is derived through the coefficient of variation of the channel gain. Our analyses and closed form results based on the used physical model are consistent with those of the literature relying on more abstract Rayleigh fading models, but offer further insights on the relationship with channel characteristics.

I. INTRODUCTION

Over the last decades, multi-antenna techniques have been identified as key technologies to improve the throughput and reliability of future communication systems. They offer a potential massive improvement of spectral efficiency over classical SISO (single-input-single-output) systems proportionally to the number of involved antennas. This promising gain has been quantified in terms of capacity in the seminal work of Telatar [START_REF] Telatar | Capacity of Multi-antenna Gaussian Channels[END_REF] and has recently been even more emphasized with the newly introduced massive MIMO paradigm [START_REF] Marzetta | Fundamentals of massive MIMO[END_REF].

Moving from SISO to MIMO, the reliability of communication systems improves tremendously. On the one hand in SISO, the signal is emitted from one single antenna and captured at the receive antenna as a sum of constructive or destructive echoes. This results in fading effects leading to a potentially very unstable signal to noise ratio (SNR) depending on the richness of the scattering environment. On the other hand in a MIMO system, with appropriate precoding, small-scale multipath fading is averaged over the multiple transmit and receive antennas. This yields a strong reduction of the received power fluctuations, hence the channel gain becomes locally deterministic essentially driven by its large-scale properties. This effect, sometimes referred to as channel hardening [START_REF] Hochwald | Multipleantenna channel hardening and its implications for rate feedback and scheduling[END_REF] has recently been given a formal definition based on the channel power fluctuations [START_REF] Ngo | No downlink pilots are needed in TDD massive MIMO[END_REF]. Indeed, studies on the stability of the SNR are essential to the practical design of MIMO systems, in particular on scheduling, rate feedback, channel coding and modulation dimensioning [START_REF] Marzetta | Fundamentals of massive MIMO[END_REF], [START_REF] Hochwald | Multipleantenna channel hardening and its implications for rate feedback and scheduling[END_REF], [START_REF] Björnson | Random access protocol for massive MIMO: Strongest-user collision resolution (SUCR)[END_REF]. From the definition in [START_REF] Ngo | No downlink pilots are needed in TDD massive MIMO[END_REF], we propose in this paper a comprehensive study on channel hardening through a statistical analysis of received power variations derived from the propagation characteristics of a generic ray-based spatial channel model. Related work. Channel hardening, measured as the channel gain variance, has recently been studied from several points of view. The authors in [START_REF] Martínez | Massive MIMO properties based on measured channels: Channel hardening, user decorrelation and channel sparsity[END_REF] used data from measurement campaigns and extracted the variance of the received power. A rigorous definition of channel hardening was then given in the seminal work [START_REF] Ngo | No downlink pilots are needed in TDD massive MIMO[END_REF] based on the asymptotic behavior of the channel gain for large antenna arrays. This definition was applied to pinhole channels, i.i.d. correlated and uncorrelated Rayleigh fading models [START_REF] Björnson | Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency[END_REF]. Contributions. Complementary to this pioneer work, we propose a non-asymptotic analysis of channel hardening, as well as new derivations of the coefficient of variation of the channel not limited to classically assumed Rayleigh fading models. Indeed, channel hardening is analyzed herein using a physically motivated ray-based channel model widely used in wave propagation. Our approach is consistent with previous studies [START_REF] Ngo | No downlink pilots are needed in TDD massive MIMO[END_REF], [START_REF] Martínez | Massive MIMO properties based on measured channels: Channel hardening, user decorrelation and channel sparsity[END_REF], but gives deeper insights on channel hardening. In particular, we managed to provide an expression of the channel hardening measure in which the contributions of the transmit and receive antenna arrays, and the propagation conditions can easily be identified, and thus interpreted. Notations. Upper case and lower case bold symbols are used for matrices and vectors. z * denotes the conjugate of z. u stands for a three-dimensional (3D) vector. .,. and a• u denote the inner product between two vectors of C N and 3D vectors, respectively. [H] p,q is the element of matrix H at row p and column q. H F , h and h p stand for the Frobenius norm, the euclidean norm and the p-norm, respectively. H H and H T denotes the conjugate transpose and the transpose matrices. H denote the normalized matrix H/||H|| F . E{.} and Var{.} denote the expectation and variance.

II. CHANNEL MODEL

We consider a narrowband MIMO system (interpretable as an OFDM subcarrier) with N t antennas at the transmitter and N r antennas at the receiver, such that y = Hx+n, with x ∈ C Nt×1 , y ∈ C Nr×1 and n ∈ C Nr×1 the vectors of transmit, receive and noise samples, respectively. H ∈ C Nr×Nt is the MIMO channel matrix, whose entries [H] i,j are the complex gains of the SISO links between transmit antenna j and receive antenna i. The capacity of the MIMO channel can be expressed as [START_REF] Telatar | Capacity of Multi-antenna Gaussian Channels[END_REF]

C = log 2 (det(I Nt +ρ Q HH H)) bps/Hz, (1)
where ρ = Pt N0 ||H|| 2 F with Q ∈ C Nt×Nt , P t and N 0 the input correlation matrix (precoding), emitted power and noise power. C is a monotonic function of the optimal received SNR ρ [START_REF] Loyka | On physically-based normalization of MIMO channel matrices[END_REF], hence H 2 F directly influences the capacity of the MIMO channel. It is then of high interest studying the spatial channel gain variations to predict the stability of the capacity.

In the sequel, we will consider that the channel matrix H is obtained from the following generic multi-path 3D ray-based model considering planar wavefronts [START_REF] Raghavan | Sublinear Capacity Scaling Laws for Sparse MIMO Channels[END_REF], [START_REF] Zwick | A stochastic multipath channel model including path directions for indoor environments[END_REF], [7, p. 485]

H(f ) = N t N r P p=1
c p e r ( u rx,p )e t ( u tx,p ) H .

(

Such channel consists of a sum of P physical paths where c p is the complex gain of path p and u tx,p (resp. u rx,p ) its direction of departure -DoD -(resp. of arrival -DoA -). In (2) e t and e r are the so-called steering vectors associated to the transmit and receive arrays. They contain the path differences of the plane wave from one antenna to another and are defined as [START_REF] Zwick | A stochastic multipath channel model including path directions for indoor environments[END_REF] e

t ( u tx,p ) = 1 √ N t e 2jπ a tx,1 • u tx,p λ ,•••,e 2jπ a tx,N t • u tx,p λ T , (3)
and similarly for e r ( u rx,p ). The steering vectors depend not only on the DoD/DoA of the impinging rays, but also on the topology of the antenna arrays. The latter are defined by the sets of vectors A tx = { a tx,j } and A rx = { a rx,j } representing the positions of the antenna elements in each array given an arbitrary reference. Such channel model has already been widely used (especially in its 2D version) [START_REF] Raghavan | Sublinear Capacity Scaling Laws for Sparse MIMO Channels[END_REF], [START_REF] Zwick | A stochastic multipath channel model including path directions for indoor environments[END_REF], verified through measurements [START_REF] Samimi | 28 GHz Millimeter-Wave Ultrawideband Small-Scale Fading Models in Wireless Channels[END_REF] for millimeter waves and studied in the context of channel estimation [START_REF] Magoarou | Parametric channel estimation for massive MIMO[END_REF]. In contrast to Rayleigh channels, it explicitly takes into account the propagation conditions and the topology of the antenna arrays.

In the perspective of the following sections,

let c = [|c 1 |,•••,|c P |]
T denote the vector consisting of the amplitudes of the rays. c 2 is the aggregated power from all rays, corresponding to large-scale fading due to path-loss and shadowing.

III. CHANNEL HARDENING

Definition. Due to the multipath behavior of propagation channels, classical SISO systems suffer from a strong fast fading phenomenon at the scale of the wavelength resulting in strong capacity fluctuations [START_REF] Telatar | Capacity of Multi-antenna Gaussian Channels[END_REF]. MIMO systems average the fading phenomenon over the antennas so that the channel gain varies much more slowly. This effect is called channel hardening. In this paper, the relative variation of the channel gain H 2 F , called coefficient of variation (CV ) is evaluated to quantify the channel hardening effect as previously introduced in [4], [START_REF] Björnson | Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency[END_REF]:

CV 2 = Var H 2 F E{ H 2 F } 2 = E H 4 F -E H 2 F 2 E{ H 2 F } 2 (4)
In ( 4) the statistical means are obtained upon the model which govern the entries of H 2 given random positions of the transmitter and the receiver. This measure was previously applied to a N t ×1 correlated Rayleigh channel model h ∼ CN (0,R) [START_REF] Ngo | No downlink pilots are needed in TDD massive MIMO[END_REF], [7, p. 231]. In that particular case, (4) becomes

CV 2 = E |h H h| 2 -Tr(R) 2 Tr(R) 2 = Tr(R 2 ) Tr(R) 2 , (5)
where the rightmost equality comes from the properties of Gaussian vectors [START_REF] Björnson | Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency[END_REF]Lemma B.14]. This result only depends on the covariance matrix R, from which the influences of antenna array topology and propagation conditions are not explicitly identified. Moreover, small-scale and large-scale phenomena are not easily separated either. In this paper, (4) is studied using a physical channel model that leads to much more interpretable results.

Assumptions on the channel model. The multipath channel model described in Section II relies on several parameters governed by some statistical laws. Our aim is to provide an analytical analysis of CV while relying on the weakest possible set of assumptions on the channel model. Hence, we will consider that:

• For each ray, gain, DoD and DoA are independent.

• arg(c p ) ∼ U[0,2π] i.i.d.

• u tx,p and u rx,p are i.i.d. with distributions D tx and D rx . The first hypothesis is widely used and simply says that no formal relation exists between the gain and the DoD/DoA of each ray. The second one raisonnably indicates that each propagated path experiences independent phase rotation without any predominant angle. The last one assumes that all the rays come from independent directions, with the same distribution (distributions D tx at the emitter, D rx at the receiver).

It has indeed been observed through several measurement campaigns that rays can be grouped into clusters [START_REF] Saleh | A statistical model for indoor multipath propagation[END_REF], [START_REF] Wu | 60-GHz Millimeter-Wave Channel Measurements and Modeling for Indoor Office Environments[END_REF]. Considering the limited angular resolution of finite-size antenna arrays, it is possible to approximate all rays of the same cluster as a unique ray without harming a lot the channel description accuracy [START_REF] Magoarou | Parametric channel estimation for massive MIMO[END_REF]. It then makes sense to assume that this last hypothesis is valid for the main DoDs and DoAs of the clusters. Simulations. A preliminary assessment of the coefficient of variation is computed through Monte-Carlo simulations of (4) using uniform linear arrays (ULA) with inter-antenna spacing of λ 2 at both the transmitter and receiver and taking a growing number of antennas. A total of P ∈ {2, 4, 5, 6} paths were randomly generated with Complex Gaussian gains c p ∼ CN (0,1), uniform DoDs u tx,p ∼ U S2 and DoAs u rx,p ∼ U S2 .

Simulation results of CV are reported in Fig. 1 as a function of the number of antennas. It is observed that all curves seem to reach an asymptote around 1/P for large N t and N r . Hence, the higher the number of physical paths, the harder the channel. The goal of the next sections is to provide further interpretation of such phenomenon by means of analytical derivations.

IV. DERIVATION OF CV 2

In this section, CV 2 is analytically analyzed from (4).

Expectation of the channel gain. From ( 2) and (3) the channel gain H 2 F = Tr(H H H) can be written as

H 2 F = N t N r p,p c * p c p γ p,p ,
where the term γ p,p is given by γ p,p = e r ( u rx,p ),e r ( u rx,p ) e t ( u tx,p ),e t ( u tx,p ) * . Using the hypothesis arg(c p ) ∼ U[0,2π] i.i.d. introduced in the channel model and γ p,p = 1, the expectation of the channel gain can further be expressed as

E H 2 F = N t N r E c 2 . ( 6
) Thus the average channel gain increases linearly with N r and N t , which is consistent with the expected beamforming gain N t and the fact that the received power linearly depends on N r . Coefficient of variation. The coefficient of variation CV is derived using the previous hypotheses and ( 6). We introduce:

E 2 (A tx ,D tx ) = E | e t ( u tx,p ),e t * ( u tx,p ) | 2 E 2 (A rx ,D rx ) = E | e r ( u rx,p ),e r ( u rx,p ) | 2 . ( 7
)
These quantities are the second moments of the inner products of the transmit/receive steering vectors associated to two distinct rays. They represent the correlation between two rays as observed by the system. They can also be interpreted as the average inability of the antenna arrays to discriminate two rays given a specific topology and ray distribution. From such definitions, and based on the derivations given in Appendix A, CV 2 can be expressed as a sum of two terms,

CV 2 =E 2 (A tx ,D tx )E 2 (A rx ,D rx ) E c 4 -c 4 4 E{ c 2 } 2 + Var c 2 E{ c 2 } 2 . (8)
Note that this result only relies on the assumptions introduced in section II. The second term can be identified as the contribution of the spatial large-scale phenomena since it simply consists in the coefficient of variation of the previously defined large-scale fading parameter c 2 of the channel. To allow local channel behavior interpretation, conditioning the statistical model by c 2 is required. It results in the cancellation of the large-scale variations contribution of CV 2 which reduces to what is called hereafter small-scale fading.

V. INTERPRETATIONS A. Large-scale fading

The contribution of large-scale fading in CV 2 is basically the coefficient of variation of the total aggregated power c 2 of the rays. To better emphasize its let us consider a simple example with independent |c p | 2 of mean µ and variance σ 2 . The resulting large scale fading term is then

Var c 2 E{ c 2 } 2 = 1 P σ µ 2 .
It clearly appears that more rays lead to reduced large-scale variations. This stems from the fact that any shadowing phenomenon is well averaged over P independent rays, hence becoming almost deterministic in rich scattering environments. This result explains the floor levels obtained for various P in our previous simulations in Section II and is consistent with the literature on correlated Rayleigh fading channels where high rank correlation matrices provide a stronger channel hardening effect than low rank ones [START_REF] Björnson | Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency[END_REF].

B. Small-scale fading

The coefficient of variation particularized with the statistical conditional model can easily be proven to be:

CV 2 c 2 = E 2 (A tx ,D tx )E 2 (A rx ,D rx )α 2 (c)
where

α 2 (c) = 1- E c| c 2 c 4 4 c 4 . (9)
The small-scale fading contribution to CV 2 thus consists of a product of the quantities defined in ( 7) that depend only on the antenna array topologies (A tx /A rx ) and ray distributions (D tx /D rx ) multiplied by a propagation conditions factor α 2 (c) that depends only on the statistics of the ray powers c. Ray correlations. This paragraph focuses on the quantity E 2 (A tx , D tx ) (the study is done only at the emitter, the obtained results being equally valid at the receiver). Eq. ( 7) yields

E 2 (A tx , D tx ) = 1 N 2 t E    Nt i=1 e 2jπ a tx,i •( u tx,p -u tx,p ) λ 2    .
A well-known situation is when the inner sum involves exponentials of independent uniformly distributed phases and hence corresponds to a random walk with N t steps of unit length. The above expectation then consists in the second moment of a Rayleigh distribution and E 2 (A tx , D tx ) = 1 Nt . A necessary condition to such a case is to have (at least) a half wavelength antenna spacing ∆d to ensure that phases are spread over [0,2π]. On the other hand, phase independences are expected to occur for asymptotically large ∆d. It is however shown hereafter that such assumption turns out to be valid for much more raisonnable value of ∆d.

Numerical evaluations of E 2 are performed versus ∆d (Fig. 2), and versus N t (Fig. 3). Uniformly distributed rays over the 3D unit sphere (D tx = D rx = U S2 ) and Uniform Linear, Circular and Planar Arrays (ULA, UCA and UPA) are considered. As a reminder, the smaller E(A tx ,D tx ) the better the channel hardening. In Fig. 2, E 2 reaches the asymptote 1/N t for all array types with ∆d = λ 2 and remains almost constant for larger ∆d. Fig. 3 shows that E 2 merely follows the 1/N t law whatever the array type. We thus conclude that the independent uniform phases situation discussed above is a sufficient model for any array topology given that ∆d ≥ λ 2 . It is therefore assumed in the sequel that, (11) The right inequality comes from the convexity of the square function. Equality is achieved when there is only one contributing ray, i.e. no multipath occurs. In that case CV 2 c 2 = 0 and the MIMO channel power is deterministic. The left part in [START_REF] Samimi | 28 GHz Millimeter-Wave Ultrawideband Small-Scale Fading Models in Wireless Channels[END_REF] is given by Hölder's inequality. Equality is achieved when there are P rays of equal power. Then, taking the expectation on each member in [START_REF] Samimi | 28 GHz Millimeter-Wave Ultrawideband Small-Scale Fading Models in Wireless Channels[END_REF] yields [START_REF] Zwick | A stochastic multipath channel model including path directions for indoor environments[END_REF].

E 2 (A tx ,U S2 )≈1/N t , E 2 (A rx ,U S2 )≈1/N r .
In contrast to the large-scale fading, more rays lead to more small-scale fluctuations. It is indeed well known that a richer scattering environment increases small-scale fading.

Comparison with the simulations. Based on the general formula given in [START_REF] Loyka | On physically-based normalization of MIMO channel matrices[END_REF], on the interpretations and evaluations of its terms, we can derive the expression of channel hardening for the illustrating simulations of Section II:

CV 2 illustration = 1 N t N r
(1-1/P )+1/P.

Simulation and approximated formula are compared in Fig. 4 in which small-scale and large-scale contributions are easily evidenced, as intuitively expected from simulations of Fig. 1.

Comparison with the Gaussian i.i.d. model. This model assumes a rich scattering environment. Using (5) with R = I:

CV 2 iid = 1 N t N r .
Using the realistic model in a rich scattering environment, the large-scale part of (8) vanishes leading to a deterministic c 2 and small-scale variations reach the upper bound of (10). This yields the limit

CV 2 P →∞ ----→ CV 2 iid ( 12
)
which is coherent with the interpretation of the model. VI. CONCLUSION In this paper, previous studies on channel hardening have been extended using a physics-based model. We have separated influences of antenna array topologies and propagation characteristics on the channel hardening phenomenon. Large-scale and small-scale contributions to channel variations have been evidenced. Essentially, this paper provides a general framework to study channel hardening using accurate propagation models.

To illustrate the overall behavior of channel hardening, this framework have been used with generic model parameters and hypotheses. The scaling laws evidenced for simpler channel models are conserved provided the antennas are spaced by at least half a wavelength. The results are consistent with state of the art and provide further insights on the influence of array topology and propagation on channel hardening. The proposed expression can easily be exploited with various propagation environments and array topologies to provide a more precise understanding of the phenomenon compared to classical channel descriptions based on Rayleigh fading models. ACKNOWLEDGMENT

This work has been performed in the framework of the Horizon 2020 project ONE5G (ICT-760809) receiving funds from the European Union. The authors would like to acknowledge the contributions of their colleagues in the project, although the views expressed in this contribution are those of the authors and do not necessarily represent the project.

APPENDIX A COEFFICIENT OF VARIATION [START_REF] Loyka | On physically-based normalization of MIMO channel matrices[END_REF] For the sake of simplicity, an intermediary matrix A is introduced. It is defined by

Fig. 1 .

 1 Fig. 1. Simulated CV 2 for growing number of rays. Asymptotes are the black dashed lines.

Fig. 2 .

 2 Fig. 2. Numerical evaluation of E(Atx, U S 2 ) for various array types and increasing antenna spacing ∆d. The values are normalized so the asymptote is 1.

Fig. 3 .

 3 Fig. 3. Numerical evaluation of E(Atx,U S 2 ) for various antenna arrays at the half wavelength. The lower, the better.

Fig. 4 .

 4 Fig. 4. Comparison between (8) and simulated CV 2 . Uniform distribution of DoDs and DoAs over the unit sphere and complex Gaussian gains.Propagation conditions. It is now interesting to point out that the propagation factor α(c) introduced in (9) is bounded by 0 ≤ α 2 (c) ≤ 1-1/P. (10) Those bounds are deduced from the following inequality:c 4 2 /P ≤ c 4 4 ≤ c 4 2 .(11)The right inequality comes from the convexity of the square function. Equality is achieved when there is only one contributing ray, i.e. no multipath occurs. In that case CV 2 c 2 = 0 and the MIMO channel power is deterministic. The left part in[START_REF] Samimi | 28 GHz Millimeter-Wave Ultrawideband Small-Scale Fading Models in Wireless Channels[END_REF] is given by Hölder's inequality. Equality is achieved when there are P rays of equal power. Then, taking the expectation on each member in (11) yields[START_REF] Zwick | A stochastic multipath channel model including path directions for indoor environments[END_REF].In contrast to the large-scale fading, more rays lead to more small-scale fluctuations. It is indeed well known that a richer scattering environment increases small-scale fading.Comparison with the simulations. Based on the general formula given in[START_REF] Loyka | On physically-based normalization of MIMO channel matrices[END_REF], on the interpretations and evaluations of its terms, we can derive the expression of channel hardening for the illustrating simulations of Section II:

2 F 4 F

 24 [A] p,p = 2|γ p,p |cos(φ p,p ) if p = p 1 if p = p with φ p,p = arg(c * p c p γ p,p ) the whole channel phase dependence. H 2 F can be written using a quadratic form with vector c and matrix A, which can be decomposed into two terms I (identity) and J HN t N r = c T Ac = c T c+c T Jc where J = A-I. E{J} = 0 so: E H (N t N r ) 2 = E c 4 +E (c T Jc) 2 .The ray independence properties yields the following weighted sum of coupled ray powersE (c T Jc) 2 = p =p E |c p | 2 |c p | 2 E [J]2 p,p . Considering i.i.d. rays, all the weights E [J] 2 p,p are identical. Using the weights notations introduced in (7) and the definition of the 4-norm yields the second order moment E H 4 F . With the expectation (6) we derive the result (8).