Guaranteeing QoS for NOMA-Enabled URLLC Based on κ–μ Shadowed Fading Model

Sixth-generation (6G) wireless communication scenarios are complex and diverse. Small-scale fading is a key part of wireless channels and its impact on performance in scenarios with time sensitivity and 6G ultrareliable and low latency communications (URLLC) quality-of-service requirements cannot be ignored. Therefore, it is necessary to accurately characterize small-scale fading when designing wireless communication systems. In this paper, we derive approximate closed form expressions for the probability density function, cumulative distribution function and moment-generating function of the postprocessing signal-to-noise ratio following the zero-forcing detector in a cell-free massive multiple-input multiple-output (CF mMIMO) system. CF mMIMO system is a nonorthogonal multiple access (NOMA) system that enables users to share all channel uses and can ensure the fairness of the communication quality experienced by different users. Our key contributions include the extension of the κ–μ shadowed fading model to a CF mMIMO system and the proposal of theoretical tools (the derived closed-form expression) to improve its mathematical tractability. By exploiting the statistical characterizations of the arrival and service processes, another important contribution is the exploitation of the upper bound of the queuing delay violation probability (UB-QDVP) over the Mellin transforms of the arrival and service processes in the proposed CF mMIMO system under the κ–μ shadowed fading model. Corroborated by extensive simulations, our analyses validate that the CF mMIMO system outperforms the orthogonal multiple access and power-domain NOMA systems and reveal the relationships among different small-scale fading types, energy efficiency, delay and the UB-QDVP, as well as the accuracy and effectiveness of the proposed theoretical tools based on the κ–μ shadowed fading model.


Background and Motivation
Due to the highly time-varying nature of wireless channels, it is difficult to design communication systems that support ultrareliable and low latency communications (URLLC) with reductions in energy consumption. For example, the explosively increasing volumes of delay-sensitive and bandwidth-intensive applications under stringent quality-of-service (QoS) requirements have required strict delay (<0.1 ms), ultrareliability (>99.99999%), and extra-high energy efficiency (EE) demands [1]. Notably, achieving higher data rates and lowering power consumption seem to be contradictory objectives. In particular, when the delay boundary approaches the channel coherence time, the transmit power required to obtain the QoS requirements becomes unbounded [2]. Although both academia and industry have paid attention to the EE of cellular networks in the past few years [3], the related works are still unable to address the increasing complexity of devices, antennas and many frequency bands in the future [4].
Orthogonal multiple access (OMA) systems (e.g., multiple-input multiple-output orthogonal frequency division multiple access (MIMO-OFDMA)) apply the orthogonal resource allocation method, in which different uplink users occupy different channel uses (CUs). This type of system distinguishes different users through different orthogonal time-frequency resource blocks. However, all users share time-frequency resources in the non-OMA (NOMA) system. As a promising technology, NOMA makes full use of the nonorthogonal nature over limited orthogonal resources so as to provide a higher transmission accuracy with lower computational complexity [5], and to meet the massive connectivity requirements in an Internet of things scenario [6,7] and even satellite-based communications [8]. For instance, power-domain NOMA (PD-NOMA), considered as the representative of single carrier NOMA technology [9], distinguishes different groups of users through the spatial domain, and different users in the same group share spatial domain resources. Under PD-NOMA systems, power control is required to meet the QoS requirements of URLLC for energy efficiency, latency and reliability [10]. However, fairness among the users in terms of their throughput is an important goal of NOMA power allocation [11]. For the upcoming sixth generation (6G) mobile networks, a beneficial embodiment of a distributed massive MIMO system is the cell-free massive MIMO (CF mMIMO) system proposed in [12]. As a type of NOMA system, the CF mMIMO system is a better solution between PD-NOMA and MIMO-OFDMA. All users share all CUs, and different users are distinguished by the spatial domain without considering the fairness experienced by uplink users. The CF mMIMO system contains a large number of randomly distributed access points (APs) to coherently serve massive mobile devices; its implementation is more feasible for practical issues than network MIMO with low overhead that exploits channel hardening phenomena and favorable propagation properties [4,13].Relative to traditional cellular mMIMO networks, the notable characteristics of CF mMIMO system are its superior macro diversity and interference suppression capabilities. Macro diversity offers orders of magnitude of coverage probability [12,14], reliability and delay optimization through suitable signal processing. Intercell interference is effectively mitigated by the user-centric approach and coordination [15,16], thereby improving the accuracy of multiuser detection. To the best of our knowledge, a CF mMIMO system that efficiently integrates the unified fading model to provide substantial performance improvements with respect to different types of small-scale fading has still not been well studied.
On the other hand, most current papers concerning the performance analysis in CF mMIMO systems [4,[17][18][19] perform analyses in the context of small-scale fading models such as the Rayleigh and Nakagami-m models, which ignore the random nature of wireless fading channels. Such classic fading models usually do not cover all types of fading scenarios encountered in practical cases. More importantly, the analyses performed for practical channel estimations depend on the channel's statistical properties, which are affected by the selected channel models [20]. Consequently, many results are not reasonably universal enough due to simplified assumptions, such as in many studies of CF mMIMO systems performance based on the Rayleigh model [12,14,21] or Rician model [22]. The design of reliable wireless communication systems is highly dependent on the accurate characterization of the wireless fading channel. However, for the complex communication scenarios in future 6G URLLC scenarios [23,24], the classic fading model is no longer applicable [25,26]. With the research on fading models in recent years, a new fading model that can losslessly portray a variety of fading types has the potential to become a unified fading model [27][28][29][30], which can be used to accurately characterize fading channels in complex future 6G URLLC scenarios. This is a critical issue regarding the theoretical analysis and performance optimization of URLLC, given that the reliability of modeling is of capital importance. Specifically, the κ-µ shadowed fading model, which unifies the κ-µ fading distribution and the η-µ fading distribution, was introduced in [27]. First, by adjusting the parameters κ, µ and m, almost all classical fading models can be unified by the κ-µ shadowed fading model [28] with similar mathematical flexibility levels. Second, this shadowed model provides much better tractability and flexibility [31]. These two features allow the κ-µ shadowed fading model to accommodate different propagation conditions, such as in practical communication environments that experience the effects of nonhomogeneous fading conditions. However, although the κ-µ shadowed fading model can characterize a variety of fading channels, adapting the κ-µ shadowed fading model to the CF mMIMO system and specifically analyzing the system performance while supporting URLLC becomes a challenging research issue.
With respect to mathematical tools, towards the exponentially growing demands of time-sensitive services, the traditional Shannon information theory based on the infinite blocklength assumption and a random coding scheme provides a poor benchmark in this scenario [32]. Finite-blocklength (FBL) information theory was proposed in [33] as a shortpacket data communication technique for guaranteeing the stringent QoS requirements of URLLC [34]. Well-known results regarding the maximum achievable data rate of FBL packets were derived in [33] as a function of the blocklength and error probability. Moreover, as suggested by Bennis et al. [35] and Arnau et al. [36], the design of URLLC should capture the tail behavior of reliability (i.e., queuing delay violation probability) and delay instead of average metrics. Because the delay and overall reliability consist of multiple components, the corresponding performance indicators should be constrained by a delay boundary and a delay violation probability boundary for URLLC [36]. In this sense, the stochastic network calculus (SNC) approach proposed in [37] provides a nonasymptotic boundary on the delay violation probability as an effective analysis tool for connecting the delay and overall reliability. Therefore, the CF mMIMO system under the κ-µ shadowed fading model is sufficient for analyzing performance in practical communication environments to meet the strict QoS requirement of URLLC. However, a major problem with this kind of performance analysis is appropriately characterizing stochastic queuing behaviors when they are integrated with the CF mMIMO system under the κ-µ shadowed fading model with SNC and FBL information theory. Furthermore, it is still an urgent issue to exploit the effect of smallscale fading characteristics on the performance of the CF mMIMO system under the κ-µ shadowed fading model to achieve NOMA-enabled URLLC with high resource efficiency.

Related Works
Our paper builds on results derived from several research areas. On the physical layer, we consider the uplink CF mMIMO system with imperfect channel state information (CSI) under the κ-µ shadowed fading model as well as FBL information theory. On top of a service description in the physical-layer fading channel model, we then investigate service processes. We determine that the effect of different types of small-scale fading on the performance in the design, modeling and analysis of several wireless communication systems cannot be neglected. The CF mMIMO system has been gradually considered in the URLLC scenario. Nasir et al. [38] first investigated a CF mMIMO system for downlink URLLC in the FBL regime and designed a special class of conjugate beamforming to achieve a lower computational complexity and improved URLLC rate performance. Elwekeil et al. [19] introduced power control schemes in a CF mMIMO system that supported URLLC applications for both traditional ground users and unmanned aerial vehicles, and they evaluated the superior rate performance of this CF mMIMO system. Due to the spatial domain, the CF mMIMO method can be a key to boosting the reliability of URLLC applications. The study in [39] compared different transmission modes in a factory automation scenario and demonstrated that the user-centric transmission mode and power control brought substantial performance improvements in terms of reliability and latency. Zhang et al. [40] integrated the CF mMIMO approach with simultaneous wireless information and power transfer for statistical delay and error-rate bounded QoS to support URLLC and boost the data rate and energy efficiency. Hence, a CF mMIMO network can be viewed as a promising system to support URLLC with a high resource efficiency and low latency.
The κ-µ shadowed fading model [27,28] was proposed as a more general yet equally tractable model to capture a wide range of propagation conditions with clear physical interpretation and good analytical properties. Furthermore, Lopez-Martinez et al. [31] derived the probability density function (PDF) and cumulative distribution function (CDF) of the κ-µ shadowed fading model with integer fading parameters. Ramirez-Espinosa et al. [25] provided a statistical characterization of the κ-µ shadowed fading model and asserted that these results could be used to derive some performance metrics of wireless communication systems over fading channels. With respect to the characterization of the κ-µ shadowed fading model, its applications have been studied in several works. ElHalawany et al. [41] carried out a performance (i.e., ergodic capacity, outage probability and average bit error rate) analysis of downlink NOMA systems subject to the κ-µ shadowed fading model. Chun et al. [42] proposed an analytic framework to evaluate the average of an arbitrary signal-to-noise-plus-interference ratio (SINR) function over the κ-µ shadowed fading model and evaluate the spectral efficiency, moments of the SINR and outage probability of heterogeneous cellular systems. The authors in [43] characterized the impacts of realistic propagation conditions on the achievable secrecy performance of MIMO systems over the κ-µ shadowed fading model. Fully absorbing the attractive nature of the κ-µ shadowed fading model is seen as one of the potential keys to characterizing the propagation media in emerging practical scenarios.
Several works have also examined SNC-based network performance analysis and theoretical boundary calculation for URLLC. Schiessl et al. used SNC to investigate the delay performance of a multiuser MIMO system with zero-forcing beamforming [44] and NOMA with joint decoding [45] under imperfect CSI and FBL. Xiao et al. [46] derived a closed-form of the upper bound of the queuing delay violation probability (UB-QDVP) in a downlink NOMA system and designed an optimal power allocation scheme for guaranteeing the delay violation probability. Furthermore, by utilizing moment generating function (MGF)-based SNC, probabilistic delay bounds for traffic dispersion, network densification and a hybrid scheme were obtained in millimeter-wave communications over Nakagami-m fading [47]. A study close to our work is the research by Zhang et al. [22], where the stochastic QoS performance in terms of both the delay and error rate metrics was obtained for a downlink CF mMIMO based system model across Rician fading in the FBL regime. However, only traditional fading channels have been considered for discussion in the related literature, while no investigation is available with respect to the performance analysis for CF mMIMO systems with the juxtaposition of shadowing and small-scale fading.

Main Contributions
For 6G URLLC requirements, we combine the potential key technology of future 6G, CF mMIMO [48], to build the uplink wireless communication system, and extend the κ-µ shadowed fading model to the CF mMIMO system, for lossless characterization of its fading channel in the complex communication scenario of 6G [23,24]. The key contributions of the paper can be summarized as follows.

•
We derive approximate closed-form expressions for statistical characteristics (the PDF, CDF and MGF) of the sum of independent and nonidentically distributed (i.n.i.d.) κ-µ shadowed random variables (RVs). The analysis is nontrivial, as we optimize the theoretical analysis performance of the κ-µ shadowed fading model and improve its mathematical tractability. • Based on the κ-µ shadowed fading model, we derive approximate closed-form expressions for the PDF, CDF and MGF of the postprocessing signal-to-noise ratio (PPSNR) after the zero-forcing detector in the proposed CF mMIMO system, and extend the κ-µ shadowed fading model to 6G NOMA wireless communication systems. • By utilizing FBL information theory, SNC and the Mellin transform on the service process, we exploit the UB-QDVP in the proposed CF mMIMO system under the κ-µ shadowed fading model. Furthermore, based on extensive simulations, we analyze the system performance with the delay and UB-QDVP indicators and validate the necessity of analyses performed under the κ-µ shadowed fading model; the CF mMIMO system outperforms the OMA system and PD-NOMA system. (For the OMA and PD-NOMA systems, there are L antennas equipped in their base station and K users. For illustration convenience, in the PD-NOMA system, we divide K users into N pairs based on the channel gains of the users (K = 2N ). Specifically, we divide K users into a group of N "strong users" and a group of N "weak users" according to their channel gains and sort the users in the descending order of their channel gains within each of the groups. The "strong user" and "weak user" with the same label are grouped in the same pair.)

Organization
The remainder of this paper is organized as follows. In Section 2, the system model for a CF mMIMO system under the κ-µ shadowed fading model is described. Section 3 introduces SNC and presents the Mellin transform over both arrival and service processes as well as the UB-QDVP. Simulation results and an analysis are given in Section 4, followed by the conclusions in Section 5.
Notations : C K×N is a complex matrix with K rows and N columns. The modulus and expectation operators are denoted by |•| and E(•), respectively. (•) T and (•) H denote the transpose and conjugate transpose, respectively, of their arguments. (•) lk denotes the element in the lth row and kth column of a matrix. CN (0, 1) denotes a complex Gaussian distribution with a mean of zero and unit variance.

System Model
We consider an uplink CF mMIMO system in which L single-antenna APs and K singleantenna users are randomly uniformly distributed in a circular area with L K, as shown in Figure 1. Some conventional fading can be obtained by adjusting the parameters of the κ-µ shadowed fading model: and m − are nonnegative real numbers. When m − → ∞, it denotes that no shadow effect is observed). All APs simultaneously serve all users at the same time-frequency resource and are connected to a baseband unit (BBU) for centralized signal processing via ideal backhaul links [19,49] (There is no delay in the ideal backhaul links between APs and the BBU [50]). The transmission from the APs to the users (downlink transmission) and the transmission from the users to the APs (uplink transmission) proceed by a time-division-duplexing operation. Each coherence interval is divided into three phases: uplink training, downlink payload data transmission, and uplink payload data transmission. In the uplink training phase, the users send pilot sequences to the APs and the BBU estimates the channel for all users. In this paper, the QoS delay for URLLC is characterized by the queuing delay (Considering ideal backhaul links and tail behaviors based on [35,36], we consider the queuing delay as the delay of the URLLC QoS requirement). Each user transmits a short packet with D information bits through N CUs, which are spread over B MHz of bandwidth and t f milliseconds of time (N = Bt f ). The channel coefficient g lk (t) between the lth AP and the kth user can be characterized as where β lk (t) is a large-scale fading coefficient that changes very slowly with time and is assumed to be constant over many coherence conditions. h lk (t) is the small-scale fading coefficient of the corresponding links, which is assumed to be static during a finite-sized time-frequency coherence interval, and it changes independently from one coherence interval to the next [51]. The channel response is assumed to follow the κ-µ shadowed fading model. The gain of small-scale fading between the kth (k ∈ [1, K]) user and the lth AP (l ∈ [1, L]) can be modeled as where µ lk is a natural number. Each multipath cluster is modeled by one term of the sum; thus, µ lk is the number of multipath clusters. Each cluster has scattered components with the same power and dominant components with a certain arbitrary power. X lk,i and Y lk,i are independent and identically distributed (i.i.d.) Gaussian RVs with means of zero and variances of σ 2 , and X lk,i + jY lk,i denotes the scattered components of the ith cluster.p lk,i andq lk,i are real numbers and the mean values of the in-phase and quadrature parts of the multipath components of the ith cluster, respectively. ξ lk is a Nakagami-m RV with shaping parameter and E ξ 2 lk = 1, where ξ 2 lk ∼ Γ(m lk , 1/m lk ); m lk represents the fluctuation degree of dominant components caused by the shadow [28]. For each cluster, the total power of the scattered components is 2σ 2 . Furthermore, the ratio of dominant components to scattered In addition, |h lk | 2 statistically follows the κ-µ shadowed distribution, whose PDF can be expressed approximately by the Gamma ; Ω lk is the mean value of |h lk | 2 . The received signal at the BBU under imperfect CSI can be denoted as in the following equation.
where p u is the average transmit power of each user. Y (p) (t) ∈ C L×n and Y (d) (t) ∈ C L×n are received pilot signal and received data signal transmitted by the users at time slot t, respectively. n andn are the length of the pilots (LoP) and the length of the data signal, respectively, which satisfy n +n = ∈ C K×n is the transmitted data signal at time slot t, obtained by as- As we mentioned above, the expression for the channel coefficient is g lk (t) = β lk (t)h lk (t), ∀k ∈ [1, K], and ∀l ∈ [1, L], which can be simplified as g lk = β lk h lk , g k (t) ≡ g k and G(t) ≡ G because the large-scale fading coefficient β lk (t) and small-scale fading coefficient h lk (t) are generally assumed to be RVs that are independent of time. Z (p) (t) and Z (d) (t) represent additive white Gaussian noise (AWGN), whose elements are i.i.d.; this can be expressed as CN (0, 1). In a practical communication system, serving users with a large number of APs requires accurate CSI to be available. The estimated channel matrixĜ obtained via the least squares estimation method which relies on the received signal with the known pilot sequence is expressed as [52] . β lk and Ω lk are the large-scale fading coefficient and mean square value of the small-scale fading coefficient, respectively.Ê = [ê 1 , . . . ,ê K ] ∈ C L×K is the channel estimation error matrix, and its elementsê k = [ê 1k , . . . , The elements ofĜ can be expressed aŝ Adopting the standard linear zero-forcing detectorÂ =Ĝ(Ĝ HĜ ) −1 ∈ C L×K , the received signal of the kth user is After completing zero-forcing detection, the PPSNR of the kth user is [53] and Unique(T ) returns the same values as those in the set T but with no repetitions.â lk is the lth entry ofâ k , andâ k is the kth column ofÂ. ηl k is the channel gain of the lth AP and the kth user and is an RV with both large-scale fading coefficient and a small-scale fading coefficient. (a) is based on the consideration that when [53]. Because the gain of small-scale fading |h lk | 2 statistically follows the κ-µ shadowed distribution, the distribution of ηl k is obtained in Lemma 1 as follows.

Lemma 1.
Based on a κ-µ shadowed fading model, the channel gain ηl k follows the Gamma distribution, which is where the mean value of ηl k , as well as the mean value of the channel gain.
Proof. Please refer to Appendix A.
Based on the κ-µ shadowed fading model, ζ k = ∑ l∈L k ηl k can be regarded to the sum of i.n.i.d. Gamma RVs. Then, the PDF, CDF and MGF of ζ k are obtained from Theorem 1.

Theorem 1.
Considering a κ-µ distribution with parameters (A k , K k ,Ω k ), the PDF can be approximately obtained by a V-order generalized Puiseux series expansion process as follows. .
We can obtain approximate parameters c v and ω v via the following formula: . . .
The CDF of ζ k is given by where γ(·, ·) is the lower incomplete Gamma function and is expressed by The MGF of ζ k is where G q,p p,q (·) is the Meijer G-function and p q = A k = A n,k = 1. According to gcd(p, q) = 1, both p and q can be 1.
Then, the MGF of ζ k can be simplified to Proof. Please refer to Appendix B.
Based on Equation (6) and Theorem 1, the statistical characteristics of the PPSNR in the CF mMIMO system are provided in the following corollary. Corollary 1. Due toγ k ≈p u ζ k , the PDF of the PPSNRγ k at the BBU is derived as The CDF ofγ k is given by The corresponding MGF ofγ k can be expressed as

Analysis of the Delay Performance
Because each packet is very short, both imperfect CSI and FBL channel coding cause transmission errors. The maximum achievable data rate is subject to transmission errors and an instantaneous signal-to-noise ratio (SNR), as described in Section 3.1. Due to transmission errors and dynamic data rates in the physical layer, the APs need to keep their packets in a buffer until successful completion, which leads to a random queuing delay, as discussed in Section 3.2. Finally, in Section 3.3, we derive the Mellin transform over the arrival and service processes and investigate the UB-QDVP.

Service Transmission in FBL Regime
Based on FBL information theory, the maximum achievable data rate (bit/CU) with an instantaneous SNR and a decoding error probability for the kth user can be obtained as follows. [22,33] where ,n is the number of CUs occupied by data transmission and ε d is the decoding error probability.
x e −y 2 2 dy is the tail distribution function of the standard normal distribution. C(γ k ) = log (1+γ k ) 2 and V(γ k ) = 1 − 1 (1+γ k ) 2 (log e 2 ) 2 are the Shannon capacity and channel dispersion, respectively. Given a decoding error probability of ε d k , the achievable data rate r k γ k , ε d k is less than zero when the PPSNRγ k is below a threshold γ 0 . As a consequence, we redefine Furthermore, the maximum achievable data rate can be simplified as follows:

Queuing Model
For our queuing system analysis, we consider a usual statistical model such as that in [37], which is widely used in the SNC methodology. Then, the cumulative arrival, achievable service, and departure processes of user k during time interval [τ, t) are defined as follows.
The instantaneous traffic arrival a k (i) models the number of bits that arrive at the queue within a discrete time slot i. The increment of service process (service rate) r k (i) is independent of time and equal to the maximum achievable data rate r k γ k , ε d k . The departure process u k (i) denotes the number of bits that arrive successfully at the destination and is subject to both the service process and the amount of data waiting in the queue. In addition, we assume that all queues are work-conserving first-come-first-served queues. Then, the queuing delay w k (t) of user k at time slot t is defined as the number of time slots required for all data that arrive before time t to depart from the transmit buffer and is expressed as

SNC in the SNR Domain
A notable aspect of SNC is that it can obtain the delay distribution and the UB-QDVP based on statistical characterizations of the arrival and service processes in terms of their Mellin transforms. However, it is still difficult to capture the statistics of the random arrival and service processes in practical scenarios. To facilitate the analysis, we convert these processes from the bit domain through the exponential domain to the SNR domain. These cumulative processes in the SNR domain are denoted by calligraphic letters [54] as follows: M X (s) = E X s−1 is the Mellin transform of any nonnegative random variable X for a parameter s ∈ R [37]. Both a k (i) and r k (i) are independent of each other at different time slots, so we use two time-independent random variables a k and r k to denote the number of arrival and service bits at either time slot, respectively. Then, the Mellin transform over the accumulated arrival process is denoted by M A k (s, τ, t) = M A k (τ,t) (s) = If the stability condition M ϕ k (1 + s)M φ k (1 − s) < 1 holds, denoting t T as a target delay (target number of time slots), the kernel function K(s, t T ) is defined as [37] K(s, Thus, the UB-QDVP denoted by ε v,UB k is [22,37] ε v,UB where s 0 = sup s :

Mellin Transform over the Arrival Process in the SNR Domain
As shown in Equation (24), the UB-QDVP is related to the arrival process on the Mellin transform in the SNR domain. Consequently, deriving the Mellin transform of a k is essential for evaluating the UB-QDVP. We consider the traffic class of (δ(s), λ(s))-bounded arrivals, whose arrival processes characterized in the SNR domain for some s > 0 on the Mellin transform are bounded by [37] M A k (s, τ, t) = M ϕ k (s) t−τ ≤ e (s−1)·(λ(s−1)·(t−τ)+δ(s−1)) .
To simplify the notation, we consider the case where δ(s) = 0 and λ(s) ≡ λ; namely, (δ(s), λ(s)) are independent of s, which is true for constant arrivals [54]. Thus, where λ is the constant number of arrival bits that arrive successfully at one slot, as well as the constant bit arrival rate.

Mellin Transform over the Service Process in the SNR Domain
Equation (24) also shows that deriving the closed-form expression for the Mellin transform over the service process in the SNR domain at user k is a key component when analyzing the UB-QDVP, which motivates the following Theorem 2.
. ε d k is the decoding error probability. γ 0 is the upper bound of the PPSNR when the maximum achievable data rate of user k is 0.

Proof. Please refer to Appendix C.
Combining and extending Theorem 2, Equations (24) and (26) yield the following corollary.

Corollary 2. The approximate UB-QDVP of the kth user is
Proof. Please refer to Appendix D.

Simulation Results and Analyses
We present numerical results to validate and evaluate the stochastic queuing behaviors and the effects of fading characteristics on performance in the proposed CF mMIMO system under the κ-µ shadowed fading model. Different from most studies [4,12,13,17,18], which explore the effects of the number of APs, the coverage radio, the number of potential users and other traditional system parameter configurations on performance in communication scenarios, we focus on fading characteristics. Throughout our simulations, we set such parameters constant. For large-scale fading β lk = u p d k d 0 −α p , u p , α p , d 0 and d lk are the constant path loss across the minimum distance, the path loss exponent, the minimum distance and the distance between the kth user and the lth AP, respectively. The key simulation parameters are listed in Table 1. The overall reliability for the kth user can be represented by the overall error probability ε k [55], which is First, we compared the maximum achievable data rates obtained under the same decoding error probability but different fading models in Figure 2, when the transmit power is 5 dBm and the distance from the user to the center of the coverage area is 200 m. It is obvious that the maximum achievable data rates corresponding to the same decoding error probability are inconsistent under different fading environments. This indicates that the effect of different fading environments cannot be ignored with respect to system performance. In particular, when designing a wireless communication system (e.g., the CF mMIMO system) for strict decoding error probability, it is crucial to design customized coding design and power control schemes by combining the fading characteristics of the wireless communication environment to adjust the service rate to meet the reliability requirements of NOMA-enabled URLLC QoS. Therefore, the κ-µ shadowed fading model accommodating different propagation conditions that are applied to design or analyze the performance of a wireless communication system is necessary and significant.  Figure 2. Relationship between the decoding error probability and the maximum achievable data rate. Figure 3 shows the fitting ability of our derived theoretical tool (Theorem 1) for the PDFs of the sum of i.n.i.d. RVs. We used the Monte Carlo method to randomly generate six independent RVs with different distributions and sum some of them. We verify that the cylindrical diagram (obtained by a simulation using the Monte Carlo method for 10 7 iterations) and the curve with an asterisk (obtained by numerical calculation (Equation (8)) have almost the same slope. This indicates that our proposed theoretical tools based on Theorem 1 can accurately characterize the statistical characteristics of the sum of multiple i.n.i.d. RVs. Furthermore, we summarize the second-order approximate parameters obtained under the different numbers of summed elements and different distributions listed in Table 2, which were obtained by using fsolve in MATLAB. From Figure 3, we also observe that the generalized Puiseux series expansion process converges very quickly, so it can effectively fit the PDF without complex approximate orders. All these properties lay the foundation that appropriately characterizing the statistical characteristics of the sum of multiple i.n.i.d. shadowed RVs is the key to unsolved problems. This again verifies the effectiveness and good analytical properties of our proposed theoretical tool (Theorem 1) in complex future communication environments.  Figure 3.
Case 3: For the proposed CF mMIMO system, Figure 4 depicts the PDFs of the PPSNR under different fading environments when the transmit power is 10 dBm and the distance from the user to the center of the coverage area is 200 m. We verify that the cylindrical diagram (obtained by simulation using Monte Carlo method for 10 7 iterations) and the curve with asterisk (obtained by numerical calculation (Corollary 1)) have almost the same slope. This indicates that our proposed theoretical tools based on Corollary 1 have superior flexibility, mathematical tractability and practicality for wireless communication system design. Furthermore, these findings also reveal the κ-µ shadowed fading model based on our proposed theoretical tools, which form a reliable and mathematically tractable channel fading model can be extended to other common future 6G wireless communication systems such as CF mMIMO systems, enabling this model as a more applicable approach to complex 6G applications in the future.    Figure 5 plots the UB-QDVPs obtained with different target delays over a κ-µ shadowed fading model for the CF mMIMO system when the transmit power is 10 dBm and the distance from the user to the center of the coverage area is 200 m. The UB-QDVP decreases as the target delay increases. For the same target delay, the users have inconsistent UB-QDVPs under different fading environments. In Figure 5c, under Nakagami-m fading (κ → 0, µ = 3, m → ∞) and a target delay of 0.05 ms, the UB-QDVP is 10 −7 . However, when the UB-QDVP is 10 −7 under One-sided Gaussian fading (κ → 0, µ = 0.5, m → ∞), the target delay exceeds 0.05 ms. For the performance evaluation, analysis and design of the system, it is necessary to consider the impacts of different fading environments on the system performance. Otherwise, the final result or conclusion obtained will differ greatly from reality, eventually leading to analysis results with less practical value and significance.  These three communication systems are the CF mMIMO system we proposed, PD-NOMA and OMA systems equipped with 200 antennas in their base station. We observe that the UB-QDVP of the PD-NOMA system is lower than that of the OMA system, and the UB-QDVP of the CF mMIMO system is lower than that of the PD-NOMA system under the same target delay when the distance from the user to the center of the coverage area is 200 m. In detail, when the target delay is 0.06 ms and the distance from the user to the center of the coverage area is 200 m, the UB-QDVPs of the PD-NOMA and OMA systems exceed 10 −7 , and only the UB-QDVP of the CF mMIMO system is lower than 10 −7 . Under a target delay of 0.07 ms, the UB-QDVPs of the CF mMIMO system and the PD-NOMA system are below 10 −7 . However, under the same parameter settings, the UB-QDVP of the OMA system is still higher than 10 −7 . Specifically, we extended the distance from the user to the center of the coverage area to 1000 m and set the target delay to 0.08 ms. Only the CF mMIMO system can achieve a UB-QDVP lower than 10 −7 , and both the PD-NOMA and the OMA systems cannot. Figure 6b plots the UB-QDVP for a user at 200 m or 1000 m from the center of the coverage area under Nakagami-m fading (κ → 0, µ = 3, m → ∞) for the CF mMIMO system, PD-NOMA system, and OMA system when the transmit power is 10 dBm. It is also clear from Figure 6b that when the user is 200 m from the center of the coverage area, the UB-QDVP of the CF mMIMO system is lower than that of the PD-NOMA system, which in turn is lower than that of the OMA system, for the same target delay. For example, when the target delay is 0.06 ms, the UB-QDVP is lower than 10 −7 for the CF mMIMO system and PD-NOMA system, and higher than 10 −7 for the OMA system when the user is 200 m away from the center of the coverage area. When the user is 1000 m away from the center of the coverage area and the target delay is 0.07 ms, only the UB-QDVP of the CF mMIMO system is lower than 10 −7 , and the UB-QDVPs of the PD-NOMA system and OMA system are higher than 10 −7 . Figure 6a,b confirm that the CF mMIMO system not only optimizes the reliability and delay performance, but also expands the coverage area for providing URLLC services in comparison with the PD-NOMA and OMA systems.  Figure 7a shows the relationship between EE dB and the UB-QDVP under different fading environments in the proposed CF mMIMO communication system when the transmit power is 10 dBm, the target delay is 0.08 ms and the distance from the user to the center of the coverage area is 200 m. As seen from Figure 7a, as the UB-QDVP becomes tighter, EE dB decreases. Under different fading environments, the EE dB of the user for the same UB-QDVP is different, and a gap is observed. The EE dB of the user under Nakagami-m fading (κ → 0, µ = 3, m → ∞) for controlling the transmit power to meet the UB-QDVP requirement (10 −7 ) is improved by 16.6% in EE dB over that of users under Rician shadowed fading (κ = 4, µ = 1, m = 3) with the same UB-QDVP. Therefore, for delay-sensitive applications, performance analyses for indicators such as delay and reliability, power control and system design, the impacts of different small-scale fading characteristics on the performance of communication systems need to be considered. The accurate characterization of small-scale fading using the κ-µ shadowed fading model may become a necessity for future wireless communication performance analysis and system design.   Figure 7b, we can see that when the target delay is 0.1 ms and the distance to the center of the coverage area is 200 m, the user achieves ε v,UB k =10 −7 with EE dB < 40 bit/J in the OMA system. Under the same conditions, the user in the PD-NOMA system can arrive at the same UB-QDVP with EE dB > 40 bit/J, which is an EE dB improvement of 5.1% over that of the OMA system. Furthermore, the user in the CF mMIMO system achieves a 2.4% increment in EE dB compared with that of the PD-NOMA system. When EE dB = 27.65 bit/J, the target delay is 0.1 ms and the distance to the center of the coverage area expands to 1000 m, the UB-QDVPs of the user in both PD-NOMA and OMA systems are greater than 10 −7 . However, the UB-QDVP of the user in the CF mMIMO system is less than 10 −7 and is 97.5% lower than that of the PD-NOMA system. In Figure 7c, when the target delay is 0.1 ms and the distance to the center of the coverage area is 200 m, the user achieves ε v,UB k =10 −7 with EE dB < 40 bit/J in the OMA system. Under the same conditions, the user in the PD-NOMA system can arrive at the same UB-QDVP with EE dB = 40 bit/J, which is an EE dB improvement of 11.1% over that of the OMA system. Furthermore, the user in the CF mMIMO system achieves a 13.9% increment in EE dB compared with that of the OMA system. When EE dB = 27.65 bit/J, the target delay is 0.1 ms and the distance to the center of the coverage area expands to 1000 m, the UB-QDVPs of the user in both PD-NOMA and OMA systems are greater than 10 −7 . However, the UB-QDVP of the user in the CF mMIMO system is less than 10 −7 and is 99.4% lower than that of the PD-NOMA system. In Figure 7d, when the target delay is 0.1 ms and the distance to the center of the coverage area is 200 m, the user achieves ε v,UB k =10 −7 with EE dB = 40 bit/J in the OMA system. Under the same conditions, the user in the PD-NOMA system can arrive at the same UB-QDVP with EE dB > 40 bit/J, which is an EE dB improvement of 4.9% over that of the OMA system. Furthermore, the user in the CF mMIMO system achieves a 6.25% increment in EE dB compared with that of the OMA system. When EE dB = 27.65 bit/J, the target delay is 0.1 ms and the distance to the center of the coverage area expands to 1000 m, the UB-QDVPs of the user in both PD-NOMA and OMA systems are greater than 10 −7 . However, the UB-QDVP of the user in the CF mMIMO system is less than 10 −7 and is 99.6% lower than that of the PD-NOMA system. Hence, the CF mMIMO approach is a prospective technology that can improve EE while ensuring that the target delay and reliability meet the URLLC QoS requirements. At the same time, it can greatly decrease the UB-QDVP and improve user reliability, thereby ensuring ultrareliable communication for all users covered.

Conclusions
In this paper, we developed an analytical model to characterize stochastic queuing behaviors and the effects of fading characteristics on performance in the proposed CF mMIMO system under the κ-µ shadowed fading model. We derived approximate closedform expressions for the PDF, CDF and MGF of the sum of i.n.i.d. κ-µ shadowed RVs and improved their mathematical tractability. Then, based on the κ-µ shadowed fading model, we derived approximate closed-form expressions for the PDF, CDF and MGF of the PPSNR after the zero-forcing detector in the proposed CF mMIMO system. The κ-µ shadowed fading model can be extended to 6G wireless communication systems, and we verified the reasonableness and feasibility of its superior theoretical analysis performance with simulations and theoretical comparisons. In particular, we applied FBL information theory and SNC to model and characterize the arrival, service and departure processes, and we derived a closed-form expression for the UB-QDVP to analyze the effectiveness and reliability of the CF mMIMO system under the κ-µ shadowed fading model. Furthermore, by applying the theoretical analysis above, we also conducted a set of simulations to validate and evaluate the UB-QDVP gaps under different systems (the CF mMIMO, PD-NOMA and OMA systems) and fading environments and the relationships between the UB-QDVPs and EE under different fading environments. Making full use of the characteristics of the κ-µ shadowed fading model can improve the system EE, reduce the error probability and better satisfy the QoS requirements of 6G URLLC. Applying the κ-µ shadowed fading model to accurately describe small-scale fading and to model and analyze unmanned aerial vehicles, underwater acoustics and other wireless communication systems in the future is a promising direction.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.

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

Appendix B. Proof of Theorem 1
In reference [57], the authors derived the PDF of the sum of i.n.i.d. generalized Gamma RVs. Therefore, we transform the PDF of ηl k into the PDF of generalized Gamma RVs, and its expression is as follows. where Dˆl k = Γ(Kl k +1 Aˆl k ) Γ(Kl k ) . Forl ∈ L k and k ∈ [1, K], Aˆl k = 1 can be obtained.
Therefore, whenl ∈ L k , k ∈ [1, K] and Aˆl k is constant, the sum ofL i.n.i.d. generalized Gamma RVs can be approximated by the generalized Puiseux series expansion. The approximation for the PDF of the sum of L i.n.i.d. κ-µ shadowed RVs z is given by I(z; u 1 , u 2 , . . . , u L−1 ) (1−u 1 ) On the basis of the generalized Puiseux series expansion process, I(z; u 1 , u 2 , . . . , u L−1 ) can be approximated by the truncation order of any integer V, namely, When V → ∞, the series can accurately approach I(z; u 1 , u 2 , . . . , u L−1 ). Moreover, the series expansion converges very quickly. Equation (A2) can be expressed as To simplify the derivation, we transform the content of the cumulative sum into the PDF of the generalized Gamma distribution as follows.
Additionally, according to the nature of the Gamma distribution, the corresponding CDF of the approximation can be estimated as The coefficients c m and ω m can be obtained by using the moment-matching method [57]. Using this approach, the moments of the approximate PDF are equal to the exact moments.
For this purpose, the following equations are formulated, in which the approximate moments appear on the left and the exact moments appear on the right.
Theorem 1 is proven.