On the Dispersive Nature of the Power Dissipated into a Lossy Half-Space Close to a Radiating Source Andrea Cozza, Member, IEEE, and Benoît Derat, Member, IEEE Abstract-In this paper, we study how the power dissipated into a lossy medium excited by a nearby antenna is affected by drifts in the electrical parameters of the lossy medium. The statistical distribution of the sensitivity of the dissipated power is determined by means of a spectral analysis of the transmission of electromagnetic energy from air into the lossy half-space. A clear link is drawn between the reactive content of the field excited by the source and the dispersiveness of the sensitivity. The case of a stratified structure is also addressed, by defining a modification factor representing the alteration of the transmissivity and of its sensitivity when a buffer layer is introduced. All of the results provided point out that, in general, the sensitivity of the total amount of power dissipated into the half-space cannot be predicted independently from a precise knowledge of the source characteristics, unless under a paraxial propagation approximation or in a far-field configuration.

Index Terms-Antennas, Lossy Media, Near Field, Plane-Wave Spectrum, Sensitivity Analysis, Specific Absorption Rate, Stratified Media.

I. INTRODUCTION

In many practical configurations antennas stand near to a lossy half-space; an example of historical importance is that of radiating sources placed over a lossy soil [START_REF] Baños | Dipole Radiation in the Presence of a Conducting Halfspace[END_REF]. Many other applications involve a similar scenario: radio-frequency hyperthermia [START_REF] Carr | Antenna: the critical element in successful medical technology[END_REF] and Specific Absorption-Rate (SAR) assessment [START_REF] Kuster | Energy absorption mechanism by biological bodies in the near field of dipole antennas above 300 MHz[END_REF] are but two examples. All of these configurations share a common concern, that of being able to estimate the amount of power that will be dissipated into the lossy medium. In the context of the present discussion, we will refer to the concept of Total Dissipate Power (TDP). Depending on the actual application, the TDP needs to be maximized or minimized, but all the same its value closely depends on the electrical properties of the concerned medium, namely the electric conductivity σ and the relative permittivity ǫ r . The basic need of a proper estimation of these parameters thus appears as a potential difficulty in the computation of the TDP; hence, it is fundamental to investigate how the TDP can be affected by a variation of the electrical parameters of the medium, be that due to a drift in time or rather to an uncertainty in their knowledge.

The primary aim of this paper is twofold: 1) to provide a model capable of predicting the sensitivity of the TDP A. Cozza is with SUPELEC, Département de Recherche en Électromagnétism, Gif-sur-Yvette, France, e-mail: andrea.cozza@supelec.fr.

B. Derat is with Field-Imaging, Meudon, France, e-mail: benoit.derat@field-imaging.com.

Manuscript received ..., 2008; revised ... with respect to the lossy half-space characteristics for a given source and 2) to prove that the TDP sensitivity is statistically dispersed with respect to the characteristics of the source, hence leading to a non-negligible uncertainty. In particular, in the context of SAR assessment, available sensitivity analyses are very limited in scope, addressing such canonic sources as dipoles [START_REF] Douglas | Enabling the use of broadband tissue equivalent liquids for specific absorption rate measurements[END_REF] and normally impinging plane waves, which are not always representative of actual ones. Although the proposed approach is based on a closed-form analysis, the results presented here do not assume a far-field configuration for the radiation source, which is a key point for their practical use.

II. TRANSMISSION THROUGH DIELECTRIC INTERFACES

The scenario here considered is briefly depicted in Fig. 1. A planar interface Σ 0 separates two half-spaces, one containing a radiating antenna and the other one a lossy homogeneous dielectric. The first half-space is assumed to be air, but any other non-magnetic medium could be considered, whereas the second one is characterized by the electric parameters σ and ǫ r , which can be summed up by means of the complex permittivity ǫ = ǫ r -jσ/(ωǫ 0 ).

Rather than working in the space domain, it is convenient to introduce a Plane-Wave Expansion (PWE), which allows to represent the electric field distribution over a plane as a superposition of plane waves or Plane-Wave Spectrum (PWS). Considering the electric field E(R) over a plane, its PWS Ê(K) is defined as [START_REF] Scott | The Spectral Domain Method in Electromagnetics[END_REF]:

Ê(K, 0) = Σ0 E(R, 0)e jK•R dR , (1)
where R = xx + y ŷ spans the entire field topography over the plane Σ 0 and K = k x x + k y ŷ is the spectral variable;

hereafter, capital hatted quantities will stand for PWSs. The rationale for this approach is that in the spectral domain it is much easier to represent the propagation of complex field topographies through stratified media. The PWS over planes at other distances z at the right of Σ 0 can be easily computed introducing the propagator P (K, z) = exp(-jk z z):

Ê(K, z) = Ê(K, 0)P (K, z) , (2)
with k z = k 2 -K 2 as the longitudinal propagation constant, and k the propagation constant for the considered medium. The norms used throughout this work are to be intended as L 2 norms applied to complex vectors, unless explicitly declared otherwise.

The PWE formalism can be defined just over the tangential components of the vectors involved; for this reason, we will decompose the electric field PWS as

Ê = Ê + Êz ẑ , (3)
having defined the tangential part of the PWS as Ê = Êx x + Êy ŷ. The overall PWS Ê can be derived by recalling that Gauss law in a homogeneous medium in absence of local sources reads as k • Ê = 0, i.e.,

Êz = - K • Ê k z , (4)
where k = K + k z ẑ is the propagation vector for a plane wave. Thus, (3) can be recast as

Ê =    1 0 0 1 - k x k z - k y k z    Ê = T Ê . (5)
Let us now consider the interface Σ 0 : for each plane wave impinging from the left side, there will be a transmitted plane wave. Imposing the continuity of the tangential components of the electric field through the air-dielectric interface, the PWS Ê inc (K) of the incident electric field over Σ - 0 (i.e., on the air side) and the PWS Ê tx (K) of the field over Σ + 0 (on the right side) are related by [START_REF] Scott | The Spectral Domain Method in Electromagnetics[END_REF]:

Ê tx (K) = Π 0 (K) Ê inc (K) , (6)
where Π 0 (K) is the spectral transmission operator for a planar interface, as recalled in Appendix I. Hereafter, the R and Kdependency of, respectively, the spatial and spectral functions will be omitted most of the time for the sake of simplicity.

The TDP is related to the transmitted field by

TDP = σ Ω E tx (R, z) 2 dV , ( 7
)
where Ω is the volume occupied by the lossy half-space. Thanks to Parseval theorem (as applied over Σ + 0 ), we can state that

TDP = σ Êtx 2 -2Im {k z } dK . (8)
We will now introduce the only two major assumptions used throughout this work: 1) the spatial bandwidth k BW of the electric field PWS, satisfies the condition k BW k 0 |χ|, where χ = √ ǫ is the dielectric contrast of the lossy medium and 2) the medium can be described as a lossy dielectric, i.e., σ ≪ ωǫ 0 ǫ r . The first requirement is necessary so that the propagator modulus |P (K, z)| be almost flat over the entire bandwidth of the source PWS, allowing the approximation k z ≃ k 0 χ; the second assumption is needed in order to simplify the final results, although closed-form expressions could be given even in a more general case. Under such conditions, [START_REF] Harville | Matrix algebra from a statistician's perspective[END_REF] becomes

TDP ≃ √ ǫ r ζ Êtx 2 dK = √ ǫ r ζ E . ( 9
)
where ζ is the free-space wave impedance, and having defined as in signal analysis the signal energy E = E tx (R) 2 L 2 as the square of the L 2 norm of the field distribution over Σ + 0 . A major issue in this approach is that whenever the dielectric interface is in the near-field region of the source, coupling mechanisms will likely ensue. For a given distance d between the source and the half-space, the field impinging on the dielectric interface may differ from the one that would be measured at the same distance in a free-space configuration. These phenomena notwithstanding, the approach here proposed can also be applied to near-field configurations, without any loss of generality. Indeed, the coupling between one source and a passive scatterer can be regarded as given by an infinite series of separate interactions. This can be referred to as the multiple-interaction paradigm, and it has been applied, among other scenarios, in the definition of probe-correction models for near-field measurement techniques [START_REF] Hansen | Spherical Near-Field Antenna Measurements[END_REF].

The basic idea is that the time-domain evolution of the coupling is naturally represented through a series of simple interactions, that can be divided into scattering, transmitting and receiving events. The actual field Ê inc impinging over the half-space is thus generally not the one that would be radiated in a free-space configuration, e.g. in a far-field scenario, but rather the superposition of all the multiple-interaction contributions, i.e. the steady-state field including the actual coupling. This approach is completely general, with the farfield scenario as a special case.

III. SENSITIVITY ANALYSIS

We proceed by deriving a sensitivity model for the TDP, as discussed in Section III-A.A; subsequently, this model will be used in order to investigate how the TDP sensitivity is linked to the PWS of the source (and ultimately to the field it radiates), proving in Section III-B.B that for a given configuration the sensitivity is, in general, not a deterministic value, but rather a random variable.

A. Perturbation model

The sensitivity of the TDP to variations of ǫ can be assessed by applying a perturbation approach to the transmission equation [START_REF] Hansen | Spherical Near-Field Antenna Measurements[END_REF]. While a perturbation of ǫ has a direct impact on the TDP, it will have a higher-order effect on the reflectivity of the medium interface, as seen from free-space. In other words, it will be hardly noticed by the source: it will be shown in Section VI that this idea is viable, as long as the dielectric contrast |χ| stays high. This conclusion is intuitively justified by the fact that the fundamental quantity in source-half-space interactions is the dielectric contrast between the two media, which is negligibly affected by the lossy medium perturbation; this phenomenon was already pointed out in [START_REF] Correia | Mobile Broadband Multimedia Networks: Techniques, Models and Tools for 4G[END_REF]. Therefore, it will be assumed that the total incident electric field E inc (R) fulfills the condition ∂E inc /∂ǫ ≪ ∂E tx /∂ǫ .

The second order Taylor expansion of (9) with respect to the complex dielectric permittivity ǫ reads:

δTDP TDP = δE E + δ √ ǫ r √ ǫ r + δ 2 E 2E + δ 2 √ ǫ r 2 √ ǫ r + δ √ ǫ r √ ǫ r δE E . (10)
with:

δ √ ǫ r √ ǫ r = 1 2 δǫ r ǫ r (11)
δ 2 √ ǫ r √ ǫ r = - 1 4ǫ r δǫ r ǫ r (12)
The TDP sensitivity depends on δE/E and δǫ r /ǫ r sensitivity terms; since the latter are trivial, we will focus our analysis on the former. Moreover, in the subsequent statistical analysis of the sensitivity, δǫ r /ǫ r introduces a mere offset, without accounting for the statistical dispersion of the TDP sensitivity.

We will detail here the proposed analysis method for the first-order term δE/E, the second-order being derived fundamentally in the same way. Knowing that the first-order differential of the norm of a complex vector v is related to the differential of a complex variable p by

d v(p) 2 = 2Re v H ∂v ∂p dp , (13)
where the apex H stands for the Hermitian transpose, we can write:

δE = 2Re δǫ ÊH tx ∂ Êtx ∂ǫ dK . (14)
Using ( 5) and ( 6), together with the hypothesis of a weakly sensitive source, the following quadratic form is obtained:

δE = 2Re δǫ Ê H tx D 1 Ê tx dK , (15)
having defined matrix D 1 as follows:

D 1 = T H D ′ 1 (16) D ′ 1 = ∂T ∂ǫ + T ∂Π 0 ∂ǫ Π -1 0 . (17)
A great part of the results dealing with the sensitivity will be expressed as functions of the transmitted PWS. The rationale for this choice is that the PWS of E tx is more regular than the one in air, so that its envelope is more easily described and reproduced. This is important for the subsequent statistical analysis. Moreover, the mathematical analysis is also simplified thanks to this formalism, since the properties of the different operators are more readily applicable. In any case, the results are not affected by this choice, since the PWS Êinc and Êtx are biunivocally related through the transmission operator.

All the derivatives involving the transmission matrices are reported in Appendix I. It is sensible to define the operator S : C 2×2 → C, hereafter referred to as the sensitivity operator:

S(A) = Ê H tx A Ê tx dK T Ê tx 2 dK . (18)
Therefore, we can now provide the following fundamental result:

δE E = 2Re {S(D 1 )δǫ} . ( 19
)
This formulation has the advantage of pointing out the fact that the sensitivity is of course dependent on the PWS of the excitation, but that in fact its behavior is controlled by a kernel, modeled through the matrix D 1 , that is common to every source. In other words, D 1 represents the transfer function between the PWS of the impinging wave and the sensitivity that its components will experience while being transmitted into the lossy medium.

The development applied to the linear term can also be applied to the quadratic term in [START_REF] Rhodes | On a fundamental principle in the theory of planar antennas[END_REF], the only difference being that it is now necessary to consider the second-order differential of the norm of a vector v:

d 2 v 2 = 2 ∂v ∂p dp 2 + 2Re v H ∂ 2 v ∂p 2 (dp) 2 , ( 20
) hence δ 2 E 2E = Re (δǫ) 2 S(D 2 ) + |δǫ| 2 S(D 3 ) , (21)
having introduced the following matrices:

D 2 = T H ∂ 2 T ∂ǫ 2 + T ∂ 2 Π 0 ∂ǫ 2 Π -1 0 + + 2 ∂T ∂ǫ ∂Π 0 ∂ǫ Π -1 0 ( 22
)
D 3 = D ′H 1 D ′ 1 . (23)
The model thus derived allows computing the TDP sensitivity for a given source, just requiring as an input the PWS of the field impinging on the lossy half-space. An interesting by-product of the proposed analysis is that the frequency dependency of the TDP sensitivity is dominated by the term δǫ = δǫ r -jδσ/(ωǫ 0 ), as seen in ( 19) and ( 21). The fact that this term contains a 1/f associated to the electrical conductivity term means that any model approximating it through a polynomial expansion is bound to need high-order terms. This is actually the case, as proved by the choice made in [START_REF] Douglas | Enabling the use of broadband tissue equivalent liquids for specific absorption rate measurements[END_REF], where a third-order frequency polynomial had been necessary. The proposed representation, being derived from a physical analysis, thus leads to a more effective definition of the sensitivity.

B. Statistical analysis

The previous results can be used for investigating the following problem: how does the TDP sensitivity depend on the PWS of the source? We will here describe the PWS Ê tx as a stationary random process and study how the energy This problem can be solved easily by means of (18), since all the quantities involved can be computed directly without any need of actual simulations or measurements. Nevertheless, (18) requires to impose a certain PWS distribution; it is possible to prove that the sensitivity actually depends just on the envelope of the PWS, since (18) involves quadratic expressions [START_REF] Harville | Matrix algebra from a statistician's perspective[END_REF] (see Appendix II). For this reason, the energynormalized spectrum envelope w(K) will be considered, as defined by:

w(K) = Ê tx 2 Ê tx 2 dK . ( 24
)
In the context of a statistical analysis, the sensitivity operator defined in (18) can be evaluated directly; the focus will be put on studying the first two statistical moments of the sensitivity for a given PWS energy-density spectrum w(K) . Since all the quantities in ( 19) and ( 21) are deterministic but for the sensitivity operator, the statistical analysis will focus on this last term.

In particular, the following approximation holds (see Appendix II) as long as the spatial bandwidth k BW of the normalized energy-density spectrum w(K) of the average source is smaller than k 0 |χ|:

S(A) ≃ λH (K) + λS (K) w(K) dK , (25)
where λH and λS are the arithmetic means of, respectively, the eigenvalues of the Hermitian and the skew-Hermitian parts Fig. 3. Eigenvalues of the Hermitian and skew-Hermitian parts of operator D 2 and D 3 (last picture) associated to the quadratic sensitivity term, at 0.9 GHz (solid line), 2.5 GHz (dash line) and 6 GHz (dotted line). The results here shown regard the cut ky = 0. of matrix A, where this last stands for any of the derivative operators D i . This result is noteworthy, since it states that the average of the sensitivity over all the possible sources depends exclusively on the properties of the lossy medium and the average spectral content of the source.

The dispersion of the TDP sensitivity around its average can be assessed by studying the eigenvalues of the three derivative operators D i . Indeed:

min λ H wdK ≤ S(A H ) ≤ max λ H wdK (26) min λ S wdK ≤ S(A S ) ≤ max λ S wdK , ( 27
)
so that, the distance between the eigenvalue pairs is a direct measure of the dispersive nature of the sensitivity.

In order to test these ideas, we will hereafter consider a planar SAR assessment configuration as a practical application, since it is an interesting example of high-contrast lossy configuration; the nominal electrical parameters of the lossy [START_REF]Human Exposure to Radio Frequency Fields From Handheld and Body-Mounted Wireless Communication Devices -Human Models, Instrumentation and Procedures, Part 2: Procedure to Determine the Specific Absorption Rate (SAR) in the Head and Body for 30 MHz to 6 GHz Handheld and Body-Mounted Devices Used in Close Proximity to the Body[END_REF].

medium are thus imposed by international standards, such as [START_REF]Human Exposure to Radio Frequency Fields From Handheld and Body-Mounted Wireless Communication Devices -Human Models, Instrumentation and Procedures, Part 2: Procedure to Determine the Specific Absorption Rate (SAR) in the Head and Body for 30 MHz to 6 GHz Handheld and Body-Mounted Devices Used in Close Proximity to the Body[END_REF], and are reported in Table I. Therefore, the associated eigenvalues are shown in Fig. 2 and3, as computed for the three half-space configurations of Table I. The spectral variable was normalized with respect to the free-space wave-number k 0 . Matrices D 1 and D 2 were expanded into their Hermitian and skew-Hermitian components, as explained in Appendix II, whereas D 3 is already Hermitian.

Figures 2 and3 show that each pair of eigenvalues coincides around the origin K = 0, i.e., for normal paraxial incidence. The L 2 -integrable discontinuities around k 0 are due to the branch singularity in the spectrum of Green's function [START_REF] Rhodes | On a fundamental principle in the theory of planar antennas[END_REF]; a similar behavior is present around k 0 |χ|, fundamentally due to the same phenomena, but mitigated by losses. Conversely, it can be noticed that the distance between each pair of eigenvalues widens with the time and space frequency. This trend gets much stronger when getting beyond k 0 , i.e., when considering reactive components of the source PWS. More specifically, the skew-Hermitian components present an increased distance with respect to the purely Hermitian ones for a given K. Looking at (19), the skew-Hermitian component operates over the variations in the conductivity σ: this means that the reactive components of the source PWS will be affected in a more variable way in response to variations of the conductivity than equal-strength modifications of the permittivity. These simple conclusions already show that the sensitivity will present an higher degree of variability for nearfield sources, richer in reactive energy, than far-field ones, where no reactive component is available; this is coherent with the common-sense perception that the TDP sensitivity is not dispersed for a far-field configuration.

In order to assess these findings, the sensitivity operators were studied in several ways: 1) equation ( 25) was used for assessing the average sensitivity, 2) a deterministic estimate was given assuming a paraxial propagation, thus considering the approximation S(D)(K) ≃ S(D)(0); and finally 3) equation (18) was directly evaluated by considering a population of ten thousand random realizations for the transmitted PWS. For this last case, we considered a Gaussian energy-density spectrum with a -3 dB bandwidth k BW = k 0 , with a zero-mean Gaussian distribution. This choice is justified in the context of SAR assessment tests for telecommunications, where sources are usually not highly directive. Table II shows the results of this analysis. It turns out that the paraxial approximation is not very effective when compared to (25); nevertheless, it allows to introduce a very simple model for the average sensitivity, as reported in Section VI. Much worse is the fact that the paraxial approximation fails to acknowledge the intrinsic variability of the sensitivity with respect to the source PWS; the sensitivity is indeed collapsed into a deterministic value, rather than a statistical distribution. Statistical distributions for the five operators are shown in Fig. 4 for the same Gaussian PWS envelope, but for a changing spectral bandwidth k BW , for a frequency of 6 GHz. It appears that the case k BW = k 0 /2, i.e., with very little reactive components, is indeed well described by the paraxial approximation; but as soon as the reactive content increases, the sensitivity moves away and spreads. This trend is in accordance with the results presented in Fig. 2 and3 and the previous considerations about the distance between the eigenvalues. It is therefore not possible to neglect the variability of the sensitivity due to the source characteristics. The spread is quantified in Table II for k BW = k 0 , showing that considering a 95 % margin of confidence implies an uncertainty of ±20 % for the sensitivity, approximating the distributions as Gaussian ones.

IV. INCLUDING THE PRESENCE OF BUFFER LAYERS

In many practical cases the lossy dielectric may not face directly the air half-space, but buffer layers of complex dielectric permittivity ǫsi and thickness s i are interposed between the two half-spaces, for i ∈ [1, N b ] where N b is the number of buffer layers. Although the previous model was developed under a half-space assumption, it can be promptly extended to a stratified configuration; in the following discussion, we will assume that all the media have negligible losses with respect to the inner medium. Therefore, for the sake of computing the TDP, the only important quantity is the PWS Ê tx of the field transmitted into the lossy dielectric. This can be computed directly in the case of a stratified structure, by cascading the relationships reported in Appendix I. In order to simplify the notations, hereafter all of the quantities related to the electric field will refer to the tangential component. In this Section we first investigate how the presence of a buffer layer affects the TDP, in particular by extending the analysis of the role of the reactive energy of the source. Then, we focus on the modification of the TDP sensitivity.

A. Impact on the TDP

The PWS Êtx of the field transmitted into the lossy medium can be regarded as having been generated in a configuration involving just two media, i.e., the right half-space made of the lossy dielectric and the left one in air, as in the original description. To this end, an equivalent source needs to be introduced, much in the same way as for a Thevenin equivalent circuit. For the case of a single additional layer as depicted in Fig. 5, we can write:

Êtx = Π 2 R -1 Π 1 e -jkzss
Êinc ,

with

R = 1 + e -2jkzss Γ 1 Γ 2 . ( 29
)
Matrices Γ 1 and Γ 2 are the reflection matrices at the two dielectric interfaces, as given in Appendix I, whereas k z0 and k zs are the longitudinal propagation constants for, respectively, the free-space and the shell media. The PWS transmitted in the non-stratified configuration was given by ( 6):

Êtx = Π 0 Êinc e -jkz0s . (30)
By comparing (28) and (30), it is possible to define an equivalent PWS Êeq impinging from the air side (with no shell) as :

Êtx = Π 0 Êeq e -jkz0s = Π 0 Ξ Êinc e -jkz0s , (31)
having introduced the correction matrix

Ξ Ξ = Π -1 0 Π 2 R -1 Π 1 e -j(kz2-kz0)s . ( 32
)
Under a weakly-sensitive source assumption, the impact of the shell on the TDP can be assessed by computing the ratio of the signal energy E s for the shell configuration and the original one E 0 , that are related to the definition of the TDP given in [START_REF]Human Exposure to Radio Frequency Fields From Handheld and Body-Mounted Wireless Communication Devices -Human Models, Instrumentation and Procedures, Part 2: Procedure to Determine the Specific Absorption Rate (SAR) in the Head and Body for 30 MHz to 6 GHz Handheld and Body-Mounted Devices Used in Close Proximity to the Body[END_REF], while expressing all the spectral quantities as functions of the transmitted PWS:

TDP s TDP 0 = E s E 0 = Π 0 ΞΠ -1 0 E tx,0 2 L 2 E tx,0 2 L 2 , ( 33
) 0 0.5 1 1.5 2 2.5 3 3.5 4 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 k x /k 0 || Π 0 Ξ Π 0 -1
|| 0.9 GHz 2.5 GHz 6.0 GHz Fig. 6. The norm of matrix Ξ for a 2 mm thick buffer with ǫs = 4. The inner lossy dielectric are tissue-equivalent liquids as described in Table I.

where E tx,0 is the PWS transmitted with no shell present. The previous expression can be bound as

E s E 0 ≤ Π 0 ΞΠ -1 0 2 w(K)dK . ( 34
)
This formulation can be usefully employed in assessing the overall effect of a buffer layer on the TDP. Indeed, the matrix Π 0 ΞΠ -1 0 accounts for the modification of the transmissivity of the PWS from the air side to the lossy medium in the presence of a buffer layer; again, its effect is weighted by the energy-normalized spectrum w(K).

The norm of Ξ is related to

Ξ ≤ Π 1 Π 2 R -1 ; (35)
the first two norms are bound to have a finite value, due to their physical meaning: they account for the transmissivity. On the other hand, R -1 may be unbounded, since it represents standing-wave phenomena. This can lead to a strong increase in the TDP as soon as a resonance can be physically instated inside the shell region. An example of the behavior of Π 0 ΞΠ -1 0 is given in Fig. 6, for the case of SAR tissue-equivalent liquids with a shell 2 mm thick, with a relative permittivity ǫ s = 4. These results clearly prove that the TDP is indeed modified by the presence of the buffer, especially (as expected) at those frequencies where its thickness is comparable with the wavelength. These conclusions support and complete the findings reported in [START_REF] Onishi | Influence of phantom shell on SAR measurement in 3-6 GHz frequency range[END_REF], [START_REF] Christ | The dependence of electromagnetic far-field absorption on body tissue composition in the frequency range from 300 MHz to 6 GHz[END_REF].

Fig. 6 shows that the shell affects the TDP in two ways: 1) in the visible region (with respect to air), the shell operates as an impedance transformer, providing a better matching between the wave-impedances in air and in the lossy dielectric; 2) in the reactive region, it allows the transmission of more reactive energy ( K > k 0 ) if 1 < ǫ s < ǫ r , since it behaves as a lens focusing PWS components towards the direction normal to the dielectric interfaces. It is this last phenomenon that gives the strongest contribution to the TDP modification. As a consequence, modification of the reactive part may have a non negligible effect on the topography of the electric field inside the lossy medium.

As in the case treated in Section III-B.B, the impact of a shell on the TDP depends on the reactive content of the impinging PWS. This same approach yields the results shown in Fig. 7 : although the average modification of the TDP is not strongly affected, the more reactive the source, the more statistically dispersed is the TDP. It is fundamental to bear in mind that the scenario here considered presents a perfect knowledge of the electrical characteristics of the lossy medium.

The effectiveness of the bound given in ( 34) is proven in Table III, where for the 6 GHz configuration, with k BW = k 0 , the maximum modification of 1.41 well represents the dispersiveness shown in Figure 7.

B. Impact on the TDP sensitivity

The equivalent-source approach can now be applied to derive the sensitivity of the TDP in the stratified configuration, yielding

∂ Êtx ∂ǫ = ∂Π 2 ∂ǫ R -1 Π 1 e -jkzss Êinc + -Π 2 e -jkzss R -1 Γ 1 ∂Γ 2 ∂ǫ R -1 e -2jkzss Π 1 Êinc . (36)
As shown in Section VI, the sensitivity of the reflection matrices is one order of magnitude smaller than that of the transmission ones. Hence, we can claim that:

∂ Êtx ∂ǫ ≃ ∂Π 2 ∂ǫ R -1 Π 1 e -jkzss Êinc . (37)
Comparing this result to the sensitivity obtained in the first place with no buffer yields:

∂ Êtx ∂ǫ ≃ ∂Π 0 ∂ǫ Ξ ′ Êinc e -jkzss , (38)
where

Ξ ′ = ∂Π 0 ∂ǫ -1 ∂Π 2 ∂ǫ R -1 Π 1 e -j(kz2-kz0)s . (39)
Plugging (31) and (38) into ( 14) provides a tool for computing the new linear sensitivities as in the case with no shell. The fundamental quantity that dominates the modification of the sensitivities is the product of the norms of the two matrices Ξ and Ξ ′ . An example of the spectral behaviour of Ξ ′ is shown in Fig. 8: as opposed to the modification of the TDP, the norm of the first derivative of the transmitted PWS is generally not strongly affected when compared to the case with no shell, especially for sources with a poor reactive content. Again, the statistical distributions related to the source PWS have been computed and are shown in Fig. 9. These results give a better insight into the modification of the linear sensitivities; indeed, the sensitivity to the conductivity (imaginary part of S) is more strongly affected than the one to the permittivity. In particular, its average value decreases for a more reactive source and it spreads over a quite larger support.

The derivation of a closed-form upper-bound is not feasible for the modification of the linear sensitivities. Nevertheless, the following bound provides some information:

|S(D 1H ) + jS(D 1S )| s |S(D 1H ) + jS(D 1S )| 0 ≤ ≤ Π 0 ΞΠ -1 0 2 Ξ ′ 2 w(K)dK , (40)
where the indices stands for "shell" and "no shell". The statistical results for the modification of the TDP and its sensitivities are summarized in Table III for k BW = k 0 , together with the value of the correction factors evaluated for a paraxial propagation and the upper bounds given in (34) and (40). As for the results shown in Table II, though the paraxial model provides a fairly good estimate for the average modifications, it is unable to account for the dispersion they induce. Moreover, it cannot explain neither the stronger impact for more reactive sources, nor the shift in the linear sensitivity related to the conductivity (see Fig. 9).

For the sake of brevity, we will not derive here the secondorder sensitivity in the case of a shell. Nevertheless, the same approach can be extended to include such analysis. We can conclude that, in the case of SAR applications, the presence of a shell has not a fundamental impact on how the TDP reacts to modifications in the permittivity ǫ r , whereas the sensitivity to the conductivity σ is more strongly spread and lowered. Moreover, the shell can indeed strongly increase the transfer of energy between the two half-spaces, as well as lead to a further uncertainty in the evaluation of the TDP for a given near-field source. 0.9 GHz 2.5 GHz 6.0 GHz Fig. 10. The maximum eigenvalue of the matrix T H T, as computed for the three SAR configurations in Table I.

V. NEGLECTING THE NORMAL COMPONENT

There have been several discussions about whether the normal component Êz be negligible or not [START_REF] Kuster | Energy absorption mechanism by biological bodies in the near field of dipole antennas above 300 MHz[END_REF]. In other words, is it possible to approximate Ê with Ê ? This same assumption is used throughout this paper. An answer to this question can be given by recalling that [START_REF] Golub | Matrix Computations[END_REF]

Ê 2 = T Ê 2 = Ê H T H T Ê ≤ ν(K) Ê 2 , (41)
where ν(K) is the maximum eigenvalue of the matrix T H T. This allows writing

1 ≤ Ê(K) 2 Ê (K) 2 ≤ ν(K) . (42)
A useful rule of thumb for deciding whether the normal component can be neglected involves again the hypothesis k BW k 0 |χ|: in this case, the normal component can be neglected. As an example validating this statement, the upper bound ν(K) is shown in Fig. 10 for several frequencies, as computed for standard tissue-equivalent liquids.

Back in the spatial domain, the impact of E z on the TDP can also be assessed; indeed:

1 ≤ TDP TDP ≤ ν(K)w(K)dK , (43)
having defined TDP as the TDP obtained by neglecting E z . The band-limitedness of the PWS means that the upper-bound for the error on the TDP for SAR configurations must be lower than 4 % for a Gaussian-envelope PWS with k BW = k 0 and 1 % for k BW = k 0 /2, whereas it can be as high as 20 % for k BW = 2k 0 .

VI. PARAXIAL APPROXIMATION

In the previous analysis, no assumption was made on the distance between the source and the dielectric interface, so that the reactive components (in air) of the PWS were considered as contributing to the TDP. Nevertheless, the propagator P (K, z) in air imposes a very strong attenuation to the reactive parts, as expressed by the attenuation-per-wavelength A dB/λ :

A dB/λ K, z λ = 40π log(e) z λ K k 0 2 -1 . ( 44
)
As an example, the PWS components over the spectral radius K = 1.1k 0 experience an attenuation of 25 dB after a propagation distance of one wavelength. In many practical situations such a distance is not enough to justify a far-field analysis; nevertheless, although the field topography has not yet assumed its far-field configuration, a distance of λ is enough to enforce a PWS with almost the same envelope w(K) that would be obtained in a truly far-field analysis. Thus, from the TDP point of view, these two configurations coincide. Recalling that the eigenvalues of all the sensitivity operators vary very slowly over the air active region (see Fig. 2 and3), they can be fairly approximated by considering their values in the origin. This means that even in a configuration that does not satisfy the far-field criteria, a paraxial analysis can be applied.

Under a paraxial approach, the transmission operator Π 0 is just a scalar function, given by

Π 0 = 2 1 + χ , (45)
to be compared with the exact expression (55). It is thus necessary to define the sensitivities of the transmission terms:

1 Π 0 ∂Π 0 ∂ǫ ≃ - 1 2χ(1 + χ) (46) 1 Π 0 ∂ 2 Π 0 ∂ǫ 2 ≃ 1 + 3χ 4χ 3 (1 + χ) 2 , (47)
and, finally, the sensitivities of the signal energy E as required by ( 10)

δE E ≃ -Re δǫ χ(1 + χ) (48) δ 2 E E ≃ 1 2 δǫ χ(1 + χ) 2 + + 1 2 Re 1 + 3χ χ 3 (1 + χ) 2 (δǫ) 2 . ( 49
)
These results can be applied in practice as long as the reactive content of a source is negligible, thus providing a simple way of predicting the sensitivity of the TDP together with [START_REF] Rhodes | On a fundamental principle in the theory of planar antennas[END_REF].

It is also interesting to consider the reflection coefficient Γ d at the air-dielectric interface as seen from air and its relative sensitivity:

Γ d = χ -1 χ + 1 (50) 1 Γ d ∂Γ d ∂ǫ ≃ 1 |χ(χ 2 -1)| . ( 51
)
This last result can be used to explain the low sensitivity of the incident field to drifts in the lossy medium. To this end, let us consider the power P R = P av (1 -|Γ in | 2 ) radiated by a lossless antenna; Γ in is here the steady-state reflection coefficient as seen from the antenna feed-line, and it is thus affected by the reflections over the dielectric interface. Considering the free-space reflection coefficient Γ a at the antenna input port, Γ in can be written as

Γ in = Γ a + αΓ d , (52)
where α is related to the transmitting properties of the antenna [START_REF] Kerns | Plane-Wave Scattering-Matrix Theory of Antennas and Antenna-Antenna Interactions[END_REF]; the finite directivity of the antenna implies that |α| < 1. These considerations lead to the following result:

1 P R ∂P R ∂ǫ = 2|αΓ in | 1 -|Γ in | 2 ∂Γ d ∂ǫ (53)
and finally to

1 P R ∂P R ∂ǫ < |Γ in | 1 -|Γ in | 2 1 |χ(1 + χ) 2 | . (54)
Therefore, the relative sensitivity of the power radiated by the antenna is approximately bounded by |Γ in |/|χ| 3 , which is at least one order of magnitude smaller than the signalenergy sensitivity. This result proves that, unless in a very near-field configuration, the power radiated by the antenna is less strongly affected by changes in the lossy half-space characteristics than the TDP, as reported in [START_REF] Correia | Mobile Broadband Multimedia Networks: Techniques, Models and Tools for 4G[END_REF].

VII. CONCLUSIONS

We have introduced a spectral approach for the analysis of the sensitivity of the TDP to drifts of the electrical properties of a lossy half-space. The definition of a sensitivity operator and of the derivative matrices has shown that their eigenvalues distribution leads to a clear understanding of complex phenomena, such as the dispersiveness of the TDP sensitivity and in particular the fact that the sensitivity to the conductivity is more critical. The same approach was extended to the case of a stratified structure, in order to investigate how a lossless shell modifies the TDP and its sensitivity.

In all these scenarios, the fundamental role played by the reactive content of the source PWS was highlighted, pointing out how it gives rise to a statistically dispersive behaviour of the TDP and its sensitivity to drifts in the electrical parameters of the lossy medium. Hence, the very idea of characterizing the sensitivity in a deterministic way, independently from the source, is not physically sound, especially for near-field sources. These results should thus lead to a better understanding of the phenomena involved in near-field configurations, such as in SAR applications.

APPENDIX I DEFINITION OF THE SPECTRAL TRANSMISSION OPERATOR AND RELATED DERIVATIVES

The transmission and the reflection operators, respectively Π and Γ, for a dielectric interface between two media are defined as [START_REF] Scott | The Spectral Domain Method in Electromagnetics[END_REF]:

Π = 2X -1 (Y 1 + Y 2 ) -1 Y 1 X (55) Γ = X -1 (Y 1 + Y 2 ) -1 (Y 1 -Y 2 )X (56)
Y i = - 1 ωµ 0 k z,i k 2 i -k 2 x -k x k y -k x k y k 2 i -k 2 y (57) X = 0 1 -1 0 , ( 58
) being k z,i = k 2 i -K 2 , i ∈ [1, 2].
By applying the derivative chain-rule, the first derivative ∂Π/∂ǫ is given by:

∂Π ∂ǫ = ∂Π
∂k 2 ∂ǫ = k 1 2 √ ǫ . (63)
In the same way, the second derivative is given by:

∂ 2 Π ∂ǫ 2 = ∂ 2 Π ∂k 2 2 ∂k 2 ∂ǫ 2 + ∂Π ∂k 2 ∂ 2 k 2 ∂ǫ 2 , ( 64
)
where

∂ 2 Π ∂k 2 2 = 2X -1 B 2 ∂Y 2 ∂k 2 B ∂Y 2 ∂k 2 + - ∂ 2 Y 2 ∂k 2 2 BY 1 X (65) ∂ 2 Y 2 ∂k 2 2 = Y 2 K 2 k 4 z,2 + 1 2 ωµ 0 k z,2 2 K 2 k 2 z,2
+ 1 (66)

∂ 2 k 2 ∂ǫ 2 = - k 1 4ǫ 3/2 . ( 67
)
Concerning the derivatives of the T operator, we get: The sensitivity operator defined in (18) operates over a matrix A ∈ C 2×2 . We consider at first the fact that any matrix can be decomposed into the sum of an Hermitian part A H and a skew-Hermitian one A S . Hence, the matrix A H is orthonormal, i.e., with eigenvalues λ Hi ∈ R, as well as jA S ; this last claim implies that the eigenvalues λ Si of A S are purely imaginary. Furthermore, thanks to their symmetry properties the diagonalization matrices are orthonormal too, so that one can write:

∂T ∂ǫ =   0 0 0 0 k x k y   1 2 k 2 1 k 3 z,2 (68)
A H = X H Λ H X H H (70) A S = X S Λ S X H S , (71)
where Λ H,S are diagonal matrices containing the eigenvalues of, respectively, the Hermitian and the skew-Hermitian component of A. Imposing v = Ê tx , we are able to write the integrand of the numerator of (18) as: In order to compute the statistical average of S(A), the following functions need to be studied

v H Av =
w i 2 = |e i | 2 m |e m | 2 dK , (77)
where e stands for e H or e S . This function represents the energy distribution of the PWS normalized to its total energy; it can thus be regarded as a function describing an envelope. Only one component of the PWS, as represented over the diagonalized basis, is considered, according to the value of the index i. It is reasonable to assume that the two components of this function are identically distributed; therefore, their averages are identical too. Hence, the energy-density spectrum w = w i , ∀i can be defined, yielding

S(A) ≃ λH (K) + λS (K) w(K) dK , (78)
where λH and λS are the arithmetic means of, respectively, the eigenvalues of matrices A H and A S . The fact that the same envelope has been used for the Hermitian and the skew-Hermitian parts is due to the unitary property of the diagonalization matrices. The energy content of the PWS is therefore not modified passing from one basis to the other.

No assumption has been made in order to come to this result; it is therefore independent of the probability distribution of the energy-density spectrum.

Fig. 1 .

 1 Fig. 1. The planar configuration used for the spectral analysis.

Fig. 2 .

 2 Fig.2. Eigenvalues of the Hermitian and skew-Hermitian parts of operator D 1 associated to the linear sensitivity term, at 0.9 GHz (solid line), 2.5 GHz (dash line) and 6 GHz (dotted line). The results here shown regard the cut ky = 0.

- 0 .Fig. 4 .

 04 Fig.[START_REF] Douglas | Enabling the use of broadband tissue equivalent liquids for specific absorption rate measurements[END_REF]. Statistical distributions of the sensitivity operators for the Hermitian and skew-Hermitian components of the three derivative matrices D i , as computed for a Gaussian-distributed PWS with a Gaussian-envelope energydensity function with -3 dB spatial bandwidth k BW /k 0 equal to 0.5, 1 and 2. The working frequency is 6 GHz and the lossy medium is specified in TableI.

Fig. 5 .

 5 Fig.5. A stratified configuration with a lossless buffer layer of thickness s.

Fig. 7 .

 7 Fig. 7. Statistical distributions of the modification factor Es/E 0 for the TDP, due to the presence of a dielectric shell for a 6 GHz SAR configuration.

Fig. 8 .

 8 Fig. 8. The norm of matrix Ξ ′ for a 2 mm thick buffer with ǫs = 4.

Fig. 9 .

 9 Fig.9. Statistical distributions of the modification of the linear sensitivity operator due to the presence of a dielectric shell for a 6 GHz SAR configuration. The real part is related to the sensitivity to the dielectric permittivity ǫr, while the imaginary part deals with the sensitivity to the conductivity σ.

and ∂ 2 T ∂ǫ 2

 22

TABLE II SENSITIVITY

 II OPERATORS FOR A GAUSSIAN-ENVELOPE PWS, WITH BANDWIDTH k 0 . ALL THE RESULTS MUST BE MULTIPLIED BY A FACTOR 10 -4 .

  MM AND ǫs = 4, FOR A GAUSSIAN-ENVELOPE PWS WITH k BW = k 0 . THE BOUND FOR THE TDP

			TDP	S(D 1H ) S(D 1S )
		•	1.008	0.87	0.80
	0.9 GHz	• 2 paraxial	0.0034 1.01	0.0066 0.92	0.02 0.86
		upper bound	1.009	(1.07)	-
		•	1.04	0.88	0.79
	2.5 GHz	• 2 paraxial	0.010 1.03	0.0074 0.92	0.034 0.88
		upper bound	1.03	(1.13)	-
		•	1.19	0.99	0.89
	6.0 GHz	• 2 paraxial	0.03 1.16	0.02 1.02	0.05 0.98
		upper bound	1.41	(1.33)	-
			TABLE III
	MODIFICATION OF THE TDP AND ITS LINEAR SENSITIVITIES IN PRESENCE
	OF A DIELECTRIC SHELL, WITH s = 2 MODIFICATION IS GIVEN BY (34), WHEREAS THE VALUES IN
	PARENTHESES CORRESPOND TO THE BOUND (40).

  e H H Λ H e H + e H S Λ S e S = Hi |e Hi | 2 + Si |e Si | 2 , (72)having introduced the representation of the vector v into the new basis given by the diagonalization matrices A H and A S as: Hi and λ Si are the eigenvalues of, respectively, matrices A H and A S .As shown in Section V, matrixT H T ≃ 1 for any source with k BW k 0 |χ|. Recalling that X H X H H = X S X H S = 1, we can state that v H T H Tv ≃ i |e Hi | 2 =

			e H = X H H v		(73)
			e S = X H S v ,		(74)
	where the terms λ i	|e Si | 2 ,	(75)
	leading to
	S(A) ≃	i	λ Hi	m	|e Hm | 2 dK |e Hi | 2	dK +
	+		λ Si		|e Si | 2
		i

i λ i λ m |e Sm | 2 dK dK . (76)


Statistics of the Performance of Time Reversal in a Lossy Reverberating Medium

Andrea Cozza * SUPELEC, Département de Recherche en Électromagnétisme, 3 rue Joliot-Curie, 91192 Gif-sur-Yvette, France.

(Dated: October 20, 2009)

It has been proposed [J. Derosny, Ph.D. Thesis, Université Paris VI, 2000] that the performance of time reversal at recreating a coherent pulse in a strongly reverberating medium is directly proportional to the number of resonant modes M actively taking part at the transmission of energy. This idea is here tested against experimental results, showing that as soon as losses are taken into account, the quality of the focused pulse is a sublinear function of M , leading to a saturation phenomenon that was previously unacknowledged. This is here proven to be caused by mutual coupling between lossy resonant modes, thanks to a statistical modal description of the transmission of signals through the medium. Closed-form relationships are proposed for the first two moments of the pulse signal-to-noise ratio, linking them to the occupied bandwidth, the number of active modes and the degree of resonance of the medium. These formulae, supported by experimental and numerical results, prove that the performance of time reversal can be affected by a strong statistical dispersion. The proposed analysis also predicts that time reversal is a self-averaging process when applied to a reverberating medium, thus allowing the use of models developed in an ensemble-average framework.

I. INTRODUCTION

The interest of time-reversal techniques has been demonstrated in the fields of acoustics as well as electromagnetics, giving rise to a host of applications as diverse as pulsed-energy focusing in complex media [1], imaging techniques [2], and selective focusing [3]. Among these, one of the most surprising features of time reversal is its ability to ensure the transmission of coherent pulses through reverberating media [1,4]. In the context of this paper, we will consider a reverberating medium as a generally inhomogeneous medium where the propagation of electromagnetic or acoustic energy is strongly constrained into a finite volume. This region of space can be identified by an ideally closed surface imposing highly reflective boundary conditions, a configuration often referred to as a cavity. The provision of a finite volume does not exclude the existence of small apertures, through which a limited amount i(t) y(t) transducers reverberating medium FIG. 1: Synoptic of a generic communication system embedded into a reverberating environment. A pulse y(t) is to be focused at the receiver location by feeding the transmitter with an ad hoc signal i(t) defined by means of time-reversal techniques. * andrea.cozza@supelec.fr of the energy can leak out of the reverberating environment leading to a small perturbation of its behaviour assimilable to an energy loss. Coherently with this scenario, the media filling the cavity, as well as its boundary surface, will be regarded as lossy, introducing power dissipation along wave propagation. Following this description, a cavity is characterized by an infinite but countable set of resonances, associated to the eigenmodes of Helmholtz equation when the reflectivesurface boundary conditions are enforced [START_REF] Morse | Methods of Theoretical Physics[END_REF].

Typically, for a pulse transmitted within such a medium, the signal received would be dominated by a long non-coherent tail, made up of a large number of echoes of the original pulse; conversely, use of time-reversal techniques allows a predominantly coherent transmission of the pulse. Applications of this property to electromagnetics range from new signal-forming schemes for telecommunications in multipath channels [START_REF] Derode | [END_REF] to the generation of high-intensity local fields for device/material testing [START_REF] Moussa | IEEE EMC Symposium[END_REF]. Although the basics for the physical interpretation of time-reversal in such context are known [4], there is no available model allowing to predict the statistical behaviour of this technique when used in a reverberating medium, and in particular how its performance depends on the relative position of the receiver-transmitter pair within the system (see Fig. 1). In particular, the analyses found in the literature focus just on the mean asymptotic performance, without giving any hint of its statistical dispersion. Furthermore, to the best of our knowledge, these models assume the system to be lossless [4].

In this paper we fill this gap by proposing a study of the performance of time reversal for more realistic scenarios, assessing how the signal-to-noise ratio (SNR) of received pulses evolves while changing, on the one hand, the positions of the transducers and, on the other hand, the properties of the medium. This is done by studying the statistics of the performance, in particular by proposing closed-form expressions for the first two moments of the SNR of the received pulse. These results should allow predicting more thoroughly the way time reversal behaves in a reverberating medium, especially thanks to the knowledge of the variance of the SNR: indeed, this is a fundamental piece of data for ensuring, within a certain confidence margin, a given performance for any po-sition of the receiver and transmitter. Moreover, by acknowledging the existence of non-negligible loss mechanisms, the proposed model predicts phenomena that were not previously highlighted, such as the fact that losses lead to a saturation of the SNR, because of the mutual coupling between resonant modes through the tails of their frequency spectra. A major point is the generality of these results, which are valid for any system obeying to a modal resonant description, be it acoustic or electromagnetic, while requiring a very limited number of simplifying assumptions.

The paper starts with the introduction of tools for assessing the quality of the received pulse, generalized to the case of lossy media by applying novel definitions with respect to [4]. These tools are then applied for the analysis of experimental results obtained in an electromagnetic reverberation chamber, pointing out some of the previously recalled phenomena that have not yet been acknowledged. A theoretical analysis based on a modal description is then proposed in Section III, leading to the first two statistical moments of the SNR; subsequently, we focus in Section IV on the asymptotic response, proving that thanks to the self-averaging properties of time reversal, its statistical description is asymptotically independent on the actual realization of the reverberating medium, and entirely describable through few global parameters. Numerical results, as well as experimental ones are presented in Section V, validating the accuracy of the proposed analysis. Finally, a simple heuristic interpretation of our findings is given in Section VI, providing a framework for intuitively understanding the reasons for the saturation of the SNR.

II. ON THE QUALITY OF RECEIVED PULSES

A. Mathematical tools for quality assessment

In this Section we are concerned by the use of time-reversal techniques for transmitting a coherent pulse to a given receiver placed in a reverberating medium. Hereafter, this will be supposed to be reciprocal. The configuration we deal with is depicted in Fig. 1: two transducers are placed within the medium, one acting as a transmitter and the other one as a receiver. Defining h(t) as the impulse response between the two transducers, it was shown in [4] that by transmitting the signal i(t) = x(-t) ⋆ h(-t), the received signal y(t) will be a fair replica of x(-t), even in a strongly reverberating medium. Due to this feature, time reversal has been proposed as a way of communicating through complex media, and in particular multipath channels, whose characteristics are well represented by reverberating media [START_REF] Holloway | [END_REF]. An example of pulse received for this setup is shown in Fig. 2, where x(t) is a cardinal sinus modulating an harmonic carrier. Indeed, it appears that the received signal is almost undistinguishable from the original one around its peak region, whereas it is affected by a stronger modification over the signal tails.

In the following, we will consider the pulse x(t) to have a spectral content comprised in the frequency range [f 1 , f 2 ], i.e., with a frequency bandwidth B T = f 2 -f 1 centered around the frequency f c = (f 1 + f 2 )/2. A total number of M reso- FIG. 2: An example of a pulse transmitted through a reverberating medium using time reversal, for the case of Q = 5000, M = 500, BT = 100 MHz, central frequency 1 GHz: the region around the peak of the pulse (a) and a wider perspective highlighting the residual noise distribution (b). The dark trace is the original pulse to be transmitted, whereas the light one is the signal actually received, affected by residual noise. The thicker line represents the equivalent noise ne(t). nant modes will be assumed to exist over this frequency range, with resonance frequencies given by the set {f k }. In a general way, a transfer function H(f ) in a reverberating medium can be expressed as a superposition of these modes, weighted by complex coefficients γ k = α k + jβ k [START_REF] Morse | Methods of Theoretical Physics[END_REF], i.e.,

H(f ) = M k=1 γ k (f )φ k (f ) , f ∈ [f 1 , f 2 ] . (1)
The response φ k (f ) of the k-th mode will be assumed, with no loss of generality, to be a Lorentzian function. By considering the main effect of losses to lead to a small perturbation of these functions, one can write [START_REF] Morse | Methods of Theoretical Physics[END_REF]:

φ k (f ) = f 2 k /Q k f 2 k (1 + j/2Q k ) 2 -f 2 , (2)
where Q k is the quality factor associated to the k-th mode, which thus has a -3 dB bandwidth B M,k = f k /Q k . The modal weights {γ k } are functions of the transmitter and receiver positions and of the spatial field distribution associated to each resonant mode.

In order to simplify the notations in the following analysis, we will consider the reference signal to be x(-t), so that the received one is rather linked to x(t). For the same reason, we will avoid delaying the time-reversed pulse, and consider a non-causal description, as this does not affect the final results. Bearing these definitions in mind, the fitness of the focused pulse can be assessed by computing the component s(t) = ρx(t) that is coherent with respect to x(t), leading to

y(t) = ρx(t) + n(t) , (3)
where n(t) is the residual noise due to the distortion of the pulse introduced by the non-flat transfer function H(f ) of the medium, with H(f ) the Fourier transform of h(t). Residual n(t) being orthogonal to x(t) by definition, ρ can be computed as

ρ = Y, X X, X = f2 f1 X|H| 2 X ⋆ df f2 f1 |X| 2 df , (4)
having applied Parseval equality, where X(f ) and Y (f ) are the Fourier spectra of the respective time-domain signals. The brackets stand for the projection operator. It is interesting to notice that because of the quadratic form in the previous result, ρ ∈ R + , so that the sign of the peak of the received pulse will always be unchanged. Defining the energy E G of a spectrum G(f ) as

E G = 2 f2 f1 |G(f )| 2 df , (5)
the energy E S of the coherent part is thus given by

E S = ρ 2 E X = 2ρ 2 f2 f1 |X| 2 df , (6)
whereas the energy E N of the residual noise is

E N = 2 f2 f1 |X| 2 |H| 4 df -E S , (7)
so that we can introduce the energy SNR Λ as

Λ = E S E N . (8)
While the energy ratio Λ will be extensively used in the rest of the paper, the ability of time reversal in transmitting coherent pulses is better assessed by means of the peak SNR Λ p as defined as follows [4]:

Λ p = s 2 (0) n 2 rms , (9)
having assumed that the peak of the received pulse occurs in t = 0. The quantity n rms is the root-mean-square (rms) value of n(t):

n 2 rms = lim T →∞ 1 T T /2 -T /2 n 2 (t)dt . (10)
Equation ( 9) thus measures how much the coherent part of the transmitted pulse stands out of the residual noise. Now, the use of the rms value is consistent only in the case of a lossless system, as done in [4], since in this case n(t) would have an infinite energy but a finite non-zero average power. Conversely, for the case of a lossy system, n rms would be equal to zero, since the noise has finite energy; as a matter of fact, the time-constant characterizing the decay of a mode with quality factor

Q k is τ k = Q k /(πf k ) = 1/(πB M,k
). We will assume that all the modes involved have the same bandwidth, and as a consequence the same time-constant τ = τ k , ∀k; this assumption is valid as long as B T /f c is sufficiently smaller than one, i.e., for configurations that cannot yet be regarded as wide-band, though not strictly narrow-band. Under such conditions, n(t) will also obey to a time-decay with constant τ ; therefore, we introduce the equivalent noise signal n e (t):

n e (t) = n 0 e -πfct/Q , ( 11
)
where Q is the average quality factor, and by imposing the same overall energy for the two noise signals, we get

n 2 0 = E N πf c Q . ( 12
)
The equivalent noise n e (t) behaves as a smoothed version of the actual noise n(t), maintaining the same overall timedecay, and thus the same average instantaneous power content. The example shown in Fig. 2 illustrates this approach. It is now possible to define the peak SNR by considering the equivalent instantaneous noise energy at the signal peak in t = 0,

Λ p = s 2 (0) n 2 e (0) = Q πf c χΛ , (13)
where χ = x 2 (0)/E X is a factor related to the shape of x(t). This simple relationship between the two SNR definitions allows focusing on the energy SNR, which is much simpler to compute in the frequency domain.

Thanks to these definitions, it is possible to predict the performance of time-reversal transmission for any pulse x(t), just by knowing the transfer function H(f ). Most remarkably, this just requires having access to the absolute values of the spectra X(f ) and H(f ).

B. Experimental investigations in a reverberation chamber

In general, the pulse SNR will depend on the relative position of the transmitter and the receiver within the system; in order to exploit the time-reversal technique for real-life applications, it is of paramount importance to be able to ensure that a given minimum SNR be respected for any transducers position, at least with respect to a certain confidence margin. In order to assess the variability of the SNR, we carried out experimental tests, by considering an electromagnetic reverberation chamber, with a fixed antenna acting as the transmitter and a linearly polarized electro-optical sensor (connected to an optical fiber) as a receiver. This last choice was imposed by the fact that a receiving antenna being moved inside the chamber would have changed its fundamental characteristics, modifying the sets {f k } and {γ k }, and as thus impairing the validity of any comparison. The very weak interaction ensured by the optical sensor means that it can be regarded as an almost ideal electric-field probe, minimizing the modification of the quantity being measured.

A total of 100 randomly chosen positions and orientations of the receiver were considered, measuring the respective transfer functions over six frequency ranges, centered at frequencies from 0.5 GHz to 3 GHz, by steps of 0.5 GHz. For each central frequency, two bandwidths B T were considered, namely 100 MHz and 200 MHz.

We assumed, for the sake of simplicity, x(t) to have a flat spectrum over the frequency range [f 1 , f 2 ]; this would be the case, e.g., for cardinal sinus pulses. The energy SNR was computed by means of ( 8), leading to the empirical statistical distributions shown in Fig. 3, and the first two statistical moments of Λ resumed in Table I, together with the average peak SNR Λ p . These results point to three important facts: 1) while the average energy SNR increases with the frequency, its progression slows down at the same time, converging on an asymptotic value close to one, 2) the actual value of Λ is strongly dispersed, depending on the position of the receiver and, 3) an increasing bandwidth has apparently little effect on the average energy SNR, whereas its dispersion decreases. To the best of our knowledge, these conclusions have never been discussed before, and no theoretical framework is available for interpreting them.

The only available model is the one proposed in [4] for a lossless configuration. It predicts a direct proportionality between the number of modes M and the average peak SNR Λ p . In order to apply this model to our experimental results, we estimated the number M w of resonant modes existing in the chamber for a given frequency range through Weyl's formula [9]

M w (f c ) ≃ 8πV c 3 f 2 c B T 1 + B 2 T 12f 2 c , ( 14
)
where V is the volume of the reverberation chamber and c is the speed of light in the medium filling it. Results obtained for the case of Supelec's chamber (V = 3.08 × 1.84 × 2.44 m 3 ) are shown in Table I. By comparing M w to the average Λ p , it is clear that their relationship is more complex, and characterized by a form of saturation of the performance, since even for large increases in M w , Λ p is barely affected.

The reasons for such a peculiar behaviour are to be sought in the existence of loss mechanisms in actual reverberation chambers; we will show in the next Section that this leads to mutual couplings between resonant modes, and ultimately to a saturation of the performance.

Concerning the statistical dispersion discussed in points 2) and 3), it has never been addressed before. It is worthwhile noticing that the orientation of the probe has little effect on the dispersion, since the cavity was over-moded for all the frequency ranges, apart for f c = 0.5 GHz. For such configurations, the field is statistically isotropic, and it presents very similar statistical properties along its three Cartesian components [10]. These results point out that apart from being able to explain the limitations of time reversal, it is of paramount importance to have a model predicting the statistical dispersion of the SNR. These are indeed the basic motivations of this paper.

fc (GHz) BT (MHz) mean Λ std Λ mean Λp Q Mw 0.

III. STATISTICAL MOMENTS OF Λ

As often done in statistical descriptions, the parameters of the model will need to be regarded as random variables. This approach is not just dictated by mathematics, but it comes with physical meaning. In particular, the resonance frequencies {f k } are indeed distributed over the bandwidth B T in a way that is hardly predictable, unless in canonical configurations. For the associated modal weights {γ k }, since describing the projection of the transducer characteristic response over the modal topographies, a modification in the position of the receiver or the transmitter leads to a modification of the {γ k }, so that a random position of the transducers implies a random set of modal weights. Moreover, the fact that the modal topographies, as well the excitation of the transducers are, in general, sign-changing functions, implies that the {α k } and {β k } should be treated as zero-mean random variables, and they will be assumed to be independent and identical distributed (iid). No further assumption will be necessary about the type of distributions.

In order to simplify the model, the φ k (f ) will be assumed to be frequency-shifted replica, with approximately the same bandwidth B M = f k /Q k , ∀k. This also implies that all the modal responses φ k (f ) have the same energy E φ . In principle, this assumption holds only when relatively narrow bands are considered, although the results shown in Section V prove that

Λ f c = 0.5 GHz f c = 1.0 GHz f c = 1.5 GHz f c = 2.0 GHz f c = 2.5 GHz f c = 3.0 GHz (a) (b) (c) (d) (e) (f
Λ (g) (h) (i) (j) (k) (l)
FIG. 3: Frequencies of occurrence for the energy SNR Λ as computed from experimental data measured over a bandwidth BT = 100 MHz (left column) and BT = 200 MHz (right column), centered around the frequencies fc = {0.5, 1, 1.5, 2, 2.5, 3} GHz (top to bottom). Each histogram was obtained from a population of 100 sample transfer functions as measured between a fixed transmitter antenna and a moving electro-optical probe.

this is not necessarily the case.

It is important to understand the physical role of the following statistical analysis, which aims at accounting for the impact of the random position of the transducers on the received pulse SNR. Indeed, equation (1), when coupled with equations ( 6)-( 8), leads to the definition of Λ as a random function, depending on the probability density functions (pdfs) of {γ k } and {f k }.

We start our analysis by considering a specific configuration for the reverberating medium, i.e., for a given set of known deterministic {f k }, whereas the {γ k } will be regarded as random variables. This scenario corresponds to the case of a single realization for the medium, while the positions of the transducers are let free to change, so that all the statistical moments will be conditional to the set {f k }. The ensemble behaviour of the SNR considering random {f k } will be studied in Section IV.

In order to simplify our analysis, but with no loss of generality in the conclusions, we will assume the modulus of X(f ) to be directly proportional to the characteristic function of the interval [f 1 , f 2 ], leading to χ = 2B T . This choice corresponds, e.g., to a cardinal sine excitation in the time-domain, modulating an harmonic carrier of frequency f c , as for the example shown in Fig. 2. Attention should be paid to the fact that the definitions of the SNR actually depends just on the modulus of X(f ), so that an infinite number of pulse shapes sharing the same spectral occupation would be characterized by the same SNR.

Following this assumption, (8) can be recast as

Λ = f2 f1 |H| 2 df 2 B T f2 f1 |H| 4 df - f2 f1 |H| 2 df 2 . ( 15
)
Thus Λ is entirely defined by the properties of the random function |H(f )| 2 over the frequency-range [f 1 , f 2 ]. In order to study the statistical properties of Λ, we introduce the auxiliary random variables W i ∈ R + , as defined as

W i = f2 f1 |H| 2i df , (16)
yielding

Λ = W 2 1 B T W 2 -W 2 1 . (17)
The rationale for introducing these auxiliary variables is that the statistical moments of Λ cannot be expressed as a direct function of the {γ k } and {f k } moments. Nevertheless, the moments of W i can be linked more easily to those of {γ k } and {f k }; an estimation of the moments of Λ can then be given by linearizing (17) around the W i ensemble averages [START_REF] Khuri | Advanced Calculus with Applications in Statistics[END_REF]. For the sake of simplicity, the following convention is introduced:

W i = E[W i |{f k }] . (18)
Applying this approach to the average of Λ conditional to a given realization {f k } yields

E[Λ|{f k }] ≃ W 2 1 B T W 2 -W 2 1 . (19)
In the same way, the conditional variance σ 2 Λ|{f k } can be approximated as

σ 2 Λ|{f k } = E[Λ 2 |{f k }] -E[Λ|{f k }] 2 ≃ J T ΣJ , ( 20
)
where Σ is the covariance matrix of the random vector [W 1 , W 2 ] T . The column vector J is the Jacobian of Λ as computed with respect to W 1 and W 2 , evaluated at (W 1 , W 2 ):

J = B T W 1 B T W 2 -W 2 1 2 2W 2 -W 1 . (21)
A higher-order estimate of Λ could be given, but the resulting expression would be quite unwieldy without delivering considerable improvement in the final accuracy. For the same reason, we just consider the first two moments of W 1 and W 2 . Thanks to the following expansion

|H(f )| 2 = M k=1 |γ k | 2 |φ k (f )| 2 + + 2 M k=1 M m=k+1 Re {γ k γ ⋆ m φ k (f )φ ⋆ m (f )} (22)
and assuming the {γ k } to be independent from the {f k }, while recalling the hypothesis of all the modes having the same energy E φ , as defined in ( 5), we can write

E[W 1 |{f k }] = 2M µ 2 E φ (23) E[W 2 |{f k }] = 2M (µ 2 2 + µ 4 )E |φ| 2 + + 16µ 2 2 M k=1 M m=k+1 E φ k φm (24)
having introduced the moments µ i

µ i = E[α i k ] = E[β i k ] . (25)
Equation ( 24) differs from (23) in a fundamental aspect, i.e., the presence of the mutual energies E φ k φm shared between each couple of modes of the system. This term can be shown to be the source of the limitations of the SNR as the modal density M/B T increases enough to lead to nonnegligible interactions between the modes.

The same type of analysis was carried out for the elements

Σ ij = E[(W i -E[W i ])(W j -E[W j ])
] of the covariance matrix, but this led to too complex expressions, especially for Σ 22 . We would rather propose approximate results, where the mutual-energy terms are neglected:

Σ ≃ 4M ν , (26)
having defined the elements of ν as

ν 11 = E 2 φ (µ 4 -µ 2 2 ) (27)
ν 12 = E φ E |φ| 2 (µ 2 µ 4 + µ 6 -2µ 3
2 ) (28)

ν 22 = E 2 |φ| 2 (µ 8 + 4µ 2 µ 6 + µ 2 4 -2µ 4 2 -4µ 2 2 µ 4 ) (29)
It is worth noting that although interactions between modes have been neglected in Σ, the Jacobian in (21) takes them into account. It will be shown in Section V that when applied as an input to (20), these expressions provide a good estimate of the variance of Λ, and as thus they are a useful tool in stating the uncertainty that affects time-reversal performances in a reverberating medium.

IV. AVERAGE ASYMPTOTIC PERFORMANCE

The formulae presented in the previous Section were derived considering a given deterministic set of resonant frequencies {f k }, and as thus ( 19) and (20) depends, in principle, on the actual realization of {f k }. In fact, this dependence subsists only in the mutual energies E φ k φm in (24).

In this Section, we consider the performance of timereversal when averaged over all the realizations of {f k }, hence related to its general trend rather than for a specific configuration, proving that under certain conditions, the statistics of a single realization are well approximated by the simpler ensemble statistics.

To this end, let us consider the ensemble average of Λ with respect to the random set {f k }; thanks to the linearization of Λ, it will suffice to carry out this averaging over the W i , leading to:

E[W 1 ] = 2M µ 2 E φ (30) E[W 2 ] = 2M (µ 2 2 + µ 4 )E |φ| 2 + + 8µ 2 2 M (M -1)E[E φ k φm ] . (31)
The double sum in ( 24) is thus simplified by introducing the average mutual energy E[E φ k φm ]. The result in ( 24) and (31) would then be identical if the following condition were satisfied:

2 M (M -1) M k=1 M m=k+1 E φ k φm ≃ E[E φ k φm ] . (32)
This requirement corresponds to assuming the system to be ergodic, approximating the ensemble average over all the realizations with the average carried out over the set of mutual energies within a single realization. The strong law of large numbers [START_REF] Leon-Garcia | Probability, Statistics, and Random Processes for Electrical Engineering[END_REF] states that, if the system is ergodic, the left hand of (32) converges in probability to the ensemble average of the mutual energy; therefore, for a sufficiently high number of active modes M one gets

E[Λ|{f k }] ≃ E[Λ] . ( 33
)
This phenomenon, often referred to as self-averaging, had already been experimentally highlighted in [START_REF] Derode | [END_REF], although in a different context, as one of the most interesting features of time reversal, and it implies that its performance in transmitting coherent pulses in a reverberating medium is asymptotically independent from the actual realization of the set of the resonance frequencies {f k }, underpinning the robustness of this technique. An example highlighting this property is given in Fig. 4, where empirical pdfs are shown for 10 different realizations {f k }, for three values of M , namely 50, 100 and 200, with f c = 2 GHz, B T = 200 MHz and Q = 1000. For each realization of {f k }, ten thousand sets of modal weights {γ k } were considered. It appears that indeed the pdfs converge toward the ensemble average as M increases, even for such small values of M .

We can now write

E[Λ] = 2M E 2 φ B T (1 + µ 4 /µ 2 2 )E |φ| 2 + 4(M -1)E[E φ k φm ] -2M E 2 φ .
(34) It is clear that the behaviour of E[Λ] could be easily predicted should the three energy terms be known. As a matter of fact, a closed-form expression can be given for E[E φ k φm ] by exchanging the order of integration:

E[E φ k φm ] = f2 f1 f2 f1 |φ k (f )| 2 p(f k )df k 2 df , (35)
having exploited the fact that the {f k } are iid random variables. Noticing that the inner integral does not depend on the actual frequency f at which it is computed, as long as

B T ≫ B M , (35) becomes E[E φ k φm ] = B T f2 f1 |φ k (f c )| 2 p(f k )df k 2 , ( 36
)
and by assuming a uniform distribution for the {f k } over the bandwidth B T , this yields

E[E φ k φm ] = E 2 φ B T . ( 37
)
The two remaining energies E φ and E |φ| 2 can also be ex-pressed in closed-form as

E φ k = f k π 2 Q k Q 2 k + 1 ≃ π 2 B M (38) E |φ k | 2 = f k π 4 Q 3 k (Q 2 k + 5) (Q 2 k + 1) 3 ≃ π 4 B M , (39)
so that (34) can be restated in a simpler form

E[Λ] = M M + (1 + µ 4 /µ 2 2 ) B T 2πB M . ( 40
)
Hence E[Λ] is linearly dependent on M at low modal density levels, whereas it converges to an asymptotic value for an higher M . Thus, in lossy reverberating media, the potential gain obtained by increasing the number of active modes (i.e., increasing the central frequency f c ) is put in jeopardy by the coupling existing between lossy resonant modes, with E[Λ] converging to a fixed value. Interestingly, this asymptotic value is simply equal to one. The physical significance of this result will be given a simplified explanation in Section VI.

Equation ( 40) is remarkably simple, and it shows that a handful of global parameters is sufficient for an accurate prediction of the quality of the received pulse. It is worth noting that the central frequency does not appear explicitly, as a consequence of the identical-mode assumption. These results also point to the fact that the most fundamental quantity for understanding the phenomena behind pulse focusing in a reverberating medium is B T /B M . This quantity will be hereafter referred to as N s , for reasons that will be made clear in Section VI, yielding

E[Λ] = M/N s M/N s + (1 + µ 4 /µ 2 2 )/(2π) . ( 41
)
This reformulation states that the average performance is entirely predicted by means of the ratio M/N s . As soon as M N s the marginal gain brought by the availability of new modes is increasingly reduced, leading to a saturation for higher M .

V. MODEL VALIDATION

In order to check the accuracy of the proposed description, we considered numerical simulations, by synthesizing random realizations of transfer functions, thanks to (1). The rationale for this approach is the possibility to closely monitor the number of modes M , their quality factor, and so on. Indeed, as recalled later in this Section, experimental validations are impaired by the impossibility to assess the exact number of modes taking part to the transmission.

Thanks to the fact that (40) is not directly dependent on the central frequency f c , but rather on the bandwidths B T and B M , the validation can be carried out at any value of f c . We set for f c = 2 GHz, with a varying bandwidth B T and several average quality factors Q. Random complex weights γ k

FIG. 5: Empirical pdfs for the energy SNR Λ, depending on the number M of active modes and the quality factor Q of the medium. The six curves presented in each picture correspond to M = {100, 250, 500, 1000, 2500, 5000}, respectively, from left to right.

are drawn accordingly to normally distributed α k and β k ; the energy SNR for the transfer functions thus obtained are subsequently computed thanks to [START_REF] Holloway | [END_REF]. Contrary to the assumption of a constant f k /Q k , the φ k (f ) were assumed to have a constant quality factor Q = Q k , ∀k. Therefore, the modal responses φ k (f ) will not be identical as assumed in the model derivation.

The first tests aimed at showing how the energy SNR is distributed as M and N s = B T /B M vary, and is more of a qualitative investigation. A bandwidth B T of 50 MHz was chosen, while four values of Q were tested, ranging from 5000 to 40000. The number of modes varied from 100 to 5000. For each set of global parameters ten thousand realizations were generated, in order to establish empirical pdfs; these are shown in Fig. 5, and illustrate quite clearly that: 1) increasing losses tend to saturate the energy SNR faster, as M increases, 2) decreasing losses slow the saturation down, but reduce the average energy SNR, as the length of the residual-noise tails increases and so does the noise energy, and 3) the SNR experiences a standard deviation that is far from negligible when compared to the average value, although, as Q increases, the dispersion appears to decrease. The trend in the simulated pdfs recalls that of the experimental ones shown in Fig. 3. Even though the energy SNR Λ decreases with Q, the peak SNR Λ p increases monotonically, since the relationship between Λ and Q is actually sublinear. This fact is to be expected intuitively, and it also confirms the trend predicted by the model proposed in [4], since

lim Q→∞ E[Λ p ] = M , (42)
again, in the case of a cardinal-sinus pulse and gaussian statistics for the modal weights. Quantitative validations were then carried out by considering f c = 2 GHz, Q = {1000, 2500, 5000, 10000, 20000}, and a varying bandwidth B T = {200, 400, 600} MHz. The number of modes spanned the values 100 to 5000, and a population of 500 random realizations {γ k } per configuration was considered, each configuration representing just a single realization of {f k }. The results thus obtained for the average value and the standard deviation of Λ are shown in Fig. 6 and Fig. 7, respectively. Fig. 6 shows that the average value is predicted within a few percent points as long as N s 200. Indeed, (40) is unable to predict any energy SNR greater than one; this is actually not due to a bad estimate of the moments of the auxiliary variables W i , which are indeed precisely estimated in all of the considered tests, since (30) and (31) do not involve any approximation. This rather points directly to the conditions that are necessary for applying the linearization in (9), implying that it is necessary for the condition N s 200 to be fulfilled. We investigated the possibility of including the Hessian term in the expansion ( 19), but we dropped this option, since it brought no tangible improvement, thus implying that the SNR as a non linear function of W i would require terms higher than quadratic ones.

The standard deviation depicts a rather different scenario. Expression (20) neglects any modal interaction in the covariance matrix Σ, but it includes them through the use of W i . For this reason, (20) is expected to underperform as soon as the modal interactions get more important, i.e., as the modal density M/B T increases. The results in Fig. 7 support these ideas, showing that (20) is a very good estimate of the standard deviation, as long as it has not yet attained its maximum value σ max . After this point, (20) is no more a valid estimate, but the actual standard deviation gets to a plateau fairly approximated by σ max . In general, this value needs to be computed numerically, but for the case of modal weights distributed as Gaussian random variables, the number of modes M σ for which the standard deviation reaches its maximum value can be approximated by

M σ ≃ 6 π N s . (43)
Knowledge of the saturation point allows extending the validity of (20) over the entire range of values of M , i.e.,

σ 2 Λ ≃ (J T ΣJ)(M ) M ≤ M σ 1.7π N s M > M σ . (44)
Therefore, the maximum standard deviation goes like 1/ √ N s , whereas it is inversely proportional to N s for M ≪ M σ . The former conclusion explains the behaviour previously highlighted, with the standard deviation decreasing when B T and/or Q increase.

These numerical validations prove the effectiveness of the asymptotic models, even for a relatively low number of modes and with no ensemble averaging in {f k }. This implies that the ergodic assumption formulated in the previous Section does indeed hold. The greatest limitation in the proposed models is the need for (9) to be well approximated by its tangent plane over the range of values spanned by the W i , requiring N s 200.

The last validation is a tentative experimental one. As recalled at the beginning of this Section, the exact number of modes excited in a reverberating system is usually not known.

Therefore a direct validation is not feasible; nevertheless, it is current practice in electrical engineering to assumed that a linearly polarized antenna placed in an electrically large reverberating chamber will excite most of the modes existing over the frequency range of emission of the antenna. As a consequence, Weyl's formula is often used as a reference. Hence, we computed the moments of the energy SNR predicted by our model, considering a number of modes equal to three fractions of the estimate M w given by Weyl's formula (14).

The results are resumed in Table II, together with the quantity N s : the range of variation of the SNR is very well identified, both for the average value and the standard deviation, and the experimental results are consistently approached when considering a number of modes close to 0.9M w . Furthermore, as expected from the numerical validation, as soon as N s decreases towards 200 the experimental average Λ goes beyond one; in this case, the model will underestimate the statistical moments.

The fact that considering the same fraction of M w over the six frequencies leads to good results, strongly reduces the odds that this accuracy be a random result; we thus consider that 90 % of the available modes were indeed effectively excited. The only exception is for f c = 0.5 GHz and B T = 200 MHz; in fact, the transmitting antenna had a cut-off frequency around 450 MHz, so that of the 200 MHz pulse to be received, it actually transmitted only three-quarters of the signal spectrum, hence exciting roughly three-quarters of the available modes. By taking into account this fact, the actual number of modes to be considered is rather 3/4 • 0.9M w ≃ 2/3M w : indeed, the results agree.

Overall, it appears that the average of Λ is hardly affected by an increase in B T . Actually, this is predictable, since both N s and M w are linearly dependent on B T , so that (41) is not modified. Conversely, the peak SNR will increase proportionally to B T . At the same time, the standard deviation is sensitive to an increasing B T . This was predicted in (44), and the reduction of a factor 1/ √ 2 subsequent to a doubling B T is indeed well confirmed by the experimental results. These findings are of the utmost importance should time-reversal techniques be used for pulse transmission.

VI. AN HEURISTIC INTERPRETATION

We will here try to give an interpretation of the reported phenomena from a more physical, yet approximate, point of view. To this end, let us recall that the maximum value attained by the peak SNR Λ p , as long as N s 200, is simply given by N s /(2π). It is thus not dependent on the actual number of modes M , but rather to a, usually, much lower quantity.

Let us look at N s from a different perspective: knowing that B M is the average bandwidth of the frequency response of each mode, N s states the maximum number of modes that could be placed one after the other over the bandwidth B T . The energy SNR corresponding to this configuration is equal to one, and it corresponds to the best efficiency time-reversal can provide in concentrating energy in the coherent part rather than in the residual noise. This fact can be used to give a simple intuitive interpretation, by introducing the idea of a number N s of available slots, to be occupied by the actual number of active modes. Although simplistic, this vision of the spectrum as a quantified space makes sense. Hence, each mode introduces a further degree of freedom only if it can be allocated to a free slot; otherwise, it will be lost, just leading to a different modal weight, as a function depending on the weights of the modes previously allocated to the same slot. Therefore, the performance of time reversal is not related to the actual number of active modes, but rather to the number of slots being used, which could thus be regarded as an effective number of modes or degrees of freedom, or M e , leading to an efficiency and, ultimately, to an energy SNR equal to M e /N s . The allocation of a mode to a specific slot being a "rare" event, this random process can be modelled by a Poisson law, with mean M/N s . The mean number of occupied slots, and thus the effective number of modes M e , is thus simply given by M e = N s (1 -e -M/Ns ) ,

and the related energy SNR

E[Λ] ≃ M e N s = 1 -e -M/Ns , (46)
highlighting the dominant role of the quantity M/N s , as previously shown in (41). We could thus dub the quantity M/N s as the modal slot occupancy: it defines completely the SNR and is sufficient for predicting the performance of time reversal in any configuration. The validity of this reasoning is proven in Fig. 8. Indeed, for a low number of modes (with respect to N s ) the results predicted by (41) and (46) correspond fairly well. For higher slot occupancies, (46) saturates faster, since this model is incapable of acknowledging the partial superposition of two modes, something that would just lead to a partial loss of a degree of freedom. In spite of this over-simplification, this approach yields results consistent with those predicted by (41), while providing a simple framework for understanding the SNR saturation phenomenon.

VII. CONCLUSIONS

This paper has addressed the main phenomena underlying the quality of pulses received by a transducer as transmitted through a reverberating medium, when using time-reversal techniques. The quality of the received pulse has been analyzed with respect to global parameters identifying the properties of the medium, according to a modal description. Having included loss mechanisms, it was proven that the sharing of energy between finite-bandwidth resonant modes is at the origin of the limitations in the SNR of the received pulse. A statistical approach has led to general results based on very few assumptions, mainly that of a sufficiently "wide-band" configuration with B T /B M 200: the developed model predicts correctly the first two statistical moments of the SNR, acknowledging its non-negligible statistical dispersion. Although mainly based on a mathematical approach, the physical meaning of these results were explained in plain terms by introducing a simplified heuristic description, proving that the SNR is in fact limited by the finite number of degrees of freedom available in a lossy reverberating system. These results should be useful for both the design of experiments and the interpretation of their results, and pertain to any type of wave propagation problem in a reverberating environment.

FIG. 4 :

 4 FIG. 4: Empirical conditional pdfs p(Λ|{f k }) for Q = 1000, fc = 2 GHz and BT = 200 MHz, obtained for an increasing number of active modes. The thicker curve stands for the ensemble average of each group of realizations.

FIG. 6 :

 6 FIG. 6: Validation tests for E[Λ]: numerical results obtained from 500 realizations (dots) and the values predicted by (40) (solid lines). All the presented results were computed for fc = 2 GHz and for Q={1000, 2000, 5000, 10000, 20000}, while the corresponding Ns are displayed.

FIG. 7 :

 7 FIG. 7:Validation tests for the standard deviation of Λ, for the same configurations as in Fig.6. The dashed lines correspond to the values predicted by (20), but considering ensemble-averaged Wi, whereas the thick solid ones highlight the model results up to the curve local maximum, for M < Mσ.

FIG. 8 :

 8 FIG.8: Comparison between the mean energy SNR, as predicted by the modal approach (41) (solid line) and the slot occupancy description (46) (dashed line).

TABLE I :

 I Statistical moments of the energy SNR Λ and average peak SNR Λp, as computed from the experimental data presented in Fig.3. The approximate number of modes Mw was computed by means of Weyl's formula (14), whereas the average quality factor Q was directly estimated from the time constant of the residual noise n(t).

	5	100 200	0.37 0.074 0.30 0.051	187 300	3900	322 644
	1.0	100 200	0.75 0.12 0.72 0.080	272 130	5700	1287 2574
	1.5	100 200	0.91 0.14 0.90 0.10	243 481	6300	2896 5792
	2.0	100 200	0.96 0.13 0.95 0.12	192 380	6300	5149 10297
	2.5	100 200	1.02 0.13 1.00 0.10	187 365	7200	8045 16090
	3.0	100 200	1.02 0.14 0.98 0.10	188 358	8700	11584 23169

TABLE II :

 II fc (GHz) BT (MHz) experiments mod. 2Mw/3 mod. 3Mw/4 mod. Mw Ns Mw/Ns Experimental validation against the results presented in TableI: mean values are given directly, while standard deviations are in parenthesis. The results computed by means of (40) and (20) were obtained considering a number of active modes equal to 2Mw/3, 3Mw/4 and Mw, due to the uncertainty on the actual M . The reliability of the estimates can be tested by checking the condition Ns 200.

	0.5	100 200	0.37 (0.074) 0.30 (0.082) 0.33 (0.083) 0.39 (0.083) 785 0.30 (0.051) 0.30 (0.058) 0.33 (0.058) 0.39 (0.058) 1570	0.41
	1.0	100 200	0.75 (0.12) 0.70 (0.097) 0.73 (0.097) 0.78 (0.097) 570 0.72 (0.080) 0.70 (0.069) 0.73 (0.069) 0.78 (0.069) 1140	2.3
	1.5	100 200	0.91 (0.14) 0.88 (0.11) 0.90 (0.095) 0.88 (0.080) 0.89 (0.080) 0.92 (0.080) 840 0.89 (0.11) 0.92 (0.11) 420	6.9
	2.0	100 200	0.96 (0.13) 0.95 (0.013) 0.95 (0.13) 0.96 (0.13) 315 0.95 (0.010) 0.95 (0.092) 0.95 (0.092) 0.96 (0.092) 630	16
	2.5	100 200	1.02 (0.13) 0.97 (0.14) 1.00 (0.12) 0.97 (0.097) 0.97 (0.097) 0.98 (0.097) 576 0.97 (0.14) 0.98 (0.14) 288	28
	3.0	100 200	1.02 (0.14) 0.98 (0.14) 0.97 (0.10) 0.98 (0.096) 0.98 (0.096) 0.98 (0.096) 580 0.98 (0.14) 0.98 (0.14) 290	40