%% BioMed_Central_Tex_Template_v1.06
%%                                      %
%  bmc_article.tex            ver: 1.06 %
%                                       %

%%IMPORTANT: do not delete the first line of this template
%%It must be present to enable the BMC Submission system to
%%recognise this template!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                     %%
%%  LaTeX template for BioMed Central  %%
%%     journal article submissions     %%
%%                                     %%
%%          <8 June 2012>              %%
%%                                     %%
%%                                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                 %%
%% For instructions on how to fill out this Tex template           %%
%% document please refer to Readme.html and the instructions for   %%
%% authors page on the biomed central website                      %%
%% http://www.biomedcentral.com/info/authors/                      %%
%%                                                                 %%
%% Please do not use \input{...} to include other tex files.       %%
%% Submit your LaTeX manuscript as one .tex document.              %%
%%                                                                 %%
%% All additional figures and files should be attached             %%
%% separately and not embedded in the \TeX\ document itself.       %%
%%                                                                 %%
%% BioMed Central currently use the MikTex distribution of         %%
%% TeX for Windows) of TeX and LaTeX.  This is available from      %%
%% http://www.miktex.org                                           %%
%%                                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% additional documentclass options:
%  [doublespacing]
%  [linenumbers]   - put the line numbers on margins

%%% loading packages, author definitions

%\documentclass[twocolumn]{bmcart}% uncomment this for twocolumn layout and comment line below
\documentclass{bmcart}

%%% Load packages
\usepackage{amsthm,amsmath}
%\RequirePackage{natbib}
%\RequirePackage[authoryear]{natbib}% uncomment this for author-year bibliography
%\RequirePackage{hyperref}
\usepackage[utf8]{inputenc} %unicode support
%\usepackage[applemac]{inputenc} %applemac support if unicode package fails
%\usepackage[latin1]{inputenc} %UNIX support if unicode package fails
\usepackage[dvips]{graphicx}
%\usepackage{graphicx}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                             %%
%%  If you wish to display your graphics for   %%
%%  your own use using includegraphic or       %%
%%  includegraphics, then comment out the      %%
%%  following two lines of code.               %%
%%  NB: These line *must* be included when     %%
%%  submitting to BMC.                         %%
%%  All figure files must be submitted as      %%
%%  separate graphics through the BMC          %%
%%  submission process, not included in the    %%
%%  submitted article.                         %%
%%                                             %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%\def\includegraphic{}
%\def\includegraphics{}



%%% Put your definitions there:
\startlocaldefs
\endlocaldefs


%%% Begin ...
\begin{document}

%%% Start of article front matter
\begin{frontmatter}

\begin{fmbox}
\dochead{Research}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the title of your article here     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{On the Performance of Cooperative Non-orthogonal Multiple Access With Energy Harvesting Over Nakagami-$m$ Fading Channels}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors here                   %%
%%                                          %%
%% Specify information, if available,       %%
%% in the form:                             %%
%%   <key>={<id1>,<id2>}                    %%
%%   <key>=                                 %%
%% Comment or delete the keys which are     %%
%% not used. Repeat \author command as much %%
%% as required.                             %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\author[
   addressref={aff1,aff2},                   % id's of addresses, e.g. {aff1,aff2}
   corref={aff1},                       % id of corresponding address, if any
   %noteref={n1},                        % id's of article notes, if any
   email={gud2001@qfnu.edu.cn}   % email address
]{\inits{JE}\fnm{Weidong} \snm{Guo}}
%\author[
%   addressref={aff1,aff3},
%   email={chaozhai@sdu.edu.cn}
%]{\inits{JRS}\fnm{Chao} \snm{Zhai}}
\author[
   addressref={aff2},
   email={zcbsunny@hrbeu.edu.cn}
]{\inits{JRS}\fnm{Lin} \snm{Sun}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors' addresses here        %%
%%                                          %%
%% Repeat \address commands as much as      %%
%% required.                                %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\address[id=aff1]{%                           % unique id
  \orgname{School of Information Science and Engineering, Shandong University}, % university, etc
  \street{72 Binhai Road},                     %
  \postcode{266237}                                % post or zip code
  \city{Qingdao},                              % city
  \cny{China}                                    % country
}
\address[id=aff2]{%
  \orgname{School of Physics and Physical Engineering, Qufu Normal University},
  \street{57 Jingxuanxi Road},
  \postcode{273165}
  \city{Qufu},
  \cny{China}
}
%\address[id=aff3]{%
%  \orgname{National Mobile Communications Research Laboratory, Southeast University},
%  \street{2 Sipailou Road},
%  \postcode{210096}
%  \city{Nanjing},
%  \cny{China}
%}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter short notes here                   %%
%%                                          %%
%% Short notes will be after addresses      %%
%% on first page.                           %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{artnotes}
%\note{Sample of title note}     % note to the article
%\note[id=n1]{Equal contributor} % note, connected to author
\end{artnotes}

\end{fmbox}% comment this for two column layout

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Abstract begins here                 %%
%%                                          %%
%% Please refer to the Instructions for     %%
%% authors on http://www.biomedcentral.com  %%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstractbox}

\begin{abstract} % abstract
In this paper, we analyze the performance of a full-duplex (FD) cooperative non-orthogonal multiple access (NOMA) model in multi-relay multiuser networks. In contrast to conventional cooperative networks, relays in the considered network have no embedded energy supply, but they can harvest energy from the radio-frequency signals broadcasted by the source to facilitate the cooperative NOMA transmission. Based on this structure, we propose a novel relay selection scheme in Nakagami-\emph{m} fading channels by considering both the channel state information (CSI) and the battery statuses of relays. Assuming each relay has finite energy storage for the energy accumulation, we adpot a finite Markov chain to capture the evolution of relay batteries and simplify the performance analysis by making some reasonable assumptions. General closed-form expressions of the outage probability and the ergodic sum-rate are derived. All the theoretical results are validated by Monte-Carlo simulations. The impacts of various system parameters, such as the number of relays and the self-interference (SI) at the involved relay, on the performance are extensively investigated. It is shown that the FD relaying achieves better performance in the moderate and high SNR regions and the usage of NOMA with FD relaying outperforms the conventional orthogonal multiple access (OMA) network in terms of terminal fairness.
\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The keywords begin here                  %%
%%                                          %%
%% Put each keyword in separate \kwd{}.     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{keyword}
\kwd{Non-orthogonal multiple access}
\kwd{Full-duplex}
\kwd{Energy harvesting}
\kwd{Markov chain}
\end{keyword}

% MSC classifications codes, if any
%\begin{keyword}[class=AMS]
%\kwd[Primary ]{}
%\kwd{}
%\kwd[; secondary ]{}
%\end{keyword}

\end{abstractbox}
%
%\end{fmbox}% uncomment this for twcolumn layout

\end{frontmatter}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Main Body begins here                %%
%%                                          %%
%% Please refer to the instructions for     %%
%% authors on:                              %%
%% http://www.biomedcentral.com/info/authors%%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%% See the Results and Discussion section   %%
%% for details on how to create sub-sections%%
%%                                          %%
%% use \cite{...} to cite references        %%
%%  \cite{koon} and                         %%
%%  \cite{oreg,khar,zvai,xjon,schn,pond}    %%
%%  \nocite{smith,marg,hunn,advi,koha,mouse}%%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%% start of article main body
% <put your article body there>

%%%%%%%%%%%%%%%%
%% Background %%
%%
\section{Introduction}
With the aim of increasing the capacity of wireless systems, non-orthogonal multiple access (NOMA) has been recognized as a promising technique for the fifth generation (5G) networks due to its superior spectral efficiency [1]. The key idea of NOMA is serving multiple users in the same time solt over the same frequency band, but with different power levels to achieve a tradeoff between system throughput and user fairness [2]. Lan et al. [3], explored the impact of error propagation of successive interference cancellation (SIC) and user velocity on the NOMA performance. Their results showed that even in the worst error propagation scenario, NOMA can outperform the conventional orthogonal multiple access (OMA) for various user mobilities. Chen et al. [4], studied NOMA over the downlink of a wireless system. Traditional minimum-mean-squared-error (MMSE) precoding matrices have been used. Considering that wireless relaying is an effective means to combat the effects of channel fading and shadowing on the transmitted signal, it is reasonable to integrate NOMA into relaying networks. In [5], a new cooperative NOMA scheme was proposed and analyzed in terms of outage probability and diversity gain. The authors in [6] studied the impact of relay selection on the performance of cooperative NOMA. In particular, a two-stage relay selection strategy is proposed, which can achieve the minimal outage probability among all the possible selection schemes, and achieve the full diversity gain. Yu et al. [7] propose a cooperative truncated automatic repeat-request (T-ARQ)-based NOMA scheme with relay selection. The numerical results show that the proposed scheme can achieve much higher throughput than the T-ARQ-based time-division multiple access schemes.

Although cooperative NOMA is capable of improving the performance for far-users by using a half-duplex relay, it results in additional bandwidth cost for the system. To avoid this issue, full duplex (FD) is a promising technique for the next generation wireless systems, because it can double the spectrum efficiency by realizing transmission and reception on the same carrier frequency, simultaneously [8]. Two main types of FD relay techniques, namely FD amplify-and-forward (AF) relaying and FD decode-and-forward (DF) relaying, have been discussed in [9]-[11]. The outage probability of FD-AF relaying were analyzed in [9], which considers the processing delay of relaying in practical scenarios. Osorio et al. [10], analyzed the outage probability of FD-AF relaying with direct link available. The authors in [11] characterized the outage performance of FD-DF relaying. It demonstrated that the optimal duplex mode can be selected according to the outage probability.

In addition to improving the spectral efficiency which is the key motivation of NOMA, another critical goal of 5G networks is to maximize the energy efficiency. Energy harvesting (EH), has recently received considerable attention as a sustainable solution to alleviate the energy shortage of wireless networks. Radio frequency (RF) signals radiated by ambient transmitters can be treated as a viable new source for EH. A fundamental limitation of EH-based wireless communications lies in the restricted transmission range. Recently, relays with EH capabilities have attracted lots of research interests, where the relay can use the energy harvested from the source signals to perform the data forwarding [12]. The relay can be wirelessly powered and the communication coverage can be expanded. In [13], a wireless cooperative network was considered, in which multiple source-destination pairs communicate with each other via an EH relay. In [14], a best cooperative mechanism was proposed  to facilitate the wireless EH and spectrum sharing, which can greatly improve both the energy and spectral efficiencies. Nevertheless, the amount of energy harvested from RF radiation is often restricted and it is desirable for relays to accumulate the harvested energy in the energy storage such as super-capacitors or rechargeable batteries [15]. In [16], a threshold-based energy storage scheme was employed at relays, the stored energy in  the battery of each relay actually forms a Markov chain over time. By investigating the properties of this Markov chain, the asymptotic average throughput was derived. In [17], to support an efficient utilization of harvested energy to improve the system throughput, a harvest-use-store relaying strategy with distributed beamforming were proposed. Zhai et al. [18], proposed a high-efficient EH-based cooperative transmission protocol using either power-splitting or time-switching technique, where both the source and the relay are capable of harvesting energy from a power beacon and then using the harvested energy for the cooperative DF or AF transmission.

In order to improve the spectral efficiency and energy efficiency, we consider the FD-NOMA relaying for the downlink transmission over Nakagami-$m$ fading channels and provide a comprehensive performance analysis. Our analysis provides more general results, as the Nakagami-$m$ fading can incorporate the most commonly used Rayleigh fading and Rice fading as special cases, which can be obtained by adjusting the fading parameter. The main contributions of this paper are summarized as follows:
\begin{itemize}
\item We consider a NOMA-based downlink FD-AF relaying network over Nakagami-$m$ fading channels. The relays have no other energy supplies, but they are equipped with chargeable battery and can harvest and store the wireless energy broadcasted by the source. Based on this structure, a new relay selection scheme in Nakagami-$m$ fading
channels is proposed considering both the channel state information (CSI) and the battery statuses of relays.
\item This model combines
NOMA, cooperative communication, FD and EH together reasonably and can improve
spectral efficiency and energy efficiency simultaneously.
\item We model the amount of energy stored in the battery of each relay as finite states and use a finite Markov chain to model the variation of energy at each relay. Both the transition probability and the steady-state probability are derived for the next performance analysis.
\item The asymptotic expressions of outage probability and the ergodic rate for both near and far users are derived. Finally, simulation results are presented to validate our theoretical analysis and verify the superiority of FD-NOMA relaying over the traditional half-duplex relaying based OMA.
\end{itemize}

The rest of this paper is organized as follows. In section 2, we introduce the system model. The relay selection scheme is illustrated in section 3. The closed-form expressions of outage probability and ergodic rate are derived in section 4. Numerical and simulations results are presented in section 5 and conclusions are drawn in section 6.

\section{Method}
Next, we describe the setting and design of the system and the relay selection rule we provide.
\subsection{System model}

In our system, a source $S$ intends to communicate to a number of potential users ${D_j}$ $(j = 1,2,...,N)$ over channels with flat fading. Multiple potential FD relays ${R_i}$ $\left( {i = 1,2,...M} \right)$ are willing to amplify and forward the signal from $S$ to the users. We assume that there exist direct links between $S$ and users. $S$ and each ${D_j}$ are configured with one antenna and operate in half-duplex (HD) mode, while relays are equipped with one transmit antenna and one receive antenna that enables a FD operation. The channel between any transmitter and any receiver is assumed to follow Nakagami-$m$ fading.  The channels pertaining to the direct link, first-hop and second-hop undergo independent identically (i. i. d.) fading and the channel coefficients are denoted as ${h_{S{D_j}}}$, ${h_{S{R_i}}}$ and ${h_{{R_i}{D_j}}}$ with parameters ${m_0}$, ${m_1}$ and ${m_2}$, respectively. The channel power gain ${G_{S{D_j}}} = {\left| {{{h}_{S{D_j}}}} \right|^2}$, ${G_{S{R_i}}} = {\left| {{{h}_{S{R_i}}}} \right|^2}$ and ${G_{{R_i}{D_j}}} = {\left| {{{h}_{{R_i}{D_j}}}} \right|^2}$ thus follow the Gamma distribution, with mean ${\Phi _{SD}}$, ${\Phi _{SR}}$ and ${\Phi _{RD}}$, respectively.

 In a general case, due to the imperfect isolation or cancellation process, FD operation may suffer from residual self-interference (SI). In order to simplify the analysis, we assume each relay with spatially separated transmit and receive antennas is fixed, and the interference channel is non-fading, so the instantaneous interference-to-noise-ratio (INR) equals the corresponding average INR. As shown later in the simulations, this assumption does not degrade the accuracy of our analysis. Without loss of generality, we assume that users' channel gains are ordered as ${\hat G_{S{D_1}}} \le {\hat G_{S{D_2}}} \le  \cdot  \cdot  \cdot  \le {\hat G_{S{D_N}}}$ (${\hat G_{S{D_j}}}$ is the sorted random variable of ${G_{S{D_j}}}$). In what follows, we focus on the analysis of the group users consisting of ${D_p}$ and ${D_q}$ (${G_{S{D_P}}} < {G_{S{D_q}}}$). The additive white Gaussian noise (AWGN) of all links each link is assumed to have zero mean and equal variance ${\sigma ^2}$.

\begin{figure} [htbp]
\centering \scalebox{0.6}{\includegraphics{1.eps}}
\caption{ A reference model for multi-relay cooperative FD-NOMA network.\label{Fig4}}
\end{figure}

In this paper, we use the conventional AF protocol for the NOMA downlink cooperative transmission and consider a normalized unit block time (i.e., $T=1$). Before the transmission, each relay checks their battery at the beginning of a transmission block and judges whether it has enough energy to forward the source information. If the relay does not have enough energy or is not selected to forward the information, it will perform EH in this time block and store the harvested energy into its battery. We assume ${P_S}$ is the transmission power of source and is sufficiently large such that the energy harvested from the noise is negligible. Thus, for a relay node $R_i$, the amount of energy harvested from the source can be expressed as
\begin{equation}\label{eqn01}
\begin{split}
{E_{{R_i}}} = \eta {P_S}{G_{S{R_i}}}T,
\end{split}
\end{equation}
where $\eta$ is the energy conversion efficiency at each relay.

For those relays with sufficient energy, they report their CSI to $S$ for the relay selection. Let ${E_i}$ denote the amount of energy in the battery of relay $R{}_i$, we define those relays with enough energy as the eligible set
\begin{equation}\label{eqn01}
\begin{split}
\phi  = \left\{ {\left. {{R_i}} \right|{E_i} \ge VC,i = 1,2,...M} \right\},
\end{split}
\end{equation}
where $C$ is the EH threshold to activate the EH circuit and $V$ is a positive integer. If we assume the selected relay for forwarding data has the same transmitting power as $S$, $V$ must satisfy
\begin{equation}\label{eqn01}
\begin{split}
\left( {V - 1} \right)C < {P_S}T \le VC.
\end{split}
\end{equation}
Among this set of relays, ${R_{\bar m}}$ is selected and the index of ${R_{\bar m}}$ can be expressed as
\begin{equation}\label{eqn01}
\begin{split}
\bar m = \arg \mathop {\max }\limits_{i:{R_i} \in \phi } \left\{ {{G_{S{R_i}}}} \right\}.
\end{split}
\end{equation}
In (4), the relay with the best channel power gain of the first hop in this set is selected for forwarding information, all other relays perform EH in this block time.  We assume the transmitting power of the selected relay (${P_{{R_{\bar m}}}}$) is the same as $S$, which means ${P_S} = {P_{{R_{\bar m}}}} = P$. Although the energy consumed of ${R_{\bar m}}$ for sending date in one block time is $PT$, which can be seen from (3), we still consider it as $VC$, because $C$ is the basic unit of the relay battery. The difference between $PT$ and $VC$ does not severely degrade the accuracy of the performance analysis.

After selection, $S$ broadcasts the signal ${x_S}[k] = \sqrt {{a_p}} {x_p}[k] + \sqrt {{a_q}} {x_q}[k]$ (${a_p}$ and ${a_q}$ are the power allocation factors and satisfy ${a_p} + {a_q} = 1$, ${a_p} > {a_q}$), while the selected relay ${R_{\bar m}}$ receives ${x_S}[k]$ and forwards the signal generated by
co-channel transmission to ${D_p}$ and ${D_q}$ in FD mode, while the corresponding ${D_p}$ and ${D_q}$ simultaneously receive the signals transmitted from $S$ and forwarded by ${R_{\bar m}}$. In the $k$th time slot, $k$=1,2,3..., the signals received at ${R_{\bar m}}$, ${D_p}$ and ${D_q}$ can be, respectively, expressed as
\begin{equation}\label{eqn01}
\begin{split}
{y_{{R_{\bar m}}}}[k] = {h_{S{R_{\bar m}}}}\sqrt P {x_S}[k] + {h_{R{R_{\bar m}}}}\sqrt P A{y_{{R_{\bar m}}}}[k - {T_d}] + {n_{{R_{\bar m}}}}[k],
\end{split}
\end{equation}
\begin{equation}\label{eqn01}
\begin{split}
{y_{{D_p}}}[k] = {h_{S{D_p}}}\sqrt P {x_S}[k] + {h_{{R_{\bar m}}{D_p}}}\sqrt P A{y_{{R_{\bar m}}}}[k - {T_d}] + {n_{{D_p}}}[k],
\end{split}
\end{equation}
\begin{equation}\label{eqn01}
\begin{split}
{y_{{D_q}}}[k] = {h_{S{D_q}}}\sqrt P {x_S}[k] + {h_{{R_{\bar m}}{D_q}}}\sqrt P A{y_{{R_{\bar m}}}}[k - {T_d}] + {n_{{D_q}}}[k],
\end{split}
\end{equation}
where $A$ is the amplification coefficient, ${h_{R{R_{\bar m}}}}$ is the residual SI between transmit antenna
and receive antenna of relay, and ${T_d}$ is the processing delay. ${n_{{R_{\bar m}}}}[k]$, ${n_{{D_p}}}[k]$ and ${n_{{D_q}}}[k]$ are independent circularly-symmetric Gaussian random variables. Without loss of generality and for the overall simplicity, we assume $\gamma  = {P \mathord{\left/
 {\vphantom {P {{\sigma ^2}}}} \right.
 \kern-\nulldelimiterspace} {{\sigma ^2}}}$ is the average signal-to-noise-ratio (SNR). We also set ${\gamma _{\bar m}} = \gamma {\rm{E}}\left( {{{\left| {{{h}_{R{R_{\bar m}}}}} \right|}^2}} \right)$ as an actual measured value, which averages out the fading effect. With unity energy transmit signals and after some simple substitutions, we can obtain [19]
\begin{equation}\label{eqn01}
\begin{split}
A= {\left( {{G_{S{R_{\bar m}}}}P + {{\left| {{{h}_{R{R_{\bar m}}}}} \right|}^2}P + {\sigma ^2}} \right)^{ - \frac{1}{2}}}.
\end{split}
\end{equation}

In this study, we assume that the two signals from $S$ and $R_{\bar m}$ are resolvable absolutely (herein the relaying link corresponding to direct link has small time delay for any transmitted signals, the same as in [20]), therefore the instantaneous signal-to-interference-plus-noise-ratios (SINRs) at ${D_p}$ from $S$ and ${R_{\bar m}}$ are, respectively, expressed by
\begin{equation}\label{eqn01}
\begin{split}
{\gamma _{S{D_p}}} = \frac{{{a_p}{{\hat G}_{S{D_p}}}\gamma }}{{{a_q}{{\hat G}_{S{D_p}}}\gamma  + 1}},
\end{split}
\end{equation}
\begin{equation}\label{eqn01}
\begin{split}
{\gamma _{{R_{\bar m}}{D_p}}} = \frac{{{a_p}{\gamma ^2}{G_{S{R_{\bar m}}}}{G_{{R_{\bar m}}{D_p}}}}}{{\gamma \left( {{G_{S{R_{\bar m}}}} + {G_{{R_{\bar m}}{D_p}}}} \right) + \gamma {G_{{R_{\bar m}}{D_p}}}\left( {{a_q}\gamma {G_{S{R_{\bar m}}}} + {\gamma _{\bar m}}} \right) + {\gamma _{\bar m}} + 1}}.
\end{split}
\end{equation}
Using the selection combining (SC), we can obtain the instantaneous SINR of ${D_p}$ as
\begin{equation}\label{eqn01}
\begin{split}
{\gamma _{{D_p}}} = \max \left( {{\gamma _{S{D_p}}},{\gamma _{{R_{\bar m}}{D_p}}}} \right).
\end{split}
\end{equation}
As for ${D_q}$, the SIC technique is applied to decode information of ${D_p}$. Particularly, ${D_q}$ first decodes the
message of ${D_p}$, then subtracts this component from the received signal to detect its own information. Therefore, the received SINRs at ${D_q}$ for the decoding of the information of ${D_p}$ are given by
\begin{equation}\label{eqn01}
\begin{split}
{\gamma _{S{D_{q \to p}}}} = \frac{{{a_p}{{\hat G}_{S{D_q}}}\gamma }}{{{a_q}{{\hat G}_{S{D_q}}}\gamma  + 1}},
\end{split}
\end{equation}
\begin{equation}\label{eqn01}
\begin{split}
{\gamma _{{R_{\bar m}}{D_{q \to p}}}} = \frac{{{a_p}{\gamma ^2}{G_{S{R_{\bar m}}}}{G_{{R_{\bar m}}{D_q}}}}}{{\gamma \left( {{G_{S{R_{\bar m}}}} + {G_{{R_{\bar m}}{D_q}}}} \right) + \gamma {G_{{R_{\bar m}}{D_q}}}\left( {{a_q}\gamma {G_{S{R_{\bar m}}}} + {\gamma _{\bar m}}} \right) + {\gamma _{\bar m}} + 1}}.
\end{split}
\end{equation}
After the SIC decoding, the instantaneous SINRs at ${D_q}$ for the decoding of the information of  its own information are expressed by
\begin{equation}\label{eqn01}
\begin{split}
{\gamma _{S{D_q}}} = {a_q}{G_{S{D_q}}}\gamma,
\end{split}
\end{equation}
\begin{equation}\label{eqn01}
\begin{split}
{\gamma _{{R_{\bar m}}{D_q}}} = \frac{{{a_q}{\gamma ^2}{G_{S{R_{\bar m}}}}{G_{{R_{\bar m}}{D_q}}}}}{{\gamma \left( {{G_{S{R_{\bar m}}}} + {G_{{R_{\bar m}}{D_q}}}} \right) + \gamma {G_{{R_{\bar m}}{D_q}}}{\gamma _{\bar m}} + {\gamma _{\bar m}} + 1}}.
\end{split}
\end{equation}
Similarly to (10), according to the SC, the instantaneous SINR of ${D_q}$ can be expressed as
\begin{equation}\label{eqn01}
\begin{split}
{\gamma _{{D_q}}} = \max \left( {{\gamma _{S{D_q}}},{\gamma _{{R_{\bar m}}{D_q}}}} \right).
\end{split}
\end{equation}

\subsection{Relay selection rule}

In the proposed scheme, each relay checks their battery at the beginning of a transmission block and judges whether it has enough energy to forward the source information. For those relays with sufficient energy, they report their CSI to $S$ for making the relay selection decision. All the other relays apart from ${R_{\bar m}}$ perform the EH in this time block. We assume each relay accumulates the harvested energy using a finite energy storage with capacity of $BC$ ($B$=1, 2, 3...). We also define $C$ as the EH threshold and $PT$ as the energy required for the date forwarding. The relay can harvest $uC$ amount of energy in one time slot and $u$ is given by
\begin{equation}\label{eqn01}
\begin{split}
uC \le {E_{{R_i}}} < \left( {u + 1} \right)C. \left({\rm{ }}u = 0,{\rm{ }}1,{\rm{ }}2.....\right)
\end{split}
\end{equation}
This assumption is closer to the practical scenario, and the evolution of the battery status of each relay can be modeled as a finite-state Markov chain. Using the transition probability matrix of this chain, we can get the steady-state probability vector which can be used for analyzing the performance.

Due to the lack of general form of steady-state probability, the above analysis is computationally intense when $M$ is large. To facilitate the computation, we propose an approximate approach based on two simplified assumptions. Firstly, we denote the energy amount of relay at the selection epoch as a random variable $Z$. The exact distribution of $Z$ is very complex. To ease the computation, we approximate $Z$ by a uniform random variable over [0, $BC$]. The adopted approximation is inspired by the fact that the amount of harvested energy in a transmission block follows the geometric distribution with parameter 1/2 [21]. In general, these conditions may not always hold. We will discuss the effectiveness of the assumption through numerical results. Secondly, we found that an arbitrary relay may be either short of enough power to participate in the relay selection or otherwise, so the evolution of energy amount of relay is captured by using two states, either active or inactive. With this simplified two-state Markov chain, a relay is in ${s_0}$ if the relay lacks sufficient energy to transmit, or in ${s_1}$ when the relay has enough energy for transmission. Next, we will explain how to obtain the transition probability matrix of the two-state Markov chain.

The transition from state ${s_0}$ to ${s_0}$ happens when a relay has no enough energy to transmit (i.e., $Z<VC$) in the current block and the accumulated energy after the EH in a time block remains below the given threshold. The corresponding transition probability is given by
\begin{equation}\label{eqn01}
\begin{split}
{p_{0,0}} = \Pr \left( {\left. {Z + {E_{{R_i}}} < VC} \right|0 \le Z < VC} \right){\rm{ = Pr}}\left( {\hat Z + {E_{{R_i}}} < VC} \right),
\end{split}
\end{equation}
where $\hat Z$ is a truncated random variable defined as
\begin{equation}\label{eqn01}
\begin{split}
\hat Z = \left\{ {\begin{array}{*{20}{c}}
{Z \qquad Z < VC}\\
{0\qquad Z \ge VC}
\end{array}}. \right.
\end{split}
\end{equation}
Since $Z$ is approximated as uniformly distributed, the probability density function (PDF) of $\hat Z$ can be obtained as
\begin{equation}\label{eqn01}
\begin{split}
{f_{\hat Z}}\left( z \right) = \frac{1}{{BC}}u(VC - z) + \left( {1 - \frac{V}{B}} \right)\delta \left( {z - VC} \right),
\end{split}
\end{equation}
where $u\left(  \cdot  \right)$ and $\delta \left(  \cdot  \right)$ denote the unit step function and the Dirac delta function, respectively. Considering the channels obey Nakagami-$m$ distribution, a general random variable ${h_\Lambda }$, $\Lambda  \in \left\{ {S{D_j},S{R_i},{R_i}{D_j}} \right\}$, is subject to a Gamma distribution, so the corresponding PDF is
\begin{equation}\label{eqn01}
\begin{split}
{f_{{G_\Lambda }}}(x) = {\left( {\frac{{{m_\Lambda }}}{{{\Phi _\Lambda }}}} \right)^{{m_\Lambda }}}\frac{{{x^{{m_\Lambda } - 1}}}}{{\Gamma ({m_\Lambda })}}\exp \left( { - \frac{{{m_\Lambda }x}}{{{\Phi _\Lambda }}}} \right),
\end{split}
\end{equation}
where $\Gamma \left(  \cdot  \right)$ is Gamma function [22]. The corresponding cumulative distribution function (CDF) is expressed by
\begin{equation}\label{eqn01}
\begin{split}
{F_{{G_\Lambda }}}\left( x \right)\; = 1 - \exp \left( { - \frac{{{m_\Lambda }x}}{{{\Phi _\Lambda }}}} \right)\sum\limits_{n = 0}^{{m_\Lambda } - 1} {\frac{{{{\left( {{{{m_\Lambda }x} \mathord{\left/
 {\vphantom {{{m_\Lambda }x} {{\Phi _\Lambda }}}} \right.
 \kern-\nulldelimiterspace} {{\Phi _\Lambda }}}} \right)}^n}}}{{n!}}}.
\end{split}
\end{equation}
If ${m_\Lambda }$ takes positive integer values, then with the help of [22, (3.351)], (18) can be calculated as
\begin{equation}\label{eqn01}
\begin{split}
{p_{0,0}} &= \int_0^{VC} {\Pr \left( {{G_{S{R_i}}} < \frac{{VC - z}}{{\eta PT}}} \right)} {f_{\hat Z}}\left( z \right)dz\\
 &= 1 - \frac{1}{{BC}}\sum\limits_{n = 0}^{{m_1} - 1} {\left( {\frac{1}{{{{{m_1}} \mathord{\left/
 {\vphantom {{{m_1}} {{{\bar E}_{{R_i}}}}}} \right.
 \kern-\nulldelimiterspace} {{{\bar E}_{{R_i}}}}}}} - \exp \left( { - \frac{{{m_1}VC}}{{{{\bar E}_{{R_i}}}}}} \right)\sum\limits_{k = 0}^n {\frac{{n!{{\left( {VC} \right)}^k}}}{{k!{{\left( {{{{m_1}} \mathord{\left/
 {\vphantom {{{m_1}} {{{\bar E}_{{R_i}}}}}} \right.
 \kern-\nulldelimiterspace} {{{\bar E}_{{R_i}}}}}} \right)}^{1 - k}}}}} } \right)},
\end{split}
\end{equation}
where ${\bar E_{{R_i}}} = {\rm{E}}\left( {{E_{{R_i}}}} \right) = P\eta {\Phi _{SR}}$.

The transition from ${s_0}$ to ${s_1}$ happens when the relay enters the EH mode in the current bock and the accumulated energy exceeds $VC$. Hence, we have
\begin{equation}\label{eqn01}
\begin{split}
{p_{0,1}} = \Pr \left( {\left. {Z + {E_{{R_i}}} \ge VC} \right|0 \le Z < VC} \right){\rm{ = Pr}}\left( {\hat Z + {E_{{R_i}}} \ge VC} \right).
\end{split}
\end{equation}
Similar to the derivation of (22), we can obtain
\begin{equation}\label{eqn01}
\begin{split}
{p_{0,1}} = \frac{1}{{BC}}\sum\limits_{n = 0}^{{m_1} - 1} {\left( {\frac{1}{{{{{m_1}} \mathord{\left/
 {\vphantom {{{m_1}} {{{\bar E}_{{R_i}}}}}} \right.
 \kern-\nulldelimiterspace} {{{\bar E}_{{R_i}}}}}}} - \exp \left( { - \frac{{{m_1}VC}}{{{{\bar E}_{{R_i}}}}}} \right)\sum\limits_{k = 0}^n {\frac{{n!{{\left( {VC} \right)}^k}}}{{k!{{\left( {{{{m_1}} \mathord{\left/
 {\vphantom {{{m_1}} {{{\bar E}_{{R_i}}}}}} \right.
 \kern-\nulldelimiterspace} {{{\bar E}_{{R_i}}}}}} \right)}^{1 - k}}}}} } \right)}.
\end{split}
\end{equation}

If the relay which is in state ${s_1}$ is selected as the best relay for the data forwarding and its harvested energy transits from state ${s_1}$ to state ${s_0}$. In this case, the transition probability ${p_{1,0}}$ is then given by
\begin{equation}\label{eqn01}
\begin{split}
{p_{1,0}}&= \Pr \left( {\left. {Z - VC < VC} \right|VC \le Z \le BC} \right){\rm{Pr}}\left( {{R_{\bar m}} = {R_i}} \right)\\
&= \Pr\left( {Z' < 2VC} \right){\rm{Pr}}\left( {{R_{\bar m}} = {R_i}} \right),
\end{split}
\end{equation}
where $Z'$ is a truncated random variable defined as
\begin{equation}\label{eqn01}
%\begin{split}
Z' = \left\{ {\begin{array}{rcl}
 Z&  &VC \le Z \le BC\\
 0&   &Z < VC
\end{array}}. \right.
%\end{split}
\end{equation}
Since $Z$ is uniformly distributed, the PDF of $Z'$ can be obtained as
\begin{equation}\label{eqn01}
\begin{split}
{f_{Z'}}\left( z \right) = \frac{1}{{BC}}\left( {u(z - BC) - u(z - VC)} \right) + \frac{V}{B}\delta \left( {z - VC} \right).
\end{split}
\end{equation}
Using this PDF, the first term on the right-hand side of (26) can be obtained as
\begin{equation}\label{eqn02}
%\begin{split}
{\rm{Pr}}\left( {Z' < 2VC} \right) = \left\{ {\begin{array}{rcl}
&{{2V} \mathord{\left/
 {\vphantom {{2V} B}} \right.
 \kern-\nulldelimiterspace} B}\qquad &2VC < BC\\
&1 \qquad  &2VC \ge BC
\end{array}}. \right.
%\end{split}
\end{equation}
As to the second term, we note that the i. i. d. fading assumption implies each relay in $\phi$ has an equal chance to be selected as the best relay. To simplify the analysis, we approximate the cardinality of $\phi$ by its mean such that
\begin{equation}\label{eqn01}
\begin{split}
{\rm{Pr}}\left( {{R_{\bar m}} = {R_i}} \right) \approx \frac{1}{{{{\bar M}_e}}} = \frac{1}{{M{v_1}}},
\end{split}
\end{equation}
where ${\bar M_e}$ is the average number of relays in $\phi$ and ${v_1}$ is the steady-state probability of state $s_1$. Combining (29) and (30), we can obtain the closed-form of ${p_{1,0}}$.

The transition probability from $s_1$ to $s_1$ can be derived in the similar manner as the previous case, so we omit the derivation here.

When the two-state Markov chain is formulated, the steady-state probability vector can be directly obtained as
\begin{equation}\label{eqn01}
\begin{split}
v = ({v_0},{v_1}) = \left( {\frac{{{p_{1,0}}}}{{{p_{0,1}} + {p_{1,0}}}},\frac{{{p_{0,1}}}}{{{p_{0,1}} + {p_{1,0}}}}} \right).
\end{split}
\end{equation}
We note that both $v_0$ and $v_1$ involve ${p_{1,0}}$, which is a function of $v_1$. By substituting ${p_{0,1}}$ and ${p_{1,0}}$ into (31), $v_1$ can be solved explicitly. Take the condition $2VC<BC$ in (29), for example, $v_1$ can be derived in closed form as
\begin{equation}\label{eqn01}
\begin{split}
{v_1} = 1 - \frac{{2V}}{{BM}}{\left( {\frac{1}{{BC}}\sum\limits_{n = 0}^{{m_i} - 1} {\left( {\frac{1}{{{{{m_1}} \mathord{\left/
 {\vphantom {{{m_1}} {{{\bar E}_{{R_i}}}}}} \right.
 \kern-\nulldelimiterspace} {{{\bar E}_{{R_i}}}}}}} - \exp \left( { - \frac{{{m_1}VC}}{{{{\bar E}_{{R_i}}}}}} \right)\sum\limits_{k = 0}^n {\frac{{n!{{\left( {VC} \right)}^k}}}{{k!{{\left( {{{{m_1}} \mathord{\left/
 {\vphantom {{{m_1}} {{{\bar E}_{{R_i}}}}}} \right.
 \kern-\nulldelimiterspace} {{{\bar E}_{{R_i}}}}}} \right)}^{1 - k}}}}} } \right)} } \right)^{ - 1}}.
\end{split}
\end{equation}
With $v_1$ at hand, we can get $\Pr \left( {\left| \phi  \right| = \Omega } \right)$ ($\Omega  = 0,{\rm{ }}1,{\rm{ }}2...$), which follows the binomial distribution with the probability mass function given as
\begin{equation}\label{eqn01}
\begin{split}
\Pr \left( {\left| \phi  \right| = \Omega } \right) = {{M}\choose {\Omega}}{\left( {{v_1}} \right)^\Omega }{\left( {1 - {v_1}} \right)^{M - \Omega }},
\end{split}
\end{equation}
where $ {a \choose b}= \frac{{a!}}{{(a - b)!b!}}$. With the close form of $\Pr \left( {\left| \phi  \right| = \Omega } \right)$, we can analyze the performance in the following part.

\section{Performance analysis}
\subsection{Outage behavior}

According to the total probability law, the outage probability of $D_p$ or $D_q$ can be expressed as
\begin{equation}\label{eqn01}
\begin{split}
P_{out}^\Gamma  = \sum\nolimits_{\Omega  = 1}^M {\Pr \left( {\left| \phi  \right| = \Omega } \right)} P_{out}^{\Gamma ,\Omega },
\end{split}
\end{equation}
where $\Gamma  = \left\{ {p,q} \right\}$. As $\Pr \left( {\left| \phi  \right| = \Omega } \right)$ can be got in (33), the next step is to get $P_{out}^{\Gamma ,\Omega }$, which will be given in the following.

The scheme would be in outage if the SINR falls below the pre-set threshold, so the outage probability of $D_p$ is given by
\begin{equation}\label{eqn01}
\begin{split}
P_{out}^{p,\Omega } &= \Pr \left( {{\gamma _{{D_p}}} < {\gamma _{thp}}} \right) = \Pr \left( {\max \left( {{\gamma _{S{D_p}}},{\gamma _{{R_{\bar m}}{D_p}}}} \right) < {\gamma _{thp}}} \right)\\
 &= \underbrace {\Pr \left( {{\gamma _{S{D_p}}} < {\gamma _{thp}}} \right)}_{{J_1}}\underbrace {\Pr \left( {{\gamma _{{R_{\bar m}}{D_p}}} < {\gamma _{thp}}} \right)}_{{J_2}},
\end{split}
\end{equation}
where ${\gamma _{thp}}$ is the threshold of $D_p$. Next, according to order statistics [23], the PDF and CDF of the sorted random variable ${\hat h_{S{D_p}}}$ are expressed, respectively, by
\begin{equation}\label{eqn01}
\begin{split}
{f_{{{\hat G}_{S{D_p}}}}}(x) = \sum\nolimits_{j = 0}^{N - p} {{{\left( { - 1} \right)}^j}} Z_p^N {{N-p}\choose {j}}{f_{{G_{S{D_j}}}}}(x){\left( {{F_{{G_{S{D_j}}}}}\left( x \right)\;} \right)^{p + j - 1}},
\end{split}
\end{equation}
\begin{equation}\label{eqn01}
\begin{split}
{F_{{{\hat G}_{S{D_p}}}}}\left( x \right)\; = \sum\nolimits_{j = 0}^{N - p} {\frac{{{{\left( { - 1} \right)}^j}Z_p^N}}{{p + j}}}{{N-p}\choose {j}}{\left( {{F_{{G_{S{D_j}}}}}\left( x \right)\;} \right)^{p + j}},
\end{split}
\end{equation}
where $ Z_p^N = {{N!} \mathord{\left/
 {\vphantom {{N!} {\left( {\left( {N - p} \right)!\left( {p - 1} \right)!} \right)}}} \right.
 \kern-\nulldelimiterspace} {\left( {\left( {N - p} \right)!\left( {p - 1} \right)!} \right)}}$. Based on the above definitions, $J_1$ is calculated by
\begin{equation}\label{eqn01}
\begin{split}
{J_1}& = \Pr \left( {\frac{{{a_p}{{\hat G}_{S{D_p}}}\gamma }}{{{a_q}{{\hat G}_{S{D_p}}}\gamma  + 1}} < {\gamma _{thp}}} \right) = \int_0^\tau  {{f_{{{\hat G}_{S{D_p}}}}}(x)dx} \\
 &= \sum\nolimits_{j = 0}^{N - p} {Z_p^N\frac{{{{\left( { - 1} \right)}^j}}}{{p + j}}} {{N-p}\choose j}{\left( {1 - \exp \left( { - \frac{{{m_0}\tau }}{{{\Phi _{SD}}}}} \right)\sum\limits_{n = 0}^{{m_0} - 1} {\frac{{{{\left( {{{{m_0}\tau } \mathord{\left/
 {\vphantom {{{m_0}\tau } {{\Phi _{SD}}}}} \right.
 \kern-\nulldelimiterspace} {{\Phi _{SD}}}}} \right)}^n}}}{{n!}}} \;} \right)^{p + j}},
\end{split}
\end{equation}
where $\tau  = {{{\gamma _{thp}}} \mathord{\left/
 {\vphantom {{{\gamma _{thp}}} {\gamma \left( {{a_p} - {a_q}{\gamma _{thp}}} \right)}}} \right.
 \kern-\nulldelimiterspace} {\gamma \left( {{a_p} - {a_q}{\gamma _{thp}}} \right)}}$. It is noteworthy that (38) exists if and only if ${\gamma _{thp}} < {{{a_p}} \mathord{\left/
 {\vphantom {{{a_p}} {{a_q}}}} \right.
 \kern-\nulldelimiterspace} {{a_q}}}$, otherwise the probability will be always be one. $J_2$ can be simplified as
\begin{equation}\label{eqn01}
\begin{split}
{J_2} =& \Pr \left( {{G_{{R_{\bar m}}{D_p}}} \le \tau } \right)\\& + \Pr \left( {{G_{{R_{\bar m}}{D_p}}} > \tau ,{h_{S{R_{\bar m}}}} \le \frac{{\tau \left( {\gamma {G_{{R_{\bar m}}{D_p}}}{\gamma _{\bar m}} + \gamma {G_{{R_{\bar m}}{D_p}}} + {\gamma _{\bar m}} + 1} \right)}}{{\gamma \left( {{G_{{R_{\bar m}}{D_p}}} - \tau } \right)}}} \right)\\
 =& \int_0^\tau  {{f_{{G_{{R_{\bar m}}{D_p}}}}}} \left( x \right)dx + \int_\tau ^\infty  {{f_{{G_{{R_{\bar m}}{D_p}}}}}} \left( y \right)dy\int_0^{\frac{{\tau \left( {\gamma y{\gamma _{\bar m}} + \gamma y + {\gamma _{\bar m}} + 1} \right)}}{{\gamma \left( {y - \tau } \right)}}} {{f_{{G_{S{R_{\bar m}}}}}}\left( x \right)} dx.
\end{split}
\end{equation}
To solve (39), we first need to calculate ${f_{{G_{S{R_{\bar m}}}}}}\left( x \right)$, which can be get according to (36) when $p=1$ and $N = \Omega $. Using ${f_{{G_{S{R_{\bar m}}}}}}\left( x \right)$ and with the help of [22, (3.471.9)], $J_2$ can be formulated as
\begin{equation}\label{eqn01}
\begin{split}
{J_2} =& 1 - \exp \left( { - \frac{{{m_2}\tau }}{{{\Phi _{RD}}}}} \right)\sum\limits_{n = 0}^{{m_2} - 1} {\frac{{{{\left( {{{{m_2}\tau } \mathord{\left/
 {\vphantom {{{m_2}\tau } {{\Phi _{RD}}}}} \right.
 \kern-\nulldelimiterspace} {{\Phi _{RD}}}}} \right)}^n}}}{{n!}}} \\
 &+ \frac{{Z_1^\Omega }}{{\Gamma ({m_2})}}{\left( {\frac{{{m_2}}}{{{\Phi _{RD}}}}} \right)^{{m_2}}}\exp \left( { - \frac{{{m_2}\tau }}{{{\Phi _{RD}}}}} \right)
 \\&\times\sum\nolimits_{j = 0}^{\Omega  - 1} {\frac{{{{\left( { - 1} \right)}^j}}}{{j + 1}}}
{{\Omega  - 1}\choose j}
\sum\nolimits_{r = 0}^j {{{\left( { - 1} \right)}^r}} \exp \left( { - \frac{{{m_1}r{u_{\bar m}}}}{{{\Phi _{SR}}}}} \right)\\
& \times \sum\limits_{{r_i} \ge 0,{r_0} + {r_1} + ... + {r_{{m_1} - 1}} = r} {\frac{{
{r \choose {{r_0},{r_1},...,{r_{{m_1} - 1}}}}
{{\left( {\frac{{{m_1}{u_{\bar m}}}}{{{\Phi _{SR}}}}} \right)}^{\sum\nolimits_{z = 0}^{{m_1} - 1} {z{r_z}} }}}}{{\prod\nolimits_{z = 0}^{{m_1} - 1} {{{\left( {z!} \right)}^{{r_{_z}}}}} }}} \\
 &\times \sum\nolimits_{{k_1} = 0}^{\sum\nolimits_{z = 0}^{{m_1} - 1} {z{r_z}} } {\sum\nolimits_{{k_2} = 0}^{{m_2} - 1} {
{{\sum\nolimits_{z = 0}^{{m_1} - 1} {z{r_z}} } \choose {k_1}}
{{{m_2} - 1} \choose {k_2}}
{d^{{k_1}}}{\tau ^{{m_2} - 1 - {k_2}}}} } \\
 &\times 2{\left( {\frac{{{m_1}r{u_{\bar m}}d{\Phi _{RD}}}}{{{\Phi _{SR}}{m_2}}}} \right)^{\frac{{{k_2} - {k_1} + 1}}{2}}}{K_{{k_2} - {k_1} + 1}}\left( {\sqrt {\frac{{2{m_1}r{u_{\bar m}}d{m_2}}}{{{\Phi _{SR}}{\Phi _{RD}}}}} } \right),
\end{split}
\end{equation}
where ${u_{\bar m}} = \tau \left( {{\gamma _{_{\bar m}}} + 1} \right)$, $d = \tau  + {1 \mathord{\left/
 {\vphantom {1 \gamma }} \right.
 \kern-\nulldelimiterspace} \gamma }$, $Z_1^\Omega  = {{\Omega !} \mathord{\left/
 {\vphantom {{\Omega !} {\left( {\Omega  - 1} \right)!}}} \right.
 \kern-\nulldelimiterspace} {\left( {\Omega  - 1} \right)!}}$, ${r \choose {{r_0},{r_1},...,{r_{{m_1} - 1}}}}
 = \frac{{r!}}{{{r_0}!,{\rm{ }}{r_1}!,{\rm{ }}...,{\rm{ }}{r_{{m_1} - 1}}!}}$, $\sum\limits_{{r_i} \ge 0,{r_0} + {r_1} + ... + {r_{{m_1} - 1}} = r}$ means the sum of all the possible values of which satisfies ${r_0} + {r_1} + ... + {r_{{m_1} - 1}} = r$ and ${K_v}\left(  \cdot  \right)$ expresses the $v$th-order modified Bessel function of the second kind. By substituting (33), (38) and (40) into (34), the outage probability of $D_p$ can be obtained.

 Next, the analysis of an outage event at $D_q$ is presented. Owing to the mechanism of SIC, the outage probability of $D_q$ is given as
\begin{equation}\label{eqn01}
\begin{split}
P_{out}^{q,\Omega } =& \Pr \left( {{\gamma _{{D_q}}} < {\gamma _{thq}}} \right)\\
 =& \Pr \left( {{\gamma _{S{D_q}}} < {\gamma _{thq}}} \right)\Pr \left( {{\gamma _{R{D_q}}} < {\gamma _{thq}}} \right)\\
 =& \underbrace {\left( {1 - \Pr \left( {{\gamma _{S{D_{q \to p}}}} \ge {\gamma _{thp}}} \right)\Pr \left( {{\gamma _{S{D_q}}} \ge {\gamma _{thq}}} \right)} \right)}_{{J_3}}\\
 &\times \underbrace {\left( {1 - \Pr \left( {{\gamma _{{R_{\bar m}}{D_{q \to p}}}} \ge {\gamma _{thp}}} \right)\Pr \left( {{\gamma _{{R_{\bar m}}{D_q}}} \ge {\gamma _{thq}}} \right)} \right)}_{{J_4}},
\end{split}
\end{equation}
where ${\gamma _{thq}}$ is the threshold of $D_q$. $J_3$ can be expanded as
\begin{equation}\label{eqn01}
\begin{split}
{J_3} = &1 - \Pr \left( {\frac{{{a_p}\gamma {{\hat G}_{S{D_q}}}}}{{{a_q}\gamma {{\hat G}_{S{D_q}}} + 1}} \ge {\gamma _{thp}}} \right)\Pr \left( {\gamma {a_q}{{\hat G}_{S{D_q}}} \ge {\gamma _{thq}}} \right)\\
 = &1 - \Pr \left( {{{\hat G}_{S{D_q}}} \ge \max \left( {\frac{{{\gamma _{thp}}}}{{\gamma \left( {{a_p} - {a_q}{\gamma _{thp}}} \right)}},\frac{{{\gamma _{thq}}}}{{\gamma {a_q}}}} \right) \buildrel \Delta \over = \beta } \right)\\
 = &\sum\nolimits_{j = 0}^{N - q} {Z_p^N\frac{{{{\left( { - 1} \right)}^j}}}{{q + j}}}
{{N - q}\choose j}{\left( {1 - \exp \left( { - \frac{{{m_0}\beta }}{{{\Phi _{SD}}}}} \right)\sum\limits_{n = 0}^{{m_0} - 1} {\frac{{{{\left( {{{{m_0}\beta } \mathord{\left/
 {\vphantom {{{m_0}\beta } {{\Phi _{SD}}}}} \right.
 \kern-\nulldelimiterspace} {{\Phi _{SD}}}}} \right)}^n}}}{{n!}}} \;} \right)^{q + j}}.
\end{split}
\end{equation}
Next, $J_4$ can be written as
\begin{equation}\label{eqn01}
\begin{split}
{J_4} =& 1 - \Pr \left( {\frac{{{a_p}{\gamma ^2}{G_{S{R_{\bar m}}}}{G_{{R_{\bar m}}{D_q}}}}}{{\gamma \left( {{G_{S{R_{\bar m}}}} + {G_{{R_{\bar m}}{D_q}}}} \right) + \gamma {G_{{R_{\bar m}}{D_q}}}\left( {{a_q}\gamma {G_{S{R_{\bar m}}}} + {\gamma _{{R_{\bar m}}}}} \right) + {\gamma _{_{\bar m}}} + 1}} \ge {\gamma _{thp}}} \right)\\
 &\times \Pr \left( {\frac{{{a_q}{\gamma ^2}{G_{S{R_{\bar m}}}}{G_{{R_{\bar m}}{D_q}}}}}{{\gamma \left( {{G_{S{R_{\bar m}}}} + {G_{{R_{\bar m}}{D_q}}}} \right) + \gamma {G_{{R_{\bar m}}{D_q}}}{\gamma _{_{\bar m}}} + {\gamma _{_{\bar m}}} + 1}} \ge {\gamma _{thq}}} \right)\\
 = &1 - \int_\beta ^\infty  {\int_{\frac{{\beta \left( {\gamma y{\gamma _{_{\bar m}}} + \gamma y + {\gamma _{_{\bar m}}} + 1} \right)}}{{\gamma \left( {y - \beta } \right)}}}^\infty  {{f_{{G_{S{R_{\bar m}}}}}}\left( x \right){f_{{G_{{R_{\bar m}}{D_q}}}}}\left( y \right)dxdy} } \\
 = &\int_0^\beta  {{f_{{G_{{R_{\bar m}}{D_q}}}}}\left( y \right)dy + } \int_\beta ^\infty  {{f_{{G_{{R_{\bar m}}{D_q}}}}}\left( y \right)dy\int_0^{\frac{{\beta \left( {\gamma y{\gamma _{_{\bar m}}} + \gamma y + {\gamma _{_{\bar m}}} + 1} \right)}}{{\gamma \left( {y - \beta } \right)}}} {{f_{{G_{S{R_{\bar m}}}}}}\left( x \right)dx} }.
\end{split}
\end{equation}
The derivation of (43) is the same as (39), so we omit it here. So far, we have derived the outage probability of $D_q$ by substituting (33), (42) and (43) into (34).
\subsection{Ergodic sum-rate}
In this section, we analyse the ergodic sum-rate for a pair of users, which is given by
\begin{equation}\label{eqn01}
\begin{split}
{R_{sum}} = R_{ave}^p + R_{ave}^q = {\rm{E}}\left( {{{\log }_2}\left( {1 + {\gamma _{{D_p}}}} \right)} \right) + {\rm{E}}\left( {{{\log }_2}\left( {1 + {\gamma _{{D_q}}}} \right)} \right).
\end{split}
\end{equation}
The rate of $D_p$ can be calculated in the following, where approximation is attained under the condition of $\gamma  \to \infty $
\begin{equation}\label{eqn01}
\begin{split}
 R_{ave}^p& = {\rm{E}}\left( {{{\log }_2}\left( {1 + {\gamma _{{D_p}}}} \right)} \right) \\
  &= {\rm{E}}\left( {{{\log }_2}\left( {1 + \max \left( {{\gamma _{S{D_p}}},{\gamma _{{R_{\bar m}}{D_p}}}} \right)} \right)} \right) \\
  &\approx {\rm{E}}\left( {{{\log }_2}\left( {1 + \max \left( {\frac{{{a_p}}}{{{a_q}}},\frac{{{a_p}\gamma {G_{S{R_{\bar m}}}}}}{{{a_q}\gamma {G_{S{R_{\bar m}}}}}}} \right)} \right)} \right) = {\log _2}\left( {\frac{{{a_p}}}{{{a_q}}}} \right). \\
\end{split}
\end{equation}
The rate of $D_q$ can be calculated as
\begin{equation}\label{eqn01}
\begin{split}
R_{ave}^q =& {\rm{E}}\left( {{{\log }_2}\left( {1 + {\gamma _{{D_q}}}} \right)} \right)\\
 =& \Pr \left( {{\gamma _{S{D_{q \to p}}}} \ge {\gamma _{thp}}} \right)\Pr \left( {{\gamma _{{R_{\bar m}}{D_{q \to p}}}} \le {\gamma _{thp}}} \right)\underbrace {{\rm{E}}\left( {{{\log }_2}\left( {1 + {\gamma _{S{D_q}}}} \right)} \right)}_{{J_5}}\\
 &+ \Pr \left( {{\gamma _{S{D_{q \to p}}}} \le {\gamma _{thp}}} \right)\Pr \left( {{\gamma _{{R_{\bar m}}{D_{q \to p}}}} \ge {\gamma _{thp}}} \right)\underbrace {{\rm{E}}\left( {{{\log }_2}\left( {1 + {\gamma _{{R_{\bar m}}{D_q}}}} \right)} \right)}_{{J_6}}\\
 &+ \Pr \left( {{\gamma _{S{D_{q \to p}}}} \ge {\gamma _{thp}}} \right)\Pr \left( {{\gamma _{{R_{\bar m}}{D_{q \to p}}}} \ge {\gamma _{thp}}} \right)\\&\times\underbrace {{\rm{E}}\left( {{{\log }_2}\left( {1 + \max \left( {{\gamma _{S{D_q}}},{\gamma _{{R_{\bar m}}{D_q}}}} \right)} \right)} \right)}_{{J_7}}.
\end{split}
\end{equation}
All the probabilities in the (46) can be obtained according to the above analysis, we only need derive $J_5$, $J_6$ and $J_7$ in the following. Firstly, $J_5$ can be written as
\begin{equation}\label{eqn01}
\begin{split}
{J_5} &= {\rm{E}}\left( {\frac{1}{{\ln 2}}\ln \left( {1 + {a_q}{{\hat G}_{S{D_q}}}\gamma } \right)} \right) \\&= \frac{1}{{\ln 2}}\int_0^\infty  {\ln \left( {1 + {a_q}\gamma x} \right){f_{{{\hat G}_{S{D_q}}}}}} \left( x \right)dx\\
 &= \frac{1}{{\ln 2}}\int_0^\infty  {\frac{{1 - {F_{{{\hat G}_{S{D_q}}}}}\left( {{x \mathord{\left/
 {\vphantom {x {{a_q}\gamma }}} \right.
 \kern-\nulldelimiterspace} {{a_q}\gamma }}} \right)}}{{1 + x}}} dx.
\end{split}
\end{equation}
Using [22, (3.351.2)] and [22, (3.353.5)], the closed form of $J_5$ can be get
\begin{equation}\label{eqn01}
\begin{split}
 {J_5} =& \frac{{Z_q^N}}{{\ln 2\Gamma ({m_0})}}\sum\nolimits_{j = 0}^{N - q} {{{\left( { - 1} \right)}^j}}
   {{N - q} \choose{j}}\sum\nolimits_{r = 0}^{q + j - 1} {}  \\
  &\times \sum\limits_{{r_i} \ge 0,{r_0} + {r_1} + ... + {r_{{m_2} - 1}} = r} {{{\left( { - 1} \right)}^r}\frac{
   {r\choose {{r_0},{r_1},...,{r_{{m_2} - 1}}} } }{{\prod\nolimits_{z = 0}^{{m_2} - 1} {{{\left( {z!} \right)}^{{r_z}}}} }}\left( {\sum\nolimits_{z = 0}^{{m_0} - 1} {z{r_z}}  + {m_0} - 1} \right)!}  \\
  &\times \sum\nolimits_{t = 1}^{\sum\nolimits_{z = 0}^{{m_0} - 1} {z{r_z}}  + {m_0} - 1} {\frac{{{\Phi _{SD}}^{ - t}}}{{t!{{\left( {{a_q}\gamma } \right)}^t}{m_0}^{ - t}{{\left( {r + 1} \right)}^{\sum\nolimits_{z = 0}^{{m_0} - 1} {z{r_z} + {m_0} - t} }}}}}  \\
  &\times \left( {{{\left( { - 1} \right)}^{t - 1}}\exp \left( {\frac{{{m_0}\left( {r + 1} \right)}}{{{a_q}\gamma {\Phi _{SD}}}}} \right){\rm{Ei}}\left( { - \frac{{{m_0}\left( {r + 1} \right)}}{{{a_q}\gamma {\Phi _{SD}}}}} \right)} \right. \\
 &\left. { + \sum\nolimits_{o = 1}^t {\left( {o - 1} \right)!{{\left( { - 1} \right)}^{t - o}}{{\left( {\frac{{{m_0}\left( {r + 1} \right)}}{{{a_q}\gamma {\Phi _{SD}}}}} \right)}^{ - o}}} } \right). \\
\end{split}
\end{equation}
where ${{\rm{E}}_{\rm{i}}}(x) = \int_x^\infty  {{t^{ - 1}}} \exp ( - t)dt$ denotes the exponential integral function [22, (5.1.1)]. Secondly, In order to calculate $J_6$, we need to obtain the asymptotic expression as
\begin{equation}\label{eqn01}
\begin{split}
{J_6} &= {\rm{E}}\left( {{{\log }_2}\left( {1 + \frac{{{a_q}{\gamma ^2}{G_{S{R_{\bar m}}}}{G_{{R_{\bar m}}{D_q}}}}}{{\gamma \left( {{G_{S{R_{\bar m}}}} + {G_{{R_{\bar m}}{D_q}}}} \right) + \gamma {G_{{R_{\bar m}}{D_q}}}{\gamma _{_{\bar m}}} + {\gamma _{_{\bar m}}} + 1}}} \right)} \right)\\
 &= {\rm{E}}\left( {{{\log }_2}\left( {1 + \frac{{{a_q}{\gamma ^2}\frac{{{G_{S{R_{\bar m}}}}}}{{{\gamma _{_{\bar m}}} + 1}}{G_{{R_{\bar m}}{D_q}}}}}{{\frac{{\gamma {G_{S{R_{\bar m}}}}}}{{{\gamma _{_{\bar m}}} + 1}} + \gamma {G_{{R_{\bar m}}{D_q}}} + 1}}} \right)} \right) \\&< {\rm{E}}\left( {{{\log }_2}\left( {1 + {a_q}\gamma \min \left( {{G_{{R_{\bar m}}{D_q}}},\frac{{{G_{S{R_{\bar m}}}}}}{{{\gamma _{_{\bar m}}} + 1}}} \right)} \right)} \right).
\end{split}
\end{equation}
Now denoting a random variable
\begin{equation}\label{eqn01}
\begin{split}
L = {a_q}\gamma \min \left( {{G_{{R_{\bar m}}{D_q}}},\frac{{{G_{S{R_{\bar m}}}}}}{{{\gamma _{_{\bar m}}} + 1}}} \right).
\end{split}
\end{equation}
Thus an upper bound of $J_6$ can be rewritten as
\begin{equation}\label{eqn01}
\begin{split}
{\rm{E}}\left( {{{\log }_2}\left( {1 + {\gamma _{{R_{\bar m}}{D_q}}}} \right)} \right) < {\rm{E}}\left( {\frac{1}{{\ln 2}}\ln \left( {1 + l} \right)} \right) = \frac{1}{{\ln 2}}\int_0^\infty  {\frac{{1 - {F_L}\left( l \right)}}{{1 + l}}dl}.
\end{split}
\end{equation}
Submitting the PDFs of ${G_{S{R_{\bar m}}}}$ and ${G_{{R_{\bar m}}{D_q}}}$ into (51) and using [22, (3.353.5)], we can get
\begin{equation}\label{eqn01}
\begin{split}
 {J_6} <& \frac{{Z_1^\Omega }}{{\ln 2\Gamma ({m_1})}}\sum\limits_{n = 0}^{{m_2} - 1} {\frac{{{{\left( {{{{m_2}} \mathord{\left/
 {\vphantom {{{m_2}} {{\Phi _{RD}}}}} \right.
 \kern-\nulldelimiterspace} {{\Phi _{RD}}}}{a_q}\gamma } \right)}^n}}}{{n!}}} \sum\nolimits_{j = 0}^{\Omega  - 1} {{{\left( { - 1} \right)}^j}}
   {{\Omega  - 1} \choose j}  \\
  &\times \sum\nolimits_{r = 0}^j {\sum\limits_{{r_i} \ge 0,{r_0} + {r_1} + ... + {r_{{m_1} - 1}} = r} {\frac{{{{\left( { - 1} \right)}^r}}}{{\prod\nolimits_{z = 0}^{{m_1} - 1} {{{\left( {z!} \right)}^{{r_z}}}} }}} }  \\
  &\times \sum\nolimits_{o = 0}^{^{\sum\nolimits_{z = 0}^{{m_1} - 1} {z{r_z}}  + {m_1} - 1}} {\frac{{\left( {\sum\nolimits_{z = 0}^{{m_1} - 1} {z{r_z}}  + {m_1} - 1} \right)!{{\left( {\left( {{\gamma _{_{\bar m}}} + 1} \right)} \right)}^o}{\Phi _{SR}}^{ - o}}}{{o!{{\left( {{a_q}\gamma } \right)}^o}{{\left( {\left( {r + 1} \right)} \right)}^{\sum\nolimits_{z = 0}^{{m_1} - 1} {z{r_z} + {m_1}}  - o}}{m_1}^{ - o}}}}  \\
  &\times \left( {{{\left( { - 1} \right)}^{n + o - 1}}\exp \left( {\frac{{{m_1}\left( {{\gamma _{_{\bar m}}} + 1} \right)\left( {r + 1} \right)}}{{{\Phi _{SR}}{a_q}\gamma }} + \frac{{{m_2}}}{{{a_q}\gamma {\Phi _{RD}}}}} \right)} \right. \\
  &\times {\rm{Ei}}\left( { - \frac{{{m_1}\left( {{\gamma _{_{\bar m}}} + 1} \right)\left( {r + 1} \right)}}{{{\Phi _{SR}}{a_q}\gamma }} - \frac{{{m_2}}}{{{a_q}\gamma {\Phi _{RD}}}}} \right) \\
 &\left. { + \sum\nolimits_{\varepsilon  = 1}^{n + o} {\left( {\varepsilon  - 1} \right)} !{{\left( { - 1} \right)}^{n + o - \varepsilon }}{{\left( {\frac{{{m_1}\left( {{\gamma _{_{\bar m}}} + 1} \right)\left( {r + 1} \right)}}{{{\Phi _{SR}}{a_q}\gamma }} + \frac{{{m_2}}}{{{a_q}\gamma {\Phi _{RD}}}}} \right)}^{ - \varepsilon }}} \right). \\
\end{split}
\end{equation}
Thirdly, $J_7$ can be rewritten as
\begin{equation}\label{eqn01}
\begin{split}
 {J_7}& = {\rm{E}}\left( {{{\log }_2}\left( {1 + \max \left( {{a_q}\gamma {G_{S{D_q}}},\frac{{{a_q}{\gamma ^2}\frac{{{G_{S{R_{\bar m}}}}}}{{{\gamma _{{R_{\bar m}}{R_{\bar m}}}} + 1}}{G_{{R_{\bar m}}{D_q}}}}}{{\gamma \frac{{{G_{S{R_{\bar m}}}}}}{{{\gamma _{{R_{\bar m}}{R_{\bar m}}}} + 1}} + \gamma {G_{{R_{\bar m}}{D_q}}} + 1}}} \right)} \right)} \right) \\
  &< {\rm{E}}\left( {{{\log }_2}\left( {1 + \max \left( {{a_q}\gamma {G_{S{D_q}}},{a_q}\gamma \min \left( {\frac{{{G_{S{R_{\bar m}}}}}}{{{\gamma _{{R_{\bar m}}{R_{\bar m}}}} + 1}},{G_{{R_{\bar m}}{D_q}}}} \right)} \right)} \right)} \right). \\
\end{split}
\end{equation}
As the derivation process of (49), we denote the random variable
\begin{equation}\label{eqn01}
\begin{split}
W= \max \left( {{a_q}\gamma {h_{S{D_q}}},{a_q}\gamma \min \left( {\frac{{{G_{S{R_{\bar m}}}}}}{{{\gamma _{{R_{\bar m}}{R_{\bar m}}}} + 1}},{G_{{R_{\bar m}}{D_q}}}} \right)} \right).
\end{split}
\end{equation}
Then, (53) can be obtained as
\begin{equation}\label{eqn01}
\begin{split}
{\rm{E}}\left( {{{\log }_2}\left( {1 + \max \left( {{\gamma _{S{D_q}}},{\gamma _{{R_{\bar m}}{D_q}}}} \right)} \right)} \right){\rm{ < E}}\left( {\frac{1}{{\ln 2}}\ln \left( {1 +W} \right)} \right) = \frac{1}{{\ln 2}}\int_0^\infty  {\frac{{1 - {F_W}\left( w \right)}}{{1 + w}}dw}.
\end{split}
\end{equation}
The derivation of (55) is the same as (49), so we omit it here. By submitting (48), (52), (55) and all the probabilities into (46), we can get the close form expression of $R_{ave}^q$. With the $R_{ave}^p$ and $R_{ave}^q$, we can get the expression of ${R_{sum}}$.

\section{Results and discussion}
In this section, Monto-Carlo simulations are performed to validate our theoretical analysis. In all the simulations, the energy conversion efficiency $\eta  = 0.5$, the fixed transmission rate of the source is 1 bit per channel use (bpcu). We set $C$ as the multiple of the source transmission energy, i.e., $C = \delta PT$, where $\delta  > 0$ is the scaling factor. Unless stated otherwise, the following sets of parameters are used: $N=3$, ${\sigma ^2} = 1$, ${m_0} = {m_1} = {m_2} = 2$, ${\Phi _{SD}} = {\Phi _{SR}} = {\Phi _{RD}} = 1$, ${a_p} = {8 \mathord{\left/
 {\vphantom {8 9}} \right.
 \kern-\nulldelimiterspace} 9}$, ${a_q} = {1 \mathord{\left/
 {\vphantom {1 9}} \right.
 \kern-\nulldelimiterspace} 9}$, ${\gamma _{thp}} = 1{\rm{dB}}$, ${\gamma _{thq}} = 2{\rm{dB}}$, $p=1$, $q=3$.

\begin{figure} [htbp]
\centering \scalebox{0.6}{\includegraphics{2.eps}}
\caption{Outage probability with  different scheme (SI=-10dB, $\delta =0.5$, $M$=6, $B$  = 10).\label{Fig4}}
\end{figure}

\subsection{Results}

\subsubsection{Numerical and simulation results of outage behavior}
 Fig. 2 shows the performance of our proposed scheme versus the transmitting SNR of source. We can see that the approximate analysis of finite Markov chain is accurate in the whole SNRs. Compared with the HD and conventional OMA relaying scheme, the outage behavior of each terminal in the FD relaying system is better, even in the existence of SI at the relay. We can also find that the diversity order for each user is equal to that in the HD cooperative NOMA system. Since the theoretical analyses agree well with the simulations in medium and high SNR ranges, we will only plot the analytical results in the remaining outage probability figures.

In Fig. 3, we investigate the impacts of SI on the outage probability. We can see that that the outage performance of the FD relaying (only take one user, i.e. $D_p$, e.g. and vice versa) system will become worse than that in the HD relaying network for high-SNR values when SI is extremely severe. But both two schemes can have better outage performance than using the conventional OMA scheme. Also, two users can be served simultaneously. Therefore, the user fairness of the FD cooperative NOMA system can be improved.
\begin{figure} [htbp]
\centering \scalebox{0.6}{\includegraphics{3.eps}}
\caption{Outage probability with different SI ($\delta  = 0.5$, $M$=6, $B$= 10).\label{Fig4}}
\end{figure}

In Fig. 4, we illustrate the impacts of $\delta$ on the performance of the proposed scheme with different SI in medium SNR conditions (${P \mathord{\left/
 {\vphantom {P {{\sigma ^2}}}} \right.
 \kern-\nulldelimiterspace} {{\sigma ^2}}}$= 20dB).  For all the curves, the trends are the same for all schemes, namely, the probability first decreases then increases as $\delta$ varies from 0 to 1. This means that, when the other parameters are determined, there must be an optimal value of $\delta$. However, the values of the inflection points are not always the same for different users and SI. We can see that for the user with better channel condition ($D_q$), the optimal value is smaller than user with worse channel condition ($D_p$), this is because $D_p$ has the worse direct link, so it needs more power to enhance the relay link. We can also see that for the FD cooperative NOMA system, the optimal value shifts to the left with the increase in the level of SI and when $\delta$ is too big, the outage behavior is nearly the same under different values of SI. More importantly, compared to the HD cooperative NOMA system, the outage behavior of $D_p$ is getting inferior since the value of SI is large. However, the optimal value of $\delta$ can easily be obtained by a one-dimensional exhaustive search, with this optimal value of $\delta$, the system can resist fading more effectively.
\begin{figure} [htbp]
\centering \scalebox{0.6}{\includegraphics{4.eps}}
\caption{ Outage probability against ($M$=6, $B=10$).\label{Fig4}}
\end{figure}

\subsubsection{Numerical and simulation results of  ergodic sum-rate}
Fig. 5 presents the rates of the two paired users varies with system SNR. Firstly, the sum-rate increases with the growth of transmitting SNR, and asymptotic curves of the ergodic sum-rate and simulation results are almost overlapping in the medium- and high-SNR regions, we will only plot the analytical results in the remaining outage rate figures. Secondly, the results clearly show that the FD cooperative NOMA system can achieve a larger sum-rate compared with the HD relaying system and conventional OMA system, but the double rate cannot be achieved owing to the effect of SI. More importantly, with the increase of system SNR, the rate of $D_p$ approaches a certain value, whereas $D_q$ can obtain improved performance. Both the two conclusions are consistent with the derivations in (45) and (46).

\begin{figure} [htbp]
\centering \scalebox{0.6}{\includegraphics{5.eps}}
\caption{Rate with different user and different scheme (SI=-10dB, $\delta =0.5$, $M$=6, $B= 10$).\label{Fig4}}
\end{figure}

In Fig. 6, we investigate the impact of battery size on the proposed finite battery scheme by varying the relay number and SI. we can observe that the performance increases as $B$ increases. However, the gain provided by a larger battery size does not increase when exceeds a certain value and this value is approximate to 6. From this figure, we can also see that the battery can not offer enough energy for forwarding the date when $B=1$, so the performances are the same for all curves and the performance gap for different SI  increases as the relay number increases.

\begin{figure} [htbp]
\centering \scalebox{0.6}{\includegraphics{6.eps}}
\caption{Sum-rate against $B$ ($\delta=0.5$).\label{Fig4}}
\end{figure}

\subsection{Discussion}
Compared with existing literatures, the channel in this paper we use is more in line with the actual communication scenario
and the model we propose combines NOMA, cooperative communication, FD and EH together reasonably and
can improve spectral efficiency and energy efficiency simultaneously. However, there are still
some shortages in our proposed scheme. The joint optimization of multiple parameters such as the transmission
power, power splitting factor and user pairing has not been mentioned in this paper, this joint optimization can
make the proposed schene more useful in practical application which we will study in the future.

\section{Conclusions}
In this paper, we have proposed a relay selection scheme for the FD-NOMA networks with EH over Nakagami-$m$ fading channels. We modeled the amount of harvested energy at each relay using a finite Markov chain and then derived the approximate closed-form expression of the outage probability and ergodic sum-rate. Simulations were carried out to verify the effectiveness of our theoretical analysis. We concluded that through carefully choosing the system parameters of the network, (e.g., EH threshold), FD-NOMA performs much better than the HD-NOMA or the conventional OMA.
%\nocite{oreg,schn,pond,smith,marg,hunn,advi,koha,mouse}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Backmatter begins here                   %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{backmatter}

\section*{Abbreviations}
FD: full-duplex; NOMA: Non-orthogonal multiple access; CSI: Channel state information;
SI: Self-interference; OMA: Orthogonal multiple access; SIC: Successive interference cancellation;
MMSE: Minimum-mean-squared-error; T-ARQ: Truncated automatic repeat-request;
AF: Amplify-and-forward; DF: Decode-and-forward; EH: Energy harvesting; RF: Radio frequency;
HD: Half-duplex; INR: Interference-to-noise ratio; AWGN: Additive white Gaussian noise;
SNR: Signal-to-noise-ratio; SINR: Signal-to-interference-plus-noise-ratio; SC: selection combining;
PDF:  Probability density function; CDF: Cumulative distribution function

\section*{Acknowledgements}
Not applicable.

\section*{Availability of data and materials}
  Not available online. Please contact corresponding author for data requests.


\section*{Funding}
This work was supported by the Shandong Provincial Natural Science Foundation,
China (Grant No. ZR2017MD019).

\section*{Author's contributions}
WG proposed the system model, analyzed the performance and drafted the paper. LS
participated in the simulations and helped to draft the manuscript. All authors
read and approved the final manuscript.

\section*{Competing interests}
  The authors declare that they have no competing interests.






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  The Bibliography                       %%
%%                                                         %%
%%  Bmc_mathpys.bst  will be used to                       %%
%%  create a .BBL file for submission.                     %%
%%  After submission of the .TEX file,                     %%
%%  you will be prompted to submit your .BBL file.         %%
%%                                                         %%
%%                                                         %%
%%  Note that the displayed Bibliography will not          %%
%%  necessarily be rendered by Latex exactly as specified  %%
%%  in the online Instructions for Authors.                %%
%%                                                         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% if your bibliography is in bibtex format, use those commands:
%\bibliographystyle{bmc-mathphys} % Style BST file (bmc-mathphys, vancouver, spbasic).
%\bibliography{bmc_article}      % Bibliography file (usually '*.bib' )
% for author-year bibliography (bmc-mathphys or spbasic)
% a) write to bib file (bmc-mathphys only)
% @settings{label, options="nameyear"}
% b) uncomment next line
%\nocite{label}

% or include bibliography directly:
% \begin{thebibliography}
% \bibitem{b1}
% \end{thebibliography}

\begin{thebibliography}{1}

\bibitem{Monks2} Y. Saito, A. Benjebbour, Y. Kishiyama, T. Nakamura, System level performance evaluation of downlink non-orthogonal multiple access (NOMA).
%\hskip 1em plus 0.5em minus 0.4em
in Proceedings of the IEEE Annual International Symposium on Personal, In door and Mobile Radio Communications (PIMRC), London, United Kingdom, 2013.

\bibitem{Monks3} Z. Ding, Z. Yang, P. Fan, H. V. Poor, On the performance of non-orthogonal multiple access in 5G systems with randomly deployed users.
%\hskip 1em plus 0.5em minus 0.4em
 IEEE Signal Process. Lett. 21(12), 1501-1505 (2014).

\bibitem{Muqattash3}Y. Lan, A. Benjebboiu, X. Chen, A. Li, H. Jiang,  Considerations on downlink non-orthogonal multiple access (NOMA) combined with closed-loop SU-MIMO. in Proceedings of the IEEE International Conference on Signal Process and Communication System (ICSPCS), Gold Coast, Australia, 2014.

\bibitem{Muqattash3}X. Chen, A. Benjebbour, Y. Lan, A. Li, H. Jiang, Impact of rank optimization on downlink non-orthogonal multiple access (NOMA) with SU-MIMO. in Proceedings of the IEEE International Conference on Communication System (ICCS), Macau, China, 2014.

\bibitem{Krunz1}Z. Ding, M. Peng, H. V. Poor, Cooperative non-orthogonal multiple access in 5G systems. IEEE Commun. Lett. 19(8),1462-1465 (2015).

\bibitem{Li6}Z. Ding, H. Dai, H. V. Poor, Relay selection for cooperative NOMA. IEEE Wireless Commun. Lett. 5(4), 416-419 (2016).

\bibitem{Li6}Z. Yu, C. Zhai, W. Ni, D. Wang, Non-orthogonal multiple access with cooperative truncated ARQ and relay selection. IEEE Access, 7, 56228-56243 (2019).


\bibitem{10}Z. Zhang, X. Chai, K. Long, A. V. Vasilakos, L. Hanzo, Full duplex techniques for 5G networks: Self-interference cancellation, protocol design, and relay selection. IEEE Commun. Mag. 53(5), 128-137 (2015).

\bibitem{10}Q. Wang, Y. Dong, X. Xu, X. Tao, Outage probability of full-duplex AF relaying with processing delay and residual self-interference. IEEE Commun. Lett. 19(5), 783-786 (2015).

\bibitem{10}D. P. M. Osorio, E. E. B. Olivo, H. Alves, J. C. S. S. Filho, M. Latva-Aho, Exploiting the direct link in full-duplex amplify-and-forward relaying networks. IEEE Signal Processing Lett. 22(10), 1766-1770 (2015).

\bibitem{10}T. Kwon, S. Lim, S. Choi, D. Hong, Optimal duplex mode for DF relay in terms of the outage probability. IEEE Trans. Veh. Technol. 59(7), 3628-3634 (2010).

\bibitem{10}A. A. Nasir, X. Zhou, S. Durrani, R. A. Kennedy, Relaying protocols for wireless energy harvesting and information processing. IEEE Trans. Wireless Commun. 12(7), 3622-3636 (2013).

\bibitem{10}Z. Ding, S. M. Perlaza, I. Esnaola, H. V. Poor, Power allocation strategies in energy harvesting wireless cooperative networks. IEEE Trans Wireless Commun. 13(2), 846-860 (2014).
\bibitem{10}H. Gao, W. Ejaz, M. Jo, Cooperative Wireless Energy Harvesting and Spectrum Sharing in 5G Networks. IEEE Access. 4, 3647-3658 (2016).

\bibitem{10}S. Sudevalayam, P. Kulkarni, Energy harvesting sensor nodes: Survey and implications. IEEE Commun. Surveys Tuts. 13(3), 443-461 (2011).

\bibitem{10}C. Huang, J. Zhang, P. Zhang, S. Cui, Threshold-based transmissions for large relay networks powered by renewable energy. in Proceedings of the IEEE Global Communications Conferences (ICC), Budapest, Hungary, 2013.

\bibitem{10}Z. Zhou, M. Peng, Z. Zhao, W. Wang, R. S. Blum, Wireless-powered cooperative communications: Power-splitting relaying with energy accumulation. IEEE J. Select. Areas Commun. 34(4), 969-982 (2016).

\bibitem{10}C. Zhai, Z. Yu, X. Wang, High-efficient cooperative relaying with wireless powered source and relay. Computer Networks, 152, 199-209 (2019).

\bibitem{10}Z. Zhang, Z. Ma, Z. Ding, M. Xiao, G. K. Karagiannidis, Full-duplex two-way and one-way relaying: average rate, outage probability, and tradeoffs. IEEE Trans. Wireless Commun. 15(6), 3920-3933 (2016).

\bibitem{10}X. Yue, Y. Liu, S. Kang,  A. Nallanathan, Z. Ding, Exploiting full/Half-duplex user relaying in NOMA systems. IEEE Trans. Commun. 66(2), 560-575 (2018).
\bibitem{10}V. Rego, Characterization of equilibrium queue length distributions in M/GI/I queues. Computers and Operations Research, 15(1), 7-17 (1988).
\bibitem{10}I. Gradshteyn, I. Ryzhik, Table of integrals, series and products. Academic Press: San Diego, CA (1994).
\bibitem{10}H. A. David, H. N. Nagaraja, Order statistics. 3rd edn, John Wiley, New York (2003).

\end{thebibliography}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Figures                       %%
%%                               %%
%% NB: this is for captions and  %%
%% Titles. All graphics must be  %%
%% submitted separately and NOT  %%
%% included in the Tex document  %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
%% Do not use \listoffigures as most will included as separate files

%\section*{Figures}
%  \begin{figure}[h!]
%  \caption{\csentence{Sample figure title.}
%      A short description of the figure content
%      should go here.}
%      \end{figure}
%
%\begin{figure}[h!]
%  \caption{\csentence{Sample figure title.}
%      Figure legend text.}
%      \end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Tables                        %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Use of \listoftables is discouraged.
%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Additional Files              %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\end{backmatter}
\end{document}
