A Phase-Space Approach for Propagating Field-Field Correlation Functions

We show that radiation from complex and inherently random but correlated wave sources can be modelled efficiently by using an approach based on the Wigner distribution function. Our method exploits the connection between correlation functions and theWigner function and admits in its simplest approximation a direct representation in terms of the evolution of ray densities in phase space. We show that next leading order corrections to the ray-tracing approximation lead to Airy-function type phase space propagators. By exploiting the exact Wigner function propagator, inherently wave-like effects such as evanescent decay or radiation from more heterogeneous sources as well as diffraction and reflections can be included and analysed. We discuss in particular the role of evanescent waves in the near-field of non-paraxial sources and give explicit expressions for the growth rate of the correlation length as function of the distance from the source. Furthermore, results for the reflection of partially coherent sources from flat mirrors are given. We focus here on electromagnetic sources at microwave frequencies and modelling efforts in the context of electromagnetic compatibility.


I. INTRODUCTION
Predicting the properties of wave fields in complex environments is an extremely challenging task of crucial importance to a wide variety of technological and engineering applications, such as vibroacoustics [1] or electromagnetic (EM) wave modelling. In particular, characterising the radiation of EM sources reliably, both in free space and within enclosures, is a longstanding research issue. In the context of electromagnetic compatibility (EMC), digital circuits and large printed circuit boards (PCB) embed thousands of electronic devices and metallic tracks and can produce fields reaching dangerous but hard-to-predict levels [2].
In this paper, we set out an approach for propagating such complex and statistically characterized wave fields exploiting Wigner distribution function (WDF) techniques. This approach has its origin in quantum mechanics [3], but has more recently found widespread attention in optics, see [4][5][6] for an overview. The WDF formalism offers a direct route to pure ray-tracing approximations in an operator implementation [1], while still capturing in its exact formulation the full wave dynamics. The formalism allows one to efficiently treat radiation from complex sources, often having a significant nondeterministic, statistical character.
The method introduced below exploits a connection between the field-field correlation function (CF) and the WDF [9,31]. Both quantities have been studied intensively in the physics and optics literature. For wave chaotic systems, Berry's conjecture postulates a universal CF equivalent to correlations in Gaussian random fields [10,11]. Non-universal corrections can be retrieved by linking the CF to the Green function of the system [12][13][14][15]. In this paper, we describe how field-field CFs can be efficiently propagated using ideas based on ray propagation in phase-space. We discuss furthermore nonparaxial effects as well as including near field effects due to evanescent wave contributions. A systematic expansion of the Wigner function propagator including next-toleading-order effects in the propagating regime leads to Airy-function integral kernels containing the ray-tracing propagator in the small wavelength limit, akin to the treatments in [31,33]. We also show that our WDF representation confirms the validity of the generalised form of the Van Cittert-Zernike (VCZT) theorem discussed in [16]. We give a natural extension of this generalised VCZT for non-paraxial sources and in the near-field region where evanescent waves play a prominent role.
We illustrate these techniques in the context of applications in EMC and related issues. Here, the system under investigation represents a high-density interconnect of integrated electronic circuits. Simulating EM field distributions in a reliable way is highly topical; in addition, the wave correlation function can be measured explicitly in this regime [21], thus providing the necessary input information for numerical simulations. The system under consideration in this paper consists of a series of parallel tracks carrying partially correlated currents, and mimics the typically very complex EM sources found on PCBs. The method has much wider application, however. In particular, when combined with fast phase-space propagation methods such as the Discrete Flow Mapping techniques developed in the context of vibro-acoustics [1,8], the proposed WDF approach offers an ideal platform for developing a universal high-frequency simulation method.

II. PHASE-SPACE REPRESENTATION OF CLASSICAL FIELDS
Radiation from simple EM sources such as antennae can be characterised deterministically through classical electrodynamics methods [18]. Even though such sources are regular and homogeneous, efficiently predicting farfield emission from the near-field pattern requires nontrivial effort if the sources are extended over many wavelengths [19]. EM sources are becoming increasingly complex, however, and the problem of radiation from digital circuits or PCBs presents even greater challenges. Modelling such sources deterministically is often infeasible due to the complexity of the structures, whose details may not even be known in practice. Each component of such a complex EM source is typically driven by unknown sets of quasi-random voltages, subject to fast transients [20]. This is due to the presence of a multitude of electronic components whose switching behaviour depends on the instantaneous operation mode of the circuit, and whose excitation signals are intrinsically random, or highly sensitive to frequency [21]. Consequently, the physical investigation of these scenarios challenges existing analytical and numerical techniques, and calls for more sophisticated modelling tools.
It is thus natural to use statistics as a language for describing the radiation from such complex sources. Specifically, we do not attempt to characterise or propagate the field itself, which is typically hard to obtain in practice, but rather its two-point CF. It has been demonstrated in [21] that corresponding measurements are feasible in the context of emission from electronic devices and PCBs. Here we describe the basic elements needed to use such measurements as input for a practical algorithm with which to predict field intensities and correlations away from the source. Initially we consider radiation into free space in Sec. IV by studying a simple model source and a more realistic source obtained from a full field simulation. In Sec. V, we describe an application to a problem with reflecting boundaries, which is a first step towards our ultimate goal of extending the method to propagation of CFs in more complex environments such as cavities and larger structures.
We start from a planar source at z = 0, parametrized by coordinates x = (x 1 , · · · x d ) with d = 1 or 2 in general, and radiating into in the half-space z > 0. We aim to predict the CF for z > 0 under the assumption than it can be measured (or otherwise modelled) near the source screen z = 0. Here, . denotes an ensemble average over different source field correlations such as a time or frequencyband average. Furthermore, ψ(x, z) denotes one of the tangential field components in the frequency domain.
The results easily extend to cross-correlation between different components.
In the past, the focus has often been on predicting the propagation of probability density functions (PDF) of waves passing through time-domain random [22] or turbulent [23] media. In our approach, the propagation itself is treated deterministically, whereas the radiation from the source is characterised statistically. This can be done, for example, by measuring the spatial field along a surface close to the source and determining the source CF by averaging the signal over time. We thereby eliminate statistical fluctuations carried by the wave fields by ensemble averaging physical observables over suitable parameters.
We now present the CF propagation rule explicitly for single field components. Polarisation effects can also be accounted for by propagating the field-field correlation tensor, which can be derived from the dyadic free-space Green's function [7,24].
The field in the region z > 0 is naturally presented in terms of the partial Fourier transform, where ψ(x, z) denotes a field component on the screen itself (z = 0) or to its right (z > 0) and k is the wave vector. The radiated fields can then be reconstructed using the evolution of this partial field. This can be calculated by using the dyadic second Green identity which, in a source-free region, becomes the dyadic version of Huygen's principle [25]. Being a convolution integral, the partial Fourier transform of the surface integral transforms to an algebraic equation. Then, the boundary conditions given by the fields sampled in the near-field region of the source can be used to eliminate the magnetic field in such an equation. The result of this procedure, restricted to the electric field components parallel to the source plane, is the following inhomogeneous plane-wave solution where For the moment, we neglect waves incident from the right and thus only describe radiation from a strong, directional source; including incoming waves at the interface can be introduced formally using the boundary integral equations according to the discussion in [26]. An example of this scenario, involving a planar reflector beyond the source, will be given in Sec. V. Here, p = (p 1 , . . . , p d ) takes the meaning of a momentum tangential to the d dimensional source plane. In the ray-dynamical limit, we may identify where the angle α describes the direction of the ray with respect to the local outward normal to the source. In this perspective, T (p) represents a generalized kinetic energy of the ray. The case |p| 2 > 1 in (3) corresponds to evanescent propagation, which does not contribute to the far-field, but may be detectable in the near field; see also the discussion in Secs. IV B, IV C. In order to represent wave fields in phase-space using canonical coordinates (x, p) parallel to the source plane, we define the WDF Upon insertion of (2) in (6), and by exploiting the inverse transformation to represent the source correlation (at z = 0) in terms of the source Wigner function W 0 (x, p), we find This provides us with a propagator of the Wigner function taking the form where the δ-function represents translational invariance in x and the corresponding conservation of momentum. Eq. (8) provides a scheme to propagate wave densities in phase-space for arbitrary sources, no matter how complex or rapidly varying. The propagation of the correlation functions themselves can subsequently be retrieved by an inverse Fourier transform of (7). That is, where x = (x A + x B )/2, and s = x B − x A . The intensity I z as function of the distance z can be retrieved using [5]

III. RAY TRACING APPROXIMATIONS
Asymptotic approximation of the propagator (8) leads to a direct propagation method for the WDF in terms of rays [6,[27][28][29]31]. We will give a derivation of this ray limit below and will also discuss more subtle wave effects such as evanescent decay into the near-field and higher order (in 1/k) wave corrections.
The simplest ray-based approximation is obtained under the assumption that the CF is quasihomogeneous at the source, that is, Γ 0 (x B , x A ) = Γ 0 (x + s/2, x − s/2) varies with x on a larger-than-wavelength scale. In that case, significant contributions to (8) are obtained only for small q and we can expand the phase difference ∆T (p, q) = T (p + q/2) − T * (p − q/2) around q = 0.
In the region |p| 2 < 1 corresponding to propagating waves, the difference ∆T receives contributions only from odd powers of q. Neglecting cubic and higher order terms we find that This is the Frobenius-Perron (FP) propagator [30] for radiation into free space and leads to the evolution [31] of the WDF in the region |p| 2 < 1. This approximation is equivalent to identifying the propagation of the WDF with the propagation of phase space densities along rays according to the evolution The paraxial approximation is obtained by using the linearized flow in the regime |p| 2 ≪ 1 [6]. Outside this regime, the full flow has asymptotes along |p| 2 = 1, the transition to evanescent propagation. Note that inserting Eq. (12) in Eq. (9) and including only the contributions from the propagating region |p| 2 < 1 leads to In the paraxial regime, that is for |p| 2 ≪ 1, and for quasihomogeneous sources, Eq. (14) retrieves the well-known van Cittert-Zernike theorem (VCZT) or generalisation thereof [16,32]. In the next section we will show how to extend the VCZT to the non-paraxial regime and including evanescent waves. Expanding ∆T further in q leads to an improved propagator which is capable of mapping less homogeneous CF's. Including cubic terms leads in the 2D case, for example, to the Airy form where δ a (u) = aAi(au), with a = 2k/(kzT ′′′ (p)) 1/3 and Ai denotes the Airy function. Note that lim |a|→∞ δ a (u) = δ(u), so the FP form is obtained in the limit of large wavenumber k as expected. Similar results have been obtained in the context of the propagation of EM waves through inhomogeneous media [33]. Further improvements over the basic FP propagation (11) are obtained by accounting for evanescent decay into the near-field, which emerges from contributions |p| 2 > 1 in Eq. (7). Since the kinetic operators in Eq. (8) now add constructively, the leading contribution is formed by the zeroth order term in the expansion of ∆T (p, q) and we obtain Improved approximations may be achieved by treating the exponent beyond leading order, but we find that (16) gives a good description of evanescent decay already as discussed in the next section.

IV. RADIATION INTO FREE SPACE
We now test the effectiveness of the FP propagator in the simple case of radiation into free space. The first example treated in Sec. IV A assumes a quasi-homogenuous source distributed according to the Gauss-Schell model. We will examine the near-field behaviour in more detail in Sec. IV B, considering in particular the limit of completely uncorrelated sources. In a second example in Sec. IV C, we consider a more complex set-up mimicking the more realistic sources expected in typical EM applications. We will restrict ourselves in these examples to 2D models (so d = 1) and characterise the behaviour of field-field correlations by focusing on the propagation of one field component along z. We select the field tangent to the source.

A. Propagation of Gauss-Schell model
Using a simple 2D model for the emission of partially coherent EM radiation, we assume a source correlation in terms of a truncated 1D Gaussian-Schell model where l is the length of the source. Here, the characteristic functions account for the finite size of the source. The quasihomogeneity condition can be expressed through demanding σ s ∼ λ ≪ σ x , where λ = 2π/k is the optical wavelength. The source WDF is then found to be For extended sources, for which l ≫ λ, and for x inside the region occupied by the source, Eq. (19) simplifies to Fig. 1 shows the WDF in phase-space at z = 0 for a spatially extended source [37]. Here, and in all other computations with the Gauss-Schell model, we work at a frequency of operation of 1 GHz corresponding to λ = 0.3 m and choose σ x = 1.0 m, σ s = 0.1 m. Fig. 2 shows the propagation of (19) as computed through the full integral operator (7) together with the propagation obtained by the Frobenius-Perron approximation (12). There is surprisingly good agreement between the exact and approximate behaviour even far from the paraxial regime. This is remarkable given that the ray tracing approximation is only valid to leading order. This constitutes a major computational advantage as the FP approximation reduces an integral equation to a simple coordinate transformation. The overall behaviour shown in Fig. 2(a) and Fig. 2(b) reflects the distribution shearing due to the geometrical ray propagation based on Eq. (13); see also [6]. The CFs can now be obtained by a back transformation according to Eq. (9) and are shown in Fig. 2(c) and Fig. 2(d).

B. Non-paraxial Van Cittert-Zernike theorem
In the following, we will focus on near-field effects for small distances from the source as a function of the source correlation parameter σ s . We are in particular interested in how the correlation length propagates in the near-field before reaching the linear VCZT regime. In the near-field limit, the WDF shows exponentially decaying evanescent components according to (16), while the WDF remains essentially unchanged for the propagating part |p| 2 < 1. This leads to a model for the WDF with source distribution (20) of the form Far enough from the source, such that evanescent components have completely decayed, while close enough that evolution in the propagating region of phase space can still be neglected, we model the WDF using Using the inverse Fourier transform, Eq. (9), we now obtain the correlation function from the WDF given by Eq. (21) or Eq. (22) acccordingly. In Fig. 3 we show the resulting near-field evolution of the correlation function, placing the midpoint x = (x A + x B )/2 = 0 at the centre of the source. One observes that in the near-field regime, the width ∆s of the correlation function is smaller than λ, but it increases rapidly towards λ as z approaches and exceeds In the presence of evanescent waves the correlation width increases rapidly before the correlation function becomes a sinc function, whose width then increases linearly according to the VCZT.
λ. The second moment of the correlation function is not defined and cannot therefore be used to define a correlation length. Instead, we define the correlation lengths to be the spacing at which the correlation has fallen by a factor 1/ √ e: Note that for a Gaussian correlation function such as assumed for the source in (17), this definition coincides with the standard variance: ∆s = σ s . With the definition adopted above, we can now obtain the correlation lengths from exact wave propagation calculations. The results are shown in Fig. 4 as a function of the distance z for different source correlation lengths σ s . The universal regime is shown as the blue dashed line. The VCZT regime starts around kz > 1.
From Eqs. (21) and (22), one can estimate the growth rate both in the near and the far field. Including nonparaxial effects, there are three different regimes: i) in the deep near field, with kz < 1 and kσ s < 1, the correlation length increases linearly with a slope that is independent of the frequency as well as of σ s and σ x ; ii) a no-growth regime with ∆s = const.
The latter regime has already been described by Cerbino for paraxial sources [16]; A plateau is observed between kz ≈ 1 and the Rayleigh range discovered in [16,17].
iii) the VCZT regime for large z with a linear growth of the correlation length according to with a slope depending on the ratio of wavelength to source dimension ℓ.
We now motivate these three regimes in more detail, beginning with case (i), which corresponds to kσ s < 1, kz < 1. The WDF W z (x, p) described by (21) then decays slowly along the p axis as |p| increases beyond the propagating region |p| 2 < 1. In the extreme nearfield the correlation function is proportional to the inverse Fourier transform of the function The correlation length defined by Eq. (23) then takes the form That is, we find in regime (i) that evanescent decay of the sub-wavelength correlations in the source dominates in such a way that there is a universal growth rate in the correlation length. The numerical value of the slope in (25) is particular to the form taken in Eq. (23) for the correlation length, but the qualitative conclusion applies more generally. The presence of evanescent waves thus leads to a rapid increase of the correlation length in the near field in this regime. This is important for sources that show fluctuations on scales smaller than the wavelength, such as in the case of a fully uncorrelated source σ s = 0, which may serve as a model for thermal sources [35]. The plateau behaviour corresponding to regime (ii) arises when z is sufficiently large that (22) describes the WDF, while kσ s ≪ 1. The correlation function is then proportional to the inverse Fourier transform of p. In this case the correlation length defined by Eq. Finally, regime (iii) applies once evolution of the phase space takes effect in the propagating region |p| 2 < 1. Assuming the quasihomogeneous case σ x ≫ l and considering first only the near-field region kσ s ≪ kz ≪ 1, then for a given midpoint x the finite size of the source reduces the support in p of the Wigner function and (26) is replaced by For simplicity consider the case x = 0. Then the correlation function obtained from the inverse Fourier transform of this function is Γ z (s) ∼ 1 π sinc ks 1 + (2z/l) 2 and the correlation length defined by (23) takes the form generalising (27). In the farfield z ≫ l, we find (where the numerical prefactor is particular to the convention (23)). Alternatively, if σ x ≫ l, then the screen length becomes unimportant and σ x provides the length scale appropriate to the source intensity. An analogous calculation then allows us instead to recover the basic form

C. Application to a complex source
The field at the source at z = 0 is often produced by a complex process such as tracks on a printed circuit board or integrated circuits in electronic devices; the radiation produced in the source region then propagates into free space. We model such a complex source here by a set of N metallic wires driven by random time-domain voltages, as illustrated in Fig. 5; a realization of the voltage s p (t) driving a pin of the bundle is reported in Fig. 6 along with its spectrum S p (f ). This represents a typical problem in EMC where one tries to obtain statistical information about an erratic signal.
The presence of a perfect electric conductor along an oblique plane makes the source radiate only into the halfspace z > 0: this mimics a configuration that is widely used in the design of printed circuit boards. We use N = 12 wires very close to each other and to the metallic plane in terms of wavelength. Therefore, it is reasonable to think of this circuit as a collection of random sources of partially coherent radiation. The exact fields emitted from such a complex structure are computed through an in-house Transmission Line Matrix (TLM) code [40]. This is a time domain method for modelling 3D electromagnetic field interactions with complex structures that may include a variety of materials. The technique is based on the equivalence between electric and magnetic fields and the voltages and currents on a network of transmission lines. After discretizing space, the fields in individual cells are modelled by transmission lines incident from each cell-face and intersecting at the cell centre forming a junction. Each of these orthogonal transmission lines allows for the propagation of electromagnetic waves. The waves are characterised by voltage and current and their associated electric and magnetic fields. In order to obtain the desired correlation functions, we sample the numerically obtained fields in a plane above the tracks at different times in order to create a suitable ensemble of uncorrelated circuit realisations. This is used as a basis for calculating field-field correlation functions and their Wigner functions both in the near-and far-field. Figs. 7 (a) and (b) show the comparison between the WDF as computed through the full-wave (TLM) simulations, and the WDF obtained by the FP approximation (12) in the far-field at z = 2.3λ. In the TLM calculation, the full time-dependent field is propagated out from the source, while in the FP approximation, the WDF obtained from the signal at the source (as shown in Fig.  6, see also Fig. 8 (a)) is propagated according to (12). There is good agreement between the behavior predicted from full-wave simulations and the FP approximation, even though the source exhibits strong inhomogeneities. Interestingly, Fig. 7 shows the same Wigner distribution shearing as in Fig. 2 following the geometrical interpretation (13) of the correlation propagation. It is worth stressing that such a Wigner function challenges the FP approximation (12), whose underlying assump-  tion is quasi-homogeneity. Note that we can always also switch to the exact transport rule (7), which is computationally more expensive than the FP approximation, but still orders of magnitudes faster than a full TLM calculation. Propagated CFs as shown in Fig. 7 (lower plots (c) and (d)) are finally obtained by applying the inverse Fourier transform (9).
Note that we also find a pronounced broad side radiation around p ≈ ±1 (corresponding to α ≈ ±π/2), and a strong asymmetry of the Wigner distribution due to the oblique metallic reflector. Those features can be captured by inspection of the WDF representation in phase-space, while they are less apparent in the propagated correlation function shown in Fig. 7 (c) and (d).
The source distribution W 0 (x, p) as obtained from the radiated signal, see Fig. 6, is shown in Fig. 8 (a). Note that the region with |p| 2 > 1 corresponds to evanescent contributions. In Fig. 8 (b)-(d), a comparison between WDFs as computed through the full-wave (TLM) simulations and those obtained using the WDF propagator incorporating evanescent contributions, Eq. (16), are shown along the line x = 0. We find that propagation beyond z = 0.1 λ results predominantly in an exponential reduction of the WDF in the region |p| 2 > 1. In the far-field, the radiation energy is restricted to the phasespace region |p| < 1, as can be seen in the WDF in Fig.  7 (a) and (b). A comparison of (TLM) simulated and (Frobenius-Perron) approximate far-field propagated energy is shown in Fig. 9. We see that the two numerical methods show qualitatively the same features, however, there are quantitative differences. We think that these deviations are due to a difference in the numerical treatment of the boundary conditions at x = ±3 m. While the FP approach has no difficulties in treating these boundaries as completely open, the TLM method needs to model this with absorbing boundary conditions. These conditions tend to be still slightly reflective, as is evident from the source distribution in Fig. 8 (a) around x = ±3, p = ±1. This comparison highlights another advantage of the Wigner function propagation method. The evaluation of the WDF of actual circuits can be done for the full electromagnetic field by using the approach described here component by component.

V. REFLECTION OF PARTIALLY CORRELATED SOURCES
Having developed a framework for the propagation of CFs in free space, we are now interested in tackling the case of reflection from planar boundaries. In particular we would like to test the FP approximation in the presence of interference. It is then interesting to solve the canonical situation depicted in Fig. 10, where a planar reflector is located at distance z = L from the source at z = 0. The reflecting boundary is here for simplicity assumed to be parallel to the source plane, indefinitely extended in thexy-plane, and made of an ideal perfect electric conductor (PEC). Therefore, for electric (TE) or magnetic (TM) fields perpendicular toẑ, the Fresnel reflection coefficient reads r (α) = −1, for all incoming angles α [38]. We again consider for simplicity only a scalar field, or a single component of the vector field, emitted from the source.

A. Theory
Consider a plane located at an arbitrary longitudinal coordinate z = D between source and detector. The field distribution in the plane consists then of two contributions: the direct wave coming from the source, and the reflected wave bouncing off the reflector back to the source, that is, where φ (p, 0) is the field at the source plane z = 0, T (p) is defined as in (3), and ∆ = L − D. The momentum space CF is formed as the product of the two fields in (30) and an ensemble average is taken as in Eq. (1). By plugging the closed-form expression (30) into the definition of the WDF (6), we find the phase space representation where the first two terms are direct and reflected contributions respectively, coming to the detector straight from the source or through the reflector, and the last two terms express the interference between direct and reflected waves with cc standing for the complex conjugate.
Following the procedure described in the previous subsection, it can be shown that direct and reflected terms in (31) can be calculated through the free-space propagation scheme in (7) and (8), with z = D and z = 2L − D respectively, while the interference terms lead to with a modified Green integral operator For the class of statistically quasi-homogeneous sources, we may again expand the exponent in (33) in a Taylor series in q, and retain only terms up to first order. This results in a Frobenius-Perron approximation of the interference terms, leading to a phase-factor of the optical length ∆ besides the Dirac's delta in (11). Adopting the same linear approximation for each term in (31) gives the updated WDF Similar expressions have been found in quantum mechanics [39] and optics [5] for two overlapping wave-functions. Again, the propagated CF can be obtained by the inverse Fourier transform (9) of (31) or (34), the latter being closely related to the free-space VCZT.

B. Numerical results
We chose again an initial correlation density distributed according to the Gauss-Schell model, Eq. (17), with corresponding source WDF shown in Fig. 1. We work as usual at a frequency of operation of 1 GHz corresponding to λ = 0.3 m and choose σ x = 1.0 m, σ s = 0.1 m.
We further suppose a metallic mirror at L = 1.8 m (6λ). The propagation of the intensity from the source to the mirror can be found by evolving the source WDF with the exact rule composed of Eqs. (7), (8) and those for the interference terms, Eqs. (31) - (33), and then inverse Fourier transforming the propagated WDF according to Eq. (9). The coherent energy I z (x) reaching the scan plane at z = D is given by Eq. (10), that is, by considering the correlation function at s = 0. Figure 11 shows the behavior of the intensity I z (x = 0) near the mirror, from D = 1.0 m to D = 1.8 m. The solid black line is computed through the full Green's integral operators (31) and (33), while the dashed red line is obtained by the Frobenius-Perron approximation (34). The oscillatory behaviour in (34) is due to the interference terms in the WDF.
In Figs. 12 and 13, we show the magnitude of the WDF and the associated CFs at a distance ∆ = 1.25λ (position A in Fig. 11) and at a distance ∆ = 2λ (position B in Fig.  11) from the mirror, respectively. While good agreement between the exact and the approximate propagation using the FP approximation is achieved at position A, a maximum in the correlation function, the same is not true at position B. Here the intensity is suppressed due to destructive interference and the magnitude of the CF is itself only of order O(1/k). To obtain the good agreement shown in Fig. 13, we need to take into account higher order corrections in the WDF propagator such as using the Airy function integral kernel, Eq. (15). The improvement when going from the leading order FP to the Airy function approximation is shown in Fig. 13 (b) to (c), which need to be compared with the exact WF Fig. 13 (a); the corresponding propagated CF is displayed in Fig. 13 (d). Only after going beyond the FP approximation in this way are we able to reconstruct the fine structure of the WDF. This finding is not surprising, but remarkable nevertheless; computing WDFs in a multi-scattering environment will encounter exactly these problems and we have shown that the Airy-function approximation -still faster than a full WDF propagation -can handle interference corrections successfully. We note that these corrections have been reported also in the "diffusive" Green function presented in [33].

VI. CONCLUSION
An exact propagator has been derived for field-field correlation functions of complex sources. It has been applied to a problem mimicking EM radiation from a complex source; extending this to other wave problems such as in vibro-acoustics or quantum mechanics is straightforward. The phase-space representation based on the Wigner function provides a useful means of physically interpreting the propagated data. It also serves as a very efficient computational technique both for an exact propagation of CFs and in terms of a ray approximation leading to the Frobenius-Perron operator. This provides a good description of the propagated data even when applied to source data that are relatively far from homogeneity. Where necessary, more heterogeneous sources can be ac-  counted for by higher-order approximations leading to an Airy propagator. This propagator proved important in the case of a planar random source emitting in presence of a planar reflector, for which we are able to reconstruct the fine structure of the phase space in presence of interference. Evanescent decay into the near field can also be accounted for using simple propagation rules. These rules has been used to investigate the effect of evanescent waves in near-field correlation functions. For source correlations exhibiting smaller-than-wavelength scales, we predicted a rapid initial increase of the correlation length (with distance from the source), before it saturates with the onset of the Van Cittert-Zernike behaviour at a distance of a wavelength. The approximations used have been validated through full-wave simulations using model sources and numerical sources exhibiting strong statistical inhomogeneities.