Propagating wave correlations in complex systems

TL;DR

This paper introduces a semiclassical method for computing wave correlation functions in complex systems, linking classical dynamics with wave fluctuations, applicable in quantum and classical wave contexts.

Contribution

It provides explicit algorithms for calculating wave correlation functions using classical dynamics and characterizes fluctuations through ensemble averages, applicable to complex quantum and wave systems.

Findings

Derived explicit expressions relating wave correlation fluctuations to classical correlations.

Applied methods to quantum maps modeling chaotic systems, demonstrating practical utility.

Showed the approach's relevance to both quantum mechanics and classical wave problems.

Abstract

We describe a novel approach for computing wave correlation functions inside finite spatial domains driven by complex and statistical sources. By exploiting semiclassical approximations, we provide explicit algorithms to calculate the local mean of these correlation functions in terms of the underlying classical dynamics. By defining appropriate ensemble averages, we show that fluctuations about the mean can be characterised in terms of classical correlations. We give in particular an explicit expression relating fluctuations of diagonal contributions to those of the full wave correlation function. The methods have a wide range of applications both in quantum mechanics and for classical wave problems such as in vibro-acoustics and electromagnetism. We apply the methods here to simple quantum systems, so-called quantum maps, which model the behaviour of generic problems on Poincar\'e sections. Although low-dimensional, these models exhibit a chaotic classical limit and share common characteristics with wave propagation in complex structures.

Figures

Equations

$$-{\nabla }^2\psi -k^2\psi ={\psi }_0.$$

$$|\psi ⟩=|{\psi }_{-}⟩+|{\psi }_{+}⟩$$

$$(1-\hat{T})|{\psi }_{-}⟩=|{\psi }_{-}^0⟩,$$

$$T(x,x^{\prime };\omega )\approx {\left(\frac{k}{2\pi \text{rmi}}\right)}^{(d-1)/2}\sum _{\stackrel{\mathrm{rays}}{x^{\prime }\to x}}D(x,x^{\prime })\exp \left(\text{rmi}kS(x,x^{\prime };\omega )-\text{rmi}\mu \frac{\pi }{2}\right)$$

$$D(x,x^{\prime })=\sqrt{\left|\det \left(\frac{{\partial }^{2}S}{\partial x\partial x^{\prime }}\right)\right|}.$$

$$\Gamma (x_1,x_2)=⟨x_1|\hat{\Gamma }|x_2⟩,$$

$${\Gamma }_0(x_1,x_2)=⟨{\psi }_0^{-}(x_1){\psi }_0^{-}(x_2{)}^{\ast }⟩,$$

$$\hat{\Gamma }=(1-\hat{T}{)}^{-1}{\hat{\Gamma }}_0(1-{\hat{T}}^{†}{)}^{-1}.$$

$$\hat{\Gamma }=\sum _{n,n^{\prime }=0}^{\infty }{\hat{T}}^n{\hat{\Gamma }}_0({\hat{T}}^{†}{)}^{n^{\prime }}=\hat{K}+\sum _{m=1}^{\infty }\left({\hat{T}}^m\hat{K}+\hat{K}({\hat{T}}^{†}{)}^m\right),$$

$$\hat{K}=\sum _{n=0}^{\infty }{\hat{K}}_n$$

$$W_{\Gamma }(x,p)=\int \mathrm{d}s{\text{rme}}^{-\text{rmi}kps}\Gamma \left(x+\frac{s}{2},x-\frac{s}{2}\right),$$

$$\Gamma (x_1,x_2)={\left(\frac{k}{2\pi }\right)}^{d-1}\int \mathrm{d}p{\text{rme}}^{\text{rmi}k(x_1-x_2)p}{W}_{\Gamma }\left(\frac{x_1+x_2}{2},p\right).$$

$$\mathcal{L}[{\rho }_n](X)={\rho }_{n-1}({\phi }^{-1}(X))=\int \mathrm{d}{X}^{\prime }\delta (X-\phi (X^{\prime })){\rho }_{n-1}(X^{\prime })$$

$$\rho =\mathcal{L}\rho +{\rho }_0\Rightarrow \rho =\frac{1}{1-\mathcal{L}}{\rho }_0=\sum _{n=0}^{\infty }{\mathcal{L}}^n{\rho }_0=\sum _{n=0}^{\infty }{\rho }_n.$$

$$⟨W_{K_n}(X)⟩\approx W_{{\Gamma }_0}\left({\phi }^{-n}(X)\right)={\mathcal{L}}^n\left[{\rho }_0\right](X)={\rho }_n(X),$$

$$⟨W_K⟩=\sum _{n=0}^{\infty }⟨W_{K_n}⟩\approx \sum _{n=0}^{\infty }{\mathcal{L}}^n[{\rho }_0]=\frac{1}{1-\mathcal{L}}{\rho }_0=\rho .$$

$$⟨W_{\Gamma }⟩=⟨W_K⟩\approx \rho .$$

$$⟨\hat{\Gamma }⟩=⟨\hat{K}⟩\approx W^{-1}[\rho ].$$

$${\mathcal{L}}_{\mu }[{\rho }_0](X)$$

$$=$$

$$\hat{T}=v\hat{U},$$

$$W_{{\Gamma }_0}(x,p)=C{\text{rme}}^{-(1-\cos 2\pi (x-x_0))/(2\pi {\sigma }_x{)}^2-(1-\cos 2\pi (p-p_0))/(2\pi {\sigma }_p{)}^2},$$

$$⟨|{\Gamma }_{ij}{|}^2{⟩}_N=\frac{1}{N^2}\text{Tr}{\hat{\Gamma }}^2.$$

$$⟨\text{Tr}{\hat{\Gamma }}^2⟩=\frac{1+v^2}{1-v^2}⟨\text{Tr}{\hat{K}}^2⟩.$$

$$⟨|{\Gamma }_{ij}{|}^2{⟩}_N=\frac{1+v^2}{1-v^2}⟨|K_{ij}{|}^2{⟩}_N$$

$$⟨|W_{\Gamma }(x,p){|}^2⟩=\frac{1+v^2}{1-v^2}⟨|W_K(x,p){|}^2⟩$$

$$\hat{K}=\sum _{n=0}^{\infty }{v}^{2n}{\hat{\Gamma }}_n=\sum _{n=0}^{\infty }{v}^{2n}{\hat{U}}^n{\hat{\Gamma }}_0{\hat{U}}^{-n}$$

$$\mathrm{Tr}({\hat{K}}^2)$$

$$=$$

$$A(n)=\mathrm{Tr}\left({\hat{\Gamma }}_0{\hat{\Gamma }}_n\right)$$

$$\mathrm{Tr}({\hat{K}}^2)=\frac{1}{1-v^4}\sum _{n=-\infty }^{\infty }{v}^{2|n|}A(n).$$

Full text

††: J. Phys. A: Math. Gen.

1 Introduction

There is a long tradition of describing statistical properties of wave fields and spectra in terms of semiclassical or high-frequency techniques and random matrix arguments, both for classical waves such as vibroacoustics [1, 2], electromagnetics [3, 4] and for quantum mechanics [5, 6, 7, 8]. Spectral and wave function statistics are universal under very general conditions and, in particular, if the system under consideration is complex with an associated ray dynamics being chaotic [9]. The statistical properties depend then only on a few generic parameters such as the mean level density, the mean wave amplitude and underlying symmetries such as time-reversal symmetry. Deviations from universality may occur due to non-chaotic dynamical features and for many wave systems in practical engineering applications, these deviations are of great practical importance. In addition, the spatial dependence of mean field values is itself of interest, particularly if absorption or other dissipative effects become relevant.

Classical wave problems – such as the vibro-acoustic response of complex built-up structures or the electromagnetic wave field inside a complex cavity – have in common that they are typically driven by complex, that is, spatially extended but finite, stochastic, and broad-band sources. In addition, they are typically dissipative and are often composed both of regular components (rectangular rooms, walls or corridors, for example) and irregular components (such as cable bundles, circuit boards, moulded vehicle parts or support beams) whose exact location, shape, and topology may be uncertain. Such problems pose computational challenges for full numerical simulation, particularly in the high frequency regime where the scale of a wavelength is small compared to the size of the structure. The inherent uncertainty that one commonly encounters in the intrinsic geometry may in any case make detailed characterisation of the response irrelevant, even if it was computable.

We therefore take an approach in this paper which uses ray propagation to predict averaged, coarse-grained features of the response, as well as providing statistical information. The averages are not sensitive to structural changes, while providing a platform for the prediction of more detailed statistical characteristics in a post-processing step. The connection between ray propagation and wave-field correlations is based here on transfer or boundary operator approaches [10]. Transfer operators allow us to present the problem of wave propagation via multiple reflection in a format which mirrors very closely the phase space mappings used for ray propagation (while in principle allowing an exact solution). This transfer operator approach means that we effectively characterise the problem in terms of a Poincaré-section representation, but we emphasise this lower-dimensional representation can ultimately be mapped to the full problem by using Green function identities (see [11]) for a detailed discussion of propagation of a correlation function from a straight boundary in this context).

Furthermore, an explicit quantitative connection with phase space densities can be made by presenting the correlation function as a Wigner function. The Wigner distribution function (WDF) approach has been developed in the context of quantum mechanics [12], but has found widespread applications also for microwaves [6] and in optics [13, 14, 15]. The method introduced below exploits a connection between the field-field correlation function (CF) and the WDF [16, 17, 18]. Both quantities have been studied intensively in the physics and optics literature. For wave chaotic systems, Berry’s conjecture postulates a universal CF equivalent to correlations in Gaussian random fields [19, 20, 3]. Non-universal corrections can be retrieved by making a link between the CF and the Green function of the system after suitable averaging [21, 22, 23, 24]. The Wigner function formalism can be used to derive the ray-tracing limit for propagating wave fields, see for example [11]. Higher-order (Airy-function) correction to a ray-tracing approach as well as treating evanescent contributions have been discussed in [25, 16, 11].

In this paper, we will use the WDF approach inside finite domains which may be regularly or irregularly shaped. Extending the treatment in [11], we will consider here multiple reflections and interference effects explicitly. We derive the ray-tracing limit in this case and discuss various limits leading to deviations from a pure ray-tracing approach. The ray-tracing component itself can be evaluated efficiently when combined with fast phase-space propagation methods such the Dynamical Energy Analysis (DEA) developed in the context of vibro-acoustics [26] and mesh-based implementation tools such as Discrete Flow Mapping (DFM) techniques [27].

We will focus here on linear and stationary (frequency domain) scalar wave problems using as reference point the wavenumber $k$ or the frequency $\omega$. When considering quantum systems, we may identify the scale $1/k$ with $\hbar$; extensions to vector wave equations is straightforward [28]. In Sec. 2, we introduce the concept of transfer operators and write correlation functions of stationary wave fields in terms of these operators. We derive relations between correlation functions and phase space densities and verify the results with the help of a simple quantum map, the kicked Harper map [29]. In Sec. 3, we extend the results to derive explicit expressions for the variance around the ensemble mean in terms of classical phase space quantities. The results are again validated numerically, here for a fully chaotic quantum map, the perturbed cat map [30].

2 Transfer operator formulation and correlation functions

2.1 The transfer operator $\hat{T}$

The wave problem under consideration is described in an operator formulation using transfer operators $\hat{T}$ defined on a $(d-1)$-dimensional manifold, or surface of section (SOS), where $d$ is the dimension of the underlying space. As a concrete example, we consider the Helmholtz equation in a $d$-dimensional region $\Omega\subset\mathbb{R}^{d}$

[TABLE]

Boundary conditions given on $\partial\Omega$ are assumed to take the form of a linear relationship between the solution $\psi$ itself and its normal derivative $\partial\psi/\partial n$. We will in general use $\partial\Omega$ as the SOS although other choices are possible. We include an inhomogeneous source term $\psi_{0}$ here, which drives the wave dynamics. An operator formulation on the boundary is naturally given in terms of boundary integral equations, which yield relations between the wave function $\psi$ and its normal derivative on the boundary. For our purposes, it is more convenient to follow [10] and decompose the boundary field into incoming and outgoing components

[TABLE]

as depicted in Fig. 1. Note that we use the ket notation for the solutions restricted to the SOS given here by $\partial\Omega$.

The solution to the inhomogeneous wave problem can now be cast in terms of incoming waves at boundaries,

[TABLE]

where $\ket{\psi_{-}^{0}}$ is the source wave field incoming on the boundary. Transfer operators can be defined for generic wave systems and for generalised SOSs other than $\partial\Omega$, see [31, 32]. An explicit construction of the exact $\hat{T}$ is, however, often not straightforward and implementations have been presented for model systems only [33, 34, 35, 10]. General expressions for $\hat{T}$ can be given explicitly in the semiclassical limit $k$ (or $\omega$) $\to\infty$ in terms of the so-called Bogomolny transfer operator [31, 36],

[TABLE]

with

[TABLE]

The sum in (4) is over all rays passing through the interior to cross from point $x^{\prime}$ to point $x$ on the $(d-1)$ dimensional SOS. The phase $S(x,x^{\prime})$ is the classical action at a fixed frequency $\omega$ of the ray and $\mu$ is a phase index arising due to boundary conditions or ray caustics; $S$ can be complex if there is damping. Note that, when considering the Helmholtz case in a domain $\Omega$, there is only one such ray trajectory from $x^{\prime}\to x$. The wave number $k$ is introduced here for convenience; in the special case of the Helmholtz equation (1), we have $S(x,x^{\prime})\equiv k\,L(x,x^{\prime})$, the optical length of the chord connecting $x^{\prime}$ to $x$.

2.2 The correlation function $\hat{\Gamma}$

We are primarily interested in classical wave problems driven by noisy, stochastic sources. In this context, the natural solution of the problem is in terms of a two-point correlation function

[TABLE]

where $\hat{\Gamma}$ can be interpreted as a density matrix of the system and $x_{1}$ and $x_{2}$ are coordinates on the SOS. In the simplest case of coherent driving, we may consider $\hat{\Gamma}=\ket{\psi^{-}}\bra{\psi^{-}}$, with $\ket{\psi^{-}}$ a solution of the form (2) for a specific source term and at a given frequency. More generally, we consider driving by stochastic sources and correspondingly density matrices obtained from, for example, frequency, spatial or time averages (when considering time dependent, in particular, non-harmonic, stochastic sources) [37]. Note that we always assume here that $\hat{T}$ itself is not explicitly time dependent, even if the driving is.

Consider a source correlation function $\Gamma_{0}$ with

[TABLE]

where $\psi_{0}^{-}(x)$ is the source wave field as defined in (2) and $\braket{.}$ denotes an ensemble average over time intervals, frequency or local spatial averaging. Note, that while the source distribution $\ket{\psi_{0}^{-}}$ show strong spatial fluctuations, the averaged quantity ${\Gamma_{0}}$ can be a fairly smooth function of both $x_{1}$ and $x_{2}$. Starting from (3), we write the system correlation function, including multiple reflections at boundaries, in terms of the source correlation function and the transfer operator $\hat{T}$. That is,

[TABLE]

Formally, one can write this in terms of the source correlation function $\Gamma_{0}$ being propagated along outgoing waves and undergoing multiple reflections at boundaries. After geometrically expanding the right hand side of Eq. (7), we write

[TABLE]

where the $n$th iteration of the operator $\hat{T}^{n}$ is related to waves undergoing $n$ reflections at $\partial\Omega$. The operator $\hat{K}$ in (8) contains the diagonal contribution of the double sum and is defined as

[TABLE]

It contains the smooth part of the correlation function as will be shown below. The terms $\hat{K}_{n}$ can be interpreted as the $n$th iteration or reflection contribution to the smooth part of the correlation function. Note that some damping is implicitly assumed here in order for these sums to converge.

We finally remark that, once the correlation function has been characterised in a boundary representation, as set out in this section, one can use Green function identities to propagate this information to the interior. See [11] for a more detailed discussion of such propagation of correlation functions from a straight boundary, including the treatment of evanescent components, for example.

2.3 Relation to classical phase space densities

The quantities $\hat{\Gamma}_{0}$, $\hat{\Gamma}$ and $\hat{K}$ can be related to phase space densities using Wigner transformation after additional averaging as demonstrated below and in A. The WDF of an operator $\hat{\Gamma}$ is defined as

[TABLE]

with back transformation given by

[TABLE]

Here, $k$ represents a wave number such as introduced in (1), (4) and $d$ is the dimension of the full, interior problem. It is shown in A, that for sufficiently smooth initial phase space distributions $W_{\Gamma_{0}}=\rho_{0}$, the averaged Wigner transform of $\hat{K}_{n}$ in (9) is given in terms of the classical flow equations, see also [11, 38].

Classical phase space densities are driven by the phase space dynamics, in our case a boundary map or more generally the Poincaré map of an SOS. The initial density $\rho_{0}=W_{\Gamma_{0}}$ can then be associated with a boundary density of rays arriving directly from a source distribution in the interior. Mapping the source ray density through subsequent reflections leads to the iterated densities $\rho_{0}\to\rho_{1}\to\cdots\rho_{n}\to\cdots$ which can be described in terms of the (linear) integral operator $\cal L$ defined in the lossless limit as [39]

[TABLE]

(see below for the treatment of losses). Here, $X=(x,p)$ denotes the collective phase space coordinates on the SOS (with $p$ the momentum conjugate to $x$) and $\varphi:X^{\prime}\to X$ is the classical map describing the flow of trajectories from the SOS back to itself after a single reflection. We use that $\varphi$ is Hamiltonian, and, in particular, phase space volume preserving. The operator $\cal L$ is also referred to as the Frobenius-Perron (FP) operator [39]. The integral representation in (12) is useful for considering effects like absorption and mode conversion as well as uncertainty, see [40]. The ray-dynamical, or classical, analogues of equations (9) are now provided by

[TABLE]

From the relation derived in A, we can now write

[TABLE]

and we obtain for the Wigner function of the full $\hat{K}$ operator to leading order

[TABLE]

The averaging $\langle\cdot\rangle$ is understood as defined in (6) and can be performed in terms of an average over an ensemble of similar systems, appropriately chosen frequency averaging or local (spatial) averaging, for example. In addition, it is assumed that the initial density $\rho_{0}=W_{\Gamma_{0}}$ is a smooth function on the scale of $\Delta x\Delta k=1$.

For the full correlation function, we find that contributions to (8) with $n\neq n^{\prime}$ are removed by averaging so that we may assert in addition that

[TABLE]

The quantities $\hat{\Gamma}$ and $\hat{K}$ thus have the same mean when taking a suitable average and this mean value is given by the classical equilibrium (phase space) density $\rho$ obtained from (15), that is,

[TABLE]

We will find in the next section, however, that fluctuations of $\hat{\Gamma}$ about this mean are much stronger than those of $\hat{K}$, particularly in the limit of weak damping.

The averaged quantities in (17) contain detailed information about the system. Powerful numerical tools have been developed for computing $\rho$, and thus indirectly for the mean of $\hat{\Gamma}$. Among those, the DEA technique together with the DFM implementation on meshes is particularly suitable for complex structures, which has been used in a range of engineering applications [26, 27, 41].

Note that, without losses, the FP operator $\cal L$ has a leading eigenvalue 1 and the solution to Eq. (15) diverges. It is thus necessary - in order to arrive at mathematically and physically meaningful solutions - to account for losses such as naturally occur in wave systems due to wall absorption or absorption in the interior of $\Omega$ or due to radiation through apertures. Absorption can be included in a FP operator formalism by introducing, for example,

[TABLE]

where we require that the damping coefficient $\mu\geq 0$ is additive along the flow, mimicking distributed losses in many real-life engineering systems. For interior damping, such as for acoustic waves in air, we typically have $\mu(X)=\mu_{0}\,L(X)$, where $L$ is the length of a trajectory segment between two reflections on the boundary starting at $X$ and $\mu_{0}\geq 0$ is a real constant.

2.4 Numerical illustration of averaged response using quantum map models

We will demonstrate the validity of the relations described in the last section with the help of a simple model in which the transfer operator is simulated by

[TABLE]

where $\hat{U}$ is a unitary operator acting on a space of dimension $N$ and the prefactor $v$, satisfying $0<v<1$, accounts for dissipation. In numerical investigations here and in following sections, we choose $\hat{U}$ to be the quantum analogue of a classical map on the unit torus (see B for further details). By choosing maps with different levels of chaotic behaviour, we aim to simulate the response of larger, complex systems governed by classical wave theories. In these models, the dimension $N$ takes on the role of the wave number $k$ as the large parameter.

In order to illustrate the average response set out in the preceding discussion, we let $\hat{U}$ take the form of a quantum kicked Harper (QKH) map in this section, see [42, 43, 44] or the review [29] (other maps are used when discussing fluctuations about the mean later). The corresponding classical kicked Harper map exhibits a predominantly chaotic dynamics but with small regular islands. Although these islands are small, their presence leads to a significant degree of tanglement of the stable and unstable manifolds as demonstrated in Figs. 2c and 2d. This gives rise to a strong dependence on initial source distributions when considering correlation functions, which makes it a good test model for studying relations between classical and wave operators.

The source $\hat{\Gamma}_{0}$ is chosen so that the associated Wigner distribution function takes the form

[TABLE]

where $C$ is a normalisation constant chosen so that $\Tr\Gamma_{0}=1$ and $\sigma_{x}$ and $\sigma_{p}$ respectively determine the variances in the $x$ and $p$ coordinates. For small enough $\sigma_{x}$ and $\sigma_{p}$ this is a nearly Gaussian source distribution in phase space centred on $X_{0}=(x_{0},p_{0})$, and is associated with a corresponding nearly Gaussian source correlation function $\Gamma_{0}(x_{1},x_{2})$.

The time evolution of the Wigner transform $W_{K_{n}}$ of $\hat{K}_{n}$ is displayed in Fig. 2 and is compared with the evolution of a classical phase space density with initial distribution according to (20). The quantum evolution follows the classical dynamics closely for small $n$ as expected from Eq. (12) and detailed in A. Once the classical dynamics foliates phase space down to phase space cells of area $1/N$, wave fluctuations take over and wash out the classical phase space structures, see Fig. 2d at $n=10$. Note that, after several interations of the map, the classical density is stretched along the unstable manifolds of the system and tracks, for example, the folds and switchbacks in them that arise from the small islands present. Some of these features are also seen in the propagated Wigner function but the details are lost.

In Fig. 3, we compare the classical (stationary) phase space distribution $\rho$ defined in (13) with the Wigner transforms of $\hat{K}$ and $\hat{\Gamma}$ for a damping value $v=0.9$. The classical and wave results coincide remarkably well for $\rho$ and $W_{K}$ even without further averaging, see Fig. 3a, while only traces of the classical phase space structure are left for $W_{\Gamma}$ in this case. When performing a further parameter average in the Harper map, the quantum fluctuations are suppressed and the classical phase space structure comes out more clearly, both for $W_{K}$ and $W_{\Gamma}$, see Fig. 3b. The fluctuations in $W_{\Gamma}$ are in both cases much larger than in $W_{K}$, with similar differences to be expected for the original objects $\hat{\Gamma}$ and $\hat{K}$. We can in fact quantify the increased fluctuations in $\hat{\Gamma}$ compared to $\hat{K}$ as a function of $v$, as will be described in the next section.

The results show that the classical phase space density gives the mean value for both quantities $W_{K}$ and $W_{\Gamma}$ as expressed in (16). It also shows how well the classical dynamics describes a wave operator such as $W_{K}$ down to a detailed description of the complex phase space structure present in a mixed system like the kicked Harper map. Note, that the wave fluctuations due to higher iterates of the map are suppressed here due to the 10% damping ($v=0.9$) introduced.

3 Fluctuations of $\hat{K}$ and $\hat{\Gamma}$ about the mean

Although $\hat{K}$ and $\hat{\Gamma}$ have the same average, the fluctuations about the mean are much greater for $\hat{\Gamma}$ than for $\hat{K}$ as illustrated in Fig. 3, an effect which is greatly amplified in the weak damping limit. In this section we quantify this effect for the simple model $\hat{T}=v\hat{U}$ introduced in Sec. 2.4 by calculating mean values for $\Tr\hat{\Gamma}^{2}$ and $\Tr\hat{K}^{2}$. Formally, these also provide us with variances of individual matrix elements $K_{ij}$ and $\Gamma_{ij}$ in any given basis, using, for example,

[TABLE]

Here $\langle\cdot\rangle_{N}$ denotes an average over the $N$-dimensional basis on the space on which $\hat{T}$ is defined.

Explicit expressions can be given for these quantities, as outlined in the remainder of this section. In particular, if we average over a parameter such as frequency (or system dimension in the quantum map model), we can demonstrate that fluctuations in $\hat{\Gamma}$ scale simply with fluctuations in $\hat{K}$, according to

[TABLE]

This expression is quite general and is equally valid in chaotic and integrable limits, for example. It is also true irrespective of whether the system has time reversal symmetry or not. It provides a direct quantitative measure of the greater fluctuations seen in $\hat{\Gamma}$, relative to those of $\hat{K}$, as $v\to 1$. Furthermore, we will see that the fluctuations of $\hat{K}$ itself can be obtained using information that is readily available from propagation of a corresponding classical density, as performed in complex structures using the DEA method, for example.

Using (21), the scaling relationship (22) provides an equivalent scaling

[TABLE]

of the relative sizes of corresponding matrix elements. It should be noted that in practice such averages over all matrix elements are dominated by diagonal elements in the weak damping limit or when the source intensity is extended over length scales much larger than a wavelength: however, it can be shown using alternative but more involved calculations, not reported here, that these results can also be applied element by element. Finally, any statement involving such traces can be expressed also as a relation between Wigner functions. For example, (22) can also be expressed in the form,

[TABLE]

using standard properties of Wigner functions. This relation provides an explicit quantitative statement of the qualitative observations made in Fig. 3.

3.1 The variance of the diagonal part $\hat{K}$

We consider $\hat{K}$ first and write Eq. (9) as

[TABLE]

with $\hat{\Gamma}_{n}=v^{-2n}\hat{K}_{n}$ representing the smooth $n$-step correlation function as in Eq. (9). Using the cyclic permutability of the trace we may write,

[TABLE]

where

[TABLE]

is the $n$-step return probability. Note that ${A}(n)={\rm Tr}\hat{\Gamma}_{0}\hat{\Gamma}_{n}={\rm Tr}\hat{\Gamma}_{n}\hat{\Gamma}_{0}={\rm Tr}\hat{\Gamma}_{0}\hat{\Gamma}_{-n}={A}(-n)$. After reordering the double sum, this may be written in the form

[TABLE]

So far the calculation is exact and completely independent of the classical limit.

The quantity ${A}(n)$ can again be related to the classical phase space dynamics using

[TABLE]

We now make the approximation used in (12) that, on averaging, the fluctuations in $A(n)$ are removed and that it may be approximated by its classical or ray-dynamical analogue, the classical phase space autocorrelation function

[TABLE]

Here, $\rho_{0}$ is the correspondingly averaged classical phase space density of the source and $\rho_{n}={\cal L}^{n}\rho_{0}$ its $n$th-reflection iterate. Then

[TABLE]

For strongly chaotic systems, the autocorrelation function decays exponentially as

[TABLE]

where $\gamma_{\rm cor}$ denotes the correlation exponent and $a_{0},a_{1}$ are constants depending on the initial density $\rho_{0}$ ($\rme^{-\gamma_{\rm cor}}$ typically being the 2nd largest eigenvalue of the FP operator.) Note that, for strong damping, it suffices to know the autocorrelation function $A(n)$ for relatively small $n$ in order to compute $\Tr(\hat{K}^{2})$. In this case one finds that the exact $A(n)$ may be well approximated by its ray-dynamical analogue $A_{\rm cl}(n)$ even without averaging and that the result above then also holds for individual maps. (Of course, the trace operation itself performs an average over basis states even if no further averaging is performed.)

3.2 The variance of the correlation function $\hat{\Gamma}$

To get an equivalent expression for the variance of the full correlation function of the stationary wave field, $\hat{\Gamma}$, we start from Eq. (8), giving

[TABLE]

from which we obtain (using cyclic permutability of the trace again)

[TABLE]

We argue next that, after averaging,

[TABLE]

due to the phase difference that one necessarily finds in contributions to $\hat{U}^{n}$ and $\hat{U}^{n^{\prime}}$ when $n\neq n^{\prime}$. On average, the sum in (31) is then dominated by terms with $n_{1}-n_{4}=n_{2}-n_{3}=n$ and, after reordering, we may write

[TABLE]

as previously reported in (22). Note that this condition has been derived based on neglecting interfering contributions which are wiped out by appropriate averaging. We have not made any further assumptions about the underlying classical dynamics and the relations (22) – (24) are universal, independent of whether the system dynamics is regular, chaotic or mixed. In particular, the variance of $\hat{\Gamma}$ always exceeds that of $\hat{K}$, only approaching the same value in the limit of strong damping, $v\to 0$.

We note that relation between wave fluctuations (as a function of energy or frequency) and classical decay rates is a well established concept and has been discussed in the context of fluctuations in scattering cross sections [45, 46] or conductance fluctuations for transport through open chaotic cavities [47], see also [48]. The relevant classical quantity is here the classical escape rate $\gamma_{\rm esc}$, which in our setting corresponds to the rate of dissipation $\gamma_{esc}=-\log v$. A generalised treatment involving the full spectrum of the FP operator is given in [49]. Our result highlights the influence of higher order eigenvalues of the FP operator on the variances in Eqs. (27) and (32) such as through the classical decay of correlation contributions, Eq. (30). The main result of this section is, however, the relation (33) relating the fluctuations in the ’smooth’ part, $\hat{K}$, to the fluctuations in the total correlation function. Together with (17) and (29), we can now relate the first and second moments of the distributions of both $\hat{\Gamma}$ and $\hat{K}$ to classical phase space observables such as the stationary phase space density $\rho$ and the classical phase space autocorrelation function $A_{\rm cl}(n)$. These quantities depend of course on the underlying classical dynamics.

3.3 Numerical illustration of fluctuations using quantum map models

The fluctuations of $\hat{K}$ and $\hat{\Gamma}$ are now illustrated using a quantum map model. As asserted in the preceding discussion, averages (27) and (32) and the corresponding variances are insensitive to the underlying symmetries and to how chaotic or integrable the system is, although the detailed distributions are not.

In this section, we use a perturbed cat map [30], in which a quantum cat map [50, 51, 52] is perturbed with a QKH map of the form used in Sec. 2.4 (see B for more details). The perturbation is large enough to break any underlying time-reversal or spatial symmetries of the quantum version of the cat map, but small enough that the overall classical dynamics is still completely chaotic. A source term $\Gamma_{0}$ is used which corresponds to the nearly Gaussian density given in (20), with $\sigma_{x}=\sigma_{p}=1/2$. It is centred on a period-one fixed point $X_{0}=(x_{0},p_{0})$ near the origin of phase space, which slows the short-time decay of the autocorrelation function $A(n)$.

Corresponding numerically computed values of $\Tr\hat{K}^{2}$ and $\Tr\hat{\Gamma}^{2}$ are shown as circles and triangles respectively in Fig. 4 as $v\to 1$. As $v$ is varied, the dimension $N$ of the quantum map is also changed, so that typical variations can be seen in individual values of $\Tr\hat{K}^{2}$ and $\Tr\hat{\Gamma}^{2}$ around the averages (27) and (32), presented in the figure as solid curves. Over the range shown, individual values of $\Tr\hat{K}^{2}$ follow the average prediction (27) quite closely. The trace operation is effectively self-averaging in this case. In contrast, although (32) is found to be a good predictor of average behaviour, individual values of $\Tr\hat{\Gamma}^{2}$ are seen to fluctuate significantly about the mean, and to an increasing extent as $v$ approaches one.

In fact, it is found in Fig. 4 that most individual values of $\Tr\hat{\Gamma}^{2}$ fall significantly below the mean (32) when $v$ is very close to one, with the average being achieved because increasingly rare individual cases arise with exceptionally large values. This can be understood simply in terms of the physics of the weakly damped resonant response of the system and is illustrated in more detail in Fig. 5. Here we present fluctuations with changing $N$ of $\Tr\hat{\Gamma}^{2}$ for the particular values $v=0.995$ and $v=0.999$ of the damping parameter, represented respectively by crosses and circles. The smallest individual values of $\Tr\hat{\Gamma}^{2}$ are insensitive to $v$ in this illustration. Physically, these correspond to instances of the system being off resonance and changing only slightly as $v$ makes the final steps towards unity. In the weak damping limit, where resonances are very narrow, this accounts for most parameter values. However, the minority of cases which are at or near resonance respond very wildly to changing $v$ and become dramatically larger as $v$ is increased. These instances of resonant response bring the average value of $\Tr\hat{\Gamma}^{2}$ up to its predicted level (32), represented in Fig. 5 by the horizontal lines.

Clearly then, the average fluctuations predicted in this paper present only a partial characterisation of the system response in the weak damping limit. A complete description demands a characterisation of the distribution of values. This lies outside the scope of the current discussion but will be reported in a future publication. The calculations in this section do provide, however, a simple characterisation of the variability of the response of the system about the mean.

4 Conclusion

Ray tracing methods often provide the only feasible means of treating wave propagation in the high-frequency limit. Recently developed DEA methods have in particular provided an effective means of implementing such phase-space transport in very large, complex structures. They provide an effective description of the local averaged response of large systems to driving by sources that are themselves also often complex and statistically characterised. However, such methods used in isolation miss important physical phenomena related to interference and to higher order effects such as diffraction.

The aim of this paper has been to establish and exploit an entirely wave-based analogue of this phase-space transport problem so that wave effects such as multi-path interference can be incorporated into large-scale simulations with complex and noisy forcing. We have argued that an effective platform for the calculation of such wave effects can be built on the relationship between correlation functions and Wigner functions that has been established in the contexts of quantum mechanics and optics. Two-point correlation functions provide an effective means of exploiting available information about spatial localisation and directionality of waves radiated from a noisy, complex source. They have proven to be an effective way of characterising EM emissions from electronic circuitry through direct measurement, see for example [37]. Correlation function propagators then provide a completely wave-based analogue of phase-space transport approaches such as DEA.

Importantly, this allows us to provide a statistical description of fluctuations due to interference in the response of a forced wave problem, using only information that is readily available from a direct phase-space simulation. In particular, the global variance of the wave correlation function can be described in terms of an autocorrelation function of propagated phase-space densities. In essence this allows us to boot-strap phase-space transport simulations to predict fluctuation about the mean of the response of the system, as well as the mean itself. The approach has been tested on simple quantum map models based on a representation of wave transport by boundary transfer operators, but using a framework that we believe will scale up effectively to much larger systems. Achieved results are relevant in the statistical characterisation of large vibro-acoustics and electromagnetic structures, including reverberation chambers operated at arbitrary frequencies.

We gratefully ackowledge support from from EPSRC under grant number EP/K019694/1, from the EU FP7 project MHiVec and from the EU Horizon 2020 network NEMF21.

Appendix A The evolution of correlation functions

We will give here a derivation of Eq. (14) valid in the semiclassical limit (see also [38]). Starting from Eq. (9), we write $\hat{K_{n}}=\hat{T}^{n}\hat{\Gamma}_{0}(\hat{T}^{n})^{\dagger}$ in $x$ representation as

[TABLE]

where we set $K_{0}=\Gamma_{0}$ for convenience. We employ the short-wavelength approximation

[TABLE]

of the operator $\hat{T}^{n}$, where the sum index $\alpha$ runs over all trajectories from $x^{\prime}$ to $x$ that encounter the SOS $n$ times, while other quantities are defined as in (4) and (5) with obvious modifications. For notational compactness, the phase factor $\mu$ is assumed to be absorbed into the definition of the amplitude here. Using the transformation

[TABLE]

and starting from the semiclassical expression (34), we obtain

[TABLE]

where $\alpha$ and $\beta$ label $n$-bounce orbits from $x_{1}^{\prime}$ to $x_{1}$ and from $x_{2}^{\prime}$ to $x_{2}$, respectively.

At this point, we need to make two crucial assumptions about the quantities of interest. First, we concentrate on the averaged response, so that orbit combinations with topologically distinct $\alpha$ and $\beta$, which arrive with significant noncancelling phase differences, are washed out. We are then left just with the diagonal contributions in which an orbit $\alpha$ coincides with an orbit $\beta$ up to a slight deformation. In the following we therefore set $\alpha=\beta$. Second, we assume that the source has a sufficiently short correlation length that propagation can be approximated using Taylor expansions of the surviving amplitude and phase contributions, truncated at the terms

[TABLE]

and

[TABLE]

Here, $p_{\alpha}^{\prime}(x,x^{\prime})=-\partial S_{\alpha}(x,x^{\prime})/\partial x^{\prime}$ and $p_{\alpha}(x,x^{\prime})=\partial S(x,x^{\prime})/\partial x$ [53] are respectively the initial and final momenta of the $n$-bounce trajectory from $x^{\prime}$ to $x$ labelled by $\alpha$. The condition that the correlation length is sufficiently short translates in the Wigner representation to the condition that the dependence on $p$ of the Wigner function is sufficiently slow. Note that even very smooth initial densities wrinkle under chaotic evolution so that rapid oscillations develop as $n$ increased. For a fixed underlying dynamics, the truncated Taylor expansions above may therefore fail for moderately large $n$. By using the derived identities below for longer times we are therefore implicitly assuming that there is also sufficient averaging over system parameters that very fine-scale classical features are smoothed out in propagated densities.

With these assumptions we obtain

[TABLE]

Using the definition of $D_{\alpha}$ given in (5), we may write

[TABLE]

where on the left we sum over all orbits arriving at $x$, labelled by initial position $x^{\prime}$ and topology $\alpha$. On the right we reformulate the same sum as a simple integration over the momentum $p$ at arrival (in which there is no need for the label $\alpha$ since $x$ and $p$ uniquely determine the orbit topology). Using this reformulation of the sum we can write

[TABLE]

with $W_{K_{0}}(x^{\prime},p^{\prime})$ denoting the Wigner transform of $K_{0}$ as defined in (10). Note that $X^{\prime}=(x^{\prime},p^{\prime})$ in (37) is now a function of $X=(x,p)$ through the relation $X=\varphi^{n}(X^{\prime})$, with $\varphi$ defining the dynamics on the SOS (see Sec. 2.2). After applying the Wigner transformation on both sides of Eq. (37) and evaluating the resulting $\delta$-function, we obtain

[TABLE]

which mirrors the action of the classical FP operator.

Appendix B Conventions for quantum maps

Here we summarise the quantum maps used in Secs. 2.4 and 3.3 to test correlation function propagation. We use maps defined on a toral phase space with unit period in each of the phase space coordinates $x$ and $p$.

We begin with the classically-defined kicked map defined by

[TABLE]

where the primes on $F$ and $G$ denote differentiation. This map can be viewed as being the result of using unit-time flow under the Hamiltonian $G(p)$, followed by unit-time flow under the Hamiltonian $F(x)$. The illustrations in Figs. 2-3 assumed $G^{\prime}(p)=a\sin 2\pi p$ and $F^{\prime}(x)=b\sin 2\pi x$ and $a=-2b=2/\pi$, in which case the map is a kicked Harper map, which will be the nomenclature we use henceforth. These maps provide generic examples with predominantly chaotic dynamics mixed with small regular islands for the chosen parameter values.

The corresponding quantum map is then defined on a Hilbert space of dimension $N=1/(2\pi\hbar)$ such that

[TABLE]

The geometry of phase space is reflected in the detailed quantisation of the position and momentum operators $\hat{x}$ and $\hat{p}$, and on the boundary conditions imposed on quantum states. We can use a position basis $\hat{x}\ket{x_{i}}=x_{i}\ket{x_{i}}$, with quantised positions $x_{i}=(i+\alpha_{x})/N$, whose index $i$ runs over $i=0,\cdots,N-1$. Alternatively, we can use a momentum basis $\hat{p}\ket{p_{i}}=p_{i}\ket{p_{i}}$ on the same index set with quantised momenta $p_{i}=(i+\alpha_{p})/N$ and related to the position basis by $\langle{x_{l}}|{p_{j}}\rangle=\rme^{2\pi N\rmi x_{l}p_{j}}/\sqrt{N}$. The shifts $\alpha_{x}$ and $\alpha_{p}$ are determined by the boundary conditions satisfied by states under translation across the torus. The direct wavefunction $\psi(x_{i})=\langle x_{i}|\psi\rangle$ satisfies $\psi(x_{i}+1)=\rme^{2\pi\rmi\alpha_{p}}\psi(x_{i})$. The momentum representation $\varphi(p_{i})=\langle p_{i}|\psi\rangle$ satisfies $\varphi(p_{i}+1)=\rme^{-2\pi\rmi\alpha_{x}}\varphi(p_{i})$. All of the maps used in this paper use either $\alpha_{x}=\alpha_{p}=0$ or $\alpha_{x}=\alpha_{p}=1/2$ and, in particular, Figs. 2-3 assumed $\alpha_{x}=\alpha_{p}=0$.

Alternatively, cat maps and their quantisations provide examples of fully chaotic, hyperbolic dynamical systems. Quantisations of the simple cat map exhibit nongeneric degeneracies but these can be removed by perturbations which retain the fully chaotic dynamics. In particular, we use a perturbed cat map

[TABLE]

where $\hat{U}_{\rm QKH}$ is a quantum kicked Harper map and $\hat{U}_{\rm C}$, defined by

[TABLE]

quantises the (unperturbed) cat map

[TABLE]

This quantisation is well-defined for even $N$ with $\alpha_{x}=\alpha_{p}=0$ and for odd $N$ with $\alpha_{x}=\alpha_{p}=1/2$. The QKH map used for the numerical illustration in Sec. 3.3 was also of this form with $G^{\prime}(p)=a\sin 2\pi p$, $F^{\prime}(x)=-b\sin 2\pi x$ and $a=b=0.1$. This choice completely eliminates symmetries in the combined map, while presenting a small enough perturbation of the cat map that the dynamics is fully chaotic.