Random Matrix Models, Double-Time Painlev\'e Equations, and Wireless Relaying

This paper gives an in-depth study of a multiple-antenna wireless communication scenario in which a weak signal received at an intermediate relay station is amplified and then forwarded to the final destination. The key quantity determining system performance is the statistical properties of the signal-to-noise ratio (SNR) \gamma\ at the destination. Under certain assumptions on the encoding structure, recent work has characterized the SNR distribution through its moment generating function, in terms of a certain Hankel determinant generated via a deformed Laguerre weight. Here, we employ two different methods to describe the Hankel determinant. First, we make use of ladder operators satisfied by orthogonal polynomials to give an exact characterization in terms of a"double-time"Painlev\'e differential equation, which reduces to Painlev\'e V under certain limits. Second, we employ Dyson's Coulomb Fluid method to derive a closed form approximation for the Hankel determinant. The two characterizations are used to derive closed-form expressions for the cumulants of \gamma, and to compute performance quantities of engineering interest.


I. INTRODUCTION
Over the past decade, multiple-input multiple-output (MIMO) antenna systems employing space-time coding have revolutionized the wireless industry. Such systems are well-known to offer substantial benefits in terms of channel capacity, as well as improved diversity and link reliability, and as such these techniques are being incorporated into a range of emerging industry standards. More recently, relaying strategies have been also proposed as a means of further improving the performance of space-time coded MIMO networks. Such methods are particularly effective for improving the data transmission quality of users who are located near the periphery of a communication cell.
Various MIMO relaying strategies have been proposed in the literature; see, e.g., Refs. 1-6, offering different trade-offs in various factors such as performance, complexity, channel estimation requirements, and feedback requirements. In Ref. 7, a communication technique was considered which employed a form of orthogonal space-time block coding (OSTBC), along with non-coherent amplify-and-forward (AF) processing at the multi-antenna relay. This method has the advantage of operating with only low complexity, requiring only linear processing at all terminals, achieving high diversity, and not requiring any short-term channel information at either the relays or the transmitter. For this system, an important performance measure-the so-called outage probability-which gives a fundamental probabilistic performance measure over random communication channels, was considered in Refs. 7 and 8, based on deriving an expression for the moment generating function of the received signal to noise ratio (SNR). The outage probability could then be computed by performing a subsequent numerical Laplace Transform inversion. The exact moment generating function results presented in Refs. 7 and 8, however, are quite complicated functions involving a particular Hankel determinant, and they yield little insight. Moreover, for all but small system configurations, the expressions are somewhat difficult to compute, particularly when implementing the Laplace inversion.
In this paper, we significantly expand and elaborate upon the results of Refs. 7 and 8 by employing analytical tools from random matrix theory and statistical physics. Technically, the challenge is to appropriately characterize the Hankel determinant which arises in the moment generating function computation. In general, the Hankel determinant takes the form, generated from the moments of a certain weight function w(x), 0 ≤ x < ∞, µ k := ∞ 0 x k w(x)dx , k = 0, 1, 2, . . ..
Note that the above is merely one of several possible representations for D n , with equivalent forms regularly found in the fields of statistical physics and random matrix theory.
For the problem at hand, we are faced with a Hankel determinant that is generated from a weight function of the form, which is a "two-time" deformation of the classical Laguerre weight, x α e −x . The parameters in this weight, satisfy the following conditions: For our problem both α and N s are integers, where α := |N R − N D | is the absolute difference between the number of relay and destination antennas, N R and N D respectively. The parameters α and N s are determined by the model, but in fact can be extended to take non-integer values such that α > −1 and N s > 0. 9 We will apply two different methods to characterize this Hankel determinant. For the second method, we derive an approximation for the Hankel determinant using the general linear statistics theorem obtained in Ref. 21, based on Dyson's Coulomb Fluid models. 15,[22][23][24][25] These results are essentially the Hankel analogue of asymptotic results for Toeplitz determinants. 26 Once again, these techniques have been used extensively, particularly in statistical physics, though only very recently have they been applied to address problems in communications and information theory. 13,27 It is important to note that in addition to the two methodologies advocated above, there exists other integrable systems approaches which can be used for characterizing Hankel determinants. For example, D n may be written in an equivalent matrix integral formulation, from which a 'deform-and-study' or isomonodromic deformation approach may be adopted. Essentially, this idea involves embedding D n into a more general theory of the τ-function, 28,29 using bilinear identities and linear Virasoro constraints. For details, see Refs. [30][31][32]. Yet another approach is to characterize D n through the use of Fredholm determinants as employed by Tracy and Widom. 10 Both these exact, non-perturbative integrable systems methods have their own specific advantages and disadvantages, and in most cases lead to non-linear ordinary differential or partial differential equations (ODE/PDEs) satisfied by the Hankel determinant. However, these equations are usually of higher order (equations of Chazy type usually appear), from which first integrals have to be found to reduce to a second order ODE. The advantage of the ladder operator method is that closed-form second order equations are directly obtained for a quantity related to the Hankel determinant, bypassing the need to find a first integral, or evaluate any multiple integral.
The rest of this paper is organized as follows. In Sections I A and I B, we present a brief discussion of the model which underpins the MIMO-AF wireless communication system of interest, and some basic measures which are employed to quantify system performance. Then, in Sections I C and I D, we go on to introduce the moment generating function and cumulants of interest, and pose the key mathematical problems to be dealt with in the remainder of the paper.
In Section II, we establish an exact finite n characterization of the Hankel determinant D n , employing the theory of orthogonal polynomials and their ladder operators.
In Section III, we introduce the Coulomb Fluid method, where we first give a brief overview of its key elements, following Refs. 25, 33, 13 and 21. In Section III B, we compute an asymptotic (large n) approximation for the Hankel determinant D n , and thus a corresponding characterization for the moment generating function of interest. In Section III C, the Coulomb Fluid representation for the moment generating function of the received SNR is shown to yield extremely accurate approximations for the error performance of OSTBC MIMO systems with AF relaying, even when the system dimensions are particularly small.
We also employ our analytical results to compute closed-form expressions for the cumulants of the received SNR, first via the Coulomb Fluid method in Section III D, and then presenting a refined analysis based on Painlevé equations in Sections IV and V.
Subsequently, we give an asymptotic characterization of the moment generating function, valid for scenarios for which the average received SNR is high, deriving key quantities of interest to communication engineers, including the so-called diversity order and array gain. These results are, once again, established via the Coulomb Fluid approximation in Section VI, and subsequently validated with the help of a Painlevé characterization in Section VII.

A. Amplify and Forward Wireless Relay Model
Here we briefly recall the background for the model, having been developed in Ref. 7. The dual hop MIMO communication system features a source, relay and destination terminal, having N s , N R and N D antennas respectively. In the process of transmitting a signal, each transmission period is divided into two time slots. In the first time slot, the source transmits to the relay. The relay then amplifies its received signal subject to an average power constraint, prior to transmitting the amplified signal to the destination terminal during the second time-slot. We assume that the source and destination terminals are sufficiently separated such that the direct link between them is negligible (i.e., all communication is done via the relay).
Let H 1 ∈ C N R ×N s and H 2 ∈ C N D ×N R represent the channel matrices between the source and relay, and the relay and destination terminals respectively. Each channel matrix is assumed to have uncorrelated elements distributed as C N (0, 1 A situation is considered where the source terminal transmits data using a method called OS-TBC encoding. 34 In this situation, groups of independent and identically distributed complex Gaussian random variables (referred to as information "symbols") s i , i = 1, . . ., N are assigned via a special codeword mapping to a row orthogonal matrix X = (x 1 , . . ., x N P ) ∈ C N s ×N P satisfying the power constraint E x k 2 =γ, where N P is the number of symbol periods used to send each codeword.
Since it takes N P symbol periods to transmit N symbols, the coding rate is then defined as The received signal matrix at the relay terminal at the end of the first time slot, Y = (y 1 , . . . , y N P ) ∈ C N R ×N P , is given by 1 The notation C N (µ, σ 2 ) represents a complex Gaussian distribution with mean µ and variance σ 2 .
where N ∈ C N R ×N P has uncorrelated entries distributed as C N (0, 1), representing normalized noise samples at the relay terminal. The relay then amplifies the signal it has received by a constant In the above,b is the total power constraint imposed at the relay, i.e., E Gy k 2 ≤b, whilstγ represents the average received SNR at the relay. The received signal R ∈ C N D ×N P at the destination terminal at the end of the second time slot is then given by where W ∈ C N D ×N P has uncorrelated entries distributed as C N (0, 1), representing normalized noise samples at the destination terminal.
Based on the above relationship, standard linear OSTBC decoding can be applied (see, e.g., Ref. 34). This results in decomposing the matrix model (9) into a set of parallel non-interacting single-input single-output relationships given bỹ where η l is distributed as C N (0, 1), with ||H|| F representing the Frobenius norm (or matrix norm) ofH.
The quantity that we are interested in, the instantaneous SNR for the lth symbol γ l , can then be written as

B. Wireless Communication Performance Measures
Here we recall some basic measures which are employed to quantify the performance of wireless communication systems. One of the most common measures is the so-called symbol error rate (SER), which quantifies the rate in which the transmitted symbols (or signals) are detected incorrectly at the receiver.
where M γ (·) is the moment generating function of the instantaneous SNR γ at the receiver, and Θ = π(M − 1)/M and g MPSK = sin 2 (π/M) are modulation-specific constants. As also presented in Ref. 36, the SER may be approximated in terms of M γ (·) but without the integral, via the following expression: In addition to the SER, another useful quantity is the so-called amount of fading (AoF), which serves to quantify the degree of fading (i.e., the level of randomness) in the wireless channel, and is expressed directly in terms of the cumulants of γ. Specifically, the AoF is defined as follows: where κ 1 and κ 2 are the mean and variance of γ respectively. As shown in Ref. 37, for example, this quantity can be used to describe the achievable capacity of the wireless link when the average SNR is low.
For our AF relaying system under consideration, it is clear that a major challenge is to characterize the moment generating function and also the cumulants of the instantaneous SNR given in (10). This is the key focus of the paper.

C. Statistical Characterization of the SNR γ
Here we introduce the moment generating function and cumulants of interest, and pose the key mathematical problems to be dealt with in the remainder of the paper.
Defining the positive integers, where n = 0 due to physical considerations, it was shown in Ref. 7 that the moment generating function of the instantaneous SNR (10) may be expressed as the multiple integral, From the above, we may then compute the l th cumulant of γ via Let We now consider the moment generating function (15) as a function of two variables, (T,t), or (s,t), since T = t/(1 + cs). We may then write the moment generating function (15) as We see that this involves the weight with the parameters T and t which is a deformation of the classical generalized Laguerre weight w (α) Note that w AF (x,t,t) = w (α) Lag (x). The multiple integral representation (21) is expressed as a ratio of Hankel determinants: where µ j (T,t) is the jth moment of the weight For N s ∈ N, the moments of the weight w AF are expressed in terms of the Kummer function of the second kind, The Hankel determinant for T = t, namely, can be regarded as a normalization constant, so that M γ (t,t) = 1, and its closed form expression is well-known (see Ref. 38).
From a simple change of variables, the cumulants can be calculated by transforming (16) using (20) as Note that an equivalent expression to the result given in (26) has been previously reported in Ref. 7. This result, whilst having numerical computation advantages when the system dimensions are small, has a number of shortcomings. Most importantly, its complexity does not lead to useful insight into the characteristics of the probability distribution of the SNR γ, and it does not clearly reveal the dependence on the system parameters. Moreover, in this form, the expression is not amenable to asymptotic analysis, e.g., as the number of antennas grow sufficiently large, and in such cases its numerical computation becomes unwieldy.
In the following, we will seek alternative simplified representations with the aim of overcoming these shortcomings. The key challenge is to characterize the Hankel determinant in (24 where z := t/T and M Shot-Noise (T, z) is the Hankel determinant generated via Of course, T , z and N s have different interpretations in this situation.
We mention here that the operator theory approach of Ref. 40

D. Alternative Characterization of the SNR γ
An alternative characterization for the moment generating function that will also prove to be useful is derived as follows. From the definitions of m and n, there are two possible choices for the parameter t, We first consider the sub-case N R ≤ N D . To this end, starting with (21) and with the change of variables, x i → nx i , i = 1, . . . , n, we obtain 2 where Note that in terms of s, T ′ can be written as Equivalently, we can write (32) as where and Remark 2. We introduce the variables T ′ and t ′ in order to account for the n-dependence of the variables T and t. This is important since the above representation will be useful for deriving a large n approximation for the moment generating function based on the Coulomb Fluid linear statistics approach in Section III.
For the sub-case N R > N D , we have from (31) that t is instead given by where β = m n − 1. Hence equations (37)- (39) and the results of Sections III-V are valid for the case N R > N D upon transforming t ′ and T ′ to t † and T † via and respectively. In this case, we may write T † in terms of the variable s as

FRAMEWORK
Using the ladder operator framework, and treating T and t as independent variables, the result below is proved in Appendix A. This gives a PDE satisfied by log M γ (T,t) in the variables T and t.
Theorem 1. Let the quantity H n (T,t) be defined through the Hankel determinant D n (T,t) as where the second equality follows from (24). Then H n (T,t) satisfies the following PDE: 3 : 3 Dropping the T and t dependence notation for the sake of brevity. where and We have that v = t and s = 1 c ( t T − 1) and consequently, Under the change of variables, the original PDE (46) becomes the following PDE in the variables (s, v): where A * 1 (H n ) and A * 2 (H n ) obtained from a corresponding transformation are and We note that in the PDE (54), the only higher order derivatives are the mixed derivatives with respect to s and v, and second order derivatives with respect to v.
It is of interest to note that the Hankel determinant in (21) or (24) has the following asymptotic forms: where the generating weight is the "first-time" deformation of the Laguerre weight In previous work by the authors, 13 the Hankel determinant D n [w dLag (·, T, α, N s )] was shown to arise in the analysis of the moment generating function of the single-user MIMO channel capacity, and was shown to involve a Painlevé V differential equation.

III. COULOMB FLUID METHOD FOR LARGE n ANALYSIS
In this section, we make use of the Coulomb Fluid method, which is particularly convenient when the size of the matrix n is large, to describe the MIMO-AF problem. The idea is to treat the eigenvalues as identically charged particles, with logarithmic repulsion, and held together by an external potential. When n, in this context the number of particles, is large, this assembly is regarded as a continuous fluid, first put forward by Dyson,[22][23][24] where the eigenvalues were can be found in Ref. 21 and will be used extensively in this paper. The main benefit of this approach, based on singular integral equations, is that it leads to relatively simple expressions for characterizing our moment generating function.
In Section III A, for completeness, we give a brief overview of the key elements of the Coulomb

A. Preliminaries of the Coulomb Fluid Method
We start by considering the ratio in the moment generating function expression (37), noting that it is of the form where and F n := − log Z n is known as the free energy.
This expression embraces the moment generating function expression (37) with appropriate selection of the functions v 0 (x) and f (x, T ′ ,t ′ ). A key motivation for writing our problem in this form is that it admits a very intuitive interpretation in terms of statistical physics, as originally observed by Dyson (see Refs. [22][23][24]. Specifically, if the eigenvalues x 1 , . . . , x n are interpreted as the positions of n identically charged particles, then Φ(x 1 , . . . , x n ) is recognized as the total energy of the repelling charged particles, which are confined by the external potential nv 0 (x). The function f (x, T ′ ,t ′ ) acts as a perturbation to the system, resulting in a modification to the external potential. The quantity F n (T ′ ,t ′ ) may be interpreted as the free energy of the system under an external perturbation f (x, T ′ ,t ′ ), with F n (t ′ ,t ′ ) the free energy of the unperturbed system.
For sufficiently large n, the system of particles, following Dyson, may be approximated as a continuous fluid where techniques of macroscopic physics and electrostatics can be applied. For large n we expect the external potential nv 0 (x) to be strong enough to overcome the logarithmic repulsion between the particles (or eigenvalues), and hence the particles or fluid will be confined within a finite interval to be determined through a minimization process. For this continuous fluid, we introduce a macroscopic density σ (x)dx, referred to as the equilibrium density. Since v 0 (x) is convex for x ∈ R, this density is supported on a single interval denoted by (a, b), to be determined later (see Ref. 33 for a detailed explanation). The equilibrium density is obtained by minimizing the free-energy functional: With Frostman's Lemma, 42 (p. 65) the minimizing σ (x)dx can be characterized through the integral equation where x ∈ [a, b] and A is the Lagrange multiplier for the normalization condition (64), which can be interpreted as the chemical potential of the fluid. 25,33 Noting that the integral equation above has a logarithmic kernel, taking a derivative with respect to x ∈ (a, b) converts it into a singular where P denotes Cauchy principal value.
Noting the structure (in n) of the left-hand side of (66), it is clear that σ (·) must take the general form: where σ 0 (x)dx is the density of the original system in the absence of any perturbation, while σ c (x, T ′ ,t ′ ) represents the deformation of σ 0 (x) caused by f (x, T ′ ,t ′ ). Furthermore, to satisfy Substituting (67) into (66), and comparing orders of n, we see that Following Ref. 33, where the choice for the solution for σ 0 has been extensively discussed based on the theory described in Ref. 43; the solution to (69) subject to the boundary condition σ 0 (a) = σ 0 (b) = 0 reads together with a supplementary condition, The solution to (70) Finally, the normalization condition The end points of the support of the density σ 0 (x), a, and b, are determined by (72)  With the above results, for sufficiently large n, we may approximate the ratio (60) as (see

Ref. 21 for more details)
where In the sequel, we will find explicit solutions for these quantities.
From here, the l-th cumulant can be extracted from the formula The objective of the subsequent analysis is to evaluate the quantities S AF 1 and S AF 2 . As we shall see, we will need to solve numerous integrals which are not readily available. Thus, to aid the reader, we have compiled these integrals in Appendix B.
We start by considering the end points of the support a and b. These are determined by equations (72) and (74). With the integral identities (B1)-(B3) in Appendix B, we obtain, after a few easy steps, which leads to The limiting density σ 0 (x) in (71) can be computed using the integral identity (B2) as which is the Marĉenko-Pastur law. 13,46,47 Meanwhile, with the aid of the integral identity (B4) in Using σ 0 (x) and invoking the integral identities (B5)-(B7) in Appendix B, gives Moreover, using σ c (x, T ′ ,t ′ ) and invoking the integral identities (B2), (B8) and (B12) in Appendix B gives (87)

C. SER Performance Measure Analysis Based on Coulomb Fluid
Combining (78)  The curves labeled "Coulomb (Exact SER)" were generated by substituting our Coulomb Fluid approximation into the exact expression of the SER (11), whilst the curves labeled "Coulomb (Approx SER)" are generated by substituting our Coulomb Fluid approximation into the approximate expression for the SER (12).
From these curves, it is remarkable that the Coulomb Fluid based approximations, derived under the assumption of large n, very accurately predict the SER even for small values of n. This is evident, for example, by examining the set of curves corresponding to the configuration (2, 3, 2), for which n = 2. In fact, even for the extreme cases with n = 1, the Coulomb Fluid curves still yield quite high accuracy.
We note that the results in Fig. 1  The key difference therein was that the analytic curves were based on substituting into either (11) or (12) an equivalent determinant form of the moment generating function to that given in (25) and (26). That representation, involving a determinant of a matrix with elements comprising sums of Kummer functions, is far more complicated than the Coulomb Fluid representation, which is a simple algebraic equation involving only elementary functions. The fact that the Coulomb Fluid representation yields accurate approximations for the SER for small as well as large numbers of antennas makes it a useful analytical tool for studying the performance of arbitrary MIMO-AF systems.

D. Coulomb Fluid Analysis of Large n Cumulants of SNR
In this subsection we use the Coulomb Fluid results to study the asymptotic cumulants of the SNR. We suppose β is fixed, and distinguish between two cases, β = 0 and β > 0.

Case 1: β = 0
Here N R = N D . Then from (83), a = 0 and b = 4, leading to and From the closed-form approximation of the moment generating function (78) the cumulants can be extracted using (81).
Alternatively, recall that the variable T ′ is related to s via where c =γ/(RN s ). Therefore we expand log M γ (t ′ /(1 +cs),t ′ ) obtained from the Coulomb Fluid formalism (88) and (89) in a Taylor series about s = 0. From (78), where the first few cumulants are and κ 6 (t ′ ) (96)

Case 2: β > 0
In this case, the end points of the support, a and b, are given by (83). Hence, going through the same process, the first five cumulants are and The accuracy of these approximations is demonstrated in Fig.2, where they are compared with simulation results generated by numerically computing the exact cumulants via Monte Carlo methods. As before, different antenna configurations are shown, as represented by the form . From these curves, the accuracy of our Coulomb Fluid based approximations is quite remarkable, even for small values of n.
We note that exact expressions for κ 1 and κ 2 (i.e., the mean and variance) were derived pre-

IV. POWER SERIES EXPANSION FOR H n (s, v)
In this section we seek power series solutions for (54), to determine the cumulants. We begin by removing the square roots of the PDE and arrive at the following lemma. and then the PDE (54) may be rewritten in an equivalent square-root-free form, as where A 1 and A 2 are given by (55) and (56) respectively. Now, recall from (51) that the function H n (s, v) is related to the moment generating function (21) through and note that the logarithm of the moment generating function has the expansion, (105), H n (s, v) has the following expansion in s, where ′ denotes differentiation with respect to v, and a j (v) is related to the jth cumulant κ j (v) by A. Analysis of κ 1 In this subsection, the computation of the mean κ 1 is presented. To this end, upon substituting (107) into the PDE (104), keeping the lowest power of s in the resulting expansion gives a highly nonlinear ODE expressed in the factored form where Ψ (1) 4 We introduce a j in order to make a clearer comparison to the Coulomb Fluid model and 7 without having to exactly specify a value for the constant c, since it falls out of the subsequent expansion. and Ψ (1) Hence, (109) is equivalent to either Ψ (1) (1) 1 = 0 has two types of solutions, the first of the form a 1 (v) = c 1 v + c 2 where c 1 and c 2 are constants to be determined. The second is more complicated but can be succinctly written as where c 1 and c 2 are integration constants, and f and g are algebraic functions (we omit these for brevity).
Whilst these are valid mathematical solutions, it is clear that by taking n large, they differ drastically from that predicted by the Coulomb Fluid method and, for small n, they differ with the solutions obtained in Ref. 7. This suggests that these are not the solutions of interest to our problem; thus, we set them aside and henceforth focus on Ψ (1) 2 = 0. It becomes evident at this point that we require initial conditions for the cumulants κ j (v)/(c j N s ), since the PDE (104) gives rise to a system of ODEs satisfied by κ j (v)/(c j N s ). However, (104) is not supplemented by any initial conditions from which one deduces the initial conditions at t ′ = 0 of the system of ODEs.
As we shall see, results from the Coulomb Fluid analysis (97)-(99) are crucial, as these will specify initial conditions at t ′ = 0. Consequently, the Coulomb Fluid results are in fact leading order approximations to the exact results. Without the Coulomb Fluid analysis, each cumulant κ j (v)/(c j N s ) would carry unknown constants of integration.
We now examine the equation Ψ (1) 2 = 0 in the large n region, which leads to an asymptotic characterization of a 1 (v) and thus κ 1 . To proceed, we scale the variable v by v = nt ′ . Also, note that α = n(m/n − 1) ≡ nβ . We find that a 1 (t ′ ) satisfies the following ODE: For large n, keeping the leading order term, we arrive at whose solutions are where C 1 is a constant of integration. We retain the solution that is bounded as t ′ → ∞, Comparing this with the corresponding result (97) obtained from the Coulomb fluid approximation, we see that κ CF 1 (0) = cN s n, which suggests a 1 (0) = n and hence C 1 = n(1 + β /2). Therefore, the first cumulant κ 1 to O(n) becomes κ Large n B. Analysis of κ 2 Having characterized κ 1 , we now turn to the variance κ 2 . To this end, equating the coefficients of the next lowest power of s to zero in the expansion of the PDE (104), results in an ODE which can be factored into the following form The ODE Ψ (1) 1 = 0, makes a reappearance, which we set aside, leaving the ODE Ψ where Ψ 1 , Ψ 2 and Ψ 3 are given by and Clearly, Ψ 1 depends on a 1 only, whilst Ψ 2 and Ψ 3 depend on both a 1 and a 2 . To extract the large n behavior of κ 2 , we replace v by nt ′ in (118) and make use of the large n formula for a 1 (t ′ ) in (116), resulting to a highly non-linear ODE satisfied by a 2 (t ′ ). To proceed further, keeping the highest powers of n, a first order equation is obtained for a 2 (t ′ ), Integrating this leads to an expression for a 2 (t ′ ) which, again yields an unknown integration constant. This constant can be determined through a comparison with the Coulomb Fluid results for κ 2 (t ′ ) in (98) at t ′ = 0, giving a 2 (0) = n.
The above analysis gives the leading order characterization of the variance, C. Beyond κ 1 and κ 2 The same procedure for computing κ 1 and κ 2 easily extends to higher cumulants. Here we give the leading order formula for κ 3 , κ 4 and κ 5 , and κ Large n where c j (s,t) := t α+ j+1 Γ(α + j + 1) which was derived in Ref. 7. For small fixed integer values of n and N s , e.g., n = 2, 3 and N s = 10, and n = 4, 5 and N s = 1, the above determinant can computed without much difficulty, from which the cumulants follow. It can be seen that (127), satisfied third order ODEs for a 1 (t) given by (111) and a 2 (t) given by (118). Similar results hold for the higher cumulants, which provide a consistency check.

V. LARGE n CORRECTIONS OF CUMULANTS OBTAINED FROM COULOMB FLUID
We have shown that the PDE (104) satisfied by may be used to generate a series of non-linear ODEs that are satisfied by the cumulants κ j (v).
Under a large n assumption, where v = nt ′ , the first few of these ODEs are approximated as first order ODEs for κ Large n l (t ′ ), whose solutions matched exactly with that obtained from the Coulomb Fluid analysis for the cumulants κ CF l (t ′ ). We give a method in this section where the non-linear ODEs generated from the PDE (104) are employed systematically to obtain "correction terms" to the Coulomb Fluid results.

A. Large n expansion of κ 1
We assume the first cumulant or κ 1 has the following large n expansion Substituting the above into (113), and setting the coefficients of n j to zero, a system of first order ODEs for e k (t ′ ) is obtained. These are are solved successively starting from e −1 (t ′ ), followed by e 0 (t ′ ) and so on.
As a result e −1 (t ′ ) is found to satisfy the ODE and has two solutions both satisfying the initial condition e 1 (0) = 1.
We retain the solution that is finite at infinity to match with that obtained from the Coulomb Fluid computation, namely, (97).
Continuing, we set the coefficient of n 3 to zero, which implies, Setting the coefficient of n 2 to zero gives rise to an ODE involving e ′ −1 (t ′ ), e ′′ −1 (t ′ ), e ′′′ −1 (t ′ ), e ′ 0 (t ′ ) and e ′ 1 (t ′ ). Simplifying, with e −1 (t ′ ) given by (132) and e ′ 0 (t ′ ) given by (133), we find that e 1 (t ′ ) satisfies: Setting the coefficient of the next lowest power of n to zero gives rise to an ODE involving e , and e ′ 2 (t ′ ). From the expression of e −1 (t ′ ) and e ′ 0 (t ′ ) given by (132) and (133) respectively, we find Continue with this process, the next three terms e 3 (t ′ ), e 4 (t ′ ) and e 5 (t ′ ) are found to satisfy the following equations: where l Solving ODEs (133)-(138) with initial conditions e k (0) = 0, we obtain, where A (1)

B. Large n Expansion of κ 2
In computing a series expansion for the variance, we proceed in a similar way as was just done for the mean. First, we substitute the expansion for a 1 (t ′ ) from (141) and into (118). Equating the coefficients of n k to zero, a system of first order ODEs for f k (t ′ ) are found.
In this case, setting the coefficient of n 12 to 0 results in a first order ODE for f −1 (t ′ ): The solution of this, with the initial condition f −1 (0) = 1, reads Continuing the process, we obtain an ODE for f 0 (t ′ ): The solution with the condition f 0 (0) = 0 reads from which it can be immediately seen that and we recover κ CF 2 (t ′ ) found previously. The ODEs satisfied by f k (t ′ ), k = 1, 2, 3, 4 are reported in Appendix C, equations (C1)-(C4).
These will allow us to compute f j (t ′ ) with appropriate initial conditions. Briefly, the idea is that to determine the ODE satisfied by the k-th correction term f k (t ′ ), for k odd, the preceding 1 2 (k + 1) ODEs satisfied by f j (t ′ ), j = −1, 1, 3 . . ., k are employed. For even k, the previous k 2 ODEs satisfied by f j (t ′ ), j = 0, 2, 4 . . ., k are employed. Going through the procedure described, we find that κ 2 (t ′ ) has the following large n expansion, where A (2) C. Beyond κ 1 and κ 2 The procedure adopted above for computing the large n expansion series for κ 1 (t ′ ) and κ 2 (t ′ ) easily extends to the higher cumulants. By way of example, here we focus on the third cumulant κ 3 . In this case, an asymptotic expansion for a 3 (t ′ ) is assumed, along with the initial conditions, g −1 (0) = 2 and g k (0) = 0 for k = 0, 1, 2, 3, . . ..
In this case, we find that the large n expansion of the third cumulant reads For ease of reference, the ODEs satisfied by g k (t ′ ), k = 1, 2, 3, 4, are placed in Appendix C, In summary, this entire section has shown that by expanding the cumulants κ 1 (t ′ ), κ 2 (t ′ ) and κ 3 (t ′ ) into asymptotic series in n, the Coulomb Fluid results are recovered as the leading order contributions in the large n scenario. It is also seen that, by examining the finite-n correction terms for each cumulant, no terms of O(n 2 ) or higher are present within the expansions.

VI. ASYMPTOTIC PERFORMANCE ANALYSIS BASED ON COULOMB FLUID
In this section, we return to the analysis of the moment generating function, and consider the high SNR scenario (i.e., asγ → ∞). To this end, we will study the Coulomb Fluid based approximation derived in Section III (to be complemented in Section VII through analysis based on Painlevé equations), where the variables T ′ and t ′ are taken to be dependent onγ, namely, Note that asγ → ∞, To obtain the desired high SNR expansion, we compute the moment generating function M γ as s → ∞ andγ → ∞. The Coulomb Fluid based representation (78), when expressed in terms of where S AF 1 and S AF 2 are given by (86) and (87) respectively. Our goal is to compute the expansion of M γ (s) as s → ∞ andγ → ∞, which turns out to be an expansion in (γs) −1 . It turns out that the cases β = 0 and β = 0 behave fundamentally differently, and as such, these are treated separately.

A. The Case of β = 0
With β = 0, we have N R = n. An easy computation shows that M γ admits the following expansion: where A 0 , A 1 , A 2 , A 3 , . . . are constants independent of s andγ. The leading exponent d is given by whilst the first few A ℓ are We have refrained from presenting A k , k ≥ 4, as these are rather long.

High SNR Analysis of the Symbol Error Rate (SER)
Based on (11), the SER of MPSK modulation can be expanded at high SNR using (165), resulting in where Considering the first order expansion, following Ref. 48, we may write where we identify as the so-called diversity order, and identify the factor Hence, using (176), G a admits the simplified form The high SNR results above are illustrated in Fig. 3. The "Simulation" curves are based on numerically evaluating the exact SER relation (11); the "Coulomb Fluid (Exact)" curves are based on substituting (164) into (11) and numerically evaluating the resulting integral; the "Coulomb Fluid (Leading term only)" curves are based on (174); whilst "Coulomb Fluid (Leading 4 terms)" curves are based on the first four terms of (172). The leading-order approximation is shown to give a reasonably good approximation at high SNR, whilst the additional accuracy obtained by including a few correction terms is also clearly evident.

High SNR Analysis of the Probability Density Function of γ
With the moment generating function expansion given above in (165), we may also readily obtain an approximation for the probability density function (PDF) of γ, denoted f γ (x), by direct Laplace Transform inversion. In particular, we obtain For the case of very largeγ, with the leading term gives an approximation for the PDF deep in the left-hand tail. Of course, with the inclusion of more terms, a more refined approximation is obtained.

B. The Case of β = 0 (with N R < N D )
For the situation where β = 0, two sub-cases arise. This first is N R < N D , for which N R = n; the second is N R > N D , for which N R = (1 + β )n.
In the following, we will focus on the first sub-case, N R < N D . It turns out however, that our results also apply for N R > N D upon transforming the quantity ν to ν * by where ν * (> 0) is interpreted as the (fixed) scaling factor betweenb andγ, i.e., The moment generating function given by (164) admits an expansion distinct from the β = 0 case which does not have fractional powers of 1/(γs), reading and In this situation, the sub-leading terms are very complicated, however, the jth term in the expansion can be written in the following form: where E j and F j also depend upon ν, N s , n and β .
In A 2 , E 2 , and F 2 are given by In A 3 , E 3 , and F 3 are given by In A 4 , E 4 , and F 4 are given by

High SNR Analysis of the Symbol Error Rate (SER)
Based on (11), the SER of MPSK modulation can be expanded at high SNR using (183) into where Note that here, in contrast to (173), the exponent r is a positive integer. As such, (195) admits the following closed-form solution: (Ref. 49, 2.513.1) As before, a first-order approximation is of key interest, giving where we identify the diversity order and the array gain The result for G d above is consistent with a result obtained via a different method in Ref. 8, whilst the expression for G a appears new. The high SNR results above, for the case β = 0, are illustrated in Fig. 4. As before, the "Simulation" curves are based on numerically evaluating the exact SER relation (11), and the "Coulomb Fluid (Exact)" curves are based on substituting (164) into (11) and numerically evaluating the resulting integral. Moreover, the "Coulomb Fluid (Leading term only)" curves are based on (197), whilst the "Coulomb Fluid (Leading 5 terms)" curves are based on the first five terms of (194). Again, the leading-order approximation is shown to give a reasonably good approximation at high SNR, whilst the additional accuracy obtained by including a few correction terms is also evident.

High SNR Analysis of the Probability Density Function of γ
As before, based on the moment generating function expansion (183), applying for β = 0, we can immediately take a Laplace inversion to obtain the following high SNR representation for the PDF of γ, We are mainly interested in the leading order term, and so we write the above as Note that despite the similarity with (180), interestingly, these results do not coincide upon taking β → 0 in (201), due to the differences in d and A 0 . This seems to indicate that the double asymptoticsγ → ∞ and β → 0 are non-commutative.

VII. CHARACTERIZING A 0 THROUGH PAINLEVÉ V
In this section, we obtain the leading term of the large s and largeγ expansion (183) from a Painlevé V differential equation, thus demonstrating the accuracy of the Coulomb Fluid approximation. We will focus in this section on the case β = 0.
Recall that the moment generating function of SNR γ, regarded a function of s and t, given by (21), reads, where K n,α is a normalization constant in (27). We also recall that c and t are given by respectively. So asγ → ∞, we note that A simple computation shows that M γ (s,t) admits the following expansion for largeγs: 1 (γs) nN s +1 nD n w dLag (·,t, α − N s , N s ) where w dLag (x,t, α − N s , N s ) is the deformation of the classical Laguerre weight, i.e., and The condition α − N s > −1 is required for the validity of the orthogonality relation with respect to the generalized Laguerre weight w Ref. 9). This, in turn, will ensure the validity of the first two terms in (204), while for the multiple integral in the third term to converge, the condition α − N s > 0 is required.
For the problem at hand, α and N s are integers, satisfying α ≥ 0 and N s > 0. Therefore α ≥ N s implies that the result for A 0 presented below is valid for α > 0.
Comparing the expansion (204) with (183), we see that the diversity order is and We see that A 0 is up to a constant the Hankel determinant which generates a particular Painlevé V and shows up in the single-user MIMO problem studied in Ref. 13.
We obtain A 0 , through a large n expansion of for α > 0. We will see that A 0 precisely matches that obtained in (185).
From Ref. 13 we learned that the logarithmic derivative of D n [w dLag (·,t, α − N s , N s )] with respect to t, satisfies the Painlevé V: where ′ denote derivative w.r.t. t.
We restrict to the case where N R < N D , for which N R = n. The case where N R > N D can be considered in a similar fashion as outlined in Section VI B. Setting α = nβ in the above equation, where β = m n − 1 is a fixed positive number, and with the change of variable an easy computation shows that Y n (ν) := H n (n/ν) satisfies ν 4 where ′ denotes derivative with respect to ν.
Assuming that the coefficient of n k is zero, we find that the equation c −2 = 0 gives us With MAPLE, the solutions of the above ODE for p −1 (ν) are found to be where C 2 is a constant of integration. 5 .
The equation c −1 = 0 is a coupled differential equation involving both p −1 (ν) and p 0 (ν) given by With p −1 (ν) given by (217), chosen to match the result from the Coulomb Fluid (185), we find that the first order ODE in p 0 (ν) has the solution, We disregard the second and third solutions for p −1 (ν); as these would lead to p 0 (ν) which do not generate the A 0 in agreement with that obtained from the Coulomb Fluid method.
Using a method similar to the cumulant analysis of Section V, we can also compute the first correction term to A 0 (i.e., the quantity q 1 (ν)) using (222), which reads Higher order corrections could also be obtained in a similar way.

VIII. CONCLUSION
In this paper, we have introduced two methods for characterizing the received SNR distribution in a certain MIMO communication system adopting AF relaying. We showed that the mathe- i.e., where h n is the square of the L 2 norm of P n (x). The Hankel determinant (1) is reduced to the following product: Our convention is to write P n (x) as P n (x) := x n + p 1 (n)x n−1 + p 2 (n)x n−2 + · · · + P n (0).
Hence, this implies that properties of the Hankel determinant may be obtained by characterizing the class of polynomials which are orthogonal with respect to w AF (x, T,t), over [0, ∞). It is clear that the coefficients of the polynomial P n (x), p i (n), will depend on T , t, α and N s ; for brevity, we do not display this dependence.
The main aim is to determine these unknown recurrence coefficients α n and β n from the given weight. Substituting (A4) into the three term recurrence relation results in where p 1 (0) := 0. Taking a telescopic sum gives Moreover, combining the orthogonality relationship (A2) with the three term recurrence relation leads to which can also be expressed in terms of the Hankel determinant D n in (A3) through since h n = D n+1 D n .

Ladder Operators, Compatibility Conditions, and Difference Equations
In the theory of Hermitian random matrices, orthogonal polynomials plays an important role, As we now show, there is a recursive algorithm that facilitates the determination of the recurrence coefficients α n and β n . This is implemented through the use of so-called "ladder operators" as well as their associated compatibility conditions. This approach can be traced back to Laguerre From the weight function w AF (x, T,t), one constructs the associated potential v(x) through and therefore, As shown in Ref. 15, a pair of ladder operators, to be satisfied by our orthogonal polynomials of interest, are expressed in terms of v(x) and are given by where x − y P n (y)P n−1 (y)w AF (y)dy, and where, for the sake of brevity, we have dropped the t, T dependence in w AF .
Moreover, there are associated fundamental compatibility conditions to be satisfied by A n (x) and B n (x), which are given by 16 These were initially derived for any polynomial v(x) (see Refs. [54][55][56], and then were shown to hold for all x ∈ C ∪ {∞} in greater generality. 16 We now combine (S 1 ) and (S 2 ) as follows. First, multiplying (S 2 ) by A n (x), it can be seen that the RHS is a first order difference, while (x − α n )A n (x) on the LHS can be replaced by B n+1 (x) + B n (x) + v ′ (x) from (S 1 ). Then, taking a telescopic sum with initial conditions leads to the useful identity The condition (S ′ 2 ) is of considerable interest, since it is intimately related to the logarithm of the Hankel determinant. In order to gain further information about the determinant, we need to find a way to reduce the sum to fixed number of quantities; for which, (S ′ 2 ) ultimately provides a way of going forward.

Remark 6. Since our v ′ (x) is a rational function of x, we see that
is also a rational function of x, which in turn implies that A n (x) and B n (x) are rational functions of x. Consequently, equating the residues of the simple and double pole at x = 0, x = −T , x = −t on both sides of the compatibility conditions (S 1 ), (S 2 ) and (S ′ 2 ), we obtain equations containing numerous n, T and t dependant quantities; which we call the "auxiliary variables" (to be introduced below). The resulting non-linear discrete equations are likely very complicated, but the main idea is to express the recurrence coefficients α n and β n in terms of these auxiliary variables, and eventually take advantage of the product representation (A3) to obtain an equation satisfied by the logarithmic derivative of the Hankel determinant. Now substituting (A14) into (A13) which define A n (x) and B n (x), and followed by integration by parts, we find where R * n (T,t) ≡ N s h n ∞ 0 w AF (y)P 2 n (y) y + t dy, r * n (T,t) ≡ N s h n−1 ∞ 0 w AF (y)P n (y)P n−1 (y) y + t dy, w AF (y)P 2 n (y) y + T dy, r n (T,t) ≡ N s h n−1 ∞ 0 w AF (y)P n (y)P n−1 (y) y + T dy, are the auxiliary variables.

Difference Equations from Compatibility Conditions
Inserting A n (x) and B n (x) into the compatibility conditions (S 1 ), (S 2 ) and (S ′ 2 ), and equating the residues as described above yields a system of 12 equations. Note that there are actually 14 equations, the compatibility conditions (S 1 ) and (S 2 ) also have O(x 0 ) terms which yield equations that lead to 0 = 0.
The compatibility conditions (S 1 ) give the following set of equations, Similarly, the compatibility condition (S 2 ) gives the following set of equations, To obtain a system of difference equations from the compatibility condition (S ′ 2 ) is somewhat more complicated, but carrying out the process, equating all respective residues in (S ′ 2 ) yields six equations. The first three are obtained by equating the residues of the double pole at x = 0, x = −t and x = −T respectively: − 2r n r * n − (2n − N s + α)r n +(2n + N s + α)r * n + n(n + α) = β n − β n (R n R * n−1 + R * n R n−1 ) while the last three are given by equating the residues of the simple pole at x = 0, x = −t and x = −T respectively: n−1 ∑ j=0 R j + (N s − 2r n )(r n − n − r * n ) + αr n T + (N s − 2r n )r * n + N s r n T − t + r n = − β n (R n + R n−1 + 2R n R n−1 )

Remark 7. It can be seen that (A27) is a combination of (A28) and (A29), and serves as a consis-
tency check. We also see that (A28) and (A29) respectively gives us the value of ∑ j R * j and ∑ j R j automatically in closed form. Hence we can also obtain any linear combination of these sums in closed form, which is a crucial step in obtaining a link between a linear combination of the logarithmic partial derivatives of the Hankel determinant, and the quantities β n , r n and r * n . This step would not be possible without (S ′ 2 ), also known as the bilinear identity.

Analysis of Non-linear System
Whilst some of the above difference equations (A18)-(A29) look rather complicated, our aim is to manipulate these equations in such a way as to give us insight into the recurrence coefficients α n and β n . Our dominant strategy is to always try to describe the recurrence coefficients α n and β n in terms of the auxiliary variables R n , r n , R * n and r * n . The first thing we can do is to sum (A18) through (A20) to get a simple expression for the recurrence coefficient α n in terms of the auxiliary variables R n and R * n : α n = T R n − tR * n + α + 2n + 1.
We can immediately see that by taking a telescopic sum of the above from j = 0 to j = n − 1, and recalling equation (A7), gives rise to Looking at (S 2 ) next, we sum equations (A21) through (A23) to get Taking a telescopic sum of the above from j = 0 to j = n − 1 and rearranging yields β n = Tr n − tr * n − p 1 (n).
We are now in a position to derive the following important lemma, which describes the recurrence coefficients α n and β n in terms of the set of auxiliary variables: Lemma 2. The quantities α n and β n are given in terms of the auxiliary variables R n , r n , R * n and r * n as α n = T R n − tR * n + α + 2n + 1, (A34) Proof. Equation (A34) is a restatement of equation (A30).
To obtain (A35), we use the equations obtained from(S ′ 2 ). We eliminate R * n−1 and R n−1 from (A24) using (A25) and (A26) respectively, and then rearrange to express β n in terms of the auxiliary quantities.

Toda Evolution
In this section, n is kept fixed while we vary two parameters in the weight function (22) ∂ t β n = β n (R * n − R * n−1 ).
Applying ∂ T and ∂ t to the orthogonality relation ∞ 0 w AF (x, T,t)P n (x)P n−1 (x)dx = 0 results in the following two relations: ∂ t p 1 (n) = −r * n .
Proof. To obtain (A45), we substitute in for α n in (A43) using (A34) to yield r n − r n+1 = ∂ T α n , Eliminating r n+1 in the above formula using (A20) gives 2r n = 1 − α n − T R n + T ∂ T R n − t∂ T R * n + N s .
We obtain (A46) in a similar method to (A45). This time, we first substitute in for α n in (A44) using (A34). We then proceed to eliminate r * n+1 using (A19), and finally eliminate α n again using (A34) to obtain (A46).
To obtain (A47), we differentiate equation (A33) with respect to T . Following this, we substitute in for ∂ T β n and ∂ T p 1 (n) using (A38) and (A40) respectively to yield T ∂ T r n − t∂ T r * n = β n (R n−1 − R n ).
We then proceed to replace β n R n−1 by r n (r n − N s )/R n using (A26), and finally we eliminate β n in favor of R n , R * n , r n and r * n using (A35) to obtain equation (A47). Similarly, we obtain (A48) by first differentiating (A33) with respect to t. Following this, we substitute in for ∂ t β n and ∂ t p 1 (n) using (A39) and (A41) respectively to yield T ∂ t r n − t∂ t r * n = β n (R * n − R * n−1 ).
We then proceed to replace β n R * n−1 by r * n (r * n − N s )/R * n using (A25), and finally we eliminate β n in favor of R n , R * n , r n and r * n using (A35) to obtain equation (A48), completing our proof.
From Lemma 4, the auxiliary variables r n and r * n appear linearly in equations (A45) and (A46) respectively. By eliminating r n and r * n from (A47) and (A48) using (A45) and (A46), we obtain a second order PDE system for the auxiliary variables R n and R * n . This is stated in the following lemma: Lemma 5. The auxiliary variables R n and R * n satisfy the following second-order PDE system: +2T (tR * n + 1) ∂ T R * n ) + 2t(T R n + 1) ∂ t R * n − 2T 2 R * n ∂ T R n + R n ∂ t R n − 2T (T − t) ∂ t R n +2R * n (T R n − tR * n + t) T R n − tR * n + 2n + α + 1 + (A54)

Toda Evolution of Hankel Determinant
Recall that the moment generating function is related to the Hankel determinant through equation (24): while the Hankel determinant is related to h n (T,t) by the relation D n (T,t) = ∏ n−1 j=0 h j (T,t). We compute the partial derivatives of the Hankel determinant by taking a telescopic sum from j = 0 to j = n − 1 of the partial derivatives of log h n , i.e. equations (A36) and (A37). Using (A3), To continue, we apply t∂ T + t∂ t to equation (A59), and find that (T ∂ T + t∂ t )H n = p 1 (n) + β n , where we have used (A42) to replace (T ∂ T + t∂ t )p 1 (n) by p 1 (n) + β n .
We see that equations (A59) and (A61) can be regarded as a set of simultaneous equations for p 1 (n) and β n . Solving for p 1 (n) and β n , we find that p 1 (n) = H n − n(n + α), and β n = (T ∂ T + t∂ t )H n − H n + n(n + α).
Before completing the proof of Theorem 1, we show here that our Hankel determinant D n (T,t) is related to the τ-function of a Toda PDE.
Substituting the definition (44) into (A62) and together with (A9), we find The above equation is a two parameter generalization of the Toda molecule equation, 57 and hence we identifyD n (T,t) as the corresponding τ-function of the two-parameter Toda equations.
A reduction may be obtained by writingD n (T,t) as Proof. We rewrite log(x + t 2 ) using (B11) and then use (B2) to obtain .
For our problem, where a and b are given by (83) and β ≥ 0, we see that the discriminant of the quadratic form is given by while the constant term is t 2 2 t 2 1 + 2t 1 (2 + β ) + β 2 > 0.
Hence we may invoke (B9) and then take the limit ε → 0 to get Substituting this back into (B13) and after some algebra, we get (B12).