#  Hilbert space fragmentation imposed real spectrum of a non-Hermitian system

Somsubhra Ghosh1, K. Sengupta1, and Indranil Paul2 1School of Physical
Sciences, Indian Association for the Cultivation of Science, Kolkata 700032,
India.
2Université Paris Cité, CNRS, Laboratoire Matériaux et Phénomènes Quantiques,
75205 Paris, France.

###### Abstract

We show that constraints imposed by strong Hilbert space fragmentation (HSF)
along with the presence of certain global symmetries can provide a sufficient
condition for the reality of eigenspectra of non-Hermitian quantum systems;
such a reality cannot be guaranteed by global symmetries alone. We demonstrate
this insight for an interacting finite fermionic Nelson-Hatano chain. We show
analytically that strong HSF and real spectrum are both consequences of the
same dynamical constraints in the limit of large interaction, provided the
system has sufficient global symmetries. The spectrum stays real for
interactions above a finite critical value, where the system encounters a
many-body exceptional point. We provide a method to detect this exceptional
point using a local equal-time correlation function.

_Introduction._ — Non-Hermitian many-body Hamiltonians are of great current
interest for their relevance to open quantum systems, and also for their novel
properties without Hermitian analog reviews_nonH ; nonhlit1 ; nonhlit2 ;
nonhlit3 ; nonhlit4 ; nonhlit5 ; nonhlit6 ; nonhlit7 ; nonhlit8 ; nonhlit9 ;
nonhlit10 ; nonhlit11 ; nonhlit12 ; nonhlit13 ; nhdyn1 ; nhdyn2 ; nhdyn3 ;
nhdyn4 ; nhdyn5 ; nhdyn6 . One such feature in the spectrum is an exceptional
point, where certain energy eigenvalues and eigenfunctions coalesce, while
across it eigenvalues can transform from being real to complex reviews_EP ;
eptop1 ; eptop2 ; eptop3 ; eptop4 ; eptop5 ; eptop6 .

The purpose of the current work is to investigate an important related
feature. Namely, why the spectra of certain non-Hermitian Hamiltonians are
entirely real in some parameter regimes. Note, this question cannot be
addressed completely by invoking global symmetry properties. For example, a
property such as pseudo-Hermitcity only guarantees that complex eigenvalues,
if they appear, come in complex conjugate pairs Mostafazadeh2002a ;
Mostafazadeh2002b . Likewise, a so-called $\mathcal{P}\mathcal{T}$-symmetric
system has completely real eigenvalues only in the regime where all the energy
eigenfunctions are also simultaneously eigenfunctions of the
$\mathcal{P}\mathcal{T}$ operator bender2007 ; zyablovsky2014 ; ozdemir2019 ,
and the question remains as to what guarantees the latter.

In this work we study a canonical non-Hermitian system, namely the interacting
fermionic Hatano-Nelson model comprising of a finite chain of spinless
fermions with non-reciprocal hopping and with nearest-neighbor interaction
nonhlit11 , shown schematically in Fig. 1(a). It is well-known that this model
has a “phase” at large enough interaction in which all the many-body energy
eigenvalues are real zhang2022 . Our goal here is to examine microscopically
what gives rise to and protects this phase.

Our work is built upon earlier results which showed that, in the limit of
infinitely large interaction, the Hermitian version of the above model shows
strong Hilbert space fragmentation, where the Fock space breaks up into
dynamically disjoint pieces whose number scales exponentially with the system
size khemani2020 ; sala2020 . This phenomena is the focus of intense research
at present since it leads to non-ergodicity and the ability to generate exotic
non-equilibrium states which cannot be accessed in equilibrium rakovszky2020 ;
yang2020 ; tomasi2019 ; frey2022 ; hsf6 ; hsf7 ; hsf8 ; hsf9 ; hsf10 ;
ghosh2023 .

Figure 1: (Color Online) (a) Schematic of a spinless fermionic chain with non-
reciprocal hopping, and nearest neighbor interaction. Filled sites are marked
red. (b), (c) show the maximum imaginary part and the real parts of the
spectrum, respectively, as a function of the interaction $V_{1}$, for $L=10$
and $\gamma=0.2J$. $V_{1}^{c}$ is the critical interaction where the system
encounters a many-body exceptional point. Below $V_{1}^{c}$ complex conjugate
pairs of eigenvalues first appear. In (c) the two eigenvalues which coalesce
at the exceptional point are delineated in blue. (d) $V_{1}^{c}$ diverges as
$\gamma\rightarrow J$, and also with system size.

Our main results are the following. First, we show analytically that,
precisely in the fragmented limit there is a many-body similarity
transformation that maps the non-Hermitian system to a Hermitian one, thereby
guaranteeing the reality of the spectrum in this limit. In fact, fragmentation
and real spectrum are both consequences of the same dynamical constraints
which emerge in the limit of infinitely large interaction. Second, we show
that the spectrum persists to be real for arbitrarily large but finite
interaction, provided the system has sufficient global symmetry protection.
These symmetries impose a hidden Hermiticity in a subspace where the reality
of the spectrum is not guaranteed by the fragmentation limit. These
observations imply that the spectrum is real for interactions above a finite
critical value, where the system encounters a many-body exceptional point.
Finally we compute a local equal-time density-density fermionic correlation
function and show that it can be used to detect the position of this
exceptional point. Overall, our work provides the first analysis of the role
of dynamical constraints in determining the spectrum of a non-Hermitian
system.

_Model._ — In terms of fermionic creation and annihilation operators
$(c^{\dagger}_{i},c_{i})$ at site $i$ the Hamiltonian is

$\mathcal{H}=\sum_{i=1}^{L}\left[(J-\gamma)c^{\dagger}_{i}c_{i+1}+(J+\gamma)c^{\dagger}_{i+1}c_{i}+V_{1}\hat{n}_{i}\hat{n}_{i+1}\right],$
(1)

where $\hat{n}_{i}\equiv c^{\dagger}_{i}c_{i}$ is the number operator at site
$i$, and $L$ is the system size, see Fig. 1(a). The non-Hermiticity of
$\mathcal{H}$ is due to the non-reciprocal hopping parameter $\gamma>0$. We
study the system at half-filling with $\sum_{i}n_{i}=L/2$ where $n_{i}$ is the
fermion occupation at site $i$, and we impose periodic (anti) boundary
condition for total particle number odd (even). This choice ensures that the
system is traslationally invariant, as discussed in the Supplementary
information (SI) SI . For open boundary condition and $\gamma<J$ the problem
is trivial because there is a one-body similarity transformation which makes
$\mathcal{H}$ Hermitian eptop3 . As discussed in detail in the SI SI ,
$\mathcal{H}$ is pseudo-Hermitian, while its global symmetries are
$\mathcal{G}=(\mathcal{P}\mathcal{C},\mathcal{R})$ with
$[\mathcal{H},\mathcal{G}]=0$, where $(\mathcal{P},\mathcal{C},\mathcal{R})$
are parity, charge conjugation and translation by one site, respectively.
Furthermore, since integrability plays no role, our results are valid even in
the presence of next-nearest neighbor interaction.

The spectral properties of $\mathcal{H}$ are summarized in Fig. 1, (b) - (d).
Panels (b) and (c) show that, for $\gamma<J$, the spectrum is real for $V_{1}$
above a critical value $V_{1}^{c}$. At $V_{1}^{c}$ one pair (or two pairs) of
eigenvalues and eigenvectors coalesce at a many-body exceptional point, and
below $V_{1}^{c}$ they become complex conjugate pairs, as dictated by the
pseudo-Hermiticity of $\mathcal{H}$. This is shown in blue in (c). Panel (d)
shows that $V_{1}^{c}$ diverges as $\gamma\rightarrow J$, and also with system
size $L$. Thus, the regime with real spectrum is relevant only for finite size
systems. Note, while these features were reported in a recent work zhang2022 ,
the link between Hilbert space fragmentation and symmetry protection with the
reality of the spectrum, which is the focus of this work, has not been
explored earlier.

_Limit of fragmentation._ — We consider the limit of large interaction,
keeping terms to linear order in $(J,\gamma)$ and ignoring those of order
$(J,\gamma)^{2}/V_{1}$ and smaller. This gives dias2000 ,

$\displaystyle\mathcal{H}\approx\mathcal{H}_{f}$
$\displaystyle=\sum_{i=1}^{L}\left[\hat{P}_{i}\left((J-\gamma)c^{\dagger}_{i}c_{i+1}+(J+\gamma)c^{\dagger}_{i+1}c_{i}\right)\hat{P}_{i}\right.$
$\displaystyle\left.+V_{1}\hat{n}_{i}\hat{n}_{i+1}\right],$ (2)

where the projector $\hat{P}_{i}\equiv 1-(\hat{n}_{i-1}-\hat{n}_{i+2})^{2}$
ensures that the hopping is constrained, and is allowed only if the process
does not change the total number of nearest neighbor occupations
$\hat{N}_{d}\equiv\sum_{i}\hat{n}_{i}\hat{n}_{i+1}$ . Thus, by suitably
increasing $V_{1}$, the spectra of $\mathcal{H}_{f}$ and $\mathcal{H}$ can be
made to coincide with arbitrary accuracy, as shown in the SI.

The Hermitian version of $\mathcal{H}_{f}$ has been shown to display strong
Hilbert space fragmentation tomasi2019 ; frey2022 . Since fragmentation does
not depend on whether the hopping mediated connectivity between the many-body
Fock states is reciprocal or not, the non-Hermitian $\mathcal{H}_{f}$ will
show the same property. Below we prove that the dynamical constraints that
give rise to fragmentation also allows the existence of a many-body similarity
transformation that maps $\mathcal{H}_{f}$ into a Hermitian form for
$\gamma<J$.

_Similarity transformation._ — The first step of the proof is to label the
many-body states. Traditionally, this is done using “spins and movers”
dias2000 ; tomasi2019 ; frey2022 . Here we take a different strategy, and we
label them by “defects”. A “particle-defect” and a “hole-defect” are two
occupied or two unoccupied nearest-neighbor sites, respectively. Due to half-
filling particle- and hole-defects appear in pair, and their numbers are
conserved, since $\left[\mathcal{H}_{f},\hat{N}_{d}\right]=0$. Thus, the
Hilbert space factorizes into sectors with eigenvalue
$N_{d}=0,1,\ldots,L/2-1$.

All dynamically frozen (i.e., zero connectivity) states have real energies.
This is the case for $N_{d}=0$, where the two defect-free wavefunctions with
occupations $|0,1,0,1,\ldots\rangle$ and $|1,0,1,0,\ldots\rangle$ have zero
energy. For $N_{d}\neq 0$ we label a defect position by the location of the
first of the two nearest-neighbor sites. Thus, any state with $N_{d}=1$ has
label $|(i)(j)\rangle$, where $i$ and $j$ are locations of the particle- and
hole-defect, respectively. Likewise, a state with $N_{d}=2$ is labeled by
$|(i_{1},i_{2})(j_{1},j_{2})\rangle$, and $N_{d}=n$ by $|(i_{1},i_{2},\ldots
i_{n})(j_{1},j_{2},\ldots j_{n})\rangle$. Since the fermions are
indistinguishable, permutations of $i$ and of $j$ imply the same state. Thus,
the states shown in Fig. 2(a) have labels $|(5)(7)\rangle$ and
$|(8)(4)\rangle$ with $N_{d}=1$, and $|(2,6)(4,8)\rangle$ and
$|(2,3)(5,8)\rangle$ with $N_{d}=2$.

Due to half-filling the defect locations follow certain rules. (a) If two
particle-defects at $i_{1}$ and $i_{2}$ are “adjacent”, then $(i_{1},i_{2})$
can only be (odd, even) or (even, odd). The same applies for two adjacent
hole-defects. Here, “adjacent” does not imply defects located right next to
one another. Two defects are adjacent if there is no third defect in between
the two while traversing either clockwise or counterclockwise. (b) If a
particle-defect at $i_{1}$ is adjacent to a hole-defect at $j_{1}$, then
$(i_{1},j_{1})$ can only be (even, even) or (odd, odd). One can verify that
the wavefunctions in Fig. 2(a) follow these rules.

Figure 2: (Color Online) (a) Examples of labeling many-body wavefunctions by
the location of “defects”. Two nearest neighbor sites form a particle- or a
hole-defect when they are both occupied or both unoccupied, respectively. The
defect position is the location of the first of the two sites from the left.
(b) and (c) are examples of connectivities for $L=10$, $N_{d}=2,1$,
respectively, see text. Solid (green) and dashed (brown) arrows denote
fermions hopping to the right and left, respectively. Reversing an arrow
direction also implies exchanging solid $\leftrightarrow$ dashed. (d) Three
examples of non-reciprocal hopping over four sites. (i) has open boundary,
while (ii) and (iii) have periodic boundary conditions. $1,r,r^{2}$, _et
cetera_ (in blue) in (b)-(d) are the scaling factors $\lambda$, such as in Eq.
3, of the wavefunctions next to them, which define the similarity
transformation, wherever possible. In (d, ii) one of the sites, indicated by a
question symbol, cannot be scaled consistently. For a closed loop a similarity
transformation is only possible if there are equal number of solid and dashed
arrows while traversing the loop clockwise or anticlockwise, as in (d, iii).

The second step is to determine the defect dynamics which, due to the
constrained hopping, obey the following rules. (i) An allowed fermion hop
changes $i$ or $j$ by $\pm$2 modulo $L$. (ii) Since second nearest neighbor
hopping is absent, two defects cannot cross each other. It can be shown that
due to rules (i) and (ii) each sector of $N_{d}$ breaks into exponentially
large number of disjoint subsectors, i.e. fragments, that scale as $e^{N_{d}}$
sala2020 .

The third step is to establish the constrained hopping induced connectivity
between the many-body wavefunctions within each non-trivial subsector. There
is no general pattern for these connectivities, and they need to be worked out
case by case, even though the proof below holds for all the connectivities. To
show few examples, Fig. 2(b) is the connectivity for $L=10$, $N_{d}=2$ with
$(i_{1},i_{2})(j_{1},j_{2})$ (odd, odd)(odd, odd), while Fig. 2(c) is for
$L=10$, $N_{d}=1$ with $(i)(j)$ (odd)(odd). The dashed and solid arrows denote
fermions hopping to the left (with amplitude $J_{1}\equiv J-\gamma$) and right
(with amplitude $J_{2}\equiv J+\gamma$), respectively. Reversing an arrow
implies the exchange $J_{1}\leftrightarrow J_{2}$. A fermion hopping left can
result either a particle-defect to move left, i.e., $i\rightarrow(i-2)$ mod
$L$, or a hole-defect to move right, i.e., $j\rightarrow(j+2)$ mod $L$. Thus,
the connectivity diagram of each sub-sector can be viewed as a single non-
interacting “particle” hopping in the abstract space of many-body
wavefunctions in a non-reciprocal manner.

The fourth and final step of the proof is to establish the existence of the
similarity transformation in each sub-sector. For pedagogical reason we first
consider few examples of non-reciprocal hopping of a single particle in a
four-site system. Fig. 2(d, i) is a linear chain with open boundary condition.
This can be mapped to a Hermitian form for $\gamma<J$ by the scaling

$|i\rangle\rightarrow\lambda_{i}|i\rangle,\quad\langle
i|\rightarrow(1/\lambda_{i})\langle i|,$ (3)

where $\lambda_{i}=1,r,r^{2},r^{3}$, for $i=1,\ldots,4$, respectively, and
$r\equiv\sqrt{J_{2}/J_{1}}$ eptop3 . However, for periodic boundary condition,
as in Fig. 2(d, ii), the transformation will not work since the new link 4-1
cannot be made Hermitian. Thus, finding similarity transformations is non-
trivial where the connections form closed loops, which is precisely our case
as seen in Fig. 2(b, c). Now consider Fig. 2(d, iii) which is also a closed
loop, but where the hops are $J_{2},J_{2},J_{1},J_{1}$ moving clockwise. In
this case, once again, a similarity mapping exist, with
$\lambda_{i}=1,r,r^{2},r$, respectively. This example illustrates the
principle that a closed loop which has equal number of $J_{1}$ and $J_{2}$
hops while traversing it clockwise or anticlockwise can be mapped to a
Hermitian form. This is because every $J_{2}$ link requires an additional
scaling of $r$ for the second site compared to the first, while a $J_{1}$ link
requires a scaling of $1/r$. This is exactly the case for the connectivities
of Fig. 2 (b, c), where the scalings associated with the wavefunctions are
marked in blue. Additional examples of such scalings are discussed in the SI
SI . We prove below that _all_ the connections of $\mathcal{H}_{f}$ are such
that each and every possible loop has this property.

Starting from any $|(i_{1},i_{2},\ldots i_{n})(j_{1},j_{2},\ldots
j_{n})\rangle$, a closed loop is obtained in three basic ways.

(1) If one or more of the site indices change as, say, $i\rightarrow
i^{\prime}\rightarrow i^{\prime\prime}$ and so on, and then reverse the path
to go back to $i$, while obeying the rules (i) and (ii). Since the reverse of
a $J_{1}$ hop is a $J_{2}$ hop, and vice versa, traversing the loop along one
direction will necessarily have equal number of $J_{1}$ and $J_{2}$ hops. The
loop (1)(7) $\rightarrow$ (1)(5) $\rightarrow$ (3)(5) $\rightarrow$ (3)(7)
$\rightarrow$ (1)(7) in Fig. 2(c) is an example.

(2) If a defect does not retrace its path, but moves across the chain,
traversing $L$ sites, and gets back to its original position using the
periodic boundary condition. However, according to rule (ii) this can happen
only if all other defects perform the same circular motion in the same
direction and regain their original positions, each having traversed $L$
sites. Since a particle-defect moving to the right is a $J_{2}$ hop, while a
hole-defect moving to the right is a $J_{1}$ hop, and since there are equal
number of particle- and hole-defects, this loop, too, will have equal number
of $J_{1}$ and $J_{2}$ hops. Starting from state (1)(3) on the left side of
Fig. 2(c) and ending again at (1)(3) on the right side of the figure is an
example of such a loop.

(3) In some cases, such as in Fig. 2(b), it is possible for the defects to
exchange positions such that $i_{1}\rightarrow i_{2}\rightarrow
i_{3}\ldots\rightarrow i_{n}\rightarrow i_{1}$, and $j_{1}\rightarrow
j_{2}\rightarrow j_{3}\ldots\rightarrow j_{n}\rightarrow j_{1}$. In this case
a loop is completed by permuting the indices, while the defects neither
retrace their paths nor complete the circle. Here the sum of the sites
traversed by all the particle-defects is $L$ and the same is true for all the
hole-defects, and they are along the same direction. Thus, again here the loop
has equal number of $J_{1}$ and $J_{2}$ hops.

This completes the proof that $\mathcal{H}_{f}$ can be mapped into a Hermitian
form for $\gamma<J$; this feature guarantees reality of eigenspectrum of
${\mathcal{H}}_{f}$ in this limit.

_Finite $V_{1}^{c}$ and symmetry protection_.— The above conclusion, however,
is not sufficient to deduce that the spectrum remains real once $V_{1}$ is
finite. To understand why, consider two eigenstates of $\mathcal{H}_{f}$ from
the same sub-sector of $\hat{N}_{d}$. Measuring energies from the average
eigenvalue, the sub-system has the structure

$\mathcal{H}_{2}=\begin{pmatrix}l&m_{1}+m_{2}\\\ m_{1}-m_{2}&-l\end{pmatrix},$

with eigenvalues $\pm\sqrt{l^{2}+m_{1}^{2}-m_{2}^{2}}$. Since
$m_{1,2}\sim\mathcal{O}((J,\gamma)^{2}/V_{1})$ or smaller, for finite $l$ the
reality of the eigenvalues is guaranteed for $V_{1}$ sufficiently large. But,
this argument fails when the two states are degenerate and $l=0$. However, in
this subspace the reality of the spectrum can still be protected by the global
symmetries $\mathcal{G}$, provided any two degenerate states are related by
$\mathcal{G}|\psi\rangle=|\phi\rangle$, since it implies that the off-diagonal
matrix elements
$\langle\langle\psi|\mathcal{H}|\phi\rangle=\langle\langle\phi|\mathcal{H}|\psi\rangle$.
Here $(|\psi\rangle,|\phi\rangle)$ and
$(|\psi\rangle\rangle,|\phi\rangle\rangle)$ are the right and left
eigenvectors, respectively, of $\mathcal{H}_{f}$ in the degenerate subspace.
Since these off-diagonal terms are also real (because it is possible to choose
the eigenvectors of $\mathcal{H}_{f}$ to be real), we conclude that in this
subspace there can be a hidden Hermiticity of $\mathcal{H}$ which is symmetry
protected. Thus, the dynamical constraints and the global symmetries together
ensure that the spectrum stay real for $V_{1}$ greater than a finite value
$V_{1}^{c}$.

As discussed in the SI, the above symmetry protection can be destroyed by a
suitable choice of boundary condition. In this case one can have complex
eigenvalues for any finite value of $V_{1}$, even though the spectrum of the
corresponding $\mathcal{H}_{f}$ is real SI .

The above argument also implies that the subspace $\mathcal{H}_{2}$ that
eventually triggers the exceptional point at $V_{1}^{c}$ cannot be symmetry
related. In this case $l\sim\sqrt{J^{2}-\gamma^{2}}/e^{cL}$ is the average
level spacing of the sub-sector, and the constant $c>0$ depends on the sub-
sector size. Empirically, we find that $m_{1}$ is at least one order of
magnitude smaller than $m_{2}$, while $m_{2}\propto
V_{1}^{-\alpha}(J^{2}-\gamma^{2})^{-\beta/2}$ where the exponents
$(\alpha,\beta)$ are $L$-dependent. This implies that $V_{1}^{c}\sim
Je^{cL/\alpha}/(1-(\gamma/J)^{2})^{(\beta+1)/(2\alpha)}$. Thus, $V_{1}^{c}$
diverges exponentially with the system size, and it diverges as a power-law
with a $L$-dependent exponent for $\gamma\rightarrow J$. These features are
illustrated in Fig. 1(d).

Note, in passing, that for certain values of $L$ the two coalescing levels at
$V_{1}^{c}$ are each doubly degenerate, so that below $V_{1}^{c}$ there are
two pairs of complex conjugate eigenvalues. This degeneracy is related to
$\mathcal{P}\mathcal{C}$ invariance SI .

Figure 3: (Color Online) (a)-(c) Time evolution of the correlation function
$\chi(t)$ for $V_{1}$ = 3J, 4.285J, and 6J, respectively, for $L=14$ and
$\gamma=0.2J$. In (b) the system is very close to the exceptional point at
$V_{1}^{c}\approx 4.2863$. (d) Variation of the relaxation timescale $\tau$
with $V_{1}$ showing a one-sided divergence as $V_{1}\rightarrow V_{1}^{c}$
from below.

_Correlation function._ — Finally, we study the time evolution of a
correlation function which, in principle, can be measured to identify the
location of the many-body exceptional point. One such example is
$\chi(t)=\langle\psi(t)|\hat{N}_{d}|\psi(t)\rangle/L$ starting from a random
Fock state $|\psi(0)\rangle=\sum_{m}c_{m}|\phi_{m}\rangle$, expanded in the
basis of the right eigenvectors $|\phi_{m}\rangle$ of $\mathcal{H}$. The time-
evolved wavefunction, suitably normalized to account for the non-Hermiticity
of $\mathcal{H}$, is

$|\psi(t)\rangle=\frac{e^{-i\mathcal{H}t/\hbar}|\psi(0)\rangle}{||e^{-i\mathcal{H}t/\hbar}|\psi(0)\rangle||}=\frac{\sum_{m}c_{m}(t)|\phi_{m}\rangle}{\sqrt{\sum_{m,n}c_{m}^{*}(t)c_{n}(t)\langle\phi_{m}|\phi_{n}\rangle}},$

where $c_{m}(t)=c_{m}e^{-i\epsilon_{m}t/\hbar}$.

Fig. 3, panels (a)-(c) show the time evolution of $\chi(t)$ for $V_{1}$ less
than, nearly equal to, and greater than $V_{1}^{c}$, respectively. For
$V_{1}<V_{1}^{c}$ the time evolution is dominated by the eigenvalue with the
largest imaginary component. Consequently, after a timescale $\tau\sim 1/{\rm
Max[Im}\,\epsilon]$, the correlation function attains a steady state value
$\chi(t\gg\tau)\sim 1/L\langle\phi_{m}^{*}|\hat{N}_{d}|\phi_{m}^{*}\rangle$,
where $|\phi_{m}^{*}\rangle$ is the eigenvector with the largest ${\rm
Im}\,\epsilon$. This implies that $\tau$ diverges as $V_{1}\rightarrow
V_{1}^{c}$ from below, as seen in Fig. 3(d). For $V_{1}\geq V_{1}^{c}$ all the
eigenvalues are real and the system quickly attains a diagonal ensemble, and
$\chi(t)$ fluctuates about an average value $\chi(t\gg\tau)\sim
1/L\sum_{m}|c_{m}|^{2}\langle\phi_{m}|\hat{N}_{d}|\phi_{m}\rangle$ rigol1 ;
rigol2 , implying that the the peak of $\tau(V_{1})$ in Fig. 3(d) is one-
sided. This peak can be used to detect the exceptional point.

_Conclusion._ — To summarize, we explained microscopically why the interacting
fermionic Hatano-Nelson model has purely real many-body spectrum for nearest
neighbor interaction $V_{1}>V_{1}^{c}$. This is a consequence of two
ingredients. First, the dynamical constraints in the infinitely large
interaction limit which also fragments the Hilbert space of the model. Second,
the global symmetries of the Hamiltonian. While the role of the global
symmetries has been widely studied, that of the first ingredient has not been
explored earlier. Thus, we reveal a deep link between the physics of
fragmentation and the property of real spectrum of an interacting non-
Hermitian system. This link is worth investigating in the future for other
interacting non-Hermitian models.

_Acknowledgement._ – The authors thank Diptiman Sen for several comments. IP
is thankful to Masudul Haque for insightful discussions. SG acknowledges the
financial support provided by CSIR, India through file
09/080(1133)/2019-EMR-I. KS thanks DST, India for support through SERB project
JCB/2021/000030.

## References

  * (1) For reviews see, e.g., I. Rotter, J. Phys. A: Math. Theor. 42, 153001 (2009); Y. Ashida, Z. Gong, and M. Ueda, Adv. Phys. 69, 249 (2021); I. Rotter and J. P. Bird, Rep. Prog. Phys. 78, 114001 (2015).
  * (2) J. Gonzalez and R. A. Molina, Phys. Rev. B 96, 045437 (2017); V. Kozii and L. Fu, arXiv:1708.05841 (unpublished); A. A. Zyuzin and A. Y. Zyuzin, Phys. Rev. B 97, 041203(R) (2018); H. Shen and L. Fu, Phys. Rev. Lett. 121, 026403 (2018); R. A. Molina and J. Gonzalez, Phys. Rev. Lett. 120, 146601 (2018);T. Yoshida, R. Peters, and N. Kawakami, Phys. Rev. B 98, 035141 (2018); J. Carlstrom and E. J. Bergholtz, Phys. Rev. A 98, 042114 (2018).
  * (3) T. M. Philip, M. R. Hirsbrunner, and M. J. Gilbert, Phys. Rev. B 98, 155430 (2018); Y. Chen and H. Zhai, Phys. Rev. B 98, 245130 (2018); K. Moors, A. A. Zyuzin, A. Y. Zyuzin, R. P. Tiwari, and T. L. Schmidt, Phys. Rev. B 99, 041116(R) (2018); R. Okugawa and T. Yokoyama, Phys. Rev. B 99, 041202(R) (2019); J. C. Budich, J. Carlstrom, F. K. Kunst, and E. J. Bergholtz, Phys. Rev. B 99, 041406(R) (2019).
  * (4) Z. Yang and J. Hu, Phys. Rev. B 99, 081102(R) (2019); T. Yoshida, R. Peters, N. Kawakami, and Y. Hatsugai, Phys. Rev. B 99, 121101(R) (2019); Y. Wu, W. Liu, J. Geng, X. Song, X. Ye, C.-K. Duan, X. Rong, and J. Du, Science 364, 878 (2019); P. San-Jose, J. Cayao, E. Prada, and R. Aguado, Sci. Rep. 6, 21427 (2016); Q.-B. Zeng, B. Zhu, S. Chen, L. You, and R. Lu, Phys. Rev. A 94, 022119 (2016); C. Li, X. Z. Zhang, G. Zhang, and Z. Song, Phys. Rev. B 97, 115436 (2018); J. Cayao and A. M. Black-Schaffer Phys. Rev. B 105, 094502 (2022); R. Arouca, J. Cayao, A. M. Black-Schaffer, arXiv:2206.15324 (unpublished).
  * (5) K. Kawabata, Y. Ashida, H. Katsura, and M. Ueda, Phys. Rev. B 98, 085116 (2018); A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez, G. A. Siviloglou, and D. N. Christodoulides, Phys. Rev. Lett. 103, 093902 (2009); C. E. Ruter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, and D. Kip, Nat. Phys. 6, 192 (2010); L. Feng, M. Ayache, J. Huang, Y.-L. Xu, M.-H. Lu, Y.-F. Chen, Y. Fainman, and A. Scherer, Science 333, 729 (2011); A. Regensburger, C. Bersch, M.-A. Miri, G. Onishchukov, D. N. Christodoulides, and U. Peschel, Nature (London) 488, 167 (2012).
  * (6) L. Feng, Y.-L. Xu, W. S. Fegadolli, M.-H. Lu, J. E. Oliveira, V. R. Almeida, Y.-F. Chen, and A. Scherer, Nat. Mater. 12, 108 (2013); C. Poli, M. Bellec, U. Kuhl, F. Mortessagne, and H. Schomerus, Nat. Commun. 6, 6710 (2015); B. Zhen, C.W. Hsu, Y. Igarashi, L. Lu, I. Kaminer, A. Pick, S.-L. Chua, J. D. Joannopoulos, and M. Solja.i., Nature (London) 525, 354 (2015); H. Zhao, S. Longhi, and L. Feng, Sci. Rep. 5, 17022 (2015); K. Ding, Z. Q. Zhang, and C. T. Chan, Phys. Rev. B 92, 235310 (2015).
  * (7) S. Weimann, M. Kremer, Y. Plotnik, Y. Lumer, S. Nolte, K. Makris, M. Segev, M. Rechtsman, and A. Szameit, Nat. Mater. 16, 433 (2017); H. Hodaei, A. U. Hassan, S. Wittek, H. Garcia-Gracia, R. El-Ganainy, D. N. Christodoulides, and M. Khajavikhan, Nature (London) 548, 187 (2017); W. Chen, K. Ozdemir, G. Zhao, J. Wiersig, and L. Yang, Nature (London) 548, 192 (2017); P. St-Jean, V. Goblot, E. Galopin, A. Lemaitre, T. Ozawa, L. Le Gratiet, I. Sagnes, J. Bloch, and A. Amo, Nat. Photonics 11, 651 (2017).
  * (8) B. Bahari, A. Ndao, F. Vallini, A. E. Amili, Y. Fainman, and B. K. Le, Science 358, 636 (2017); J. Wang, H. Y. Dong, Q. Y. Shi, W. Wang, and K. H. Fung, Phys. Rev. B 97, 014428 (2018); H. Zhou, C. Peng, Y. Yoon, C.W. Hsu, K. A. Nelson, L. Fu, J. D. Joannopoulos, M. Solja.i., and B. Zhen, Science 359, 1009 (2018); M. Parto, S.Wittek, H. Hodaei, G. Harari, M. A. Bandres, J. Ren, M. C. Rechtsman, M. Segev, D. N. Christodoulides, and M. Khajavikhan, Phys. Rev. Lett. 120, 113901 (2018); H. Zhao, P. Miao, M. H. Teimourpour, S. Malzard, R. El-Ganainy, H. Schomerus, and L. Feng, Nat. Commun. 9, 981 (2018).
  * (9) G. Harari, M. A. Bandres, Y. Lumer, M. C. Rechtsman, Y. D. Chong, M. Khajavikhan, D. N. Christodoulides, and M. Segev, Science 359, 1230 (2018); M. A. Bandres, S. Wittek, G. Harari, M. Parto, J. Ren, M. Segev, D. N. Christodoulides, and M. Khajavikhan, Science 359, 1231 (2018); M. Pan, H. Zhao, P. Miao, S. Longhi, and L. Feng, Nat. Commun. 9, 1308 (2018); L. Jin and Z. Song, Phys. Rev. Lett. 121, 073901 (2018); S. Malzard and H. Schomerus, Phys. Rev. A 98, 033807 (2018); Z. Oztas and C. Yuce, Phys. Rev. A 98, 042104 (2018).
  * (10) M. Kremer, T. Biesenthal, L. J. Maczewsky, M. Heinrich, R. Thomale, and A. Szameit, Nat. Commun. 10, 435 (2019); K. Y. Bliokh, D. Leykam, M. Lein, and F. Nori, Nat. Commun. 10, 580 (2019); S. Wang, B. Hou, W. Lu, Y. Chen, Z. Zhang, and C. Chan, Nat. Commun. 10, 832 (2019); S. Chen,W. Zhang, B. Yang, T.Wu, and X. Zhang, Sci. Rep. 9, 5551 (2019); T. E. Lee and C.-K. Chan, Phys. Rev. X 4, 041001 (2014); Y. Xu, S.-T. Wang, and L.-M. Duan, Phys. Rev. Lett. 118, 045701 (2017); Y. Ashida, S. Furukawa, and M. Ueda, Nat. Commun. 8, 15791 (2017); Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Phys. Rev. X 8, 031079 (2018); M. Nakagawa, N. Kawakami, and M. Ueda, Phys. Rev. Lett. 121, 203001 (2018); K. Takata and M. Notomi, Phys. Rev. Lett. 121, 213902 (2018); L. Pan, S. Chen, and X. Cui, Phys. Rev. A 99, 011601(R) (2019).
  * (11) J. Li, A. K. Harter, J. Liu, L. de Melo, Y. N. Joglekar, and L. Luo, Nat. Commun. 10, 855 (2019); T. Liu, Y.-R. Zhang, Q. Ai, Z. Gong, K. Kawabata, M. Ueda, and F. Nori, Phys. Rev. Lett. 122, 076801 (2019); M. S. Rudner and L. S. Levitov, Phys. Rev. Lett. 102, 065703 (2009); J. M. Zeuner, M. C. Rechtsman, Y. Plotnik, Y. Lumer, S. Nolte, M. S. Rudner, M. Segev, and A. Szameit, Phys. Rev. Lett. 115, 040402 (2015); K. Mochizuki, D. Kim, and H. Obuse, Phys. Rev. A 93, 062116 (2016); L. Xiao, X. Zhan, Z. Bian, K. Wang, X. Zhang, X. Wang, J. Li, K. Mochizuki, D. Kim, N. Kawakami et al., Nat. Phys. 13, 1117 (2017).
  * (12) N. Hatano and D. R. Nelson, Phys. Rev. Lett. 77, 570 (1996); N. Hatano and D. R. Nelson, Phys. Rev. B 56, 8651 (1997); N. Hatano and D. R. Nelson, Phys. Rev. B 58, 8384 (1998).
  * (13) J. A. S. Lourenco, R. L. Eneias, and R. G. Pereira, Phys. Rev. B 98, 085126 (2018); E. I. Rosenthal, N. K. Ehrlich, M. S. Rudner, A. P. Higginbotham, and K.W. Lehnert, Phys. Rev. B 97, 220301(R) (2018); M. Wang, L. Ye, J. Christensen, and Z. Liu, Phys. Rev. Lett. 120, 246601 (2018).
  * (14) M. Ezawa, Phys. Rev. B 99, 121411(R) (2019); M. Ezawa, Phys. Rev. B 99, 201411(R) (2019); M. Ezawa, Phys. Rev. B 100, 045407 (2019).
  * (15) L. Zhou, Q.-h.Wang, H.Wang, and J. Gong, Phys. Rev. A 98, 022129 (2018); L. Zhou and Q. Du, New J. Phys. 23, 063041 (2021); B. Zhu, Y. Ke, H. Zhong, and C. Lee, Phys. Rev. Research 2, 023043 (2020); L. Zhou and J. Gong, Phys. Rev. B 98, 205417 (2018); L. Zhou, Phys. Rev. B 100, 184314 (2019); L. Zhou, Y. Gu, and J. Gong, Phys. Rev. B 103, L041404 (2021).
  * (16) L. Zhou and W. Han, Phys. Rev. B 106, 054307 (2022); C-H Liu, H. Hu, and S. Chen, Phys. Rev. B 105, 214305 (2022); L. Zhou, R. W. Bomantara, and S. Wu, SciPost Phys. 13, 015 (2022).
  * (17) S. Zamani, R. Jafari, and A. Langari, Phys. Rev. B 102, 144306 (2020); R. Jafari and A. Akbari, Phys. Rev. A 103, 012204 (2021); K. Yang, L. Zhou,W. Ma, X. Kong, P.Wang, X. Qin, X. Rong, Y.Wang, F. Shi, J. Gong, and J. Du, Phys. Rev. B 100, 085308 (2019); D. Chowdhury, A. Banerjee, and A. Narayan Phys. Rev. A 103, L051101 (2021).
  * (18) P. He and Z-H Huang, Phys. Rev. A 102, 062201 (2020); S. Longhi, J. Phys. A: Math. Theor. 50, 505201 (2017).
  * (19) X. Turkeshi and M. Schiro, arXiv:2201.09895 (unpublished); T. Banerjee and K. Sengupta, Phys. Rev. B 107, 155117 (2023).
  * (20) J. Ren, P. Hanggi, and B. Li, Phys. Rev. Lett. 104, 170601 (2010); J. Ren, S. Liu, and B. Li Phys. Rev. Lett. 108, 210603 (2012); H. Xu, D. Mason, L. Jiang and J. G. E. Harris, Nature 537, 80 (2016); Z. Wang, J. Chen, and J. Ren, Phys. Rev. E 106, L032102 (2020);L. J. Fernández-Alcazar, R. Kononchuk, H. Li, and T. Kottos, Phys. Rev. Lett. 126, 204101 (2021).
  * (21) for reviews see, e.g., E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Rev. Mod. Phys. 93, 015005 (2021); W. D. Heiss, J. Phys. A: Math. Theor. 45, 444016 (2012); M. Müller and I. Rotter, J. Phys. A: Math. Theor. 41, 244018 (2008).
  * (22) Y. C. Hu and T. L. Hughes, Phys. Rev. B 84, 153101 (2011); K. Esaki, M. Sato, K. Hasebe, and M. Kohmoto, Phys. Rev. B 84, 205128 (2011); T. E. Lee, Phys. Rev. Lett. 116, 133903 (2016); D. Leykam, K. Y. Bliokh, C. Huang, Y. D. Chong, and F. Nori, Phys. Rev. Lett. 118, 040401 (2017); V. M. Martinez Alvarez, J. E. Barrios Vargas, and L. E. F. Foa Torres, Phys. Rev. B 97, 121401(R) (2018); Y. Xiong, J. Phys. Commun. 2, 035043 (2018); H. Shen, B. Zhen, and L. Fu, Phys. Rev. Lett. 120, 146402 (2018).
  * (23) C. Yuce, Phys. Rev. A 97, 042118 (2018); C. Yin, H. Jiang, L. Li, R. Lu, and S. Chen, Phys. Rev. A 97, 052115 (2018); C. Yuce, Phys. Rev. A 98, 012111 (2018); F. K. Kunst, E. Edvardsson, J. C. Budich, and E. J. Bergholtz, Phys. Rev. Lett. 121, 026808 (2018); S. Yao, F. Song, and Z. Wang, Phys. Rev. Lett. 121, 136802 (2018); K. Kawabata, K. Shiozaki, and M. Ueda, Phys. Rev. B 98, 165148 (2018); C. Yuce and Z. Oztas, Sci. Rep. 8, 17416 (2018).
  * (24) S. Yao and Z. Wang, Phys. Rev. Lett. 121, 086803 (2018)
  * (25) K. Kawabata, S. Higashikawa, Z. Gong, Y. Ashida, and M. Ueda, Nat. Commun. 10, 297 (2019); L. Jin and Z. Song, Phys. Rev. B 99, 081103(R) (2019); H. Wang, J. Ruan, and H. Zhang, Phys. Rev. B 99, 075130 (2019); D. S. Borgnia, A. J. Kruchkov, and R.-J. Slager, Phys. Rev. Lett. 124, 056802 (2020); Z. Ozcakmakli Turker and C. Yuce, Phys. Rev. A 99, 022127 (2019); E. Edvardsson, F. K. Kunst, and E. J. Bergholtz, Phys. Rev. B 99, 081302(R) (2019).
  * (26) C.-H. Liu, H. Jiang, and S. Chen, Phys. Rev. B 99, 125103 (2019); C. H. Lee and R. Thomale, Phys. Rev. B 99, 201103(R) (2019); F. K. Kunst and V. Dwivedi, Phys. Rev. B 99, 245116 (2019); K. Yokomizo and S. Murakami, Phys. Rev. Lett. 123 066404 (2019).
  * (27) R. Nehra, and D. Roy, Phys. Rev. B 105, 195407 (2022); K. Kawabata, K. Shiozaki, and S. Ryu, Phys. Rev. B 105, 165137 (2022); K. Yang, D. Varjas, E. J. Bergholtz, S. Morampudi, and F. Wilczek, arXiv:2202.04435 (unpublished).
  * (28) A. Mostafazadeh, J. Math. Phys. 43, 205 (2002).
  * (29) A. Mostafazadeh, J. Math. Phys. 43, 2814 (2002).
  * (30) C. M Bender Rep. Prog. Phys. 70, 947 (2007).
  * (31) A. A. Zyablovsky, A. P. Vinogradov, A. A. Pukhov, A. V. Dorofeenko, A. A. Lisyansky, Phys.-Uspekhi 57, 1063 (2014).
  * (32) S. K. Özdemir, S. Rotter, F. Nori, and L. Yang, Nat. Mater. 18, 783 (2019).
  * (33) S.-B. Zhang, M. M. Denner, T. Bzdušek, M. A. Sentef, and T. Neupert, Phys. Rev. B 106, L121102 (2022).
  * (34) V. Khemani, M. Hermele and R. Nandkishore, Phys. Rev. B 101, 174204 (2020).
  * (35) P. Sala, T. Rakovszky, R. Verresen, M. Knap and F. Pollmann, Phys. Rev. X 10, 011047 (2020).
  * (36) T. Rakovszky, P. Sala, R. Verresen, M. Knap and F. Pollmann, Phys. Rev. B 101, 125126 (2020).
  * (37) Z.-C. Yang, F. Liu, A. V. Gorshkov and T. Iadecola, Phys. Rev. Lett. 124, 207602 (2020).
  * (38) G. De Tomasi, D. Hetterich, P. Sala, and F. Pollmann, Phys. Rev. B 100, 214313(2019).
  * (39) P. Frey, L. Hackl, and S. Rachel, Phys. Rev. B 106, L220301 (2022).
  * (40) S. Moudgalya and O. I. Motrunich, Phys. Rev. X 12, 011050 (2022); D. T. Stephen, O. Hart, and R. M. Nandkishore, arXiv:2209.03966 (unpublished); D. Hahn, P. A. McClarty, D. J. Luitz, SciPost Phys. 11, 074 (2021); N. Regnault and B. A. Bernevig, arXiv:2210.08019 (unpublished); D. Vu, K. Huang, X. Li, and S. Das Sarma, Phys. Rev. Lett. 128, 146601 (2022).
  * (41) T. Kohlert, S. Scherg, P. Sala, F. Pollmann, B. H. Madhusudhana, I. Bloch, and M. Aidelsburger, arXiv:2106.15586 (unpublished).
  * (42) B. Mukherjee, D. Banerjee, K. Sengupta, and A. Sen, Phys. Rev. B 104, 155117 (2021); P. Brighi, M. Ljubotina, and M. Serbyn, arXiv:2210.5607 (unpublished).
  * (43) J. Lehmann, P. Sala, F. Pollmann, and T. Rakovszky, arXiv:2208.12260 (unpublished).
  * (44) A. Chattopadhyay, B. Mukherjee, K. Sengupta, and A. Sen, arXiv:2208.13800 (unpublished).
  * (45) S. Ghosh, I. Paul, and K. Sengupta, Phys. Rev. Lett. 130, 120401 (2023).
  * (46) See supplementary information for more details.
  * (47) R. G. Dias, Phys. Rev. B 62, 7791 (2000).
  * (48) M. Rigol, V. Dunjko, M. Olshanii, Nature 452, 854 (2008).
  * (49) E. Khatami, G. Pupillo, M. Srednicki, M. Rigol, Phys. Rev. Lett. 111, 050403 (2013).


# Optimal Power Allocation for HARQ Schemes over Time-Correlated Nakagami-m
Fading Channels

Zheng Shi1 1Department of Electrical and Computer Engineering, University of
Macau, Macau
2Department of Electrical and Electronic Engineering, The University of Hong
Kong, Hong Kong  Shaodan Ma1 1Department of Electrical and Computer
Engineering, University of Macau, Macau
2Department of Electrical and Electronic Engineering, The University of Hong
Kong, Hong Kong  Fen Hou1 1Department of Electrical and Computer Engineering,
University of Macau, Macau
2Department of Electrical and Electronic Engineering, The University of Hong
Kong, Hong Kong  Kam-Weng Tam1 1Department of Electrical and Computer
Engineering, University of Macau, Macau
2Department of Electrical and Electronic Engineering, The University of Hong
Kong, Hong Kong  and Yik-Chung Wu2 1Department of Electrical and Computer
Engineering, University of Macau, Macau
2Department of Electrical and Electronic Engineering, The University of Hong
Kong, Hong Kong

###### Abstract

This paper investigates the problem of power allocation for hybrid automatic
repeat request (HARQ) schemes over time-correlated Nakagami-m fading channels
under outage constraint. The presence of time correlation complicates the
power allocation problem due to the involvement of multiple correlated fading
channels. Under a general time-correlated Nakagami-m fading channel with
exponential correlation, outage probabilities for three widely adopted HARQ
schemes, including Type I HARQ, HARQ with chase combining (HARQ-CC) and HARQ
with incremental redundancy (HARQ-IR), are first derived. With these results,
power allocation schemes are proposed to minimize the average total
transmission power with guaranteed outage performance. Simulation results
demonstrate the accuracy of our outage analysis and the effectiveness of our
proposed power allocation schemes. It is shown that our proposed power
allocation schemes can achieve significant power savings when compared with
fixed power allocation. Moreover, under practical low outage constraint, the
power efficiency is further improved when the time correlation is reduced
and/or the fading order is increased.

###### Index Terms:

Time-correlated Nakagami-m fading, hybrid automatic repeat request, power
allocation.

##  I Introduction

Hybrid automatic repeat request (HARQ) is a powerful transmission protocol to
combat the detrimental effects of channel fading and noise due to its
combination of automatic repeat request and forward error correction.
Generally, there are three types of HARQ schemes, including Type I HARQ, HARQ
with chase combining (HARQ-CC) and HARQ with incremental redundancy (HARQ-IR).
For Type I HARQ, the erroneously received packets are discarded and only the
most recently received packet is used for decoding. Since the failed packet
may still contain some useful information, it can be exploited for performance
enhancement and the other two HARQ schemes are thus designed for this purpose.
They combine the erroneously received packets with subsequently received
packets for joint decoding to improve the performance. Their difference lies
in whether the same set of coded bits are transmitted in each HARQ round.
Specifically, for HARQ-CC, the same coded sequence is repetitively transmitted
in each HARQ round and maximal-ratio-combining (MRC) is employed to combine
all the received packets to recover the message, whereas HARQ-IR transmits
different sets of coded bits in each retransmission and code combining is
adopted for joint decoding.

Power allocation for HARQ schemes has attracted considerable research
attention recently. Most of prior works consider either quasi-static fading
channels [1, 2, 3] or fast fading channels [4, 5, 6]. To be specific, in [1],
an optimal power allocation scheme is proposed to minimize the average total
transmission power of HARQ-CC over quasi-static fading channels, where the
channel response remains constant during multiple HARQ rounds. Similar to [1],
outage-limited power allocation is investigated for HARQ-CC and HARQ-IR
schemes in both continuous and bursting communication systems in [2].
Considering the same quasi-static fading channels, power allocation is
investigated in [3]. A backward sequential calculation method is developed to
find the optimum power allocation. On the other hand, some of prior literature
considers fast fading channels, where channel responses vary independently
among multiple transmissions. For example, in [4], power allocation is
discussed for HARQ-IR enabled distributed cooperative beamforming system,
where the source and the relay have fixed transmission power in each HARQ
round. Another power allocation scheme is proposed for HARQ-CC over
independent Rayleigh fading channels in [5]. By reformulating the power
allocation problem as a geometric programming problem and using dominant-term
approximation, the optimal solution is found efficiently. The same approach is
further extended to the power allocation for HARQ-enabled incremental MIMO
systems in [6].

Apart from quasi-static and fast fading channels, another frequently
experienced channel is time-correlated fading channel [7, 8], which usually
occurs when the transceiver has low-to-medium mobility. Under time correlated
fading channels, power allocation becomes much more challenging due to the
involvement of multiple correlated random variables and there are few
solutions if any in the literature. In this paper, we investigate power
allocation for HARQ schemes over time-correlated Nakagami-m fading channels. A
general multivariate Nakagami-m distribution with exponential correlation is
adopted to model time-correlated fading channels. The outage probabilities and
their asymptotic expressions are first derived for three HARQ schemes, i.e.,
Type I HARQ, HARQ-CC and HARQ-IR. These analytical results then enable the
optimal power allocation to minimize the average total transmission power with
guaranteed outage performance. Closed-form optimal solutions are found based
on the asymptotic outage expressions. Finally, these theoretical results are
validated through simulations. It is found that our proposed power allocation
schemes can achieve significant power savings when compared with fixed power
allocation. Moreover, under practical low outage constraint, the power
efficiency is further improved with the reduction of time correlation and the
increase of fading order.

The remainder of this paper is organized as follows. In Section II, system
model is given and outage analysis is conducted for three HARQ schemes.
Section III generalizes the problem of outage-limited power allocation for
three HARQ schemes, and optimal solutions are proposed in closed forms. In
Section IV, numerical results are presented and discussed to demonstrate the
efficiency of our proposed power allocation schemes. Finally, Section V
concludes this paper.

##  II System Model and Outage Analysis

A point-to-point HARQ enabled system operating over time-correlated Nakagami-m
block-fading channels is considered in this paper. Following the HARQ
protocol, $L$ maximal transmissions are allowed for each single message. The
received signal $y_{l}$ in the $l$th HARQ round is written as

$y_{l}=\sqrt{P_{l}}h_{l}x_{l}+\eta_{l},\quad 0\leq l\leq L,$ (1)

where $x_{l}$ denotes the transmitted signal with unit power in the $l$th HARQ
round, $P_{l}$ refers to the transmit power in the $l$th HARQ round,
$\eta_{l}$ represents the complex Gaussian white noise with zero mean and unit
variance, i.e., $\eta_{l}\sim{\mathcal{CN}}(0,1)$, and $h_{l}$ is the channel
coefficient in the $l$th HARQ round. Unlike prior literature, time correlated
Nakagami-m fading channels are considered. More precisely, the joint
distribution of channel amplitudes $|{\bf{h}}_{L}|=[|h_{1}|,\cdots,|h_{L}|]$
are modeled as a multivariate Nakagami-m distribution with exponential
correlation [7, 8, 9], whose joint probability density function (PDF) is given
by

${f_{|{\bf{h}}_{L}|}}\left({{z_{1}},\cdots,{z_{L}}}\right)=\int_{t=0}^{\infty}{\frac{{{t^{m-1}}}}{{\Gamma\left(m\right)}}{{\rm{e}}^{-t}}}\times\\\
\prod\limits_{l=1}^{L}{\frac{{2{z_{l}}^{2m-1}}}{{\Gamma\left(m\right){{\left({\frac{{{\Omega_{l}}\left({1-{\rho^{2\left({l+\delta-1}\right)}}}\right)}}{m}}\right)}^{m}}}}{e^{-\frac{{m{z_{l}}^{2}}}{{{\Omega_{l}}\left({1-{\rho^{2\left({l+\delta-1}\right)}}}\right)}}}}}\times\\\
{e^{-\frac{{{\rho^{2\left({l+\delta-1}\right)}}t}}{{1-{\rho^{2\left({l+\delta-1}\right)}}}}}}_{0}{F_{1}}\left({;m;\frac{{m{z_{l}}^{2}{\rho^{2\left({l+\delta-1}\right)}}t}}{{{\Omega_{l}}{{\left({1-{\rho^{2\left({l+\delta-1}\right)}}}\right)}^{2}}}}}\right)dt,\rho\neq
1,$ (2)

where $\rho$ and $\delta$ denote the time correlation and the channel feedback
delay, $m$ denotes the fading order, ${{\Omega_{l}}}$ is defined as the
average power of $h_{l}$, i.e., ${{\Omega_{l}}}={\rm E}\\{|h_{l}|^{2}\\}$,
$\Gamma(\cdot)$ denotes Gamma function and ${}_{0}{F_{1}}(\cdot)$ denotes the
confluent hypergeometric limit function [10, Eq. 9.14.1].

The system performance is fully characterized by outage probability, which is
defined as the probability that the message cannot be successfully decoded,
i.e., the mutual information is smaller than the target transmission rate $R$
bps/Hz. For different HARQ schemes, the outage probability over time-
correlated Nakagami-m fading channels are analyzed as follows.

###  II-A Outage Probability of Type I HARQ

For Type I HARQ, only the most recently received packet is employed for
recovering the message. The outage probability $p_{out,l}^{\rm I}$ after $l$
transmissions can be formulated as

$p_{out,l}^{\rm{I}}=\Pr\left({{I_{1}}<R,\cdots,{I_{l}}<R}\right)\\\
={F_{|{{\bf{h}}_{l}}|}}\left({\left|{{h_{1}}}\right|<\sqrt{\frac{{{2^{R}}-1}}{{{P_{1}}}}},\cdots,\left|{{h_{l}}}\right|<\sqrt{\frac{{{2^{R}}-1}}{{{P_{l}}}}}}\right),$
(3)

where
${I_{\iota}}={\log_{2}}\left({1+{P_{\iota}}{{\left|{{h_{\iota}}}\right|}^{2}}}\right)$
denotes the mutual information in the $\iota$th transmission, and
${F_{|{{\bf{h}}_{l}}|}}(\cdot)$ denotes the joint cumulative distribution
function (CDF) with respect to $|{{\bf{h}}_{l}}|$, which can be derived in the
following lemma.

######  Lemma 1.

The joint CDF ${F_{|{{\bf{h}}_{l}}|}}(y_{1},\cdots,y_{l})$ can be written as a
weighted sum of joint CDF of $l$ independent Nakagami RVs ${\bf A}_{\bf n}$
with parameters
$(m+n_{\iota},{{{\Omega_{\iota}}\left({1-{\rho^{2\left({\iota+\delta-1}\right)}}}\right)}}(m+n_{\iota})/m)$,
where ${\bf n}=[n_{1},\cdots,n_{l}]$ and $0\leq\iota\leq l$. Precisely,

${F_{|{{\bf{h}}_{l}}|}}\left({{y_{1}},\cdots,{y_{l}}}\right)=\sum\limits_{{n_{1}},\cdots,{n_{l}}=0}^{\infty}{{W_{\bf{n}}}{F_{{{\bf{A}}_{\bf{n}}}}}\left({{y_{1}},\cdots,{y_{l}}}\right)},$
(4)

where the coefficient $W_{\bf n}$ is given by

${W_{\bf{n}}}=\frac{{\Gamma\left({m+\sum\limits_{\iota=1}^{l}{{n_{\iota}}}}\right)}}{{\Gamma\left(m\right){{\left({1+\sum\limits_{\iota=1}^{l}{{\omega_{\iota}}}}\right)}^{m}}}}\prod\limits_{\iota=1}^{l}{\frac{1}{{{n_{\iota}}!}}{{\left({\frac{{{\omega_{\iota}}}}{{1+\sum\limits_{\iota=1}^{l}{{\omega_{\iota}}}}}}\right)}^{{n_{\iota}}}}}$
(5)

and satisfies
$\sum\limits_{{n_{1}},\cdots,{n_{l}}=0}^{\infty}{{W_{\bf{n}}}}=1$,
${\omega_{\iota}}{\rm{=}}\frac{{{\rho^{2\left({\iota+\delta-1}\right)}}}}{{1-{\rho^{2\left({\iota+\delta-1}\right)}}}}$,
and the joint CDF with respect to ${\bf A}_{\bf n}$,
${F_{{{\bf{A}}_{\bf{n}}}}}\left({{y_{1}},\cdots,{y_{l}}}\right)$, is
explicitly expressed as

${F_{{{\bf{A}}_{\bf{n}}}}}\left({{y_{1}},\cdots,{y_{l}}}\right)=\prod\limits_{\iota=1}^{l}{\frac{{\Upsilon\left({m+{n_{\iota}},\frac{{m{y_{\iota}}^{2}}}{{{\Omega_{\iota}}\left({1-{\rho^{2\left({\iota+\delta-1}\right)}}}\right)}}}\right)}}{{\Gamma\left({m+{n_{\iota}}}\right)}}}$
(6)

with $\Upsilon(\cdot,\cdot)$ being the lower incomplete Gamma function.

###### Proof.

The result directly follows from (2) and the series expansion of
${}_{0}{F_{1}}(\cdot)$ [10, Eq. 9.14.1]. ∎

With Lemma 1, the outage probability of Type I HARQ can be obtained as

$p_{out,l}^{\rm{I}}=\sum\limits_{{n_{1}},\cdots,{n_{l}}=0}^{\infty}{{W_{\bf{n}}}{F_{{{\bf{A}}_{\bf{n}}}}}\left({\sqrt{\frac{{{2^{R}}-1}}{{{P_{1}}}}},\cdots,\sqrt{\frac{{{2^{R}}-1}}{{{P_{l}}}}}}\right)}.$
(7)

In practice, the outage probability can be computed by truncating the infinite
series in (7). Herein, an efficient truncation method is proposed as

$\tilde{p}_{out,l}^{\rm{I}}=\sum\limits_{t=0}^{N}{\sum\limits_{{n_{1}}{\rm{+}}\cdots{\rm{+}}{n_{l}}=t}^{\infty}{{W_{\bf{n}}}{F_{{{\bf{A}}_{\bf{n}}}}}\left({\sqrt{\frac{{{2^{R}}-1}}{{{P_{1}}}}},\cdots,\sqrt{\frac{{{2^{R}}-1}}{{{P_{l}}}}}}\right)}},$
(8)

where $N$ defines the truncation order. It can be proved that the truncation
error exponentially decreases with $N$. The proof is omitted here due to space
limit.

Under high SNR, the outage probability can be asymptotically derived as shown
in following theorem.

######  Theorem 1.

Under high SNR regime, i.e., $P_{1},\cdots,P_{l}\to\infty$, the outage
probability $p_{out,l}^{\rm{I}}$ is written as

$p_{out,l}^{\rm{I}}=\frac{{{m^{ml}}\ell\left(l,\rho\right){{\left({{2^{R}}-1}\right)}^{lm}}}}{{{\Gamma^{l}}\left({m+1}\right)\prod\limits_{\iota=1}^{l}{{\Omega_{\iota}}^{m}{P_{\iota}}^{m}}}},$
(9)

where
$\ell\left(l,\rho\right)={\left({\left({1+\sum\limits_{\iota=1}^{l}{\frac{{{\rho^{2\left({\iota+\delta-1}\right)}}}}{{1-{\rho^{2\left({\iota+\delta-1}\right)}}}}}}\right)\prod\limits_{\iota=1}^{l}{\left({1-{\rho^{2\left({\iota+\delta-1}\right)}}}\right)}}\right)^{-m}}$.

###### Proof.

By using the series expansion of $\Upsilon(\cdot,\cdot)$ [10, Eq. 8.354.1] and
omitting the higher order infinitesimal of
${{{{\prod\limits_{\iota=1}^{l}{{P_{\iota}}}^{-m}}}}}$, the outage probability
(7) can be asymptotically expressed as (9). ∎

###  II-B Outage Probability of HARQ-CC

In HARQ-CC scheme, all the previously received packets are combined through
MRC for decoding. The outage probability after $l$ HARQ rounds is thus written
as

$\displaystyle p_{out,l}^{CC}$
$\displaystyle=\Pr\left({{{\log}_{2}}\left({1+\sum\limits_{\iota=1}^{l}{{P_{\iota}}{{\left|{{h_{\iota}}}\right|}^{2}}}}\right)<R}\right)$
$\displaystyle=\Pr\left({Y_{l}\triangleq\sum\limits_{\iota=1}^{l}{{P_{\iota}}{{\left|{{h_{\iota}}}\right|}^{2}}}<{2^{R}}-1}\right)={F_{{Y_{l}}}}\left(2^{R}-1\right).$
(10)

where ${F_{{Y_{l}}}}\left(\cdot\right)$ denotes the CDF of $Y_{l}$. After
deriving the CDF ${F_{{Y_{l}}}}\left(\cdot\right)$ using the method of moment
generating function (MGF), the outage probability $p_{out,l}^{CC}$ is derived
in the following theorem.

######  Theorem 2.

The outage probability for HARQ-CC scheme $p_{out,l}^{CC}$ can be obtained as

$p_{out,l}^{CC}=1+\frac{{{m^{ml}}\ell\left(l,\rho\right)}}{{\prod\nolimits_{\iota=1}^{l}{{\Omega_{\iota}}^{m}{P_{\iota}}^{m}}}}\times\\\
\sum\limits_{\kappa=1}^{\cal
K}{\sum\limits_{\varsigma=1}^{m{q_{\kappa}}}{\frac{{{\Phi_{\kappa\varsigma}}\left({-{\lambda_{\kappa}}}\right)}}{{\left({m{q_{\kappa}}-\varsigma}\right)!\left({\varsigma-1}\right)!}}{(2^{R}-1)^{m{q_{\kappa}}-\varsigma}}{e^{-{\lambda_{\kappa}}(2^{R}-1)}}}}$
(11)

where $\lambda_{1},\cdots,\lambda_{\mathcal{K}}$ define $\mathcal{K}$ distinct
poles of the MGF of $Y_{l}$ with multiplicities
$q_{1},\cdots,q_{\mathcal{K}}$, respectively,
$\sum\nolimits_{\kappa=1}^{\mathcal{K}}{{q_{\kappa}}}=l$, and
${\Phi_{\kappa\varsigma}}\left(s\right)=\frac{{{d^{\varsigma-1}}}}{{d{s^{\varsigma-1}}}}\left({{s^{-1}}\prod\limits_{\tau=1,\tau\neq\kappa}^{\cal
K}{{{\left({s{\rm{+}}{\lambda_{\tau}}}\right)}^{-m{q_{\tau}}}}}}\right)$.
Under high SNR regime, the outage probability $p_{out,l}^{CC}$ can also be
expressed asymptotically as

$p_{out,l}^{CC}=\frac{{{m^{ml}}\ell\left(l,\rho\right){{\left({{2^{R}}-1}\right)}^{ml}}}}{{\Gamma\left({ml+1}\right)\prod\nolimits_{\iota=1}^{l}{{\Omega_{\iota}}^{m}{P_{\iota}}^{m}}}}.$
(12)

###### Proof.

Please see Appendix A. ∎

###  II-C Outage Probability of HARQ-IR

Different from Type I HARQ and HARQ-CC, HARQ-IR accumulates mutual information
in all previous HARQ rounds for decoding. From information theoretical
perspective, an outage happens when the accumulated mutual information is less
than the target transmission rate $R$. Thus the outage probability after $l$
HARQ rounds is formulated as

$p_{out,l}^{IR}=\Pr\left({\sum\limits_{\iota=1}^{l}{{{\log}_{2}}\left({1+{P_{\iota}}{{\left|{{h_{\iota}}}\right|}^{2}}}\right)}<R}\right).$
(13)

Due to the time correlation among $h_{l}$, it is intractable to find closed-
form expression for (13). Instead, a lower bound of $p_{out,l}^{IR}$ is
adopted to characterize the outage probability of HARQ-IR. By using Jensen’s
inequality, $p_{out,l}^{IR}$ is lower bounded as

$\displaystyle p_{out,l}^{IR}$
$\displaystyle\geq\Pr\left({{{\log}_{2}}\left({1+\frac{1}{l}\sum\limits_{\iota=1}^{l}{{P_{\iota}}{{\left|{{h_{\iota}}}\right|}^{2}}}}\right)<\frac{R}{l}}\right)$
$\displaystyle={F_{{Y_{l}}}}\left({l\left({{2^{R/l}}-1}\right)}\right)\triangleq{p_{out,l}^{IR,lower}}.$
(14)

With the CDF ${F_{{Y_{l}}}}\left(\cdot\right)$ derived in Theorem 2, the lower
bound ${p_{out,l}^{IR,lower}}$ and its asymptotic expression can be derived in
the following theorem.

######  Theorem 3.

The lower bound of the outage probability $p_{out,l}^{IR,lower}$ can be
obtained as

${p_{out,l}^{IR,lower}}=1+\frac{{{m^{ml}}\ell\left(l,\rho\right)}}{{\prod\nolimits_{\iota=1}^{l}{{\Omega_{\iota}}^{m}{P_{\iota}}^{m}}}}\times\\\
\sum\limits_{\kappa=1}^{\cal
K}{\sum\limits_{\varsigma=1}^{m{q_{\kappa}}}{\frac{{{\Phi_{\kappa\varsigma}}\left({-{\lambda_{\kappa}}}\right)}}{{\left({m{q_{\kappa}}-\varsigma}\right)!\left({\varsigma-1}\right)!}}{(l({{2^{R/l}}-1}))^{m{q_{\kappa}}-\varsigma}}{e^{-{\lambda_{\kappa}}(l({{2^{R/l}}-1}))}}}}.$
(15)

Under high SNR regime, ${p_{out,l}^{IR,lower}}$ is further simplified as

$p_{out,l}^{IR,lower}=\frac{{{(lm)^{ml}}\ell\left(l,\rho\right){{\left({{2^{R/l}}-1}\right)}^{ml}}}}{{\Gamma\left({ml+1}\right)\prod\nolimits_{\iota=1}^{l}{{\Omega_{\iota}}^{m}{P_{\iota}}^{m}}}}.$
(16)

##  III Optimal Power allocation

In this section, the problem of power allocation is studied for the three HARQ
schemes. Generally, the average total transmission power for HARQ is defined
as ${\bar{P}=\sum\nolimits_{l=1}^{L}{{P_{l}}{p_{out,l-1}}}}$ [5]. Here
${p_{out,l}}$ refers to the outage probability after $l$ transmissions and it
unifies the cases of ${p_{out,l}^{I}}$, ${p_{out,l}^{CC}}$ and
${p_{out,l}^{IR,lower}}$. When power efficiency is concerned with certain
performance requirement, the transmission power among multiple HARQ rounds
should be properly designed to minimize the total transmission power while
guaranteeing the performance. The power allocation problem can be formulated
as

$\begin{array}[]{*{20}{cl}}{\mathop{\min}\limits_{{P_{1}},{P_{2}},\cdots{P_{L}}}}&{\bar{P}=\sum\limits_{l=1}^{L}{{P_{l}}{p_{out,l-1}}}}\\\
{{\rm{s}}{\rm{.t}}{\rm{.}}}&{{P_{l}}\geq 0,0\leq l\leq L}\\\
{}\hfil&{{p_{out,L}}\leq\varepsilon}\\\ \end{array},$ (17)

where $\varepsilon$ represents the outage tolerance.

Due to the complicated expressions of the exact outage probabilities given in
(7), (11) and (15), it is impossible to find closed-form power allocation
solutions directly. However, interior-point methods can be exploited to
numerically solve the problem (17). Meanwhile, based on the asymptotic
expressions of the outage probabilities, an efficient power allocation scheme
can be found as follows.

Notice that the asymptotic outage probabilities in (9), (12) and (16) can be
unified as

${p_{out,l}}\simeq\frac{{{\phi_{l}}}}{{{{\left({\prod\limits_{k=1}^{l}{{P_{k}}}}\right)}^{m}}}},\,0\leq
l\leq L,$ (18)

where $\phi_{l}$ depends on HARQ schemes, more precisely,

${\phi_{l}}=\left\\{{\begin{array}[]{*{20}{cl}}{\frac{{{m^{ml}}\ell\left({l,\rho}\right){{\left({{2^{R}}-1}\right)}^{ml}}}}{{{\Gamma^{l}}\left({m+1}\right)\prod\nolimits_{\iota=1}^{l}{{\Omega_{\iota}}^{m}}}},}&{{\textrm{Type}}\;{\textrm{I}};}\\\
{\frac{{{m^{ml}}\ell\left({l,\rho}\right){{\left({{2^{R}}-1}\right)}^{ml}}}}{{\Gamma\left({ml+1}\right)\prod\nolimits_{\iota=1}^{l}{{\Omega_{\iota}}^{m}}}},}&{{\textrm{HARQ-
CC}};}\\\
{\frac{{{{\left({ml}\right)}^{ml}}\ell\left({l,\rho}\right){{\left({{2^{R/l}}-1}\right)}^{ml}}}}{{\Gamma\left({ml+1}\right)\prod\nolimits_{\iota=1}^{l}{{\Omega_{\iota}}^{m}}}},}&{{\textrm{HARQ-
IR}}.}\end{array}}\right.$ (19)

Substituting (18) into (17), the Lagrangian of the optimization problem (17)
is formed as

$\mathfrak{L}\left({{P_{1}},\cdots,{P_{L}},\mu,{\nu_{1}},\cdots,{\nu_{L}}}\right)=\sum\limits_{l=1}^{L}{{P_{l}}\frac{{{\phi_{l-1}}}}{{{{\left({\prod\limits_{k=1}^{l-1}{{P_{k}}}}\right)}^{m}}}}}\\\
+\mu\left({\frac{{{\phi_{L}}}}{{{{\left({\prod\limits_{k=1}^{L}{{P_{k}}}}\right)}^{m}}}}-\varepsilon}\right)-\sum\limits_{l=1}^{L}{{\nu_{l}}{P_{l}}},$
(20)

where $\mu,\nu_{1},\cdots,\nu_{L}$ are the Lagrangian multipliers of the
constraints in the problem (17). Since the Karush-Khun-Tucker (KKT) conditions
are necessary for an optimal solution, we have

${\left.{\frac{{\partial\mathfrak{L}}}{{\partial{P_{n}}}}}\right|_{\left({{P_{1}^{*}},\cdots,{P_{L}^{*}},{\mu^{*}},{\nu_{1}}^{*},\cdots,{\nu_{L}}^{*}}\right)}}=0,$
(21)
${\mu^{*}}\left({\frac{{{\phi_{L}}}}{{{{\left({\prod\limits_{k=1}^{L}{P_{k}^{*}}}\right)}^{m}}}}-\varepsilon}\right)=0,$
(22) ${\nu_{l}}^{*}{P_{l}}^{*}=0,$ (23)

where $\mu^{*},{\nu_{1}}^{*},\cdots,{\nu_{L}}^{*},{P_{l}}^{*}$ denote the
optimal Lagrangian multipliers and the optimal power allocation, respectively.
Based on the KKT conditions (21)-(23), the optimal power allocation solution
to (17) could be found in closed form as follows.

######  Theorem 4.

The optimal solution to the problem (17) is uniquely given as

$\displaystyle P_{L}^{*}$
$\displaystyle={\left({\frac{{{\phi_{L}}\prod\limits_{k=2}^{L}{{{\left({\left({m+1}\right)\frac{{{\phi_{k-1}}}}{{{\phi_{k-2}}}}}\right)}^{\frac{1}{{{{\left({m+1}\right)}^{k-1}}}}}}}}}{{{\phi_{L-1}}{{\left({m+1}\right)}^{L-1}}\varepsilon}}}\right)^{\frac{{{{\left({m+1}\right)}^{L-1}}}}{{{{\left({m+1}\right)}^{L}}-1}}}},$
(24)
$P_{n}^{*}=\prod\limits_{k=n+1}^{L}{{{\left({\left({m+1}\right)\frac{{{\phi_{k-1}}}}{{{\phi_{k-2}}}}}\right)}^{\frac{1}{{{{\left({m+1}\right)}^{k-n}}}}}}}{P_{L}^{*}}^{\frac{1}{{{{\left({m+1}\right)}^{L-n}}}}},\\\
\quad{\textrm{f}or}\quad 1\leq n\leq L-1.$ (25)

Moreover, the minimal average total transmission power $\bar{P}^{*}$ is

$\displaystyle\bar{P}^{*}$
$\displaystyle=\frac{{\varepsilon\left({P_{L}}{{}^{*}}\right)^{{m+1}}{\phi_{L-1}}}}{{{\phi_{L}}}}\frac{{{{{\left({m+1}\right)}^{L}}-1}}}{m}.$
(26)

###### Proof.

Please see Appendix B. ∎

##  IV Numerical Results and Discussions

In this section, the analytical results are verified through simulations, and
our proposed power allocation (PPA) scheme is compared to the fixed power
allocation (FPA) scheme in [4]. Notice that for FPA, the problem (17) is
solved by adding the constraint $P_{1}=\cdots=P_{L}$. In the sequel, we take
systems with $\Omega_{1}=\cdots=\Omega_{l}=1$, $\delta=1$s and $R=2$bps/Hz as
examples.

###  IV-A Comparison of PPA and FPA

The minimal total transmission powers $\bar{P}^{*}$ of the PPA and FPA schemes
are compared in Fig. 1. The outage-limited systems with various HARQ and
parameters as $L=2$, $m=2$ and $\rho=0.5$ are considered. Clearly, the results
of PPA using asymptotic outage expressions (PPA-A) agree well with that of PPA
using exact outage expressions (PPA-E) under low outage constraint
$\varepsilon$. It is also readily found that PPA-A performs better than FPA
under low outage constraint $\varepsilon$, and their performance gap
significantly increases when $\varepsilon$ decreases. Moreover, it reveals
that HARQ-IR is superior to Type I HARQ and HARQ-CC in terms of power
efficiency.

Figure 1: Comparison of the proposed power allocation with the fixed power
allocation.

###  IV-B Impacts of Time Correlation

Since the performance of PPA-A is asymptotically equivalent to that of PPA-E
under low $\varepsilon$, PPA-A is adopted to test the impact of time
correlation on power allocation. Fig. 2 plots the minimal total transmission
power $\bar{P}^{*}$ against time correlation coefficient $\rho$ by setting
parameters as $m=2$ and $\varepsilon=10^{-6}$. It can be easily seen that the
increase of time correlation $\rho$ would lead to the increase of the minimal
total transmission power $\bar{P}^{*}$ for both $L=2$ and $L=4$. It means that
time correlation has negative effect on power efficiency under low outage
constraint.

Figure 2: Impact of time correlation.

###  IV-C Impacts of Fading Order

Fig. 3 depicts the impact of fading order $m$ on power allocation by setting
$\rho=0.5$ and $\varepsilon=10^{-6}$. Clearly, the minimal total transmission
power $\bar{P}^{*}$ decreases with the increase of fading order $m$. In fact,
higher fading order leads to higher diversity introduced by the channel, thus
reducing the power consumption given a certain outage performance constraint.

Figure 3: Impact of fading order.

##  V Conclusions

Outage-limited power allocation for various HARQ schemes operating over time-
correlated Nakagami-m fading channels has been investigated in this paper.
After deriving the outage probabilities and their asymptotic expressions, an
optimal power allocation solution has been proposed in closed form. It has
been demonstrated that our proposed solution can achieve significant power
saving when compared to the fixed power allocation solution. The superiority
of the proposed optimal solution in terms of power consumption is further
enhanced when the channel time correlation is reduced and/or the fading order
increases.

##  VI Acknowledgements

This work was supported in part by the Research Committee of University of
Macau under grants: MYRG101(Y1-L3)-FST13-MSD, MYRG2014-00146-FST and
MYRG2016-00146-FST, in part by the Macau Science and Technology Development
Fund under grant FDCT 091/2015/A3, and in part by the National Natural Science
Foundation of China under grant No.61601524.

##  Appendix A Proof of Theorem 2

The moment generation function (MGF) with respect to $Y_{l}$ can be written as

${M_{{Y_{l}}}}\left(s\right)={\rm{E}}\left({{e^{s{Y_{l}}}}}\right)\\\
=\int\limits_{0}^{\infty}{\cdots\int\limits_{0}^{\infty}{{e^{s\sum\limits_{\iota=1}^{l}{{P_{\iota}}{x_{\iota}}^{2}}}}{f_{|{{\bf{h}}_{l}}|}}\left({{x_{1}},\cdots,{x_{l}}}\right)d{x_{1}}\cdots{dx_{l}}}}.$
(27)

Plugging (2) into (27), after some algebraic manipulations, it follows that

$\displaystyle{M_{{Y_{l}}}}\left(s\right)$
$\displaystyle={{{\left(\begin{array}[]{l}\prod\limits_{\iota=1}^{l}{\left({1-s{P_{\iota}}{\Omega_{\iota}}\left({1-{\rho^{2\left({\iota+\delta-1}\right)}}}\right)/m}\right)}\\\
\times\left({1+\sum\limits_{\iota=1}^{l}{\frac{{s{P_{\iota}}{\Omega_{\iota}}{\rho^{2\left({\iota+\delta-1}\right)}}/m}}{{s{P_{\iota}}{\Omega_{\iota}}\left({1-{\rho^{2\left({\iota+\delta-1}\right)}}}\right)/m-1}}}}\right)\end{array}\right)}^{-m}}}$
(30)
$\displaystyle=\frac{{{m^{ml}}\ell\left(l,\rho\right)}}{{\prod\limits_{\iota=1}^{l}{{\Omega_{\iota}}^{m}{P_{\iota}}^{m}}}}\prod\limits_{\kappa=1}^{\cal
K}{{{\left({{\lambda_{\kappa}}-s}\right)}^{-m{q_{\kappa}}}}},$ (31)

where $\lambda_{1},\cdots,\lambda_{\mathcal{K}}$ define $\mathcal{K}$ distinct
poles of ${M_{{Y_{l}}}}\left(s\right)$ with multiplicities
$q_{1},\cdots,q_{\mathcal{K}}$, respectively, and
$\sum\limits_{\kappa=1}^{\mathcal{K}}{{q_{\kappa}}}=l$. After some tedious
manipulations, ${M_{{Y_{l}}}}\left(s\right)$ can be simplified as
${M_{{Y_{l}}}}\left(s\right)=\det{\left({{\bf{I}}-s{{\bf{F}}^{1/2}{\bf{E}}{\bf{F}}^{1/2}}}\right)^{-m}}$,
where the notation $\det(\cdot)$ refers to the determinant, $\bf I$ represents
an identity matrix, $\bf F$ is an $l\times l$ diagonal matrix with diagonal
entries as $\\{\Omega_{\iota}P_{\iota}/m\\}_{\iota=1}^{l}$, and $\bf E$ is an
$l\times l$ positive definite matrix given by

${\bf{E}}=\left[{\begin{array}[]{*{20}{c}}1&{{\rho^{2\delta+1}}}&\cdots&{{\rho^{2\delta+l-1}}}\\\
{{\rho^{2\delta+1}}}&1&\cdots&{{\rho^{2\delta+l}}}\\\
\vdots&\vdots&\ddots&\vdots\\\
{{\rho^{2\delta+l-1}}}&{{\rho^{2\delta+l}}}&\cdots&1\end{array}}\right].$ (32)

Since $1/\lambda_{1},\cdots,1/\lambda_{\mathcal{K}}$ are the eigenvalues of
the positive definite matrix ${{\bf{F}}^{1/2}{\bf{E}}{\bf{F}}^{1/2}}$, we have
$\lambda_{1},\cdots,\lambda_{\mathcal{K}}>0$. By using inverse Laplace
transform and its integration property [11, Eq.9.109], the CDF with respect to
$Y_{l}$ is derived as

$\displaystyle{F_{{Y_{l}}}}\left(y\right)=\mathcal{L}^{-1}\left\\{{{M_{{Y_{l}}}}\left({-s}\right)}\right\\}\left(y\right)=\frac{1}{{2\pi
j}}\int\limits_{a-j\infty}^{a+j\infty}{\frac{{{M_{{Y_{l}}}}\left({-s}\right)}}{s}{e^{sy}}ds},$
(33)

By using [12, Eq. 5.2.21], (33) can be calculated as

${F_{{Y_{l}}}}\left(y\right)=\frac{{{m^{ml}}\ell\left(l,\rho\right)}}{{\prod\limits_{\iota=1}^{l}{{\Omega_{\iota}}^{m}{P_{\iota}}^{m}}}}\frac{1}{{2\pi
j}}\int\limits_{a-j\infty}^{a+j\infty}{\frac{{{e^{sy}}}}{{s\prod\limits_{\kappa=1}^{\cal
K}{{{\left({s{\rm{+}}{\lambda_{\kappa}}}\right)}^{m{q_{\kappa}}}}}}}ds}\\\
=1+\frac{{{m^{ml}}\ell\left(l,\rho\right)}}{{\prod\limits_{\iota=1}^{l}{{\Omega_{\iota}}^{m}{P_{\iota}}^{m}}}}\sum\limits_{\kappa=1}^{\cal
K}{\sum\limits_{\varsigma=1}^{m{q_{\kappa}}}{\frac{{{\Phi_{\kappa\varsigma}}\left({-{\lambda_{\kappa}}}\right)}}{{\left({m{q_{\kappa}}-\varsigma}\right)!\left({\varsigma-1}\right)!}}{y^{m{q_{\kappa}}-\varsigma}}{e^{-{\lambda_{\kappa}}y}}}},$
(34)

where
${\Phi_{\kappa\varsigma}}\left(s\right)=\frac{{{d^{\varsigma-1}}}}{{d{s^{\varsigma-1}}}}\left({{s^{-1}}\prod\limits_{\tau=1,\tau\neq\kappa}^{\cal
K}{{{\left({s{\rm{+}}{\lambda_{\tau}}}\right)}^{-m{q_{\tau}}}}}}\right)$.
Therefore, by using (34) and together with $y=2^{R}-1$, (11) in Theorem 2 is
proved.

By using the expansion of Maclaurin series, (34) can be further rewritten as

${F_{{Y_{l}}}}\left(y\right)=\sum\limits_{n=0}^{\infty}{\frac{{{F_{{Y_{l}}}}^{\left(n\right)}\left(0\right)}}{{n!}}{y^{n}}}.$
(35)

Since

${F_{{Y_{l}}}}^{\left(n\right)}\left(0\right)=\frac{{{m^{ml}}\ell\left(l,\rho\right)}}{{\prod\limits_{\iota=1}^{l}{{\Omega_{\iota}}^{m}{P_{\iota}}^{m}}}}\frac{1}{{2\pi
j}}\int\limits_{a-j\infty}^{a+j\infty}{\frac{{{s^{n-1}}}}{{\prod\limits_{\kappa=1}^{\cal
K}{{{\left({s{\rm{+}}{\lambda_{\kappa}}}\right)}^{m{q_{\kappa}}}}}}}ds},\\\
0\leq n\leq ml,$ (36)

it follows by using the initial-value theorems of Laplace transform [11, Eq.
9.5.10] that
${F_{{Y_{l}}}}^{\left(1\right)}\left(0\right)=\cdots={F_{{Y_{l}}}}^{\left(ml-1\right)}\left(0\right)=0$,
and
${F_{{Y_{l}}}}^{\left(ml\right)}\left(0\right)=\frac{{{m^{ml}}\ell\left(l,\rho\right)}}{{\prod\nolimits_{\iota=1}^{l}{{\Omega_{\iota}}^{m}{P_{\iota}}^{m}}}}$.
Moreover, it can be proved that ${F_{{Y_{l}}}}^{\left(n\right)}\left(0\right)$
is a higher order term of $\prod\nolimits_{\iota=1}^{l}{{P_{\iota}}^{-m}}$ for
$n>ml$. The proof is omitted due to space limit. Thus it yields

${F_{{Y_{l}}}}\left(y\right)=\frac{{{m^{ml}}\ell\left(l,\rho\right){y^{ml}}}}{{\Gamma\left({ml+1}\right)\prod\limits_{\iota=1}^{l}{{\Omega_{\iota}}^{m}{P_{\iota}}^{m}}}}+o\left({\prod\limits_{\iota=1}^{l}{{P_{\iota}}^{-m}}}\right).$
(37)

Hence under high SNR, i.e., $P_{\iota}\to\infty$, (12) can be derived by using
(37) together with $y=2^{R}-1$. Thus Theorem 2 is proved.

##  Appendix B Proof of Theorem 4

Clearly from (18), since $P_{l}^{*}\neq 0$, we have ${\nu_{l}}=0$. Therefore,
after some algebraic manipulations, (21) can be rewritten as

${\left.{\frac{{\partial\mathfrak{L}}}{{\partial{P_{n}}}}}\right|_{\left({{P_{1}^{*}},\cdots,{P_{L}^{*}},{\mu^{*}}}\right)}}=\frac{{{\phi_{n-1}}}}{{{{\left({\prod\limits_{k=1}^{n-1}{P_{k}^{*}}}\right)}^{m}}}}\\\
-\frac{m}{{{P_{n}^{*}}}}\sum\limits_{l=n+1}^{L}{{P_{l}}^{*}\frac{{{\phi_{l-1}}}}{{{{\left({\prod\limits_{k=1}^{l-1}{P_{k}^{*}}}\right)}^{m}}}}}-\frac{m}{{{P_{n}^{*}}}}{\mu^{*}}\frac{{{\phi_{L}}}}{{{{\left({\prod\limits_{k=1}^{L}{P_{k}^{*}}}\right)}^{m}}}}=0.$
(38)

Together with
${\left.{\frac{{\partial\mathfrak{L}}}{{\partial{P_{n-1}}}}}\right|_{\left({{P_{1}^{*}},\cdots,{P_{L}^{*}},{\mu^{*}}}\right)}}=0$,
(38) reduces to

$P_{n}^{*}=\left({m+1}\right){P_{n+1}^{*}}\frac{{{\phi_{n}}}}{{{\phi_{n-1}}{{\left({P_{n}^{*}}\right)}^{m}}}}.$
(39)

Now from (39), we can derive $P_{n}^{*}$ recursively as (25). Regarding to
$P_{L}^{*}$, by letting $n=L$ in (38), we have

$\frac{{{\phi_{L-1}}}}{{{{\left({\prod\limits_{k=1}^{L-1}{P_{k}^{*}}}\right)}^{m}}}}=m{\mu^{*}}\frac{{{\phi_{L}}}}{{{P_{L}^{*}}{{\left({\prod\limits_{k=1}^{L}{P_{k}^{*}}}\right)}^{m}}}}\Rightarrow{\mu^{*}}{\rm{=}}\frac{{{P_{L}^{*}}^{m+1}{\phi_{L-1}}}}{{m{\phi_{L}}}}.$
(40)

Recalling that $P_{L}\neq 0$, thus ${\mu^{*}}\neq 0$. According to (22), we
have

$\frac{{{\phi_{L}}}}{{{{\left({\prod\limits_{k=1}^{L}{P_{k}^{*}}}\right)}^{m}}}}-\varepsilon=0\Rightarrow{\left({\prod\limits_{n=1}^{L}{{P_{n}^{*}}}}\right)^{m}}=\frac{{{\phi_{L}}}}{\varepsilon}.$
(41)

Substituting (25) into (41) yields (24). Moreover, by using (39), it follows
that

${P_{n}^{*}}=\frac{{{{\left({m+1}\right)}^{L-n}}{\phi_{L-1}}}}{{{\phi_{n-1}}}}\frac{{{P_{L}^{*}}}}{{{{\left({{P_{n}^{*}}\cdots
P_{L-1}^{*}}\right)}^{m}}}}.$ (42)

## References

  * [1] W. Su, S. Lee, D. Pados, J. D. Matyjas _et al._ , “Optimal power assignment for minimizing the average total transmission power in hybrid-ARQ Rayleigh fading links,” _IEEE Trans. Commun._ , vol. 59, no. 7, pp. 1867–1877, Jul. 2011.
  * [2] B. Makki, A. Graell i Amat, and T. Eriksson, “Green communication via power-optimized HARQ protocols,” _IEEE Trans. Veh. Technol._ , vol. 63, no. 1, pp. 161–177, Jan. 2013.
  * [3] G. Wang, J. Wu, and Y. R. Zheng, “Optimum energy and spectral efficient transmissions for delay-constrained hybrid ARQ systems,” _IEEE Trans. Veh. Technol._ , vol. PP, no. 99, pp. 1–1, Aug. 2014.
  * [4] J. Choi, W. Xing, D. To, Y. Wu, and S. Xu, “On the energy efficiency of a relaying protocol with HARQ-IR and distributed cooperative beamforming,” _IEEE Trans. Wireless Commun._ , vol. 12, no. 2, pp. 769–781, Feb. 2013\.
  * [5] T. V. Chaitanya and E. G. Larsson, “Optimal power allocation for hybrid ARQ with chase combining in i.i.d. Rayleigh fading channels,” _IEEE Trans. Commun._ , vol. 61, no. 5, pp. 1835–1846, May 2013.
  * [6] T. Chaitanya and T. Le-Ngoc, “Energy-efficient adaptive power allocation for incremental MIMO systems,” _IEEE Trans. Veh. Technol._ , vol. PP, no. 99, pp. 1–1, Mar. 2015.
  * [7] S. M. Kim, W. Choi, T. W. Ban, and D. K. Sung, “Optimal rate adaptation for hybrid ARQ in time-correlated Rayleigh fading channels,” _IEEE Trans. Wireless Commun._ , vol. 10, no. 3, pp. 968–979, Mar. 2011.
  * [8] H. Jin, C. Cho, N.-O. Song, and D. K. Sung, “Optimal rate selection for persistent scheduling with HARQ in time-correlated Nakagami-m fading channels,” _IEEE Trans. Wireless Commun._ , vol. 10, no. 2, pp. 637–647, Feb. 2011.
  * [9] N. C. Beaulieu and K. T. Hemachandra, “Novel simple representations for Gaussian class multivariate distributions with generalized correlation,” _IEEE Trans. Inf. Theory_ , vol. 57, no. 12, pp. 8072–8083, Dec. 2011.
  * [10] I. S. Gradshteyn, I. M. Ryzhik, A. Jeffrey, D. Zwillinger, and S. Technica, _Table of integrals, series, and products_. Academic press New York, 1965, vol. 6.
  * [11] A. V. Oppenheim, A. S. Willsky, and S. H. Nawab, _Signals &Amp; Systems (2Nd Ed.)_. Upper Saddle River, NJ, USA: Prentice-Hall, Inc., 1996.
  * [12] H. Bateman, “Tables of integral transforms,” _California Institute of Technology Bateman Manuscript Project, New York: McGraw-Hill, 1954, edited by Erdelyi, Arthur_ , vol. 1, 1954.