*Article* **Throughput Optimization for NOMA Cognitive Relay Network with RF Energy Harvesting Based on Improved Bat Algorithm**

**Yi Luo 1, Chenyang Wu 1, Yi Leng 2,\*, Nüshan Huang 1, Lingxi Mao <sup>1</sup> and Junhao Tang <sup>1</sup>**


**Abstract:** Due to the shortcomings of the standard bat algorithm (BA) for multi-parameter optimization, an improved bat algorithm is proposed. The benchmark function test shows that the proposed algorithm has better realization of high-dimensional function optimization by introducing multiple flight modes, adopting adaptive strategy based on group trend, and employing loudness mutation flight selection strategy based on Brownian motion. Aiming at the characteristics of complex networks structure and multiple design variables of energy harvesting non-orthogonal multiple access cognitive relay networks (EH-NOMA-CRNs), we utilize the proposed hybrid strategy improved bat algorithm (HSIBA) to optimize the performance of EH-NOMA-CRNs. At first, we construct a novel two-hop underlay power beacon assisted EH-NOMA-CRN, and derive the closed-form expressions of secondary network's outage probability and throughput. Then, the secondary network performance optimization is formulated as the throughput maximation problem with regard to EH ratio and power allocation factors. Subsequently, the HSIBA is employed to optimize the above parameters. Numerical results show that the proposed HSIBA can achieve optimization to the constructed EH-NOMA-CRN with faster convergence speed and higher stability.

**Keywords:** bat algorithm; hybrid strategy; energy harvesting; NOMA; cognitive relay network

**MSC:** 94A05

#### **1. Introduction**

Non-orthogonal multiple access (NOMA) is an effective method to increase frequency efficiency. It has drawn great attention for its promising applications in the fields of 5G, 6G and Internet of things (IoT) networks [1–3]. Especially in the military scenarios of unmanned aerial vehicles (UAVs) battlefield situation awareness [4], covert military operations and tactical area communications [5], UAV-NOMA relaying communications [6], and hybrid satellite-UAV networks [7], the NOMA technology has been widely used. The core concept of NOMA is that when the non-orthogonal transmission is adopted at the transmitter, the interference signals are actively introduced, and the correct demodulation is achieved at the receiver by applying serial interference cancellation (SIC) technique. It can be seen that reliable interference cancellation techniques and strategies [8,9] are the basis for realizing NOMA. Although the use of SIC technique will increase the complexity of the receiver, it can improve the spectrum efficiency.

In order to further improve the utilization efficiency of the authorized spectrum and the expand the coverage of networks, many scholars have combined NOMA technique with cognitive relay network (CRN) to build and optimize some novel NOMA-CRN models. In [10], the authors achieved in-band full-duplex and two-way cognitive relaying transmission in a cooperative NOMA system, and applied the successive inner approximation technique to maximize the energy efficiency. In [11], the error rate performance of an underlay NOMA-CRN with partial relay selection (PRS) scheme was studied, and the optimum

**Citation:** Luo, Y.; Wu, C.; Leng, Y.; Huang, N.; Mao, L.; Tang, J. Throughput Optimization for NOMA Cognitive Relay Network with RF Energy Harvesting Based on Improved Bat Algorithm. *Mathematics* **2022**, *10*, 4357. https://doi.org/10.3390/ math10224357

Academic Editor: Oliviu Matei

Received: 22 August 2022 Accepted: 11 November 2022 Published: 19 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

power coefficients were optimized to minimize the average bit error rate union bound by using numerical methods. In [12], the physical-layer security of NOMA-CRN with relay cooperation was investigated, and the particle swarm optimization (PSO) algorithm was adopted to optimize the power allocation to maximize the secrecy sum-rate. In [13], the authors utilized the NOMA-CRN for device to device (D2D) transmission, and proposed a deep neural networks framework for optimizing resource allocation to maximize the sum rate. In [14], the NOMA-CRN with a single primary transmitter/receiver (PT/PR) and PRS scheme was constructed, and the PSO algorithm was employed for determining the jointly optimal power allocation factors to realize the maximization of throughput.

In addition to extend the lifetime and improve the energy efficiency of energy-constrained relay networks using NOMA technique, the radio frequency energy harvesting (RF-EH) technique has been introduced into NOMA relay network to construct the new energy harvesting NOMA relay networks (EH-NOMA-RNs). Ref. [15] considered a cooperative EH-NOMA-RN, and adopted a one-dimensional search algorithm for optimizing power allocation coefficient to attain the maximization of weighted sum rate. Ref. [16] used the EH-NOMA-RN with interfering signal for IoT systems, and adopted the Golden section search method to maximize the sum-throughput. Ref. [17] applied EH-NOMA-RN system to layered video multicast transmission, and used the initial feasible point search algorithm to minimize the base station's average transmission power. Ref. [18] constructed a full-duplex EH-NOMA-RN with cooperative relaying, and made use of the alternating optimization technique to achieve the outage probability (OP) minimization and throughput maximization of the system. Similar to [16], ref. [19] also applied Golden section search algorithm for optimization of a power beacon (PB)-assisted EH-NOMA-RN IoT-based system to minimize OP and maximize sum-throughput of the system.

Recently, based on the research of NOMA-CRNs and EH-NOMA-RNs, integrating RF-EH technique into NOMA-CRN to form the brand-new energy harvesting NOMA cognitive relay networks (EH-NOMA-CRNs) and solving the corresponding networks' performance optimization problems have gradually become a novel research direction. In [20], the authors employed a two-loop procedure using one-dimensional search to solve the minimum transmission power problem of a multiple-input single-output EH-NOMA-CRN with non-linear EH model and robust beamforming. In [21], two-level bisection search algorithms were developed to realize the maximal secondary throughput of EH-NOMA-CRN by optimizing resource allocation coefficients. In [22], a joint optimization algorithm based on one-dimensional search was proposed to maximize the energy efficiency of EH-NOMA-CRN with a discrete EH scheme. In [23], the physical layer security of EH-NOMA-CRN with multi-input multi-output two-way relaying was investigated, and a joint path-following-based optimization algorithm was put to use for maximizing the sum achievable secrecy rate. Later, ref. [24] designed two iterative optimization algorithms based on a deep neural network frame for the ergodic capacity prediction of IoT PB assisted EH-NOMA-CRN to minimize OP users and maximizing system throughput.

According to the existing researches, it is known as follows:


and then convex optimization is performed. However, the transformation process is more complicated. (iii) Algorithms based on deep learning [13,24] are introduced into system performance optimization. Nevertheless, the deep neural networks framework requires stronger computing power, which is not very suitable for energy-constrained EH-NOMA-CRNs' nodes with lower storage and computing capacity. (iv) The metaheuristic algorithm based on PSO is initially used to solve the system optimization problem [12,14]. Although the PSO algorithm is simple, its performance is poor. It is worth studying to design a new high-performance meta-heuristic algorithm suitable for performance optimization of EH-NOMA-CRN.


To the best of our knowledge, there is no existing paper studying meta-heuristic algorithms other than PSO algorithm for system performance optimization of EH-NOMA-CRNs, which motivates us to write this treatise. Compared with other existing swarm intelligence algorithms, Bat algorithm (BA) has many advantages, such as implementation simplicity, less control parameters and excellent global ability. Since the BA was proposed, its improvement has been carried out continuously. Especially in the last two years, some novel improved BA algorithms have been proposed successively, such as the enhanced Levy flight bat algorithm (ELBA) [25], and the improved BA with extremal optimization algorithm (IBA-EO) [26]. Therefore, how to further enhance the global search ability and convergence rate of improved BA for coping with different engineering application scenarios is still a problem worth studying. Like other meta-heuristic algorithms, BA has some inherent deficiencies. Due to a single flight mode that has limited ability to search the solution space, the accuracy of the algorithm is insufficient, and it is easy to fall into the local optimal solution in high-dimensional search space. In this paper, a novel improved BA has been proposed.

Explicitly, the major contributions of our paper are summarized as follows:


The remainder of our paper is organized as follows. Section 2 introduces the standard BA and HSIBA, and verifies the performance of HSIBA. Section 3 describes the EH-NOMA-CRN model, and derives the closed-form solutions of SN's OP and throughput. In Section 4, the HSIBA is adopted to optimize the throughput of SN, and simulation results are shown. Then, the conclusions are presented in Section 5.

#### **2. Hybrid Strategy Improved Bat Algorithm**

#### *2.1. Standard Bat Algorithm*

The bat algorithm (BA) proposed by Xin-She Yang is a novel metaheuristic algorithm [27], which adopts the echolocation of microbats. The standard BA is based on the following three idealized rules:


Originally, the *i*-th bat is randomly given a frequency obeying uniform distribution from *fmin* to *fmax*. In the *t*-th iteration, the new position *x<sup>t</sup> <sup>i</sup>* and velocity *<sup>v</sup><sup>t</sup> <sup>i</sup>* of the *i*-th bat can be respectively calculated as

$$f\_i = f\_{min} + (f\_{max} - f\_{min})\beta \tag{1}$$

$$v\_i^t = v\_i^{t-1} + \left(\mathbf{x}\_i^{t-1} - \mathbf{x}^\*\right) f\_i \tag{2}$$

$$\mathbf{x}\_i^t = \mathbf{x}\_i^{t-1} + \mathbf{v}\_i^t \tag{3}$$

where *β* ∈ [0, 1] is a random vector obeying uniform distribution, and *x*<sup>∗</sup> is the current global best location (solution) determined after comparing the locations of all *n* bats.

During the local search process, once a solution is chosen from the current best solutions, a new solution for each bat can be generated by

$$
\varkappa\_{\text{new}} = \varkappa\_{\text{old}} + \varepsilon A^{\text{t}} \tag{4}
$$

where *<sup>ε</sup>* <sup>∈</sup> [−1, 1] is a random scaling factor, and *<sup>A</sup><sup>t</sup>* is the mean loudness of all *<sup>n</sup>* bats during the *t*-th iteration.

Furthermore, in the (*t*+1)-th iteration, the *i*-th bat's loudness *At*+<sup>1</sup> *<sup>i</sup>* and rate of pulse emission *rt*+<sup>1</sup> *<sup>i</sup>* can be respectively updated as

$$A\_i^{t+1} = \theta A\_i^t \tag{5}$$

$$r\_i^{t+1} = r\_i^0 \left(1 - \mathbf{e}^{-qt}\right) \tag{6}$$

where *θ* and *ϕ* are constants. For simplicity, *θ* and *ϕ* are respectively set as 0.8 and 0.9 in our paper. In order to clearly describe the optimization process of standard BA, the pseudocode of standard BA is given in Algorithm 1.

#### *2.2. Algorithm Improvement Based on Hybrid Strategy*

The standard BA has a fast convergence speed, which can adjust and balance the local search and global search in the optimization process by gradually increasing the pulse transmission frequency *ri* and reducing the pulse loudness *Ai*. However, the accuracy of the algorithm is insufficient due to a single flight mode that has limited ability to search the solution space. In our paper, we adopt adaptive strategy based on group trend, add several different flight modes, and introduce loudness mutation factor to select between the proposed modes as the bat flight mode during different search stages. This method can greatly improve the accuracy and astringency of the hybrid strategy improved bat algorithm.

#### **Algorithm 1** Pseudocode of standard BA


#### 2.2.1. Adaptive Strategy Based on Group Trend

For the standard BA, bat individuals fly freely in a wide area independently of each other in the process of solving. Due to the lack of group awareness to adapt to the current situation, it is difficult to balance the distance of each flight in the iterative process, and it is easy to miss the optimal solution, which is not conducive to local optimization. For example, it is assumed that the *i*-th bat is selected for local optimization at the global optimization in the *t*-th iteration. For the (*t*+1)-th iteration, if its position update is still carried out according to a large step size, it is easy to make it out of the optimal position, which is not conducive to local optimization. Therefore, HSIBA introduces an adaptive strategy based on group trend to change flight weight *ω* which can update the bat's flight weight *ω* according to the current search situation. To describe the group trend, *B<sup>t</sup> <sup>i</sup>* and *<sup>O</sup><sup>t</sup>* are respectively defined as

$$B\_i^t = \begin{cases} 1 & if\left(fit(\mathbf{x}\_i^{t-1}) < fit(\mathbf{x}\_i^t)\right) \\ 0 & otherwise \end{cases} \tag{7}$$

where *fit*(·) denotes fitness operation.

$$O^t = \frac{\sum\_{i=1}^{I} B\_i^{t-1}}{I} \tag{8}$$

When there are more bat individuals flying in the better direction in the population, the value of *O<sup>t</sup>* is larger. The adaptive flight weight *ω* can be calculated as

$$
\omega = \omega\_{\rm min} + (\omega\_{\rm max} - \omega\_{\rm min})\mathbf{O}^t \tag{9}
$$

When there are more individuals flying in the favorable direction in the population, the value of *ω* is larger, so that bats can adopt a larger velocity weight to expand the search range and prevent the algorithm from falling into the local optimal solution; when most individuals in the population are difficult to find a better solution, the value of *ω* is smaller to avoid missing the optimal solution, which will effectively speed up the convergence speed of the algorithm.

#### 2.2.2. Introducing Multiple Flight Modes

To change the single flight mode such as Equation (3) adopted in the standard BA, four flight modes are introduced in HSIBA. The four flight modes are adjusted by adaptive velocity weights *ω*, and selected by loudness mutation factor *β<sup>i</sup>* (see Section 2.2.3) and random number *rand* with uniform distribution of 0∼1.

(1). Adaptive weight bit flight mode

The position update Equation (3) in the standard BA is changed to Equation (10), and the bat's flying speed is adjusted by adaptive flight weight.

$$\mathbf{x}\_{i}^{t+1} = \mathbf{x}\_{i}^{t} + \boldsymbol{\omega} \cdot \mathbf{v}\_{i}^{t+1} \tag{10}$$

When |*βi*| < 0.1 and *rand* < 0.4, the bats adopt the adaptive weight bit flight mode and update their positions by using Equation (10).

(2). Position exchange variant flight mode

In the standard BA, all bats only take the optimal individual of the current population as the flight direction, and there is no information exchange between bat individuals. Therefore, in the HSIBA, position exchange variant flight is innovatively introduced to enhance information exchange among bat individuals and obtain more information about the feasible solution space. In position-exchange variant flight mode, the bats' positions are updated as follows

$$\mathbf{x}\_{i}^{t+1} = \boldsymbol{\omega} \times \left| 2 \times rand \times \mathbf{x}^{\*} - \mathbf{x}\_{I}^{t} \right| \times \left( \beta\_{i}^{t} + rand \right), \quad i = 1, \tag{11}$$

$$\mathbf{x}\_{i}^{t+1} = \boldsymbol{\omega} \times \left| \begin{array}{c} 2 \times \quad rand \ \times \mathbf{x}^{\*} - \mathbf{x}\_{i-1}^{t+1} \mid \times \left( \beta\_{i}^{t} + \quad rand \ \right) , \quad i \ge 2 \end{array} \right. \tag{12}$$

The first bat performs a variant flight based on the *I*-th bat's position of the *t*-th iteration, and obtains the position of the (*t*+1)-th iteration through Equation (11), and the *i*-th bat performs a variant flight based on the (*i*-1)-th bat's position of the (*t*+1)-th iteration, and obtains the position of the (*t*+1)-th iteration through Equation (12). When |*βi*| ≥ 0.1 and *rand* < 0.6, the bats apply the position exchange variant flight

mode and update their positions by using Equations (11) and (2).

(3). Rotary flight mode

In [28], the African vulture's rotating flight foraging model was proposed and simulated. Inspired by this model, the rotary flight mode is introduced to our HSIBA, which makes the bats rotate around the optimal individual of the population to expand the search range and improve the ability of jumping out the local optimal solution. The rotation flight vectors are respectively expressed as

$$\phi\_1 = \mathbf{x}^\* \times \left(\frac{rand \times \mathbf{x}\_i^t}{2\pi}\right) \times \cos\left(\mathbf{x}\_i^t\right) \tag{13}$$

$$\phi\_2 = \mathbf{x}^\* \times \left(\frac{rand \times \mathbf{x}\_i^t}{2\pi}\right) \times \sin\left(\mathbf{x}\_i^t\right) \tag{14}$$

and the *i*-th bats' position is updated as

$$\mathbf{x}\_{i}^{t+1} = \mathbf{x}^\* - \boldsymbol{\omega} \cdot (\boldsymbol{\phi}\_1 + \boldsymbol{\phi}\_2) \tag{15}$$

When |*βi*| ≥ 0.1 and *rand* ≥ 0.6, the bats employ the rotary flight mode and update their positions by using Equation (15).

(4). Levy flight mode

In order to ensure that every bat can catch food in the process of predation, it is sure that each bat does not fly alone. If there is no companion around the *i*-th bat, it adopts the random walk scheme to update the step. The *i*-th bat can apply Levy flight [29] mode, its position can be updated as

$$\mathbf{x}\_{i}^{t+1} = \mathbf{x}^\* - \left| \mathbf{x}^\* - \mathbf{x}\_i^t \right| \times \boldsymbol{\omega} \times L \boldsymbol{x} \mathbf{y}(d) \tag{16}$$

where *d* is the dimension of the *i*-th bat's position vector, and Levy function can be calculated as

$$Levy(\mathbf{x}) = 0.01 \times \frac{\mu\_1 \times \varepsilon}{|\mu\_2|^{\frac{1}{\Phi}}} \tag{17}$$

where *μ*<sup>1</sup> and *μ*<sup>2</sup> are random numbers on the interval of 0∼1, *ψ* is a constant, and *ε* can be computed as

$$\varepsilon = \left( \frac{\Gamma(1+\psi) \times \sin\left(\frac{\pi\beta}{2}\right)}{\Gamma\left(\frac{1+\psi}{2}\right) \times \psi \times 2^{\frac{\psi-1}{2}}} \right)^{\frac{1}{\psi}} \tag{18}$$

where Γ(*x*) = (*x* − 1)!

When |*βi*| < 0.1 and *rand* ≥ 0.4, the bats use Levy flight mode and update their positions by applying Equation (16).

#### 2.2.3. Loudness Mutation Flight Selection Strategy Based on Brownian Motion

This strategy is inspired by the physical phenomenon that the random motion of particles decreases with the temperature reduction in Brownian motion. In the HSIBA, the decreasing loudness *A* is understood as the decreasing temperature, and the loudness mutation factor *β* is the intensity of random motion that decreases with the decrease of temperature. During the *t*-th iteration, the loudness mutation factor of the *i*-th bat can be expressed as

$$
\beta\_i^t = A\_i^t \times (randn - 1) \tag{19}
$$

where *randn* is a random number that obeys a normal distribution with a mean of 0 and a variance of 1. In the search process, the pulse loudness *Ai* gradually decreases, and the loudness mutation factor *β<sup>i</sup>* dynamically converges to zero.

In the HSIBA, the flight mode employed by the *i*-th bat is decided by *βi*. When |*βi*| ≥ 0.1, the algorithm is in the early stage of the search. The bat tends to fly with more intense motion, and update its position by Equations (11), (12) and (15), which can achieve faster convergence speed; When |*βi*| < 0.1, the algorithm is at the end of the search. The bat tends to fly in a flight mode with small intensity and higher degree of randomness, and update its position through Equations (10) and (16), which can improve the ability of jumping out of the local optimal solution.

The improved bat algorithm fully takes into account the adjustment of the flight mode of bats in different search stages, which effectively improves the shortage of single flight mode used by standard BA in the optimization process. In order to conveniently understand the advantages of the HSIBA, its pseudocode and optimization flow chart are respectively shown in Algorithm 2 and Figure 1.

#### *2.3. Algorithm Performance Verification*

In order to verify the performance of HSIBA, sixteen frequently-used benchmark functions are selected as test functions. The test functions are shown in Table 1, including single-peak and multi-peak functions. At the same time, advanced African vulture optimization algorithm (AVOA) [28], the wild horse optimizer (WHO) algorithm [30], the arithmetic optimization algorithm (AOA) [31], hunter-prey optimization (HPO) algorithm [32], enhanced Lévy flight bat algorithm (ELBA) [26] and standard BA are selected for comparative analysis. The parameter settings of above algorithms and HSIBA are shown in Table 2. The number of all algorithm populations tested is 30, the maximum number of iterations is 100, and each test function runs 30 times independently, the test results are shown in Tables 3 and 4, and the iteration diagram and total time for each algorithm are

shown in Figures 2 and 3 . (The algorithms have been tested and executed using the Matlab 9.0 (R2016a) on a laptop computer, which runs Windows 10 Enterprise 64-bit with an AMD Ryzen 7 4800H, Radeon Graphics 2.90 GHz processor, and 16.00 GB RAM).


As shown in Tables 3 and 4, the test results and the performance rank of seven algorithms under 16 test functions are given. By comparing the test results in Tables 3 and 4, it can be seen that the optimization performance of HSIBA is far better than standard BA on single-peak and multi-peak test functions, and the HSIBA also has higher accuracy and stability than other advanced optimization algorithms. Moreover, by comparing the test results of all algorithms in 30 and 100 dimensions, the HSIBA has superior performance in solving high-dimensional optimization problems. Therefore, HSIBA is highly capable of maintaining balance in exploration and exploitation against large-scale issues.

In order to further compare the performance of different optimization algorithms, as shown in Figure 2 the iterative curves of the five optimization algorithms under 8 test functions are given. By the comparison of the five algorithm iterations curves, it can be seen that the HSIBA achieves better results with a minimum number of iterations in all test functions. Consequently, when solving the same problem, HSIBA can find a better solution quickly, so as to save computing time and power. Meanwhile, it can be seen that the fitness calculated by HSIBA decreases steadily with the increase of the number of iterations, which can effectively avoid the abrupt change of the convergence speed and falling into the local optimal value.

Figure 3 illustrates the computing time of seven algorithms. The results show that HSIBA does not significantly increase the complexity of the operation, but can improve the convergence speed and accuracy. The main reason is that the HSIBA adopts advanced search strategies and a variety of flight modes for selection.

**Figure 1.** Optimization flow chart of HSIBA.




**Table 2.** Parameter settings of optimization algorithms for comparison and evaluation of the HSIBA.

**Figure 2.** Iteration diagram of different optimization algorithms. (**a**) F1. (**b**) F4. (**c**) F5. (**d**) F6. (**e**) F9. (**f**) F10. (**g**) F11. (**h**) F12.


**Table 3.** Test results (*d* = 30).

**Figure 3.** Computing time of different algorithms.


**Table 4.** Test results (*d* = 100).

#### **3. Throughput Optimization for EH-NOMA-CRN**

#### *3.1. System Model*

As shown in Figure 4, we consider an underlay PB-assisted EH-NOMA-CRN, where primary network (PN) includes a set of *N* PR nodes, *Qn*, *n* ∈ {1, 2, ···, *N*}, and SN consists of a secondary source node S that employs the power domain NOMA technique to communicate with two secondary destination nodes D1 and D2 with the help of a secondary decode-and-forward (DF) relay node R. Similar to [11,12,24], It is assumed that the interference from PN to SN is neglected. It is also assumed that all nodes except the PB node W in our system model are equipped with a single antenna, and all SN nodes working in half-duplex mode only use the energy harvested from RF signals of W installed *M* antennas. All wireless channels are assumed to undergo quasi-static independent Rayleigh flat fading where each channel coefficient keeps constant during a frame but varies independently between different frames.

Let *hj*,*<sup>k</sup>* denote the fading channel coefficient corresponding to the link from node *j* to *k*, where *j* ∈ {S,R,W}, *k* ∈ {Q*n*,R, D1, D2}, and *j* = *k*. Accordingly, The channel power gain *hj*,*<sup>k</sup>* 2 is exponentially distributed with expectation value <sup>1</sup> *λj*,*<sup>k</sup>* ; for example, *<sup>λ</sup>j*,*<sup>k</sup>* <sup>=</sup> *<sup>d</sup><sup>ξ</sup> j*,*k*, where *dj*,*<sup>k</sup>* is the distance between node *j* and *k*, and *ξ* is the path loss factor. Moreover, it is assumed that all D*<sup>n</sup>* are closely located in one center point. Therefore, *h*S,Q*<sup>n</sup>* <sup>2</sup> and *h*R,Q*<sup>n</sup>* 2 are independent identically distributed random variables, respectively, i.e., *λ*S,Q*<sup>n</sup>* = *λ*S,Q, and *λ*R,Q*<sup>n</sup>* = *λ*R,Q.

Similar to [33], the three-phase time division broadcasting protocol with precise synchronization is adopted in our system model. Node *l*, *l* ∈ {S,R} simultaneously harvest energy from RF signals of W for a duration of *αT* at the beginning of every frame, where

*T* is the period of one frame, and 0 < *α* < 1 is EH ratio. During the EH phase, the energy harvested by node l can be written as

$$E\_l = \sum\_{m=1}^{M} \eta \left( P\_{\mathbf{t}} \left| h\_{m,l} \right|^2 \right) \alpha T \tag{20}$$

where 0 < *η* < 1 is the energy conversion efficiency, *P*<sup>t</sup> is the transmit power of single W's antenna, and *hm*,*<sup>l</sup>* <sup>2</sup> is the channel power gain of link from the W's *<sup>m</sup>*-th antenna to node *<sup>l</sup>*, *m* ∈ {1, 2, ···, *M*}, *λm*,*<sup>l</sup>* = *λ*W,*l*, respectively. Subsequently, the phase (1 − *α*)*T* is equally divided into two time slots for the two-hop data transmission of SN. During the phase of data transmission, the RF-EH circuit of node *l* is turned off, and the transceivers of R, D1 and D2 are turned on. We assume that regardless of energy harvested by node *l* in each frame, they can be stored in its configured storage device (e.g., a supercapacitor) and can be immediately applied for subsequent data transmission. Furthermore, the storage device of node *l* lacks energy management function and has obvious leakage of electricity, which causes the node *l*'s residual energy at the end of each frame to be completely leaked [24,33].

**Figure 4.** Network system model.

In underlay paradigm, the instantaneous transmit power of node *l* must be strictly constrained such that the interference induced by node *l* remains below the threshold *P*I, i.e., the peak interference power that Q*<sup>n</sup>* can tolerate. It is also assumed that the circuit energy consumption of node *l* is ignored, i.e., the energy harvested by node *l* is only used for data transmission. Consequently, the transmit power of node *i* can be set as

$$P\_l = \min\left(\frac{2E\_l}{(1-\alpha)T}, \frac{P\_l}{Y\_l}\right) = \min\left(\rho P\_l Z\_l, \frac{P\_l}{Y\_l}\right) \tag{21}$$

where *ρ* = <sup>2</sup>*ηα* <sup>1</sup>−*<sup>α</sup>* ,*Zl* <sup>=</sup> *<sup>M</sup>* ∑ *m*=1 *hm*,*<sup>l</sup>* 2 , and *Yl* <sup>=</sup> max *<sup>n</sup>*=1,2,···,*<sup>N</sup> hl*,*<sup>n</sup>* 2 

Let *S*1(*t*) and *S*2(*t*) respectively denote the data signals transmitted to D1 and D2 by S at time *t*. S produces the NOMA signal *Ss*(*t*) by allocating distinct power coefficients *a*<sup>1</sup> and *a*<sup>2</sup> for D1 and D2 respectively, i.e., the power *P*s*a*1and *P*s*a*<sup>2</sup> are respectively allocated to D1 and D2 at S. Moreover, it is assumed that D1 is near node and D2 is far node. Due to the principle that lower power is allocated to near node D1 at S, *a*<sup>2</sup> > *a*<sup>1</sup> and *a*<sup>1</sup> + *a*<sup>2</sup> = 1. During the first time slot of every frame, R receives the signal *S*(*t*) transmitted by S. Since the far node D2's data *S*2(*t*) with large power are firstly decoded and the near node D1 's data *S*1(*t*) are secondly decoded using successfully perfect SIC [22] by R, the corresponding signal-to-noise ratio (SNR) to *S*2(*t*) and *S*1(*t*) at R can be respectively given by

$$\gamma\_{\mathbb{R},S\_2} = \frac{a\_2 P\_S X\_1}{a\_1 P\_S X\_1 + \sigma^2} \tag{22}$$

and

$$
\gamma\_{\mathbf{R}, S\_1} = \frac{a\_1 P\_S X\_1}{\sigma^2} \tag{23}
$$

where *<sup>σ</sup>*<sup>2</sup> is the variance of additive complex white Gaussian noise, and *<sup>X</sup>*<sup>1</sup> <sup>=</sup> <sup>|</sup>*h*S,*R*<sup>|</sup> 2 .

During the second time slot of every frame, R adopts superposition coding to make the new NOMA signal *SR*(*t*) and send it to D1 and D2. Let *b*<sup>1</sup> and *b*<sup>2</sup> be the power allocation coefficients at R for D1 and D2, respectively. Since D1 is closer to R than D2, higher power is allocated to D2's data in *SR*(*t*) at R as well, i.e., *b*<sup>2</sup> > *b*<sup>1</sup> and *b*<sup>1</sup> + *b*<sup>2</sup> = 1. Now the far node D2 decodes its data *S*2(*t*) from the signals forwarded by R, and the corresponding SNR for *S*2(*t*) at D2 can be expressed as

$$\gamma\_{\rm D\_2,S\_2} = \frac{b\_2 P\_R X\_2}{b\_1 P\_R X\_2 + \sigma^2} \tag{24}$$

where *X*<sup>2</sup> = *hR*,*D*<sup>2</sup> 2

At the near node D1, the far node D2's data *S*2(*t*) with larger power are first decoded and later D1's data *S*1(*t*) are decoded by achieving successfully perfect SIC [22]. Hence, the corresponding SINRs for *S*2(*t*) and *S*1(*t*) at D1 can be respectively given as

$$\gamma\_{D\_1, S\_2} = \frac{b\_2 P\_R X\_3}{b\_1 P\_R X\_3 + \sigma^2} \tag{25}$$

and

$$
\gamma\_{D\_1, S\_1} = \frac{b\_1 P\_R X\_3}{\sigma^2} \tag{26}
$$

where *X*<sup>3</sup> = *hR*,*D*<sup>1</sup> 2

#### *3.2. Outage Probability Analysis of SN*

In this section, the outage probability experienced by D1 and D2 will be studied. Let *R*<sup>1</sup> and *R*<sup>2</sup> respectively denote the target rates for D1 and D2. Since two consecutive time slots are needed to realize the communication from S to D1 and D2 in every frame, the achieved rates are halved. Hence, the target SINRs for successful decoding of *S*1(*t*) and *S*2(*t*) are respectively expressed as *γ*<sup>1</sup> = 2 <sup>2</sup>*R*<sup>1</sup> <sup>1</sup>−*<sup>α</sup>* <sup>−</sup> 1 and *<sup>γ</sup>*<sup>2</sup> <sup>=</sup> <sup>2</sup> <sup>2</sup>*R*<sup>2</sup> <sup>1</sup>−*<sup>α</sup>* <sup>−</sup> 1.

#### 3.2.1. OP for D1

To achieve the reliable transmission of *S*1(*t*) from S to D1, both *S*1(*t*) and *S*2(*t*) should be successfully decoded by R and D1. Accordingly, the OP for D1 is defined as

$$P\_{\text{out},1} = 1 - \Pr\{\gamma\_{\mathbb{R},\mathbb{S}\_2} \ge \gamma\_{\mathbb{2}} \gamma\_{\mathbb{R},\mathbb{S}\_1} \ge \gamma\_{1\prime} \gamma\_{\mathbb{D}\_1,\mathbb{S}\_2} \ge \gamma\_{\mathbb{2}} \gamma\_{\mathbb{D}\_1,\mathbb{S}\_1} \ge \gamma\_1\} \tag{27}$$

**Proposition 1.** *Assuming* 0 < *a*<sup>1</sup> < <sup>1</sup> <sup>1</sup>+*γ*<sup>2</sup> *and* <sup>0</sup> <sup>&</sup>lt; *<sup>b</sup>*<sup>1</sup> <sup>&</sup>lt; <sup>1</sup> 1+*γ*<sup>2</sup> *, OP expression for* D1 *is calculated as*

$$\begin{array}{l} P\_{out,1} = 1 - \left(\sum\_{n=0}^{N} \binom{N}{n} (-1)^n \frac{2\lambda\_{S,R}}{\Gamma(M)} \mathbb{C}\_1 \stackrel{M}{\mathbb{Z}} \mathbb{C}\_2 \stackrel{M}{\mathbb{Z}}^{-1} K\_M \left(2\sqrt{\mathbb{C}\_1 \mathbb{C}\_2} \right) \right) \\ \times \left(\sum\_{n=0}^{N} \binom{N}{n} (-1)^n \frac{2\lambda\_{R,D\_1}}{\Gamma(M)} \mathbb{C}\_3 \stackrel{M}{\mathbb{Z}} \mathbb{C}\_4 \stackrel{M}{\mathbb{Z}}^{-1} K\_M \left(2\sqrt{\mathbb{C}\_3 \mathbb{C}\_4} \right) \right) \end{array} \tag{28}$$

*where* Γ(·) *and Kϑ*(·) *are respectively gamma function and the modified Bessel function of second kind with order ϑ, and C*1*, C*2*, C*3*, and C*<sup>4</sup> *are defined in Appendix A.*

**Proof.** See Appendix A.

3.2.2. OP for D2

For completing successful transmission of *S*2(*t*) from S to D2, at first both *S*1(*t*) and *S*2(*t*) need to be reliably decoded by R; subsequently, *S*2(*t*) in *SR*(*t*) transmitted by R should be successfully decoded at D2. Accordingly, OP for D2 is determined by

$$P\_{\text{out},2} = 1 - \Pr\{\gamma\_{\mathbb{R},\mathbb{S}\_2} \ge \gamma\_{2\prime}\gamma\_{\mathbb{R},\mathbb{S}\_1} \ge \gamma\_{1\prime}\gamma\_{D\_2,\mathbb{S}\_2} \ge \gamma\_{2\prime}\}\tag{29}$$

**Proposition 2.** *Assuming* 0 < *a*<sup>1</sup> < <sup>1</sup> <sup>1</sup>+*γ*<sup>2</sup> *and* <sup>0</sup> <sup>&</sup>lt; *<sup>b</sup>*<sup>1</sup> <sup>&</sup>lt; <sup>1</sup> 1+*γ*<sup>2</sup> *, OP expression for* D2 *is given by*

$$\begin{array}{l} P\_{out,2} = 1 - \left(\sum\_{n=0}^{N} \binom{N}{n} (-1)^n \frac{2\lambda\_{\overline{5},\mathsf{R}}}{\Gamma(M)} \mathbb{C}\_1 \stackrel{\scriptstyle \mathsf{M}}{\ } \mathbb{C}\_2 \stackrel{\scriptstyle \mathsf{M}}{\ } -1 \,\mathsf{K}\_M (2\sqrt{\mathsf{C}\_1 \mathsf{C}\_2})\right) \\ \times \left(\sum\_{n=0}^{N} \binom{N}{n} (-1)^n \frac{2\lambda\_{\overline{5},\mathsf{D}\_2}}{\Gamma(M)} \mathbb{C}\_5 \stackrel{\scriptstyle \mathsf{M}}{\ } \mathbb{C}\_6 \stackrel{\scriptstyle \mathsf{M}}{\ } -1 \,\mathsf{K}\_M \left(2\sqrt{\mathsf{C}\_5 \mathsf{C}\_6}\right)\right) \end{array} \tag{30}$$

*where C*<sup>5</sup> *and C*<sup>6</sup> *are defined in Appendix B.*

#### **Proof.** See Appendix B.

#### *3.3. Analysis and Optimization of SN's Throughput*

In this section, to achieve the maximization of SN's throughput, the jointly optimal EH ratio (*α*∗) and power allocation factors at S *a*∗ <sup>1</sup>, *a*<sup>∗</sup> 2 and R *b*∗ <sup>1</sup> , *b*<sup>∗</sup> 2 have to be determined. Considering delay-limited transmission mode, where the fixed target rates at D1 and D2 are respectively *R*<sup>1</sup> and *R*2, the throughput of SN is expressed as [14]

$$\pi = R\_1(1 - P\_{out,1}) + R\_2(1 - P\_{out,2})\tag{31}$$

where Pout1 is the outage probability of SN1 and Pout2 is the outage probability of SN2.

For the purpose of maximize the throughput of SN, the resource allocation problem with respect to EH ratio *α* and power allocation factors *a*1,*a*2,*b*<sup>1</sup> and *b*<sup>2</sup> can be formulated as

$$\begin{array}{c} \max\_{\left(a, a\_1, a\_2, b\_1, b\_2\right) = \left(a^\*, a\_1^\*, a\_2^\*, b\_1^\*, b\_2^\*\right)} \tau \\ \text{s.t.} \\ a\_1 + a\_2 = 1, b\_1 + b\_2 = 1 \\ 0 < a < 1, a\_1 < a\_2, b\_1 < b\_2 \end{array} \tag{32}$$

To our best knowledge, the convexity of the objective function is hard to be estimated due to the complex expression of SN's throughput. Consequently, the throughput optimization problem is difficult to be settled in general methods, such as convex optimization. For resolving the SN's throughput optimization problem, a meta-heuristic algorithm based on HSIBA is applied in our paper.

#### **4. Optimization Results and Performance Analysis**

In this section, to maximize the throughput of SN, we apply the HSIBA to optimize the EH ratio and power allocation factors of the EH-NOMA-CRN. Without any loss of generality, the network system environment follows by the system model described in Section 3. The nodes S, R, W, D1, D2 and Q*<sup>n</sup>* are located at (−1, 0), (0, 0), (0, −1), (0.5, 0.5), (1, 0) and (1, 1), respectively. Unless otherwise stated, the key simulation parameters are listed in Table 5.

In Figures 5 and 6, we plot unoptimized SN's throughput and the optimized SN's throughput by the proposed HSIBA as a function of *P*t(*P*I) with different *N*(*R*<sup>1</sup> and *R*2), respectively. Moreover, we also compare the HSIBA with the exhaustive searching algorithm. The results show that: (1) the throughput increases with the increase of *P*t(*P*I) and tends to be saturated gradually. That is because SN nodes' transmission power increases with the increase of *P*t(*P*I), which leads to the throughput raise. Nevertheless, due to the limitation of *P*I(*P*t), SN nodes' transmission power, OP and throughput of the SN all tend to

saturation; and (2) In Figure 5, for given *P*t, the probability that interference channel obtains higher power gain increases with increase of *N*, which leads to throughput increasing; and (3) In Figure 6, for given *P*I, different values of *R*<sup>1</sup> and *R*<sup>2</sup> correspond to different SN's throughput; and (4) From Figures 5 and 6, it is distinctly observed that the optimized throughput values are significantly higher than those unoptimized values either for given *P*<sup>t</sup> or *P*I; and (5) It also can be clearly observed that the optimal throughput curves by proposed HSIBA extremely approximate the optimal curves by the greedy search algorithm. However, the latter has to search the optimal solution in a three-dimensional continuous space generated by *α*, *a*<sup>1</sup> and *b*1, which results in a much higher computational complexity than the former.

**Table 5.** Simulation parameters and values.


**Figure 5.** Throughput versus *P*t for various *N*.

Figure 7 shows the relationship between optimal SN's throughput value search of five different optimization algorithms and the number of iteration *t*. It can be observed that: (1) The four kinds of optimization algorithms except AOA can all effectively optimize the proposed EH-NOMA-CRN and successively search the optimal SN's throughput value; (2) For HSIBA, HPO, AOVA, ELBA and standard BA, they need 8, 21, 31, 33 and 42 iterations

to converge, respectively. Thus it can be seen that the HSIBA has better convergence speed and higher stability.

**Figure 6.** Throughput versus *P*<sup>I</sup> for various *R*<sup>1</sup> and *R*2.

**Figure 7.** Optimal throughput search versus *t*.

#### **5. Conclusions**

In this paper, we propose a novel HSIBA with adaptive strategy and multiple flight modes selection, and respectively derive the closed-form solutions of SN's OP and throughput. On the base of the throughput analysis, we further formulate the throughput maximation problem with regard to EH ratio and power allocation factors. In addition, due to the complexity of throughput solution, we employ the HSIBA to jointly optimize EH ratio and power allocation factors to maximize the SN's throughput. Through comparing with other advanced meta-heuristic algorithms, we confirm the correctness and effectiveness of proposed HSIBA by a large number of benchmark functions' test and numerical results. Based on the network parameters setting and grasp of channel state information, the proposed

HSIBA can achieve joint optimization to EH-NOMA-CRN based on three parameters with faster convergence speed and higher stability.

**Author Contributions:** Conceptualization: Y.L. (Yi Luo) and C.W.; experimentation and data analysis: N.H. and L.M.; writing—original draft preparation: Y.L. (Yi Luo), N.H. and C.W.; writing—review and editing: C.W., Y.L. (Yi Leng), N.H. and J.T.; funding acquisition: Y.L. (Yi Leng), C.W. and N.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by Airforce of People's Liberation Army, Government of China under grant No. PGGC-2021-003, the Hunan Provincial Science and Technology Project Foundation under Grant 2018TP1018, and Hunan Normal University Undergraduates Innovative Experiment Project and Entrepreneurship Program (Grant No. 2022180 and No. 2022183).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

In this Appendix A, we offer a proof of Proposition 1. To simplify the derivation process, we firstly provide the probability density function (PDF) and cumulative distribution function (CDF) of *Yl*, *Zl*, and *Ul* = *Pl* B *σ*<sup>2</sup> .

Note that *Yl* is the maximum of *N* independent exponential random variables (RVs) with the same expectation value 1 B *λ*W,*l*, and *Zl* is the sum of M independent exponential RVs with the same expectation value 1 B *λl*,Q . Therefore, the CDF and PDF of *Yl* and *Zl* can be respectively computed as

$$F\_{Y\_l}(y\_l) = \left(1 - \mathbf{e}^{-\lambda\_{l,Q}y\_l}\right)^N = \sum\_{n=0}^N \binom{N}{n} (-1)^n \mathbf{e}^{-n\lambda\_{l,Q}y\_l} \tag{A1}$$

$$f\_{\mathbf{Y}\_l}(y\_l) = N\lambda\_{l,\mathbf{Q}} \sum\_{n=0}^{N-1} \binom{N-1}{n} (-1)^n \mathbf{e}^{-(n+1)\lambda\_{l,\mathbf{Q}}y\_l} \tag{A2}$$

$$F\_{\mathbb{Z}\_l}(z\_l) = 1 - \frac{\Gamma(M\_\prime z\_l \lambda\_{\mathbb{W},l})}{\Gamma(M)} \tag{A3}$$

and

$$f\_{Z\_l}(z\_l) = \frac{\lambda\_{W,l}^M z\_l^{M-1} \mathbf{e}^{-\lambda\_{W,l} z\_l}}{\Gamma(M)} \tag{A4}$$

where Γ(.) is the upper incomplete gamma function.

Subsequently, we can calculate the CDF of *Ul* as follows

$$\begin{split} F\_{lI\_l}(u\_l) &= 1 - \Pr\left\{ Z\_l > \frac{u\_l \varphi^2}{\rho P\_l} \right\} \Pr\left\{ Y\_l < \frac{P\_l}{u\_l \varphi^2} \right\} \\ &= 1 - \sum\_{n=0}^N \binom{N}{n} (-1)^n \frac{\Gamma\left( M, \frac{\lambda\_{W\_l} u\_l \varphi^2}{\rho P\_l} \right)}{\Gamma(M)} \mathbf{e}^{-\frac{n \lambda\_{l,Q} P\_l}{u\_l \varphi^2}} \end{split} \tag{A5}$$

Then, Equation (27) can be reformulated as

$$P\_{\rm out,1} = 1 - \underbrace{\Pr\left\{\frac{a\_2 \text{l} \mathbf{I}\_S \mathbf{X}\_1}{a\_1 \text{l} \mathbf{I}\_S \mathbf{X}\_1 + 1} \ge \gamma\_{2\prime} a\_1 \text{l} \mathbf{I}\_S \mathbf{X}\_1 \ge \gamma\_1 \right\}}\_{A\_{00}} \underbrace{\Pr\left\{\frac{b\_2 \text{l} \mathbf{I}\_R \mathbf{X}\_3}{b\_1 \text{l} \mathbf{I}\_R \mathbf{X}\_3 + 1} \ge \gamma\_{2\prime} b\_1 \text{l} \mathbf{I}\_R \mathbf{X}\_3 \ge \gamma\_1 \right\}}\_{A\_{01}} \tag{A6}$$

By conditioning *A*<sup>00</sup> in Equation (A6) on *X*<sup>1</sup> and taking the expected value of the results over the distribution of *X*1, we can obtain

$$\begin{split} A\_{00} &= \Pr\{lL\_{\mathbb{S}}X\_{\mathbb{I}} \ge H \mid X\_{\mathbb{I}} = x\_{\mathbb{I}}\} \\ &= 1 - \int\_{0}^{\infty} F\_{\text{LI}\_{\mathbb{S}}} \left( \frac{H}{x\_{\mathbb{I}}} \right) \lambda\_{\mathbb{S},\mathbb{R}} \mathbf{e}^{-\lambda\_{\mathbb{S},\mathbb{R}} x\_{\mathbb{I}}} \mathbf{d} x\_{\mathbb{I}} \end{split} \tag{A7}$$

where *H* = max *<sup>γ</sup>*<sup>2</sup> *a*2−*a*1*γ*<sup>2</sup> , *γ*1 *a*1 

Applying Equation (6.453) in [34], *A*<sup>00</sup> can be expressed as

$$A\_{00} = \sum\_{n=0}^{N} \binom{N}{n} (-1)^n \frac{2\lambda\_{\text{S},R}}{\Gamma(M)} \mathcal{C}\_1 \,^{\frac{M}{2}} \mathcal{C}\_2 \,^{\frac{M}{2}-1} K\_M \left(2\sqrt{\mathcal{C}\_1 \mathcal{C}\_2}\right) \tag{A8}$$

where *<sup>C</sup>*<sup>1</sup> <sup>=</sup> *<sup>λ</sup>W*,*SHσ*<sup>2</sup> *<sup>ρ</sup>Pt* and *<sup>C</sup>*<sup>2</sup> <sup>=</sup> *<sup>n</sup>λ*S,*QPI <sup>H</sup>σ*<sup>2</sup> + *<sup>λ</sup>*S,R. In the same manner, the term *A*<sup>01</sup> in Equation (A10) can be obtained as

$$A\_{01} = \sum\_{n=0}^{N} \binom{N}{n} (-1)^n \frac{2\lambda\_{\rm R,D\_1}}{\Gamma(M)} \mathcal{C}\_3^{\frac{M}{2}} \mathcal{C}\_4^{\frac{M}{2}-1} K\_M \left( 2\sqrt{\mathcal{C}\_3 \mathcal{C}\_4} \right) \tag{A9}$$

where *<sup>C</sup>*<sup>3</sup> <sup>=</sup> *<sup>λ</sup>*W,*RGσ*<sup>2</sup> *<sup>ρ</sup>Pt* , *<sup>C</sup>*<sup>4</sup> <sup>=</sup> *<sup>n</sup>λ*R,*QP*<sup>I</sup> *<sup>G</sup>σ*<sup>2</sup> <sup>+</sup> *<sup>λ</sup>R*,*D*<sup>1</sup> , and *<sup>G</sup>* <sup>=</sup> max *<sup>γ</sup>*<sup>2</sup> *b*2−*b*1*γ*<sup>2</sup> , *γ*1 *b*1 

The OP for D1 can be calculated by applying the results obtained from Equations (A8) and (A9).

#### **Appendix B**

In this Appendix B, we provide a proof of Proposition 2. Equation (29) can be rewritten as follows

$$P\_{\text{out},2} = 1 - \underbrace{\Pr\{\gamma\_{\text{R},S\_2} \ge \gamma\_{2\text{-}}\gamma\_{\text{R},S\_1} \ge \gamma\_1\}}\_{A\_{00}}\underbrace{\Pr\{\gamma\_{D\_2,S\_2} \ge \gamma\_2\}}\_{B\_{00}}\tag{A10}$$

The term *A*<sup>00</sup> of Equation (A10) has been calculated in Appendix A, and the term *B*<sup>00</sup> of Equation (A10) can be computed as

$$B\_{00} = 1 - \int\_0^\infty F\_{\rm II\_R} \left(\frac{V}{X\_2}\right) \lambda\_{\rm R,D\_2} \mathbf{e}^{-\lambda\_{\rm R,D\_2} x\_2} \mathbf{d}x\_2 \tag{A11}$$

where *V* = *<sup>γ</sup>*<sup>2</sup> *b*2−*b*1*γ*<sup>2</sup> . Applying Equation (6.453) in [34], *B*<sup>00</sup> can be presented as

$$B\_{00} = \sum\_{n=0}^{N} \binom{N}{n} (-1)^n \frac{2\lambda\_{\rm R,D\_2}}{\Gamma(M)} \mathbb{C}\_5 \overset{\Delta}{\leftarrow} \mathbb{C}\_6 \overset{\Delta}{\leftarrow} {\mathop{\rm L}}\_M \Big( 2\sqrt{\mathbb{C}\_5 \mathbb{C}\_6} \Big) \tag{A12}$$

where *<sup>C</sup>*<sup>5</sup> <sup>=</sup> *<sup>λ</sup>W*,*RVσ*<sup>2</sup> *<sup>ρ</sup>Pt* and *<sup>C</sup>*<sup>6</sup> <sup>=</sup> *<sup>n</sup>λ*R,*QPI <sup>V</sup>σ*<sup>2</sup> + *<sup>λ</sup>*R,*D*<sup>2</sup> .

The OP for D2 can be calculated by applying the results obtained from Equations (A8) and (A12).

#### **References**

