An Efficient Plane-Waves Superposition Method for Improved Spatial Correlation in Simulated Reverberation Chambers

Accurate Rayleigh propagation environments within an ideal reverberation chamber (RC) are synthesized computationally using a superposition of plane-waves (PWs). Randomness in the electro-magnetic (EM) field distribution is achieved by synthesizing many field realizations featuring slant-polarized PWs with randomly defined transverse field components. These PWs propagate along a large number of predefined and uniformly distributed directions determined only once via an efficient spiraling sampling scheme over the unit sphere. Comparisons with theoretical Rayleigh field statistics indicate that the proposed method yields random EM field distributions that in several cases follow more closely the field spatial correlation properties, for comparable numerical burden, than other available numerical techniques. Finally, a formula relating the number of PWs to the size of a desired RC analysis volume was provided.


I. INTRODUCTION
Reverberation chambers (RCs) are widely employed as a standard test equipment in electromagnetic compatibility (EMC) assessments [1]. RCs present several advantages compared to other EMC testing facilities because nearly ideal random field statistics can be synthesized through so-called mode stirring, achieved using mechanical tuners, spatial diversity, or frequency stirring [2]. RCs random field properties make them appealing not only for electronic equipment immunity and emission tests, but also to study potential biological effects of exposures to electromagnetic (EM) fields [3], [4], [5], [6].
Different approaches have been recently employed to simulate the EM environment within a low-loss, highly reverberant chamber or cavity. Deterministic approaches, involving the simulation of EM field distributions inside RCs, have employed Ray-Tracing [7] or image theory [8] techniques, the circuital transmission line matrix (TLM) method [9], as well as full-wave techniques such as the Method of Moments (MoM) [10], the finite-difference time-domain (FDTD) technique [11], and the finite element method (FEM) [12]. On the other hand, statistical approaches based on the Monte Carlo (MC) superposition of plane-waves (PWs) have been proposed in [13], [14], [15], [16], [17], and [18] to synthesize Rayleigh random field distributions within an AV. In general, full-wave deterministic models result in extremely large VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ numbers of unknowns to model RCs walls (integral equation methods) or the volume within the walls (differential equation methods), while MC-based methods appear better suited to reproduce the desired statistical field characteristics within the AV with a modest computational burden. Several MC-based approaches have been proposed. In [13] and [14], PWs with randomly determined propagation directions, polarization, and phase, but constant magnitudes, were assumed, whereas in [17] and [18], the PWs magnitudes were randomized as well. Rather than randomizing the PWs impinging directions, a different approach proposed in [15] and [16] invoked spectrum sampling theory to derive a predefined set of directions, randomizing only the complex amplitudes of elliptically-polarized PWs. In this way, computational efficiency was notably improved since the PWs directions were common to each random field realization, hence the EUT response to each individual PW had to be determined only once.
In this paper, the Rayleigh field synthesis in RCs is also carried out using PWs superposition based on a predefined set of directions. However, compared with West et al. [15], a different parametric sampling method over the unit sphere is employed to derive the set of PWs propagation directions. Moreover, unitary-magnitude slanted-polarized PW is associated to each impinging PW direction. Comparisons with the theoretical spatial correlation properties characterizing an ideal Rayleigh environment indicate that the random fields generated through the herein proposed method reproduce them well over larger distances, resulting in roughly twice as large AVs for a comparable computational workload. Consequently, the proposed method may be usefully adopted in the design of state-of-the-art bioeffects studies, where the AV must be maximized to enable the analysis of concurrent exposures of large animal cohorts, such as those of rodents.

II. MODELS AND METHODS
The field statistics within an ideal RC can be rigorously modeled through PWs superposition, as described in [19]. Some key findings are briefly summarized hereinafter. Upon assuming and suppressing a time-harmonic dependence e jωt , ω being the radian frequency, the E-field within an RC is expressed through its uniform PW spectrum integral as where k and r are the wave number and the field point vector, respectively. Within the unit sphere illustrated in Fig. 1, the former can be related to impinging PWs elevation and azimuth angles (θ, φ) as k = k(xsinθcosφ +ŷsinθsinφ + zcosθ), where k = 2π/λ, λ being the free-space wavelength. Finally, denoting withθ andφ the spherical reference frame unit-vectors, the angular spectrum F in (1) is given by The complex amplitudes F θ and F φ are completely defined by the statistically independent distributions of their respective real and imaginary parts [20]. Hill's theory [21] invokes the Central Limit Theorem to show that Rayleigh fields Cartesian components are normally distributed, and to derive the closed-form field spatial correlation properties summarized in II-C.
Once a suitable spherical sampling scheme is implemented to closely approximate uniformly distributed PWs propagation directions, as detailed in the following, the random Rayleigh field computational generation becomes straightforward as Eq. (1) naturally embodies PWs superposition.

A. DEFINITION OF THE PWs IMPINGING DIRECTIONS
Exact uniform spherical sampling with an arbitrary number of points is very cumbersome, hence impractical in many applications. Indeed, several methods have been proposed to approximate an ideal uniform distribution [22], [23], [24], [25] and obtain a set of D points, corresponding to D propagation directions of the impinging PWs. The choice of the method may depend on the criterion adopted to quantify the sampling uniformity. For instance, in [15], PWs propagating directions coincided with the zeros of the L − th degree Legendre polynomial in elevation (θ) and over 2L uniformly spaced points in azimuth (φ) [22], [23], leading to D = 2L 2 impinging directions across a regular grid, as shown in Fig. 1. On the other hand, sampling schemes along spherical spirals have been implemented in [24] and [25] to yield less regular yet visibly more uniformly distributed samples on the sphere, as clearly illustrated in Fig. 1.
In this work, the angles defining the PWs propagation directions were derived from a uniform spherical spiral sampling, as suggested in [24] and [25]. Deviations from those schemes can be found in the distribution of points on the spiral that, while still uniform, includes the opposite poles, as shown in Fig. 1(a), as well as in the whole number of turns m ∈ N, resulting in the spherical spiral being expressed as φ = 2mθ, r = 1. As noted in [25], spiral schemes do not satisfy the antipodal symmetry desired in some applications. Nevertheless, the above choices, together with uniform sampling along the spherical spiral, lead to certain symmetries, in particular with respect to the following transformation: The proposed spherical spiral sampling method hinges upon determining the distance δ between neighboring points along the curve. Good sampling uniformity in both φ and θ is achieved when δ approximates the elevation gap between adjacent turns (π/m), simultaneously enforcing uniform segmentation along the curve path. With this approach, however, the total number of sampling points (i.e., PWs directions D) grows approximately as m 2 , not allowing enough granularity. Instead, to choose any desired number of sampling points, the initial estimate of δ can be rather linked to a seed number P by computing δ = π/ √ (P − 1)/2, where now m = √ (P − 1)/2 . In this way, a desired number of sampling points D can be attained for a certain P, which is found starting from the initial guess P = πD/2 , typically with just a single iteration.
The spiral sampling approach leads to a greater uniformity of the sphere surface sampling and provides flexibility in choosing an arbitrary number of sample points, unlike in [15], where the number of sampled points is limited to 2L 2 , with L an integer number. Finally, the spiral sampling approach provides almost uniform projections of the sample points along all three main cardinal directions, in particular along the z-axis, as shown in Fig. 1(b). This figure reports histograms of the Cartesian coordinates for a set of D = 200 points sampled with both methods.

B. PLANE-WAVE SUPERPOSITION ALGORITHM
The proposed PWs superposition algorithm to generate stochastic Rayleigh field realizations is described by the flowchart in Fig. 2. Upon defining the main input parameters, such as the desired number of field realizations (R) and PWs (N = 2D following the formalism in [15], where each field component PW counts individually), the aforementioned spherical sampling algorithm is employed to define the set of PWs propagation directions to be employed in generating the field realizations. Hence, the generic slant-polarized PW electric field vector is expressed as: where i = 1, . . . , D. As the PWs propagation directions are determined once, randomness is introduced through the definition of its field components, expressed as follows where k (i) ,θ (i) andφ (i) are mutually orthogonal. Corresponding random slanted polarization states are attained upon enforcing in-phase θ and φ components as follows where M φ and α (i) are random real variables. While α (i) is simply obtained as an independent, uniformly distributed stochastic variable in [0, 2π), the remaining variables are defined within the innermost cycle of the algorithm (see Fig. 2 so that each slant-polarized PW features unit field magnitude. This step is repeated for each of the D slant-polarized PWs, which are then superposed at all desired points within the AV (last flowchart block in Fig. 2), thereby yielding a unique Rayleigh field realization. Repeating R times this procedure leads to a set of Rayleigh field realizations, suitable to carry out the following ensemble statistics.

C. IDEAL RAYLEIGH-FIELD SPATIAL STATISTICS
To quantify the goodness of the proposed PWs superposition algorithm, theoretical Rayleigh field spatial statistics are illustrated and later used as a measure of the synthesized field fidelity to an ideal RC. In particular, the spatial autocorrelation of the electric field vector and those of its components have been tracked. The former features the following expression embodying a noteworthy spatial isotropy resulting from its behavior being independent of direction [20], [21]: where d = |r 1 − r 2 |, with r 1 and r 2 being arbitrary observation points within the AV. However, the same does not hold true for the field components [20], [21]. For example, the spatial autocorrelations of their real or the imaginary part of the z-component depends on whether the field points r 1 , r 2 lie in the xy-plane or along the z-direction. Focusing on the real part, without loss of generality, in the xy-plane we have ρ Re(E z,xy ) (r 1 , r 2 ) = Re(E z (r 1 ))Re(E z (r 2 )) Re(E z (r 1 )) 2 Re(E z (r 2 )) 2 where r 1 −r 2 =ξ d, withξ ·ẑ = 0, while along the z-direction ρ Re(E z,z ) (r 1 , r 2 ) = Re(E z (r 1 ))Re(E z (r 2 )) Re(E z (r 1 )) 2 Re(E z (r 2 )) 2 where now r 1 − r 2 =ẑd. For the sake of brevity, the autocorrelation functions in Eqs. (11) and (12) will be henceforth referred to as ρ E z,xy and ρ E z,z , respectively. Finally, another key statistical property of ideal Rayleigh fields, stemming from the real and imaginary parts of the corresponding Cartesian field components being normally distributed, features the chi-squared probability distribution function (PDF) with 6 degrees of freedom (χ 2 6 ) [19] for the squared E-field magnitude (|E(r)| 2 ).

III. RESULTS
In the following, the theoretical spatial autocorrelation function in Eqs. (10)- (12), as well as the aforementioned χ 2 6 PDF, are adopted as reference goodness metrics for the quality of the realized Rayleigh field distributions. These metrics are compared with the reference method presented in [15], which represent the relevant benchmark to further assess the accuracy and computational burden of the proposed approach.

A. METHOD COMPARISON
A large number (R = 5000) of Rayleigh field realizations, involving the superposition of N = 3600 PWs (L = 30), were used to compute the |E(r)| 2 PDF for both the proposed and reference method. Fig. 3 illustrates how the respective normalized PDFs fit quite well the expected χ 2 6 distribution. The application of Kolmogorov-Smirnov (KS) [26] and Anderson-Darling (AD) [27] goodness-of-fit (GOF) tests indicated that in all the evaluated cases the null hypothesis (i.e. that the samples come from the specified distribution) fails to be rejected (with a significance level α = 0.05), thus both samples can be reasonably assumed to come from a χ 2 6 distribution. Moreover, the statistics from the two GOF tests computed for both data samples showed that the difference in terms of how well they fit the χ 2 6 distribution can be considered to be not statistically significant. More precisely, the KS statistics computed on the analysed samples corresponds to 0.0067 and 0.0084 for the reference and proposed method, respectively, while the AD statistics resulted in 0.2066 and 0.3160, respectively.
As done in [15], spatial autocorrelations were derived numerically upon choosing a reference starting point r 1 and considering a set of applicable points r 2 within the computational domain using the same number (R = 5000) of random realizations but a smaller number (N = 400 , L = 10) of PWs than in the previous comparison so as to ''stress'' both methods already within a few wavelengths. Specifically, the AV spanned a cube with 25λ long edges, with all distances used for correlation estimates taken from a corner point towards random directions in 3D space for ρ E and on the  xy-plane for ρ E z,xy , whereas uniformly spaced sampling points along the z-axis were selected for ρ E z,z . An example of the 3D distribution of evaluation points is given in Fig. 4, while the autocorrelation results are presented in Fig. 5. The numerically estimated correlation functions for ρ E z,z , based on field points r 1 , r 2 laying on a line, behave smoothly (as seen in Fig. 5 (c)). However, those for ρ E and ρ E z,xy exhibit rather noisy behaviors starting just a few wavelengths away from r 1 . For this reason, upper and lower envelopes are reported in the respective plots in Fig. 5 (a) and (b) to better appreciate the differences between the proposed and reference methods.
All three plots in Fig. 5 show that the proposed approach yields Rayleigh field realizations with improved spatial correlations than the reference method, in terms of maximum covered AV for a set number of PWs. In particular, the reference method exhibits the aforementioned noisy behavior sooner (∼ 2.6λ vs. ∼ 3.5λ), and in more pronounced fashion than the proposed method (see Fig. 5 (a) and (b)). Moreover, regarding ρ E z,z the proposed approach yields an excellent agreement with the theoretical behavior across the entire domain, while the reference method suffers from visible inaccuracies already between two and three wavelengths (see Fig. 5(c)). This good agreement may be explained by the greater uniformity of the z−projections of the spiral spherical sampling points, as shown in Fig. 1(b).

B. ACCURACY ASSESSMENT
Quantitative comparisons between the proposed and reference methods, in terms of spatial autocorrelations, were carried out based on the Pearson linear correlation coefficient VOLUME 10, 2022 between the computed spatial correlation curve and the theoretical one. This will be referred to as r i = r ρ i ,ρ i , where ρ is the theoretical correlation curve,ρ is the computed one, and the index i specifies the field or field component spatial correlation. This similarity measure is particularly suitable for the case under investigation, since both the computed and the theoretical curves are evaluated on the same set of distance points, hence a perfect match would yield the Pearson coefficient upper bound of r i = 1.
Furthermore, the Pearson coefficient can be computed for increasing distances, yielding clear indications of where the agreement with the theoretical behaviors begin to fall apart. Upon defining the distance-dependent correlation coefficient r i (x) as the Pearson coefficient between computed and theoretical curves up to a distance x, the maximum accuracy distanced is defined as the largest distance such that r i (d) ≥r i , wherer i is a predefined threshold value. Fig. 6 (a), (b) and (c) show the distance-dependent Pearson coefficient of the three tested spatial correlation curves ρ E , ρ E z,xy and ρ E z,z for N = 400 and R = 5000. As expected, distance-dependent correlations of both methods drop at increasing separation distances after an initial range where r i ∼ 1, with those of the proposed method decaying less rapidly and after a longer distance. In particular, Pearson correlation of the proposed ρ E z,z curve stays closer to 1 for about the whole spanned distance, while it falls abruptly around ∼ 3λ for the reference method, as shown in Fig. 5(c).
An approximated theoretical evaluation of the separation distance at which inaccuracy begins is given in [15], where a formula for the maximum AV radial dimension d rad is expressed in terms of the number of PWs impinging directions D. Recalling that N = 2D in [15], such a formula yields Based on numerical estimates of the linear correlation coefficient versus distance for ρ E , carried out for an increasing number of PWs (up to 14400, i.e. for L = 60), upon setting d rad =d, the threshold correlation value providing the best fit to (13) was found to ber E = 0.998. Hence, carrying out analogous linear correlation estimates for the proposed method, such a threshold value was adopted to arrive at a corresponding relationship, expressed as with the best fit observed for γ ≈ 0.80. Hence, as long as γ < 1, the proposed method yields larger AVs than the reference one for a given number N of PWs. As illustrated in Fig. 7, the number of required PWs to attain the predefined accuracy up to a desired AV size is consistently smaller for the proposed method. In fact, for increasing N , the ratio between the respective AVs is about γ 3 ≈ 0.5, implying roughly twice the volume when using the proposed method.

C. COMPUTATIONAL RUN TIMES
In order to assess the numerical burden, a comparative analysis of the computational times has been conducted. Table 1 reports the run times needed for each method to generate Rayleigh fields with a certain accuracy. As a first step, the reference method was used to generate Rayleigh field distributions with increasing number of PWs (N ), obtaining the maximum AV radial dimension values (atr E = 0.998) reported in the second row. Then, the proposed method was employed to generate fields with the same maximum AV and accuracy, which in this case requires less PWs (see Fig. 7). The computational set up made use of a Matlab R2019b (The MathWorks Inc.) code executed on a pair of NVIDIA® GeForce GTX 1080 Ti GPUs, for the analysis carried out up to d/λ = 25.
This benchmark test shows that the proposed method performs notably better for the same desired maximum AV (d rad = d rad ), especially for increasing AV size. In fact, taking the ratio between Eqs. (13) and (14) for the same, sufficiently large d rad , shows that the ratio between PWs (thus, between computational times, which depend linearly on N ) approaches 1/γ 2 ≈ 1.56 (i.e., one-third or more). This agrees well with the recorded computational times reported in the last row of Table 1, where the ratio between computational times obtained with the proposed and the reference method (speed-up factor) is reported.

IV. CONCLUDING REMARKS
A novel computational method to generate Rayleigh field distributions in three-dimensional space has been illustrated, based on the classical description of the field inside an ideal reverberation chamber in terms of PW superposition. An efficient and highly uniform sampling of the unit sphere has been implemented to define the PW propagation directions. Slanted linearly-polarized PWs with randomly defined components have been then associated to each direction and superposed within a predefined computational domain. Due to the PW analytical definition, the computational domain can be sampled with arbitrary resolution and employed as impressed sources in the solution of scattering problems.
The performance of the proposed approach has been quantified in terms of ensemble-average and spatial correlation field properties, comparing those of the synthesized fields with theoretical closed-form expressions, as well as those produced by a reference method described in [15]. An extensive benchmark analysis documented how the proposed approach provides an advantageous trade-off in terms of computational workload when compared with the reference method. In fact, it is shown that the proposed method requires significantly lower computational times for the same accuracy, especially for increasing AV size. Finally, starting from a threshold for the linear correlation coefficient between theoretical and computed Rayleigh field spatial correlation functions, a relationship between the number of PWs and the AV size was derived. Such a relationship indicated that, as the number of PWs grows larger, the proposed method yields almost twice the volume, compared with prior approaches.
Based on the accuracy and efficiency characteristics illustrated in the foregoing, the proposed method can be usefully employed in the design of state-of-the-art bioeffects studies, where the RC analysis volume must be maximized to enable the concurrent exposures of large animal cohorts.