© The Author(s), Published by EDP Sciences, 2025

1 Introduction

Reverberation has a long history in acoustics. Even before scientific times, the word was used in Medieval French to describe the repercussion of voice in a room [1]. As early as 1835, Dr. Reid realised “that any difficulty in the communication of sound in large rooms arises generally from the interruption of sound produced by a prolonged reverberation” ([2], pp. 14–16). He later went on to renovate the House of Commons in London, where he successfully applied his experimental findings ([3], pp. 19, 145).

Most acousticians, however, consider that reverberation history started with W.C. Sabine's study of the lecture hall in the Fogg Art Museum of Harvard, where he conducted at the end of the 19th century extensive experiments to assess reverberation and sound absorption [4]. This led him subsequently to derive his famous formula, which still remains the basis of Architectural Acoustics.

Scientists quickly realised, however, that the Sabine formula was not a panacea because it gives a finite reverberation time for total absorption at all boundaries. Moreover, experimentally measured absorption coefficients can exceed one, that is, absorption can be more than total. Solving the absorption problem, as it was then called, became a challenge for acoustical engineers.

Since Sabine's pioneering days, many acousticians have proposed improved reverberation formulas. In a recent Internet publication [5], no less than 11 reverberation formulas are listed with references: Sabine [6], p. 117), Eyring [6], p. 116), Millington-Sette [7], Eq. (7.58)), Zhang [8], Kuttruff ([7], Eq. (7.64)), Arau-Puchades [9], Fitzroy-Eyring ([10], p. 176), Fitzroy-Sabine ([10], p. 176), Fitzroy-Kuttruff [11], Fitzroy-Neubauer ([12], p. 295), and Kang-Orlowski [13]. There are probably more formulas in the literature. Common to most of these formulas is the lack of a proper underlying theory that justifies them on rigorous mathematical and physical principles.

In 1992, the author proposed a generalized Sabine formula that developed reverberation time over a power series of the reflection coefficient on the boundaries of an ergodic enclosure [14]. Building on the work of Joyce [15] and Kuttruff [6], pp. 103–116), this generalized formula considered the long-term distribution of reflections along one single ray as a Lévy process, for which Paul Lévy ([16], pp. 180–186) developed the complete characterization. As a consequence, the decay process reduces to the analytical continuation to negative real values of the characteristic function of a Lévy process.

The present investigation was started in the hope of obtaining formulas piloted by dimension. When this idea turned out to be unsuccessful, the profile of the probability distribution for short free paths was investigated, because it is well known from Fourier transform theory that the asymptotic behavior at large times is piloted by the short frequencies – here, the short paths. Indeed, this probability is proportional to path length for very short free paths, therefore vanishes for null paths, contrarily to most published distribution of free paths, for example in Kuttruff [6]. It was thus assumed that the Sabine formula, that relies on a strictly positive probability for null paths (see Sect. 2.1), could not be the most general reverberation formula and had to be replaced by a revised reverberation formula, the formulation of which would depend on the standard deviation of the free paths. This is indeed the path followed in this paper.

The present paper, therefore, revisits the author's development of 1992 in the case of 2D and 3D circular and rectangular ergodic enclosures that, according to Joyce [17], follow Lambert law reflections. Circular enclosures were selected because they were treated by Joyce [17]; and rectangular enclosures because of their relevance in Building Acoustics (most rooms are rectangular), but also because of a renewed interest from practitioners in reverberation formulas for rectangular rooms, as experienced from students’ internships. As explained by Kuttruff [6], reverberation is piloted by the first two statistical moments of the free path distribution between successive reflections. Analytical formulations of these two moments have been derived for ergodic 2D and 3D circular enclosures [17], as well as 2D and 3D rectangular enclosures [18, 19]. Here, these analytical results are presented only, leaving the full derivations to Appendices A and B and to the Supplementary material submitted with the paper, because of their length and complexity.

The paper starts with a presentation of the most relevant publications on reverberation theory, including the author's reverberation developments of 1992, which are compared to Lévy's original expression for a characteristic function and to Kuttruff's formula [6]. It also revisits Kuttruff's formula and shows that Kuttruff's original integration [20] leads to an expression that does not break down for large absorptions coefficients. It then turns to the analytical expressions of the mean free paths and mean quadratic free paths for 2D and 3D ergodic enclosures, both circular and rectangular, and to the corresponding formulas that can then be derived from Lévy's expression. The case of high mean absorption coefficients is also discussed, when Kuttruff's formula breaks down, as well as the case of different absorption coefficients on different boundaries and the case of angle dependent absorption, before briefly analysing how the theory breaks down in non-ergodic enclosures, such as rectangular enclosures with purely specular reflections.

Finally, the implications of this “revised theory” for absorption measurements in reverberation chambers and for auditorium design are briefly discussed.

2 State of the art

Among the many publications devoted to reverberation theory, those of Schroeder, Joyce and Kuttruff stand out. Note that in the following, the expression “free flight time” is restricted to the time span between two successive reflections along one trajectory; whereas the expression “waiting time” indicates the time span for several reflections. The free flight path is equivalent to the free flight time, but for a factor equal to the speed of sound c.

2.1 Schroeder, Joyce and Kuttruff

Schroeder was probably the first to realize that reverberation theory essentially depends on the probability distribution of the number of reflections along rays [21]. He concluded that “Sabine's decay formula […] emerges here not as an approximation to the Eyring formula for live room, but as the result of certain assumptions concerning the distribution of the number of reflections”. However, he went on to consider instead the distribution of decay rates, a position he held for many years [14].

Joyce published two papers on reverberation theory. In the first one [15], he proved that Sabine's theory, that is, an exponential decay, is only obtained in ergodic enclosures: it is then sufficient to follow one ray to obtain the reflection statistics, as all rays are equivalent. In the second one [17], he showed that Lambert reflection law at the boundaries is a sufficient – but not necessary – condition for obtaining ergodicity: direct computation of the mean free path recovers the theoretical value of 4V/S in the 3D case.

Following Schroeder's first idea, Kuttruff focussed on reflection distribution along one single trajectory. In a first Letter [20], he considered that the probability distribution P(t|n) of waiting times t along a ray undergoing n reflections tends toward a Gaussian distribution for large n; inverting the roles of the variables, he then considered that the same expression expresses the probability distribution P(n|t) of reflections along a ray during the time span t. He then obtained, but for a different factor, the expected value of Rⁿ as a function of time, where R is the mean energy-based reflection coefficient at the boundaries:

$$\overline{R^n}=\frac{exp\left(-t\frac{\sqrt{1-2{\gamma }^{2}ln(R)}-1}{\overline{t}{\gamma }^2}\right)}{\sqrt{1-2{\gamma }^{2}ln(R)}}$$(1)

where $\overline{t}$ is the mean free flight time between two successive reflections and ${\gamma }^2=\frac{{\sigma }^2}{{\overline{t}}^2}$ the relative variance, that is, the variance of the free flight times divided by the squared mean free flight time in order to obtain an adimensional quantity. In his next paper [22], Kuttruff considered a Γ-distribution for the free flight times between successive reflections (n=1):

$$P(t|1)=\frac{1}{\mathrm{\Gamma }\left(\frac{1}{{\gamma }^2}\right)}{\left(\frac{t}{{\gamma }^2\overline{t}}\right)}^{\frac{1}{{\gamma }^2}-1}exp\left(-\frac{t}{{\gamma }^2\overline{t}}\right)$$(2)

which leads to a different expression for the expected value of Rⁿ:

$$\overline{R^n}=exp\left(-t\frac{1-R^{{\gamma }^2}}{\overline{t}{\gamma }^2}\right)·$$(3)

Note that for t≈0 (short paths), P(t|1) in equation (2) is vanishing if γ²<1, very large if γ²>1, and roughly equal to 1 if γ²=1, the last case corresponding to Sabine's theory (see Eq. (3)). Hence the statement that the Sabine formula relies strictly on a non-zero probability for null paths in Section 1.

Kuttruff did not derive equation (3), but thanks to his well known radiation equation for energy transfer between diffusing boundaries, showed that the equivalent absorption coefficient in the case of uniform absorption on the boundaries is equal to:

$${\alpha }^{*}=\frac{1}{{\gamma }^2}\left(1-R^{{\gamma }^2}\right)$$(4)

which he developed in powers of lnR, stopping at the second order.

In his subsequent book [6], Kuttruff reverted to his first approach, which he simplified in order to recover the Taylor expansion of the equivalent absorption coefficient, up to the second order in lnR. From it, he obtained an approximation of the expected value of Rⁿ as:

$$\overline{R^n}=exp\left(t\frac{lnR}{\overline{t}}\left[1+\frac{{\gamma }^2}{2}lnR\right]\right)·$$(5)

2.2 Generalized reverberation formula

In [14], the author took an other approach. Instead of considering the distribution function of free flight times, he considered the distribution of reflections along one single trajectory in an ergodic enclosure since, due to the ergodic hypothesis, all rays are equivalent. Mixing, that states that two nearly parallel rays finally become totally uncorrelated, further ensures the independence of events that take place in successive instants. In other words, the process is infinitely divisible, that is, it is a Lévy process. As a consequence, as explained in Section 1, the author's generalized reverberation formula from 1992 was based on the expression of the characteristic function of a Lévy process, which generalizes Poisson processes. Quoting Kuttruff [6], the author then considered the probability distribution P(n|t) of reflections along a ray during the time span t. According to Lévy [16], pp. 180–186), the logarithm of its characteristic function, better known as the cumulant function, is expressed by:

$$\begin{array}{cc}\hfill \psi \left(z,t\right)& =ln\left({\displaystyle \sum _{n=0}^{\infty }}{e}^{\mathit{izn}}P(n|t)\right)\hfill \\ \hfill & =m(t)\mathit{iz}-g^2(t)\frac{z^2}{2}\hfill \\ \hfill & +{\displaystyle {\int }_{-\infty }^{+\infty }}\left(e^{\mathit{izu}}-1-\frac{\mathit{izu}}{1+u^2}\right)\text{d}N(t,u)\hfill \end{array}$$(6)

where m(t), g(t) and N(t, u) are continuous functions of t, null for t=0; and where N(t, u) is a non-decreasing function of u in each interval ]−∞, 0[ and ]0, ∞[, null at infinity, and satisfying for all t: ∫u²dN(t, u) finite on all finite intervals. In Lévy's original text, u represents the amplitude of a jump; in the present case, it can be assimilated to an “order of reflection” – see Section 2.3. From that expression and symmetry considerations due to ergodicity and mixing, since time averages and ensemble averages are equal and rays are totally uncorrelated after some time, m(t), g(t), and N(t, u) must be proportional to t, which leads to the inference that the decay rate takes the form:

$$\begin{array}{cc}\hfill \delta (t)& =-\frac{\partial \psi \left(-ilnR,t\right)}{\partial t}\hfill \\ \hfill & =-mlnR-g^2\frac{{ln}^{2}R}{2}\hfill \\ \hfill & -{\displaystyle {\int }_{-\infty }^{+\infty }}\left(e^{ulnR}-1-\frac{ulnR}{1+u^2}\right)\text{d}N(u)\hfill \end{array}$$(7)

where iz has been replaced by lnR, R still being the mean energy-based reflection coefficient. In [14], the author further assumed that g is null and u always positive, which led to the following formula for reverberation time:

$$\text{RT}=\frac{13.8}{-mlnR+{\sum }_u\left(1-R^u\right)N(u)}$$(8)

when u takes discrete positive values only. Note that the Sabine formula is obtained for the particular choice m=0 and a single u=1, which is independent of the reflection coefficient, as noted by Schroeder [21]; and the Eyring formula for the particular choice m≠0 and u=0, also independent of the reflection coefficient. Also note that Kuttruff's formula is equivalent to retaining the two first terms only in the cumulant function, assuming that n(u)=0 for all u, leading to a different expression of the reverberation time:

$$\text{RT}=\frac{13.8}{-mlnR\left(1+\frac{{\gamma }^2}{2}lnR\right)}$$(9)

with ${\gamma }^2=\frac{g}{m}$, which diverges for large values of R since lnR is always negative; and that he hinted at the general term in $\left(1-R^u\right)$ in his 1976 paper [22]. But his original expression for the expected value of Rⁿ (Eq. (1)) leads to a different expression for the reverberation time:

$$\text{RT}=\frac{13.8\overline{t}{\gamma }^2}{\sqrt{1-2{\gamma }^{2}ln(R)}-1}$$(10)

where $\overline{t}=\frac{1}{m}$ still is the mean free flight time between two successive reflections, which can be reduced to the general form of equation (8) since:

$$\begin{array}{cc}\hfill & -\frac{1}{\overline{t}{\gamma }^2}\left[\sqrt{1-2{\gamma }^{2}ln(R)}-1\right]\hfill \\ \hfill & ={\displaystyle {\int }_0^{+\infty }}\left(R^u-1\right)\frac{e^{-\frac{u}{2{\gamma }^2}}}{\sqrt{2\pi }\overline{t}\gamma u^{\frac{3}{2}}}\text{d}u.\hfill \end{array}$$(11)

The derivation of equation (11) is carried out in Appendix A.

2.3 Deriving terms in formula

The different terms in the reverberation time formula can be derived from the standard deviation of path fluctuations. Indeed, when m(t) and g(t) are both null and u takes only one value, P(n|t) takes the explicit form:

$$P_u(n|t)=e^{-\mathit{tN}(u)}\frac{{\left(\mathit{tN}(u)\right)}^n}{n!}·$$(12)

Here, u corresponds to an “amplitude” – or order – of reflection, which can be different from unity, n is the number of reflections of order u during time span t, so that the effective number of reflections, in the usual sense of the term, is equal to nu. From this expression, the probability of waiting time t, that is, the time to wait for observing the nth reflection of order u, is given by:

$$P_u(t|n)=e^{-\mathit{tN}(u)}\frac{{\left(\mathit{tN}(u)\right)}^{n-1}}{(n-1)!}N(u)$$(13)

since there were only (n−1) reflections during time t, the nth one occurring between instants t and t + dt. The mean value and the variance of this waiting time are therefore given by:

$$\langle t\rangle =\frac{n}{N(u)},{\sigma }_u^2(n)=\frac{n}{N^2(u)}·$$(14)

Taking into account the effective number of reflections, the mean waiting time $\overline{t}$ and its variance σ² between effective reflections are given by:

$$\begin{array}{cc}\hfill \overline{t}& =\frac{\langle t\rangle }{\mathit{nu}}=\frac{1}{\mathit{uN}(u)},\hfill \\ \hfill {\sigma }^2& =\frac{{\sigma }_u^2(n)}{\mathit{nu}}=\frac{1}{{\mathit{uN}}^2(u)}=\frac{\overline{t}}{N(u)},\text{and}\hfill \\ \hfill {\gamma }^2& =\frac{{\sigma }^2}{{\overline{t}}^2}=u\hfill \end{array}$$(15)

where γ² is the relative variance. In the case of stadiums, which are squares prolonged by two half-discs on two opposite sides, numerical simulations from Mortessagne et al. [23, 24] lead to N(u)=111.3 for u=0.576. Considering that the value of u is very close to $\frac{1}{2}$, the author refined the analysis by considering a mean value m, attested by further data sent by Mortessage and Legrand. It consists of modifying the effective number of reflections to take into account the influence of m during the mean waiting time 〈t〉: the effective number of reflections becomes equal to $\mathit{nu}+m\langle t\rangle =n\left(u+\frac{m}{N(u)}\right)$, so that the mean waiting time and variances between effective reflections become:

$$\begin{array}{cc}\hfill \overline{t}& =\frac{\langle t\rangle }{n\left(u+\frac{m}{N(u)}\right)}=\frac{1}{m+\mathit{uN}(u)}\hfill \\ \hfill {\sigma }^2& =\frac{{\sigma }_u^2(n)}{n\left(u+\frac{m}{N(u)}\right)}=\frac{\overline{t}}{N(u)},\hfill \\ \hfill {\gamma }^2& =u+\frac{m}{N(u)}·\hfill \end{array}$$(16)

In other words, the variance keeps the same expression as before, but not the relative variance. The data from Mortessagne et al. [23, 24] now lead to N(u)=111.3 and m=8.46 for u=0.5. Note that the Sabine formula corresponds to u=1, and that Mortessagne et al. [23, 24] compared their numerical simulations to Kuttruff's formula (9).

Table 1 presents the reverberation times calculated for different values of the mean absorption coefficient α=(1−R). They were computed with the data from Mortessagne et al. [23, 24], and compared with Sabine's, Eyring's, and Kuttruff's reverberation times. The latest one diverges for α>0.981, that is, for the last column of Table 1. The data correspond to a mean waiting time of 0.0156 s, or a mean free path of 5.3 m, that is, a stadium built on a square with 4.85 m long sides. The speed of sound was taken as 340 m/s. Note that Eyring's reverberation time RT_E corresponds to keeping the first term only in the denominator of Kuttruff's formula.

The last line but one in Table 1 corresponds to Kuttruff's revised formula (10). As for the last line, it corresponds to $u=\frac{{\ell }_{\mathit{TH}}}{l_{min}+{\ell }_{\mathit{TH}}}$, with $N(u)=\frac{l_{min}+{\ell }_{\mathit{TH}}}{\overline{t}{\ell }_{\mathit{TH}}}$, as will be explained in Section 3.3.1 in connection with path distribution in parallelepipeds.

Figure 1 presents the same data in graphic form. As can be seen in Table 1 and Figure 1, the last five reverberation formulas are equivalent within 5% accuracy up to α=0.5: they bundle together and remain half-way between the Sabine and the Eyring formulas that can be considered as bounds for reverberation times, as explained by Joyce [15]. Curves “uTH” correspond to equation (8) with m=0 and only one value of $u={\gamma }_{\mathit{TH}}^2=0.576$, with $N(u)=\frac{1}{\overline{t}{\gamma }_{\mathit{TH}}^2}$; and curves “uMFP” to the last line of Table 1. For larger absorptions, they spread more, especially for Kuttruff's formula that starts diverging above α>0.7 when it gives values that increase with absorption.

Figure 1 Reverberation times for the stadium. But for Sabine's and Eyring's formulas, all other reverberation time formulas give indistinguishable reverberation times.

Note that RT_u=0.576, or curve uTH, suffers from the same drawback as the Sabine formula: they give positive reverberation times for full absorption. Also note that Mortessagne et al. considered purely specular reflection in their stadium, which is ergodic and mixing by shape, hence demonstrating that the Lambert reflection law at the boundaries is not a necessary condition for obtaining ergodicity, as stated in Section 2.1.

3 Ergodic 2D and 3D enclosures

We now turn on to circular and rectangular enclosures. Complex analytical integration recovers the usual values for the mean free paths – πS/P for the 2D case and 4V/S for the 3D case, where V is the volume, S the surface, and P the perimeter of the enclosure – as demonstrated by Joyce [17]. In fact, the mean free paths were calculated as a verification of the integration procedures. Note that, in order to help comparison between different shapes of enclosures, the mean free path has been kept to the same value of 5.3 m as in the previous stadium.

3.1 Circular enclosures

In his paper on the effect of surface roughness on reverberation time, Joyce [17] explicitly developed the case of 3D spherical enclosures, and gave qualitative results for 2D circular enclosures. He was interested in the transition from random reflection – Lambert law – to specular reflection. Here, we only use the results for random reflection.

In the 3D case, the mean free path is equal to $\overline{\ell }=\frac{4}{3}R$, where R is the radius of the sphere, and the mean quadratic path to $\overline{{\ell }^2}=2R^2$, which gives a variance ${\sigma }^2=\frac{2}{9}{R}^2$, or a relative variance of ${\gamma }^2=\frac{\overline{{\ell }^2}-{\overline{\ell }}^2}{{\overline{\ell }}^2}=1/8=0.125$. In the 2D case, these values are respectively $\overline{\ell }=\frac{\pi }{2}R$, $\overline{{\ell }^2}=\frac{8}{3}{R}^2$, ${\sigma }^2=\left(\frac{8}{3}-\frac{{\pi }^2}{4}\right)R^2$, or a relative variance of ${\gamma }^2=\frac{\overline{l^2}-{\overline{l}}^2}{{\overline{l}}^2}=0.0807\approx 1/12.4$. The calculations are straightforward and are reproduced in Appendix B. Note that mean free paths and mean quadratic paths must be converted to time with the help of the speed of sound c, taken as 340 m/s in the calculations.

Table 2 gives the same reverberation times as Table 1, except for the last line which is not relevant in this case, calculated for the disc for different values of the mean absorption; and Table 3 the same data for the sphere. In both cases, the mean free path was kept to 5.3 m, as in Table 1. Thus, Sabine's and Eyring's reverberation times are the same as for the stadium (Tab. 1). However, the values of u have to be adapted to each case: for the disc, u=0.0807 leads to N(u)=794.3; but choosing u=1/13=0.0769 leads to the same value N(u)=794.3 with m=0.0608. For the sphere, u=0.125 leads to N(u)=512.8; and choosing u=1/9=0.111 leads to the same value N(u)=512.8 with m=0.0608.

Figure 2 presents the same data as Tables 2 and 3 in graphic form. Notice again that all reverberation times take values between the Sabine and the Eyring formulas, which again can be considered as bounds for reverberation times. In these cases of low relative variances, the last four formulas remain very close to each other up to α=0.9. Note that Kuttruff's formula only diverges for R=(1−α)<1.72 × 10⁻¹¹ for the disc, and R<1.125 × 10⁻⁷ for the sphere, which are almost zero. Therefore; Kuttruff's reverberation time is not discernible from the others, but for Sabine's reverberation time.

Figure 2 Reverberation times for uniform absorption coefficient in circular enclosures. Top: 2D disk; bottom: 3D sphere.

As pointed out by Joyce [17], the advantages of the disc and the sphere are that computing mean free paths and mean quadratic paths is straightforward. On the other hand, their disadvantages are that the relative variance γ² cannot be varied. This is why rectangular enclosures are now considered.

3.2 2D rectangular enclosures

The cases of 2D rectangular enclosures require more complex integration. Basically, a two-step integration procedure is used for 2D rectangular enclosures, first on directions for rays starting on a running point on one side, taking into account a Lambert distribution of direction (cosine law), then on the running point along the side. The procedure must be repeated on the four sides of the rectangle. This procedure recovers the expected mean free paths $\overline{\ell }=\frac{\pi \mathit{Lh}}{2(L+h)}$, where L and h are respectively the width and the height of the enclosure. The mean quadratic path takes a much more complex expression – see Supplementary material:

$$\begin{array}{cc}\hfill \overline{{\ell }^2}=& \frac{2}{3(L+h)}\left[{\left(L^2+h^2\right)}^{\frac{3}{2}}-\left(L^3+h^3\right)\right]\hfill \\ \hfill & +\frac{h^2}{L+h}\{Lln\left(\frac{L-h+\sqrt{L^2+h^2}}{L+h-\sqrt{L^2+h^2}}\right)\hfill \\ \hfill & -\frac{L^2}{\sqrt{L^2+h^2}-h}+2h\}\hfill \\ \hfill & +\frac{L^2}{L+h}\{hln\left(\frac{h-L+\sqrt{L^2+h^2}}{h+L-\sqrt{L^2+h^2}}\right)\hfill \\ \hfill & -\frac{h^2}{\sqrt{L^2+h^2}-L}+2L\}.\hfill \end{array}$$

The full computation of $\overline{{\ell }^2}$ diverges as h²ln(2L/h) for large values of the ratio L/h; in other words, the very long paths dominate for elongated enclosures.

Tables 4 and 5 present the same reverberation times as in Table 1, respectively for a square (γ²=0.205), and a rectangle of ratio 1:10 (γ²=0.573). Note that much larger γ² values are obtained than in the case of the disc and the sphere. In fact, any ratio can be achieved, provided the rectangle is elongated enough. Note also that Kuttruff's formula diverges for R<5.8×10⁻⁵ in the case of the square, and α>0.969 in the case of the rectangle of ratio 1:10.

Figure 3 presents some of the reverberation times in graphic form. Once more, all reverberation times take values between the Sabine and the Eyring formulas, which can be considered as bounds for reverberation times. Kuttruff's formula is not plotted for the square, as it is essentially not distinguishable from equation (8) with u=γ²=0.205 and m=0. It is however plotted in the rectangle of ratio 1:10 since it starts diverging for α>0.6.

Figure 3 Reverberation times for uniform absorption coefficient in square (top) and rectangle of dimensions 1:10 (bottom).

In these cases of larger relative variances, the differences between the last four reverberation formulae are much larger. They reach more than 10% for α=0.9.

3.3 3D rectangular enclosures

For 3D parallelepipeds, Kuttruff [6, 20] has carried out Monte-Carlo simulations for some special cases of relative enclosure dimensions. The numerical mean free paths are very close to the theoretical values; and γ² varies between 0.344 for a cube to 0.646 for a flat enclosure of ratio 1:10:10, which are roughly twice as large values as for the corresponding 2D rectangular enclosures. Tables 6–8 present the same reverberation times as in Table 1, respectively for the cube, ratio 1:1:10 (γ²=0.421), and ratio 1:10:10. In these cases, Kuttruff's formula diverges for α larger than 0.997, 0.991, and 0.955 for respectively the cube and the enclosure of dimensions 1:1:10 and 1:10:10.

Added to Tables 6–8 are the reverberation times computed with CATT-Acoustic V.9 [25], using the interactive reverberation time algorithm, because the process closely follows the theoretical assumptions of following one ray through many reflections. Diffusion was set at the maximum of 100%, so that reflections follow Lambert law as required for ergodic enclosures – see Section 2.1. Simulations are in good agreement with the different improved formulas; but computation becomes more and more unstable when absorption increases, and finally breaks down for very large absorptions (α>0.8). In the case of the flat enclosure, it even starts diverging at α=0.8 and above.

3.3.1 Path distribution

Kuttruff [6, 20] has extensively studied the case of ergodic rectangular parallelepipeds and evaluated the relative variances of their free flight paths by means of Monte-Carlo simulations. It is however possible to analytically compute the mean free flight path and the mean quadratic free flight paths of a parallelepiped of dimensions L·l·h through a four-step integration procedure: first on directions for rays starting on a running point on rectangular slice of the parallelepiped, taking into account a Lambert distribution of direction (cosine law); then on the running point along the rectangular slice; then along the direction perpendicular to the slice; and finally rotating the result around azimuths. All results are analytical, except for one integral in the last step. This procedure does recover the expected mean free paths $\overline{\ell }=\frac{4\mathit{Llh}}{2(\mathit{Ll}+\mathit{lh}+\mathit{hL})}$, where L, l, and h are respectively the length, the width, and the height of the enclosure. The mean quadratic path takes a lengthy expression – see Supplementary material:

$$\begin{array}{ccc}\hfill \overline{{\ell }^2}& \hfill =\hfill & \frac{1}{\pi S}\{2L^2{l}^{2}ln\frac{\sqrt{L^2+l^2}}{\sqrt{L^2+l^2+h^2}}\hfill \\ \hfill & \hfill \hfill & +\frac{L^4}{3}ln\frac{L\sqrt{L^2+l^2+h^2}}{\sqrt{L^2+l^2}\sqrt{L^2+h^2}}\hfill \\ \hfill & \hfill \hfill & +\frac{4}{3}{L}^{3}larctan\frac{l}{L}+\frac{4}{3}{\mathit{Ll}}^{3}arctan\frac{L}{l}\hfill \\ \hfill & \hfill \hfill & -\frac{4}{3}(h^2+L^2{)}^{\frac{3}{2}}larctan\frac{l}{\sqrt{L^2+h^2}}\hfill \\ \hfill & \hfill \hfill & -4h^2\mathit{Ll}[ln\frac{h}{\sqrt{L^2+h^2}+L}arctan\frac{l}{L}\hfill \\ \hfill & \hfill \hfill & +ln\frac{h}{\sqrt{l^2+h^2}+l}arctan\frac{L}{l}]\hfill \\ \hfill & \hfill \hfill & -2h^2\mathit{Ll}{\displaystyle \sum _{n=1}^{+\infty }}(\left[1-{\left(\frac{\sqrt{l^2+h^2}-l}{\sqrt{l^2+h^2}+l}\right)}^n\right]\hfill \\ \hfill & \hfill \hfill & -(-1{)}^n\left[1-{\left(\frac{\sqrt{L^2+h^2}-L}{\sqrt{L^2+h^2}+L}\right)}^n\right])\hfill \\ \hfill & \hfill \hfill & \times \frac{sin\left(2narctan\frac{l}{L}\right)}{n^2}\}+\mathit{cp}\hfill \end{array}$$

where S=2(Ll+lh+hL) and cp consists in all the circular permutations of L, l and h in the expression that precedes it. Combined with the mean free path ${\ell }_{\mathit{TH}}=\overline{\ell }=\frac{4V}{S}$, where V=Llh, the theoretical relative variance ${\gamma }_{\mathit{TH}}^2$ is easily computed. It is presented in Table 9 for a selection of relative dimensions together with the results from Kuttruff's simulations [6], labelled with index MC.

Figure 4 presents the path distributions as published by Kuttruff [6, 20] for three different parallelepipeds, and compares them with the Gamma distribution (Eq. (2)) when ${\gamma }^2={\gamma }_{\mathit{TH}}^2$ taken from Table 9. It can be seen on Figure 4 that the true distributions present jumps when the free paths become larger than some characteristic values, which correspond well with the maxima of the theoretical distributions. Symmetry considerations show that the jumps happen when the free paths are equal to the smallest dimensions l_min of the parallelepipeds. Therefore Figure 4 also present the theoretical distribution for ${\gamma }^2=\frac{{\ell }_{\mathit{TH}}}{l_{min}+{\ell }_{\mathit{TH}}}$ (broken line), since the maxima of the distributions occur for $\frac{\ell }{{\ell }_{\mathit{TH}}}=\frac{1}{{\gamma }^2}-1$ as the free paths ℓ are normalized by the mean free paths $\overline{\ell }={\ell }_{\mathit{TH}}$ in Figure 4. The maxima of these modified distributions are no longer aligned with the jumps of the true distributions, but the agreement improves for very disproportionate parallelepipeds where one dimensions is much smaller than the other two.

Figure 4 Free path distributions for a few shapes of parallelepipeds, compared with arbitrarily scaled theoretical distributions (after [20]). Continuous lines: ${\gamma }^2={\gamma }_{\mathit{TH}}^2$; broken lines: ${\gamma }^2=\frac{{\ell }_{\mathit{TH}}}{l_{min}+{\ell }_{\mathit{TH}}}$.

3.3.2 Uniform boundary absorption

Figure 5 presents reverberation times calculated for three different shapes of parallelepipeds as function of uniform absorption coefficient on all boundaries. Remember that curves “uTH” correspond to equation (8) with m=0 and only one value of $u={\gamma }_{\mathit{TH}}^2$, with $N(u)=\frac{1}{\overline{t}{\gamma }_{\mathit{TH}}^2}$; and curves “uMFP” to equation (8) with m=0 and only one value of $u=\frac{{\ell }_{\mathit{TH}}}{l_{min}+{\ell }_{\mathit{TH}}}$, with $N(u)=\frac{l_{min}+{\ell }_{\mathit{TH}}}{\overline{t}{\ell }_{\mathit{TH}}}$. Despite values of u that differ with up to 16% according to Table 9, the corresponding reverberation times only marginally differ over the whole range of absorptions. Note however that reverberation times according to Kuttruff's formula diverge for absorption coefficients larger than 0.8, as expected from equation (9) with γ²=0.646.

Figure 5 Reverberation times for uniform absorption coefficient in three shapes of parallelepipeds. Top: dimensions 1:10:10; middle: dimensions 1:10:10; bottom: dimensions 1:1:1.

Kuttruff's modified reverberation time (Eq. (11)) is not plotted on Figure 5 since it follows curves “uTH” and “uMFP” except for full absorption where its value is 0, just like Eyring's reverberation time. Note that all the new reverberation times lay between Sabine's and Eyring's, as noticed before and explained by Joyce [15].

4 Discussion

4.1 Influence of relative variance

In the previous examples, care was taken that the mean free path and the mean waiting time were kept to the same values. As a consequence, Sabine's and Eyring's reverberation times take the same values for all enclosures. The other reverberation times vary, as they depend on the relative variance γ². However, the dispersions are less than 5% within one same enclosure for the usual range of absorptions (α≤0.5), and increase with γ² for large absorption coefficients, where they can reach more that 20%.

In general, computing the relative variance is not an easy task, as exemplified in Sections 3.2 and 3.3. However, Section 3.3.2 offers an easy way to approximate it, and shows that even in the cases where the approximation is 16% off as in cubes, reverberation times computed with this approximation are comparable to those obtained with the theoretical values. As a consequence, the corresponding reverberation times appear as last lines in Tables 1–5, as last line but on in Tables 6–8, but not for the disc (Tab. 2) and the sphere (Tab. 3) which do not present a sudden increase in the probability distributions of free path above some threshold. They also appear as crosses in Figures 1 and 3. In all cases but for the square, for which the approximation is more than 100% off, reverberation times thus computed are comparable with those computed with the theoretical relative variance.

This approximation makes it also possible to compute asymptotic values for the relative variance γ², since mean free paths remain finite even for very elongated enclosures. This asymptotic value is $\frac{\pi }{\pi +2}=0.611$ for rectangular enclosures, $\frac{2}{3}=0.667$ for flat enclosures, and $\frac{1}{2}=0.5$ for long enclosures. Nevertheless, the exact relative variances diverge as lnL for very disproportionate enclosures, where L is the largest dimension of the enclosure, as is shown for 2D rectangles in Section 3.2; but it probably has no practical value, as explained below in Section 4.3.

Note that Kuttruff's reverberation time diverges for very small values of the reflection coefficient R, that is to say large absorptions. But for practical applications, this only happens for mean absorption coefficient larger than 0.9. Only disproportionate 2D enclosures, with γ²>0.86, reach this divergence for α<0.9; for a rectangular 2D enclosure, it corresponds to a dimension ratio larger than 45:1, which indeed has no practical use.

4.2 Validity for large absorption

The different reverberation formulas presented in this paper are only valid for ergodic enclosures, as they are all based on the asymptotic distribution of reflections along one single ray. Indeed, only one ray needs be considered in ergodic enclosures, as any randomly chosen starting position evolves in a ray that eventually comes infinitesimally near any other combination of position and direction. As the trajectories of the rays do not depend on the value of the reflection coefficient, the distribution of reflections is independent of absorption. As a consequence, the reverberation formulas are valid for all absorption values.

It should be noted, however, that the formulas make use of the asymptotic distribution of reflections. In the case of almost complete absorption on the boundaries, the energy decays very quickly along any ray and may become almost equal to zero before the asymptotic distribution is reached. In other words, reverberation formulas, although still valid in principle, cannot be applied to very large absorptions. This seriously reduces the importance of the divergence of Kuttruff's formula discussed in the previous section.

4.3 Non-uniform boundary absorption

The case of non-uniform absorption derives from the previous discussion: as long as the asymptotic distribution of reflections is attained before energy along the ray has decayed almost to zero, the different formulas studied in this paper are valid. One must only attribute a certain probability to each absorption coefficient, and when absorption is not angle dependent, the obvious choice is the proportion of boundary surface corresponding to each absorption coefficient. With this choice, a mean absorption coefficient is readily obtained, with which the different formulas can be computed.

More precisely, in case of different absorption coefficients on different boundaries, trajectories will undergo several reflections on the other boundaries before hitting the same boundary again. The general form of equation (8) holds, with m=0 and $N_j(u)=\frac{1}{\overline{t_j}{\gamma }_j^2}$, that is:

$$\text{RT}=\frac{13.8\overline{t_j}{\gamma }_i^2}{{\sum }_j\left(1-R_j^{{\gamma }_j^2}\right)}$$(17)

where $\overline{t_j}$ and ${\gamma }_j^2$ may take different values on different boundaries. However, consistency with the uniform case leads to the same value of ${\gamma }_j^2$ on all boundaries, ${\gamma }_j^2={\gamma }_{\mathit{TH}}^2$, that is, depending only on the shape of the enclosure; and to $\overline{t_j}$ inversely proportional to S_j, the area of boundary j, so that ${\sum }_j\frac{1}{\overline{t_j}}=\frac{1}{\overline{t}}=\frac{c}{{\ell }_{\mathit{TH}}}$ where ${\ell }_{\mathit{TH}}=\frac{4V}{S}$ is the mean free path. As a consequence, the reverberation time is now given by:

$$\text{RT}=\frac{13.8·4V{\gamma }_{\mathit{TH}}^2}{{\sum }_j{\mathit{cS}}_j\left(1-R_j^{{\gamma }_{\mathit{TH}}^2}\right)}·$$(18)

Now, if absorption is angle dependent, the quantity $R_j^{{\gamma }_{\mathit{TH}}^2}$ must be replaced in equation (18) by its average $\overline{R_j^{{\gamma }_{\mathit{TH}}^2}}$ over all angles of incidence. This comes from the fact that, in ergodic enclosures, rays reach almost any position with almost any incidence. Thus, in the long run, all incidences happen on any part of the boundary. In other words, the equivalent absorption coefficient of equation (4) becomes:

$${\alpha }^{*}=\frac{1}{S}{\displaystyle \sum _j}\frac{S_j}{{\gamma }_{\mathit{TH}}^2}\left(1-\overline{R_j^{{\gamma }_{\mathit{TH}}^2}}\right).$$(19)

Note that equation (18) tends toward Millington-Sette's formula when ${\gamma }_{\mathit{TH}}^2⟶0$.

Also note that in disproportionate enclosures, sound energy decreases along the long dimensions [26, 27, 28]. As a consequence, less energy reaches the far end surfaces, reducing their efficiency in absorbing sound. As a consequence, equation (18) most probably does not represent reality. Reference [26] shows that Eyring's formula probably is more appropriate for elongated enclosures.

4.4 E-type vs. S-type reverberation times

From equation (18) and Section 2.2, it is obvious that the most general reverberation formula is:

$$\text{RT}=\frac{13.8}{{\sum }_j\left[-m_{j}ln{R}_j+\frac{1}{{\gamma }_j^2\overline{t_j}}\left(1-R_j^{{\gamma }_j^2}\right)\right]}$$(20)

where the ${\gamma }_j^2$ are no longer automatically equal to ${\gamma }_{\mathit{TH}}^2$, or equivalently for the corresponding cumulant function:

$$\psi \left(z,t\right)=t{\displaystyle \sum _j}\left[m_j\mathit{iz}+{\displaystyle {\int }_0^{+\infty }}\left(e^{{\mathit{izu}}_j}-1\right)\text{d}{N}_j(u_j)\right]$$(21)

since the u_j must all be positive. According to Lévy [16], terms in m_j correspond to reflections at regular time intervals, for example in 1D enclosures, which are not likely to take place in 3D enclosures; and terms in N_j(u_j) correspond to generalized Poisson processes, with expression:

$$P(k|t)=\frac{1}{\mathrm{\Gamma }\left(\frac{1}{{\gamma }_j^2}\right)}{\left(\frac{t}{{\gamma }_j^2\overline{t}}\right)}^{\frac{k}{{\gamma }_j^2}-1}exp\left(-\frac{t}{{\gamma }_j^2\overline{t}}\right)$$(22)

where $\frac{1}{{\gamma }_j^2\overline{t}}=\text{d}{N}_j(u_j)$. Now, it is obvious from Table 9 that $\frac{k}{{\gamma }^2}$ hardly ever takes an integer value. Therefore, equation (22) must be understood as an asymptotic expression of the distribution of the number of reflections. It can however be interpolated to all values of waiting times t, just as the Eyring formula interpolates the discrete jumps of energy that occur at regular time intervals to a continuous function.

One can therefore distinguish the reverberation time formulas that behave like the Eyring formula and give a vanishing reverberation time as soon as one surface is totally absorbing: they deliver “E-type” reverberation times. They usually contain terms in lnR_j, that is, all m_j are not null. But Kuttruff's original formula equation (10) also reach zero for total reflection and gives an E-type reverberation time. Formulas with all m_j=0 and a discrete number of u_j, on the other hand, give reverberation times that remain strictly positive even for total absorption on the boundaries. They behave like the Sabine formula, obtained with equation (8) when m=0 and u=1: therefore, they deliver “S-type” reverberation times.

It should be noted that the two types of reverberation times correspond to different aspects of the reverberation process. The case of one totally absorbing surface makes it clear, with the help of equation (20), that E-type reverberation times correspond to the late part of the decay, when encounter on the totally absorbing surface instantaneously annihilates the sound energy along a trajectory. S-type reverberation times, on the other hand, correspond to the initial part of the decay: they describe the probability of survival of trajectories as time elapses. As a consequence, the debate on the most accurate reverberation time formula, which focussed the attention of researchers in the first half of last century, looses most of its relevance; but S-type reverberation times are probably more appropriate for practical applications.

5 Non-ergodic enclosures

First of all, it should be clear to the reader that non ergodic enclosures hardly exist. Unevenness on the boundaries creates diffusion, as is experimented in real rooms. Therefore, the cases of non-ergodic spheres and non-ergodic rectangular parallelepipeds remain scholarly exercises.

In non-ergodic spheres and rectangular parallelepipeds, that is, with pure specular reflections at the boundaries, reverberation time depends on the direction of the path as it fixes the time intervals between successive reflections. Besides, in a sphere, all successive free paths have the same length, defined by the initial directions of the rays. As a consequence, such enclosures are equivalent to infinitely many coupled rooms with vanishing coupling coefficients. As proved in ([29], Appendix D), the overall reverberation time steadily increases along the decay, depending on the relative energy in each direction, and stabilizes in the long run at its largest value, usually the reverberation time corresponding to the largest dimension of the enclosure. In real life, this stabilization occurs rather quickly since diffusion occurs at the edges of the enclosure, simply due to the fact that reflection intensity is discontinuous at edges. If this is not the case, the enclosure probably is composed of coupled spaces.

5.1 Statistics of reflections

As mentioned above, reflections of any ray in spheres occur at regular intervals fixed by the initial direction of the ray. As for rectangular parallelepiped enclosures, for any ray direction, the reflections on each pair of parallel walls also occur at regular intervals. In other words, in both cases for each direction, the total number n of reflections during time t is proportional to the distance traveled ct. It is then possible to introduce the reduced variable x=n/ct. After a sufficiently long time, the distribution of this reduced variable becomes independent of instant t, that is to say that the choice of the observation instant fixes the maximum value of the error thus committed.

On the other hand, the analysis of the number of reflections in each direction shows that the reduced variable x remains above a strictly positive value, X, which is none other than the inverse of the largest dimension of the enclosure; it is therefore effectively reached by x. Note that the upper bound for x is usually +∞ as paths are vanishingly short along the edges of a rectangular enclosure, or even in a sphere when reflecting points become very close.

The temporal independence of the probability distribution of the reduced variable is one of the particularities of spherical and rectangular parallelepiped enclosures with specular reflections at the boundaries. In the ergodic random enclosures studied in Section 2.2, the number of reflections satisfies a generalized Poisson distribution, incompatible with the temporal independence of the distribution of the reduced variable. More generally, the existence of a time-independent distribution for the reduced variable is incompatible with a Lévy process for which not only the mean, but also all terms in the cumulant function Ψ(z, t), are proportional to the path length ct (see Eqs. (6) and (7)).

5.2 Impulse response

The introduction of the distribution p(x) of the reduced variable x=n/ct makes it possible to evaluate the moments of the impulse response constituted by the superposition of arrivals, that is, rays that have traveled through the enclosure and reach the receiver at some time. We consider that the source emits a very short positive crenel of duration T and height $1/\sqrt{\mathit{cT}}$ immediately followed by a similar but negative crenel. In this case, the average value of the impulse response is always equal to the sum of the average values of each arrival; it is therefore zero for a very brief signal. As for the mean square value, it is proportional to the mean value of the quantity Rⁿ, where R still is the energy-based reflection coefficient. We then obtain:

$$\langle h^2(t)\rangle =\frac{\mathit{ct}}{2\pi V}{\displaystyle {\int }_X^{+\infty }}{R}^{n}p(x)\text{d}x.$$(23)

By setting z=−ctlnR, this integral takes the form of the characteristic function of the distribution p(x), with the difference that z is a pure imaginary variable for a characteristic function, while z is real here. By analogy with a characteristic function, however, we pose:

$${\displaystyle {\int }_X^{+\infty }}{e}^{-\mathit{zx}}p(x)\text{d}x=\mathrm{\Phi }(z)=e^{+\mathrm{\Psi }(z)}.$$(24)

There follows from equation (23) that the shape of the reverberation curve for spherical or rectangular parallelepiped enclosures will be deduced from the properties of the characteristic function Φ(z), or more precisely, from the cumulant function Ψ(z)=lnΦ(z).

5.3 Properties of cumulant function Ψ(z)

• i)

  Ψ(z) is a decreasing function of z. This is proved by computing:

  $${\mathrm{\Psi }}^{\prime }(z)=\frac{\mathrm{\Phi }\prime (z)}{\mathrm{\Phi }(z)}=\frac{{\int }_X^{+\infty }(-x)e^{-\mathit{zx}}p(x)\text{d}x}{{\int }_X^B{e}^{-\mathit{zx}}p(x)\text{d}x}$$(25)

  and introducing the distribution:

  $$p_z(x)=\frac{e^{-\mathit{zx}}p(x)}{{\int }_X^{+\infty }{e}^{-\mathit{zx}}p(x)\text{d}x}·$$(26)

  It is a positive distribution, of sum equal to unity. It is therefore a probability distribution, and [−Ψ′(z)] is none other than the average value $\overline{x_z}$ of the distribution p_z(x). As x is bounded below by X, Ψ′(z) is also bounded above by −X; Ψ′(z) is therefore negative. We deduce from this that Ψ(z) is indeed a decreasing function of z.

• ii)

  Ψ′(z) is an increasing function of z. According to expressions (25) and (26):

  $$\begin{array}{cc}\hfill {\mathrm{\Psi }}^{\prime \prime }(z)& ={\displaystyle {\int }_X^{+\infty }}(-x{)}^2{e}^{-\mathit{zx}}p(x)\text{d}x\hfill \\ \hfill & -{\left\{{\displaystyle {\int }_X^{+\infty }}(-x)e^{-\mathit{zx}}p(x)\text{d}x\right\}}^2\hfill \\ \hfill & =\overline{x_z^2}-{\overline{x}}^2={\sigma }_z^2\hfill \end{array}$$

  ${\mathrm{\Psi }}^{\prime \prime }(z)$ is none other than the variance of the distribution p_z(x), a quantity always positive. Therefore Ψ′(z) is an increasing function of z.

• iii)

  lim_z→∞p_z(x)=δ(x−X). Expression (26) can be written for every x as:

  $$p_z(x)=\frac{e^{-z(x-X)}p(x)}{{\int }_X^{+\infty }{e}^{-z(x-X)}p(x)\text{d}x}·$$

Now, for any x distinct from X, the quantity e^−z(x−X)p(x) can be made as small as desired by choosing z sufficiently large; while for x=X, it remains equal to p(X) for any value of z. On the other hand, the distribution p(x) is right continuous at point X (this is a property of spherical and rectangular parallelepiped enclosures). So the integral:

$${\displaystyle {\int }_X^{+\infty }}{e}^{-z(x-X)}p(x)\text{d}x$$(27)

remains of the order of p(X)/z when z becomes very large, because p(X) is finite and positive – note that this is not the case for circular enclosures, hence their exclusion. As this quantity decreases with z less quickly than e^−z(x−X)p(x) for x>X, we deduce that, for x>X, the probability p_z(x) can be made as small as desired by choosing z sufficiently large.

Now, by definition, for any value of z:

$${\displaystyle {\int }_X^{+\infty }}{p}_z(x)\text{d}x=1.$$

This allows us to affirm:

$$\underset{z\to \infty }{lim}{p}_z(x)=\delta (x-X).$$

As a result, we also obtain:

$$\underset{z\to \infty }{lim}{\mathrm{\Psi }}^{\prime }(z)=-X\hspace{0.5em}\hspace{0.5em},\hspace{0.5em}\hspace{0.5em}\underset{z\to \infty }{lim}{\mathrm{\Psi }}^{\prime \prime }(z)=-0.$$

Note that, in the case of angle dependent absorption, R is a function of x because angles of incidence remain constant along any ray. However, the same analysis is valid when setting z=ct and replacing x by −xlnR(x), so that X now is the minimum of −xlnR(x).

5.4 Reverberation of rectangular enclosures

By reporting the properties of the characteristic function in equation (23) of energy in the enclosure, we deduce that after a sufficiently long time, the energy decreases with the logarithmic decrement Xc|lnR|. As the value X depends only on the largest dimension of the enclosure, it is this dimension which imposes the decay of the late sound.

However, we cannot define a reverberation time, because the shape of the energy equation after a very long time depends on the complete estimation of the characteristic function (Eq. (24)). Now this characteristic function can be estimated using expression (27). We deduce from this that, for z sufficiently large:

$$\mathrm{\Phi }(z)={\displaystyle {\int }_X^{+\infty }}{e}^{-\mathit{zx}}p(x)\text{d}x\approx e^{-\mathit{zX}}\frac{p(X)}{z}·$$

In terms of reverberation, energy equation (23) becomes:

$$\begin{array}{cc}\hfill \langle h^2(t)\rangle & =\frac{\mathit{ct}}{2\pi V}\frac{p(X)}{c|lnR|}\frac{e^{-(\mathit{Xc}|lnR|)t}}{t}\hfill \\ \hfill & =\frac{1}{2\pi V}\frac{p(X)}{|lnR|}{e}^{-(\mathit{Xc}|lnR|)t}·\hfill \end{array}$$

Energy therefore decreases exponentially, but only after an infinite time, which makes it difficult to define a reverberation time.

In case of angle dependent reflection coefficients, equation (23) simply becomes:

$$\langle h^2(t)\rangle =\frac{1}{2\pi V}p(X)e^{-\mathit{Xct}}$$(28)

since $X=min\left\{x|lnR(x)|\right\}$.

6 Conclusion

This paper has presented the most general expression for reverberation in ergodic enclosures, valid for all values of absorption, even almost total absorption. It depends on the mean free path, but also on the relative variance of the free paths, which is a function of the shape of the enclosure. For rectangular parallelepipeds, Kuttruff's conjecture that the relative variance remains smaller than one is always verified in practical cases, even though it becomes infinite for extremely large flat enclosures.

Assuming that the asymptotic distribution of reflections along one single ray controls reverberation, which is the case for ergodic enclosures [14, 17], the generic decay rate for such enclosures is presented in Section 2.2. It is then simplified according to different assumptions, leading to the Sabine, the Eyring, and Kuttruff's formulas, and the corresponding reverberation times are presented for a few enclosures with simple shape, both in 2D and 3D cases. Good news are the small variations between shape and dimensions, as the piloting factors are the mean waiting time and its relative variance. Thus, reverberation does not depend on dimension.

Note that Kuttruff basically developed the same ideas in the 1970s. However, he lacked a theory to justify the different probability distributions he studied at that time. Compared to his achievements, the major improvement brought by the present study precisely is a general theory of reverberation, based on rigorous mathematical and physical principles.

This paper also proposes a simple mean to approximate the relative variance of a rectangular enclosure, using the ratio of it smallest dimension to the mean free path (see Sect. 3.3.2). This is a relatively large approximation, up to 16% for cubic enclosures, and even is wrong for extremely large flat enclosures; but Figure 5 shows that it nevertheless allows an accurate approximation of the reverberation time up to total absorption on the boundaries. This approximation can probably be extrapolated to non-rectangular enclosures, as shown in the case of the stadium.

The cases of non uniform and angle dependent absorptions were also investigated, leading to an expressions of the equivalent absorption coefficient (Eq. (19)) that generalizes known expressions to angle dependency.

The consequences of the general expression are twofold. Firstly, computer algorithms which apply extrapolations of reverberation tails from free path statistics need be modified and evaluate mean quadratic free paths beside the mean free paths. Secondly, there is a need, as already noted by Schroeder 60 years ago [21], for experimentally investigating whether absorption measurements in reverberation chambers can be improved by assessing the relative variance of the free paths, and the influence of diffusers on this relative variance, as Kuttruff started to investigate [30]. The same holds for auditoriums, where further investigations are needed to check the relevance of the present theory. Also note that reverberation times in non-ergodic enclosures increase with time and eventually converge to a finite value, driven by the longest free paths with the smallest absorption at their extremities.

In conclusion, the present study proves once more that the Sabine and the Eyring reverberation formulas are first order approximations that do not take into account the shape of the enclosure. It is however possible to improve these formulas, taking into account the relative variance of the free paths, or equivalently, of the waiting time between two successive reflections, leading to improved equivalent absorption coefficients. All improved formulas are compatible with the generic model of a Lévy process for reflection distribution along any ray; and they deliver very similar results as long as the enclosure is not too disproportionate, and as long as some diffusion occurs on its boundaries. It is therefore not possible to recommend one specific formula, either for reverberation chambers or for auditoriums, but only that both the Sabine and the Eyring formulas should be applied with caution, and with the awareness that they are approximations only.

Supplementary material

Supplement to Revising reverberation Access here

Acknowledgments

The author wish to thank the anonymous reviewers who signaled reference [21] that had escaped his attention, and whose challenging questions lead to serious improvements, including angle dependency.

Conflicts of interest

The author declares no conflict of interest.

Data availability statement

The research data associated with this article are included in the supplementary material of this article.

Appendix A Kuttruff's revised formula

A.1 General expression for level decay

According to Section 2.2, reverberation is described by the probability distribution of an infinitely divisible probability law. Lévy [16] has shown that such a law is defined by the canonical cumulant function:

$$\begin{array}{cc}\hfill \psi \left(z,t\right)& =\mathit{mtiz}-g^{2}t\frac{z^2}{2}\hfill \\ \hfill & +t{\displaystyle {\int }_{-\infty }^{+\infty }}\left(e^{\mathit{izu}}-1-\frac{\mathit{izu}}{1+u^2}\right)n(u)\text{d}u\hfill \end{array}$$

where N(u)=∫n(u)du is a non-decreasing function of u with the following properties:

• •

  it takes finite values for u=±∞ (can be null),

• •

  it can be discontinuous in u=0,

• •

  ∫u²n(u)du is finite on all finite intervals.

In that case, sound level decreases as:

$$\begin{array}{cc}\hfill L(R,t)& =\psi \left(-iln(\text{R}),t\right)\hfill \\ \hfill & =\mathit{mt}lnR)+g^{2}t\frac{{ln}^{2}R}{2}\hfill \\ \hfill & +t{\displaystyle {\int }_{-\infty }^{+\infty }}\left(R^u-1-\frac{ulnR}{1+u^2}\right)n(u)\text{d}u\hfill \end{array}$$

where R is the mean reflection coefficient on the boundaries.

A.2 Modified Kuttruff's level decay

Kuttruff has considered the probability distribution p(n|t) of reflections along a ray during the time span t. For this, he starts with the probability distribution p(t|n) of travelling times t spanning n reflections. The central theorem of probability leads him to postulate a Gaussian probability distribution for p(t|n), namely:

$$\begin{array}{cc}\hfill p(t|n)=\frac{1}{\sqrt{2\pi n}\gamma \overline{t}}exp\left(-\frac{(t-n\overline{t}{)}^2}{2n{\gamma }^2{\overline{t}}^2}\right)& \hfill \end{array}$$

where $\overline{t}=\frac{\overline{\ell }}{c}$ is the mean free waiting time between successive reflections – c being the speed of sound – and γ² the corresponding reduced variance. This distribution is also proportional to p(n|t) by modification of the factor:

$$\begin{array}{cc}\hfill p(n|t)=\frac{1}{\sqrt{2\pi n}\gamma }exp\left(-\frac{(t-n\overline{t}{)}^2}{2n{\gamma }^2{\overline{t}}^2}\right)& \hfill \end{array}$$

and with cumulant function ψ(z, t) equal to the logarithm of the characteristic function:

$$\begin{array}{cc}\hfill \varphi (z,t)& ={\displaystyle {\int }_0^{+\infty }}p(n|t)exp(\mathit{inz})\text{d}n\hfill \\ \hfill & ={\displaystyle {\int }_0^{+\infty }}\frac{exp\left(-\frac{(t-n\overline{t}{)}^2}{2n{\gamma }^2{\overline{t}}^2}\right)}{\sqrt{2\pi n}\gamma }exp(\mathit{inz})\text{d}n\hfill \end{array}$$

where the discrete number of reflections has been approximated by a continuous variable. The integration is carried out by rewriting the exponential term as:

$$\begin{array}{cc}\hfill & exp\left(-\frac{(t-n\overline{t}{)}^2}{2n{\gamma }^2{\overline{t}}^2}\right)exp(\mathit{inz})\hfill \\ \hfill & =exp\left(-\frac{(t-n\overline{t}{)}^2}{2n{\gamma }^2{\overline{t}}^2}+\mathit{inz}\right)\hfill \\ \hfill & =exp\left(-\frac{t^2-2\mathit{tn}\overline{t}+n^2{\overline{t}}^2-2{\mathit{in}}^2{\gamma }^2{\overline{t}}^{2}z}{2n{\gamma }^2{\overline{t}}^2}\right)\hfill \\ \hfill & =exp\left(-\frac{t^2-2\mathit{tn}\overline{t}+n^2{\overline{t}}^2(1-2i{\gamma }^{2}z)}{2n{\gamma }^2{\overline{t}}^2}\right)\hfill \\ \hfill & =exp\left(+\frac{t}{{\gamma }^2\overline{t}}\right)exp\left(-\frac{t^2}{2n{\gamma }^2{\overline{t}}^2}-\frac{n(1-2i{\gamma }^{2}z)}{2{\gamma }^2}\right)·\hfill \end{array}$$

Thus:

$$\varphi (z,t)=\frac{exp\left(+\frac{t}{{\gamma }^2\overline{t}}\right)}{\sqrt{2\pi }\gamma }{\displaystyle {\int }_0^{+\infty }}{n}^{\nu -1}exp\left(-\beta n^p-\delta n^{-p}\right)\text{d}n$$

with $\beta =\frac{1-2i{\gamma }^{2}z}{2{\gamma }^2}$, $\delta =\frac{t^2}{2{\gamma }^2{\overline{t}}^2}$, p=1, and $\nu =\frac{1}{2}$. The integral is then equal to ([31], Eq. (3.478.4), p. 342):

$$\begin{array}{cc}\hfill & \frac{2}{p}{\left(\frac{\delta }{\beta }\right)}^{\frac{\nu }{2p}}{K}_{\frac{\nu }{p}}\left(2\sqrt{\beta \delta }\right)\hfill \\ \hfill & =2{\left(\frac{t^2}{{\overline{t}}^2(1-2i{\gamma }^{2}z)}\right)}^{\frac{1}{4}}{K}_{\frac{1}{2}}\left(\frac{t\sqrt{1-2i{\gamma }^{2}z}}{{\gamma }^2\overline{t}}\right)\hfill \\ \hfill & =2{\left(\frac{t}{\overline{t}\sqrt{1-2i{\gamma }^{2}z}}\right)}^{\frac{1}{2}}\sqrt{\frac{\pi }{2}\frac{{\gamma }^2\overline{t}}{t\sqrt{1-2i{\gamma }^{2}z}}}\hfill \\ \hfill & \times exp\left(-\frac{t\sqrt{1-2i{\gamma }^{2}z}}{{\gamma }^2\overline{t}}\right)\hfill \\ \hfill & =\frac{\sqrt{2\pi }\gamma }{\sqrt{1-2i{\gamma }^{2}z}}exp\left(-\frac{t\sqrt{1-2i{\gamma }^{2}z}}{{\gamma }^2\overline{t}}\right)\hfill \end{array}$$

since $K_{\frac{1}{2}}\left(z\right)=\sqrt{\frac{\pi }{2z}}exp(-z)$. Thus:

$$\varphi (z,t)=\frac{exp\left(-\frac{t\left[\sqrt{1-2i{\gamma }^{2}z}-1\right]}{{\gamma }^2\overline{t}}\right)}{\sqrt{1-2i{\gamma }^{2}z}}·$$

In other words, the cumulant function of the modified Kuttruff distribution is:

$$\psi \left(z,t\right)=-t\frac{\sqrt{1-2{\gamma }^2\mathit{iz}}-1}{\overline{t}{\gamma }^2}-ln\sqrt{1-2i{\gamma }^{2}z}$$

and the level decay function is expressed by (see Sect. A.1):

$$\begin{array}{cc}\hfill L(R,t)& =\psi \left(-ilnR,t\right)\hfill \\ \hfill & =-t\frac{\sqrt{1-2{\gamma }^{2}lnR}-1}{\overline{t}{\gamma }^2}-\frac{1}{2}ln\left(1-2{\gamma }^{2}lnR\right).\hfill \end{array}$$

A.3 Fourier decomposition of the cumulant function

In order to retrieve the canonical form of the cumulant function, one needs to carry out the inverse Fourier transform of Kuttruff's modified cumulant function:

$$\begin{array}{cc}\hfill \mathrm{\Psi }\left(u,t\right)& ={\mathrm{\Psi }}_1\left(u,t\right)+{\mathrm{\Psi }}_2\left(u,t\right)\hfill \\ \hfill & =-\frac{t}{\overline{t}{\gamma }^2}{\displaystyle {\int }_{-\infty }^{+\infty }}\left(\sqrt{1-2{\gamma }^2\mathit{iz}}-1\right)\frac{e^{-\mathit{izu}}}{2\pi }\text{d}z\hfill \\ \hfill & -{\displaystyle {\int }_{-\infty }^{+\infty }}ln\sqrt{1-2i{\gamma }^{2}z}\frac{e^{-\mathit{izu}}}{2\pi }\text{d}z\hfill \end{array}$$

Figure A.1 Integration paths for Fourier decomposition. Solid arrow: original path; broken arrows: equivalent integration path.

where u>0. Introducing the new variable v defined by $\sqrt{1-2{\gamma }^2\mathit{iz}}=-\mathit{iv}$, that is:

$$1-2{\gamma }^2\mathit{iz}=-v^2,-\mathit{iz}=-\frac{v^2+1}{2{\gamma }^2},i\text{d}z=\frac{v\text{d}v}{{\gamma }^2}$$

one obtains:

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_1\left(u,t\right)& =+\frac{t}{2\pi \overline{t}{\gamma }^2}{\displaystyle {\int }_{-\infty }^{+\infty }}\left(\mathit{iv}+1\right)e^{-\frac{v^2+1}{2{\gamma }^2}u}\frac{v\text{d}v}{i{\gamma }^2}\hfill \\ \hfill & =\frac{t{e}^{-\frac{u}{2{\gamma }^2}}}{\sqrt{2\pi }\overline{t}{\gamma }^3\sqrt{u}}{\displaystyle {\int }_{-\infty }^{+\infty }}\frac{\left(v^2-\mathit{iv}\right)e^{-\frac{v^{2}u}{2{\gamma }^2}}}{\sqrt{2\pi }\frac{\gamma }{\sqrt{u}}}\text{d}v\hfill \end{array}$$

and:

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_2\left(u,t\right)& =-\frac{1}{2\pi }{\displaystyle {\int }_{-\infty }^{+\infty }}ln\left(-\mathit{iv}\right)e^{-\frac{v^2+1}{2{\gamma }^2}u}\frac{v\text{d}v}{i{\gamma }^2}\hfill \\ \hfill & =-\frac{e^{-\frac{u}{2{\gamma }^2}}}{2i\pi }{\displaystyle {\int }_{-\infty }^{+\infty }}ln\left(-\mathit{iv}\right)e^{-\frac{v^2}{2{\gamma }^2}u}\frac{v\text{d}v}{{\gamma }^2}\hfill \end{array}$$

The integration path for $\mathrm{\Psi }\left(u,t\right)$ is displayed as a solid arrow in Figure A.1.. Since the integrand is a holomorphic function, an equivalent path is given by the broken arrows. Along the circular arrows, |v| takes constant values, so that the integrals vanish on the arcs at infinity, except in the immediate vicinity of the original path where they diverge. As a consequence, integrating along the axis delivers Hadamard finite value of $\mathrm{\Psi }\left(u,t\right)$.

One therefore obtains for ${\mathrm{\Psi }}_1\left(u,t\right)$:

$$\begin{array}{cc}\hfill \mathfrak{H}{\mathrm{\Psi }}_1\left(u,t\right)& =\frac{t{e}^{-\frac{u}{2{\gamma }^2}}}{\sqrt{2\pi }\overline{t}{\gamma }^3\sqrt{u}}\left[\frac{{\gamma }^2}{u}-i·0\right]\hfill \\ \hfill & =\frac{t{e}^{-\frac{u}{2{\gamma }^2}}}{\sqrt{2\pi }\overline{t}\gamma u^{\frac{3}{2}}}=t·n(u)>0\hfill \end{array}$$

and for ${\mathrm{\Psi }}_2\left(u,t\right)$:

$$\begin{array}{cc}\hfill \begin{array}{cc}\mathfrak{H}{\mathrm{\Psi }}_2\left(u,t\right)\hfill & =-\frac{e^{-\frac{u}{2{\gamma }^2}}}{2i\pi }\{{\displaystyle {\int }_{-\infty }^0}ln\left(i|v|\right)e^{-\frac{v^2}{2{\gamma }^2}u}\frac{v\text{d}v}{{\gamma }^2}\hfill \\ \hfill & +{\displaystyle {\int }_0^{+\infty }}ln\left(-\mathit{iv}\right)e^{-\frac{v^2}{2{\gamma }^2}u}\frac{v\text{d}v}{{\gamma }^2}\}\hfill \\ \hfill & =-\frac{e^{-\frac{u}{2{\gamma }^2}}}{2i\pi }\{\frac{i\pi }{2}{\displaystyle {\int }_{+\infty }^0}{e}^{-\mathit{xu}}\text{d}x-\frac{i\pi }{2}{\displaystyle {\int }_0^{+\infty }}{e}^{-\mathit{xu}}\text{d}x\hfill \\ \hfill & +{\displaystyle {\int }_{-\infty }^{+\infty }}ln\left(|v|\right)e^{-\frac{v^2}{2{\gamma }^2}u}\frac{v\text{d}v}{{\gamma }^2}\}\hfill \end{array}& \hfill \end{array}$$

where $x=\frac{v^2}{2{\gamma }^2}$. The last integral is equal to 0 since the integrand is an odd function of v. Therefore:

$$\begin{array}{cc}\hfill \mathfrak{H}{\mathrm{\Psi }}_2\left(u,t\right)& =-\frac{e^{-\frac{u}{2{\gamma }^2}}}{2i\pi }\left\{\frac{i\pi }{2}{\left[\frac{e^{-\mathit{xu}}}{-u}\right]}_{+\infty }^0-\frac{i\pi }{2}{\left[\frac{e^{-\mathit{xu}}}{-u}\right]}_0^{+\infty }\right\}\hfill \\ \hfill & =-\frac{e^{-\frac{u}{2{\gamma }^2}}}{2i\pi }\left[-2\frac{i\pi }{2u}\right]=\frac{1}{2}\frac{e^{-\frac{u}{2{\gamma }^2}}}{u}.\hfill \end{array}$$

Altogether, one obtains:

$$\begin{array}{cc}\hfill \mathfrak{H}\mathrm{\Psi }\left(u,t\right)=\frac{t{e}^{-\frac{u}{2{\gamma }^2}}}{\sqrt{2\pi }\overline{t}\gamma u^{\frac{3}{2}}}+\frac{1}{2}\frac{e^{-\frac{u}{2{\gamma }^2}}}{u}=t·n(u)+n_0(u)>0.& \hfill \end{array}$$

As a consequence, N(u) is a increasing function of u between 0 and +∞. Note that $\int u^{2}n(u)\text{d}u=\int u^2\frac{e^{-\frac{u}{2{\gamma }^2}}}{\sqrt{2\pi }\overline{t}\gamma u^{\frac{3}{2}}}\text{d}u=\int \frac{u^{\frac{1}{2}}{e}^{-\frac{u}{2{\gamma }^2}}}{\sqrt{2\pi }\overline{t}\gamma }\text{d}u=\frac{2}{\overline{t}}\int \frac{v^2{e}^{-\frac{v^2}{2{\gamma }^2}}}{\sqrt{2\pi }\gamma }\text{d}v$ which is smaller than $\frac{1}{\overline{t}}{\int }_{-\infty }^{+\infty }\frac{v^2{e}^{-\frac{v^2}{2{\gamma }^2}}}{\sqrt{2\pi }\gamma }\text{d}v=\frac{{\gamma }^2}{\overline{t}}$, that is, finite for all finite intervals of integration, as expected. The same applies to $\int u^2{n}_0(u)\text{d}u=\int u^2\frac{e^{-\frac{u}{2{\gamma }^2}}}{2u}\text{d}u=\frac{1}{2}\int u{e}^{-\frac{u}{2{\gamma }^2}}\text{d}u$ which is smaller than $\frac{1}{2}{\int }_0^{+\infty }u{e}^{-\frac{u}{2{\gamma }^2}}\text{d}u=2{\gamma }^4$.

A.4 Verification

Section A.3 only gave the Hadamard finite value of $\mathrm{\Psi }\left(u,t\right)$. Therefore one must evaluate the canonical integral form of the cumulant function, using the previous expression for n(u), in order to check the convergence and obtain the missing terms. One obtains for ${\psi }_1\left(z,t\right)$, since u is always positive:

$$\begin{array}{ccc}\hfill & \hfill \hfill & t{\displaystyle {\int }_{\varepsilon }^{+\infty }}\left(e^{\mathit{izu}}-1-\frac{\mathit{izu}}{1+u^2}\right)n(u)\text{d}u\hfill \\ \hfill & \hfill =\hfill & \frac{t}{\sqrt{2\pi }\overline{t}\gamma }{\displaystyle {\int }_{\varepsilon }^{+\infty }}\left(e^{\mathit{izu}}-1-\frac{\mathit{izu}}{1+u^2}\right)\frac{e^{-\frac{u}{2{\gamma }^2}}}{u^{\frac{3}{2}}}\text{d}u\hfill \\ \hfill & \hfill =\hfill & -\frac{2t}{\sqrt{2\pi }\overline{t}\gamma }\{{\left[\left(e^{\mathit{izu}}-1-\frac{\mathit{izu}}{1+u^2}\right)\frac{e^{-\frac{u}{2{\gamma }^2}}}{u^{\frac{1}{2}}}\right]}_{\varepsilon }^{+\infty }\hfill \\ \hfill & \hfill \hfill & +\frac{1}{2{\gamma }^2}{\displaystyle {\int }_{\varepsilon }^{+\infty }}\left(e^{\mathit{izu}}-1-\frac{\mathit{izu}}{1+u^2}\right)\frac{e^{-\frac{u}{2{\gamma }^2}}}{u^{\frac{1}{2}}}\text{d}u\hfill \\ \hfill & \hfill \hfill & -\mathit{iz}{\displaystyle {\int }_{\varepsilon }^{+\infty }}\left(e^{\mathit{izu}}-\frac{1-u^2}{(1+u^2{)}^2}\right)\frac{e^{-\frac{u}{2{\gamma }^2}}}{u^{\frac{1}{2}}}\text{d}u\}\hfill \\ \hfill & \hfill =\hfill & -\frac{2t}{\sqrt{2\pi }\overline{t}\gamma }\{-\left(e^{\mathit{iz}\varepsilon }-1-\frac{\mathit{iz}\varepsilon }{1+{\varepsilon }^2}\right)\frac{e^{-\frac{\varepsilon }{2{\gamma }^2}}}{{\varepsilon }^{\frac{1}{2}}}\hfill \\ \hfill & \hfill \hfill & +\frac{1}{2{\gamma }^2}{\displaystyle {\int }_{\sqrt{\varepsilon }}^{+\infty }}\left(e^{{\mathit{izv}}^2}-1-\frac{{\mathit{izv}}^2}{1+v^4}\right)\frac{e^{-\frac{v^2}{2{\gamma }^2}}}{v}2v\text{d}v\hfill \\ \hfill & \hfill \hfill & -2\mathit{iz}{\displaystyle {\int }_{\sqrt{\varepsilon }}^{+\infty }}\left(e^{{\mathit{izv}}^2}-\frac{1-v^4}{(1+v^4{)}^2}\right)e^{-\frac{v^2}{2{\gamma }^2}}\text{d}v\}\hfill \\ \hfill & \hfill =\hfill & \frac{2t}{\sqrt{2\pi }\overline{t}\gamma }\left(\mathit{iz}\varepsilon +O(\mathit{iz}\varepsilon {)}^2-\frac{\mathit{iz}\varepsilon }{1+{\varepsilon }^2}\right)\frac{e^{-\frac{\varepsilon }{2{\gamma }^2}}}{\sqrt{\varepsilon }}\hfill \\ \hfill & \hfill \hfill & -\frac{2t}{\overline{t}{\gamma }^2}\left(1-2{\gamma }^2\mathit{iz}\right){\displaystyle {\int }_{\sqrt{\varepsilon }}^{+\infty }}\frac{e^{-\frac{v^2}{2{\gamma }^2}\left(1-2{\gamma }^2\mathit{iz}\right)}}{\sqrt{2\pi }\gamma }\text{d}v\hfill \\ \hfill & \hfill \hfill & +\frac{2t}{\overline{t}{\gamma }^2}{\displaystyle {\int }_{\sqrt{\varepsilon }}^{+\infty }}\frac{e^{-\frac{v^2}{2{\gamma }^2}}}{\sqrt{2\pi }\gamma }\text{d}v+\frac{2\mathit{tiz}}{\overline{t}{\gamma }^2}\hfill \\ \hfill & \hfill \hfill & \times {\displaystyle {\int }_{\sqrt{\varepsilon }}^{+\infty }}\left(\frac{v^2}{1+v^4}-\frac{2{\gamma }^2(1-v^4)}{(1+v^4{)}^2}\right)e^{-\frac{v^2}{2{\gamma }^2}}\text{d}v\hfill \\ \hfill & \hfill =\hfill & \frac{2\mathit{tiz}}{\overline{t}{\gamma }^2}\left[1-\frac{1}{1+{\varepsilon }^2}+O(\varepsilon )\right]\sqrt{\varepsilon }{e}^{-\frac{\varepsilon }{2{\gamma }^2}}\hfill \\ \hfill & \hfill \hfill & -\frac{t}{\overline{t}{\gamma }^2}[\sqrt{1-2{\gamma }^2\mathit{iz}}\mathit{Erfc}\left(\sqrt{\varepsilon (1-2{\gamma }^2\mathit{iz})}\right)\hfill \\ \hfill & \hfill \hfill & -\mathit{Erfc}\left(\sqrt{\varepsilon }\right)]+\frac{2\mathit{tiz}}{\overline{t}{\gamma }^2}F(\varepsilon ).\hfill \end{array}$$