# Simulation of diffraction and scattering using the Wigner distribution function

Emilie Pietersoone, Jean Michel Létang, Simon Rit, Max Langer

<figure><img src="image_1.png" /></figure>

To cite this version:

Emilie Pietersoone, Jean Michel Létang, Simon Rit, Max Langer. Simulation of diffraction and scattering using the Wigner distribution function. Optics Letters, 2024, 49 (19), pp.5431. 10.1364/ol.523608. hal-04704935

HAL Id: hal-04704935

https://hal.science/hal-04704935v1

Submitted on 21 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Simulation of diffraction and scattering using the Wigner distribution function

EMILIE PIETERSOONE,<fnref n="1" /> ID JEAN MICHEL LÉTANG,<fnref n="2" /> ID SIMON RIT,<fnref n="2" /> ID AND MAX LANGER<fnref n="1" />,*, ID

Received 19 March 2024; revised 10 April 2024; accepted 23 April 2024; posted 23 April 2024; published 20 September 2024

We present a new, to the best of our knowledge, method to simulate diffraction images accounting for both coherent and incoherent effects, based on the Wigner distribution function of the exit wave. This permits the simulation of wave and particle effects simultaneously and simulates images photon by photon. It is motivated by artifacts observed in x ray phase-contrast images after phase retrieval, present as noise in the low spatial frequency range, which can make analysis of such images challenging. Classical simulations have so far not been able to reproduce these artifacts. We hypothesize that these artifacts are due to incoherent scatter present in the images, hence the interest in developing a simulator that permits the simulation of both diffraction and incoherent scattering. Here, we give a first demonstration of the method by simulating the Gaussian double-slit experiment. We demonstrate the capability of combining diffraction and incoherent scattering, as well as simulating images for any propagation distance.

Published by Optica Publishing Group under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

https://doi.org/10.1364/OL.523608

Conventional x ray imaging relies on the attenuation of x rays for contrast. X ray phase-contrast imaging (XPCI), on the other hand, achieves contrast from the phase shift of the beam by the sample. The simplest way to achieve phase contrast is to let the beam propagate in free space after interaction with the sample, which is known as in-line phase contrast or propagation-based imaging. The main interest in XPCI is the increased contrast in soft materials, or between materials with similar electron densities. In-line phase contrast also arises in high-resolution imaging using x ray optics [1,2]. The main inconvenience is that the phase shift is not directly accessible in the detector plane. The phase contrast contains entangled information from both the attenuation and the phase shift. Phase retrieval is a process in which the phase (and potentially the attenuation) must be reconstructed from the images.

One problem that persists in phase retrieval is the occurrence of noise in the low spatial frequency range, giving rise to cloud-like artifacts in the reconstructed images. Several algorithms

have been proposed to address this problem on the reconstruction side [3–6], but a clear understanding of the origin of these artifacts has not been established. While it is well known that phase contrast is less sensitive to low spatial frequencies in the phase shift, thus making phase retrieval sensitive to noise in this range, this does not fully explain the qualitative, cloud-like appearance of the noise in the reconstructed images [4,7]. Indeed, current simulation techniques are unable to reproduce these kinds of artifacts. Instead, reconstructing the phase from simulated images yields a noise with the appearance of low-pass-filtered noise [8–11]. Here we hypothesize that the observed low-frequency artifact stems from incoherent scattering in the sample. The aim is therefore to combine phase contrast, due to diffraction, a wave effect, and incoherent scattering, a particle effect, into one simulation.

To address this problem, we propose a new method to simulate diffraction and incoherent scattering simultaneously. It is based on the Wigner distribution function (WDF) of the wave exiting the sample. As an initial demonstration, we simulate a Gaussian double-slit experiment, as well as a modified double-slit experiment that includes scattering (Fig. 1). We demonstrate that the method works for both near and far-field.

The usual model when the propagation distance D is relatively short is the Fresnel diffraction. If the exit wave from an object is given as $\Psi(x) = \exp[-B(x)] \exp[i\varphi(x)]$, the contrast $I_D(x)$ on the detector can be written as

$$I_D(x) \approx \left| \int \Psi(x) \exp \left( i \frac{\pi}{\lambda D} |x|^2 \right) dx \right|^2 \quad (1)$$

When the propagation distance is longer than a certain distance (Fraunhofer distance), the contrast is modeled by the Fraunhofer diffraction, which is essentially the amplitude of the Fourier transform of the exit wave:

$$I_D(x) \propto \left| \mathcal{F}\{\Psi\} \left( \frac{x}{\lambda D} \right) \right|^2 \quad (2)$$

Depending on the situation, either Eq. (1) or Eq. (2) is implemented numerically. Noise is added a posteriori to the ideal diffraction images. While this approach is simple to implement, it cannot account for the effects that we want to study: reflection, refraction, and incoherent scattering.

Monte Carlo (MC) imaging simulation methods simulate images in a photon-by-photon manner by generating trajectories through an object. The simulation tracks each photon by

<footnote n="1">¹Univ. Grenoble Alpes, CNRS, UMR 5225 VetAgro Sup, Grenoble INP, TIMC, Grenoble F-38000, France</footnote>

<footnote n="2">²INSA Lyon, Université Claude Bernard Lyon 1, CNRS, Inserm, CREATIS UMR 5220 U1294, Lyon F-69373, France</footnote>

<footnote>*max.langer@univ-grenoble-alpes.fr</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Numerical experiments and their corresponding Wigner distribution function (WDF). (a) Double-slit experiment is modeled by two Gaussians. (b) Analytic WDF and its marginals, corresponding to the intensity in the spatial and angular domains, respectively. (c) Modified double-slit experiment with a scattering element in one of the slits. (d) Corresponding WDF with the amplitude in the slit with the scatterer set to 0.35.</figcaption></figure>

first defining the photon's initial properties and then propagating it through the material until a scattering or absorption event is stochastically determined, based on material properties. Depending on the type of scattering event, the photon's energy and direction are updated accordingly. This procedure is continued until the photon is absorbed, exits the sample, or the number of scattering events reaches a set limit.

Several approaches combining phase contrast and incoherent scattering have recently been proposed. These are reviewed in [10–15]. They implement either a two-part or a combination of wave optics for propagation and MC simulation for scattering. Bartl et al. [16] introduced a method that combines a simulation for diffraction and a simulation for scattering. The behavior of x ray photons are defined as both particles and waves, and it introduces realistic noise. Peter et al. [17] introduced a simulation framework that accounts for both the particle and wave characteristics of x rays in a split approach. It involves a three-step simulation, rendering the framework flexible and suitable for various methods of phase-sensitive x ray imaging. Tessarini et al. [18] introduced a simulation framework, which relies on semiclassical properties of particles, enabling the integration of both particle and wave-like effects. Langer et al. [19] incorporated an MC process to simulate the refraction, total reflection, absorption, and scattering of x rays alongside wave optics within the GATE framework.

Here, we seek a way to account for both coherent and incoherent effects in the same simulation. To do this, we propose to use

the WDF to calculate the interference in the object plane rather than in the detector plane, as is done in the classical methods. This allows to convert the exit wave into position and momentum and thus a particle-type process that can be combined with incoherent scattering.

The WDF has previously been applied to optics to different ends [20,21]. For our purposes, the WDF can be written as

$$W_{\Psi}(x, \theta) = \int \Psi(x+x') \Psi^*(x-x') \exp(-i2\pi x'\theta) dx'. \quad (3)$$

The WDF is a quasi-probability density function, meaning it represents the joint probability distribution of positions x and angles θ, but it can take negative values, representing destructive interference in the diffracted wave. This is clearly seen in the WDF of two Gaussians (Fig. 1) where the WDF has an oscillatory interference component between the two Gaussians. The marginal along x corresponds to the intensity of the Fourier transform of Ψ(x),

$$|\mathcal{F}\{\Psi\}(\theta)|^2 = \int W_{\Psi}(x, \theta)dx, \quad (4)$$

which is the intensity past the Fraunhofer distance. The marginal along θ gives the intensity at the exit of the slits

$$I_0(x) = |\Psi(x)|^2 = \int W_{\Psi}(x, \theta)d\theta. \quad (5)$$

Inspired by the negative particle formulation of quantum physics [22], the idea is then to generate particle trajectories from the WDF by generating positions and corresponding angles probabilistically. To do this, we first need the probability of finding a photon in the exit plane taking into account the interference term. Since, in the projection over θ, the positive values will cancel out the negative [23], we instead project the modulus of the WDF:

$$S_{\Psi}(x) = \int |W_{\Psi}(x, \theta)| d\theta. \quad (6)$$

While the interference term is calculated in the object plane, the sign of the generated photon must be kept so that positive and negative photons are allowed to cancel each other out in the detector plane [22].

To account for scattering in a simple way, we model it by introducing a block of material of a certain thickness and scattering coefficient, which is set to yield a certain scattering probability p_s within the support of the block and is for this demonstration considered thin (one scattering interaction at most, and all scattering occurs in the plane of the slits). If a particle is scattered, its scattering angle is calculated by probability sampling from the differential cross section (DCS) of the chosen material. In a more complete implementation, the scattering would be simulated using the classical Monte Carlo particle transport.

Based on this, we propose the following algorithm: A position x_n is drawn at random with uniform probability. A random value ε_s ∈ [0, 1] is drawn, and if ε_s < p_s, the photon is considered to be scattered. The scattering angle is generated from the DCS, the particle is then ray-traced in a straight line to the detector, and one count is registered in the corresponding detector element. If the particle is not scattered, another random number ε_d is generated. If ε_d < S_Ψ(x_n), we consider that the photon in x_n is diffracted; otherwise, it is considered absorbed or otherwise annihilated. We then draw a diffraction angle θ_n by sampling from the probability distribution P_θ(x_n) = |W_Ψ(x_n, θ)| / Σ |W_Ψ(x_n, θ)|, and the
<figure><img src="image_4.png" /><figcaption>Fig. 2. Spatial probability with interference [S<sub>Ψ</sub>(x), Eq. (6)] gives the probability of finding a photon in a certain point in the slit plane taking into account the interference term. (a) Equal amplitudes. (b) Right slit amplitude set to 0.35 as in Fig. 1(b).</figcaption></figure>

photon is ray-traced in a straight line to the detector. The sign of the particle is retrieved from W<sub>Ψ</sub>(x<sub>n</sub>, θ<sub>n</sub>). If the sign is negative, a negative potential is incremented to account for future destructive interference. If the sign is positive, it either decrements the negative potential in that pixel or, if the potential was zero, registers a count on the detector. Note that scattered photons are not compared with the negative potential and are therefore always counted as hits on the detector.

As a first demonstration of the proposed method, a Gaussian double-slit (GDS) experiment is simulated. This is chosen since the GDS has an analytical solution for the WDF, which for now enables us to bypass the most challenging numerical difficulties in this initial demonstration. Using Gaussian functions reduces the oscillatory nature of the interference term. This is done mainly for clarity; the demonstration works also for the standard double slit. To include scattering, a scattering element is placed in one of the slits [Fig. 1(c)], with uniform scattering probability. The exit wave from the GDS can be written as [24,25]

$$ \Psi_{\text{GDS}}(x) = B_1 e^{-\frac{(x-a)^2}{\Lambda^2}} + B_2 e^{-\frac{(x+a)^2}{\Lambda^2}}, \quad (7) $$

where we have added factors B<sub>1</sub> and B<sub>2</sub> to account for the attenuation caused by the scattering element, a is the position of the Gaussians (symmetrically around 0), and A their standard deviation. The WDF of Ψ<sub>GDS</sub>(x) can then be written as

$$ W_{\Psi_{\text{GDS}}}(x, \theta) = B_1 e^{-\frac{(x-a)^2}{\Lambda^2}} e^{-\frac{(a)^2}{\Lambda^2}A^2} + B_2 e^{-\frac{(x+a)^2}{\Lambda^2}} e^{-\frac{(a)^2}{\Lambda^2}A^2} + B_1 B_2 2e^{-\frac{a^2}{\Lambda^2}} e^{-\frac{(a)^2}{\Lambda^2}A^2} \cos\left(\frac{2a\theta}{\lambda}\right), \quad (8) $$

where we have omitted the scaling factor that would cause the function to be a probability density (this aspect is actually never used here). This function is shown in Fig. 1 in two configurations along with its marginals.

In this experiment, the slit positions were set to a = 0.6 µm, the standard deviation to A = 60 nm, and the wavelength to λ = 1 nm. The detector had the size 20 µm and had 601 pixels, and the propagation distance was 250 µm.

First, the simulation was run without a scattering element. S<sub>Ψ</sub>(x) for this configuration is shown in Fig. 2(a). Results for increasing the number of photons are shown in Fig. 3. Photon-by-photon buildup of the diffraction pattern is shown in Visualization 1.

Further, this simulation was run over a range of propagation distances to demonstrate the capability to simulate the diffraction pattern regardless of the propagation distance. The evolution of the intensity with distance is shown in Fig. 4. As expected, the contrast evolves in the near-field; then at a certain distance, the

<figure><img src="image_5.png" /><figcaption>Fig. 3. Evolution of the intensity with 601 propagation distances over D ∈ [1 µm, 250 µm] recorded on a detector of size 20 µm. As expected, close to the slits, no interference is visible. With sufficient propagation, interference fringes become visible and evolve with the propagation distance (near-field). Beyond a certain distance, the pattern stops evolving qualitatively and only spreads more in space.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Simulations of the GDS for different photon counts at D = 250 µm: (a) 10<sup>3</sup>, (b) 10<sup>4</sup>, (c) 10<sup>5</sup>, (d) 10<sup>6</sup> photons.</figcaption></figure>

far-field pattern emerges and remains stationary except for a scaling.

To test the combination of coherent and incoherent effects, a scattering element was introduced in one of the slits. To account for the attenuation due to the scattering, the amplitude B<sub>2</sub> in Eqs. (7) and (8) was adjusted accordingly. The effect on the WDF and its marginals is shown in Fig. 1(b), and the effect on S<sub>Ψ</sub>(x) is shown in Fig. 2. To calculate scattering coefficients and scattering angles, we used the Rayleigh cross section for Si (Z = 14, ρ = 2.393 g.cm<sup>-3</sup>). The thickness of the material was varied to yield different scattering probabilities. The DCS is used as the scattering angular distribution and is sampled
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Simulations of the modified GDS framework for $10^6$ photons at $D = 250$ µm for different scattering probabilities: (a) 25% (theoretical double-slit pattern shown in dashed green), (b) 50%, (c) 75%, (d) 95% (theoretical single-slit pattern shown in dashed green).</figcaption></figure>

in $\theta \in [-\frac{\pi}{2}, \frac{\pi}{2}]$ using 9001 sampling points. Simulation with different scattering probabilities showed that when increasing the attenuation in one slit due to scattering, the diffraction pattern approaches the single-slit pattern, as expected (Fig. 5).

To conclude, we presented a new framework to simulate diffraction and scattering simultaneously, based on the WDF and the signed particle formulation of quantum physics. The intended application is to simulate XPCI with the purpose of recreating artifacts in experimental images that are currently not possible to simulate with existing methods. The proposed approach offers several interesting aspects, namely, that it can combine the simulation of coherent and incoherent effects simultaneously, it is independent of the propagation distance, and it builds up the images photon by photon. While the original motivation for this work was to combine coherent and incoherent effects in one simulation, the last aspect is especially interesting in view of recent developments of photon counting detectors. Recently, the double-slit experiment has been demonstrated with x rays using a photon counting detector [26], giving further weight to this argument.

The main limitation with this initial demonstration is that the WDF is known analytically. Calculating the WDF numerically is notoriously difficult. This is due to its large size, effectively doubling the dimensionality of the signal, its generally highly oscillatory nature [21], and the difficulties to correctly discretize it [27]. This choice was made to permit this demonstration while temporarily avoiding these difficulties. The method has several redeeming features, however making its numerical and practical implementation seem feasible. The WDF never needs to be calculated in its entirety. It can be calculated separately point by point in space, first to generate $S_\psi(x)$ once at the beginning of the simulation and then to generate $P_\theta(x_n)$ for each particle. This can be done in a parallel manner for each particle, e.g., as a part of a particle transport code. Imposing regularity on the exit wave, or localizing interference and the sample accordingly, could control the oscillatory nature of the WDF.

While the motivation for this work is to precisely simulate XPCI, taking into account both coherent and incoherent effects to address specific artifacts, the current demonstration relies

on the analytic calculation of the WDF, which is only possible in certain defined cases. While the first numerical results are encouraging, future work will focus on the appropriate numerical implementation of the WDF in this context to permit a practical implementation of an x ray phase-contrast simulator, as well as experimental validation of the method. Once these difficulties are overcome, the photon-by-photon aspect of the method has the potential to yield realistic simulations in low-dose XPCI scenarios. The simultaneous simulation of scattering and phase contrast would allow addressing more precisely the low-frequency noise problem in phase retrieval, e.g., through providing data for machine learning-based algorithms.

**Funding.** Agence Nationale de la Recherche (ANR-11-LABX-0004, ANR-11-LABX-0063, ANR-19-P3IA-0003).

**Acknowledgment.** We thank Emmanuel Brun, STROBE, Université Grenoble Alpes, for fruitful discussion and suggestions.

**Disclosures.** The authors declare no conflicts of interest.

**Data availability.** Data underlying the results presented in this paper are available in Visualization 1.

## REFERENCES

1. R. Mokso, P. Cloetens, E. Maire, et al., Appl. Phys. Lett. 90, 144104 (2007).

2. M. Langer, A. Pacureanu, H. Suhonen, et al., PLoS One 7, e35691 (2012).

3. D. Paganin, S. C. Mayo, T. E. Gureyev, et al., J. Microsc. 206, 33 (2002).

4. M. Langer, P. Cloetens, and F. Peyrin, IEEE Trans. on Image Process. 19, 2428 (2010).

5. M. Langer, P. Cloetens, B. Hesse, et al., Phil. Trans. R. Soc. A. 372, 20130129 (2014).

6. M. Langer, P. Cloetens, A. Pacureanu, et al., Opt. Lett. 37, 2151 (2012).

7. J.-P. Suuronen, B. Hesse, M. Langer, et al., J. Synchrotron Radiat. 29, 843 (2022).

8. M. Langer, P. Cloetens, J.-P. Guigay, et al., Med. Phys. 35, 4556 (2008).

9. D. Paganin, A. Barty, P. McMahon, et al., J. Microsc. 214, 51 (2004).

10. E. Pietersoone, J. M. Létang, S. Rit, et al., Instruments 8, 1 (2024).

11. L. Quénot, E. Brun, J. M. Létang, et al., Phys. Med. Biol. 66, 175027 (2021).

12. G.-H. Chen, J. Zambelli, N. Bevins, et al., Current Medical Imaging 6, 90 (2010).

13. A. Bravin, P. Coan, and P. Suortti, Phys. Med. Biol. 58, R1 (2013).

14. S. Berujon, R. Cojocaru, P. Piault, et al., J. Synchrotron Radiat. 27, 284 (2020).

15. L. Quénot, S. Bohic, and E. Brun, Appl. Sci. 12, 9539 (2022).

16. P. Bartl, J. Durst, W. Haas, et al., in *2009 IEEE Nuclear Science Symposium Conference Record (NSS/MIC)* (2009), pp. 3577–3580.

17. S. Peter, P. Modregger, M. K. Fix, et al., J. Synchrotron Radiat. 21, 613 (2014).

18. S. Tessarini, M. K. Fix, P. Manser, et al., Sci. Rep. 12, 2485 (2022).

19. M. Langer, Z. Cen, S. Rit, et al., Opt. Express 28, 14522 (2020).

20. M. Mout, M. Wick, F. Bociort, et al., Opt. Eng. 57, 014106 (2018).

21. M. A. Alonso, Adv. Opt. Photonics 3, 272 (2011).

22. J. M. Sellier, J. Comput. Phys. 297, 254 (2015).

23. M. S. Bartlett, Math. Proc. Cambridge Philos. Soc. 41, 71 (1945).

24. T. Cuypers, R. Horstmeyer, S. B. Oh, et al., in *2011 IEEE International Conference on Computational Photography (ICCP)* (2011), pp. 1–9.

25. J. Cervenka, R. Kosik, and M. Nedjalkov, J. Comput. Electron. 20, 2104 (2021).

26. T. Gureyev, C. Hall, B. Arhatari, et al., “Young double-slit interference with single hard x-ray photons,” arXiv, arXiv:2402.07377 (2024).

27. G. Björk, A. B. Klimov, and L. L. Sánchez-Soto, *Chapter 7: The Discrete Wigner Function*, Vol. 51 of Progress in Optics (Elsevier, 2008), pp. 469–516.