Monte-Carlo Based Numerical Dosimetry in Reverberation Chamber Exposure Systems Employed for in-Vivo Rodent Bioassays

A Monte-Carlo based computational approach for the statistical characterization of the whole-body specific absorption rate (wbSAR) variability in large cohorts of rodents exposed to radio-frequency (RF) energy in reverberation chambers (RCs) is applied to adult male rat exposures illustrative of those in a US National Toxicology Program (NTP) cancer bioassay. A large number of 3D electromagnetic field realizations fulfilling Rayleigh fading properties were generated within an electrically-large volume representative of an ideal RC, yielding granular wbSAR distributions for an ensemble of 96 homogeneous rodent models with random mass distribution, postures, positions and orientations. Two case studies were addressed: a “momentary exposure” with each rat fixed in posture, position and orientation, and a “day-long exposure” in which the position, orientation and posture were varied randomly for each subsequent Rayleigh field realization. Over 500 and 2500 field realizations or “snapshots”, respectively, the rats’ instantaneous wbSARs, as well as their individual time-averaged wbSARs, were found to be well fit by lognormal distributions. The large variability in instantaneous wbSARs in the cohort was due in part to the inherent Rayleigh field variability in RCs (70-80%) and in part to weight, posture and position variations (20-30%), while the effect of cage location was found to be small over day-long exposures. Averaging the exposure over field realizations substantially reduces the range of wbSARs in the cohort. Hence, when RF-induced thermal effects are studied, the relevant exposure metric (wbSAR averaged over appropriate times) features a narrower range than instantaneous wbSAR, which is the relevant metric in studies dealing with non-thermal effects. Compared to previous studies, the present approach was found to be computationally more efficient enabling thus a Monte-Carlo analysis by varying concurrently the incident field and the animals posture, position, and orientation. In practice, it can inform the choice of wbSAR targets in rodent bioassay, allowing to identify possible dose-effect trends while avoiding undue thermal stress.


I. INTRODUCTION
Radiofrequency (RF) reverberation chambers (RCs) offer the opportunity to conduct large scale bioassays to investigate potential effects of RF exposures on large groups of The associate editor coordinating the review of this manuscript and approving it for publication was Zhongyi Guo .
animals [1]. RCs are electrically large, shielded enclosures in which a random non-stationary electromagnetic (EM) field is produced by varying the EM mode structure in the enclosure, typically by using a rotating ''mode stirrer'' [2], [3], [4].
RF energy distribution in RCs is characterized by ensemble statistics across field realizations [5], [6], which are relevant when RCs are used for bioeffects studies on animals. Notably, RCs designed to emulate Rayleigh fading (i.e., non-line-ofsight RF propagation) appear to offer an attractive exposure setting for large animal cohorts based on the expectation that theoretical spatially-uniform, time-averaged field magnitudes can be approximated in practice.
To date, the largest lifetime RF exposure bioassays, designed and conducted by the National Toxicology Program (NTP) [7], [8] of the US National Institute of Environmental Health Sciences (NIEHS), used large numbers of rodents (rats and mice) exposed from gestation for up to two years after birth. Based on pilot studies that determined lethal wholebody specific absorption rates (wbSAR) and the dependence of body temperature on wbSAR [9], distinct target exposure levels were defined for rats and mice. The RF dosimetry underpinning the NTP bioassays was studied both experimentally [10] and computationally [11].
Subsequent to the NTP cancer bioassay, some government agencies, such as the Federal Office for Radiation Protection (Bundesamt für Strahlenschutz, BfS) in Germany [12] and the Australian Radiation Protection and Nuclear Safety Agency (ARPANSA) [13], together with the International Commission on Non-Ionizing Radiation Protection (ICNIRP) [14], have recommended further research to address the reproducibility of the NTP study findings. Therefore, while the NTP announced follow-up studies that still employ RCs [15], a joint Korean-Japanese research team is conducting a replication of the NTP study, limited to male rats exposed at 900 MHz CDMA-modulated signals [16], [17].
In the design of in-vivo animal studies using RCs, numerical dosimetry provides key guidance in determining the SAR dependence on the incident field and to assess the wbSAR uniformity in exposed animals. Several groups have therefore performed numerical characterization on RCs [11], [18], [19], [20], [21], [22], [23], [24]. However, to limit the computational burden, most studies were restricted to a single [11] or a small number of concurrently exposed animals [18], [19], [20], [21], [22], [23]. Recently, the above mentioned Korean-Japanese team tackled for the first time the concurrent exposure of up to 75 rats, even though a true Rayleigh environment was loosely approximated using only 52 planewaves (PWs) [24]. In fact, the large RC working volume in [24] required about 1400 PWs for an accurate Rayleigh field representation (see Eq. (14) in [25]).
In this study, an efficient PW superposition method [25] has been employed to generate a near-ideal Rayleigh fading environment in an electrically-large RC analysis volume (AV) hosting a whole cohort of 96 animals. In this way, a Monte-Carlo approach that used many random inputs is proposed for the first time to provide a priori quantitative wbSAR statistics for the bioassay design (e.g., to define target exposure levels). Such an approach can also be employed a posteriori to recover key statistical parameters of the cohorts' actual exposure histories based on experimental observation (e.g., body mass logs). The proposed approach may also be broadened to study the variation in RF exposure to internal organs (which is not presently considered).

II. MODELS AND METHODS
The present approach presumes that an ideal Rayleigh propagation environment is synthesized within a sufficiently large RC volume for the expected rodent cohort arrangement. Hence, the Rayleigh field modeling can be performed regardless of the actual RC geometry. To this end, a method yielding accurate Rayleigh field distributions within a predefined AV at the study frequency (900 MHz) was developed in [25] and is presently used to examine the statistical distribution of wbSAR over a large cohort of homogeneous rat models.
The Rayleigh field distribution evolves over time in a periodic manner as the RC mode stirrer system completes a full revolution, typically over a few minutes, repeating indefinitely. In the results section, two case studies illustrate the implementation of the proposed approach in the statistical characterization of wbSAR corresponding to short-term, or ''momentary'', exposures of the rat cohort, as well as to longer term, or ''day-long'', exposure sessions that typically last several hours [7], [8]. Figure 1 illustrates the Monte-Carlo approach, comprising three main steps. In the first one, the whole cohort of simplified, homogeneous rat models is determined via a random assignment of their weight, position, orientation and posture. In the second step, a random Rayleigh field realization is generated using the approach in [25] and enforced as a Huygens box equivalent source in the Finite-Difference Time-Domain (FDTD) simulation environment. In the last step, the FDTD simulation yields the wbSAR for each rat model, which is collected for post-processing.

A. IDEAL RAYLEIGH ENVIRONMENT SOURCE MODEL
The electromagnetic field of a ''well stirred'' RC (sufficient to randomize the fields within the AV) is described through the superposition of an infinite number of PWs [2]. Rayleigh fading is characterized by broad swings of the electric and magnetic power densities, proportional to the local |E| 2 and |H| 2 respectively, with identical statistics characterizing both. Each in-phase and in-quadrature field component sampled at an arbitrary observation point within the AV, across random field realizations, is then an independent random variable. The distribution of these samples follows a normal ( ) probability distribution function (PDF) with zero mean (µ = 0) and population variance σ 2 Denoting with E 2 0 the squared electric field magnitude |E| 2 ensemble mean (i.e., across all random field realizations) at that AV point, then σ 2 = E 2 0 /6 (σ 2 = E 2 0 /6η 2 0 for the magnetic field, η 0 being the free-space wave impedance) [5]. Hence, the corresponding PDFs for |E| 2 and |H| 2 are chisquared with 6 degrees of freedom, denoted χ 2 6 : representing a special case of the distribution where α = n/2 and β = 2σ 2 are the shape and range parameters, and n = 6 are the degrees of freedom. The mean, standard deviation (SD), and mean-normalized (i.e., relative) standard deviation (R-SD) are αβ, β √ α and 1/ √ α respectively. As α = 3, the R-SD of the distribution in |E| 2 will be ≈ 58%, hence, Rayleigh fading will lead to temporal variability of the rats exposure. The resulting wbSAR distribution in the exposed rat population may differ from that for |E| 2 since the SAR distribution in a rat's body stems from the spatially distributed coupling with the incident Rayleigh field, thus mitigating deep point-wise fades.

B. IMPLEMENTED RAYLEIGH ENVIRONMENT
In the present approach, the field in ideal Rayleigh fading environments is approximated using a semi-analytical method [25] borrowing from West et al. [26], where a novel algorithm was developed to superpose an arbitrarily large number (N ) of equal-amplitude PWs with randomly slanted linear polarizations and propagation directions, obtained via the uniform spherical spiral sampling scheme illustrated in Fig. 2. Compared with prior literature, this approach leads to better agreement with the ideal Rayleigh fields' ensemble statistical properties [5], [6]. In the case studies illustrated in Section III, the method in [25] was expanded by adding a second randomly slanted PW for each impinging direction, resulting in N = 5000 elliptically-polarized PWs.
The present approach differs from that adopted in [11], where a method mainly based on 12 linearly-polarized PWs propagating along the cardinal directions was employed to illuminate a single rodent model rather than an entire cohort. While a single anatomical rodent model provided the opportunity to characterize both wbSAR and organ SAR (oSAR) at various stages of growth, the present approach inherently addressed whole cohort exposures in an ideal RC field, providing individual as well as collective wbSAR statistics based on posable, homogeneous rat models. In this case, the synthesis of an almost ideal RC environment over an electrically-large AV capable to host a whole cohort is of paramount importance. This is well accomplished in our study, as documented in Appendix A, while using 12 PWs would not realize a true Rayleigh environment over such a large AV (see Eq. (14) in [25]). This limitation also applies to the 52 PW approach in [24], since ∼1400 PWs would be needed [25] for their AV (1.5 m × 1.5 m × 1.3 m) used to expose the whole cohort of adult anatomical rat models. A collection of M volumetric Rayleigh field realizations (M = 500 and M = 2500 in the ''momentary'' and ''daylong'' case studies, respectively) with 1/20 th of a wavelength step (≈1.67 cm at 900 MHz) were synthesized from N elliptically polarized PWs, over a uniform grid formed by 109 × 86 × 107 AV points. Equivalent Huygens sources over the AV bounding box (approximate dimensions 1.80 m × 1.42 m × 1.77 m) were then derived upon sampling the AV boundary fields (see Fig. 3). By enforcing the incident Rayleigh fields, through the Huygens box, the mutual coupling between rats and their RF losses are accounted in the FDTD simulations.
Appendix A illustrates the validation of the synthesized Rayleigh field distributions, together with examples of volumetric field distributions for an empty RC, as well as for an RC loaded with the homogeneous rat models described in the following. The validation of the Huygens source implementation in the FDTD solver (Sim4Life v6.0, Zurich MedTech, CH) is illustrated in Appendix B.

C. RAT MODELS AND COHORT EXPOSURES
The rat models and their disposition within the AV were shaped after the salient parameters of the NTP bioassay [7].

1) COHORT EXPOSURE
Two cage racks holding up to 56 rats each were placed side by side in the NTP study [7]. Each rack formed an incomplete 6 × 2 × 5 lattice where four rats at most were placed at the bottom layer (maximum of 112 rats per RC). The cohort modelled in this study instead consisted of R = 96 rat models within a uniform lattice resembling a single large cage rack holding a regular 6 × 4 × 4 cage array along x, y, and z, as shown in Fig. 3.
The plastic cages are essentially transparent at 900 MHz, while the metal tubing was not included in the simulations to preserve as close an ideal Rayleigh incident field as possible. The rat models were positioned on the floor of individual cage volumes having the dimensions of the cages used in the NTP rats bioassay (referring to Fig. 4: 23.5 cm × 26 cm × 21 cm in x, y, z), while approximate horizontal (5 cm in x and 18 cm in y) and vertical (10 cm in z) distances between cages were inferred from RC photographs in the NTP reports [7]. TABLE 1. WbSAR statistics (mean, SD, R-SD and ratio between wbSAR and incident mean square E-field) for the anatomical rat model (Sim4Life, Zurich MedTech, CH) and the simplified rat models computed across M=1000 Rayleigh field realizations. Incident mean square E-field (in absence of the rat models): 5000 (V/m) 2 .

2) SIMPLIFIED RAT MODELS
To collect representative wbSAR statistics for the rodent cohort, we opted for simplified, posed lossy homogeneous rat models and validated them against an anatomical rat model. As seen in Fig. 5, rats were modelled without tail starting from an assembly of four parametric ellipsoidal half-wedges, with properties taken from [11]: density ρ = 1000 kg/m 3 , conductivity σ = 0.95 S/m and relative permittivity ε r = 40.
A posable anatomically realistic adult male rat model (www.itis.ethz.ch/animals), scaled up uniformly along x, y, and z from 567 g to 647 g, provided a reference to define four postures: resting straight, elongated straight, asleep curved clockwise and asleep curved counter-clockwise. The rats drinking posture is most similar to the elongated model posture [10], [11]. The straight models were approximated first by combinations of ellipsoidal half-wedges and subsequently deformed to obtain the ''asleep'' curved models, mirroring each other (see Fig. 5).
WbSAR values were determined from FDTD simulations, employing the 900 MHz equivalent Huygens sources derived from the set of Rayleigh field realizations in order to obtain voxel SAR, post-processed to yield wbSARs. Field sensors placed in each rat model were monitored to ensure simulation convergence (−50 dB). The variable grid feature available in Sim4Life was employed to discretize all rat models in cubic 3 × 3 × 3 mm 3 voxels, allowing the grid to relax up to 1.5 cm in the surrounding air volume. The FDTD grid discretization caused slight deviations from the rat models nominal weights. Consequently, the voxel-model masses were employed to calculate wbSAR in the FDTD solver.
Based on FDTD simulations involving M = 1000 Rayleigh field realizations, the anatomical rat model (both with and without the tail) provided validation benchmarks for the homogeneous rat models in terms of dosimetric exposure. Table 1 provides a comparative summary of the wbSAR means and SDs, indicating close agreement between the realistic and simplified rat models. The ratio of wbSAR and mean square of the incident E-field is a key parameter to control wbSAR during exposure sessions. The corresponding values reported in Table 1 are in line with [11] and [24]. VOLUME 11, 2023 FIGURE 5. The simplified geometrical rat models (resting straight, elongated straight, asleep curved clockwise and asleep curved counter-clockwise) compared with the respective posed anatomical rat models. The simplified rat models were color coded as shown.

3) MASS DISTRIBUTION OF THE RAT MODELS
The four homogeneous models in Fig. 5 weight about 647 g (see Table 1) before any scaling needed to attain the assigned weights, based on the observed weights means [7, Tables 17 and 44] and R-SD [27, p. 36] of the adult male rats (months 13-24), about 650 g and 11.6 %, respectively. Thus, each one of the 96 homogeneous rat models in the exposed cohort was assigned a mass value from a raised-cosine mass distribution and then scaled uniformly along x, y, and z accordingly. The finite support of the raised-cosine PDF avoids outliers. Figure 6 reports the histogram of the randomly generated mass set (mean 653.5 g, SD 74.8 g) used in the case studies illustrated in the following.

D. NUMERICAL DOSIMETRY AND STATISTICAL ANALYSIS 1) CASE STUDIES
Individual and collective wbSAR statistics were evaluated for two exemplary cases, formulated to model ''momentary'' and ''day-long'' exposure scenarios. These case studies employed M = 500 and M = 2500 Rayleigh field realizations, respectively. For each field realization, an FDTD simulation provided a wbSAR value for each rat model. In both scenarios, each cage would host the same rat across field realizations; hence a weight would be permanently associated to each cage and the rat model hosted therein. Key elements of these case studies are summarized in Table 2.
The ''momentary'' exposure scenario involved field realizations emulating Rayleigh fading over short-enough time spans for the exposed rats to be assumed motionless. As discussed in Section III, this scenario allowed singling out the wbSAR variability associated with the Rayleigh fading. Each rat model was selected from one of the four postured versions, deemed to be equiprobable, and then scaled uniformly to attain the weight associated to the assigned cage. Upon assigning a random orientation angle about a vertical axis from a uniform distribution in [0,2π), the rat model was positioned on the cage floor according to uniform rectangular displacement distributions in x and y, making sure it would fit within the assigned cage volume. In this scenario, the overall geometry was kept fixed across successive FDTD simulations, differing only for the Huygens source employed. Well converging wbSAR statistics were achieved after 500 field realizations (see Appendix C).
In the ''day-long'' exposure scenario, the rat model posture, position and orientation were randomly assigned for each Rayleigh field realizations. This scenario simulated the cohort exposure over extended time spans (e.g. halfday exposure sessions), collecting wbSAR ''snapshots'' at large enough time intervals that would reasonably allow all the animals to have moved around within their cages, going to sleep or staying active. Since rodents are long sleepers, we included two sleeping postures and two active ones. In this case, 2500 realizations yielded converging wbSAR statistics.

2) STATISTICAL ANALYSIS
In the adopted PW-based Rayleigh fading modelling, each snapshot may be associated to a fixed position of hypo-thetical RC mode stirrers. As the stirrers rotate, a distinct incident field distribution is realized at each time instant.
Accordingly, the wbSARs calculated in a single FDTD simulation represent a ''snapshot'' in time of the evolving instantaneous wbSARs in the rodents. Instead, time-averaged wbSARs can be estimated by averaging the wbSAR in each rat over a sufficiently large set of snapshots.
Instantaneous and time-averaged wbSAR may play different roles in the context of an animal bioassay. When seeking potential thermal effects, time-averaged wbSAR is the exposure metric that better correlates with the thermal load imposed on the animals. Conversely, when the research focus is on non-thermal bioelectromagnetic interactions, the instantaneous wbSAR better represents the range of internal fields within the animals' bodies.
The present approach, built on a Monte-Carlo stochastic simulation model yielding the individual rats' instantaneous wbSARs over arbitrarily large collections of incident field snapshots, inherently suits both endeavours by providing cumulative distribution functions (CDFs) and PDFs for instantaneous and time-averaged wbSARs (see Section III).

III. RESULTS
In this Section, the main contributing factors to the wbSAR variability under ideal Rayleigh fading are inferred on the basis of the case study results, which are described first. The wbSAR data is based on the incident Rayleigh fields collection described in Appendix A, yielding spatiallyuniform mean square E-field (E 2 rms ) levels of ∼5000 (V/m) 2 .
The distribution of instantaneous wbSAR values (single snapshots) for the whole cohort (500 wbSAR values for each of the 96 rat models) is depicted in Fig. 7. As shown, the histogram is fitted very well by a lognormal PDF having same mean (0.2265 W/kg) and R-SD (37.4 %), given by whereμ andσ are respectively the arithmetic mean and the SD of the logarithms of adimensional wbSAR samples.
The fitted parameters and their confidence intervals (CIs) are reported in Table 3. The lognormal wbSAR distribution is less skewed than the |E| 2 distribution (dashed line in Fig. 7, depicting a χ 2 6 PDF having same mean and SD of the lognormal PDF), confirming the expectation that the spatial averaging involved in the wbSAR determination would smooth out the Rayleigh field deep fades and overshoots.
The ratios between mean wbSAR and incident mean square E-field in Table 3 are lower for both case studies than for the isolated rat models (Table 1), under statistically VOLUME 11, 2023 FIGURE 9. WbSAR for a single rat over M=500 field realizations in the ''momentary'' scenario. (a) WbSAR in that rat in the 500 realizations; (b) distribution of wbSAR in that rat over the M=500 field realizations. Vertical lines represent the mean (0.2265 W/kg) over the 500 snapshots, taken to represent the time-averaged wbSAR for that rat. Incident mean square E-field: 5000 (V/m) 2 .
alike incident fields (i.e., near-ideal Rayleigh environment). These ratios are in line with those in [24,Table 4], where the whole cohort dosimetry was carried out taking into account the actual RC environment (metallic water supply system).
In Fig. 8, the lognormal distribution parameters (μ,σ ) were fitted to the whole wbSAR data set (96 rats * 500 snapshots for each rat) using the Matlab distributionFitter application (Matlab 9.2, The MathWorks, Inc., Natick, MA, US). Clearly, the distribution of the computed wbSAR values almost precisely follows a lognormal distribution, with negligible errors on mean and SD. The 5%-95% probability range in Fig. 8 corresponds to about 0.26 W/kg wbSAR span, roughly as wide as the distribution mean (0.2265 W/kg).
The 500 wbSAR values for a single representative rat model is displayed in Fig. 9(a). The mean wbSAR for that rat is 0.2338 W/kg and R-SD is 37.4%, similar to that for the cohort at large. The corresponding wbSARs distribution in  While the 96 rats in the momentary exposure scenario were fixed in position, the total RF power absorbed in the cohort varied with different field realizations due to variations in incident field distributions. The total dissipated RF power in the 96 rats was 14.1 W with ≈4.4% R-SD across 500 Rayleigh field realizations, which is in line with ≈3.1% R-SD observed over the different snapshots for the spatiallyaveraged incident |E| 2 in the AV.
Exposure uniformity across cages is also of interest. Thus, cage locations were grouped based on the number of cage side views obstructed by neighboring cages resulting in 8 cages at the ''rack'' corners (three unobstructed side views), 32 edge cages (two), 40 face cages (one), and 16 inner cages (none). As shown in Table 4, about 20% spread in the mean wbSAR was observed among cage locations, with the animals with more unobstructed views (fewer adjacent cages) having higher mean wbSARs than those that are completely surrounded by other cages.

2) ''DAY-LONG'' EXPOSURE SCENARIO
To allow meaningful comparisons with the previous scenario, the same set of mass values (see Fig. 6) was associated to each cage. However, as the rats were moved around at each next field snapshot, the FDTD grid was consequently updated, introducing a bit of noise in the rat model mass (R-SD ≈0.5%), but with negligible impact on wbSAR statistics. Additionally, the total RF power dissipated in the cohort had very similar statistics to the previous case study (14.1 W mean with ≈4.8% R-SD).
An individual CDF was determined for each of the 96 rat models. In Fig. 10, this CDF family is plotted together with the histograms relative to the rat models associated to the CDFs bounding the curve ensemble. The CDF for the whole data set is reported as well, and a normalized abscissa axis is defined about the median of this CDF to mark the inter-group bounds (for a two-fold target wbSAR progression [7]) at √ 2 (upper) and 1/ √ 2 (lower). The inset histograms provide granular exposure ranges and distributions for individual rats.   Table 5). Dashed lines delimit the 5-95% percentile range (≈0.065 W/kg) about the mean (0.2261 W/kg). Incident mean square E-field: 5000 (V/m) 2 .
The rat corresponding to the flattest CDF would be exposed within a higher wbSAR group range for more than a third of the time, while the rat associated with the steepest CDF would be exposed within a lower wbSAR group range for more than a fourth of the time. CDFs grouped about the cohort average imply actual exposures falling outside the intended group range for about one-third of the time.
As found in the ''momentary'' exposure'' case, the instantaneous wbSARs distribution (2500 values for each of 96 rat models) is fitted very well by a lognormal PDF having the same mean (0.2261 W/kg) and R-SD (37.8 %).
A comparison between the cohort wbSAR statistics for the two scenarios, and the relative lognormal fitting parameters, is summarized in Table 3. The substantial overlap of the respective CIs makes the PDFs statistically indistinguishable. In Fig. 11, the distribution of the time-averaged wbSARs (mean across M =2500 field realizations) for each rat model is also shown to be well approximated by a lognormal PDF, constituting another useful finding as it readily enables the modeling of the cohort time-averaged wbSARs. As expected, the 5-95% percentile range corresponds to about 0.065 W/kg span, much narrower than for the instantaneous wbSAR.
In Fig. 12, histograms of the wbSAR means and medians for the 96 rat models show a narrower range of variation compared with the instantaneous wbSAR (see Fig. 9(b)), with the individual time-averaged wbSARs spanning between 0.188 W/kg and 0.28 W/kg (83% to 124% of the collective mean wbSAR).
Finally, correlation between time-averaged wbSAR and rat mass may be of interest. A negative Pearson correlation coefficient (−0.72) was observed between the rats mass values and the relative wbSARs. However, positive correlation (0.67) with the absorbed RF power indicated that heavier rats receive larger amounts of RF-induced heat than lighter ones, albeit their larger mass lowered the wbSAR.
Due to larger volume to surface ratio, heavier rats may release excess heat less efficiently, unless heat release through the tails can play an equalizing role.
For what concerns the individual rats wbSAR histograms, distributions similar to the one reported in Fig. 9(b) for the ''momentary'' case study were observed in this case as well. However, the observed spread in the mean wbSAR among cage locations was much lower (≈5%) compared with ''momentary'' exposures, as shown in Table 4. This suggests that prolonged exposures, providing the rats an opportunity to move around, mitigate the influence of cage location (at least under ideal Rayleigh fading conditions).

B. MAIN CONTRIBUTIONS TO VARIABILITY IN WBSAR
Given the marginal role of cage location on wbSAR variability, particularly in ''day-long'' scenarios (see above), two other factors assume a prominent role: variability in the fields incident on each rat (due to Rayleigh fading), and variability VOLUME 11, 2023 TABLE 5. Fitted lognormal distribution parameters of 96 (per rat) mean wbSAR data sets for the ''momentary'' and ''day-long'' case studies. The fitted lognormal PDF in the ''day long'' scenario is displayed in Fig. 11. Incident mean square E-field: 5000 (V/m) 2 .
in their mass, orientation, and posture. While orientation and posture are difficult to track over years long studies, the rat weights histories are usually documented (e.g., [7], [8]), enabling correlation estimates with wbSAR.
The following discussion will distinguish between the distribution of instantaneous wbSARs in the 96 rats and the individual rat time-averaged wbSARs (i.e. averaged over field realizations). Averaging wbSAR distributions for each rat separately can help to discriminate between these uncertainty sources by averaging out effects of field variation. This is particularly true for the ''momentary scenario'' where the rats remained motionless for all 500 field realizations, but it seems to be the case also for the ''day-long'' scenario given the similarity between the two distributions in wbSAR (Table 3). Indeed, comparing Tables 3 and 5, the distributions in the ''per rat'' time-averaged wbSARs feature SDs about one-quarter those of the lognormal distributions that fit the cohort results. This suggests that variability in the incident fields on each animal accounts for ≈75% of the overall instantaneous wbSAR variability. The wbSAR distribution similarities between the ''day-long'' and the ''momentary'' scenarios suggested the same conclusion. From Table 5, the fitting lognormal PDFs for the ''momentary'' and the ''daylong'' exposure scenarios are once again statistically alike given the substantial overlap between the respective CIs.

1) EFFECT OF THE RAT MODEL WEIGHT
Negative correlation was observed between individual timeaveraged wbSARs and the rats body mass, as clearly shown in Fig. 13, for the ''day-long'' exposure scenario. The correlation with the 96 rats medians is approximately the same (−0.70) as for the means (−0.72). Hence, the spread of body mass within the cohort will produce significant but opposite spread in mean and median wbSAR.

2) EFFECT OF TIME-AVERAGING
The effects of time-averaging can be retrieved in the present model by averaging wbSARs over a growing number of field  snapshots. This is conveniently done in Matlab using the function movmean, applied separately to each of the 96 sets of 2500 wbSAR results, leading to mean and SD histograms for each rat model. The 96 SD distributions (one for each rat) are reported in Fig. 14 for an increasing number of averaging window lengths. As expected, the SD progressively decreases with increasing window length since time averaging reduces the exposure variability.
These simulations show that the instantaneous exposure (i.e., the rats' wbSARs calculated for each field snapshot) is highly variable, mainly due to Rayleigh fading, which is an inherent feature of any RC. However, the long-term exposure (presently modeled in terms of averages over large numbers of Rayleigh field realizations) reflects a combination of variability in time-averaged fields throughout the RC (which is small in the present model, as illustrated in Appendix A) as well as variability in characteristics of the animals (posture, position,etc.). This variability will strongly affect the tails of the distribution in the exposure, which would have to be evaluated for each experiment. Statistical characteristics of the incident Rayleigh field (peak square E-field) derived through repeated samples from the gamma distribution with parameters given in Fig. A4, and of the ''day long'' wbSAR distributions derived applying the movmean fuction to the instantaneous wbSARs. Table 6 summarizes the statistical characteristics of the |E| 2 and wbSAR distributions.

C. SUMMARY OF FIELD AND WBSAR DISTRUBUTIONS
The distribution of |E| 2 represents the distribution of averages of M values of repeated samples of 1,003,018 values (the AV point grid) from the gamma distribution having the parameters given in Fig. A.4.
The distributions of wbSAR were derived applying the movmean function to all 240,000 wbSAR evaluations in the ''day long'' case study (96 rats over 2500 field snapshots), considered individually (M=1) or for averages of 500 or 2500 wbSAR evaluations for each of 96 rats.
The wbSAR statistical characteristics in Table 6 can be taken into account in the design of rodent bioassay, as discussed below.

IV. DISCUSSION AND FUTURE RESEARCH
The design of well-performing RCs is not a trivial task [3], [4], [28], and the field homogeneity and variability in any RC to be used for bioeffects studies needs to be thoroughly characterized [10], [17].
The performance of state-of-the-art RCs in terms of the R-SD of the fields averaged over time (i.e., stirrer states) is specified by IEC 61000-4-21:2011 standard [29]. In order for the field within the working volume (WV; defined in [29]) to be considered uniform, the standardized limit (R-SD not exceeding 3 dB above 400 MHz) still allows wide swings of the mean squared field magnitude within the WV. For example, the NTP rats bioassay report mentions ±2.5 dB incident field swings [7].
For the simulated RC field distributions employed in this study, the R-SD of the ''time-averaged'' incident field magnitude (i.e., the M=500 or M=2500 entries in Table 6) is negligible, with relative standard deviations of about 1% for M=2500. This field uniformity is far below the specification limit of 3 dB for time-averaged field variability allowed by the IEC 61000-4-21:2011 standard at 900 MHz. Even the field distribution in a single snapshot complies with IEC 61000-4-21. Hence, the implemented Rayleigh environment is an excellent approximation of an ideal RC. This high level of field uniformity, averaged over time, would be difficult or impossible to accomplish in a real RC.
Indeed, even in this case, the variability in |E| 2 implies that at any instant, a sizable fraction of the animals will be exposed at considerably higher or lower levels than the mean for the cohort. Unless wbSAR levels are already so high to cause thermal stress to the animals, short-term fluctuations in exposure may be inconsequential for studying thermal effects characterized by slow response times. However, they might be significant when non-thermal effects are of interest.
The proposed RC field modelling approach intentionally excludes perturbing effects due to RC walls and mode stirrers. Moreover, the present implementation does not include cages and cage racks, although they are essentially transparent to RF fields and may be included in further refinements. The metal water supply pipes used in the NTP racks and the cage-attached drinking fixtures [7] could also be modelled reliably unless the watering system is grounded to the RC walls. Therefore, these simplifications indicate that the results presented in the foregoing constitute an upper bound of the achievable exposure system performances with actual RCs. The proposed approach also assumes that none of the cages are in dominant (direct) propagation paths from the sources, another best case assumption since departing from ideal Rayleigh fields, e.g., due to the onset of dominant propagation paths, will produce systematically higher exposures in some caged animals.

A. RC EXPOSURE IN THERMAL STUDIES
Unforeseen thermal stress would likely be unacceptable in most bioassays, introducing an inextricable confounding element in the interpretation of biological outcomes. The present approach may be useful to indicate when thermal stress cannot be ruled out so that tracking stress biomarkers in the animal populations may be planned. In such cases, given the aforementioned wbSAR dependency on body weight and the fact that female and male rat weight distributions may differ (e.g., [7], [27]), tracking sex specific biomarkers over their life spans may be appropriate as well.

B. RC EXPOSURE IN NON-THERMAL STUDIES
When employing an RC exposure systems in studies focusing on low-level non-thermal bioeffects (e.g., in-vivo cellular level field interactions), the instantaneous wbSAR and oSAR distributions become then of the essence. In these cases, the present approach can be employed at the study design stage to quantify the instantaneous SAR range of variability (e.g., deriving collective and individual instantaneous wbSAR histograms like those in Figs. 7, 10) and set the separation between target wbSAR levels for different exposure groups accordingly.
Another potentially important aspect to address in lowlevel studies concerns the modulated RF signal waveforms VOLUME 11, 2023 employed. For example, 3G/4G signal waveform envelopes are characterized by a large peak-to-average power ratio (PAPR), which will stochastically combine with the Rayleigh field own fluctuations and increase the instantaneous wbSAR R-SD. In these cases, further separation between wbSAR target levels may increase the chances of quantifying strengtheffect relationships.

C. FUTURE RESEARCH
Frequently, rodent bioassays are carried out at wbSAR levels that approach or exceed 4 W/kg (e.g., [24], [30]), which is an acknowledged threshold for the onset of thermal effects. In fact, current safety standards [31], [32] define wbSAR limits upon applying reduction factors from 4 W/kg. In these cases, where the thermal load in the animals should be controlled, the present approach may be improved by using anatomically realistic rodent models and then including appropriate thermo-physiological modeling to yield an expected range of temperature in the animals. Such an analysis would provide statistical information on oSAR and organ specific temperature distributions.
Another refinement of the present approach would entail the intentional synthesis of less-than-ideal Rayleigh fields, to provide wbSAR variability expectations for actual, stateof-the-art RC design. Conceptually, this can be accomplished by degrading the accuracy of the Rayleigh field synthesis technique in [25], driving the time-averaged energy densities to gradually become less uniform across the AV. This would allow analyzing wbSAR/oSAR and body/organ temperature variability when the RC field uniformity is non-ideal but still within the standard requirements [29]. Employing less than ideal Rayleigh fields would also provide insights on the impact of empty cages due to mortality, especially useful to model exposures during the latter part of a bioassay.

V. CONCLUDING REMARKS
Reverberation chambers offer efficient means to expose groups of small animals to RF energy. However, even in RCs that meet standardized specifications the exposure is highly variable across space and time and has to be statistically quantified. Accomplishing this through numerical dosimetry involved a large computational burden often requiring approximations and simplifications [11], [18], [19], [20], [21], [22], [23], [24]. Hence, deterministic approaches based either on a small set of configurations (whole-cohort) or a reduced computational domain (single animal) have been typically adopted so far.
In this work, a novel stochastic and computationally efficient approach to quantify the wbSAR variability across rodents in large cohorts concurrently exposed to RF fields in near-ideal RCs has been proposed for the first time using a large and random series of inputs. The main contributions include aspects of the computational methodology, the statistical analysis design, and a number of findings that can be readily used by researchers in modelling rodent exposures in RCs.
Near-ideal Rayleigh field realizations were efficiently implemented through equivalent sources in the FDTD solver to carry out an extensive Monte-Carlo characterization of individual and collective, instantaneous and time-averaged wbSARs. This was preceded by the definition and validation of lossy, posed homogeneous rodent models suited to yield wbSAR statistics for arbitrarily large rodent cohorts within an efficient computational framework.
The case study results allowed quantifying the effect of the cage location in two case studies modelled after the parameters of an important study [7]. Furthermore, it was found that lognormal PDFs fit both instantaneous and time-averaged wbSAR histograms of the 96-rat cohort, providing a reliable tool for researchers to model expected wbSAR ranges. The results illustrate that instantaneous wbSAR variability was due to inherent RC field variability (70-80%) and in part to weight, posture and position (20-30%). Instead, the cage location was found to have a small effect over day-long exposures. Finally, the results were also employed to determine the correlation between rodent weight and time-averaged wbSAR, as well as documenting wbSAR time-averaging effect under ideal Rayleigh fading. However, if the actual experimental setup (RC field inhomogeneities, cages, racks, watering system, etc.) causes significant changes of the statistics across position (i.e., cage location), then the weight, posture and position dependencies may not necessarily represent the main variables, besides the RC field itself, influencing the wbSAR variability.
The proposed approach is fully parametric and accounts for key stochastic variables in order to characterize statistically key dosimetric observables, such as the wbSAR and even organ-averaged SAR. It was shown to yield very granular individual and collective wbSAR statistics (e.g., histograms and CDFs). When needed, it can be supplemented by experimental verification and subsequent modelling refinements. In this way, it can be purposefully employed to gather insightful guidance in setting target exposure levels that may enable discerning dose-effect relationships while avoiding confounders such as the onset of thermal stress in the animals. These are critical aspects in the design of RF exposure bioassays, which often yield statistically weak findings or no findings at all.

APPENDIX A. RAYLEIGH FIELD SYNTHESIS VALIDATION
Applying the Rayleigh field synthesis approach described in [25], M = 2500 volumetric realizations were computed across a uniform grid of 109 × 86 × 107 AV points, with ≈1.67 cm spacing between grid points. This resulted in a set of 1,003,018 |E| 2 and |H| 2 values for each of the M Rayleigh field realizations used in the wbSAR simulations.
An exemplary distribution of the electric and magnetic field magnitudes across an AV plane is provided in Fig. A.1. As expected, the spatial field behavior features the typical peaks and valleys occurring in classic non-line-of-sight propagation scenarios causing the Rayleigh fading. The positions of these peaks and valleys vary randomly across different field realizations.
Corresponding distributions of the electric and magnetic field magnitudes in the presence of elongated rat models, across the same AV cut-plane are presented in Fig. A.2, illustrating the strong coupling between the incident Rayleigh field and the lossy rat models. Pronounced coupling with the magnetic field indicates the onset of eddy currents, thus confirming the SAR enhancement findings from the NTP bioassay dosimetry studies at 900 MHz [10], [11], [27]. The SAR distribution, shown in Fig. A.3, indeed displays visible spatial correlation with the magnetic field.
The histogram of the ensemble (i.e., time-averaged) mean |E| 2 is reported in Fig. A.5 for M =2500 field realizations, showing a narrow distribution with R-SD ≈1.2% about the anticipated value (10 4 V 2 /m 2 ). The R-SD decreases with Rayleigh field realizations (M ), as seen upon applying the movmean Matlab function to the |E| 2 dataset, yielding the |E| 2 mean, SD, and R-SD shown in Table A.1.

APPENDIX B. HUYGENS BOX SOURCE VALIDATION
The validation of the FDTD Huygens box source implementation is illustrated with reference to Fig. B.1, where the percentage difference of the squared E-field magnitude between the Matlab-generated and the FDTD-generated Rayleigh fields indicates an average difference of 0.37% with 1.23% SD and 15.7% maximum error (occurring about a field null) across a cut-plane in the AV. The validation of the Huygens source implementation was expanded through the inclusion of two lossy spheres (15 cm radii, σ = 0.95 S/m and ε r = 40) in the AV, as illustrated in Fig. B.2. In this case, we compared the simulation where 500 random PWs were generated concurrently FIGURE B.1. Percentage difference between |E| 2 values derived from the Rayleigh field generation method proposed in [19] and from the FDTD implementation of the corresponding equivalent Huygens source within an empty AV. (a) False-color representation of the percentage difference across an AV cut-plane (axes ticks indicate the FDTD grid points); (b) corresponding percentage error distribution compared with a normal distribution featuring the same mean and variance.
at 900 MHz within the AV, and that where the corresponding equivalent Huygens source was impressed. The error analysis was conducted separately for the whole domain (mean error 0.34%, SD 0.37%, maximum error 10%) and for the lossy spheres (mean error 0.28%, SD 0.27%, maximum error 1.8%), indicating that introducing losses in the AV improves the Huygens source accuracy.

APPENDIX C. WBSAR CONVERGENCE STUDY
The minimum number of realizations required for the two exposure scenarios described in Table 2 was determined tracking the R-SD of the wbSAR distributions across the anatomically detailed and simplified rat models described in  Table 1. The analysis involved up to M = 1000 Rayleigh field realizations. As illustrated in Fig. C.1, the wbSAR running R-SD substantially converged after 500 field realizations for all animal models, yielding a criterion to set the minimum M .

ACKNOWLEDGMENT
This work was supported by the Mobile & Wireless Forum (MWF), which had no influence over the execution and findings of this study. VOLUME 11, 2023 FIGURE C.1. Ratio between the running wbSAR SD and the running wbSAR mean for the homogeneous rat and the anatomical models. R-SD values calculated over M=1000 Rayleigh field realizations are reported in Table 1.