Evaluation of Cell-Free Millimeter-Wave Massive MIMO Systems Based on Site-Specific Ray Tracing Simulations

Ultra densification of the number of antennas combined with the use of large bandwidths in the millimeter wave (mmWave) spectrum is considered one of the main methodologies to achieve the quality of service requirements for future generations of wireless communications. Massive multiple-input multiple-output (mMIMO) cell-free (CF) systems have a large number of access points (APs) distributed in the coverage area, serving simultaneously a smaller number of users consuming the same time-frequency resources. In order to support realistic CF networks designs, this work proposes a performance analysis based on ray tracing simulations. The propagation modeling considers reflection, diffraction, diffuse scattering, atmospheric molecular absorption and foliage and rainfall losses. CF networks with APs equipped with multiple antennas operating in the 26 GHz, 38 GHz and 73 GHz bands are evaluated. From the simulation results, the communication channel is characterized and parameterized. Different performance analysis of CF networks are performed, based on downlink spectral efficiency. In addition, the performance of CF networks in rainy environments is evaluated, and it has been observed that this architecture promotes resistance to the effect of rain attenuation.


I. INTRODUCTION
Massive multiple-input multiple-output (mMIMO) systems, ultra dense networks and millimeter-wave (mmWave) communications are indicated as fundamental technologies to meet the quality of service (QoS) requirements of future wireless systems [1]. Using a large number of antennas to simultaneously serve a small number of users consuming the same time-frequency resources, mMIMO systems are The associate editor coordinating the review of this manuscript and approving it for publication was Olutayo O. Oyerinde . capable of providing high spectral and energy efficiencies using low-complexity signal processing [2]. The cell-free (CF) networks combine traditional mMIMO features, such as time-division duplexing (TDD) operation and simple linear precoding, in a non-centralized architecture in which a large number of distributed access points (APs) are connected to a central processing unit (CPU) to serve a smaller number of users [2]. From the user's perspective, this distributed architecture presents no cellular boundaries, and a homogeneous coverage is feasible. Several works show that cell-free networks are able to outperform colocated systems in terms of spectral and energy efficiency, while ensuring uniform QoS for most network users [3]- [6].
Communications in the mmWave range, from 30 GHz to 300 GHz, provide a large available bandwidth for broadband applications [7]. However, communications in this frequency range suffer from severe propagation losses in free space, in addition to atmospheric molecular absorption, and are also susceptible to intense rain loss [7]. Many works on short-range outdoor mmWave communications are focused on the 28 GHz, 38 GHz and 73 GHz bands [8]- [10]. Specifically in Brazil, the National Telecommunications Agency (Anatel) evaluates the use of the 26 GHz band for broadband applications in 5G networks [11]. Works on CF networks are mostly focused on sub-6 GHz communications [12]- [14]. However, the combination of the advantages of CF networks and mmWave communications has been investigated in [15]- [18]. In general, the published results depend on statistical channel models that, despite being able to specify several characteristics of the multipath behavior [19], [20], generalize certain properties of the communication channel, which can hide the particularities of specific environments.
Focusing on CF network designs in real environments, this work proposes a hybrid methodology, which combines the deterministic ray tracing method with a small-scale random fading model. The site-specific three-dimensional ray tracing model is based on the shooting-and-bouncing rays (SBR) method and is used to determine the channel large-scale parameters [21]- [23]. To characterize the propagation, the effects of reflection, diffraction and diffuse scattering are considered. In addition, the impacts of atmospheric molecular absorption, vegetation and rainfall losses are included.
In order to expand the results presented in [24], CF networks with APs equipped with multiple antennas, operating at 26 GHz, 38 GHz and 73 GHz, are evaluated by ray tracing simulations. The main objective of this work is to present design and performance evaluations of realistic CF networks using a channel model that incorporates the particularities of real propagation channels. Regarding the contributions of the research, the authors present performance analysis of CF networks based on ray tracing simulation results. To the best of the authors' knowledge, this is the first evaluation of cellfree networks, in which the APs are equipped with multiple antennas, applying the ray tracing method. In addition, an unprecedented analysis of the influence of rain attenuation on CF networks is provided.

A. PAPER ORGANIZATION
This paper is organized as follows. Section II presents the adopted CF system channel model, followed by a brief discussion of the methodology to allocate symbols in the TDD operation, in addition to the descriptions of the uplink training and downlink data transmission phases. Section III exposes the main characteristics of the ray tracing model applied in the simulations, including descriptions of the propagation phenomena, environment model and channel model. Section IV presents the numerical results, which include the characterization of the 26 GHz, 38 GHz and 73 GHz channels and the performance analysis of the CF networks and Section V presents the conclusions about the work.

B. NOTATION
Boldface uppercase and lowercase symbols denote matrices and column vectors, respectively. The Euclidean norm of a vector is denoted by · . The N × N identity matrix is denoted by I N . The superscripts (·) * , (·) T and (·) H stand for the conjugate, the transpose and conjugate-transpose, respectively. In turn, z ∼ CN (µ, ) denotes a circularly symmetric complex Gaussian vector with mean vector µ and covariance matrix . The expected value operator is written as E [·]. The Fourier transform F(y) of a function f (x) is described as F(y) = F{f (x)} x→y . Finally, 1 S (x) denotes the indicator function of an arbitrary set S, in which 1

II. CELL-FREE SYSTEM ANALYSIS
A. SYSTEM MODEL The CF mMIMO system considered in this work is illustrated in Fig. 1. The system architecture consists of M APs, each one equipped with L antennas, randomly located along the considered sites, which serve K single-antenna users applying the same time-frequency resources [25]. The channel gain vector between the m-th AP and the k-th user, denoted as g m,k ∈ C L×1 , is written as [13], [25] in which β m,k denotes the large-scale fading, that incorporates the path losses and shadowing effects,h m,k is the line-of-sight (LoS) components deterministic vector and h m,k ∼ CN (0, I L ) is the fading vector with components assumed mutually uncorrelated. The parameters β m,k , F m,k andh m,k , as well the links LoS conditions, are calculated using the ray tracing method, as exposed in Sections III-A and III-C. Considering (1), the channel gain vector can be represented by a sum of a deterministic vectorḡ m,k (result of ray tracing simulations) and a random vectorg m,k , respectively expressed byḡ resulting in a channel gain vector distributed as g m,k ∼ CN (ḡ m,k , α m,k I L ), in which α m,k = β m,k /(F m,k + 1) [25].

B. TDD OPERATION
Assuming that the channel gain is approximately invariant within the time and frequency coherence intervals, the TDD operation must be allocated in a block with a number of samples consistent with the imposed limitations. Each link is subject to different coherence times and bands, therefore, the TDD block must be allocated with a number of samples that maintains all transmissions within the coherence intervals. Thus, denoting the coherence band and the coherence time of the channel between the m-th AP and the k-th user by B m,k and T m,k , respectively; the number of samples in the TDD scheme is equal to τ c = min(B m,k × T m,k ). This assignment ensures that all transmissions on the network are made within the coherence intervals. By fixing the speed of a user's mobility, the channel coherence time is proportional to the wavelength, which reduces the value of τ c as the carrier frequency increases.

C. UPLINK TRAINING PHASE
During the training phase, all K users send pilot sequences of length τ p < τ c to all APs in the network. Denoting the pilot sequence transmitted by the k-th user as √ τ p ϕ k ∈ C τ p ×1 , in which ϕ k 2 = 1 and ϕ H k ϕ k ∈ {0, 1}, k = k, the received L × τ p pilot matrix at the m-th AP is expressed by [25] in which P p is the uplink transmitted power and N p;m ∼ CN (0, σ 2 AP I L ) is the noise matrix at the m-th AP. To estimate the channel, consider y m,k , defined as the projection of Y p;m onto ϕ k , expressed as which can be rewritten as y m,k =ȳ m,k +ỹ m,k , in which y m,k andỹ m,k are the deterministic and random parts of y m,k , respectively described bȳ and withñ p;m,k = N p;m ϕ k denoting the equivalent noise vector with distribution CN (0, σ 2 AP I L ). Based onỹ m,k = y m,k −ȳ m,k , each AP estimates its respective channel vectors adopting a Bayesian minimum mean square error (MMSE) estimation defined by [13] and it is possible to demonstrate thatĝ m,k =ḡ m,k + c m,kỹm,k , in which c m,k is expressed as In (9), P k is a subset of indexes of the respective users who reuse the k-th pilot sequence, i.e., P k = {k ∈ {1, · · · , K }|ϕ H k ϕ k = 1}. Therefore, the estimated channel vectors areĝ m,k ∼ CN (ḡ m,k , γ m,k I L ) distributed and the estimation error, defined as m,k g m,k −ĝ m,k , is distributed If τ p ≥ K , it is possible to associate K mutually orthogonal pilot sequences for the users and, in this case, the pilot contamination in (5) is null, given that ϕ H k ϕ k = 0, ∀ k = k ∈ {1, · · · , K } or, equivalently, P k = {k}. However, the value τ p is limited by the coherence interval, which can make the condition τ p ≥ K unfeasible, making it necessary to reuse pilot sequences between different users.

D. DOWNLINK DATA TRANSMISSION AND PERFORMANCE ANALYSIS
In the downlink data transmission phase, the M APs apply the so-estimated channel vectors to precode the intended data to each user. Adopting the conjugate beamforming, the transmitted signal by the m-th AP is given by [25] in which P d is the downlink transmitted power, η m,k are the power control coefficients and q k ∼ CN (0, 1) is the data symbol intended to the k-th user [14], [25]. The simplest power assignment strategy is adopted, in which the APs uniformly distribute the power to all users satisfying the condition E x m 2 = P d , resulting in Assuming that the channel is reciprocal, i.e., the downlink and uplink channel coefficients are numerically equivalent, and that each user terminal only has knowledge of its own statistics, the received signal at the k-th user is expressed by [26] in which n u;k is the CN (0, σ 2 u ) additive noise at the k-th user, DS k is the desired signal, BU k is the beamforming uncertainty and UI k,k is the interference caused by the k -th user in the k-th user link. These quantities are respectively described by and Treating the sum n eff;k = BU k ·q k + K k =k UI k,k ·q k +n u;k as the effective noise and assuming that q k is independent and uncorrelated of DS k , BU k and UI k,k , the signal-to-noise ratio (SNR) experienced by the k-th user can be expressed as in (17), as shown at the bottom of the page. This expression is proved in Appendix A. In this equation, m,k,k = γ m,k Thus, the SE experienced by the k-th user, denoted as S k , is found replacing (17) in [25]

III. RAY TRACING MODELING A. RAY TRACING PROPAGATION MODEL
The characterization of electromagnetic propagation by the ray tracing SBR method follows three main steps: transmission, tracing and reception [21], [27]. In the transmission step, multiple rays are launched uniformly along the angular sphere. In the tracing stage, the rays' lengths are incremented and, at each incremental step, the algorithm evaluates if the ray is in the vicinity of an object [28]. In these cases, an intercept test is performed and, if so, the ray is reflected according to Snell's law. New sources of rays can be generated in these intercepts, which are applied to characterize the phenomena of diffraction and diffuse scattering. Finally, in the reception stage, a sphere is defined around the observation points (receiving points) and if a ray intercepts this sphere, its contribution is considered in the link. The diameter of the reception sphere is determined in order to minimize double receptions of the same wavefront.
In this work, the diffraction is characterized based on the uniform theory of diffraction (UTD) [29]- [32], while the diffuse scattering is modeled according to the Beckmann-Kirchhoff theory [33]- [35]. In the simulations, at most, second-order diffractions are considered, following the limitations imposed by the UTD in the transition regions [32]. In turn, in the characterization of diffuse scattering on rough surfaces, the methodology of subdivision of the illuminated regions into tiles that have areas that correspond to the far-field condition is adopted [36], [37].
Thus, an arbitrary received ray contributes in the link with a normalized power given by in which λ is the wavelength, d is the total path length and ρ p (θ a , φ a ) is the polarization mismatch loss between the receiver antenna polarization and the received electric field, which has angles of arrival in azimuth and in elevation θ a and φ a , respectively [38]. In turn, the parameters i and ρ s;i are the reflection losses, based on Fresnel coefficients, and the loss factor in the specular direction due to diffuse scattering on a surface with effective roughness σ h , respectively. On the other hand, ζ k is the diffuse scattering loss in an arbitrary direction, as described by the Beckmann-Kirchhoff theory [33], [34]. Finally, the value D j in (19) is the diffraction loss, calculated according to the UTD [29], [30]. The terms VOLUME 10, 2022 A re;i , A ds;k and A d;j express the spreading factors related to reflection, diffuse scattering and diffraction, respectively. In (19), the values ρ m , ρ r and ρ f are three absorption loss factors (expressed in dB): frequency dependent atmospheric molecular absorption (ρ m ), rain attenuation (ρ r ) and foliage losses (ρ f ). This work considers the molecular losses caused by dry air and water vapor absorptions, that are calculated according to the numerical model presented in the Recommendation ITU-R P.676-10 [39]. In turn, the rain attenuation, expressed in dB, is calculated as function of the rain precipitation rate C r (mm/hour) by the following expression [40], [41] in which (c 1 , c 2 ) are parameters that depend on the frequency, wave polarization and incidence direction to the rainfall and d eff is the effective distance in the rainy environment. Finally, the foliage losses ρ f are calculated by applying the radiative energy transfer (RET) theory [42], [43]. A vegetation medium is geometrically modeled as a polyhedron and ρ f depends on the foliage depth, which is the length of a ray inside this medium. Due to the lack of data on the scattering characteristics of the local vegetation, in the simulations, the parameters presented in [43] were applied.

B. ENVIRONMENT CHARACTERIZATION
The propagation environment considered in the simulations is composed of objects, modeled as polyhedrons, which represent houses, buildings and trees. This environment is based on coastal neighborhoods in the city of João Pessoa, capital of the state of Paraíba, Brazil. A view of the propagation environment model is shown in Fig. 2.
The calculated average height of buildings and housing density are 7,9 m and 26.8%, respectively. Due the lack of information about the heights of the treetops, it was considered that all trees are medium sized with random top heights between 5 m and 10 m. The objects in the environment are classified as house, building or vegetation, and each class has specific parameters. Objects representing vegetation are predefined and have specific scattering properties, as described in [43]. On the   other hand, objects with heights less than 7 m are classified as houses, otherwise, are considered buildings. In the case of houses, the walls and ceilings planes materials are brick and roofing tiles, respectively. For buildings, all planes are modeled as concrete surfaces. It was assumed that the ground plane is ideally smooth (σ h = 0 mm) and composed of asphalt throughout the environment. Table 1 presents the values of electrical permittivity (ε r ) and effective surface roughness (σ h ) of each considered material [44]- [47].
In the simulations, three square sites of area L s × L s were chosen. The respective central coordinates of each site are listed in Table 2. To define the positioning of users, the entire propagation environment is subdivided into N mr square microregions of area L mr × L mr , in which users are positioned at their centers with receiving equipment at a height h rx . Therefore, the reception points define a reception grid, with each element indexed by the pair (i, j), as illustrated in Fig. 3. Each user is associated with an instantaneous velocity vector v i,j that has a random direction and norm v i,j = v u , in which v u is the user scalar velocity. It is assumed that the velocity vectors do not change their direction during the reception of all multipath components and that there is no significant position variation within a coherence time.
Given that only outdoor environments are simulated, reception points contained in the limits of an object are excluded. Therefore, a site with N mr microregions has N v valid positions, in which N v < N mr . The center coordinates of the M APs and K users are randomly associated to M × K N v valid coordinates. To increase the amount of data generated per simulation, several user groups (UG) were considered [24]. Each UG is composed by K users and there is no dependency between members of different UGs. The total number of UGs in each simulation is calculated by Each reception point is mapped to a user in its respective UG. That is, there is a one-to-one relationship between the indexes (i, j) and (k, u), in which k ∈ {1, · · · , K } is the user index and u ∈ {1, · · · , N UG } is the UG index. To simplify the notation and avoid misunderstanding, in the channel and CF systems analysis, the index notation (i, j) was reduced to k (related to a specific UG).

C. RAY TRACING CHANNEL MODEL
The ray tracing simulations consider M APs, each one equipped with L antennas, transmitting rays that are captured by the receiving points in the reception grid composed of N UG UGs with K users each. To avoid excessive computational processing, it is assumed that the channels relative to the L antennas of the same AP share the same largescale parameters, differing only in the small-scale fading components. From this assumption, it is not necessary to transmit rays from the L coordinates of the antennas for the same AP. Thus, it is sufficient to apply the ray tracing only from the centroid of each AP antenna array. Therefore, consider an uniform linear antenna array, as presented in Fig. 4, where the position of the l-th antenna element in the m-th AP is given by in which c m is the array centroid position vector, t m is the array orientation vector ( t m = 1) and d ant is the minimum spacing between the antenna elements. The ray tracing simulations are performed transmitting rays only from the centroid coordinates c m . Assuming that all antenna elements have LoS conditions if there is visibility between the centroid and a given user and that the separation distance between the AP and the user is much greater than the length of the array, the elements of the LoS vectorh m,k are calculated by in whichφ m,k is the angle between t m and the received ray.
Assuming that the received signal in the link between the m-th AP with the k-th user is composed by the superposition of N r;m,k rays, the double-directional channel impulse response (CIR) is expressed by [7] H m,k (t, τ, , ) (23) in which t, τ , and are the absolute time, delay, angle of arrival and angle of departure domains; and in turn, p r;m,k,n and τ m,k,n are the normalized power, as expressed in (19), and delay of the n-th received ray, respectively. In other hand, m,k,n = [θ a;m,k,n , φ a;m,k,n ] T and m,k,n = [θ d;m,k,n , φ d;m,k,n ] T are the angle of arrival and angle of departure vectors, respectively. The terms θ d;m,k,n , θ a;m,k,n are azimuth angles and the terms φ d;m,k,n , φ a;m,k,n are elevation angles. Following in (23), ϕ m,k,n (t) is the phase component of the n-th received ray, expressed by [48] ϕ m,k,n (t) = 2π[(f + ν m,k,n )τ m,k,n − ν m,k,n t], (24) in which f is the carrier frequency and ν m,k,n is the Doppler shift, caused by the incidence of a ray with arrival direction vector denoted by k m,k,n , which is calculated by All parameters involved in the CIR are computed using the ray tracing method. Applying the CIR expressed in (23), the channel autocorrelation function, denoted as R m,k (ς, τ, , ), is calculated by a time average, in which where ς denotes an absolute time shift. Therefore, the joint power spectrum of the channel P m,k (ν, τ, , ), with ν VOLUME 10, 2022 denoting the Doppler shift domain, is calculated through the Fourier transform of R m,k (ς, τ, , ), resulting in ,k,n ). (27) The marginal power spectra can be found by integrating (27) in the remaining domains and respective supports. From (27), the large-scale fading is calculated by The distance-dependent path losses, expressed in dB, are parameterized by applying the close-in free-space (CI) reference model, described by [7] PL CI (f , d tx ; κ) = PL FS (f , d 0 ) + 10κ log 10 (d tx ), (29) in which d tx > d 0 is an arbitrary distance between the transmitter and the receiver, κ is the propagation exponent and PL FS (f , d 0 ) is the free-space path loss for a reference distance d 0 , which is considered d 0 = 1 m. The parameter κ is calculated by a linear least squares regression. Using the data obtained in simulation and the model in (29), the shadow fading (expressed in dB) is calculated by with d m,k denoting the distance between the m-th AP and the k-th user. In channel statistical models, shadow fading is often modeled by a zero-mean Gaussian random variable with standard deviation σ χ . However, in this work, the shadow fading of each individual link is calculated by means of (30) and is parameterized by its sample standard deviation. As evidenced in the channel model described in (1), the Rice factor F m,k is calculated as the ratio between the LoS component power, denoted asp m,k , and the remaining scattered power, i.e., In urban environments, the probability of a given user being in a LoS region is inversely proportional to d tx . This probability can be parameterized as [7] p LoS (d tx ) = min in which d 1 and d 2 are numerically fit parameters that can be found by minimizing the mean squared error between the simulation data and the theoretical expression. In the considered ray tracing simulation, the LoS condition of a given link is defined geometrically and do not depend on the carrier frequency.
To obtain the temporal and spectral dispersions, the RMS delay spread (RMSDS) σ τ ;m,k and RMS Doppler spread (RMSDoS) σ ν;m,k are respectively calculated as follows: σ τ ;m,k = 1 β m,k N r;m,k n=1 p r;m,k,n (τ m,k,n − µ τ ;m,k ) 2 , (33) σ ν;m,k = 1 β m,k N r;m,k n=1 p r;m,k,n (ν m,k,n − µ ν;m,k ) 2 , (34) in which µ τ ;m,k and µ ν;m,k are the mean delay and the mean Doppler shift in the respective link. As known, the RMSDS and RMSDoS are inversely proportional to the channel coherence band B m,k and coherence time T m,k [49], respectively. Analytical relationships between coherence parameters and spread parameters do not exist, therefore, the following approximations are considered: B m,k ≈ 1/50σ τ ;m,k and T m,k ≈ 1/σ ν;m,k [49]. These coherence parameters are used to define the length of the TDD block, as discussed in Section II-B.
Finally, to measure the dispersion levels in the angular domains, four different spread parameters can be calculated: RMS azimuth spread of departure (RMSASD) σ θ ;d , RMS azimuth spread of arrival (RMSASA) σ θ ;a , RMS elevation spread of departure (RMSESD) σ φ;d and RMS elevation spread of arrival (RMSESA) σ φ;a . These four parameters are calculated by the general expression 1 [50]    , (35) in which ϑ m,k,n ∈ {θ d;m,k,n , θ a;m,k,n , φ d;m,k,n , φ a;m,k,n } is an azimuth or elevation angle, of departure or arrival.

IV. NUMERICAL RESULTS
In this section, the results of channel characterization and CF networks performance obtained from the ray tracing simulations are presented. The simulations are performed on 26 GHz, 38 GHz and 73 GHz channels, considering square sites with an area equal to 400 m×400 m (L s = 400 m). The interference from neighboring sites are not considered, i.e., each simulation is performed separately. Each user equipment (UE) is positioned at a height of h rx = 1.5 m, while each AP centroid is at a height h AP = 15 m. A low mobility scenario is considered, in which users are modeled as pedestrians with a scalar velocity equal to v u = 1.2 m/s. In the simulations, the angular difference between adjacent rays was approximately 0.135 • . Each ray path is limited to five reflections and the CIRs are composed of a maximum 1 Considering the values of β m,k as expressed in (28) and denoting the absolute value in (35) by I = 1 β m,k N r;m,k n=1 p r;m,k,n exp(jϑ m,k,n ) , the triangular inequality for complex numbers guarantees that I ≤ 1 and, consequently, −2 log(I ) ≥ 0. Therefore, the angular spread as expressed in (35) is necessarily a real number (σ An;m,k ∈ R). of 50 received rays. A total of N I = 10 interactions were considered in each simulation. The evaluated CF system is composed by a network of at most M = 30 random located APs, representing a scenario with an AP density of approximately 188 APs/km 2 , serving a maximum of K = 10 users. Each AP is equipped with a vertically aligned uniform linear antenna array with elements spacing of d ant = λ/2 and transmits vertically polarized rays with an average power P d = 200 mW. In turn, users transmit pilots with power equal to P p = 100 mW. Without loss of generality, it is assumed that both APs and users equipment are affected by additive white Gaussian noises (AWGN) with a power spectral density of -174 dBm/Hz. Furthermore, the transmitters and receivers have noise figures of 6 dB. The parameters used to configure the simulations are summarized in Table 3.

A. CHANNEL MODELING RESULTS
In this section, the channel modeling results are presented, which are obtained from the ray tracing simulations at the three different sites. Large-scale parameters statistics were calculated using these simulation results, with channels being distinguished by different carrier frequencies (26 GHz,38 GHz and 73 GHz). These channel parameters and statistics are summarized in Table 4.
In LoS links, the path losses are parameterized with propagation exponent κ = 2 with a shadow fading standard deviation of approximately σ χ = 2.0 dB at all considered carriers. In turn, the frequency that imposes higher path losses on NLoS links is 73 GHz, with a propagation exponent at κ = 3.8, attached to a free space loss of 69.7 dB, combined to shadow fading, which has standard deviation σ χ = 12.2 dB. In Fig. 5 are presented the path loss results, associated with their respective CI models, in addition to the NLoS links shadow fading histogram, superimposed by the respective Gaussian probability density function, obtained from the simulations results at 73 GHz.
As explained, the probability that a user is inside an LoS region is calculated based on the link geometry, thus it does not depend on the carrier frequency. It has been observed that the parameters that best fit the model in (32), for the   Table 4, are used to calculate the channel coherence parameters and thus define the value of τ c . As expected, due to the inherent decrease in coherence time, the value of τ c decreases with the respective increase in frequency, with a minimum value of 100 samples in 73 GHz channels at Site 2. The maximum observed value of τ c was obtained from the simulations at 26 GHz channels at Site 2, adding up to 367 samples.
Analyzing the RMSASA and RMSASD statistics presented in Table 4, it can be seen that their sample averages in LoS and NLoS links differ by at most 24.1 • . Furthermore, on average, the signal power is more scattered in the angleof-arrival domains, with maximum value of 33.3 • observed in NLoS links at 38 GHz. Fig. 7 shows examples of power spectra results in the azimuthal domains at 26 GHz, resulting in a RMSASA of σ θ;a = 28.6 • and a RMSASD of σ θ;d = 24.5 • . On the other hand, the RMSESD results show that, on average, the CIRs are composed by groups of transmitted rays with a relatively low elevation angle variation, with maximum mean value of 1.5 • on 73 GHz NLoS channels. Similar to azimuthal scattering, the average values of the  RMSESA are higher than the respective average values of the RMSESD, with a maximum value equal to σ φ;a = 6.7 • in NLoS channels with carrier at 38 GHz. Table 5 presents the average SE values, as a function of L. These results were obtained considering M = 30, K = 10 and τ p = 5 samples. In these conditions, it can be seen that the average SE value is relatively stable with the frequency variation, having only small decreases with increasing frequency. Among the considered values, there is an increase of approximately 0.8 (bits/s/Hz) in the average SE value when L is duplicated. Fig. 8 shows the curves of the empirical cumulative probability functions (CDFs) of the SE for the respective carriers, considering M = 30, K = 10, τ p = 5 samples and L ∈ {8, 16, 32}. The small variation of the SE values shown in Fig. 8 is a result of the relative stability of path losses, shadow fading and Rice factor in the considered carriers, as shown in Table 4.  Due to the dispersed spatial distribution of APs in the network, it is possible that a set of antennas does not effectively participate in the communications. To quantify this effect, the number of effective APs (M eff ) is calculated. The aforementioned quantity is defined as the minimum number of APs that contribute at least 95% of the link power [4]. Table 6 shows the statistics of M eff for the different carriers considering τ p ∈ {5, 20}. Initially, it is observed that the increase of τ p raises the average value of M eff , in which, for 26 GHz channels, the increasing from τ p = 5 samples to  The expression in (17) tends to be independent of downlink transmission power in high power regimes (P d → ∞). To quantify the transmission power required for the average SES to reach its limit, the saturation powerP d is defined as the value that satisfies

B. CELL-FREE PERFORMANCE RESULTS
Curves are presented in Fig. 9 for the average SE as function of   SE increases with the number of transmitting antennas in the network. From this result, some scenarios can be analyzed. For example, on 26 GHz channels, to achieve an average SE of 3 bits/s/Hz, the CF network requires 25 APs with eight antennas, nine APs with 16 antennas or five APs with 32 antennas. Networks composed by 13 APs with 32 antennas or 22 APs equipped with 16 antennas provides the same average SE value at 4 bits/s/Hz. Fig. 11 presents the average SE curves as a function of the number of users K . As expected, the average SE decreases with the increase of K due to the growth of interference. From these results, it can be seen that the number of users served by the CF network limits the average SE. For example, to achieve an average SE higher than 5 bits/s/Hz, networks with 30 APs can only simultaneously serve a maximum of four users when L = 8 and six users when L = 16.
As a last result, the influence of the rain attenuation on the average SE is evaluated. As explained in Section III-A, rain scenarios cause additional losses to the link, influencing the SNR observed by a given user. For example, in heavy rain scenarios, with precipitation rate C r = 50 mm/hour, the specific attenuations in the 26 GHz, 38 GHz and 73 GHz bands are 7.7 dB/km, 11.8 dB/km and 17.4 dB/km, respectively [40]. Denoting the average SE as a function of the rainfall precipitation rate C r byS(C r ), the rain influence on the average SE is measured by whereS (0) is the average SE in a scenario without rain (C r = 0 mm/hour). The curves of this measure are presented in Fig. 12, in which the factor C r is varied from light rain (C r < 2.5 mm/hour) to heavy rain scenarios (C r > 50 mm/hour). Initially, it can be observed that there is a range of rainfall rate, between 0 mm/hour and a constantČ r , in which the average SE slightly increases. In this range, rain attenuation reduces the contribution of interference from distant APs and consequently increases the relative contribution of nearby APs, increasing the SNR. However, when C r >Č r , the rain attenuation is sufficient to reduce the SNR even in the group of effective APs close to a given user. The constanť C r depends on the frequency and number of antennas on each AP. TheČ r values corresponding to the CF networks operating at 26 GHz extrapolate the limits of the simulations, that is, on this carrier,Č r > 100 mm/hour. On channels operating at 38 GHz,Č r = 48 mm/hour (L = 8) anď C r = 64 mm/hour (L = 32) were calculated.
The effect of the rain attenuation is more severe in networks operating at 73 GHz, withČ r = 8 mm/hour andČ r = 16 mm/hour. Despite this, it can be observed that the rain influence on the average SE is relatively low even in very intense rainfall scenarios, with a value of -0.3 dB (73 GHz and L = 8) of loss in a rainfall with rate of 100 mm/hour. Using the meteorological analysis method presented in [51], it is estimated that, in the coastal region of João Pessoa, a rainfall with rate higher than 100 mm/hour lasting for, at most, one minute has a probability of occurrence of 0.01% (1 in 10,000 cases). Therefore, the occurrence of precipitation capable of drastically reducing the average SE of the analyzed CF networks is unlikely.

V. CONCLUSION
This work presents the analysis of CF mMIMO networks, considering APs equipped with multiple antennas communicating in the mmWave spectrum, based on site-specific ray tracing simulations in order to guide real and applicable system designs. The communication channel in the 26 GHz, 38 GHz and 73 GHz bands, for three different sites, was characterized and parameterized. According to the simulation results, the average SE grows by approximately 0.8 bit/s/Hz when the number of antennas on the APs is doubled. When APs are equipped with 32 antennas, the CF network can provide an average SE higher than 5 bits/s/Hz. In high power regime, the average SE can exceed 15 bits/s/Hz. On the other hand, it was observed that, at most, 56.6% of the AP network effectively contribute to communications. In conclusion, the CF networks show some resilience to the effect of attenuation by rain. In the most restricted case (heavy rain with carrier at 73 GHz), a maximum decrease of -0.3 dB in the average SE was observed. Furthermore, according to the results, the increase in the number of antennas in the APs makes the network more resistant to the influence of rain.

APPENDIX A PROOF OF THE EQUATION (17)
The expression in (17) is obtained from [4], [12], [52] The derivation of the spectral efficiency expression depends on the calculation of the terms DS k , E |BU k | 2 and E UI k,k 2 .

A. CALCULATION OF DS k
Assuming that the MMSE estimation error vector m,k is independent of the estimated channel vectorĝ m,k and knowing that E m,k = 0, the component of the desired signal is expressed by √ η m,k η n,k E ξ m,k ξ * n,k . (40) Given that the variables ξ m,k and ξ n,k are mutually independent for m = n, it is possible rewrite (40) as The quantity E ξ m,k ξ * m,k in (41) Replacing (42), (43) and (44) in (41) To solve the variance term in (46), the auxiliary random variablesμ m,k,k andμ m,k,k are defined as in which E μ m,k,k = 0, E μ m,k,k = 0 and E μ m,k,k μ * m,k,k = 0. Applying the auxiliary random variables is found that In turn, the term E μ m,k,k 2 in (50) (51), (52) and (54) in (50), it is obtained η m,k α m,k (γ m,k F m,k + γ m,k + α m,k F m,k ). (55) Finally, from the substitution of (47) and (55) into (46), it is found that η m,k α m,k (γ m,k F m,k + γ m,k + α m,k F m,k ). (56)