Asymptotic Outage Analysis of HARQ-IR over Time-Correlated Nakagami-$m$ Fading Channels

TL;DR

This paper derives an exact and asymptotic outage probability analysis for HARQ-IR over time-correlated Nakagami-m fading channels, revealing how power, correlation, and rate affect performance.

Contribution

It introduces a novel outage analysis method using conditional Mellin transform for correlated Nakagami fading, providing explicit asymptotic expressions and insights for system optimization.

Findings

Exact outage probability expressed as a weighted sum over independent Nakagami channels.

Asymptotic outage probability quantifies impacts of power, correlation, and rate.

Analytical results enable low-complexity optimal power and rate adaptation.

Abstract

In this paper, outage performance of hybrid automatic repeat request with incremental redundancy (HARQ-IR) is analyzed. Unlike prior analyses, time-correlated Nakagami-$m$ fading channel is considered. The outage analysis thus involves the probability distribution analysis of a product of multiple correlated shifted Gamma random variables and is more challenging than prior analyses. Based on the finding of the conditional independence of the received signal-to-noise ratios (SNRs), the outage probability is exactly derived by using conditional Mellin transform. Specifically, the outage probability of HARQ-IR under time-correlated Nakagami-$m$ fading channels can be written as a weighted sum of outage probabilities of HARQ-IR over independent Nakagami fading channels, where the weightings are determined by a negative multinomial distribution. This result enables not only an efficient truncation approximation of the outage probability with uniform convergence but also asymptotic outage analysis to further extract clear insights which have never been discovered for HARQ-IR even under fast fading channels. The asymptotic outage probability is then derived in a simple form which clearly quantifies the impacts of transmit powers, channel time correlation and information transmission rate. It is proved that the asymptotic outage probability is an inverse power function of the product of transmission powers in all HARQ rounds, an increasing function of the channel time correlation coefficients, and a monotonically increasing and convex function of information transmission rate. The simple expression of the asymptotic result enables optimal power allocation and optimal rate selection of HARQ-IR with low complexity. Finally, numerical results are provided to verify our analytical results and justify the application of the asymptotic result for optimal system design.

Equations

$${\mathbf{y}}_k=h_k{\mathbf{x}}_k+{\mathbf{n}}_k,$$

$$|{h_{k}}|=\sqrt{\frac{{{\sigma_{k}}^{2}}}{m}}\left\|{\sqrt{1-{\lambda_{k}}^{2}}\left({\begin{array}[]{*{20}{c}}{{\vartheta_{k,1}}}\\ \vdots\\ {{\vartheta_{k,m}}}\end{array}}\right)+{\lambda_{k}}\left({\begin{array}[]{*{20}{c}}{{\vartheta_{0,1}}}\\ \vdots\\ {{\vartheta_{0,m}}}\end{array}}\right)}\right\|,$$

$$\frac{\mathrm{E}\left({\left|h_l\right|}^2{\left|h_k\right|}^2\right)-\mathrm{E}\left({\left|h_l\right|}^2\right)\mathrm{E}\left({\left|h_k\right|}^2\right)}{\sqrt{\mathrm{Var}\left({\left|h_l\right|}^2\right)\mathrm{Var}\left({\left|h_k\right|}^2\right)}}={{\lambda }_l}^2{{\lambda }_k}^2,1\le l\ne k\le K,$$

$${\gamma }_k=\frac{P_k{\left|h_k\right|}^2}{{\mathcal{N}}_0},$$

$$p_{out,K}=\mathrm{Pr}\left(I_K<\mathcal{R}\right)=F_{I_K}(\mathcal{R}).$$

$$I_K={\sum }_{k=1}^K{\log }_2\left(1+{\gamma }_k\right).$$

$$p_{out,K}=\mathrm{Pr}\left(G_K\triangleq {\prod }_{k=1}^K\left(1+{\gamma }_k\right)<2^{\mathcal{R}}\right)=F_{G_K}\left(2^{\mathcal{R}}\right),$$

$$f_{{\gamma }_k|T}\left(x_k|t\right)$$

$$f_{G_K|T}\left(x|t\right)$$

$$f_{{\mathcal{A}}_{\mathbf{\ell }}}(x)$$

$$\displaystyle=\frac{1}{{\prod\nolimits_{k=1}^{K}{{\Omega_{k}}}}}Y_{0,K}^{K,0}\left[{\left.{\begin{array}[]{*{20}{c}}-\\ {\left({0,1,\frac{1}{{{\Omega_{1}}}},m+{\ell_{1}}}\right),\cdots,\left({0,1,\frac{1}{{{\Omega_{K}}}},m+{\ell_{K}}}\right)}\end{array}}\right|\frac{x}{{\prod\nolimits_{k=1}^{K}{{\Omega_{k}}}}}}\right],$$

$$Y_{p,q}^{m,n}\left[{\left.{\begin{array}[]{*{20}{c}}{\left({{a_{1}},{\alpha_{1}},{A_{1}},{\varphi_{1}}}\right),\cdots,\left({{a_{p}},{\alpha_{p}},{A_{p}},{\varphi_{p}}}\right)}\\ {\left({{b_{1}},{\beta_{1}},{B_{1}},{\phi_{1}}}\right),\cdots,\left({{b_{q}},{\beta_{q}},{B_{q}},{\phi_{q}}}\right)}\end{array}}\right|x}\right]=\frac{1}{{2\pi{\rm i}}}\int_{\mathcal{L}}{M_{p,q}^{m,n}\left[s\right]{x^{-s}}ds},$$

$$\begin{gathered} M_{p,q}^{m,n}\left[s\right]=\frac{{\prod }_{j=1}^m{B_j}^{{\varphi }_j+b_j+{\beta }_{j}s-1}\Psi \left({\varphi }_j,{\varphi }_j+b_j+{\beta }_{j}s;B_j\right)}{{\prod }_{i=n+1}^p{A_i}^{{\phi }_i+a_i+{\alpha }_{i}s-1}\Psi \left({\phi }_i,{\phi }_i+a_i+{\alpha }_{i}s;A_i\right)}\times \\ \frac{{\prod }_{i=1}^n{A_i}^{{\phi }_i-a_i-{\alpha }_{i}s}\Psi \left({\phi }_i,{\phi }_i+1-a_i-{\alpha }_{i}s;A_i\right)}{{\prod }_{j=m+1}^q{B_j}^{{\varphi }_j-b_j-{\beta }_{j}s}\Psi \left({\varphi }_j,{\varphi }_j+1-b_j-{\beta }_{j}s;B_j\right)}. \end{gathered}$$

$$f_{G_K}\left(x\right)={\mathrm{E}}_T\left\{f_{G_K|T}\left(x|t\right)\right\}={\int }_0^{\infty }{f}_{G_K|T}\left(x|t\right)f_T\left(t\right)dt,$$

$$f_T\left(t\right)=\frac{t^{m-1}}{\Gamma \left(m\right)}{e}^{-t},t\ge 0.$$

$$f_{G_K}\left(x\right)$$

$$={\sum }_{{\ell }_1,\cdots ,{\ell }_K=0}^{\infty }{W}_{\mathbf{\ell }}{f}_{{\mathcal{A}}_{\mathbf{\ell }}}(x),$$

$$W_{\mathbf{\ell }}\triangleq \frac{\Gamma \left(m+{\sum }_{k=1}^K{\ell }_k\right)}{\Gamma \left(m\right)}{\left(1+\sum _{k=1}^K\frac{{{\lambda }_k}^2}{1-{{\lambda }_k}^2}\right)}^{-m}\prod _{k=1}^K\frac{{w_k}^{{\ell }_k}}{{\ell }_k!},\mathbf{\ell }=({\ell }_1,\cdots ,{\ell }_K)\in {{\mathbb{N}}_0}^K,$$

$$F_{G_K}\left(x\right)={\int }_0^x{f}_{G_K}\left(t\right)dt={\sum }_{{\ell }_1,\cdots ,{\ell }_K=0}^{\infty }{W}_{\mathbf{\ell }}{\int }_0^x{f}_{{\mathcal{A}}_{\mathbf{\ell }}}(t)dt={\sum }_{{\ell }_1,\cdots ,{\ell }_K=0}^{\infty }{W}_{\mathbf{\ell }}{F}_{{\mathcal{A}}_{\mathbf{\ell }}}(x),$$

$${F_{{{{\mathcal{A}}}_{\boldsymbol{\ell}}}}}(x)=Y_{1,K+1}^{K,1}\left[{\left.{\begin{array}[]{*{20}{c}}{\left({1,1,0,1}\right)}\\ {\left({1,1,\frac{1}{{{\Omega_{1}}}},m+{\ell_{1}}}\right),\cdots,\left({1,1,\frac{1}{{{\Omega_{K}}}},m+{\ell_{K}}}\right),\left({0,1,0,1}\right)}\end{array}}\right|\frac{x}{{\prod\nolimits_{k=1}^{K}{{\Omega_{k}}}}}}\right].$$

$$p_{out,K}=F_G\left(2^{\mathcal{R}}\right)={\sum }_{{\ell }_1,\cdots ,{\ell }_K=0}^{\infty }{W}_{\mathbf{\ell }}{F}_{{\mathcal{A}}_{\mathbf{\ell }}}(2^{\mathcal{R}}).$$

$$\left|h_{{\mathcal{A}}_{\mathbf{\ell }},k}\right|\sim \mathrm{Nakagami}\left(m+{\ell }_k,\frac{\left(m+{\ell }_k\right)(1-{{\lambda }_k}^2){{\sigma }_k}^2}{m}\right),1\le k\le K.$$

$$p_{out,K}\approx {\sum }_{{\sum }_{k=1}^K{\ell }_k\le N}{W}_{\mathbf{\ell }}{F}_{{\mathcal{A}}_{\mathbf{\ell }}}\left(2^{\mathcal{R}}\right)\triangleq {\tilde{p}}_{out,K}^N.$$

$$\nabla (N)=p_{out,K}-{\tilde{p}}_{out,K}^N={\sum }_{{\sum }_{k=1}^K{\ell }_k\ge N+1}{W}_{\mathbf{\ell }}{F}_{{\mathcal{A}}_{\mathbf{\ell }}}\left(2^{\mathcal{R}}\right)\ge 0.$$

$$\nabla (N)$$

$$\left(\frac{P_1}{{\mathcal{N}}_0},\frac{P_2}{{\mathcal{N}}_0},\cdots ,\frac{P_K}{{\mathcal{N}}_0}\right)=\frac{P_{total}}{{\mathcal{N}}_0}\left(\frac{P_1}{P_{total}},\frac{P_2}{P_{total}},\cdots ,\frac{P_K}{P_{total}}\right)=\gamma \mathbf{\theta },$$

$$\begin{gathered} F_{{\mathcal{A}}_{\mathbf{\ell }}}(x)={\gamma }^{-d_{{\mathcal{A}}_{\mathbf{\ell }}}}{\prod }_{k=1}^K\frac{1}{\Gamma \left(m+{\ell }_k\right)}{\left(\frac{m}{{\theta }_k{{\sigma }_k}^2\left(1-{{\lambda }_k}^2\right)}\right)}^{m+{\ell }_k} \\ \times {\sum }_{n_1,\cdots ,n_K=0}^{\infty }{\prod }_{k=1}^K\frac{1}{n_k!}{\left(-\frac{m}{\gamma {\theta }_k{{\sigma }_k}^2\left(1-{{\lambda }_k}^2\right)}\right)}^{n_k}{g}_{\mathbf{n}+\mathbf{\ell }}\left(x\right), \end{gathered}$$

$$g_{\mathbf{\ell }}\left(x\right)$$

$$\displaystyle=\prod\limits_{k=1}^{K}{\Gamma\left({m+{\ell_{k}}}\right)}G_{K+1,K+1}^{0,K+1}\left({\left.{\begin{array}[]{*{20}{c}}{1,1+{\ell_{1}}+m,\cdots,1+{\ell_{K}}+m}\\ {1,\cdots,1,0}\end{array}}\right|x}\right),$$

$$F_{{\mathcal{A}}_{\mathbf{\ell }}}(x)={\gamma }^{-d_{{\mathcal{A}}_{\mathbf{\ell }}}}\prod _{k=1}^K\frac{1}{\Gamma \left(m+{\ell }_k\right)}{\left(\frac{m}{{\theta }_k{{\sigma }_k}^2\left(1-{{\lambda }_k}^2\right)}\right)}^{m+{\ell }_k}{g}_{\mathbf{\ell }}\left(x\right)+o\left({\gamma }^{-d_{{\mathcal{A}}_{\mathbf{\ell }}}}\right),$$

$$\frac{F_{{\mathcal{A}}_{\mathbf{\ell }}}\left(x\right)}{F_{{\mathcal{A}}_0}\left(x\right)}=o\left(1\right),\mathbf{\ell }\in {{\mathbb{N}}_0}^K\mathrm{and}\mathbf{\ell }\ne 0.$$

$$\underset{\gamma \to \infty }{\lim }\frac{F_{{\mathcal{A}}_{\mathbf{\ell }}}\left(x\right)}{F_{{\mathcal{A}}_0}\left(x\right)}$$

Full text

I Introduction

Recently, hybrid automatic repeat request (HARQ) technique has found wide applications in wireless communications due to its high potential for reliable transmissions. It has been proved from an information-theoretic view in [1] that HARQ with incremental redundancy (HARQ-IR) can achieve the ergodic capacity in Gaussian collision channels. It also provides superior performance over other types of HARQ since extra coding gain is obtained through code combining. Thus this paper focuses on the analysis of HARQ-IR. As shown in [2], the most fundamental metric to evaluate the performance of HARQ-IR is outage probability and its analysis essentially turns to determine the cumulative distribution function (CDF) of accumulated mutual information.

In prior literature, various methods have been proposed to derive the CDF of accumulated mutual information for HARQ-IR under either quasi-static [3, 4] or fast fading channels [5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. To be more specific, for quasi-static fading channels where the same channel realization is experienced by the transmitted signal in each HARQ round, the CDF of accumulated mutual information is easy to be derived because of the simplicity of handling only one single random variable (RV). Hereby, [3] and [4] have conducted average rate analysis and power optimization for HARQ-IR, respectively. Unfortunately, the analytical results under quasi-static fading channels are only applicable to low mobility environment. In high mobility environment, the transmitted signals among all HARQ rounds would experience fast fading channels, where the channel responses vary independently from one transmission to another. Under fast fading channels, the derivation of the CDF of accumulated mutual information essentially turns to determine the distribution of the product of multiple independent shifted RVs. There are several approaches proposed to tackle this problem in the literature. For example, in [5, 6, 7], Log-normal approximation is proposed based on central limit theorem (CLT). In [8, 9, 10], a lower bound and an upper bound of the CDF are derived based on Jensen’s inequality and Minkowski inequality, respectively. To calculate the exact CDF, Mellin transform and multi-fold convolution have been applied in [11, 12] and [13] respectively. Specifically, based on Mellin transform, [11] derives the exact CDF in terms of generalized Fox’s H function. The analytical results are then applied for outage analysis in [12]. In [13], noticing that the accumulated mutual information is a sum of independent RVs under fast fading channels, the exact CDF is thus derived by using multi-fold convolution. Moreover, by associating the analytical results in [13] with matrix exponential distribution through Laplace transform, [14] approximates the CDF of accumulated mutual information to a matrix exponential distribution. Unfortunately, the exact results in [11, 12, 13] are too complicated to provide meaningful insights due to difficulties in handling generalized Fox’s H function and multi-fold integrals, while the approximated results in [5, 6, 7] and [14] involve certain numerical calculations which also hinder the extraction of clear insights.

The analysis of HARQ-IR becomes more challenging while considering time-correlated fading channels, which usually occur in low-to-medium mobility environment, because a product of multiple shifted correlated RVs is involved in the CDF of accumulated mutual information. To the best of our knowledge, there are only few approximation approaches available to analyze the performance of HARQ-IR over time-correlated fading channels, that is, Log-normal approximation [15], polynomial fitting technique [16] and inverse moment matching method [17]. Unfortunately, the Log-normal approximation in [15] is inaccurate when the fading channels have medium-to-high time correlation. Although the results in [16] and [17] usually can achieve a good approximation of the CDF of accumulated mutual information, their tightness is proved in mean square error (MSE) sense. It does not necessarily imply that the approximation error always approaches to zero and the occurrence of inaccuracy is found especially under low outage probability/hign signal-to-noise ratio (SNR). In addition, the approximated results in [16, 17] are still complicated with little insights, thus limiting their applications in practical system design.

In this paper, asymptotic outage analysis is conducted to thoroughly investigate the performance of HARQ-IR over time-correlated Nakagami-$m$ fading channels with meaningful insights. Based on our finding that the received SNRs in multiple HARQ transmissions are conditionally independent given a certain RV, the CDF of accumulated mutual information can be exactly derived by using conditional Mellin transform. With the result, the outage probability can be derived as a weighted sum of outage probabilities of HARQ-IR over independent Nakagami fading channels, where the weightings are determined by a negative multinomial distribution. A truncation approximation with uniform convergence is then proposed to ease the computation of the exact outage probability. Notice that the uniform convergence is stronger than the MSE convergence in [16, 17] and thus offers higher accuracy in the approximation. The analytical result in this paper also enables the asymptotic outage analysis, which has never been conducted for HARQ-IR even under fast fading channels. More specifically, the asymptotic outage probability is derived in a simple form which clearly quantifies the impacts of transmission powers in multiple transmissions, channel time correlation and information transmission rate. It is also proved that the asymptotic outage probability is an inverse power function of the product of transmission powers in all HARQ rounds, an increasing function of the channel time correlation coefficients, and a monotonically increasing and convex function of information transmission rate. The simple expression of the asymptotic result enables optimal power allocation and optimal rate selection of HARQ-IR through plenty of readily available optimization techniques with low complexity. Finally, numerical results are provided to verify our analytical results and justify the application of the asymptotic result for optimal system design.

The remainder of this paper is organized as follows. Section II introduces the system model and formulates the outage probability. The probability distribution of the product of multiple correlated shifted RVs which is necessary for outage analysis is derived and the exact outage probability is given in Section III. Asymptotic outage probability is then derived in a simple form and meaningful insights are shown in Section IV. In Section V, numerical results are presented for validation and discussion. Finally, Section VI concludes this paper.

II System Model and Outage Formulation

A point-to-point HARQ-IR enabled system operating over time-correlated Nakagami-$m$ fading channels is considered in this paper. The details of the system are introduced as follows.

II-A HARQ-IR Protocol

Following HARQ-IR protocol, the source first encodes every $b$-bits information message into $K$ codewords each with length of $L$, where $K$ denotes the maximum allowable number of transmissions for each message. As a result, the initial coding rate is $c=b/L$ [18]. The $K$ codewords will be sequentially transmitted to the destination until the message is successfully decoded. In each transmission, the previously received codewords are combined with the most recently received codeword for joint decoding. If successful, an acknowledgement (ACK) message is fed back from the destination to the source and the source then moves to the transmission of the next information message. Otherwise, a negative acknowledgement (NACK) message is fed back from the destination to the source and the source transmits the next codeword until the maximum number of transmissions $K$ is reached. Similarly to the analyses in the literature [12, 13, 14], error-free feedback channels are assumed here, that is, all feedback messages can be successfully decoded.

II-B Time-Correlated Nakagami-$m$ Fading Channels

Denote ${\bf x}_{k}$ as the $k$th codeword with length of $L$. It is transmitted over a block fading channel, i.e., each symbol of ${\bf x}_{k}$ experiences an identical channel realization during the $k$th transmission. Therefore, the signal received in the $k$th transmission is given by

[TABLE]

where ${\bf n}_{k}$ denotes a complex additive white Gaussian noise (AWGN) vector with zero mean vector and covariance matrix $\mathcal{N}_{0}{\bf I}_{L}$, i.e., ${\bf n}_{k}\sim{\cal CN}(0,\mathcal{N}_{0}{{\bf I}_{L}})$, ${\bf I}_{L}$ represents an $L\times L$ identity matrix, and $h_{k}$ denotes the block fading channel coefficient in the $k$th transmission.

Notice that Nakagami-$m$ fading is a general channel model whose parameters can be adjusted to fit a variety of empirical measurements [19], e.g., it covers one-sided Gaussian channel and Rayleigh fading as special cases by setting $m=\frac{1}{2}$ and $m=1$ respectively. It is thus considered here. Unlike most of the prior analyses, time-correlated fading channels are considered, that is, the channel coefficients among the $K$ transmissions are correlated. The time correlation usually occurs in low-to-medium mobility environment [20]. Herein, a widely used Nakagami-$m$ fading channel model with generalized correlation is adopted and the channel magnitude $|h_{k}|$ is written as [21, 22, 23]

[TABLE]

where $\left\|\cdot\right\|$ represents ${\ell}^{2}$ norm, ${\vartheta_{k,l}}$ and ${\vartheta_{0,l}}$ are independent circularly-symmetric complex Gaussian random variables (RVs) with zero mean and unit variance, i.e., ${\vartheta_{0,l}},{\vartheta_{k,l}}\sim\mathcal{CN}\left({0,1}\right)$, $m$ represents the fading order that indicates the severity of fading channels, ${\sigma_{k}}^{2}$ denotes Nakagami spread and is equal to the expectation of the squared channel magnitude, i.e., ${\rm E}\{|h_{k}|^{2}\}={\sigma_{k}}^{2}$, and ${\boldsymbol{\lambda}}=\left({{\lambda_{1}},{\lambda_{2}},\cdots,{\lambda_{K}}}\right)$ specifies the generalized time correlation among fading channels of all HARQ transmissions [21]. Under this model, the magnitude of channel coefficient ${\left|{{h_{k}}}\right|}$ follows a Nakagami-$m$ distribution, i.e., $\left|{{h_{k}}}\right|\sim{\rm Nakagami}(m,{\sigma_{k}}^{2})$, and the channel magnitudes ${\left|{{h_{l}}}\right|}$ and ${\left|{{h_{k}}}\right|}$ are correlated. Moreover, the cross correlation coefficient between the squared channel magnitudes ${\left|{{h_{l}}}\right|^{2}}$ and ${\left|{{h_{k}}}\right|^{2}}$ is determined by the time correlation coefficients $\boldsymbol{\lambda}$ as [21]

[TABLE]

where ${\rm Var}(\cdot)$ denotes the variance of the random variable in the brackets. Without loss of generality, the time correlation coefficients are assumed to be non-negative, i.e., $0\leq\lambda_{1},\cdots,\lambda_{K}\leq 1$. Clearly from (2), when ${\boldsymbol{\lambda}}={\bf 1}_{K}$, the time-correlated fading channel reduces to quasi-static fading channel with $h_{1}=h_{2}=\cdots=h_{K}$, while when ${\boldsymbol{\lambda}}={\bf 0}_{K}$, the time-correlated fading channel reduces to a fast fading channel where the channel coefficients $h_{1}$, $h_{2}$, $\cdots$, $h_{K}$ are mutually independent. In other words, the Nakagami-$m$ fading channel with generalized time correlation in (2) includes quasi-static fading channel and fast fading channel as its special cases. Here ${\bf 1}_{K}$ and ${\bf 0}_{K}$ denote an all-ones vector and a null vector, each with length $K$, respectively. Unless otherwise indicated, subscript $K$ is omitted in the sequel.

According to (1), the received signal-to-noise ratio (SNR) in the $k$th transmission is given by

[TABLE]

where $P_{k}$ denotes the transmitted signal power in the $k$th transmission. Since $\left|{{h_{k}}}\right|$ is Nakagami-$m$ distributed, i.e., $\left|{{h_{k}}}\right|\sim{\rm Nakagami}(m,{\sigma_{k}}^{2})$, it is readily proved that ${\gamma_{k}}$ complies with Gamma distribution, i.e., ${\gamma_{k}}\sim\mathcal{G}\left({m,\frac{{P_{k}}{\sigma_{k}}^{2}}{m\mathcal{N}_{0}}}\right)$. Due to the time correlation among the channel coefficients as given in (3), the SNRs $\boldsymbol{\gamma}=(\gamma_{1},\gamma_{2},\cdots,\gamma_{K})$ are correlated Gamma RVs.

II-C Outage Formulation

Outage probability has been proved as the most fundamental performance metric of HARQ schemes [1]. For HARQ-IR, the outage probability is directly determined by the CDF of accumulated mutual information. Specifically, assuming that information-theoretic capacity achieving channel coding is adopted for HARQ-IR, an outage event happens when the information cannot be successfully decoded after $K$ transmissions, i.e., the accumulated mutual information $I_{K}$ is below the information transmission rate $\mathcal{R}$. Notice that the information transmission rate $\mathcal{R}$ depends on the coding rate $c$, modulation scheme and symbol transmission rate. The outage probability after $K$ transmissions is thus written as

[TABLE]

where $F_{{I_{K}}}(\cdot)$ denotes the CDF of ${I_{K}}$. In HARQ-IR protocol, decoding at each transmission is performed based on the combined codewords received in the previously transmissions and the current transmission and the accumulated mutual information after $K$ transmissions is given by

[TABLE]

Accordingly, the outage probability ${p_{out,K}}$ becomes

[TABLE]

where $F_{G_{K}}(\cdot)$ denotes the CDF of $G_{K}$. Noticing that $(\gamma_{1},\gamma_{2},\cdots,\gamma_{K})$ are correlated Gamma RVs, the derivation of outage probability ${p_{out,K}}$ essentially turns to determining the CDF of the product of multiple correlated shifted-Gamma RVs, i.e., ${F_{G_{K}}}\left({{x}}\right)$. Due to the presence of time correlation and Gamma distribution, the outage analysis is more challenging than those for quasi-static [3, 4] or fast fading channels [5, 6, 7, 8, 9, 10, 24, 11, 12, 13, 14].

It is worth noting that the CDF of the product of multiple shifted and correlated RVs has numerous applications in wireless communications and the outage formulation in (5) is also applicable to parallel transmission systems and orthogonal frequency-division multiplexing (OFDM) systems [25, 26]. Specifically, in OFDM systems, coded signals are transmitted through multiple subcarriers. Due to close frequency spacing of the subcarriers and limited Doppler spread, the fading channels on multiple subcarriers are usually correlated [27]. The outage analysis in OFDM systems is thus also reduced to the analysis of the CDF of the product of multiple shifted and correlated RVs in (7).

III Probability Distribution of the Product of Multiple Correlated Shifted-Gamma RVs

As aforementioned, the derivation of outage probability ${p_{out,K}}$ essentially turns to determining the CDF of the product of multiple shifted and correlated Gamma RVs. It is mathematically difficult due to the presence of time correlation. In the literature, only a few approximated results are available in [15, 16, 17]. They are complicated without clear insights and thus are not favorable for system design. Moreover, their accuracy can not be guaranteed under certain scenarios. In this paper, we will derive the exact CDF of the product of multiple shifted and correlated Gamma RVs, based on which a truncation approximation with uniform convergence can be proposed and asymptotic outage probability can be derived in a simple form with clear insights.

III-A Exact CDF

From (2) and (4), we can find that the SNRs $\{\gamma_{1},\gamma_{2},\cdots,\gamma_{K}\}$ are conditionally independent, given the complex Gaussian RVs ${\vartheta_{0,1}},\cdots,{\vartheta_{0,m}}$. It has been proved in [28, Theorem 1.3.4] that $\{\gamma_{1},\gamma_{2},\cdots,\gamma_{K}\}$ follow independent noncentral chi-squared distributions with $2m$ degrees-of-freedom when conditioned on $T\triangleq\sum\nolimits_{l=1}^{m}{{{\left|{\vartheta_{0,l}}\right|}^{2}}}$. The conditional probability density function (PDF) is given in the following lemma.

Lemma 1**.**

[28, Theorem 1.3.4]** The conditional PDF of the SNR $\gamma_{k}$ given $T=t$ is written as

[TABLE]

where ${u_{k}}=\frac{{P_{k}}{\sigma_{k}}^{2}}{m\mathcal{N}_{0}}$, $\Omega_{k}=\frac{{P_{k}}{\sigma_{k}}^{2}\left({1-{\lambda_{k}}^{2}}\right)}{m\mathcal{N}_{0}}$, $\Gamma(\cdot)$ and ${}_{0}F_{1}(\cdot)$ denote Gamma function and the confluent hypergeometric limit function [29, Eq.16.2.1], respectively.

It should be mentioned that (8) is not applicable to quasi-static fading channels, i.e., $\boldsymbol{\lambda}=\bf 1$, because of the occurrence of singularity. Due to its speciality, the outage analysis of HARQ-IR over quasi-static fading channels will be discussed separately.

Noticing that Mellin transform111 The Mellin transform with respect to a function $f(x)$ is defined as $\left\{{\mathcal{M}f}\right\}\left(s\right)=\int\nolimits_{0}^{\infty}{{x^{s-1}}f\left(x\right)dx}\triangleq\tilde{f}\left(s\right)$, and the associated inverse Mellin transform is given by $f\left(x\right)=\frac{1}{{2\pi{\rm{i}}}}\int\nolimits_{c-{\rm{i}}\infty}^{c+{\rm{i}}\infty}{{x^{-s}}\tilde{f}\left(s\right)dx}$ [30]. is a powerful mathematical tool to obtain the probability distribution of a product of multiple independent RVs [24, 11, 12], the conditional independence of SNRs $\boldsymbol{\gamma}$ given $T$ inspires us to derive the conditional PDF of $G_{K}$ in closed-form by using conditional Mellin transform. Specifically, with the conditional independence, the Mellin transform of the conditional PDF of $G_{K}$ can be written as a product of the Mellin transforms corresponding to the conditional PDFs of the shifted SNRs. With this special property, as proved in Appendix A, the conditional PDF of $G_{K}$ given $T=t$, ${f_{\left.G_{K}\right|T}}(x|t)$, can be derived as

[TABLE]

where ${\boldsymbol{\ell}}$ is a vector with $K$ elements as ${\boldsymbol{\ell}}=(\ell_{1},\cdots,\ell_{K})$ and ${{f_{{{{\cal A}}_{\boldsymbol{\ell}}}}}(x)}$ is given by

[TABLE]

with ${\rm i}=\sqrt{-1}$, $\Psi\left({\alpha,\gamma;z}\right)=\frac{1}{{\Gamma\left(\alpha\right)}}\int\nolimits_{0}^{\infty}{{e^{-zt}}{t^{\alpha-1}}{{\left({1+t}\right)}^{\gamma-\alpha-1}}dt}$ denoting Tricomi’s confluent hypergeometric function [31, Eq.9.211.4], and $Y_{p,q}^{m,n}[\cdot]$ denoting the generalized Fox’s H function defined by Mellin-Barnes integral as [11, 12]

[TABLE]

where $\mathcal{L}$ is a Mellin-Barnes contour in the complex $s$-plane running from $c-{\rm i}\infty$ to $c+{\rm i}\infty$, $c\in\mathbb{R}$ and $M_{p,q}^{m,n}\left[s\right]$ is written as

[TABLE]

Notice that the generalized Fox’s H function does have two special properties as shown in Appendix B which can simplify our mathematical derivations, and an efficient MATHEMATICA® implementation of the generalized Fox’s H function in (13) can be found in [11]. Interestingly from (III-A), $f_{{\mathcal{A}}_{\boldsymbol{\ell}}}(x)$ is the inverse Mellin transform of a product function $\prod\nolimits_{k=1}^{K}{\frac{{\Psi\left({m+{\ell_{k}},s+m+{\ell_{k}};\frac{1}{{{\Omega_{k}}}}}\right)}}{{{{\left({{\Omega_{k}}}\right)}^{m+{\ell_{k}}}}}}}$. Here the multiplier ${\frac{{\Psi\left({m+{\ell_{k}},s+m+{\ell_{k}};\frac{1}{{{\Omega_{k}}}}}\right)}}{{{{\left({{\Omega_{k}}}\right)}^{m+{\ell_{k}}}}}}}$ in fact is the Mellin transform of the PDF of a shifted-Gamma RV $(1+R_{{\boldsymbol{\ell}},k})$, where $R_{{\boldsymbol{\ell}},k}\sim\mathcal{G}(m+\ell_{k},{\Omega_{k}})$. Therefore, ${f_{{{{\mathcal{A}}}_{\boldsymbol{\ell}}}}}(x)$ can be regarded as the PDF of the product of $K$ independent shifted-Gamma RVs i.e., $\mathcal{A}_{\boldsymbol{\ell}}=\prod\nolimits_{k=1}^{K}(1+R_{{\boldsymbol{\ell}},k})$.

With (9), the PDF of $G_{K}$ can be eventually obtained by integrating the conditional probability ${f_{\left.G_{K}\right|T}}\left({\left.x\right|t}\right)$ over the distribution of $T$, such that

[TABLE]

where ${f_{T}}\left(t\right)$ denotes the PDF of $T$. With respect to the distribution of $T\triangleq\sum\nolimits_{l=1}^{m}{{{\left|{\vartheta_{0,l}}\right|}^{2}}}$, since $\vartheta_{0,1},\cdots,\vartheta_{0,m}$ are independent and identically distributed (i.i.d.) complex Gaussian RVs with zero mean and unit variance, it can be proved that $\left|\vartheta_{0,1}\right|^{2},\cdots,\left|\vartheta_{0,m}\right|^{2}$ are i.i.d. Gamma RVs, i.e., ${{\left|{\vartheta_{0,l}}\right|}^{2}}\sim\mathcal{G}\left({1,1}\right)$. As a sum of $m$ independent Gamma RVs with identical scale parameter, the random variable $T$ then follows a Gamma distribution with the PDF of [32, p289]

[TABLE]

Plugging (9) and (16) into (15), it follows by using [31, Eq.3.381.4] that

[TABLE]

where ${{\boldsymbol{\ell}}}$ is a vector with $K$ elements as ${{\boldsymbol{\ell}}}=(\ell_{1},\cdots,\ell_{K})$ and the weighting ${W_{\boldsymbol{\ell}}}$ is given as

[TABLE]

${w_{k}}=\frac{{{\lambda_{k}}^{2}}}{{1-{\lambda_{k}}^{2}}}{\left({1+\sum\nolimits_{l=1}^{K}{\frac{{{\lambda_{l}}^{2}}}{{1-{\lambda_{l}}^{2}}}}}\right)^{-1}}$ and ${{\mathbb{N}}_{0}}^{K}$ denotes $K$-ary Cartesian power of natural number set. By comparing (18) with [33, Eq.13.8-1], it is found that $\left\{{W_{\boldsymbol{\ell}}},{{{\boldsymbol{\ell}}}\in{{\mathbb{N}}_{0}}^{K}}\right\}$ are probabilities of a negative multinomial distributed vector RV ${\boldsymbol{\ell}}$, i.e., ${\boldsymbol{\ell}}\sim{\rm NM}(m,\bf w)$, where ${\bf w}=(w_{1},\cdots,w_{K})$. Clearly, we have $\sum\nolimits_{{\boldsymbol{\ell}}\in{{\mathbb{N}}_{0}}^{K}}{{W_{\boldsymbol{\ell}}}}=1$. As such, ${f_{G_{K}}}\left(x\right)$ is expressed as a weighted sum of the PDFs corresponding to $\left\{\mathcal{A}_{\boldsymbol{\ell}},{{{\boldsymbol{\ell}}}\in{{\mathbb{N}}_{0}}^{K}}\right\}$.

Based on (III-A), the CDF of the product of multiple correlated shifted-Gamma RVs ${F_{G_{K}}}\left(x\right)$ can be derived as shown in the following theorem.

Theorem 1**.**

The CDF of $G_{K}=\prod\nolimits_{k=1}^{K}{\left({1+{\gamma_{k}}}\right)}$ is given by

[TABLE]

where ${{F_{{{{\mathcal{A}}}_{\boldsymbol{\ell}}}}}(x)}$ denotes the CDF of the product of independent shifted-Gamma RVs $\mathcal{A}_{\boldsymbol{\ell}}=\prod\nolimits_{k=1}^{K}(1+R_{{\boldsymbol{\ell}},k})$ with $R_{{\boldsymbol{\ell}},k}\sim\mathcal{G}(m+\ell_{k},{\Omega_{k}})$ and ${{\boldsymbol{\ell}}}=(\ell_{1},\cdots,\ell_{K})$. Specifically, ${{F_{{{{\mathcal{A}}}_{\boldsymbol{\ell}}}}}(x)}$ can be expressed in terms of the generalized Fox’s H function as

[TABLE]

Proof.

Please see Appendix B.

∎

Therefore, the outage probability ${p_{out,K}}$ can be obtained by substituting (19) into (7), such that

[TABLE]

Notice that $\mathcal{A}_{\boldsymbol{\ell}}=\prod\nolimits_{k=1}^{K}(1+R_{{\boldsymbol{\ell}},k})$, $R_{{\boldsymbol{\ell}},k}$ follows Gamma distribution as $R_{{\boldsymbol{\ell}},k}\sim\mathcal{G}(m+\ell_{k},{\Omega_{k}})$ and $\Omega_{k}=\frac{{P_{k}}{\sigma_{k}}^{2}\left({1-{\lambda_{k}}^{2}}\right)}{m\mathcal{N}_{0}}$. Clearly, the random variable $R_{{\boldsymbol{\ell}},k}$ can be factorized as $R_{{\boldsymbol{\ell}},k}=\frac{P_{k}\left|{{h_{{{\mathcal{A}}_{\boldsymbol{\ell}}},k}}}\right|^{2}}{\mathcal{N}_{0}}$ where $\left|{{h_{{{\mathcal{A}}_{\boldsymbol{\ell}}},k}}}\right|$ follows Nakagami distribution as

[TABLE]

Therefore, ${F_{{{{\cal A}}_{\boldsymbol{\ell}}}}}\left(2^{\mathcal{R}}\right)$ can be regarded as the outage probability of HARQ-IR after $K$ transmissions over independent Nakagami-$m$ fading channels, where $R_{{\boldsymbol{\ell}},k}$ denotes the received SNR in the $k$th HARQ transmission and ${{h_{{{\mathcal{A}}_{\boldsymbol{\ell}}},k}}}$ denotes the Nakagami-$m$ fading channel coefficient in the $k$th transmission. Consequently, the outage probability ${p_{out,K}}$ can be rephrased as a weighted sum of outage probabilities of HARQ-IR over independent Nakagami fading channels where the weightings ${W_{\boldsymbol{\ell}}}$ are probabilities of the negative multinomial distribution ${\rm NM}(m,\bf w)$.

It is worth mentioning that the outage result in (21) is applicable to fast fading channels. Under fast fading channels, the time correlation coefficients are equal to zero, i.e., $\boldsymbol{\lambda}=\bf 0$. Putting it into (18), we have $W_{\bf{0}}=1$ and $W_{\boldsymbol{\ell}}=0$ for all ${\boldsymbol{\ell}}\neq{\bf 0}$. Then the outage probability ${p_{out,K}}$ reduces to ${F_{{{\mathcal{A}}_{\bf 0}}}}({2^{R}})$.

III-B Truncation Approximation

Considering that $p_{out,K}$ is represented by an infinite series in (21), it is impossible to compute the exact value by adding an infinite number of terms up. To enable its computation, it is natural to truncate $p_{out,K}$ into a finite series. Towards this end, an effective truncation approach is proposed. Specifically, $p_{out,K}$ in (21) is approximated with a truncation order $N$ as

[TABLE]

The truncation error $\nabla(N)$ is characterized by the difference between $p_{out,K}$ and $\tilde{p}_{out,K}^{N}$, such that

[TABLE]

Clearly from (24), $\nabla(N)$ is a monotonically decreasing function of $N$ and satisfies $\lim_{N\rightarrow\infty}\nabla(N)=0$. Moreover, since ${{F_{{{\cal A}_{\boldsymbol{\ell}}}}}\left(2^{\cal R}\right)}\leq 1$ and ${{W_{\boldsymbol{\ell}}}}\leq 1$, the truncation error $\nabla(N)$ is in fact uniformly upper bounded by ${\sum\nolimits_{{\sum\nolimits_{k=1}^{K}\ell_{k}\geq N+1}}{{W_{\boldsymbol{\ell}}}}}$ which is irrespective of the rate $\cal R$ and converges to zero when $N$ increases. Thus strictly speaking, the truncation approximation in (23) actually admits a uniform convergence[34, p147]222We say that $f_{n}(x)$ converges uniformly to $f(x)$ on its domain $\cal D$ if, given $\varepsilon$, there exists an integer $N(\varepsilon)$ independent of $x$, such that $|f_{n}(x)-f(x)|\leq\varepsilon$ for all $x\in\cal D$ whenever $n\geq N(\varepsilon)$[34, p147].. This uniform convergence guarantees the truncation approximation of outage probability with high accuracy for any $\cal R$. It is different from the convergence in MSE in [16, 17] which does not necessarily imply that the approximation error approaches to zero for arbitrary $\cal R$ [35, p86]. Therefore our truncation approximation is expected to perform better than those in [16, 17], which will be further demonstrated in Section V.

As proved in Appendix C, we further notice that the truncation error $\nabla(N)$ is bounded as

[TABLE]

where ${{F_{{{\cal A}_{\boldsymbol{\ell}}},N}^{\rm max}}\left(2^{\cal R}\right)}=\mathop{\max}\nolimits_{\sum\nolimits_{k=1}^{K}\ell_{k}=N+1}\left({{F_{{{\cal A}_{\boldsymbol{\ell}}}}}\left(2^{\cal R}\right)}\right)$, ${\xi\left(N\right)}$ is a decreasing function of $N$ as $\xi\left(N\right)={\left({\sum\nolimits_{k=1}^{K}{{w_{k}}}}\right)^{N+1}}\frac{{{{\left(m\right)}_{N+1}}}}{{\left({N+1}\right)!}}{}_{2}{F_{1}}\left({m+N+1,1;N+2;\sum\nolimits_{k=1}^{K}{{w_{k}}}}\right)$, the notation $(\cdot)_{n}$ stands for Pochhammer symbol, and ${}_{2}F_{1}(a,b;c;z)=\sum\nolimits_{s=0}^{\infty}{\frac{{{{\left(a\right)}_{s}}{{\left(b\right)}_{s}}}}{{\Gamma\left({c+s}\right)s!}}{z^{s}}}$ denotes hypergeometric function [29, Eq.15.1.1]. For Rayleigh fading channels, i.e., $m=1$, the term $\xi\left(N\right)$ reduces to $\xi\left(N\right)={W_{\bf{0}}}^{-1}{\left({\sum\nolimits_{k=1}^{K}{{w_{k}}}}\right)^{N+1}}$ by using ${}_{2}{F_{1}}\left({N+2,1;N+2;\sum\nolimits_{k=1}^{K}{{w_{k}}}}\right)={W_{\bf{0}}}^{-1}$ in [36, Eq.1.40]. Together with ${{F_{{{\cal A}_{\boldsymbol{\ell}}},N}^{\rm max}}\left(2^{\cal R}\right)}\leq 1$, the truncation error under Rayleigh fading channels is thus further bounded as $\nabla(N)\leq\mathcal{B}_{u}\leq{\left({\sum\nolimits_{k=1}^{K}{{w_{k}}}}\right)^{N+1}}$. Since ${\sum\nolimits_{k=1}^{K}{{w_{k}}}}<1$, the truncation error under Rayleigh fading channels exponentially decays with $N$. It implies that the convergence speed of our truncation approach is fast and also justifies the effectiveness of our truncation approach.

In practice, the truncation order plays an important role in improving the approximation accuracy and an efficient selection of the truncation order is necessary. Specifically, for a given accuracy requirement $\varepsilon$, the minimal truncation order $N(\varepsilon)$ can be determined by setting the upper bound $\mathcal{B}_{u}$ in (25) less than the maximum tolerable approximation error $\varepsilon$, such that $N(\varepsilon)=\mathop{\min}\left\{\left.N\geq 0\right|\mathcal{B}_{u}\leq\varepsilon\right\}$.

The impact of truncation order on the approximation accuracy is then examined and the results are shown in Fig. 1. Clearly, the truncated outage probability $\tilde{p}_{out,K}^{N}$ converges quite fast to the exact outage probability $p_{out,K}$ and the truncation order of $N=3$ is sufficient for a very good approximation of the outage probability. Plugging $N=3$ together with ${\cal R}=1$bps/Hz into (25), the upper bounds $\mathcal{B}_{u}$ of four cases from I to IV can be calculated as $7.0*10^{-5}$, $7.4*10^{-6}$, $1.4*10^{-4}$ and $3.7*10^{-8}$, respectively. They are very small and thus verify the accuracy of our truncation approach. By comparing the upper bound $\mathcal{B}_{u}$ of Case I with that of Case II, it can be observed that the increase of $K$ decreases $\mathcal{B}_{u}$ from $7.0*10^{-5}$ to $7.4*10^{-6}$, which means that the increase of $K$ would be beneficial to the reduction of the truncation error and thus improve the accuracy of our truncation approach. This can be roughly explained as follows. From the upper bound $\mathcal{B}_{u}$ in (25), we can see that increasing $K$ has twofold impact on the upper bound $\mathcal{B}_{u}$. Specifically, the increase of $K$ reduces the first component ${W_{\bf{0}}}{{F_{{{\cal A}_{\boldsymbol{\ell}}},N}^{\rm max}}\left(2^{\cal R}\right)}$ of $\mathcal{B}_{u}$ on one hand, while increases the second component $\xi\left(N\right)$ on the other hand. Since the first component dominates particularly in a low-to-medium outage region which is the most concerned region for practical applications, the upper bound $\mathcal{B}_{u}$ can be viewed as a decreasing function of $K$.

IV Asymptotic Outage Analysis

IV-A Asymptotic Outage Probability

Although the outage probability of HARQ-IR over time-correlated Nakagami fading channels can be exactly derived in closed-form as shown in (21) and can be approximated as (23), they are still complex and hinder the extraction of meaningful insights. Fortunately, after analyzing the CDF ${F_{{{{\cal A}}_{\boldsymbol{\ell}}}}}\left(x\right)$, we can find some special properties of ${F_{{{{\cal A}}_{\boldsymbol{\ell}}}}}\left(x\right)$, which enable us to further conduct asymptotic outage analysis so that the expression of outage probability can be simplified with clear insights. To proceed with the analysis, we rewrite the transmission SNRs as

[TABLE]

where ${P_{total}}={P_{1}}+\cdots+{P_{K}}$, $\gamma\triangleq\frac{{{P_{total}}}}{{{{\cal N}_{0}}}}$ denotes the total transmission SNR and $\boldsymbol{\theta}=({\theta_{1}},\cdots,{\theta_{K}})\triangleq\left({\frac{{{P_{1}}}}{{{P_{total}}}},\cdots,\frac{{{P_{K}}}}{{{P_{total}}}}}\right)$ represents the power allocation vector. With this definition of the transmission SNR, the CDF ${F_{{{{\cal A}}_{\boldsymbol{\ell}}}}}\left(x\right)$ in (20) can be rewritten as shown in the following theorem.

Theorem 2**.**

The CDF ${F_{{{{\mathcal{A}}}_{\boldsymbol{\ell}}}}}(x)$ can be rewritten as

[TABLE]

where $d_{{{\cal A}_{\boldsymbol{\ell}}}}=mK+\sum\nolimits_{k=1}^{K}{{{\ell_{k}}}}$, ${\bf n}=(n_{1},\cdots,n_{K})$, ${{\boldsymbol{\ell}}}=(\ell_{1},\cdots,\ell_{K})$ and ${g_{{\boldsymbol{\ell}}}}\left(x\right)$ is given by

[TABLE]

*with $G_{p,q}^{m,n}\left({\cdot|x}\right)$ denoting Meijer-G function [31, Eq.9.301]. *

Proof.

Please see Appendix D. ∎

From Theorem 2, when the transmission SNR is high, i.e., $\gamma\to\infty$, the CDF ${F_{{{{\mathcal{A}}}_{\boldsymbol{\ell}}}}}(x)$ can be further expressed as

[TABLE]

where $o(\cdot)$ refers to the little-O notation, and $f(\gamma)\in o(\phi(\gamma))$ provided that $\lim\limits_{\gamma\to\infty}f(\gamma)/\phi(\gamma)=0$. Based on (32), the following property of the CDF ${F_{{{{\mathcal{A}}}_{\boldsymbol{\ell}}}}}(x)$ holds since $d_{{{\cal A}_{\boldsymbol{\ell}}}}=mK+\sum\nolimits_{k=1}^{K}{{{\ell_{k}}}}>d_{{{\cal A}_{\bf{0}}}}=mK$ for ${\boldsymbol{\ell}\neq\bf{0}}$.

Lemma 2**.**

As $\gamma\to\infty$, the ratio of ${{F_{{{\mathcal{A}}_{\boldsymbol{\ell}}}}}\left(x\right)}$ to ${{F_{{{\mathcal{A}}_{\bf{0}}}}}\left(x\right)}$ satisfies

[TABLE]

Proof.

From (32), we have

[TABLE]

where the last equality holds because of ${d_{{{\cal A}_{\bf 0}}}-{d_{{{\cal A}_{\boldsymbol{\ell}}}}}}<0$ when $\boldsymbol{\ell}\neq\bf 0$. Then by using the little-O notation, (33) directly follows.∎

By applying Lemma 2 into (21), as $\gamma\to\infty$, the outage probability ${p_{out,K}}$ can be rewritten as

[TABLE]

where the last equality holds since $W_{\boldsymbol{\ell}}$ is irrelevant to the transmission SNR $\gamma$ and $\sum\nolimits_{{\ell_{1}}+\cdots+{\ell_{K}}>0}{{W_{\boldsymbol{\ell}}}}=1-W_{\bf 0}<1$. Now putting (32) into (IV-A) and neglecting the high order terms $o\left(1\right)$ and $o\left({{\gamma^{-{d_{{\mathcal{A}_{\boldsymbol{\ell}}}}}}}}\right)$ when $\gamma\to\infty$, the outage probability can be asymptotically approximated as

[TABLE]

Substituting (18) into (36), under high SNR regime, i.e., as $\gamma\to\infty$, the asymptotic outage probability ${p_{out\_{asy},K}}$ can be factorized as

[TABLE]

where $\ell\left({{\boldsymbol{\lambda}},K}\right)$ is explicitly given as

[TABLE]

$\zeta(\boldsymbol{\theta})$ and $\varrho(\boldsymbol{\lambda})$ quantify the impacts of transmission power allocation and channel time correlation and thus are regarded as power allocation impact factor and time correlation impact factor, respectively, and $C(\mathcal{R})$ is termed as coding and modulation gain since it depends on the information transmission rate $\mathcal{R}$ which eventually is determined by the coding rate, modulation scheme and symbol transmission rate.

IV-B Discussions

Clearly from (37), the outage performance is determined by the number of transmissions, transmission power allocation, channel time correlation and information transmission rate. Their impacts will be thoroughly investigated through the analysis of the terms $\zeta(\boldsymbol{\theta})$, $\varrho(\boldsymbol{\lambda})$ and $C(\mathcal{R})$. In addition, outage probability of HARQ-IR under quasi-static fading channels is particularly discussed as a special case of time-correlated fading channels.

IV-B1 Diversity Order

Diversity order indicates the number of degrees of freedom in communication systems. Roughly speaking, it is equivalent to the number of independently faded paths that a transmitted signal experiences. Specifically, it is defined as the slope of the outage probability against the transmission SNR on a log-log scale as [37]

[TABLE]

Putting (37) into (39), the diversity order of HARQ-IR over time-correlated Nakagami-$m$ fading channels directly follows as $d=mK$. Noticing that a Nakagami-$m$ fading channel can be viewed as a set of $m$ parallel independent Rayleigh fading channels, the maximum number of independent fading channels in a HARQ-IR system with a maximum number of $K$ transmissions is $mK$. In other words, channel time correlation would not degrade the diversity order and full diversity can be achieved by HARQ-IR even under time-correlated fading channels. This result is consistent with those in [16, 17], and demonstrates the correctness of our asymptotic outage analysis.

IV-B2 Power Allocation Impact Factor $\zeta(\boldsymbol{\theta})$

It is clear from (37) that power allocation impact factor $\zeta(\boldsymbol{\theta})$ characterizes the impact of power allocation on outage probability. Specifically, $\zeta(\boldsymbol{\theta})$ is an inverse power function of the the product of power allocation factors ${\prod\nolimits_{k=1}^{K}{{\theta_{k}}}}$. It would decrease as the product of power allocation factors $\boldsymbol{\theta}$ increases, which eventually results in the improvement of outage performance. With the definition of $\theta_{k}=\frac{{{P_{k}}}}{{{P_{total}}}}$, we also can conclude that the asymptotic outage probability is an inverse power function of the product of the transmission powers in all HARQ rounds. Notice that this clear quantitative relationship between the power allocation factors and the outage probability hasn’t been discovered even under independent fading channels. With this quantitative relationship, optimal power allocation to achieve various objectives is enabled. Taking the energy-limited outage minimization as an example, the power allocation problem can be formulated as

[TABLE]

where ${\sum\nolimits_{k=1}^{K}{{P_{k}}{p_{out,k-1}}}}$ refers to the average transmission energy normalized to the codeword length [4], $p_{out,0}=1$ and $P_{\rm given}$ is the average energy constraint. By substituting $P_{k}={\theta_{k}}P_{total}$ and using the asymptotic outage probability ${p_{out\_{asy},k}}$ as the approximation of $p_{out,k}$, the optimization problem is reduced to the maximization of the product of power allocation factors ${\prod\nolimits_{k=1}^{K}{{\theta_{k}}}}$ subject to certain constraints. Accordingly, (40) can be easily converted into the generalized power optimization problem in [38, Eq.6]. Karush-Khun-Tucker conditions can then be applied to derive the optimal solution in closed-form as shown in [38, Eq.15].

IV-B3 Time Correlation Impact Factor $\varrho(\boldsymbol{\lambda})$

The time correlation impact factor $\varrho(\boldsymbol{\lambda})$ quantifies the impact of channel time correlation on the outage probability. Notice that when only one transmission is allowed, i.e., $K=1$, no time correlation is involved and we always have $\ell\left({{\boldsymbol{\lambda}},K}\right)=1$ and $\varrho(\boldsymbol{\lambda})=1$. When $K>1$, the time correlation impact factor $\varrho(\boldsymbol{\lambda})$ has one property as shown in the following lemma. For notational convenience, we define a partial ordering for two vectors ${\bf x}=(x_{1},x_{2},\cdots,x_{K}),{\bf y}=(y_{1},y_{2},\cdots,y_{K})\in\mathbb{R}^{K}$ as ${\bf x}\preceq\bf y$ if $x_{i}\leq y_{i},i=1,2,\cdots,K$.

Lemma 3**.**

Given $K>1$ and two time correlation vectors $\boldsymbol{\lambda}_{1}\preceq\boldsymbol{\lambda}_{2}$, we have

[TABLE]

[TABLE]

The left equalities in (41)-(42) hold if and only if $\boldsymbol{\lambda}_{1}=\bf 0$, while the right equalities in (41)-(42) hold if and only if $\boldsymbol{\lambda}_{1}=\boldsymbol{\lambda}_{2}$ .

Proof.

Please see Appendix E. ∎

It can be concluded from (37) and Lemma 3 that although the time correlation does not affect the diversity order, the increase of time correlation coefficients would cause the increases of $\varrho(\boldsymbol{\lambda})$ and the outage probability, thus resulting in the degradation of outage performance.

IV-B4 Coding and Modulation Gain $C(\mathcal{R})$

As defined in (37), $C(\mathcal{R})=\left({g_{{\bf{0}}}}\left(2^{\cal R}\right)\right)^{-\frac{1}{d}}$. For a given outage probability, the increase of $C(\mathcal{R})$ can result in the reduction of the SNR $\gamma$. In other words, $C(\mathcal{R})$ can quantify the amount of SNR reduction for a given outage probability under certain coding and modulation scheme. It is thus termed as coding and modulation gain [19]. After analyzing the function ${g_{{\bf{0}}}}\left(2^{\mathcal{R}}\right)$, we have the following property of the coding and modulation gain.

Lemma 4**.**

The function ${g_{{\bf{0}}}}\left(2^{\mathcal{R}}\right)$ is a monotonically increasing function of the information transmission rate $\mathcal{R}$, and is convex with respect to $\mathcal{R}$ for any fading order $m\geq 1$, and thus the coding and modulation gain $C(\mathcal{R})$ is a monotonically decreasing function of the information transmission rate $\mathcal{R}$.

Proof.

Please see Appendix F. ∎

It follows from (37) and Lemma 4 that the asymptotic outage probability is a monotonically increasing and convex function of the information transmission rate $\mathcal{R}$ when $m\geq 1$. In order to achieve a desired performance, the information transmission rate should be properly chosen. Owing to the simple analytical expression in (37), the optimal rate can be easily found. Taking the maximization of the long term average throughput (LTAT) given an allowable outage constraint $\epsilon$ as an example, the rate selection problem can be formulated as

[TABLE]

where $\bar{\mathcal{T}}$ denotes the LTAT. Clearly, the numerator ${\mathcal{R}\left({1-{p_{out,K}}}\right)}$ and the denominator ${\sum\nolimits_{k=0}^{K-1}{{p_{out,k}}}}$ in (43) are concave and convex with respect to $\mathcal{R}$, respectively, while the feasible region $\left\{\mathcal{R}\in\mathbb{R}_{+}:{p_{out,K}}\leq\epsilon\right\}$ is a convex set, when $m\geq 1$. As shown in [39], the optimization problem in (43) is a concave fractional programming problem and the globally optimal solution can be easily found using the techniques proposed in [39].

IV-B5 Quasi-Static Fading channels

As aforementioned, the probability distribution in Section III and the asymptotic outage probability in (37) are not applicable to the case of quasi-static fading channels with ${\boldsymbol{\lambda}}={\bf 1}$. Particularly, under quasi-static fading channels, the channel coefficients are constant among multiple transmissions, i.e., $|h_{1}|=\cdots=|h_{K}|\triangleq|h|\sim{\rm Nakagami}(m,{\sigma}^{2})$. Assuming constant transmission powers $P_{1}=\cdots=P_{K}\triangleq P$, that is, $\theta_{1}=\cdots=\theta_{K}\triangleq\theta$, the outage probability ${p_{out,K}}$ under quasi-static fading channels is readily obtained as

[TABLE]

By applying [31, Eq.8.354.1] into (44), the outage probability can be rewritten as

[TABLE]

Then the asymptotic outage probability under this quasi-static fading channels can be written as ${p_{out\_{asy},K}}=\zeta(\boldsymbol{\theta})\left(C(\mathcal{R}\right)\gamma)^{-m}$, where the power allocation impact factor is given as $\zeta(\boldsymbol{\theta})={\left({\frac{m}{{{\theta{\sigma}^{2}}}}}\right)^{m}}$, while the coding and modulation gain becomes $C(\mathcal{R})=\frac{{{{\left({\Gamma\left({m+1}\right)}\right)}^{1/m}}}}{{{2^{{\mathcal{R}}/K}}-1}}$ and the diversity order reduces to $d=m$ due to the full correlation of fading channels. It means that no time diversity can be achieved from multiple transmissions under quasi-static fading channels.

V Numerical Results And Optimal Design

The analytical results are now verified and optimal system design is discussed in this section. For illustration, we take systems with equicorrelated channels (i.e., $\boldsymbol{\lambda}_{\rm eq}=(\rho,\cdots,\rho)$) [40, 41, 42] and unit Nakagami spread ${\sigma_{1}}^{2}=\cdots={\sigma_{K}}^{2}=1$ as examples, unless otherwise indicated. In the following numerical analysis, the exact outage probability is approximated by (23) with the truncation order set as $N=3$ and the involved generalized Fox’s H function is efficiently calculated by the MATHEMATICA® program in [11] with a properly chosen Mellin-Barnes contour.

V-A Performance Evaluation

In Fig. 2, the outage probability $p_{out,K}$ is plotted against the transmission SNR $\gamma$ by setting $K=4$ and $\mathcal{R}=4$bps/Hz. The approximated results based on polynomial fitting technique [16] and inverse moment matching method [17] are also presented for comparison. Notice that polynomial fitting technique in [16] is proposed for HARQ-IR under Rayleigh fading channels and thus its result only for Rayleigh fading (i.e., $m=1$) is shown in Fig. 2. It is clear that the approximated results (23) coincide with the simulation results, while the asymptotic results (37) approach to the approximated/simulation results under high SNR regime. However, under the considered scenarios, neither [16] nor [17] can provide an accurate approximation under high SNR regime due to the fact that those proposed methods can only guarantee their convergence in MSE. In addition, the outage probabilities decrease with the increase of the transmission SNR $\gamma$ and the decreasing rate becomes larger when the fading order $m$ increases, because the diversity order is $mK$ and $p_{out\_asy,K}$ is directly proportional to $\gamma^{-mK}$. As expected, channel time correlation has a detrimental impact on outage probability. For a given fading order $m$, the outage probability curves corresponding to different correlations become parallel as $\gamma$ increases due to the same diversity order which is irrelevant to the correlation coefficient. These numerical results thus validate the results in Section IV.

To further investigate the effect of time correlation, Fig. 3 depicts time correlation impact factors $\varrho(\boldsymbol{\lambda})$ under two different channel correlation models, i.e., equal correlation and exponential correlation. For exponential correlation model, the correlation coefficients follow as $\boldsymbol{\lambda}_{\rm exp}=(\rho,\rho^{2},\cdots,\rho^{K})$ [41, 20]. It can be seen that the time correlation impact factor $\varrho(\boldsymbol{\lambda})$ increases as $\rho$ increases, thus causing the degradation of the outage performance. The curves become steeper with the increase of $\rho$, which means the outage performance degradation would become more significant when time correlation is high. By comparing the time correlation impact factors $\varrho(\boldsymbol{\lambda}_{\rm eq})$ and $\varrho(\boldsymbol{\lambda}_{\rm exp})$, we can find that $\boldsymbol{\lambda}_{\rm exp}\preceq\boldsymbol{\lambda}_{\rm eq}$ and $\varrho(\boldsymbol{\lambda}_{\rm exp})\leq\varrho(\boldsymbol{\lambda}_{\rm eq})$, which is consistent with Lemma 3. In addition, it is found that the increase of the maximum number of transmissions $K$ will lead to the increase of the time correlation impact factor $\varrho(\boldsymbol{\lambda})$ no matter under equal correlation model or under exponential correlation model. This is because that $\varrho(\boldsymbol{\lambda}_{K})=\varrho((\boldsymbol{\lambda}_{K},0))$ and $(\boldsymbol{\lambda}_{K},0)\preceq(\boldsymbol{\lambda}_{K},\lambda_{K+1})$. It follows from Lemma 3 that $\varrho(\boldsymbol{\lambda}_{K})\leq\varrho((\boldsymbol{\lambda}_{K},\lambda_{K+1}))$ and the equality holds if and only if $\lambda_{K+1}=0$. However, when $K$ increases, the increase of the time correlation impact factor will be offset by the decrease of the term $\left(C(\mathcal{R}\right)\gamma)^{-mK}$ and thus would not cause the degradation of the outage performance.

Finally, Fig. 4 illustrates the impacts of the information transmission rate $\cal R$ and the fading order $m$ on the coding and modulation gain $C(\cal R)$ under both time-correlated fading channels and quasi-static fading channels. Clearly, the coding and modulation gain $C(\cal R)$ decreases with the increase of the rate $\mathcal{R}$, which has already been proved in Lemma 4. Additionally, the coding and modulation gain over quasi-static fading channels is superior to the gain over time-correlated fading channels. The reason behind this is that multiple HARQ transmissions under quasi-static fading channels can be viewed as one single transmission with lower coding rate due to the same channel realization experienced, thus leading to the improvement of coding and modulation gain. However, such coding and modulation gain improvement is achieved in sacrifice of time diversity, that is, the diversity order of quasi-static fading channels reduces to $m$. Moreover, Fig. 4 shows that the coding and modulation gain $C(\mathcal{R})$ could benefit from the increase of fading order $m$.

V-B Optimal System Design

With the simple expression, the asymptotic result in (37) would facilitate optimal system design for HARQ-IR with low complexity333Specifically, the asymptotic outage probability (37) can be easily computed by a few arithmetical operations. Its computational complexity is much lower than the simulation results which require the calculations on a large number of channel realizations. For example, in order to obtain the outage probability of $10^{-5}$, more than $10^{6}$ channel realizations are required to be simulated, which definitely causes very high computational complexity.. Here optimal power allocation in (40) and optimal rate selection in (43) are particularly investigated as examples.

The optimal power allocation (OPA) schemes in (40) designed based on approximated outage probability (23) and asymptotic outage probability (37) are first compared with an optimal equal power allocation (OEPA) scheme where the outage probability is minimized as (40) with additional equal power allocation constraint (i.e., $P_{1}=\cdots=P_{K}$). The achieved optimal outage probabilities $p_{out,K}^{*}$ for a system with $K=2$, $\mathcal{R}=2$bps/Hz and $\rho=0.5$ are shown in Fig. 5. It is clear that the OPA solutions found based on approximated outage probability and asymptotic outage probability agree well and can lead to similar outage probability, which justifies the adoption of the asymptotic outage probability for system design. Notice that the OPA solution based on asymptotic outage probability can be easily found in closed-form as shown in [38] and its computational complexity is significantly reduced compared to the optimization based on approximated outage probability which requires exhaustive search with a large number of numerical computations. Moreover, OPA scheme performs better than OEPA since the transmission powers in all HARQ rounds are optimized.

Fig. 6 then illustrates the maximal LTAT $\bar{\mathcal{T}}$ achieved through optimal rate selection in (43) by setting $m=1$ and $K=4$. Similarly, the results based on the approximated outage probability and the asymptotic outage probability match well, further validating our asymptotic results. It is also observed that the maximal LTAT increases when the outage constraint $\epsilon$ increases and/or the channel time correlation reduces.

VI Conclusions

In this paper, asymptotic outage analysis has been conducted to thoroughly investigate the impacts of channel time correlation, transmit powers and information transmission rate on the performance of HARQ-IR over time correlated Nakagami fading channels. Clear insights have been discovered. Particularly, it has been revealed that the power allocation impact factor is an inverse power function of the product of power allocation factors, the time correlation impact factor is an increasing function of the channel time correlation coefficients, while the coding and modulation gain is a monotonically decreasing function of the information transmission rate. Therefore, high product of power allocation factors, low channel time correlation and low information transmission rate are favorable for improving outage performance. The simple form and special properties of asymptotic outage probability would effectively facilitate the optimal system design to achieve various objectives, e.g., optimal power allocation to minimize outage probability and optimal rate selection to maximize the LTAT.

Appendix A The derivation of ${f_{\left.G_{K}\right|T}}(x|t)$

The Mellin transform with respect to the conditional PDF of $G_{K}$ given $T=t$, ${f_{\left.G_{K}\right|T}}(x|t)$, can be written as

[TABLE]

Due to the independence of SNRs $\boldsymbol{\gamma}$ given $T=t$ and with the definition of $G_{K}\triangleq\prod\nolimits_{k=1}^{K}{\left({1+{\gamma_{k}}}\right)}$, $\phi\left({\left.s\right|t}\right)$ can be rewritten as

[TABLE]

Herein, with the conditional PDF ${{f_{\left.{{\gamma_{k}}}\right|T}}\left({\left.{{x_{k}}}\right|t}\right)}$ in (8), ${\rm E}\left\{{\left.{{{\left({1+{\gamma_{k}}}\right)}^{s-1}}}\right|t}\right\}$ can be derived as

[TABLE]

By using the series expansion of the confluent hypergeometric limit function [29, Eq.16.2.1], (48) is further derived as

[TABLE]

where $\Psi\left(\cdot\right)$ denotes Tricomi’s confluent hypergeometric function [31, Eq.9.211.4]. Plugging (A) into (47) together with some algebraic manipulations yields

[TABLE]

Meanwhile, by using inverse Mellin transform, the conditional PDF $f_{\left.G_{K}\right|T}(x|t)$ can be written as [30]

[TABLE]

where ${\rm i}=\sqrt{-1}$. Putting (50) into (51) finally leads to (9).

Appendix B Proof of Theorem 1

Based on (III-A), the CDF of the product of multiple correlated shifted-Gamma RVs ${F_{G_{K}}}\left(x\right)$ directly follows as (19) where the CDF of $\mathcal{A}_{\boldsymbol{\ell}}$ can be written based on (III-A) as

[TABLE]

Regarding to the generalized Fox’s H function in (52), we have the following properties.

Property 1**.**

[TABLE]

Property 2**.**

[TABLE]

Notice that Property 1 is similar to the property of Meijer-G function in [31, Eq.9.31.5] and the property of generalized upper incomplete Fox’s H functions in [24, Eq.A.10], while Property 2 is similar to the property of Meijer-G function in [31, Eq.9.31.2] and the property of generalized upper incomplete Fox’s H functions in [24, Eq.A.9]. They can thus be proved using similar approaches in [24] and [31].

By using Property 1, (52) can be derived as

[TABLE]

Then applying Property 2 into (55) yields

[TABLE]

With the definition of generalized Fox’s H function [11, 12], (56) can be further derived as

[TABLE]

where $M_{K,0}^{0,K}\left(s\right)=\frac{1}{{\prod\nolimits_{j=1}^{K}{\Xi\left({0,1,\frac{1}{{{\Omega_{j}}}},m+{\ell_{j}}}\right)}}}$ and $\Xi\left({a,\alpha,A,\varphi}\right)={A^{\varphi+a+\alpha s-1}}\Psi\left({\varphi,\varphi+a+\alpha s;A}\right)$. The function $\Xi\left({a,\alpha,A,\varphi}\right)$ can be rewritten by using one important property of Tricomi’s confluent hypergeometric function $\Psi\left({\alpha,\gamma;z}\right)$ in [31, Eq.9.210.2] as

[TABLE]

where ${}_{1}F_{1}(\cdot)$ represents the confluent hypergeometric function [31, Eq.9.210.1]. Thus it follows that

[TABLE]

By using (59) together with the definition of generalized Fox’s H function, (B) can be further rewritten as

[TABLE]

Finally, applying Property 2 to (60) yields (20).

Appendix C Proof of (25)

Following from (24), the truncation error $\nabla(N)$ is upper bounded as

[TABLE]

where ${{F_{{{\cal A}_{\boldsymbol{\ell}}},N}^{\rm max}}\left(2^{\cal R}\right)}=\mathop{\max}\nolimits_{\sum\nolimits_{k=1}^{K}\ell_{k}\geq N+1}\left({{F_{{{\cal A}_{\boldsymbol{\ell}}}}}\left(2^{\cal R}\right)}\right)$ and $(\cdot)_{n}$ denotes Pochhammer symbol. Herein, it is readily found that $\mathop{\max}\nolimits_{\sum\nolimits_{k=1}^{K}\ell_{k}=N_{1}}\left({{F_{{{\cal A}_{\boldsymbol{\ell}}}}}\left(2^{\cal R}\right)}\right)\geq\mathop{\max}\nolimits_{\sum\nolimits_{k=1}^{K}\ell_{k}=N_{2}}\left({{F_{{{\cal A}_{\boldsymbol{\ell}}}}}\left(2^{\cal R}\right)}\right)$ if $N_{1}<N_{2}$ by the physical interpretation of ${{F_{{{\cal A}_{\boldsymbol{\ell}}}}\left(2^{\cal R}\right)}}$ in Section III. Thus we have ${{F_{{{\cal A}_{\boldsymbol{\ell}}},N}^{\rm max}}\left(2^{\cal R}\right)}=\mathop{\max}\nolimits_{\sum\nolimits_{k=1}^{K}\ell_{k}=N+1}\left({{F_{{{\cal A}_{\boldsymbol{\ell}}}}}\left(2^{\cal R}\right)}\right)$. By means of multinomial theorem [43, Eq.26.4.9], it yields

[TABLE]

Reformulating the summation term in the right hand side of (62) in terms of hypergeometric function ${}_{2}F_{1}(\cdot)$ [29,Eq.15.1.1], the truncation error $\nabla(N)$ is consequently upper bounded as (25).

Appendix D Proof of Theorem 2

Recalling from Theorem 1 that ${F_{{{{\mathcal{A}}}_{\boldsymbol{\ell}}}}}(x)$ is the CDF of the product of multiple independent shifted-Gamma RVs, i.e., $\mathcal{A}_{\boldsymbol{\ell}}=\prod\nolimits_{k=1}^{K}(1+R_{{\boldsymbol{\ell}},k})$ and $R_{{\boldsymbol{\ell}},k}\sim\mathcal{G}(m+\ell_{k},{\Omega_{k}})$, it can be expressed as

[TABLE]

By using Maclaurin series of exponential function, ${F_{{{{\mathcal{A}}}_{\boldsymbol{\ell}}}}}(x)$ in (63) can then be derived as

[TABLE]

where ${\bf n}=(n_{1},\cdots,n_{K})$ and ${g_{{\boldsymbol{\ell}}}}\left(x\right)$ is given by (28). After putting $\Omega_{k}=\frac{{P_{k}}{\sigma_{k}}^{2}\left({1-{\lambda_{k}}^{2}}\right)}{m\mathcal{N}_{0}}$ and (26) into (D), the CDF expression in (27) directly follows.

On the other hand, ${F_{{{{\cal A}}_{\boldsymbol{\ell}}}}}(x)$ can also be expressed in the form of Mellin-Barnes integral from (B) as

[TABLE]

where $c_{2}=-c_{1}$ and the fundamental strip444$\left\{{\mathcal{M}{{F_{{{{\cal A}}_{\boldsymbol{\ell}}}}}}}\right\}\left(s\right)$ exists for any complex number $s$ in the fundamental strip. of $\left\{{\mathcal{M}{{F_{{{{\cal A}}_{\boldsymbol{\ell}}}}}}}\right\}\left(s\right)$ implies $c_{2}\in(-\infty,0)$ because ${F_{{\mathcal{A}_{\boldsymbol{\ell}}}}}\left(x\right)$ admits[44, p400]

[TABLE]

where $\mathcal{O}(\cdot)$ denotes the big O notation and the first equation of (66) holds because of ${F_{{\mathcal{A}_{\boldsymbol{\ell}}}}}\left(x\right)=0$ for $x\leq 1$. Thus $c_{1}\in(0,\infty)$.

By adopting [31, Eq.9.210.2] into (D), it yields

[TABLE]

where ${}_{1}{F_{1}}\left({\alpha,\beta;\frac{1}{{{\Omega_{k}}}}}\right)$ can be expanded as

[TABLE]

Clearly, as $\gamma$ approaches to infinity, the dominant term in (68) is $1$.

Herein, we assume $c_{1}>m+{\rm max}\{\boldsymbol{\ell}\}$ because $c_{1}$ could be any point in $(0,\infty)$. Then by substituting (68) into (67), we have

[TABLE]

Putting (26) into (D), ${{F_{{{{\cal A}}_{\boldsymbol{\ell}}}}}(x)}$ can be finally written as

[TABLE]

Comparing (74) with (27) and setting the coefficients corresponding to ${\gamma^{-d_{{{\cal A}_{\boldsymbol{\ell}}}}}}$ equal, we finally have (31).

Appendix E Proof of Lemma 3

Defining ${\boldsymbol{\lambda}_{1}}=\left({{\lambda_{1}},\cdots,{\lambda_{K}}}\right)$ and given $K>1$, we first consider the simplest case where $\boldsymbol{\lambda}_{1}\preceq\boldsymbol{\lambda}_{2}$ and ${\boldsymbol{\lambda}_{2}}=\left({{\lambda_{1}}+{\Delta_{1}},\cdots,{\lambda_{K}}}\right)$ with ${\Delta_{1}}\geq 0$. With (38), $\ell\left({{\boldsymbol{\lambda}_{2}},K}\right)$ can be written as

[TABLE]

Since ${\Delta_{1}}\geq 0$, we have

[TABLE]

where the equality holds if and only if ${\Delta_{1}}=0$. Now gradually adding non-negative increments $\Delta_{i}$ to the other elements $\lambda_{i}$, $i=2,\cdots,K$, in ${\boldsymbol{\lambda}_{2}}$ and using the similar approach as (E), we can finally have $\ell\left({{\boldsymbol{\lambda}_{2}},K}\right)\leq\ell\left({{\boldsymbol{\lambda}_{1}},K}\right)$ for general ${\boldsymbol{\lambda}_{2}}=\left({{\lambda_{1}}+{\Delta_{1}},\cdots,{\lambda_{i}}+{\Delta_{i}},\cdots,{\lambda_{K}}+{\Delta_{K}}}\right)$ with ${\Delta_{i}}\geq 0$ where the equality holds if and only if ${\Delta_{i}}=0$. Clearly, when $K>1$, we have $\ell\left({{\boldsymbol{\lambda}_{1}},K}\right)\leq\ell\left({{\bf 0},K}\right)=1$, where the equality holds if and only if ${\boldsymbol{\lambda}_{1}}={\bf 0}$. Then the proof completes.

Appendix F Proof of Lemma 4

It is readily found from (28) that ${g_{\boldsymbol{\ell}}}\left(2^{\mathcal{R}}\right)$ is a monotonically increasing function of $\mathcal{R}$. Therefore the first derivative of ${g_{\boldsymbol{\ell}}}\left(2^{\mathcal{R}}\right)$ with respect to $\mathcal{R}$ is greater than [math], i.e., ${g_{\boldsymbol{\ell}}}^{\prime}\left(2^{\mathcal{R}}\right)>0$. By applying the property of derivatives of Laplace transform into (31) [30], it follows that

[TABLE]

where $\Theta=\ln 2$. Similarly, the second derivative of ${g_{\boldsymbol{\ell}}}\left(2^{\mathcal{R}}\right)$ with respect to $\mathcal{R}$ is given by

[TABLE]

Without loss of generality, we assume that $\ell_{1}={\rm max}\{\boldsymbol{\ell}\}$. By rewriting $s=(s-m-\ell_{1})+(m+\ell_{1})$, the integral in (79) can be derived as

[TABLE]

Putting (80) into (79) along with (78), it follows that

[TABLE]

where ${{\boldsymbol{\ell}}_{/1}}=\left({{\ell_{2}},{\ell_{3}},\cdots,{\ell_{K}}}\right)$ and $\bar{{\boldsymbol{\ell}}}=\left({{\ell_{1}}-1,{\ell_{2}},{\ell_{3}},\cdots,{\ell_{K}}}\right)$. Thus it proves that ${{g_{\boldsymbol{\ell}}}\left(2^{\cal R}\right)}$ is a convex function of $\mathcal{R}$ if ${\rm max}{\{\boldsymbol{\ell}\}}+m\geq 1$. Then the proof directly follows.