Wave Motion Uncertainty quantification for phase-space boundary integral models of ray propagation

Vibrational and acoustic energy distributions of wave fields in the high-frequency regime are often modeled using flow transport equations. This study concerns the case when the flow of rays or non-interacting particles is driven by an uncertain force or velocity field and the dynamics are determined only up to a degree of uncertainty. A boundary integralequationdescriptionofwaveenergyflowalonguncertaintrajectoriesinfinitetwo- dimensional domains is presented, which is based on the truncated normal distribution, and interpolates between a deterministic and a completely random description of the tra- jectory propagation. The properties of the Gaussian probability density function appearing in the model are applied to derive expressions for the variance of a propagated initial Gaussiandensityintheweaknoisecase.Numericalexperimentsareperformedtoillustrate these findings and to study the properties of the stationary density, which is obtained in the limit of infinitely many reflections at the boundary. linearly propagation speed c and the parameters σ s and σ p considered here. A similar analysis is also performed for an uncertain boundary location.


Introduction
Noise and vibration simulations for mechanical structures are commonly based on numerical solutions of linear wave equations. Popular methods include finite element methods, finite volume methods, boundary element methods or a variety of spectral methods. However, numerically approximating the solutions of wave equations directly has two potential drawbacks. Firstly, the number of degrees of freedom in the numerical model needs to be increased as the frequency is increased to maintain a fixed accuracy level. If the local wavelengths are significantly smaller than the dimensions of the physical system, then this requirement can become computationally prohibitive. Secondly, uncertainties play a more important role at higher frequencies when the structural modes become sufficiently dense that they can switch positions due to small structural differences within standard manufacturing tolerances. When the resonance shifts become larger than the mean distance between the structural modes then wave based models of a vibrating system deviate from experimental data, no matter how accurate the numerical wave solver. These issues can be addressed by instead using statistical methods, such as Statistical Energy Analysis (SEA) [1], to predict averaged energy distributions as well as higher order statistical moments via random matrix theory [2]. However, SEA only provides a coarse description of the modeled structure since constant energy levels are assumed throughout relatively large substructures and the validity of such a model can therefore be difficult to ascertain a-priori [3].
The vibrational or acoustic energy distributions predicted in SEA can be modeled with a higher spatial resolution by considering high-frequency wave energy transport in terms of a geometrical ray description [4]. Wave effects such as interference and diffraction are not included in the model, which is recast in terms of tracking densities of rays or particles in phase-space. The computational cost of directly tracing rays or beams from a source to a receiver scales with the number of paths that must be modeled. For complex domains including many reflections this can lead to unfeasibly large computational models. For example, applications in room acoustics are typically limited to (at most) second order reflections [5]. In builtup mechanical structures it may also be necessary to include mode conversion and refraction in the model. Here, indirect methods based on propagating ray densities (instead of the rays themselves) through phase-space such as the so-called Dynamical Energy Analysis (DEA) method [4,[6][7][8] can provide a more practical alternative. The deterministic propagation of ray densities through phase-space is governed by the Liouville equation, the solution of which can be expressed in terms the so-called Frobenius-Perron operator [9]. In general, this linear integral operator transports ray densities along the trajectories of a dynamical system and is often termed the transfer operator. These ideas have also found their way into the literature in computer graphics [10] and room acoustics [11] amongst others, where the corresponding transfer operator equation is more commonly known as the rendering equation.
In this work we consider wave energy transport through uncertain two-dimensional domains in the high-frequency ray limit, which leads to a stochastic velocity field driving the energy transport. We discuss the propagation of these uncertainties using stochastically smoothed transfer operators [9,[12][13][14]. These operators propagate densities along stochastic ray trajectories, where both the arrival position of the ray and the nature of its reflection at the boundary will be considered as uncertain. Physically, these uncertainties may be due to an uncertain boundary geometry, including rough boundary scattering, or propagation through media with random inhomogeneities. This gives rise to both uncertain ray trajectories and uncertainties in the reflected direction, leading to probabilistic positions and momenta. Typically, the mean transported position and reflected direction relate to those of the underlying deterministic ray flow with specular reflections.
Stochastic transfer operators have also been employed more widely for modeling probability density distributions ρ as they are transported along trajectories X in phase-space [15]. Assuming X (t) is a stochastic process defining the phase-space coordinates of the trajectory at time t, then it can be described via a stochastic dynamical system of the form dX (t) = V(X (t))dt + σ dW (t), (1) where σ > 0, W is a phase-space Wiener process and V is a (deterministic) vector field. The solutions of this dynamical system X are related to the probability density ρ via For sufficiently smooth vector fields V, then conservation of probability leads to a continuity equation of the form known as the Fokker-Planck or forward Kolmogorov equation. The solution of this equation in R 2 (with phase-space R 4 ) may be represented using a Wiener integral along a small time increment δt of the form [16] L t which is sometimes referred to as the Fokker-Planck operator [13]. We note that the proper definition ofẊ (τ ) is a delicate issue that we avoid here by using a reformulation of the model for discrete flow maps, described later. The classical method for computing phase-space densities transported via deterministic transfer operators dates back to the early sixties and is known as Ulam's method [17]. In this approach the phase-space is subdivided into cells, and local approximations of the transition rates between these cells are applied. More recently, modified Ulam-type methods have been devised for the case of stochastic transfer operators, see for example Refs. [15,18]. Alternatively, periodic orbit techniques have been developed for stochastic transfer operators since the late nineties [12,[19][20][21][22]. Initial studies focused on determining spectral properties of the Fokker-Planck operator for Langevin flows in the weak noise limit. More recent work has considered higher dimensional cases [13] and the estimation of stationary distributions [23]. On the other hand, one may also consider the modeling of uncertainties in stochastic dynamical systems propagated via a deterministic transfer operator [24]. The equivalence of such an approach to stochastic transfer operator based models is well-known, see for example Ref. [21].
In this work we consider the case where V describes a Hamiltonian vector field, and hence the divergence operator in the Fokker-Planck equation (2) reduces to the Poisson bracket ∇ · (Vρ) = {ρ, H}, where H is the associated Hamiltonian. Ray tracing can be interpreted in the sense of modeling rays as trajectories of particles following Hamiltonian dynamics, and so this special case includes short wavelength asymptotic approximations to high-frequency waves modeled via the evolution of a probability density (or wave energy density if ∫ ρ(Y , t)dY is not necessarily equal to unity) along ray trajectories. We will focus on stationary solutions of the Fokker-Planck equation satisfying where the corresponding wave problem is now formulated in the frequency-domain, instead of the time-domain [7]. Eq. (4) will therefore form our underlying basic model for the long-time stochastic evolution of densities along ray trajectories through two-dimensional domains. We will also consider the generalization to the case when σ is replaced by a vector σ, which specifies differing levels of diffusion in position and momentum space. The aim of this work is to uncover how a source probability density evolves through uncertain two-dimensional domains under the action of the stochastic evolution operator, and thus to quantify uncertainties in solutions of Eq. (4). In order to track the density through reflections at the boundary of the domain, we reformulate the Fokker-Planck operator (3) to evolve a discrete map ϕ, instead of a continuous time trajectory evolution X (t). The map ϕ corresponds to discrete values of X (t) at times when the trajectory is positioned on the boundary of the domain Ω ⊂ R 2 . Note that this reformulation leads to a reduction in dimensionality; the boundary phase-space coordinates in position and momentum are both one dimension lower than before. The trajectory evolution X (t) at the discrete set of boundary intersection times may be written a stochastic boundary map ϕ σ (X ) = ϕ(X) + X ε , where X ε is a vector of independent random variables. In particular, the probability density function (PDF) f σ from which these random variables are drawn appears as the kernel of our modified stochastic transfer operator [14]. In contrast with the Gaussian PDF appearing in (3), the support of the PDFs f σ considered here are, in general, bounded. In the simplest case the PDF could describe a uniformly distributed probability of location and tangential momentum, which leads to a non-parametric SEA-type approach as described above. In this study we will instead base our model on the truncated normal distribution, which includes the uniform distribution model as a limiting case. We will focus, in particular, on the weak noise case with near deterministic dynamics when the truncated Gaussian PDF approximates a Dirac delta distribution. Finally, we consider the influence of damping on the dynamics, and ultimately the level of noise present in the stationary density ρ.

A boundary integral description of stochastic propagation
Let Ω denote a finite two-dimensional domain with an associated speed of propagation c. Consider phase-space coordinates (r, p) ∈ Ω × R 2 , where r is the spatial position and p is the momentum or slowness vector. When c is a constant, the HamiltonianĤ = c|p| = 1 describes trajectories within Ω moving in straight-line paths between reflections due to hitting the boundary Γ = ∂Ω. We write the phase-space coordinates on Γ as X = (s, p), where s is an arc-length parameter for Γ and p = c −1 sin(θ ) is the tangential component of the momentum vector p at the point s. Here θ ∈ (−π /2, π/2) is the angle between the trajectory leaving the boundary at s and the normal vector to Γ at s. We will restrict to convex polygonal domains Ω to avoid additional complications due to incorporating visibility functions and curved edges.
The stochastic propagation of a phase-space density ρ between intersections with the boundary Γ is described by an operator of the form [14] L σ ρ(X) = Here, Q = Γ × (−c −1 , c −1 ) denotes the phase-space on the boundary and ϕ : Q → Q defines the boundary flow map mentioned earlier, which maps a vector in Q to another vector in Q . The boundary map ϕ(X ′ ) describes a deterministic next intersects Γ at position ϕ s (X ′ ) and the momentum ϕ p (X ′ ) corresponds to a specular reflection at ϕ s (X ′ ). The kernel of the boundary integral operator (5) is given by a probability density function (PDF) f σ such that where σ = (σ s , σ p ) is a parameter set controlling the shape of f σ as a function of position s and momentum p. The physical meaning of the term f σ (X − ϕ(X ′ )) is that it gives the probability of a particular arrival position and tangential momentum X relative to the deterministic boundary map ϕ. For example, if f σ is a PDF with mean zero, then the mean arrival position and tangential momentum X is ϕ(X ′ ) and so, on average, trajectories will travel across Ω in a straight-line path and undergo specular reflections at the boundary. On the other hand, if f σ is a uniform distribution then the trajectories are equally likely to arrive at any position on the boundary within the support of the PDF and undergo Lambertian reflection when they reach Γ .
We choose f σ to be an uncorrelated bivariate truncated normal distribution as detailed in [14] and sketched below. The truncation limits arise from the interpretation of the evolution given by the operator in Eq. (5) as originating from a stochastic boundary map ϕ σ (X ′ ) = ϕ(X ′ ) + X ε with added noise X ε = (s ε , p ε ). For X ∈ Q given, we have to ensure that ϕ(X ′ ) = X − X ε is still in the range of the deterministic map ϕ; this gives rise to restrictions on the possible values of X ε and thus on the domain of f σ . Denoting these truncated ranges by (X − , X + ) where X ± = (s ± , p ± ), the PDF f σ will have support on X ε ∈ (s − , s + ) × (p − , p + ) only, as depicted in Fig. 1 for the spatial limits s ± . In the momentum coordinate we have that . A cut-off function of the form χ (X−,X+) (X ε ) = χ (s−,s+) (s ε )χ (p−,p+) (p ε ) for restricting the support of f σ to (X − , X + ) can be defined using step functions as detailed in Ref. [14]. An uncorrelated bivariate truncated normal distribution may then be expressed as where the normalization defined through ψ σs and ψ σp is given as and ψ σp is defined analogously. The scaling ensures that the PDF satisfies condition (6) for the truncated sampling ranges specified through χ (X−,X+) . The mean and variance of f σ differ, in general, from that of the underlying normal distribution, that is, from 0 = (0, 0) and σ = (σ s , σ p ), respectively. As σ s → 0 and σ p → 0, then one obtains a deterministic propagation model as f σ becomes an increasingly peaked Gaussian, which tends to a two-dimensional Dirac delta in the limiting case.
That is, the ray will follow a straight-line path prescribed by the initial vector X ′ when σ s → 0 and will reflect specularly upon reaching Γ when σ p → 0. On the other hand, in the limit as σ s → ∞ and σ p → ∞ then f σ becomes a uniform distribution [14]. That is, for σ s → ∞ the ray will travel to all points on Γ \ E(s ′ ) with equal probability, where E(s ′ ) is the polygonal edge containing s ′ as shown in Fig. 1. The limit σ p → ∞ corresponds to the ray undergoing Lambertian reflection upon reaching Γ . We note that the connection to the Lambertian reflection law can be explicitly derived by noting that a uniform distribution over p means that each possible p ∈ (−c −1 , c −1 ) must have probability c/2. If we change the variable from p to θ using p = c −1 sin θ, then we obtain which is Lambert's (cosine law) distribution cos (θ) /2 for the probability of each reflected direction θ. The transition from σ p = 0 to σ p → ∞ therefore corresponds to the classical model in room acoustics for interpolating from pure specular reflection to Lambertian reflection. The mean and variance of the truncated normal distribution may be calculated from the PDF (7) using the standard formulae. The variance of the bivariate distribution will tend to σ in the weak noise limit as σ → 0. In the next section, we begin by deriving a model for quantifying the uncertainties propagated by the operator (5) in this weak noise limit.

Uncertainty quantification
The stochastic propagation model introduced above may be used to describe a number of types of parametric uncertainty arising in a high-frequency wave model, including uncertain sources, geometry and material parameters. An additional paper on this topic is currently being prepared [25], where it is shown that in the weak noise asymptotic regime one can derive direct connections between the variance for a normally distributed and linearly varying noise term in the propagation speed c and the parameters σ s and σ p considered here. A similar analysis is also performed for an uncertain boundary location. Some introductory work in this direction is reported in Ref. [26]. However, in this work we instead study the effect of the stochastic propagation model on the statistical properties of a prescribed initial probability density ρ 0 . We make use of the fact that in the weak noise regime, the noisy trajectory map may be linearized about the corresponding deterministic trajectories, using the noise as the expansion variable [13,22]. An alternative approach for the strong noise case is to instead analyze non-parametric uncertainties in the sense of SEA. In this case a uniformly distributed density will be mapped to another uniformly distributed density and the variance depends only on the sizes of the domain and range of the noisy flow map ϕ σ .

Weak noise approximation of the variance evolution
Consider the propagation of a density from the phase-space boundary coordinate X ′ =X ′ + X ′ ε , via the noisy flow map ϕ σ , to the phase-space boundary coordinate X =X + X ε . Here the over-bar notation relates to the deterministic (or mean) trajectory coordinate and the ε subscript is used to indicate the noise term in each case. Let the truncation limits for X ε be denoted X ± as before, and for X ′ ε we correspondingly denote the truncation limits by X ′ ± . The deterministic trajectory evolution is then described by ϕ(X ′ ) =X, and taking the Taylor expansion of ϕ(X ′ ) aboutX ′ leads to ε , which is valid in the weak noise limit. The stochastic transfer operator (5) may therefore be approximated as the convolution Note that for a flow between two prescribed edges of a polygonal domain Ω, the deterministic trajectory flow has the property that the final direction is independent of the initial position, and hence we have J 21 = 0 in the matrix J above.
In this work we consider boundary source terms of the form for s belonging to specified edges of a polygonal domain Ω of cumulative length A, and ρ 0 (s, p) = 0 otherwise. This type of source term was introduced in Ref. [14] and for small σ 0 , it corresponds to a unit boundary density propagating (on average) with momentum p = 0, perpendicular to the boundary (the extension to other mean propagation directions is straightforward). For large σ 0 it corresponds to random propagation from the specified part of the boundary. Applying the weak noise stochastic transfer operator (10) to ρ 0 given by (11) in the weak noise limit leads to a convolution of two Gaussians in the momentum coordinate as follows: Here we have used that J 21 = 0 and the scaling factors tend to 1 and the integration limits can be extended to ±∞ in the weak noise limit; the rapid decay of the Gaussian away from the deterministic trajectory means that the extended regions of the integral are essentially zero. The integration with respect to s ′ ε can be evaluated in terms of the error function and on substitution of the limits simplifies to |J 11 | −1 . Hence which is a Gaussian in the momentum variable with variance term J 2 We observe that for the choice of initial density ρ 0 given by (11), the variance term σ 2 s for the stochastic propagation of the position variable does not enter the above expression for L σ ρ 0 in the weak noise regime. However, the above analysis can in principle be generalized to (truncated) Gaussian boundary source terms in both the position and momentum variables, and even to sources and propagation operators including correlations-see Ref. [13]. In the most general case, a bivariate Gaussian source distribution with covariance matrix Σ 0 is mapped, under the action of a weak noise stochastic transfer operator with a bivariate Gaussian kernel and covariance matrix Σ, to a bivariate normal distribution with covariance JΣ 0 J T + Σ [13]. We have therefore shown how the variance of the initial probability density ρ 0 evolves as it is propagated by the operator L σ in the weak noise regime. After applying sufficiently many iterates of L σ , the variance will converge only if J 2 22 < 1, in which case the limiting value is σ 2 p /(1 − J 2 22 ) and is independent of σ 2 0 . Note that if J 2 22 is close to, but still less than 1, then the weak noise assumption will eventually be violated. This is also the case if J 2 where the weak noise theory would predict incorrectly that the variance tends to infinity (linearly or exponentially, respectively) as the number of iterations of L σ is increased. For the polygonal domains considered here we have J 22 = − cos(θ )/cos(θ ′ ) wherē θ ′ ∈ (−π /2, π /2) is the angle between the trajectory specified byX ′ = (s ′ ,p ′ ) and the normal vector to Γ ats ′ , and analogously forθ . For a general polygonal domain Ω, the situation can become complex with some regions of phase-space where the variance is converging, and other regions where it is diverging as the number of iterations of L σ is increased. In fact, for every trajectory with a converging variance there exists the same trajectory, which if traversed in the opposite direction, will diverge. Whether a particular trajectory actually enters the model in the weak noise regime will depend on both the source term and the geometry of Ω. The predictions for the variance of the propagated density detailed in this section will be tested numerically in Section 4.2.

Long-time behavior and the stationary density
The long-time behavior of the stochastic propagation model will be studied by considering the sum where L j σ denotes the jth iterate of L σ . Taking the limit as n → ∞ in (14) leads to a stationary density ρ = ρ ∞ as described by Eq. (4). However, this sum only converges if damping is present, in which case the stationary density may be computed using the standard Neumann series result and hence by solving the second-kind integral equation Note that here L σ ρ 0 no longer corresponds to a probability density since the damped version of the operator L σ no longer conserves probability/energy. For the case when the weak noise assumption is only violated in L j σ ρ 0 when j is sufficiently large, then the level of noise contributing to the stationary density from later iterates could be controlled by damping. We will apply damping in the model via a multiplicative factor of the form exp(−µd(s, s ′ )) inside the integral defining L σ (5), where d(s, s ′ ) denotes the Euclidean distance between the boundary points s ′ and s and µ ⩾ 0 is the attenuation per unit (ray) length. In Section 4.2, we will investigate numerically the extent to which a given damping parameter µ will diminish the contributions to the stationary density ρ from the later iterates of L σ , relative to the early iterates where the weak noise approximation may be valid.

Implementation and results
In this section we describe the discretization of the operator L σ and then calculate the statistical properties of a truncated Gaussian probability density (11) after it has been transported under the action of L σ .

Discretization
The operator L σ is bounded and compact for σ s , σ p > 0 [27], unlike the corresponding deterministic operator (Frobenius-Perron operator) studied, for example, in Ref. [7] and obtained in the limit σ s , σ p → 0. As a consequence we may apply any standard method for the discretization of the integral operator L σ , such as the Nyström, collocation or Galerkin methods. In contrast, for a direct discretization of the Frobenius-Perron operator one is restricted to methods that regularize the operator by imposing it in a weak form and so typically Galerkin methods are employed [7].
The discretization of the operator L σ is performed using piecewise collocation in the position variable s and the Nyström method in the momentum variable p as described in Ref. [14]. However, a notable improvement on the methods presented in [14] is achieved by using appropriately sub-divided Clenshaw-Curtis quadrature for the Nyström method, which leads to spectral convergence of the approximation scheme in p as reported in Ref. [28]. The choice of the piecewise constant collocation method in space has the major advantage that the boundary integrals required for the evaluation of L σ ρ may be carried out analytically. Note also that for the polygonal domains considered in this work one cannot, in general, expect a high degree of regularity in the spatial dependence of the density ρ.
The first step towards implementing the discretization scheme described above is to divide Γ into n boundary elements. We then separate the spatial and momentum variable dependence of the density ρ and project onto a piecewise constant spatial basis as follows: Here ρ j are a set of directionally dependent expansion functions to be determined and b j are a set of piecewise constant spatial basis functions; i.e. b j (s) = 1 if s lies on the jth boundary element and zero elsewhere. Substituting (16) into (5) we where e j = [α j , β j ] ⊂ Γ denotes the jth boundary element. Note that the (spatial) integral with respect to s ′ appearing in (17) can be solved analytically in terms of the error function erf; this remains tractable if we include an additional damping factor of the form exp(−µd(ϕ s (X ′ ), s ′ )). In what follows, we denote this spatial integral as S j µ , where for the undamped case when µ = 0. Note that the integration limits lie at the end points of (or inside) the jth boundary element e j ; if the limit is inside the element then it corresponds to the pre-image of a vertex of the domain for rays with direction p ′ .
We now change the variable in the momentum integration from p = c −1 sin(θ ) to the direction angle θ and define collocation points s i , i = 1, . . . , n, where s i is chosen as the centroid of the ith boundary element. Applying a numerical integration rule with nodes p k and weights w k , for k = 0, . . . , K , then the combined collocation and Nyström method discretization of L σ ρ is given by where i = 1, . . . , n, κ = 0, . . . , K and f σp is given by (7). As mentioned earlier, the quadrature rules employed for the simulations in the next section are based on the Clenshaw-Curtis rule applied to a number of sub-divisions of the direction integral. The subdivisions are performed at momenta p ′ where the spatial integral S j µ (s, p ′ ) can become non-smooth.
Explicitly, these momenta correspond to directions where the ray arrives at a corner of the polygonal domain Ω, or at one of the collocation points. Note therefore that the number of sub-divisions in the Nyström method quadrature rule both grows as the spatial mesh is refined, and depends on which element j contains the starting point of the trajectory. For the numerical experiments reported in the next section, we simplify the implementation by fixing the number of subdivisions for all elements to be the maximum required for any particular j = 1, . . . , n.

Numerical results
In this section we evaluate L σ ρ 0 for the case when Ω is a rectangular domain and estimate its variance in both the strong and weak noise regimes. We then move on to consider damped problems and properties of the stationary density. The reason for choosing this relatively simple domain is that it is straightforward to keep track of the deterministic trajectory paths. This property helps when interpreting the results in the weak noise case, particularly after multiple iterations of L σ . We choose the rectangular domain to be Ω = {r = (r 1 , r 2 ) ∈ R 2 such that 0 < r 1 < 0.75, 0 < r 2 < 0.25} and impose a non-zero source boundary density ρ 0 given by (11) along the left-hand edge at r 1 = 0 only. We choose the propagation speed as c = 1 meaning that both the position and momentum variables have the same range. Hence the values of the parameters σ s and σ p , governing the level of noise in the position and momentum of the flow map ϕ σ , are directly comparable. The predictions in Section 3 for the limiting cases of weak and strong noise will be compared against direct calculations on the computed densities ρ(s, p) using the standard formulae for the meanρ(s) and variance v(s) given byρ (s) = ∫ 1 −1 p ρ(s, p)dp ∫ 1 −1 ρ(s, p)dp respectively. We note that the normalization factors (denominators) have been introduced to generalize to the case when ρ(s, p) is not necessarily a probability density in the momentum variable as a function of s.
For the particular example of the rectangle described above, the underlying deterministic propagation model corresponds simply to rays leaving the left edge perpendicular to the boundary, then arriving at the right edge and reflecting through π rad, and so on. Before proceeding to the computational results for this example, we summarize the corresponding main results from Section 3 (with no damping µ = 0) for clarity as follows: • In the weak noise case σ 0 , σ s , σ p ≪ 1, an initial density ρ 0 of the form (11) on the left-hand edge with approximate variance σ 2 0 will map (after one reflection) to a density L σ ρ 0 on the right-hand edge with approximate variance v(s) = σ 2 0 + σ 2 p , since cos(θ ) = cos(θ ′ ) = 1 and so J 2 22 = 1, |J 11 | = 1 -see Eq. (13). • The weak noise theory predicts that the limiting value of the variance after infinitely many reflections will diverge, since the explicit expression for the limit is σ 2 p /(1 − J 2 22 ) and J 2 22 = 1. We therefore expect that the weak noise predictions will become invalid after a finite number of reflections.
• In the strong noise case with large σ p > 1, the variance of L σ ρ 0 (after one reflection) on the right-hand edge will be v(s) = 1/3. As the reflection order increases, the density will spread over all edges with (normalized) variance v(s) = 1/3.

Statistical properties of L σ ρ 0
Initially we consider the undamped case (µ = 0) and analyze the effect of different choices of the parameters σ 0 , σ s and σ p on the density L σ ρ 0 . The left sub-plot of Fig. 2 shows the computed boundary phase-space density ρ 1 = ρ 0 + L σ ρ 0 in the weak noise case with σ 2 p = σ 2 0 = 1e−4 and σ s = 0. (Note that the exact spatial integral S j µ , j = 1, . . . , n defined in the previous section still converges in the case σ s = 0.) For this plot we have used Eq. (19) to compute L σ ρ 0 with a total of 1025 quadrature points for the Nyström method in the momentum variable (that is, K = 1024) and n = 88 spatial boundary elements. The right sub-plot of Fig. 2 shows the numbering of the spatial collocation points around the rectangular boundary Γ . The space axis in Fig. 2 has been divided by a white line separating ρ 0 on the left side of the rectangle Γ (elements 78 to 88), and L σ ρ 0 on the other three sides of Γ (elements 1 to 77). Elements 34 to 44 correspond to the right hand side of the rectangle Γ as shown in the right sub-plot of Fig. 2, and in the weak noise regime it is essentially only these elements on which the density L σ ρ 0 has support. Note that the sharp change of solution behavior between different edges is a consequence of the parameter choice σ s = 0. The weak noise behavior of the model is therefore as expected; a Gaussian with a sharp peak at the p = 0 (perpendicular to the boundary) on the left-hand edge of Γ is transported by L σ to another Gaussian on the opposite edge with a slightly broadened distribution, again peaking at p = 0. We note that physically, the behavior shown in Fig. 2 is close to that of a specular reflection at the right-hand edge of the domain; however, a pure specular reflection would appear as an infinite Dirac delta peak at p = 0 from collocation points 34 to 44.
In Fig. 3 we investigate the effect of the parameters σ p , σ s and σ 0 on the variance of L σ ρ 0 , computed using (20). Each of the three rows of the plot correspond to one of the parameters being varied whilst the other two are fixed sufficiently small for the weak noise theory to be applicable. The default values for the fixed parameters are those from Fig. 2: σ 2 s = 0 and σ 2 p = σ 2 0 = 1e−4. The left column shows the variance computed at the collocation point (0.75, 0.125), which is located at the center of the right hand edge, compared to the theoretical prediction for the weak noise regime (v(s) = σ 2 0 + σ 2 p ) described above. The right column shows the plots of the boundary phase-space density ρ 1 = ρ 0 + L σ ρ 0 when the parameter being varied in that particular row is set to its largest value. These plots provide a direct comparison with the weak noise case shown in Fig. 2.
In the first row of Fig. 3 we observe that as σ 2 p is increased, the computed variance follows the theory from the weak noise regime reasonably closely until almost the point where the weak noise theory plot intersects the plot of the variance for the corresponding uniform distribution, which has value 1/3 as described above. For large σ 2 p , the computed variance then corresponds to this uniform distribution. These observations are consistent with the right sub-plot, which shows that the Gaussian source ρ 0 (above the white line) has been transported by L σ to a uniform distribution on the opposite edge (elements 34 to 44) with height 2 as expected (recall that the momentum coordinate range is 2). Note that the color axis has been truncated to more clearly show L σ ρ 0 ; the maximum value of ρ 0 is still approximately 160 as shown in Fig. 2. The physical interpretation of the right sub-plot is a Lambertian reflection (9) along the right hand edge of the domain since all momenta p are (approximately) equally probable for collocation points 34 to 44. Furthermore, although the computations in the left plot are performed for just one point in space, the computed variance is independent of the point chosen on the right hand edge due to the highly directive nature of ρ 0 ; the points on the right edge of Γ receive a highly dominant contribution to the ray density from the corresponding point on the left edge with the same r 2 coordinate. For this same reason, the computed density L σ ρ 0 is a probability density (only) at the collocation points along the right hand edge of Γ .
The second row of Fig. 3 shows that the weak noise theory is not accurate for σ 2 s ⩾ 0.01. In this case, since σ 0 and σ p are fixed, the weak noise theory predicts a constant variance of σ 2 0 + σ 2 p = 2e−4. When σ s is increased, the computed variance eventually reaches a steady value of around 0.015 where the density is spread uniformly across all spatial collocation points on the three receiving edges of Γ , as shown below the horizontal white line in the right sub-plot. The distribution L σ ρ 0 remains localized in momentum with a lower peak and larger variance than the initial density ρ 0 . As before, the color axis has been truncated to more clearly show L σ ρ 0 . Owing to the even spreading of the initial density ρ 0 , which has support on an edge of length 0.25, to the computed density L σ ρ 0 , which has support on edges of cumulative length 1.75, then L σ ρ 0 is Motion ( not a probability density on the receiver edges in the strong noise limit and instead integrates to 1/7 (the ratio of the initial edge length to the cumulative receiver edge lengths).
In the third row of Fig. 3 we observe that the weak noise theory is accurate for σ 2 0 ⩽ 1e−4. As σ 2 0 is increased, there is a transition region until the uniform distribution estimate for the large noise limit is valid. This transition region appears to be centered approximately at the point where the weak noise theory plot intersects the plot of the variance for the corresponding uniform distribution. The variance of the uniform distribution here is smaller than before since the range of propagation angles at (0.75, 0.125) has been reduced as shown in the right plot at collocation point number 39. This reduction is because the source direction is approximately random, but the reflection at the right-hand edge is approximately specular and so the range of propagation angles relates to the range of admissible incoming angles at (0.75, 0.125). The range can be approximated using basic geometry (see Fig. 4) since for a uniformly distributed source on the left edge of Γ , the rays that can reach the point (0.75, 0.125) at the center of the right edge will arrive at angles within the range (− tan −1 (1/6), tan −1 (1/6)). The reflection angles will therefore also lie in the same range due to symmetry, and because σ 2 p is small, this will also approximately be the range for the propagation angles from (0.75, 0.125). The variance for a uniform distribution in the momentum variable is therefore sin 2 (tan −1 (1/6))/3 = 111 −1 . Again, the normalization introduced in (20) plays a role; the right sub-plot shows that the uniform distribution has height 2 (the same height as ρ 0 ) and so it will not integrate to one with respect to momentum since the range of directions has been reduced from (−π /2, π/2). The range of the uniform distribution is therefore geometry dependent and so both the mean and variance depend on the position on Γ , as shown in Fig. 5. The mean and variance values shown in this plot can also be verified, as above, using basic geometric arguments to determine the range of uniformly distributed tangential momentum p from a given collocation point on Γ .

Higher order reflections and the stationary density
We now consider how the variance of the initial density (11) evolves under multiple iterations of the operator L σ . We will then introduce damping and explore links to the variance of the stationary density, again computed via the formulae (20). Throughout this section we fix σ 2 s = 0 and σ 2 0 = 1e−4. onward. The large values on these edges in the weak noise case can be explained since the most probable trajectories to reach these edges do so at directions which are close to tangential to the edge; the underlying deterministic propagation in the model is between the left and right edges in the direction perpendicular to these edges. Hence for the upper and lower edges of the rectangle, the associated density in the weak noise limit will be bi-modal at the limiting tangential directions ±π/2 rad, but with mean 0. The decay of the computed variance (on the upper and lower edges) for later iterates therefore actually corresponds to a build up of noise in the model.  In the other extreme case when σ 2 p = 100, we notice that the computed variance reaches the steady value of 1/3 after only one iteration and takes this value for all collocation points and all future iterates. Recall that the value of 1/3 corresponds to the variance of the uniform distribution on the momentum coordinate range (−1, 1). In the upper-left corner of the lowerright sub-plot one can notice the initial density variance of 1e−4 for iteration 0 and collocation points 78 to 88 along the left hand edge of the rectangle. For the intermediate cases one notices that the steady variance is becoming more uniform (and closer to the value 1/3 overall) as σ 2 p is increased. When σ 2 p = 1e−2, the noise builds up quickly enough to violate the weak noise theory after only a few iterates and the steady state is achieved from around iteration j = 50 onwards. For σ 2 p = 1, the steady equilibrium is reached after fewer than 10 iterations of L σ , and closely resembles the uniform distribution observed when σ 2 p = 100. Note that the transition between different steady equilibria observed here is a consequence of the simple integrable dynamics of the underlying deterministic model. For irregular domains we will instead observe ray chaos [4,29] and the steady equilibrium will always be a uniform distribution with variance 1/3, even for a deterministic propagation model.
We now consider the properties of the stationary density ρ, which we recall is given by the limit as n → ∞ of ρ n = ∑ n j=0 L j σ ρ 0 . Note that the computed variance of L n σ ρ 0 closely resembles the computed variance of the sum ρ n (14) provided that n is taken sufficiently large; the repeated addition of the steady equilibrium distribution eventually leads to a distribution ρ n with approximately the same shape as the equilibrium distribution L n σ ρ 0 . The results in Fig. 6 therefore suggest that noise will always build up in the model and hence in the stationary density ρ, leading to a distribution very different from those described by the weak noise theory outlined in Section 3. However, this may not necessarily be the case when we introduce damping into the model, which is necessary for computing the stationary density ρ. We now consider the influence of damping on the properties of the stationary density for the case σ 2 p = σ 2 0 = 1e−4 and σ 2 s = 0. In particular, we investigate whether the influence of the later and more noisy iterates in the sum (14) can be controlled.
The upper sub-plot of Fig. 7 shows the computed variance of the stationary density ρ obtained by solving Eq. (15) for different values of the viscous damping parameter µ. This plot should be read in conjunction with the lower plot, which shows the total density present at each collocation point, and therefore which parts of the boundary Γ contain the dominant D.J. Chappell, G. Tanner / Wave Motion ( ) - contributions to the ray density ρ. The variance of the undamped cumulative density ρ 1500 is also shown on the upper subplot for reference; note that it bears a strong similarity the profile shown in the upper-left plot of Fig. 6 once the steady equilibrium has been reached. Fig. 7 shows that for smaller values of the damping parameter, the density ρ is concentrated along the left and right edges of Γ and that as the damping is increased to µ = 6.4 or 12.8, only the source density ρ 0 remains dominant with total density 1 and variance 1e−4 as expected. The computed variances along the left and right edges of Γ for µ ⩾ 0.4 are around 2e−3 or lower and start to enter the range where the weak noise theory can be applied (see Fig. 3). This result is also demonstrated qualitatively in Fig. 8, which shows the stationary density for µ = 0.4. The left sub-plot shows a similar density distribution to the weak noise case shown in Fig. 2, and is dominated by peaks at 0 rad on the left and right edges of Γ . The right sub-plot shows the same result but with the color-axis maximum truncated at 2. This allows one to observe the behavior of the density on the upper and lower edges of Γ and explains the relatively large variances computed here. In particular, one has a bi-modal distribution with peaks at momenta ±1, and as discussed earlier, this is the expected behavior in the weak noise case.
Finally, we apply the weak noise theory from Section 3 to provide approximate upper bounds for the variance of the stationary density. We note that in the weak noise regime and for a sufficiently large choice of damping parameter µ, the stationary density along the left and right edges of Γ is a weighted superposition of Gaussians of increasing width and decaying height. This is a consequence of the fact that the trajectories for the corresponding deterministic problem all have equal length, meaning that the damping term would merely apply a constant decay factor that does not alter the shape of the distribution. The computed variance will therefore approximately be bounded above by the variance of the broadest Gaussian contributing to the superposition before damping renders such contributions insignificant. For example, consider the variance along the right-hand edge of Γ with the damping parameter µ = 3.2. Here, only odd numbered iterates j of L j σ ρ 0 will contribute to the superposition and the damping will lead to a decay factor of approximately exp(−1.5µ) for the iterate j + 2 relative to the iterate j, as well as a decay factor of exp(−0.75µ) for the first iterate relative to the initial density. Hence, as a result of damping, the contribution to the stationary density will decrease by in excess of three orders of magnitude relative to the initial density after only j = 3 iterates. The weak noise theory then tells us (in the undamped case) that L 3 σ ρ 0 is a Gaussian with variance σ 2 0 + 3σ 2 p = 4e−4, which can be observed as approximately the maximum variance of the stationary density for µ = 3.2 along the right-hand edge shown in Fig. 7 damping factors, but it is important that the damping leads to sufficient decay of the contributions to the superposition before the noise builds up sufficiently to violate the weak noise theory.

Discussion: extension to more complex geometries
The results presented in the previous two subsections for ray propagation in a simple rectangular domain provide a proof of concept for the uncertainty quantification theory discussed in Section 3. Whilst this example is useful for understanding how the uncertainty changes as it propagates in a simple intuitive situation, it does naturally lead to the question of how widely these ideas can be applied in general. The theory in Section 3 can, in principle, be directly applied to any convex polygon as considered in Section 2 and so no additional restrictions are added on the range of applicability at this stage. However, the weak noise theory requires knowledge of the mean incoming and outgoing ray directions at each boundary position where one wishes to estimate the variance. Since these mean directions correspond to the underlying deterministic boundary map ϕ, then any practical application of the weak noise theory described in Section 3 relies on a full knowledge of the ray trajectories ϕ at each reflection. The method therefore has a similar limitation to traditional ray tracing and will D.J. Chappell, G. Tanner / Wave Motion ( ) - become unfeasible for high reflection orders in complex geometries, when the number of distinct ray paths to trace becomes intractably large. However, as discussed in the previous section, one typically has a build up of noise in the system as the reflection order increases and hence the range of practical applicability of the weak noise theory also coincides with the typical validity range of the theory in the first place. Hence, the wider application of the weak noise theory should be feasible provided that the focus is on similar cases to those studied here; that is the properties of the density ρ 1 after one reflection as considered in Section 4.2.1, or the properties of the stationary density ρ in the presence of strong enough damping to ensure that only the early reflections dominate its behavior.

Conclusions
An integral equation model for the stochastic propagation of densities through phase-space in two-dimensional domains has been presented. The integral kernel is chosen as an uncorrelated bivariate truncated normal distribution, which has the property of interpolating between deterministic ray tracing and random propagation according to the level of noise prescribed in the model. We have shown that in the weak noise limit the variance of an initial Gaussian probability density will be propagated to another Gaussian, and derived an explicit expression for its variance. In the strong noise limit we have shown how the variance corresponds to that of a uniform distribution, and have demonstrated the applicability of both of these theories for the limiting cases numerically. Finally, we have studied the properties of the stationary density in the limit of infinitely many reflections. We found that the weak noise theory could provide upper bounds on the variance of the stationary density, provided that the damping parameter was taken to be large enough to effectively remove the noise that builds up in the system after many reflections.