On the Secure Performance of Intelligent Reflecting Surface-Assisted HARQ Systems

This paper analyzes the physical layer security performance of hybrid automatic repeat request (HARQ) systems with the assistance of an intelligent reflecting surface (IRS) and aims to reveal the primary factors that enhance PLS. First, closed-form expressions for the connection outage probability (COP) and secrecy outage probability (SOP) in HARQ with chase combining (HARQ-CC) are acquired using the generalized-K (KG) distribution. Then, these two critical metrics are derived while adopting HARQ with incremental redundancy (HARQ-IR), resorting to the mixture gamma (MG) distribution and the Mellin transform. Diversity and coding gain are also addressed through an asymptotic analysis of the COP and SOP. Finally, an evaluation of the numerical results demonstrates that a greater gain in the main channel and the wiretap channel can be produced by increasing the number of meta-surfaces rather than increasing the maximum transmission number, except for the higher signal-to-noise (SNR) region of HARQ-IR where the latter is preferred. This finding provides a significant guidance for the joint configuration of IRS and HARQ to achieve secure communication.


Introduction
Rapid developments in wireless communications have made intelligent reflecting surface (IRS), sometimes referred to as reconfigurable intelligent surface (RIS), an important candidate technology for 6G due to its flexible control and superior performance. IRS can be deployed in different locations with extremely low costs compared to base stations (BSs) or small cells. Smart radio environments (SREs) can also be achieved by adjusting the phase, frequency or other parameters of the reflected signal using a soft-defined architecture, which effectively enhances signal coverage and network transmission performance [1,2].
Significant research has been conducted to speed up the application of IRS, including studies on channel modeling, performance analysis and resource allocation. Generally, IRS channel modeling consists of two components; those are the transmitter-IRS channel and the IRS-receiver channel. A typical approach is to uniformly model the transmitter-IRSreceiver process as a cascade channel (i.e., the product of a two-part channel fading [3][4][5]). Basar et al. proposed a phase adjusting scheme aiming to improve the received signalto-noise ratio (SNR) [6]. The central limit theorem (CLT) was applied so as to assess the probability density function (PDF) for this SNR when the amount of meta-surfaces (N) was great, to further derive a symbol error rate (SER) under different modulation conditions. However, in some scenarios, it is not necessary to allocate multiple meta-surfaces to a single user, and thus CLT is not suitable when N is small. Yang et al. noted a cascade channel can be accurately modeled with the generalized-K (KG) distribution [7]. An analysis conducted by Zhang et al. further demonstrated the influence of limited phase shifts on achievable performance levels [8]. Outage probability and spectral efficiency have also been deduced in multiple-element IRS-assisted systems [9]. In addition, Boulogeorgos et al. compared the outage probability and ergodic capacity of IRS-assisted and amplify-and-forward (AF) relaying wireless systems [10]. Huang et al. maximized spectral and energy efficiency through transmitting power allocation and meta-surface phase shifts [11]. Guo and Yang et al. optimized the achievable rate through dynamic passive beamforming in an IRSassisted system while considering orthogonal frequency division multiplexing (OFDM) and multi-user channels [12,13].
Physical layer security (PLS), which ensures secure performance through a wireless propagation environment, has received increased attention recently. This is in part due to IRS systems that can flexibly configure a wireless propagation channel with the assistance of IRS. Connection outage probability (COP), secrecy outage probability (SOP) and secrecy rate (i.e., secrecy capacity) are commonly used performance optimization metrics for these systems. COP mainly describes the likelihood that the legitimate receiver (Bob) cannot decode transmitted codewords, and SOP gives the possibility that the eavesdropper (Eve) cannot be confused by secrecy redundancy after the k-th transmission. They indicate the reliability and security, respectively, while secrecy rate shows the overall secure performance [14,15]. For example, Yang et al. developed a closed-form solution of SOP and identified a corresponding asymptotic value in an IRS-assisted communication system using CLT [16]. Khoshafa et al. proposed an SOP for an IRS-assisted D2D system involving Meijer's G-function [17]. Shen et al. investigated secrecy rate maximization for multiantenna IRS systems by combined optimum design of source transmission covariance and IRS phase shift matrices [18]. This joint optimization was also achieved by Cui et al., who solved the problem through beamforming to enhance the secrecy rate [19] while taking into account a complex scenario where the wiretap channel exceeded the main channel. Similarly, Hong et al. applied a covariance matrix of artificial noise (AN) to the joint optimization of a block coordinate descent (BCD) algorithm [20]. This was done to maximize the secrecy rate when the transmitting power limit and unit modulus of the IRS phase shifts were constrained. A fast solution was proposed to solve these constrained optimization problems with a modified genetic algorithm [21]. Yu et al. introduced a robust secure transmission scheme, which considered imperfect channel state information (CSI) in the wiretap channel [22]. This system was designed to improve the sum-rate capacity with constrained highest information leakage for possible eavesdroppers. Wang et al. proposed a joint beamforming and jamming strategy that improved security in a more realistic scenario, in which the CSI was absent [23]. Gu et al. expressed the probabilities of a nonzero secrecy capacity and an ergodic secrecy capacity using closed-form expressions while considering the distribution of eavesdropper locations [24].
As a classical error-control technology, HARQ simultaneously ensures reliability and effectiveness using re-transmission and combining mechanisms. Previous studies have suggested that HARQ is also capable of improving physical layer security. Typical HARQ schemes primarily include chase combining (CC) and incremental redundancy (IR). In the case of HARQ-CC, decoding errors in the receiver trigger the re-transmission of same redundancy versions. Soft information from initial transmissions and re-transmissions is then combined, producing gains in the SNR. Unlike HARQ-CC, HARQ-IR constructs re-transmitted data using various redundancy versions, thus achieving higher coding gain. This approach is not only more efficient, it is also more complex and improvements to its design have been the subject of several recent studies [25][26][27]. In addition, Tang et al. gave the discussion of secrecy outage probability and secrecy throughput in secure HARQ transmissions [28]. Mheich et al. designed a rate adaptive secure HARQ scheme [29], and Guan et al. addressed the reliability-security trade-off and accurate secrecy throughput [30]. It has been demonstrated that an eavesdropper cannot obtain the same diversity gain as a legitimate user, since re-transmissions rely solely on the latter user [31]. As another import metric, effective secrecy throughput (EST) with multiple eavesdroppers has been optimized in secure HARQ transmissions [32]. Additionally, Park proposed Kalman combining-based iterative detection and decoding schemes for multiple-input multiple-output (MIMO) systems with HARQ [33], and Wu et al. achieved a tradeoff between timeliness and reliability for HARQ-CC aided non-orthogonal multiple-access (NOMA) systems [34].
HARQ plays an important role in error-control mechanisms for IRS-assisted systems, as its diversity attributes will inevitably impact reliability and security. Cao et al. addressed the outage rate and progressiveness of Type-1 and HARQ-CC for multi-RIS systems [35]. However, a similar study has yet to be conducted for the more efficient HARQ-IR scheme. Moreover, the performance of PLS in IRS-assisted HARQ systems has yet to be reported. As such, this study first analyzes secure IRS-assisted HARQ performance. Additional contributions of this work include the following: • Closed-form expressions for COP and SOP are derived in IRS-assisted HARQ-CC systems while considering the differences in received SNR between the main channel and the wiretap channel. These two SNR distributions can be accurately approximated by the sum of KG distributed random variables (RVs) with varying parameters. The remainder of this paper is structured as follows. Section 2 describes the overall system model as well as secure transmission without HARQ. Section 3 derives closedform expressions for COP and SOP in IRS-assisted systems with HARQ-CC and HARQ-IR, respectively. Section 4 provides asymptotic values for COP and SOP, along with their corresponding diversity gain and coding gain. Numerical and simulated results are discussed in Section 5. Section 6 concisely concludes our study.
Notation: K v (·) represents the modified Bessel function of the second kind and order v, Γ(·) denotes the gamma function, E[·] indicates the expectation operator, G m,n p,q [·] is Meijer's G-function [36], (·) T indicates the transpose of a matrix, and R and Z represent the sets of real numbers and integers, respectively.

System Model
We take into account an IRS-assisted secure HARQ transmission model, as depicted in Figure 1. In this scenario, an obstacle is positioned between the transmitter (Alice) and the legitimate receiver (Bob). The transmitted signal reaches Bob after being reflected by the IRS while the passive eavesdropper (Eve) intercepts signals from the IRS. These two components of the cascade channel (i.e., from Alice to the IRS and from the IRS to Bob/Eve) are both assumed to be Rayleigh fading channels. The Wyner secrecy code [37] was used for secrecy encoding, with a code rate and secrecy redundancy rate of R B and R E , respectively. In the event of erroneous decoding by Bob, he will transmit a negative acknowledgment (NAK) message to Alice, requesting a re-transmission. Otherwise, an acknowledgment (ACK) message is provided for the new transmission, where a re-transmitted code word is constructed using redundancy versions based on HARQ-CC/IR protocols. Re-transmission ceases once a maximum transmission number K is reached. In this process, Eve achieves limited gain through multiple signal versions, since re-transmission relies only on decoding results from Bob.

Secure Transmissions Without HARQ
We first analyze secure transmissions without re-transmission. In this process, power scaling laws [38] were employed to fix the transmitting power for a single meta-surface path to P. The corresponding transmitted symbol x then satisfies E |x| 2 = P, where channel fading items for the i-th meta-surface from Alice to the IRS and from the IRS to Bob are represented by h i and g i , respectively. These paths constitute a cascade main channel with both symmetric complex Gaussian distribution (i.e., h i ∼ CN (0, 1) and g i ∼ CN (0, 1)). In the following sections, the shifted IRS phase is denoted by φ i , n B is additive Gaussian white noise (AWGN) for the whole cascade main channel and n B ∼ CN (0, σ 2 B ). The magnitude of the received signal in the main channel can then be expressed as These two channel fading terms can be represented using complex signals as h i = α i e −jθ i and g j = β i e −jψ i , where α i and β i are the corresponding magnitude terms and θ i and ψ i are the corresponding phase terms. The received SNR of the main channel is This SNR is maximized when whereγ B denotes the average received SNR of the main channel and X i = α i β i . Since the cascade channel is double Rayleigh distributed, X i is also generalized-K (KG) distributed.
In addition, is accurately approximated by a squared KG distribution [7,39]. Hence, the probability density function (PDF) of R 1 can be obtained by where k m and m m are shaping parameters, k m = −b+ where Ω m is the mean power of R 1 which can be obtained by its second moment denoted by µ R 1 (2). The calculation of a, b, c and µ R 1 (n) can be found in [39]. The PDF of γ B can be expressed as whereΞ m = k m m m γ B Ω m . The closed-form expression of the corresponding cumulative distribution function (CDF) can be acquired from In the case of the wiretap channel, which is the cascade channel constructed from components of Alice-IRS and IRS-Eve, channel fading and phase adjustments from Alice to IRS are consistent with the main channel. In other words, h i and φ i are defined as shown above. The other part of channel fading from IRS to Eve is indicated by g i and thus g i ∼ CN (0, 1), where n E is the AWGN of the whole cascade wiretap channel and n E ∼ CN (0, σ 2 E ). The received signal in the wiretap channel is then given by After adjusting φ i given in the main channel, h i e jφ i and h i follow the same distribution due to the characteristics of a cyclic symmetric complex Gaussian distribution. Assuming h * i = h i e jφ i , the received SNR of the wiretap channel is given by The PDF of the cascaded channel amplitude is known as a KG distributed model [40] and the corresponding CDF follows a squared KG distribution. Hence, the PDF of γ E is given by is a modified Bessel function of the second kind (order v) and k w = N and m w = 1 are shaping parameters. In addition,

Outage Probability of IRS-Assisted Secure HARQ
In this section, we discuss the outage probability (OP) performance of IRS-assisted secure HARQ communication systems, involving connection outage probability (COP) and secrecy outage probability (SOP). These two metrics are derived in HARQ-CC and HARQ-IR, respectively, and are used to evaluate the critical secure performance.

OP of IRS-Assisted Secure HARQ-CC
The HARQ-CC case is first analyzed, in which a connection outage occurs if less mutual information of the main channel is accumulated compared with R B . This corresponding COP is given by where I CC B (K) denotes the accumulated mutual information after the K-th transmission. When adopting maximum ratio combining (MRC), I CC B (K) = log 2 1 + ∑ K k=1 γ B,k , where γ B,k is the received SNR of the k-th transmission in the main channel. This term exhibits a KG distribution and its PDF can be expressed by Equation (5).
The sum of multiple KG-distributed RVs can be accurately approximated by a squared KG distribution [41]. The corresponding PDF is thus where the fading parameters are k M = Kk m and m M = Km m . Additionally, where Ω M = KΩ m . After integrating the above formula, the CDF of Υ B,K can be obtained by The COP of HARQ-CC can be acquired by instituting F Υ B,K (γ) in Equation (12). The probability of the k-th transmission is then given by where M denotes the actual transmission number and P CC e (0) = 1. Hence, the average transmission number of HARQ-CC is Since a secrecy outage occurs when the mutual information is greater than the secrecy redundancy rate in the wiretap channel, the SOP of HARQ-CC after the k-th transmission can be expressed as where I CC E (K) represents the wiretap channel mutual information accumulated by HARQ-CC after the K-th transmission and γ E,k denotes the received SNR at Eve during the k-th transmission, with a distribution given by Equation (9). We assume Υ E,K = ∑ K k=1 γ E,k and it has been already demonstrated that the sum of the aforementioned KG distributed RVs still follows a KG distribution with relevant parameters. The corresponding PDF can be represented by where k W = KN and m W = K are shaping parameters andΞ W = √ KΞ w . The CDF of Υ E,K is then given by Since the actual transmission number M is solely dependent on the main channel, the SOP of HARQ-CC can be expressed as where F Υ B,k (γ) and F Υ E,K (γ) are defined by Equations (14) and (19), respectively.

OP of IRS-Assisted Secure HARQ-IR
Unlike HARQ-CC, HARQ-IR forms code words using different redundancy versions, rather than the same re-transmitted data. After soft-combining, mutual information (not SNR) is directly accumulated. For the main channel, the COP of HARQ-IR is thus where accumulated mutual information after the K-th transmissions is given by We can assume Z B,K = ∏ K k=1 Y B,k and Y B,k = 1 + γ B,k , where Y B,k is a shifted squared KG distributed RV and γ B,k is given in Equation (5). Thus where the CDF of Z B,K is denoted by F Z B,K (·) and γ B,k follows a squared KG distribution. However, the PDF and CDF of Z B,K are difficult to determine. Therefore, we resort to applying the mixture gamma (MG) distribution [42] to express the PDF of γ B,k as where . The amount of summation terms in this distribution is denoted by L, while ω i and t i are the weight and abscissas factors in the Gaussian-Laguerre integration [43]. The k m , m m andΞ m were defined in Equation (5). Furthermore, Y B,k can be expressed as a combination of the shifted Gamma distributed RVs, with a PDF given by A Mellin transform of f Y B,k (y) can then be used to derive the distribution of Z B,K as where Ψ(α, γ; z) = 1

Analysis of Diversity and Coding Gain
In this section, the asymptotic COP and SOP of both HARQ-CC and HARQ-IR are derived at higher SNR. The corresponding diversity and coding gain are then derived from these expressions, which reveal critical factors of secure performance.

Gain of IRS-Assisted Secure HARQ-CC
We first consider the diversity gain and coding gain of IRS-assisted secure HARQ-CC. Equations (12) and (14) can be used to reformulate the COP of HARQ-CC as where Applying the asymptotic series expansion of Meijer's G-function (provided in Appendix B) whenγ B → ∞ gives that the COP can be approximated as It is worth noting that, if κ 2,l = κ 2,k , we can define Γ(κ 2,l − κ 2,k ) = Γ( ), where indicates a minor error term involved to meet the requirements of asymptotic series expansion. According to a definition in [46] The diversity gain and coding gain in the wiretap channel of HARQ-CC can then be obtained by defining the connection outage probability of the wiretap channel as

Equation (19) then becomes
where κ 3 = 1 and κ 4 =MN,M, 0. Whenγ B → ∞ andγ E → ∞, this impliesM 1 and κ 4 N, 1, 0. Applying a similar method in the main channel produces an outage probability of The diversity gain and coding gain in the HARQ-CC wiretap channel can be de- While the coding gain is

Gain of IRS-Assisted Secure HARQ-IR
We next consider the diversity gain and coding gain of IRS-assisted secure HARQ-IR. Whenγ B → ∞, the COP of HARQ-IR can be represented as where γ B,k =γ B R 2 k and R k denotes the signal magnitude of the k-th transmission. The PDF of Y B,K = ∏ K 1 R k can be determined using a moment generation function (MGF)-based method [47] by In other words, A corresponding integration implies the CDF of Z K is given by Equation (45) suggests the COP of the main channel can be expressed as where κ 5 = 1 and κ 6 = k m , m m , · · · , k m , m m 2K items , 0. As stated in Appendix B, whenγ B → ∞, The diversity gain and coding gain in the HARQ-IR main channel are denoted by G IR d,m and G IR c,m , respectively, with the former expressed as Assuming κ 6,k * = min(k w , m w ) = κ 2,k * , the coding gain can be represented as Since the actual transmission number of HARQ-IR isM = 1 when γ B → ∞, the same gain can be acquired as HARQ-CC. The diversity gain in the wiretap channel is then and the coding gain in the wiretap channel is (54)

Numerical Results
In this section, through Monte Carlo simulations and numerical results obtained from the analytical expressions above, we analyze the outage probability performance of IRS assisted secure HARQ, using the model shown in Figure 1. The main secure performance metrics, including COP, SOP and their asymptotic value of HARQ-CC and HARQ-IR, are primarily discussed with a given constraint of maximum transmission number. Figure 2 shows COP curves plotted againstγ B of a secure HARQ-CC system assisted by an IRS, demonstrating the influence of different re-transmission numbers and metasurfaces. The parameters were set asγ B = −10 ∼ 20 dB, R B = 5, R E = 2, with five pairs of curves configured as different combinations of K = 1, 2, 4 and N = 1, 2, 4. The solid lines and asterisks represent the theoretical and simulated COP, respectively. It can be observed that the differences between these two curves are trivial, verifying the theoretical accuracy of Equations (12) and (14). In addition, the COP decreases monotonically with the increase ofγ B , as higher quality of the main channel resulted in fewer connection outages. The COP also decreased significantly with increasing K, due to an accumulated SNR for HARQ-CC. The increase in N allowed the IRS to provide more channels through meta-surfaces, thereby reducing connection outage. Furthermore, it is also evident that increasing N had more of an effect on the COP of HARQ-CC than increasing K (e.g., the results of N = 2, K = 1 outperforms N = 1, K = 2 and N = 4, K = 2 outperforms N = 2, K = 4). This is because chase combining adopts simple repeated transmission with a limited gain.  Figure 3, which depicts SOP performance of an IRS-assisted HARQ-CC system, demonstrates the influence of different transmission numbers and meta-surface quantities. SOP curves were plotted against different values ofγ E , includingγ B = 10 dB, γ E = −10 ∼ 20 dB, R B = 5, R E = 2 and five additional pairs of curves configured as given in Figure 2. It is evident the theoretical SOP (solid line) values agree well with the simulated results (asterisks). In addition, increases inγ E are strongly correlated with an increase in SOP, which indicates that higher quality in the wiretap channel worsens secrecy outage. Furthermore, increases in K were seen to cause SOP increases, though the increment of N = 2, K = 4 (relative to N = 2, K = 2) was less than that of N = 2, K = 2 (relative to N = 2, K = 1). This is because the COP is higher during the first transmission than at the second round. In other words, more re-transmissions occur for lower gain in the wiretap channel. In addition, the SOP deteriorates significantly with increasing N (e.g., the results for N = 4, K = 2 are higher than for N = 2, K = 2). This is caused by more assisted channels implemented through meta-surfaces for both the main channel and wiretap channel.   Figure 2. This analysis of the COP was also verified by theoretical and simulated results, denoted by the solid lines and asterisks, respectively. The COP in Equation (29) decreases monotonically asγ B increases, indicating higher quality in the main channel with fewer connection outages. The COP decreased significantly as K increased, indicating the coding gain in HARQ-IR can be determined for better connection performance. The IRS then provides more channels as N increases, thereby reducing the COP. It is worth noting the effects of increasing N and increasing K differed for HARQ-IR under various conditions. For instance, the COP for N = 1, K = 2 outperformed that of N = 2, K = 1 and N = 2, K = 4 outperformed N = 4, K = 2. As such, the maximum transmission number should be increased at higher SNR, while meta-surface quantities should be increased at lower SNR. This is because connections benefit more from channel diversity than from coding in the case of lower SNRs and vice versa, as verified in Figure 5.  COP curves of the HARQ system assisted by IRS are plotted with their asymptotic values in Figure 5, for a case whereγ B is large. Parameters were set asγ B = 0 ∼ 60 dB, R B = 10 and R E = 5. The solid and dotted lines represent theoretical and asymptotic HARQ-CC and HARQ-IR COP values for N = 2, K = 4 and N = 4, K = 2, respectively. The simulated values have already been verified and are not repeated here for clarity. It is evident that whenγ B → ∞, each theoretical COP value comes first at lower SNR and tends to an asymptotic limit. The COP of N = 4, K = 2 is then better than that of N = 2, K = 4 when HARQ-CC is adopted. Thus, the gain resulting from an increase in the number of meta-surfaces is larger than that caused by increasing the maximum transmission number, which is consistent with our analysis from Figure 2. Furthermore, the COP for N = 2, K = 4 is superior to that of N = 4, K = 2 when HARQ-IR is adopted. In other words, increasing the maximum transmission number is more effective than increasing the number of meta-surfaces, which is concordant with the high-SNR results depicted in Figure 4. Figure 6 plots SOP performance of an IRS-assisted HARQ-IR system as a function ofγ E , N and K. The simulated parameters are the same as those of the SOP simulation for HARQ-CC as shown in Figure 3. The difference between theoretical and simulated results is slight, while the former is given in Equation (35). As shown, the SOP continues to increase withγ E and K, though the increment of N = 2, K = 4 relative to N = 2, K = 2 is lower than that of N = 2, K = 2 relative to N = 2, K = 1. The observed gain in the wiretap channel also decreased more sharply than that of HARQ-CC, due to stronger coding gain for HARQ-IR. In addition, the SOP further increased as N increased, though the increment decreased with increasing N. For instance, the SOP of N = 2, K = 2 increased significantly relative to N = 1, K = 2, while that of N = 4, K = 2 changed little compared to N = 2, K = 2. This was due to a decrease in COP, a decrease in the actual transmission number and corresponding limited IRS-assisted gain.  provides an asymptotic SOP analysis of the IRS-assisted HARQ system. Curves are included for the 1-SOP and its asymptotic values in HARQ-CC and HARQ-IR for increased clarity. Among these curves,γ B = 30 dB,γ E = 0 ∼ 60 dB and all other parameters are the same as in Figure 5. The solid and dashed lines represent the theoretical and asymptotic SOP values, respectively. Each 1-SOP theoretical value tends to an asymptotic limit whenγ E → ∞, which supports our presented analysis. In addition, the 1-SOP curves for HARQ-CC and HARQ-IR tended to be consistent for a given N configuration, such as N = 2, K = 2 and N = 2, K = 4 or N = 4, K = 2 and N = 4, K = 4.
The asymptotic SOP values were also the same, which further verifies our analysis of SOP (i.e., the asymptotic values are only dependent on N with given R E ).

Conclusions
In this paper, we derived closed-form COP and SOP expressions of HARQ-CC and HARQ-IR using the KG distribution, the MG distribution and the Mellin transform. The asymptotic value analysis of the COP and SOP was performed to determine the diversity and coding gain. An analysis of numerical and simulated results suggested that expanding the number of meta-surfaces had a significant impact on the main channel for all SNR regions in HARQ-CC and lower-SNR regions in HARQ-IR while increasing the maximum transmission number critically affected the main channel for higher-SNR regions in HARQ-IR. Additionally, only the number of meta-surfaces mainly affected the wiretap channel for both HARQ-CC and HARQ-IR. On this basis, corresponding strategies could be further designed to effectively improve physical layer security performance for an IRS-assisted HARQ system. Substituting Equation (27) into Equation (A1) gives f Z B,K (z) = 1 2πi c+i∞ c−i∞ ∑ r 1 ···r L r 1 +···+r L =K K r 1 · · · r L L ∏ j=1 α m,j Γ(β m,j )Ψ(β m,j , β m,j + s; ζ m,j ) r j z −s ds = ∑ r 1 ···r L r 1 +···+r L =K K r 1 · · · r L L ∏ j=1 α m,j Γ(β m,j ) r j 1 2πi c+i∞ c−i∞ L ∏ j=1 Ψ(β m,j , β m,j + s; ζ m,j ) r j z −s ds (A2) For convenience of presentation, we define Ψ(β m,j , β m,j + s; ζ m,j ) r j z −s ds.