Diffraction of Wigner functions

We describe the contribution of diffractive orbits to semiclassical approximations of Wigner function propagators. These contributions are based on diffractively scattered rays used in the geometrical theory of diffraction (GTD). They provide an extension of well-established approximations of Wigner-function propagators based on rays that propagate by specular reflection and refraction. The wider aim of this approach is to allow for diffractive mechanisms to be accounted for in Eulerian approaches to ray-tracing simulations. Such approaches propagate densities of rays rather than follow rays individually. They promise to be a more efficient means of performing ray-tracing simulations in complex environments with applications in, for example, planning of wireless signal coverage for mobile communication networks.


Introduction
Accounting for diffraction is a vital part of planning radio coverage in large-scale, complex environments. Scattering of radio waves from building edges, for example, allows radio signals to enter shadow regions of purely geometrical propagation and is central to the realisation of usable wireless reception in obstacle-rich, urban environments and in building interiors. Modelling of radio coverage is best tackled using ray tracing methods [1]: the aim of this paper is to allow diffractive mechanisms to be more efficiently incorporated into such phase-space simulations by calculating diffractive contributions to propagators of Wigner functions.
Diffractive mechanisms are routinely incorporated into ray-tracing simulations of radio coverage [2] using developments of the geometrical theory of diffraction (GTD) [3] and the uniform theory of diffraction (UTD) [4]. Their contribution to electromagnetic wave paths is important to key aspects of the performance of digital communications, such as path loss [5] and power and delay spread [6]. However, in typical implementations, these require searches to find orbits beginning and ending at particular locations and their rapid proliferation with length restricts calculations in practice to a handful of reflections. This is essentially a Lagrangian approach and implies computationally intensive methods. The efficiency of such approaches can be increased by using parallel implementations of ray tracing [7,8] and improved modelling of shadow loss from non-planar geometries [9].
In this paper we provide a means of incorporating diffractive rays into an Eulerian approach, that is, an approach where the propagation of power is modelled in terms of a phase-space density rather than by following individual rays. This can lead to efficiency gains in multiplescattering scenarios, such as modelling indoor wireless coverage where losses are typically low and multiple reflections dominate. Numerically implemented Eulerian approaches based on the propagation of phase-space densities, such as the dynamical energy analysis (DEA) method [10], offer a promising means of extending ray-tracing simulations to complex, multi-reflection environments. The ability to include diffractive mechanisms in such approaches is vital in many potential applications, such as for modern mobile communication networks [11]. The actual implementation is, however, far from obvious as it requires the inclusion of wavelength-sized boundary effects in an Eulerian transport equation on phase space.
The direct wave analogue of a phase-space ray density is a field-field correlation function or, in quantum-mechanical language, a density operator. A direct quantitative comparison is provided by transforming such correlation functions to Wigner functions [12], which are pseudo-densities on phase space that may be compared directly to phase-space densities calculated 'classically' or by ray tracing. Our interest here is in approximating a propagator for Wigner functions, which relates the Wigner function characterising the field-field correlation function in the receiver region to a Wigner function characterising the source (see, for example, [13][14][15][16]). Full electromagnetic wave calculations in multiple scattering environments typically show strong fluctuations due to the interference between multiple paths of propagation, which may be approximated by summing over ray pairs satisfying appropriate boundary conditions [13,17]. However, implementing such schemes in full in a complex environment is a computationally challenging task. In extending such calculations to account for diffractive orbits, we are therefore content in the first instance to focus on averaged quantities, in which multi-path interference effects are suppressed. In classical ray simulations, such averaged densities may be approximated by direct ray simulations rather than by pairs of rays.
Even in the context of predicting averaged correlation functions, it is important to account for interference between diffractive orbits and nearby specular orbits in order to satisfy global symmetries, such as flux conservation in dissipationless systems. The leading effect of including diffractive orbits in a propagation of phase space densities is to add density localised along a manifold of orbits originating on the diffractive structure (such as an aperture, an edge or a corner). This added density is positive. In order not to violate energy conservation, there must be a compensating loss elsewhere in phase space. This occurs along forward scattering directions where diffractive and specular orbits merge, and treatment of this loss is accounted for by interference between diffractive and specular orbits. We show that overall flux conservation is guaranteed, where appropriate, by an optical theorem applied to diffractive scattering.
In many diffraction problems, the GTD approximation for diffractively scattered waves breaks down in the forward direction, and must be replaced there by the UTD [4]. The previously mentioned analysis of overall flux conservation must therefore be performed in the context of UTD in the most general setting. In this paper we consider the special case of diffraction problems where a simpler GTD approximation remains valid in the forward scattering direction. Examples where this is the case include scattering by a small obstacle or diffraction through a small aperture [18]. This allows us to establish the most important features of diffracted phase space densities in the simplest setting. The most general case where UTD is required in the forward scattering direction will be treated elsewhere.
The paper is organised as follows. In section 2 we establish the background and notation assumed in the rest of the paper. In particular we summarise the boundary representation used to write the wave solutions and their phase-space counterparts. We also describe how GTD is accounted for in this notation. In section 3 we describe how the optical theorem is used to establish overall flux conservation, where appropriate. In section 4 we present the detailed calculation of Wigner-function propagators and in particular present coarse-grained approximations that can be incorporated into current approaches to the simulation of phase space densities. An explicit numerical example is given in section 5 and conclusions are presented in section 6.

Background and notation
In this section we establish the notation and setting assumed in the rest of the paper. We use a scalar model throughout, even though an important application is to electromagnetic waves, as this allows us to set out the most important features with the least cumbersome notation. We also treat propagation through uniform media, although this assumption is easily relaxed. We thus consider a wavefunction Ψ (x) satisfying the driven Helmholtz equation in a domain Ω. We assume that Ω is (d + 1)-dimensional so that d denotes the dimension of its boundary Σ. Our schematic and numerical illustrations are for d = 1, although the underlying formalism applies also to the more physically relevant case d = 2. The driving term F(x) is intended to represent internal sources in Ω. Boundary driving is also possible through the application of inhomogeneous boundary conditions. We choose to use a boundary representation of both wave and phase space solutions. This leads to conservation rules being stated in terms of power flux rather than energy density and provides, we believe, a simpler description of the underlying formalism. In the case of phase space simulation, representation of solutions as a flux density across surfaces is also the setting used by methods such as DEA which this paper is intended to augment.

Boundary representations of ray densities
We begin by setting the stage for a description of the wave problem in terms of an underlying ray dynamics. The corresponding phase space amounts in a boundary representation to the use of a surface of section, for which we now set out the notation. In much of the discussion to follow, the surface of section is formed by restricting rays to the boundary Σ of Ω. However, it is useful to frame the discussion so that surfaces of section may also be formed by restricting rays to other kinds of section, cutting through the interior of Ω for example: this more general scenario is used in section 3.1 to fix a coordinate frame around the diffractor itself.
We let p denote a (d + 1)-dimensional unit vector in the direction of a given plane wave or associated ray, so that the corresponding wave vector is k = kp, and denote by p the projection of p onto the boundary used to define the surface of section. For d = 1 we can simply let p = sin χ, where χ denotes the angle of incidence of a ray arriving at or leaving the surface. Let s denote position coordinates on the surface so that ds represents arclength for d = 1 or surface area for d = 2.
We let ρ ± (s, p) denote a density on the 2d-dimensional phase space of the surface of section so that J ray ± (s) = ρ ± (s, p)dp (2) represents a local current crossing the boundary at s. Here we label with ± the sense in which the corresponding rays cross the surface and a total current J ray (s) is then defined so that Phase-space simulations are intended in this work to calculate ρ ± (s, p). It is important to distinguish them from a density (x, p) on the full 2(d + 1)-dimensional phase space, which provides an energy density rather than a current when integrated over momentum coordinates. These are related by near a surface of section, where (s, p) are the coordinates on the surface of section of a ray reaching (x, p) and σ ∈ {+, −} labels the crossing direction defined by p. A straightforward geometrical integration verifies that the spatial energy density defined by (x, p) can, when evaluated on the boundary itself, be related to the surface densities ρ ± (s, p) by Note that the geometrical factor √ (1 − p 2 ) seen in the last integral here is also evident in boundary representations of the wave solution described in section 2.2.
Denote by (s, p) = ϕ(s , p ) the surface-of-section map which traces a ray leaving the surface of section at (s , p ) until it returns to the surface of section at (s, p). In simple cavity problems ϕ carries a ray across a chord, but we may also allow for scattering from internal obstacles or bending of rays in inhomogeneous media between surfaces of section, as shown schematically in figure 1. Note that if there is scattering from obstacles, we intend ϕ here to account only for specular reflection-diffractively scattered rays are treated separately in later sections. In a typical scenario, the system may be driven by sources either on the boundary itself or in the interior region between boundaries and we denote by ρ 0 the density returning straight to the boundary from these sources. Use the label + for a density of rays leaving the surface towards the interior and − for the density of rays returning to the boundary. Then a steady-state boundary solution satisfies a consistency equation of the form where Lρ(s, p) = G FP (s, p, s , p )ρ(s , p )ds dp Greek symbols are used to label orbits which cross the domain directly (such as orbit α) or which bounce specularly from internal obstacles (such as orbit β). Roman letters label orbits which scatter diffractively (such as orbit a).
is a Frobenius-Perron operator, defined by the kernel and [10] in which Q(s , p ) accounts for damping (0 Q(s , p ) 1). There is a second relation between ρ − and ρ + fixed by the nature of scattering from Σ. For example, when rays are reflected specularly from it, without loss.

Boundary representations of wave solutions
We will now present the solutions of the full-wave problem in a boundary representation, mirroring the setup of section 2.1. Here we follow the formalism established in the context of semiclassical approximation in [19]. Corresponding exact formulations have been suggested for general quantum potentials in [20][21][22][23] and for cavity problems in [24,25], but to present the main conclusions reached in this paper it is sufficient to consider the semiclassical setting of [19]. We assume that the total wave field Ψ(x) can be decomposed at a boundary into components approaching and leaving it as follows Here ψ − (s) and ψ + (s) respectively denote boundary functions representing wave components arriving at and leaving the boundary, andP is an operator corresponding in semiclassical approximation to the normal component of ray momentum p. This normal momentum oper-atorP can be defined outside of semiclassical approximation in the approaches of [23,24], although these are rather cumbersome for present purposes: it suffices in this paper to assert that, at leading order in semiclassical approximation,P can be represented by the classical ray symbol In the propagating region of phase space (p 2 < 1), and letting χ denote the angle of incidence of a ray hitting the boundary, we can also writê P ∼ cos χ.
Note that this geometrical factor has previously appeared in (5) relating flux to energy density. There is a corresponding relation between the intensities defined by Ψ and by ψ ± . The current crossing the boundary at s is defined here as where ∂/∂n denotes a normal derivative. This current is scaled so that it can be approximated at leading order, and neglecting evanescent contributions, as This provides a wave analogue of (2) and (3) and note that the normalisation of each of the incoming and outgoing wave components then defines a total flux (approaching or leaving the boundary) rather than the more conventional total energy obtained by normalising Ψ. In quantum-mechanical notation, we identify the density operator as being the wave analogue of the boundary densities ρ ± (s, p), whereaŝ is the wave analogue of the full density (x, p), where the outer angular brackets allow for averaging in the case of incoherent wave fields. Concretely, we will work with the two-point correlation functions defined by Note that such field-field correlation functions are an entirely natural way of measuring and characterising emissions from noisy EM sources [26,27] and in the characterization of the propagation channel between large multi-antenna systems, which is important to the design of optimal massive multiple-input-multiple-output (MIMO) systems [28]. The propagation of such correlation functions is treated at a fundamental level in [29], for example, and is extended by using the correspondence with phase space densities to arbitrary multi-reflective environments.

Transfer operators and semiclassical propagation of field amplitudes
Semiclassical propagation of field amplitudes in terms of specular and diffractively scattered rays is well established. Here we set out the main features and notation assumed, in preparation for extending their use to propagate correlation functions and Wigner functions in the rest of the paper.
The basis for all semiclassical approximation in this paper is a representation of the Green function connecting points x and x, defined here so that It is approximated semiclassically as a sum over orbits starting at x and ending at x, of the form [19,[30][31][32] where γ labels topologically distinct paths from x to x, L γ (x, x ) is the corresponding optical path length and the amplitudes A γ (x, x ) are determined by local stability properties of each path. The paths γ may be accounted for entirely by specular reflection or can include diffractively scattered segments: further details about how these are respectively treated are given in section 2.4. Using the Green function to propagate the solution through the interior of Ω leads to selfconsistency conditions on the boundary Σ in the form of a boundary integral equation. This is written here in the form of a transfer operator relating the wave ψ + leaving the boundary to the wave ψ − returning to it and is of the form where f 0 (s) represents a wave component arriving directly at the boundary from a source term derived from the forcing on the right-hand side of (1). Neglecting evanescent contributions, that is, assuming all rays included have p 2 < 1, the kernel T(s, s ) is approximated using the following modification of (12) where χ γ is the angle of incidence of orbit γ as it leaves the surface of section and χ γ is its angle of incidence as it returns to it. The cosine factors here arise because of the normalmomentum prefactor in (7) and result in the transfer operator satisfying unitarity properties under semiclassical approximation [19,24]. The unitarity condition is closely linked to flux conservation. Let us write (13) formally as In a closed, lossless system, and in the absence of internal forcing (so f 0 (s) = 0), the total flux leaving a closed boundary Σ must be balanced by the flux returning to it. If expression (8) for the local current were exact, this would imply that ψ + |ψ + = ψ − |ψ − and therefore that the operatorT is unitary: In fact, it is possible to present exact formulations of the transfer operator for which this is the case (see [23]), but this necessitates a treatment of evanescent components that requires considerable additional technical details: the basis of boundary functions used to achieve this is not obvious or intuitive. Since evanescent waves do not play a significant role in the applications that motivate this work, it is sufficient here to assert that (15) holds semiclassically on the subspace of boundary functions without evanescent components, as established originally in [19].

Specular versus diffractive orbits in the Green function
Let us return to the fundamental semiclassical approximation (12) and state more explicitly how diffractive orbits are accounted for and what the corresponding approximation for the transfer operator is. We label the classical orbits using Greek indices γ = α, β, · · · and diffractive orbits by Roman indices γ = a, b, · · · as laid out in figure 1. We now split up the corresponding contributions to the Green function [30][31][32], using the notation The explicit form of G spec (x, x ) is not needed here, see [19]; we will instead specify the corresponding contribution to the transfer operator in section 2.5. The form taken by diffractive contributions depends on the nature of the diffractor itself and its dimension. In order to minimise notational complexity we restrict our attention here to the case where the diffractor is effectively zero-dimensional (such as scattering by small obstacles or apertures or the tips of edges or wedges in 2D problems, cones in 3D, etc). We also focus on orbit contributions going through a single diffractive event. Then the corresponding contribution to the Green function can be written as [30][31][32] where x a denotes the position of the diffractor on the path of orbit a and D(p out a , p in a ) is a diffraction coefficient, depending on the direction p in a of the scattered ray as it arrives at the diffractor and the direction p out a as it leaves it. Before and after this diffraction event, the orbit undergoes specular propagation described by Green function contributions G spec a (x a , x ) and G spec a (x, x a ). These are of the same form as the first set of summands in (16), but indexed here by the diffractive orbit label a, which specifies also the specular segments before and after the diffraction event itself.

Specular versus diffractive orbits in the transfer operator
There is a corresponding decomposition for the transfer operator as defined in (13), of the form, We begin by presenting the purely specular part, which is simpler to specify for the transfer operator. This of the form [19] T spec (s, where and L α (s, s ) is the path length of an orbit leaving the boundary at position s and returning to it at position s, which is just a chord length in the absence of internal obstacles. The Maslov index μ α accounts for focusing around the ray path if it encounters internal obstacles: μ α = 0 if the path is a simple chord. This takes a standard Van Vleck form used in semiclassical approximations of unitary operators [33].
In order to describe T diff (s, s ), we begin by placing the diffractor at position s a on a virtual surface of section Σ 0 , shown schematically in figure 2; in part (a) of figure 2, Σ 0 contains a scatterer, in the absence of which rays pass directly through it, while in part (b), Σ 0 contains an aperture, without which rays reflect specularly from it. In each case it is useful to define a productT (0) spec =T outT in : Σ → Σ of operatorsT in : Σ → Σ 0 mapping Σ to the virtual section Σ 0 andT out : Σ 0 → Σ mapping the virtual section Σ 0 back to Σ, following either direct transmission through Σ 0 or specular reflection from it, as appropriate. Note that each ofT in andT out can be approximated using a kernel of the form (18), albeit including additional reflection phases for the scenario in figure 2(b). OperatorT (0) spec then collects all the contributions to the kernel T spec (s, s ) from orbits that collide with Σ 0 exactly once: the superscript (0) is used to distinguish it from the complete transfer operator, which may include orbits that that undergo more than one collision with Σ 0 before returning to Σ or that return to Σ without colliding with Σ 0 at all.
Rewriting (17) in terms of the transfer operator formalism using (14) results in the kernel where χ in a and χ out a are, respectively, the angles of arrival and departure in Σ 0 of the diffractively scattered ray a. The factors T in a (s a , s ) and T out a (s, s a ), respectively, denote the contributions to the kernels of operatorsT in andT out of the specular components of diffractive orbit a before and after the diffractive event.

Renormalised transfer operators, unitarity and the optical theorem
Although most applications of diffractive wave propagation are lossy in practice, it is useful to establish formally the symmetries that arise from power conservation in the lossless case. In this section we convert the optical theorem [34] to the notation used by us to describe diffraction in transfer operators: this will be important to us in interpreting diffracted Wigner functions in later sections.
The transfer operator formalism, and the intended application of the final results in this paper to ray-tracing simulations, can be adapted for arbitrary surfaces of section. However, although the semiclassical calculations underlying these results are quite general, the details of quantisation used in the full-wave calculations depend on features such as the topology of Σ: for example, if Σ is closed then momentum p is discretely quantised. In order to minimise notational complexity, we set out the formal structure in this section while assuming that Σ is of infinite extent, so that Σ is replaced by R d (or multiple copies of R d , as in figure 3). This also greatly simplifies the discussion of Wigner-Weyl representation in later sections.
In particular we assume in bra-ket notation the position basis |s and the momentum basis |p to be normalised so that which mirrors standard quantum-mechanical notation except that we have replaced → 1/k to emphasise the intended application to classical wave problems.
In the detailed calculations in this section, the surface of section Σ consists of two disjoint components Σ 1 and Σ 2 , on either side of the diffractor, which is itself contained within a surface of section Σ 0 , as illustrated in figure 3. We use barred coordinates (s ,p ) and (s,p) to locate where rays arrive on Σ 0 from an initial condition (s , p ) on Σ, or where they leave Σ 0 towards a final location (s, p) on Σ, respectively, as shown in figure 3. Note that these barred coordinates can also be regarded as a defining a comoving coordinate system on Σ, in which rays stay at fixed coordinates as Σ 1 and Σ 2 are moved further from the diffractor. The diffracted Wignerfunction that is the object of interest in this paper is far easier to describe when expressed in this comoving coordinate system, and this is the approach taken in the rest of the paper. The formal means of achieving this is set out in section 3.1.

Renormalisation
The combined transfer operator T tot =T spec +T diff accounts for both specular and diffractive propagation. For the simpler geometry in figure 3, there are no orbits returning to Σ 1 and Σ 2 without first colliding with Σ 0 and we can therefore simplifyT spec ≈T (0) spec =T outT in in the notation of section 2.5. We then define the rescaled accounts for diffractive propagation in a representation which removes dependence on the choice of section Σ used to define incoming and outgoing waves. Although exact formulations of the transfer-operator approach are possible, [20][21][22][23][24][25] the operatorD is in practice only ever calculated within the context of semiclassical approximation in this work and the identities above are much easier to describe in that context. Then it suffices to establish the identities above by evaluating operator products within stationary-phase approximation [19]. The resulting approximation forD is motivated next and stated explicitly at the end of this subsection.
Pre and post multiplication byT † in andT † out in (21) undoes the propagation from Σ to Σ 0 and then back to Σ in (20). However, the resulting kernel D(s,s ) = s|D|s is singular, being strongly localised arounds = s a =s , and is therefore awkward to approximate semiclassically (see appendix A). For this reason we work instead with the momentum representation D(p,p ) = p|D|p .
Evaluation of the operator product in (21) is performed in stationary phase approximation and making use of the formal unitarity ofT in andT out , leading to where p in a and p out a are respectively the direction vectors of rays arriving at Σ 0 with momentum coordinatep and leaving it with momentum coordinatep (see figure 3). Note that a detailed derivation of this result uses a stationary-phase approximation in which the diffraction coefficient D(p out a , p in a ) is assumed to be a slow function of its arguments.
By using dp cos χ in a = dΩ in in (22), where dΩ in denotes solid angle of unit vector p in a for d = 2 and dΩ in = dχ in a for d = 1, this can alternatively be expressed Alternatively, starting fromŜ †Ŝ ≈Î would lead to where dΩ out is defined analogously. Note that these last identities must hold exactly if they are to begin an asymptotic expansion for the unitarity ofŜ. We thus recover the optical theorem [34] for diffractive scattering.

Diffractive propagation of density operators and Wigner functions
In this section we assume again the scenario illustrated in figure 3 and work entirely with the renormalised operatorŜ =Î + iD, expressed in the previous section in terms of comoving coordinates (s,p). Because we work exclusively in terms of these comoving coordinates in this section, we drop the bars on them for notational simplicity as there is no need to distinguish them from the original coordinates on Σ 1 and Σ 2 . The aim of this section is to convert the geometric description of diffraction of waves described in sections 2.4 and 2.5 into an explicit scheme for the propagation of phasespace densities across diffracting obstacles. To achieve this we construct a corresponding propagator of Wigner functions [13][14][15][16], so that we can formulate a propagator of phase space densities by averaging it. Note that other representations in phase space of incoherent wave fields are available, such as the Husimi distribution [35]. However, effective propagation then requires that we work with higher-dimensional functions, such as, for example, off-diagonal matrix elements in a coherent-state basis. It is for this reason that a Wigner representation is used in the propagation step: if desired, other representations such as the Husimi distribution can be calculated from the Wigner function once the propagation step is done.
Before describing the technical details behind this, we begin in the next subsection with a summary of the main results achieved, with justification to follow in later subsections.

Summary of results derived in this section
An incident field, generally incoherent, is first characterised by converting its two-point correlation function into an incident Wigner function W in (s, p) [14,16,17]: a formal definition is given in section 4.2. The Wigner function of the diffractively scattered field is then separated into a specular part W spec (s, p) and a diffractive part W diff (s, p), each obtained by applying a corresponding Wigner-function propagator: Using the renormalised conventions in section 3.1, the specular part has the trivial propagator at leading order (see section 4.3). When unfolded back into the farfield coordinates, this reproduces the standard ray tracing evolution by the Frobenius-Perron propagator (6) [16]. The principal result derived in this section is that the diffractive part can, in a coarse-grained approximation that averages interference fringes, be replaced by an averaged kernel in the form This states that the part of the incident field carried by rays arriving at the diffractor location s a are scattered so that they re-emerge from the same location but with rearranged direction, described by the partial kernel Y diff (p, p ). This partial kernel can in turn be written The two contributions in brackets here derive respectively from orbit contributions that are entirely diffractive in nature and from mixed combinations of diffractive and specular orbits. Using the optical theorem guarantees that this combined diffractive propagation contributes no net flux to the scattered field: Ḡ diff (s, p, s , p )ds dp = 0 (26) (see section 3.2). In other words, flux taken out of the stripe s = s a by the mixed specular-diffractive contributions (second contribution in brackets in (25)) is redistributed into the rest of phase space by the purely diffractive contributions (first contribution in brackets in (25)), so that there is conservation of flux overall. The discussion in the following subsections can also be used as the basis for a fuller description of the diffractively scattered Wigner function, in which detailed interference fringes are accounted for. However, the simplified, coarse-grained description above provides the main ingredients required by the use cases we envisage, in which diffractive redirection of energy into shadow regions is to be accounted for in ray-tracing simulations.

Formalities and useful identities
We begin by setting out notation for the formal propagation of density operators and their representations as Wigner functions. Consider a density operatorρ in propagated according tô withŜ =Î + iD: this is consistent with the conventions of section 2.2 if we letρ in =T inρ+T † in andρ out =T † outρ−T out . We can formally write the propagator as a superoperator of the form Each of the contributions here is of the formÛ ⊗V † . We now set out how such transformations of operators (byÛ ⊗V † ) are represented as propagators of Weyl symbols and Wigner functions: in subsequent subsections we specialise these results to each of the contributions in (27) in turn. We use the formal approach of [17], which is centred on the algebra of operator familieŝ R(z) andĤ(z), that respectively quantise inversion of phase space through a point z = (s, p) and translation along a phase space vector z = (s, p). The inversion operatorR(z) can be defined by its action on a wave function as follows, (where we use x as the argument for the wavefunction to distinguish it from the arguments of R(z)). The Heisenberg operatorsĤ(z) are defined bŷ A range of useful identities that are satisfied by these operators are set out in appendix B. In the main text we quote only the important identitŷ where denotes the symplectic form defining area in phase space. Within this formalism, the Weyl symbol of an operatorÂ can be defined by the relation where Tr denotes a trace operation. The complementary functioñ is related to it by a Fourier transform (see (B.5) and surrounding discussion in appendix B). Note that we regard a Wigner function as being simply the Weyl symbol of the corresponding density operator. It is shown in appendix B that the transformation of Wigner functions corresponding tô U ⊗V † is described by a propagator of the form This representation of G(z, z ) associates the Wigner propagator with cycles on the following diagram: For example, going through the operators inside the trace operation in (32) from right to left, R(z) first maps z 1 to z 2 , the time-reversed operatorV † then maps z 2 to z 2 ,R(z ) maps z 2 to z 1 and finallyÛ completes the cycle by mapping z 1 back to the starting coordinates z 1 . Note that R(z) mapping z 1 to z 2 means that z is a midpoint of z 1 and z 2 and z is similarly a midpoint of z 1 and z 2 : such midpoint conditions are common in calculations of Weyl symbols and Wigner functions. The form given in (32) is useful because it allows us to go quickly to semiclassical and other identities using the formalism of trace formulas. These are exploited in the following subsections to evaluate Wigner-function propagators for the individual superoperators in (27).

Specular-specular Wigner propagators
It is clear that the propagator forŜ αα should be but still instructive to see how this emerges from (32) and using the identities in appendix B.

Diffractive-diffractive Wigner propagators: exact identities
The next easiest contribution to treat is from the diffractive-diffractive superoperatorŜ aa , whose Wigner propagator is of the form In practice we propose to use coarse-graining (see following subsections) to get simplyimplemented realisations of this, but we can also find more detailed information by exploiting identities in appendix B. For example, using (B.3) we deduce where W DD † (z) is the Weyl symbol ofDD † , according to (30). Similarly where W D † D (z ) is the Weyl symbol ofD †D , in initial coordinates z . Note that a consequence of the unitarity ofŜ =Î + iD (which is straightforwardly established within semiclassical approximation [19]) is thatDD † =D †D and W D † D and W D † D are then the same function, albeit arising here with different arguments.

Diffractive-diffractive Wigner propagators: coarse-grained approximations
Because the operatorD is spatially localised (on the manifold Λ according to the discussion in appendix A), this should be approximable in a coarse-grained sense by an operator of the form That is, the propagatorḠ aa (z, z ) approximates the action of G aa (z, z ) when acting on densities that vary over scales that are large compared to those of typical interference fringes. The function Y aa (p, p ) that achieves this is obtained by integrating G aa (z, z ) over its spatial arguments: By applying identity (B.6) to each of these spatial integrations we get Note also that the coarse-grained diffractive-diffractive propagator gives a necessarily positive contribution to the outgoing density. To maintain overall flux conservation, this must be counterbalanced by the action of mixed (diffractive-specular and specular-diffractive) propagation, which is considered next.
Note that we can also express this as whereW D (ζ) denotes a Fourier transform of the Weyl symbol ofD, according to (31).
As with the diffractive-diffractive propagator, this requires some coarse-graining to be easily implemented in practice, but exact identities are also easy to deduce here. For example, from (B.5) (with −ζ replacing z in the identity) and then using (30) to define the Weyl symbol ofD. (Actually, the same result can be obtained by directly integrating the top line of (36) and invoking (B.3)). Similarly, Finally, each of these integrals has been obtained by integrating over ζ while respectively holding z and z fixed. A third alternative is to fix while integrating over ζ and to regard as defining a function ofz and ζ. Then, since we can also write One further useful identity before embarking on approximations is that the specular-diffractive propagator is the complex conjugate of its diffractive-specular counterpart, that is, Therefore, the combined propagator is real and obtained by just taking an imaginary part of a single trace calculation. More explicit forms for this propagator can be found in appendix C; see in particular (C.5).

Diffractive-specular and specular-diffractive Wigner propagators: coarse-grained results
As for the diffractive-diffractive propagator, we make progress by asserting that G aα (z, z ) is approximately given-in a coarse-grained sense-by an operator of the form The justification for making this replacement is more subtle than for the diffractive-diffractive contribution, however, and is discussed in more detail in appendix C. Once we accept the coarse-grained ansatz above, the function Y aα (p, p ) is once again obtained by integrating G aα (z, z ) over its spatial arguments By applying identity (B.6) twice to the top line of (36) we then deduce that We get an analogous identity for the momentum part of G αa (z, z ) and we then deduce that is a coarse-grained approximation for the net contribution of mixed diffractive-specular contributions to the Wigner propagator. Note that this can also be written as so that this part of the propagator can be interpreted as subtracting density from a narrow column in phase space over s a (assuming Im (D(p, p)) > 0 as implied by the optical theorem) from the purely specular part of the propagator given in (34).

Net contribution of diffractive components
Note that combining all diffractive contributions gives a coarse-grained propagator where Y diff (p, p ) has already been defined in (25). This completes the derivation of the coarsegrained diffractive Wigner function propagator asserted in section 4.1.
The action ofḠ diff (z, z ) on a density incident on the diffractor is now easy to interpret-and to implement in practice in the context of the DEA phase-space simulation method [10]. In DEA a phase space density is represented using a mesh for spatial dependence and a polynomial basis in momentum. Spatial elements are typically large compared to the local wavelength, so use of the coarse-grained propagator is appropriate. We then transform the part of the incident density over the spatial element containing the diffractor, which is a function of momentum projected onto the polynomial basis in use. The operator with kernel Y aa (p, p ) is then used to map this incoming function of momentum into an outgoing function of momentum, which is in turn re-inserted into the outgoing density over the same spatial element. We must project the kernel Y(p, p ) onto the polynomial basis used for momentum to achieve this, but this is a reasonably straightforward operation and should be much faster than other typical calculations necessary for DEA matrix formation.
The explicit form given for Y diff (p, p ) in (25), and its derivation in this section, assumes that GTD remains valid for forward scattering. In more general diffraction scenarios, such as edge or wedge diffraction, specular and diffractive contributions must be treated jointly in the forward scattering direction using UTD. This more general case will be the subject of a future publication. Here we point out that in this more general case the structure given in (40) will remain valid, except that Y diff (p, p ) will take a different form than given in (25). Except for this change in momentum kernel, the essential way in which the method introduces diffraction to phase pspace propagation, such as described in the previous paragraph, will remain the same.

A numerical example
The calculations are illustrated in this section using a concrete numerical example. We examine scattering from a circular obstacle in a two-dimensional domain (with Dirichlet boundary conditions) of a source density corresponding in phase space to which is confined to an interval of length L = 2b in s and to the propagating region p 2 < 1 of momentum space. Note that in practical terms this source Wigner function is obtained by solving a wave problem in which the inhomogeneous Dirichlet boundary conditions lim (x 1 ,x 2 )→(s 1 ,s 2 ) are imposed on the density operator defined in (11), along the surface-of-section component Σ 1 of the geometry illustrated in figure 3. Filtering out the evanescent components, which contribute negligibly to far-field scattering, and conjugation by the square root of normal momentum operatorP then leads at leading order to the Wigner function above. This density is described in terms of regular coordinates on Σ, before the renormalisation process described in section 3.1: to use the results described in section 4, it will be necessary to convert from the comoving coordinates assumed there. It is sent towards a circular obstacle of radius R that is a distance d from the incoming surface of section Σ 1 and the forward-scattered density is recorded as it passes through surface of section Σ 2 that is distance d from the obstacle on the other side. In figures 4 and 5, these distances are chosen so that kb = kd = 150 and kR = 3. Note that in the coarse-grained approximation (24) and (25), kR need not itself be small as long as R is smaller than the length scale used to resolve the phase space density (mesh elements are typically many wavelengths wide in DEA) [10]: in this example, the obstacle is approximately one wavelength in diameter.
This source density (on Σ 1 ) and the corresponding forward-propagated Wigner function (restricted to Σ 2 ) are shown in figure 4, parts (a) and (b) respectively. The two main features are (i) the shearing of the propagated Wigner function that occurs as a natural result of free-space propagation of classical densities [16] and (ii) the dark gap running through the lit region. This gap arises from the removal of forwards flux after scattering from the obstacle and should be explained by (24) and (25), following conversion from comoving coordinates. Note that the coarse-grained calculation aims to predict the overall weight of this gap, but does not attempt to recreate the detailed interference fringes that are visible around it.
To see the diffractive contributions directly, we show in parts (c) and (d) of figure 4 the backwards scattered Wigner function and the difference in the forwards direction between the propagated Wigner densities with the obstacle present and with the obstacle removed.  ) characterising the source. It shows only the propagating part so is localised in p 2 < 1 and is also confined within a finite interval in physical space. In part (b) is shown the total Wigner function emerging from the other side of a circular obstacle. Note the dark gap running through the lit region, which corresponds to forwards flux removed by the obstacle. In parts (c) and (d) we show diffractive contributions in both the backwards scattering direction (in (c)) and the forwards scattering direction (in (d)). The latter is net negative (coloured blue) but is balanced by the positive contribution to the backwards scattering direction, according to (26). On a coarse-grained level, these plots look entirely consistent with (24) and (25), once the difference in coordinate systems is taken into account: the integrated flux in the forwards direction is negative (negative values of the Wigner function are plotted in blue), but this is balanced by the net positive flux in the backwards scattering direction. In both plots, the diffractively scattered Wigner function is localised around a manifold of orbits originating from s = s a on Σ 0 . We verify this quantitatively by comparing the predicted current in physical space, defined as where W diff (s, p) is defined in analogy with (23). This is seen from (24) and (25), after from some further manipulation and conversion from comoving coordinates, to take the form where and we have specialised the general result to the 2D case where dΩ in = dχ in a . The arguments (s (χ in ), p (χ in )) of the source Wigner function denote the initial coordinates of rays arriving at the scatterer with angle χ in . A comparison of this result against a full-wave calculation is given in figure 5, showing good agreement. As expected the net current is negative near the forward direction (part (a)) where the mixed terms of section 4.7 contribute. In the backward scattering direction (part (b)), there is only a contribution from the completely diffractive terms from section 4.5, which is intrinsically positive.
We conclude this section by exemplifying how corresponding calculations may be exploited in the context of mobile communication modelling. Multi-antenna techniques are often adopted in wireless communications at microwave, mmWave, and optics communications. An important question in this context is to estimate optimal rates of information exchange between transmitter and receiver. Answers to this question have been given in terms of coupled phasespace volumes for communication in free space [36]. In extending such theory to more complex propagation spaces, reflection, diffraction and diffusion will play a key role in the estimation of fundamental bounds of MIMO systems. The formalism being developed here will allow integration of diffraction within DEA simulations of propagation environments hosting arbitrarily complex transmitting and receiving antennas.

Conclusions
We have described how diffractively scattered orbits may be incorporated into Eulerian descriptions of ray propagation, by calculating corresponding contributions to a diffractively scattered Wigner function. In applications this approach has an advantage over traditional ray-tracing simulations of bypassing the typically exponential increase of ray number with length, which is a feature of both specularly and diffractively propagated ray orbits.
Concrete illustrations have been provided in the context of wave scattering from a small obstacle in two dimensions, but the theoretical results apply also in three dimensions and the underlying formalism is intended to extend to other diffraction scenarios. The formalism applies to incident wave fields that are incoherent, characterised by a field-field correlation function rather than a simple wave amplitude, which provide a natural wave analogue of ray densities in phase space through the Wigner function. Although the general results given in section 4 can in principle be used to predict detailed interference fringes in the diffractively scattered Wigner function, we have emphasised a simpler, coarse-grained propagation rule given in (24). This is a key result of the paper. The coarse-grained propagator predicts in broad terms where diffractively scattered wave energy is sent, without the complication of describing detailed interference patterns which average to zero-this is in keeping with the underlying goal of introducing diffractive contributions into ray-tracing simulations based on propagating phase space densities.
An important feature of this diffractive correction to ray-traced densities is that it conserves overall flux, where appropriate. Flux conservation is seen as a direct consequence of the optical theorem applied to the diffraction formalism. It guarantees exact flux conservation even for the coarse-grained propagator in (24).
A significant limitation of the results presented so far is that they are based on GTD approximations. In more general scenarios, such as edge or corner diffraction, an explicit description of the forward-diffracted density requires us to use the UTD. This will modify the kernel given in (25) in the forward direction, where specular and diffractive orbits merge, but the main qualitative features will persist. Extending the detailed calculations to problems where UTD is required will be the subject of future work.
In other words, the operatorD is thus localised at the position of the diffractor, which is what we expect intuitively.
We also note that expressing the density as a 2d-form (p,p )dp ∧ dp = k 2π 2d D(p out , p in ) 2 4k 2 dΩ out ∧ dΩ in lets us evaluate it simply in terms of the angles of arrival and departure.

Appendix B. Identities for the evaluation of Wigner functions and their propagators
In this appendix we list key identities satisfied by the reflection and translation operators defined in section 4 and used there to derive central results for diffractive Wigner-function propagators. The use of the algebra of these operators, along with the formalism of doubled phase space exploited in appendix A, to describe Wigner functions and Weyl symbols, has previously been described in [17] and is used here with slightly altered notation. The reflection operator can be formally written From these formal identities we can deduce (28). We should also note the completeness relations is here obtained as the Weyl symbol (in z ) of the operator (k/π) dV †R (z)Û.

Appendix C. Justification of spatial localisation assumed for coarse-grained propagators
The localisation assumed in (39) is not immediately obvious. Formal identities such as (36) let us establish easily that G aα (z, z ) is localised around s = s , but it is less clear that s and s should be localised around s a . In fact we need averaging to achieve this. To establish these features, it is useful to visualise the trace formula (32) as a diagram, here for the special caseÛ = iD andV =Î: This diagram shows a cycle z 1 → z 1 → z 2 = z 2 → z 1 → · · · contributing to a trace formula for (36), where the leg z 1 → z 1 is effected by iD and so constrains s 1 = s a = s 1 and the leg z 2 → z 2 is effected byÎ and so, in particular, constrains s 2 = s 2 but not s 2 = s a or s a = s 2 . Since z and z are respectively the midpoints of (z 1 , z 2 ) and (z 1 , z 2 ), we then similarly constrain s = s but not s = s a or s a = s .
To explain the further localisation around s = s a and s a = s assumed in (39), we must examine the trace in (36) in more detail. The action in the periodic orbit above is the area of a triangular region formed by the vertices (z 1 , z 2 = z 2 , z 1 ), which is easily seen to be of the form A(z, z ) = 2(p − p ) · (s − s a ). (

C.2)
This action provides the dominant phase oscillation of G aα (z, z ). If we perform averaging of the diffracted fields, G aα (z, z ) then washes out except near the lines p = p or s = s a where this action vanishes. This is at the root of the additional delta functions G aα (z, z ) ∝ δ(s − s a )δ(p − p ) in the coarse-grained propagator. We can justify this assertion more explicitly by evaluating the trace in (36) using momentum representation, to get e −ikΩ(ζ,z ) Tr DĤ † (ζ) = e 2ik(s ·p−p ·s) R d dp R d dp D(p , p ) p |Ĥ(ζ)|p * . (C. 3) The Heisenberg operator matrix elements can be shown to take the form p |Ĥ(ζ)|p = e 2ik(p−p )·(s−s )−2ik(s−s )·p δ(p − p − 2(p − p )) and following insertion into the integral above, followed by a change of variable, we can express where the approximation refers to coarse graining in u = s − s . Note that the phase factor here does indeed reproduce the triangle area predicted in (C.2).