Average SER analysis of two‐hop WP DF relay system under κ−μ shadowed fading

Correspondence Kalpana Dhaka, Department of Electronics and Electrical Engineering, Indian Institute of Technology Guwahati, Guwahati 781039, India. Email: kalpana.dhaka@iitg.ac.in Abstract This paper investigates the performance of a two-hop decode-and-forward relaying system when the source node is energy-constrained and it harnesses energy using the radio frequency signal radiated by the relay node. Data transmission at the energy-constrained source node is enabled by the energy harvested at the node. The analysis is presented for two modulation schemes: coherent M -ary phase-shift keying and orthogonal non-coherent M -ary frequency-shift keying. Analytical and asymptotic expressions for the average symbol error rate are derived for the considered modulation schemes under κ−μ shadowed fading. The results are plotted to investigate the system’s performance with variation in modulation order, fading parameters, and relay location.

to the destination. The other works available in the literature consider a general scenario of a three-node relay system [3,9,[20][21][22][23][24]. The basic idea is to enable end-to-end communication when low power devices are energy-constrained. One specific application of such three-node relay systems can be in medical implants, where a low power source device (such as pacemeter, glucose sensor etc.) is implanted in a patient's body. The implanted device sends vital information to a doctor (destination) with the help of a nearby located hub. The hub acts as a relay between the implanted device and the device at the doctor's end. Also, it provides wireless energy to the body implant. Furthermore, these three-node WP relay systems and their variants having multiple nodes can find applications in systems including Internet of Things (IoT) entities, wireless body area network (WBAN), device-to-device (D2D) communications, machine-to-machine communications, and uplink cellular transmission.
In the literature, a variety of fading models are considered to characterize the wireless medium. Rayleigh fading model is a usual consideration for scattering environment and Rician fading is more practical in case of communication with Lineof-Sight (LoS) component. The generalized Nakagami-m channel model reduces to Rayleigh for fading parameter m = 1 and approximates to Rician fading for m > 1. Rayleigh, Rician, and Nakagami-m distributions model homogeneous fading environments [25]. In practice, non-homogeneous generalized fading models, including, − fading, − fading, and − shadowed fading provide a better fit to the experimental data [25][26][27]. Using experimental setups in [26] and [27], it is demonstrated that − and − shadowed fading models provide a more accurate and reliable fit to WBAN and D2D communication environments, respectively. The experimental results also suggest that the generalized fading models when compared to the basic homogeneous fading models are more suitable for short-range communications. Furthermore, using the state-ofart technologies the applicability of WP systems is limited to short-range communications. Hence, it is more appropriate to analyse WP systems' behaviour under the consideration of the generalized fading models. The listed generalized fading models reduce to the basic ones on substituting the corresponding fading parameter as shown in [28].
Generalized − , − , and − shadowed fading models are considered in [29][30][31][32][33] to examine the performance of WP relay systems with EH at the relay node(s). To the best of our knowledge, no work on WP relay systems with EH at source node under consideration of the generalized fading models is reported in the literature. This paper considers a three-node WP relay system with a source node wirelessly powered by the relay node and links are affected by − shadowed fading. The endto-end average SERs are deduced considering two modulation schemes, namely, M -ary phase-shift keying (M -PSK) modulation with coherent detection and orthogonal M -ary frequencyshift keying (M -FSK) modulation with non-coherent detection.
In the literature, probability density function (PDF) and moment generating function (MGF)-based approaches are mainly considered for averaging the instantaneous SER expressions. MGF is used to get moments of the corresponding random variable. It can also be applied to simplify the complicated analysis. The MGF-based approach is used for M -PSK modulated data, whereas the PDF-based approach is considered for the orthogonal M -FSK modulation scheme by this study.

Contributions
The work presented in this paper is an extension of our previous work [3], where the system is considered to be under Nakgamim fading. In the present paper, to provide a more accurate and reliable analysis, we consider that the system is under − shadowed fading. The main contributions of this paper include (i) consideration of a three-node WP relay system with EH at the source node under generalized − fading, (ii) derivation for the corresponding average SERs when data is M -PSK and orthogonal M -FSK modulated, (iii) high signal-to-noise ratio (SNR) approximation of the average SERs to get simplified expressions, and (iv) examining the performance through numerical results with variation in different system parameters. In addition, another major contribution of this paper is PDF and MGF expressions for the product of two independent and non-identically distributed − shadowed random variables. The derived average SER expressions are novel and can be useful in investigating the system's behaviour for a wide range of channel conditions by varying the fading parameters. Furthermore, this work enriches the theoretical aspects by filling the gaps in this domain of research.

Organization
This paper is organized into five sections. Section 2 provides a description of the system considered. The communicationtheoretic approach is used to deduce performance of the system for M -ary modulation schemes in Section 3. The average SER expressions obtained in Section 3 are plotted in Section 4 for varying parameters affecting performance of the system. Concluding remarks based on the numerical results plotted in Section 4 are stated in Section 5.

SYSTEM MODEL
We consider a WP two-hop relay system containing an energyconstrained source node S , a relay node R, and a destination node D. All three nodes operate in half-duplex mode and feature a single antenna at each node. In odd time slots, node S transmits the data to node R, which broadcasts processed data in the following even slot. We consider node R employs fixed DF protocol to process received data, that is, information detected at the node is forwarded over the channel. In even slots, the data transmitted by the relay R is received at nodes S and D, which utilize the signals for EH and data decoding, respectively. The harvested energy is consumed by the source S for transmitting the data in subsequent (odd) slot. This cycle continues until all the data at node S are communicated to node D. We consider that energy harvested at source S is processed using the harvest-use approach, that is the energy harvested in a slot is completely consumed and it cannot be stored beyond the current slot duration [1]. In addition, we consider that the circuitry for transmitting and receiving information, and processing energy at a node consumes negligible power [34]. The average SER of the system is analysed for M -PSK and orthogonal M -FSK modulated symbols. At the receiving ends, coherent detection and non-coherent detection are performed for M -PSK and orthogonal M -FSK modulation, respectively. Each symbol has the support of [0, T s ) and is transmitted with equal a priori probability. We define the constellation  = {S 1 , … , S M }, where S 1 , … , S M are the constellation points for the M symbols [3]. The inter-node links are independent and − shadowed faded. The baseband equivalent of received RF signals at node R in the odd slot and at nodes S and D in even slot are given by and respectively. x is symbol transmitted in the odd slot with power P S andx is the detected symbol at node R which is then forwarded in even slot with power P R , where x ∈  andx ∈ . Symbols x andx have unit energy. h is the complex channel gain of ∈ {SR, RS, RD} link, (d ) is path loss of link occurred at distance d meters with path loss exponent (typically, ∈ (2, 6) for wireless channels). The noise components n R and n j are uncorrelated and zero-mean complex Gaussian with power spectral density (PSD) N 0 watt per Hertz (W/Hz). Here, energy harvested at node S is directly related to the transmission power P R at node R in even time slot of the previous transmission. Therefore at the start of communication, node R can convey an arbitrary RF signal to harvest energy at node S . In case the noise component is negligible, the energy harvested at node S can be quantified using (2) and it can be approximated as [5] where is efficiency of the energy conversion unit. Now E S = P S T s , thus using (3) the transmission power at node S is given by

PDF and MGF of instantaneous SNRs
In this subsection, we derive the PDF and MGF of the instantaneous received SNRs SR and RD for source-relay and relaydestination links, respectively. The received SNRs are the ratios of the signal energy and noise energy at the relay and destination nodes. Using (1), (2), and (4), the instantaneous received SNRs are given as . The corresponding average SNRs arē is mean power of the channel gain h . In general, the uplink and downlink channels in time-division duplexing (TDD) are assumed to be reciprocal. However, owing to non-reciprocity problems, the channel conditions can vary for different slots [35]. We consider the channel conditions of RS and SR links in consecutive time slots can vary independently and, hence, the channel gains h s and the path loss exponents RS and SR are regarded as different. The PDF of power variate for − shadowed faded link gains are given by [36,Equation (4)] where , , m are fading parameters of link , = E[|h | 2 ] is the mean power of |h |, and 1 F 1 (⋅) represents confluent hypergeometric function (38). Parameter is ratio of the total power of the LOS components and the total power of non-LOS components, parameter corresponds to the number of clusters, and parameter m captures the shadowing effect of the fading environment. Statistical characterization of − shadowed fading is elaborated in [36].
Deriving the average SER using the PDF in (6) requires dealing with the integration of integrands including confluent hypergeometric function. This makes it difficult to obtain closed-form expressions of average SERs. To simplify the derivation, we use the series form representation of the PDF given in [37]. In [37], the series form representation of confluent hypergeometric function (38) is used to get the PDF as where c and p are given by c = (1 + ) and respectively, with n = 1 for = RS , n = 2 for = SR, and n = 3 for = RD. (x) l = Γ(x + l )∕Γ(l ) represents Pochhammer symbol. The infinite series PDF given in (7) can be approximated to a series containing finite number of terms, provided the approximation error is below a desired value. For instance, the approximated PDF with N n terms is given by and the corresponding approximation error can be expressed as It is shown in Appendix A.2 that the series given in (8) converges with an increase in N n . Equation (A.17) can be utilized to numerically evaluate the desired number of summation terms. Next, the SNR SR in (5a) is a scaled version of the product of two independent random variables |h RS | 2 and |h SR | 2 . Its CDF (cumulative distribution function) can be obtained as In the above relation, we have substituted |h RS | 2 = t without actually deviating from the final result. Hence, the PDF of SR can be obtained by taking the first-order derivative of (10) with respect to SR as Now, on substituting (8) in (11) and interchanging the order of summations and integration, we have where z = l 2 − l 1 + SR − RS . Using (39), the PDF in (12) can be expressed as where a SR = RS SR (1 + RS )(1 + SR ) and K r (⋅) is the rthorder modified Bessel's function of the second kind. Here, a SR is obtained using the expression for c given in (7). The SNR RD is a scaled version of |h RD | 2 . Hence, using the property of the random variables, the PDF of RD is obtained using (5b) and (8) as where a RD = RD (1 + RD ). Next, the MGF of the instantaneous SNR is given by [5,Equation (17)] Hence, the MGF of SR can be obtained by substituting (13) in (15). On interchanging the order of summations and integration the MGF is given by where z ′ = (l 1 + RS + l 2 + SR )∕2 and ′ = (l 2 + SR − l 1 − RS ). Substituting SR = t 2 and using (A.3), the MGF in (16) is simplified as where W s,t (⋅) is Whittaker's function. Similarly, on putting (14) in (15) and using (A.4), the MGF of RD is written as

PERFORMANCE ANALYSIS
The end-to-end average SER of the system can be given by [3, Equation (19)] where P e, is the average SER and P for = {1, … , M − 1} are paired error probability when a symbol is mapped to one of the remaining M − 1 symbols [38] for link.

3.1
Analytical end-to-end average SER

M-PSK modulation
Using [3,Equations (20) and (21)], the error terms P e, and P in (19) can be represented in integral form as , and M (s) is the MGF of . On putting (20) and (21) in (19) and using MGFs (17) and (18), the end-to-end average SER for M -PSK modulation scheme under − shadowed fading can be evaluated. Numerical evaluation of integrals in the expression is done using MAT-LAB.

Orthogonal M-FSK modulation
The average SER in (19) reduces to [3, Equation (27)] for orthogonal M -FSK modulated data, that is Equation (22) requires the computation of P e, , which can be done by averaging the conditional SER [3, Equation (28)]. The error term P e,SR can be analysed by taking expectation of the conditional SER (A.5) with respect to the PDF in (17). The simplified form of expression for P e,SR can be obtained using (A.3) and following the steps used in deriving (17), we get On averaging (A.5) using the PDF of RD (18) and (A.4), the error term P e,RD in (22) is given by On substituting (23) and (24) in (22), the end-to-end average SER for orthogonal M -FSK modulation scheme under − shadowed fading can be obtained.
Note that the deduced exact average SER expressions in this subsection involve either computation of Whittaker's function or its integration, which is complex and time-consuming. In order to simplify expressions for the average SER, high SNR approximations are obtained.

Asymptotic end-to-end average SER
In order to simplify the expressions of average SER at high SNRs asymptotic approximations are obtained. The average SER in (19) can be approximated as where P ∞ e , P ∞ e,SR , and P ∞ e,RD are high SNR approximations of P e , P e,SR , and P e,RD , respectively. The approximations for the two modulation schemes are deduced in this subsection.

M-PSK modulation
To obtain P ∞ e in (25), an approximate expression of the PDF f SR ( SR ) is used to get P ∞ e,SR . The PDF in (13) has a series containing K r (⋅), which can be approximated using (43) and respectively, where  1 and  2 are sets of values for l 1 and l 2 when r = 0, sets  1 ′ and  2 ′ correspond to the values of l 1 and l 2 for |r| > 0, a = ( RS + l 1 + SR + l 2 )∕2, and b = ( RS + l 1 + SR + l 2 − | SR + l 2 − RS − l 1 |)∕2. The PDFs in (26) and (27) can be unified as where e is Euler's number, r = a , 1 = 0, and 2 = 1 for case (a), and r = b , 1 = 1, and 2 = 0 for case (b). The representation { n ,  ′ n } suggests for the selection of cases (a) and (b) depending on the values of l 1 and l 2 .
The high SNR approximation of the average SER can be found on replacing (30) and (32) in (25). Thus where We observe that the integrals ,  , and  in the asymptotic average SER (33) can be simplified in terms of beta and digamma functions, which are inbuilt functions in mathematical software and can also be represented efficiently in series form.
There is no need to compute time-consuming numerical integration as in (20) and (21) for the exact analysis. This implies the asymptotic analysis is considerably simplified in terms of computational efficiency.

Orthogonal M-FSK modulation
For orthogonal M -FSK modulated data, P ∞ e,SR is derived by taking expectation of the conditional SER (A.4) using the PDF in (28) and (A.6) as ) .
For v̄R D ∕((v + 1) RD ) ≫ 1, (24) is reduced to On putting (34) and (35) in (25), the high SNR approximation of the end-to-end average SER is given by where (l 1 , l 2 ; M ) We further observe that the asymptotic average SER (36) does not contain Whittaker's function as in (23). Hence, the average SER for orthogonal M -FSK modulated data can be computed efficiently using (36).

Diversity order (DO)
The DO for any communication system is described as the order of improvement in the performance with an increase in SNR. Generally, it is quantified as the negative slope of the average SER curve at high SNR on the log-log axis. The DO of the considered system can be analysed on substituting the approximate expressions of the average SER is given by (33) and (36) in DO = − lim →∞ ln P ∞ e ∕ ln̄ [5]. Since the expressions for asymptotic average SERs is a sum of two error terms P ∞ e,SR and P ∞ e,RD , the DO is dominated by the one with the smallest order of̄(= P R T s ∕N 0 ) in the denominator. Hence, we have RD ∕ ln̄. Furthermore, the error factors P ∞ e,SR and P ∞ e,RD are in series form. Only the first term (corresponding to l n = 0) of these two has the smallest power of in the denominator. Hence, only these first two terms would add to the resultant DO. Considering ln(̄) is sufficiently larger than the other factors in P ∞ e,SR and after doing some algebraic manipulations, we get wherēr = {( SR + RS )∕2, = 1 for case (a) and = 0 for case (b). At significantly high SNRs, the system possesses DO = min{̄r , RD }.

NUMERICAL RESULTS
We present plots for the average SER of the system described in Section 2. The suitable ranges of various parameters considered for the analysis. It is evident that − shadowed fading provides flexibility to analyse the system's performance for a wide range of channel conditions attributed by the fading parameters , , and m . In this section, we choose specific values of these fading parameters to visualize their impact on the performance. However [27,Tabel I], and other experimental results could be of more help while considering a specific scenario. Furthermore, the transmission power of low-power devices can be in the range of micro Watts ( W) and high as a few Watts. The noise PSD is usually considered in the range of −70 to -120 dBm and transmission bandwidth falls in the range of MHz. The path loss exponents for wireless communication can be considered in the range of 2-6. The results are obtained considering power transmitted at node R is P R ∈ [−10, 30] dBm (P R ∈ [100 W, 1 W]), noise PSD N 0 = −70 dBm (N 0 = 10 −10 W/Hz), and symbol bandwidth B (in MHz) (commonly used in wireless systems). Assuming B = 1 MHz (symbol duration T s (= 1∕B) is in range of microseconds), the noise power is 10 −4 W. Furthermore, we consider path loss exponent of each link is = 3 and distance d = 1 m, unless otherwise stated. The variance of channel gains is taken as unity, that is, = 1. The average SER expression for M -PSK modulation scheme is deduced as (19), where the terms are defined in (20) and (21) with MGFs in (17) and (18). In the case of orthogonal M -FSK modulation scheme, the average SER is obtained on putting (23) and (24) in (22). The factors P e,SR and P e,RD in the average SER expression are represented in a series form which converges faster with an increase in the number of terms. It is observed through numerical evaluation in MATLAB that for an increase in ∕m , more number of terms are required in (8) to get close to the actual PDF. Typically, for the fading parameters considered in this paper, one would get the desired PDF using (8) when the number of terms is around 100. Furthermore, the average SER expressions are obtained on evaluating the integrals involving the PDF (8). It is observed that the required number of summation terms reduces after each integration. Let P 1 = (N 1 N 2 ) and N 3 terms are required for the computation of P e,SR and P e,RD , respectively. The required number of terms in the series varies with change in the value of fading parameters , , and m . Also, for M -PSK and orthogonal M -FSK modulations different values of P 1 and N 3 are required to attain the same level of accuracy. Without going into details of the computation for the required number of terms depending on various parameters, analytical results are obtained which ensure the accuracy up to the fifth decimal point with respect to the simulation results. For M -PSK and orthogonal M -FSK, (17) and (23), respectively, are twofold summations. The number of terms P 1 is equally divided for the two summations, that is, the upper limit of each summation is N 1 = N 2 = P 1 ∕2 = 10. Using (26) and (27), asymptotic average SERs are also evaluated and plotted to check accuracy of the approximation. In Figures 1 and 2, plots of the average SER versus relay transmission power P R are presented for the two modulation schemes. The plots are shown for different modulation order and the fading parameters. In Figure 1, the results are plotted for parameters = 1, = 1, and m = 1 and the parameters = 1, = 2, and m = 3 are considered in Figure 2, ∈ {SR, RS, RD}. The analytical results of the two modulation schemes closely match with the simulation results for all modulation orders. Hence, our analysis of the system is validated. Moreover, it is found that the asymptotic results at high P R (directly related to SNR) provide a good approximation. as expected M -FSK outperforms M -PSK for higher modulation order. This observation is more distinguishable for the larger value of the fading parameters as seen in Figure 2.
Note that the terms P ∞ e,SR and P ∞ e,RD in the asymptotic average SER expressions should be between zero and unity (∈ [0, 1]). However, due to the approximation involved in deriving these expressions, they may violate the axiom by taking values that are negative or more than unity at low SNRs and/or with variation in different parameters. Specifically, at low SNRs, the term P ∞ e,SR can be negative, whereas a positive value greater than unity can be observed for P ∞ e,RD . The term P ∞ e,SR results in negative values when the factors − 2 ( ( r ) +  ) and − 2 ( ( r ) in (30) and (34), respectively dominate over the other factors. Moreover, at low SNRs, the factors ln(g 0̄SR ∕(4ea SR )) and ln(̄S R ∕(4e( + 1)a SR )) in (30) and (34), respectively can be negative, which further adds in overall negative value of P ∞ e,SR . On the other side, the term P ∞ e,RD is always positive but it can be greater than unity depending on the values of SNR and other parameters due to the factors (a RD ∕(g 0̄RD )) RD +l 3 and (( + 1)a RD ∕(̄R D )) RD +l 3 in (32) and (35), respectively. Thus, the asymptotic average SER expression obtained on summing the terms P ∞ e,SR and P ∞ e,RD for the considered modulation schemes at low SNRs may be less than zero or greater than unity, depending on which one of the two terms dominates over the other. The same can be seen in Figures 1 and  2. In Figure 1, the resultant summation is found to be greater than unity at low SNRs, this implies that P ∞ e,RD dominates over P ∞ e,SR . As the average SER becomes greater than unity, it cannot be shown in the figure having y-axis range between zero and unity. The irregular shape of the curve at medium SNRs is due to the transition from low to high SNR region. Similarly, in Figure 2, as P ∞ e,SR dominates over P ∞ e,RD at low SNRs, the resultant average SER is negative. The negative values cannot be shown in the figure as it corresponds to positive values only.
In Figure 3, we plot average SER versus source-to-relay distance of 4-PSK for variation in , ∈ {SR, RS, RD}. We con-  Figure 3(a), we consider fading parameters RD = 2, = 2, m = 2 and plot average SER for varying e ∈ {1, 2, 3}. We observe that with an increase in e , the relay location for optimal performance shifts closer to the destination for fixed RD . In Figure 3(b), plots are shown for RD ∈ {1, 2, 3} when other fading parameters are constant, that is, e = 2, = 2, m = 2. The relay location for optimal performance moves closer to the source with gain in RD . Thus we conclude that the optimal relay location is closer to the node with poor quality in order to facilitate reliable communication between the nodes. Alternatively, we can also say that (i) for d SR < 1 (relay is located closer to the source than destination) increment in RD has a better impact on performance than with gain in e and (ii) for d SR > 1 (relay is located closer to destination than source) the effect of e on performance is more than that due to RD . These observations hold for all modulation orders of the considered modulation schemes for the cases (a) and (b).
In Figure 3, we compared the performance with variation in e and RD . Next, the effects of change in other fading parameters for the SR link are compared in Figure 4. The assumptions regarding node placement and fading parameters are the same as viewed in Figure 3. In Figure 4(a), plots showing variation in average SER with e and m e are shown for fixed , RD , and m RD . We observe that raise in e has more advantages in terms of performance rather than those due to increment in m e . Similarly, in Figure 4(b), the performance is compared for varying e and e while considering constant values of m , RD , and m RD . The performance for increment in e dominates over that due to an increase in e . Moreover, on comparing Figures 4(a) and (b), it can be observed that performance gain is more for increment in m e than e . These results also hold for variation in fading parameters of RD link. The analysis is valid for different modulation orders of the considered schemes. Furthermore, a comparison of the considered WP three-node relay system can be done with [39] and [5] in terms of the optimal relay location. The work in [39] considers average SER analysis for a conventional three-node relay system. The term conventional implies that all the three nodes have sufficient amount of energy for the end-to-end communication. The work in [5] considers the scenario when the source node provides wireless power to the energy-constrained relay node. The optimal relay location in [39] lies between the source and the destination nodes, while it is close to the source node in [5]. Compared to these two cases, the optimal relay location with EH at the source node lies between the source and destination nodes with an inclination towards the node having poor link quality.
In Figure 5, the average SER is plotted against the relay transmission power P R for different values of the fading parameters to visualize their effect on the DO of the system. The negative slope of the plots is observed with variation in SR = RS = e and RD for = 1 and m = 1, ∈ {SR, RS, RD} to determine the DO. We observe that DO for curves 1, 2, and 3 is 1, while it is 2 for curves 4 and 5. These results validate the correctness of the derived analytical expression for DO in (37).

CONCLUSION
We have deduced the average SER of a DF relaying system with source node WP through RF signal broadcasted at the relay node. M -PSK and orthogonal M -FSK modulation schemes are considered. The system is assumed to be affected under − shadowed fading. The obtained average SER expressions are validated with simulation results and hence can be utilized to investigate the system's performance with variation in modulation order/scheme, link quality, and relay location. Through numerical results, we observed that (i) M -PSK performs better than orthogonal M -FSK for M ≤ 4 and the later outperforms for M ≥ 8 and (ii) the optimal relay location is in between the source and destination nodes with an inclination towards the node with poor link quality which increases on the decrement in its quality. Among fading parameters , , and m of each link, the increment in has the best impact on the performance improvement, while increment in has the least impact. Simplified expressions of the average SER for the considered modulation schemes are obtained at high SNRs using asymptotic approximations which are later used to obtain the DO of the system.   ) . (A.5) • An approximation of K r (t ) for t → 0 can be given using [44, Equations (9.6.6), (9.6.8) and (9.6.9)] as Γ(a + n) ) .

A.2 Convergence of the PDF in (8)
The approximate expression for the error in (9) can be rewritten as Using the relations for c and p ,l n given in (7), the approximation of error in (A.12) is expanded (after dropping the subscript and substituting |h | 2 = t ) as In (A.14), it is observed that convergence of N is mainly limited by the terms  1 (N ) and  2 (N ). It is shown below that these two terms converge with an increase in N . First, we proceed with the term  1 (N ), which is expanded using relation (a) n = Γ(n + a)∕Γ(a) and (A.10) as  1 (N ) = aN m Γ( + N + 1) ) z N +1 , (A.15) where a = Γ( )∕(Γ(m) 2 ) and b = N (N + 1). Ignoring the term o(|N | −2 ) and applying first-order differentiation on both sides of (A.15) with respect to N , then using (A.11) followed by rearranging the terms, we get ( + N + 1) i(i + + N + 1) − 1 ( + N + 1) )) .