Propagation of partially coherent radiation using Wigner functions
Boaz Nash, Nicholas Goldring, Jonathan Edelen, Stephen Webb, Rafael Celestre
PPropagation of partially coherent radiation using Wigner functions
Boaz Nash, ∗ Nicholas Goldring, Jonathan Edelen, and Stephen Webb
RadiaSoft LLC, 3380 Mitchell Lane, Boulder CO 80301
Rafael Celestre
ESRF - The European Synchrotron (Dated: September 18, 2020)Undulator radiation from synchrotron light sources must be transported down a beamline fromthe source to the sample. A partially coherent photon beam may be represented in phase space usinga Wigner function, and its transport may use some similar techniques as is familiar in particle beamtransport. We describe this process in the case that the beamline is composed of linear focusing anddefocusing sections as well as apertures. We present a compact representation of the beamline mapinvolving linear transformations and convolutions. We create a 1:1 imaging system (4f system) witha single slit on the image plane and observe the radiation downstream to it. We propagate a Gaussianbeam and undulator radiation down this sample beamline, drawing parameters from current andfuture ultra low emittance light sources. We derive an analytic expression for the partially coherentGaussian case including passage through a single slit aperture. We benchmark the Wigner functioncalculation against the analytical expression and a partially coherent calculation in the SynchrotronRadiation Workshop (SRW) code.
I. INTRODUCTION
Synchrotron radiation sources, either storage ring orFEL-based, require optical beamlines to transport theradiation to the experimental sample. As performanceof these sources is being pushed to lower emittance andhigher coherence, the beamline performance and model-ing must be accordingly improved [1].Models of radiation transport through the beamlineelements exist in a hierarchy of levels of accuracy andcomplexity. At the simplest level, one can use analyticalformulae to propagate beam sizes, divergences and co-herence lengths through an idealized beamline (see [2]).At the next level is the geometric optics description us-ing a ray-tracing approach (e.g. SHADOW [3]). At ahigher level of complexity, one may use a physical opticsapproach which requires wavefront propagation software(e.g. SRW [4]). The wavefront propagation allows theinclusion of diffraction effects in coherent optics, but ismore computationally intensive than the ray-tracing ap-proach.To go beyond fully coherent optics, accurate represen-tation of radiation requires the model include partial co-herence resulting from considerations of statistical optics[5, 6]. For synchrotron radiation, this partial coherenceresults from the finite electron beam size causing ran-domness in the phase of the emitted radiation. In thewavefront propagation method, this may be taken intoaccount by propagating multiple initial wavefronts, ei-ther via a sampling of the phase space of initial electronbeam [7] or by a coherent mode decomposition [8–10].Another approach to treating partially coherent radi-ation involves the use of Wigner functions. The Wigner ∗ [email protected]; function formalism for synchrotron radiation was pio-neered by K.-J. Kim [11]. The Wigner function was orig-inally developed in quantum statistical mechanics [12–14], as an alternative phase-space representation to thedensity operator. In an optical context, the Wigner func-tion may represent the types of systems described in sta-tistical optics, where multiple wavefronts are simultane-ously present with random phase relations between them.The properties of Wigner functions and the relation be-tween the quantum mechanics and optics context are de-scribed by Bazarov et al. [15].Although Wigner functions for fully and partially co-herent synchrotron radiation have been computed, theyhave not been widely used in the propagation downbeamlines. In this paper, we demonstrate propagation ofthe fully and partially coherent Wigner functions througha simplified beamline. We compute beamline maps thatmay be applied to any initial condition Wigner Functionwithout having to be recomputed. We limit ourselves tolinear maps, (so-called ABCD matrix, in the optics lit-erature [16, 17]) under which the Wigner function trans-forms in a straightforward manner. We will also considerphysical apertures, focusing on the case of single slits,an important element in most x-ray beamlines. We as-sume separable radiation such that we may work with2D Wigner functions [15, 18]. The extension to higherdimensionality is straightforward. Including non-linearelements, such as spherical aberrations in a lens, will bea topic of future work.As the size of the electron beam increases, the Wignerfunction becomes dominated by the Gaussian electronbeam and the coherence decreases. To understand thistransition, we consider Gaussian initial Wigner functions.We are able to derive an analytic expression for thediffraction of a partially coherent Gaussian through a sin-gle slit. We may thus validate our algorithms for lineartransport and passage through an aperture. In addition, a r X i v : . [ phy s i c s . acc - ph ] S e p we may compare the analytical Gaussian result to thecase of undulator radiation and gain understanding asto when the Gaussian result may be adequate in simula-tion. Finally, we validate our undulator radiation trans-port using our simplified Wigner function map methodto a partially coherent SRW calculation. II. X-RAY BEAMLINE MODELING
The goal of X-ray beamline modeling is to take aninitial radiation field and propagate that field througha complex series of elements to accurately predict whatthe wavefront will look like when it reaches the sample.There are several methods of approaching this: we surveyphysical optics in the Appendix A, for example. We startwith a summary of the properties of the radiation Wignerfunction. Next, we provide a brief discussion of lineargeometric optics including apertures. Finally, we showhow to evolve the Wigner function under the action of amatrix-aperture beamline.Consider an optical ray and attribute to it a wave-length λ following a trajectory starting at position s = 0where the radiation is created and ending at S = L at theend of the beamline. This trajectory will in general notbe straight due to reflections off of mirrors and gratingsand passing through other optical elements. At each po-sition s , along the trajectory, we assign transverse phasespace coordinates, (cid:126)z , (cid:126)z = xθ x yθ y (1)where x and y are transverse coordinates, θ x and θ y arecorresponding angles with θ x = dxds and θ y = dyds . A. Radiation Wigner functions
We represent the radiation along the trajectory bymeans of a Wigner function W ( (cid:126)z ). Many of theproperties of Wigner functions have been reviewed byBazarov [15]. We mention several of them so that ourtreatment here is self-contained. First, the Wigner func-tion is normalized: (cid:90) d(cid:126)z W ( (cid:126)z ) = 1 (2) Note that we will ignore the longitudinal phase space coordi-nates in this work. A 6D phase space treatment would includethe longitudinal position and relative wavelength deviation. Inparticular we ignore pulse length effects and assume a monochro-matic beam.
We adopt this normalization for clarity of presentationand close connection to the corresponding quantum me-chanical formalism. The Wigner function is related tothe brightness (or brilliance) function B by means of anoverall factor of the radiation flux φ : B ( (cid:126)z ; s ) = φ ( s ) W ( (cid:126)z ; s ) (3)Thus, as the radiation moves along the beamline, pro-gressing in s and passes through absorbing elements (suchas apertures that we consider here), the flux φ will reduce,but the normalization of W ( (cid:126)z ; s ) will remain constant.Now, suppose we know the Wigner function at s = 0, W ( x, θ x , y, θ y ). We will assume that W is separable;that is that W obeys W ( x, θ x , y, θ y ) = W x ( x, θ x ) W y ( y, θ y ) (4)In the examples we consider, the Wigner functionremains separable throughout the beamline, and thuswe may consider propagation of the components sepa-rately. We thus refer to W x ( x, θ x ) or W y ( y, θ y ) simplyby W ( x, θ ).Now, suppose that W ( x, θ ) represents fully coherentradiation. Then, there exists an electric field E ( x ) suchthat W ( x, θ ) = 1 λ (cid:90) ∞−∞ E ∗ (cid:18) x − φ (cid:19) E (cid:18) x + φ (cid:19) e − πiλ φθ dφ. (5)We may write this equation in an operator form: W ( x, θ ) = W ( E ( x )) (6)where we refer to W as the Wigner transform operator, or W ( x, θ ) as the Wigner function associated with E ( x ). Inthe case of fully coherent radiation, the electric field maybe reconstructed from the Wigner function as follows [15] E ∗ ( x ) E (0) = 1 λ (cid:90) ∞−∞ W (cid:16) x , θ (cid:17) e πiλ xθ dθ. (7)Now, in the case where W is partially coherent, theredoes not exist a single, well-defined wavefront associatedwith the Wigner function. Rather, there exists a wholesequence of fields E j ( x ) with j = 1 . . . ∞ . The Wignerfunction is given as the (infinite) sum of the Wignertransforms associated with the E j : W ( x, θ ) = (cid:88) j W ( E j ( x )) (8)The decomposition of a given partially coherent Wignerfunction into a set of underlying modes E j is not gener-ally unique. The process of finding an astute choice of Note that for the case of undulator radiation, this condition onlyapproximately holds. On resonance, the condition is satisfied towithin a few percent over a wide range of emittance values. SeeFig. 6 in [18]. such modes is known as “coherent mode decomposition.”The modes are often taken to be orthornormal (see sec-tion 4.7 in [19]). However, in the case of synchrotron ra-diation, one may also consider the single electron modesas a form of coherent mode decomposition also satisfyingEqn. (8).From the Wigner function, one may compute the quan-tity, µ , known as the degree of coherence via the followingexpression. µ = λ (cid:90) W ( x, θ ) dxdθ (9)In the fully coherent case where W may be derivedfrom an electric field, we find µ = 1. For the partiallycoherent case, µ < B. Linear Geometric Optics
Along an X-ray beamline, there are optical elementswhich the radiation will interact with. We include thesein our model by including a varying optical path dif-ference, that is, a phase and amplitude modulation, asa function of the Cartesian coordinates ( x, y ; s ). Thisoptical path length difference is often a function of theindex of refraction n ( x, y ; s ) for transmission elements.In addition we will include physical apertures, which al-low radiation within a certain transverse region to passunimpeded, and absorb all radiation outside of that re-gion. The aperture elements may be described by trans-fer functions t ( x, y ) which describe the region where raysmay pass, and where they are absorbed. In fact, wemay allow more general aperture elements with valuesbetween 0 and 1 in which the intensity of the ray may bereduced but not fully absorbed.Now consider another ray, starting at a different ini-tial condition (cid:126)Z . The evolution of this ray down thebeamline may be described by the action of the followingHamiltonian H ( x, θ x , y, θ y ; s ) = − (cid:113) n ( x, y ; s ) − θ x − θ y . (10)with θ x and θ y playing the role of momenta, and withposition along the trajectory s as independent variable. n ( x, y ; s ) is the local index of refraction that the radia-tion is passing through. The result of this Hamiltonianformulation for geometric optics is that the offset ray willfollow Hamilton’s equations:˙ Z i = J ij ∂H∂Z j (11)with summation over the repeated index j implied andthe dot representing dds . For 4D phase space (in the casethat we ignore variation in the z and δ coordinates), the matrix J is given by J = − − (12)If we allow arbitrary index of refraction n ( x, y ) in ourmodel beamline, then the equations of motion will benon-linear. For the purposes of this paper, we will restrictto the approximation that the index of refraction variesquadratically, leading to linear equations of motion forthe ray tracing. In particular, as our model for the indexof refraction n ( x, y ; s ), we will assume it to be constantalong the optical axis and then to fall off quadraticallyin the transverse directions. Thus, we parametrize it asfollows : n ( x, y ; s ) = n ( s ) − κ x ( s ) x − κ y ( s ) y . (13)Expanding the Hamiltonian for small angles, to quadraticorder, we find H ( x, y, θ x , θ y ; s ) ≈ − n ( s ) + θ x + θ y n + κ x ( s ) x + κ y ( s ) y . (14)As mentioned, this Hamiltonian will lead to linear equa-tions of motion. The solution may thus be expressed inmatrix form as (cid:126)Z ( s ) = M ( s ) (cid:126)Z (15)One may solve these equations and produce a linear mapfor a given beamline section. All the optical elements(excluding the apertures which we deal with separatelyare thus captured in the transfer matrix M ( s ) varyingalong the beamline. Because the ray tracing dynamicsare derived from a Hamiltonian, we are assured that theresulting transfer matrix is symplectic. That is: M ( s ) T JM ( s ) = J (16)for all s along the beamline. The matrix J is given inEqn. (12).Although we have formulated this section in terms ofa Hamiltonian theory to bring out some of the formalproperties of the propagation, the transfer matrix mayalso be computed for realistic beamlines using ray tracingsoftware. See [20] for an example of this calculation fora KB mirror system.In the next section, we describe the evolution of theWigner function under the matrix M . This covers boththe fully coherent and partially coherent case. In thefully coherent case, a different formalism is possible forthe propagation: that of the linear canonical transform(LCT). We outline this in Appendix A. Note that we leave out a coupling term in the Hamiltonian pro-portional to xy so that the separable condition is satisfied. C. Partially coherent propagation with Wignerfunction
The evolution equation for the Wigner function is givenas follows [13] ∂W ( x, θ x , y, θ y ; s ) ∂s = [ W, H ] (cid:63) , (17)where the Moyal bracket is defined for arbitrary phasespace functions f and g as[ f, g ] (cid:63) = 1 iλ ( f (cid:63) g − g (cid:63) f ) (18)and the Moyal star is given by (cid:63) = e iλ (cid:16) ←− ∂ x −→ ∂ θ −←− ∂ θ −→ ∂ x (cid:17) (19)with the arrows representing action of the derivative, ei-ther to the left or right, depending on arrow orientation.Fortunately, in the case of a quadratic Hamiltonian,evolution of the Wigner function is much simpler andmore intuitive. The Moyal bracket reduces to the Poissonbracket giving classical evolution (again using the quan-tum/classical mechanics analogy). One finds that themotion in phase space is a linear transformation. Theseconsiderations allow us to formulate our approach. Inparticular, consider a beamline where the geometric op-tics is defined by a transfer matrix M acting on the phasespace vector (cid:126)z : (cid:126)z f = M(cid:126)z i (20)The Wigner function evolves along this beamline accord-ing to W f ( (cid:126)z ) = W i ( M(cid:126)z ) . (21)We may describe this transformation with the opera-tor, U M , defined as U M ( W ( (cid:126)z )) = W ( M(cid:126)z ) (22)By performing a change of variables, one may showthat the degree of coherence µ is conserved under lineartransport U M . That is, µ ( W ( (cid:126)z )) = µ ( W ( M(cid:126)z )) (23)The degree of coherence is not conserved after passingthrough an aperture, which we now describe.We would now like to consider the way in which Wignerfunctions are impacted by physical apertures. As de-scribed by Bazarov, for the electric field, the effect of theaperture is given by E s (cid:48) ( (cid:126)x ) = E s ( (cid:126)x ) t ( (cid:126)x ) (24)where t ( (cid:126)x ) is the transmission function of the aperture. In terms of the Wigner function, the action of the aper-ture is given by the partial convolution (in the angularvariable) of the aperture Wigner function: W s (cid:48) = W s ∗ θ W t ≡ A t W s (25)The aperture Wigner function W t is given by apply-ing the Wigner transform to the aperture transmissionfunction. That is, we apply Eqn (5) where the aperturetransmission function t ( x ) plays the role of the electricfield.Since an aperture results in absorption of radiation,the normalization of the Wigner function would changeafter passing through. As given in Eqn. (3), the Wignerfunction is related to the brightness by a factor of thetotal flux. For simplicity, we will ignore the changingvalue of flux along the beamline and consider Wignerfunction to be consistently normalized throughout. Thuswe normalize the transmission function according to Eqn.(B2) and thus W ( t j ) will be normalized according to Eqn.(B7). Because of the different normalization, the Wignertransform of a transmission function could be consideredas a sort of filter that acts on the incoming radiationWigner Distribution Function (WDF). Thus, part of thework of propagation of WDFs through beamlines includ-ing apertures involves the calculation of the Wigner filterfunctions. We will give an example of this function forthe single slit aperture in a later section. A substantialpart of the understanding about diffraction effects can begleaned by examination of these Wigner filter functions. D. Matrix-aperture beamlines
Consider the beamline schematic as shown in Fig. 1consisting of an undulator source with an electron beamand subsequent sections which may be described by ma-trices, M j , and apertures with transmission function, t j ( x ). An electron beam with distribution f e ( (cid:126)z ) passesthrough an undulator producing synchrotron radiation.Let E ( (cid:126)x ) be the electric field produced by a single elec-tron as it appears at the center of the undulator. We maynow construct the multi-electron Wigner function for theundulator radiation as will be described in Section IV-A(see Eqn. (44)).We consider an axial ray coming from the center ofthe undulator and proceeding along the beamline untilthe final observation plane located at position s n . Wenote that the optical axis described by this axial rayis not necessarily a straight line. In particular, mirrorswill cause angular deviations from the central trajectoryand at each point, the transverse coordinates are rela-tive to the direction of the central ray. Along the beam-line, there are apertures located at positions s through s n − which are represented by transmission functions t through t n − . For the purpose of our simplified beam-line, non-linear aberrations will be ignored and we will FIG. 1. Matrix aperture beamline schematic. assume that the transport between apertures j − j may be represented by the matrix M j .Following Eqn. (22), we find an operator for this beam-line section given by U M j that simply transports thephase space by applying the matrix M j . We have thusdefined the operator for propagation through sections oflinear transport, U M j , that may include mirrors, lenses,and other elements when remaining close to the opticalaxis. Likewise, the apertures may also be represented byoperators, A t j , as given by Eqn. (25). The operator forthe entire simplified beamline may then be given by O BL = U M n A t n − U M n − ... A t U M A t U M (26)Then the fully coherent and partially coherent simula-tions can be simply written as W se,n = O BL W se (27) W me,n = O BL W me (28)respectively. III. EXAMPLES – GAUSSIAN RADIATION
We consider Gaussian radiation (GR) and undulatorradiation (UR) propagating through a simple matrix-aperture beamline doing a 1:1 imaging of the source witha horizontal slit at the image plane and observing the ra-diaion downstream to it. A schematic for this beamlineis depicted in Figure 2 and the corresponding source andbeamline parameters are provided in Table I.We first illustrate many of the useful properties of radi-ation Wigner Function Distribution (WFD) propagationwith the example of Gaussian radiation. The one-to-oneimaging section preserves the coherence properties and This assumes that that there is no beam cropping and that theimaging system can resolve the radiation source.
TABLE I. Numerical Simulation Parameters
Undulator radiation
Length, L u λ u B k E λ Electron beam (APS-U)
Energy, E I
200 mAHorizontal emittance, (cid:15) x σ x µ mHorizontal RMS divergence, σ x (cid:48) µ radVertical emittance, (cid:15) y σ y µ mVertical RMS divergence, σ y (cid:48) µ rad Beamline parameters
Distance to lens, L f µ mFinal drift length, L d thus the first element to consider is the single slit aper-ture. Propagation through the aperture and subsequentpropagation through free space allows us to observe theimpact of decreasing coherence on the WFD and corre-sponding intensity pattern. As the divergence becomessufficiently large, the oscillations in θ in the apertureWigner filter are washed out, and coherent diffractioneffects are seen to be destroyed. The GR case has theadvantage of exploiting the analytical expressions pre-viously derived for benchmarking our numerical WFDtransport methods.We then explore the more complex case of UR prop-agating through this beamline. Changing the electronbeam emittance allows us to control the degree of par-tial coherence. At small electron beam emittances, theWFD is dominated by the single electron WFD, but asthe electron beam emittance increases, this structure be-comes less relevant and the WFD asymptotically reducesto the GR case. In order to confirm the accuracy of ourcalculations, we have set up the same beamline model FIG. 2. Simple matrix-aperture beamline example. in SRW. Due to the separability of the horizontal andvertical components, we are able to demonstrate a directcomparison between SRW and our own partially coherentWFD based propagation results.Although our WFD based methods are substantiallyless time-consuming than full multi-electron partially co-herent calculations, significant computational issues re-main present regarding memory and runtime. Due tothe complexities in the radiation and aperture WFD pat-terns, sufficient resolution is required to obtain accurateresults. Often times this runs into conflict with memoryand runtime requirements for the calculations. For thisreason, we report the grid sizes ( n x number of points inposition and n θ points in angle) used in each numericalcalculation.We first consider a Gaussian beam propagatingthrough the simple matrix aperture beamline where wecan compare our numerical results to the analytical ex-pression given by Eqn. (40). A. Analytic calculation for Gaussian radiation withsingle slit aperture
In this section, we derive an analytic expression forthe propagation of a Gaussian Wigner function througha single slit aperture of width a . We are able to findan analytic expression for the radiation immediately af-ter the slit as well as the radiation after having driftedsome distance beyond. The aperture is described by thetransmission function t ( x ) = rect (cid:16) xa (cid:17) (29)with rect( x ) = , | x | < | x | = , | x | > (30) The corresponding Wigner function is given by [21] W ss ( x, θ ) = rect (cid:16) xa (cid:17) (cid:2) πθλ ( − | x | + a ) (cid:3) θ (31) ≡ θQ ) θ (32) Q = 2 π a − | x | λ (33)Plots of the single slit aperture transmission functionand corresponding WFD are displayed in Figure 3.Let us consider a partially coherent Gaussian Wignerfunction given by W ( x, θ ) = 12 πσ x σ θ e − x σ x − θ σ θ (34)Let us write the relationship between the position andangular spreads as σ x σ θ = m λ π (35)where m is the beam quality factor known in the laseroptics literature . In the case that m = 1, the Wignerfunction represents a coherent wavefront . For m > θ variable. That is: W s (cid:48) = W s ∗ θ W ss (36)= rect (cid:16) xa (cid:17) πσ x σ θ e − x σ x I ( θ, Q ) (37)with I ( θ, Q ) = (cid:90) ∞−∞ e − ( θ − τ )22 σ θ sin( τ Q ) τ dτ (38) The beam quality factor is typically represented by a capital M . We have used a lower-case m here so as not to confuse thisquantity with the transfer matrix, M We note that for non-Gaussian wavefronts, in the coherent case,the product σ x σ θ exceeds λ π . For undulator radiation, one hasapproximately λ π . See e.g. [22–25] . FIG. 3. Single slit aperture (a) transmission function and (b) corresponding Wigner filter function.
This integral may be performed (see Appendix D for de- tails) with the result I ( θ, Q ) = πe − ˆ θ Im (cid:40) i erf (cid:32) ˆ Q + i ˆ θ √ (cid:33) (cid:41) (39)where Im represents taking the imaginary part of theargument and ˆ Q = Qσ θ and ˆ θ = θ/σ θ . Writing it all outexplicitly to see the x and θ dependence, and adding thedrift following the aperture, we have: W s (cid:48) ( x, θ ) = 1 σ x σ θ rect (cid:18) x − L d θa (cid:19) e − ( x − Ldθ )22 σ x − θ σ θ Im (cid:40) i erf (cid:32) ( a − | x − L d θ | ) σ θ λ + i θσ θ √ (cid:33) (cid:41) (40)where L d is the drift length following the single slit aper-ture.From Eqn. (9), and using the Gaussian Eqn. (34),we find the expression for degree of coherence, µ , beforepassing through the slit to be µ = 1 m (41)Figure 4 shows the WFDs resulting from the convolu-tion of GR with the aperture Wigner filter function dis-played in Figure 3 (b). This calculation was performed ona discrete grid of size n x = 4000 by n θ = 4000. Becausethe analytic expressions are available in the GR case, weare free to increase the resolution to a very fine levelwithout incurring excessively large runtimes and mem-ory demand. These results demonstrate the diminishingcoherence effects for a fixed beam size and increasing di-vergences corresponding to increasing m values.Figure 5 shows the intuitive strength of the WFDmethod by illustrating the mechanism of diffraction inwhich the oscillations in θ give rise to interference effectsfollowing the drift. Before the drift, the oscillations of the Wigner functions in θ will cancel when performing aprojection, but following the drift, they result in a largedip at the center of the intensity distribution which canbe seen in Figure 5 (b). In Figure 5 (b), results are shownfor increasing values of m , showing how this interferenceeffect disappears with decreasing coherence. In order todemonstrate the validity of our numerical computations,Figure 6 compares the spatial projections of the driftedGR WFD at m = 3 calculated numerically and analyt-ically. Because the numerical result must be redepositedonto the initial grid, one tends to truncate the calculationwithin a smaller range. The WFD was computed numer-ically by first constructing a discrete Gaussian along withdiscrete representation of the single slit aperture WFD.These two functions were then convolved in θ to computethe result following the aperture. Finally, this WFD wasdrifted for 10 cm by application of the drift transfer ma-trix and redeposition upon the original grid using themethod as first reported in [26].Figure 7 displays the dependence of degree of coher-ence on increasing m value both before and after the sin-gle slit aperture. We show agreement between numericaland analytical (Eqn. (41)) calculations for this quantity FIG. 4. WFD plots for Gaussian beam after passing through single slit aperture. Beam size, σ x , has been fixed and divergencevaries according to the parameter m in Eqn. (35). The radiation wavelength λ = 3 .
98 ˚A. (a) m = 3, σ θ = 5.90 µ rad (b) m = 5, σ θ = 9.83 µ rad (c) m = 10, σ θ = 19.66 µ rad (d) m = 15, σ θ = 29.49 µ rad (e) projection on spatial axis for cases (a) -(d). before the aperture. We note that following propaga-tion through the aperture, the degree of coherence hasincreased which is to be expected since a more coherentsubset of the radiation has been selected by the aperture.In addition, we confirm that the degree of coherence isinvariant under the final free space propagation which ispredicted by Eqn. (23). IV. EXAMPLES – UNDULATOR RADIATIONA. Partially coherent undulator radiation
As an example of transporting a partially coherentWigner function, we will consider the case of undulatorradiation resulting from a beam of electrons in a syn-chrotron light source.A single electron with initial phase space coordinates (cid:126)z passing through the undulator will produce a coherentwavefront E u ( x ). The electron will be drawn from a dis-tribution of electrons, f e ( (cid:126)z ), that are circulating in theelectron storage ring. This distribution generally takes ona Gaussian form resulting from an equilibrium betweenthe damping and diffusion effects from synchrotron ra-diation [27, 28]. Due to differences in the longitudinalcoordinates of the emitting electrons, the radiated wave-fronts will add incoherently and produce partially coher-ent radiation.We assume that the radiation will satisfy W ( (cid:126)z ; (cid:126)z ) = W ( E ( (cid:126)x ; (cid:126)z )) = W ( (cid:126)z − (cid:126)z ) (42) where W ( (cid:126)z ) = W ( E ( (cid:126)x ; (cid:126)z = (cid:126) W me may be related to the single electron Wignerfunction via convolution with the electron beam distri-bution: W me ( (cid:126)z ) = W se ( (cid:126)z ) ∗ f e ( (cid:126)z ) (44)Next we note that the convolution of the single elec-tron Wigner function with the electron beam distributionmay be postponed until the first aperture by applying thefollowing identity. U M W se ( (cid:126)z ) ∗ f e ( (cid:126)z ) = W se ( M (cid:126)z ) ∗ f e ( (cid:126)z ) (45)where the transfer matrix, M , represents the propaga-tion from the source to the first aperture. Recall the dis-cussion in section II-D for the definition of the operator U M . This identity is known as the “emittance convolu-tion theorem” attributed to K.-J. Kim (see e.g. discus-sion in [29] and references therein). This theorem allowsus to use coherent optics until the first aperture wherewe then need to construct the partially coherent Wignerfunction via convolution with the electron beam phasespace distribution that has been propagated to the sameposition via the transfer matrix of the first section M .The fact that apertures are represented by only a partialconvolution prevents us from extending this identity andapplying the convolution beyond the first aperture. FIG. 5. (a) WFD from Figure 4 for m = 3 after drifting 10 cm. (b) Spatial projections of drifted WFDs for m = 3, 5, 10and 15. FIG. 6. Comparison of analytically and numerically calculated spatial projections of drifted WFDs for m = 3.FIG. 7. Degree of coherence with varying m values before and after the aperture. As an example of UR, we consider the undulator and electron beam with parameters defined in Table I.0The software package Synchrotron Radiation Workshop(SRW) [30] is used for the initial wavefront calculationthat will be used to construct the corresponding WFD.We compute the radiation at the first optical element andprogress it through the lens and drift such that we achievethe one-to-one focusing yielding the radiation as it wouldappear at the center of the undulator. Figure 8 (a) and(b) display the real and imaginary parts of the electricfield and in (c) and (d) we have projected these fieldsonto the horizontal axis and normalized them accordingto Eqn. (B2). Note that since our planar undulator hasa vertical magnetic field, we have selected the dominantpolarization component which is horizontal.The UR WFD is constructed from the electric fieldprojections and is displayed in Figure 9 (a). Figure 9(b) displays the electron beam distribution correspond-ing to the APS Upgrade parameters given in Table I.The convolution of the single electron UR WFD and theelectron beam distribution results in the multi-electronWFD given in Figure 9 (c). In this case, we can see thatwe are near the diffraction limit as the electron beam sizeis of the same order as the UR single electron radiationbeam size. The finite emittance result shows propertiesof both the underlying UR WFD and the Gaussian elec-tron beam distribution . This WFD is then propagatedthrough the example beamline and the results from bothour WFD transport method and SRW are detailed in thefollowing section. B. Comparison with multi-electron SRWsimulation
We now report the results from transporting the URthrough the single slit aperture and subsequent drift.The multi-electron WFD is computed by convolutionwith the electron beam and propagation through the sin-gle slit aperture is performed by means of convolutionin θ with the aperture WFD. The free space drift is per-formed as done in the Gaussian case with the same lineartransport algorithm used in the preceeding section. Wenote that this transport algorithm may be applied to anyWFD.This numerical result was benchmarked against a par-tially coherent SRW calculation wherein the same ex-ample beamline was used. The multi-electron SRW cal-culation requires running the coherent calculation manytimes for different macro-electrons. In this case of 42pmemittance, we found we could achieve convergence with The SRW source code is maintained by Oleg Chubar and is ava-ialble at https://github.com/ochubar/srw We note that this result is dependent on the undulator radiationenergy used of 3.115 keV. Higher energy radiation will take up asmaller phase space footprint, and the electron beam contribu-tion will be correspondingly larger. It is thus easier to reach thediffraction limit of high coherence for lower energy radiation, fora fixed electron beam size and divergence. m value of 1.33.Recall Eqn. (35). Figure 11 provides a comparison be-tween the propagated Gaussian and Undulator radiationfor the cases of fully and partially coherent radiation. Forthe Gaussian radiation, the second moments have beentaken from adding in quadrature the undulator σ x and σ θ with those of the electron beam. We note that for thefully coherent case (single electron), there is not perfectagreement between the two calculations, with the Gaus-sain case having a larger central dip from diffraction. Inthe case of the 42 pm emittance, however, the Gaussianresult quite adequately reproduces the more complex un-dulator radiation calculation. V. CONCLUSION
We have described a new, unified method fortransporting coherent and partially coherent radiationthrough a beamline of linear optical elements and aper-tures. This approach relies on transporting the Wignerdistribution of the radiation wavefront in a manner akinto single-particle tracking in particle accelerators. In con-trast with physical optics modeling, the formalism is thesame for fully coherent and partially coherent radiation,and has the same computational time and complexity.The matrices in the matrix-aperture beamline we havedefined may be computed with ray tracing, in the gen-eral case, or analytically for simplified beamline models.We thus show how ray tracing may be unified with waveoptics, at least within this linear approximation.We have demonstrated this approach by transportingboth a gaussian wave front and an undulator radiationwavefront through drift space and a single slit aperture.In the case of the gaussian wave front, we find excellentagreement between an analytical result, the Wigner dis-tribution approach, and a the well-established SRW phys-ical optics code. For the undulator radiation, we haveexcellent agreement between the Wigner distribution ap-proach and multi-electron SRW calculations. Finally, wehave seen that for large enough emittance (even as small1
FIG. 8. Undulator radiation at first harmonic energy of 3.115 keV. (a) Real part of electric field. (b) Imaginary part of electricfield. (c) Horizontal projection of real electric field. (d) Horizontal projection of imaginary electric field.FIG. 9. Undulator radiation multi-electron WFD for varying emittance (cid:15) : (a) zero emittance (fully coherent) (b) 4 pm (c) 40pm (d) 4 nm. as 42pm horizontal emittance), the Gaussian results ap-proximate the undulator radiation to a high degree ofaccuracy.Future work will extend this technique to nonlinear op-tical elements, which will allow us to account for opticalaberrations. We have briefly described the formalism fornonlinear optical elements in Appendix C. Extension tothe non-linear case will raise new issues, for future ex-ploration. See e.g. [31] for some understanding of the potential resulting complications.
VI. ACKNOWLEDGMENTS
The authors would like to thank A. Wojdyla for valu-able discussion and for alerting us to the literature onLinear Canonical Transforms, including reference [32].This material is based upon work supported by the U.S.2
FIG. 10. (a) Fully coherent UR WFD following propagation through single slit aperture and 10 cm drift. (b) Comparisonbetween WFD method and SRW calculation for fully coherent, single electron (se), and partially coherent, multi-electron (me),UR spatial projections.FIG. 11. Spatial projections of Gaussian and undulator radiation propagated through single slit aperture and 10 cm drift viaWFD method. In the legend, se denotes single electron (fully coherent) and me denotes multi-electron (partially coherent).For the latter, the photon beam has been convolved with the APS-U electron beam given in Table I.
Department of Energy, Office of Science, Office of BasicEnergy Sciences, under Award Numbers DE-SC0020593and DE-SC0018571.
Appendix A: Coherent Wave Optics
Here we describe, for completeness, the approach usedfor propagating a radiation wavefront through a beam-line using physical optics techniques. The propagation ofthis wavefront through an optical beamline is the goal ofphysical optics software such as SRW and XRT [33, 34].In these codes, there exist a variety of numerical im-plementations of propagators to transport the wavefrontthrough free space, lenses, mirrors, gratings, and otheroptical elements.We begin with a complex electric field E ( x, y, s = 0)at the entrance of the beamline. Optical transportbeamlines in synchrotron radiation facilities can often berepresented with the scalar paraxial optics approxima- tion [35, 36]. For example, the near-field Fresnel integralfor free-space propagation can be written as a convolu-tion: E ( x, y, s ) = E ( x, y ) ∗ h ( x, y, s ) (A1)with h ( x, y, s ) = e iks iλs e i k z ( x + y ) (A2)A thin lens may be traversed via E ( x, y, s (cid:48) ) = e − i k f ( x + y ) E ( x, y, s ) (A3)with k = πλ , n is the index of refraction, and f is thefocal length [37].The combination of drifts, lenses, and focussing mir-rors together can be combined to create a symplectictransport matrix M for the geometric ray optics as givenin Eqn. (15). Knowing M ( s ), one may propagate the3wavefront through the linear beamline down the channelvia Linear Canonical Transformation (LCT) [32]. Explic-itly for 4D phase space, the transfer matrix is written inthe form M = (cid:18) A BC D (cid:19) (A4)The transformed electric field E f ( (cid:126)x ) is given by E f ( (cid:126)x ) = 1 (cid:112) det( iB ) (cid:90) e iπ(cid:126)u T M (cid:126)u E i ( (cid:126)x ) d(cid:126)x i (A5)with (cid:126)u = (cid:0) (cid:126)x f (cid:126)x i (cid:1) (A6)and M = (cid:18) DB − − B − − B − B − A (cid:19) (A7)where the subscripts f and i represent initial and final.We point out here, that in addition to use of analyt-ical expressions for beamline elements to determine theHamiltonian, and thus find the transfer matrix M ( s ), onemay also use a ray tracing code, set up the beamline, andby tracking a series of rays offset from the central trajec-tory, derive the transport matrix numerically along thebeamline.The effect of the apertures, represented by the transferfunctions t j ( x, y ) on the wavefront simply by multiplica-tion: E ( x, y ; s (cid:48) ) = t j ( x, y ) E ( x, y ; s ) (A8) Appendix B: Normalization of wavefronts andWigner functions
In this paper, we will assume electric fields which sat-isfy the separability condition E ( x, y ; s ) = E E x ( x ; s ) E y ( y ; s ) , (B1)where E is a constant with units of electric field.We normalize the separate electric field components in1-D such that (cid:90) ∞−∞ E ∗ ( x ) E ( x ) dx = 1 , (B2) (cid:90) ∞−∞ E ∗ ( θ ) E ( θ ) dθ = 1 . (B3)Following Bazarov [15], we have normalized the elec-tric field in the same way as wave functions are normal-ized in quantum mechanics. The second moments of thefield distribution in coordinate and angular representa-tions may now be calculated as < x > = (cid:90) ∞−∞ x E ∗ ( x ) E ( x ) dx, (B4) < θ > = (cid:90) ∞−∞ θ E ∗ ( θ ) E ( θ ) dθ. (B5) We now introduce the Wigner function defined fromthe electric field, E ( x ), as follows W ( x, θ ) = 1 λ (cid:90) ∞−∞ E ∗ (cid:18) x − φ (cid:19) E (cid:18) x + φ (cid:19) e − πiλ φθ dφ, (B6)where W ( x, θ ) will be normalized as (cid:90) ∞−∞ (cid:90) ∞−∞ W ( x, θ ) dxdθ = 1 . (B7)The Wigner function can be thought of as a probabilitydistribution in phase space except for the fact that it maybecome negative. The second moments are given simplyas < x > = (cid:90) ∞−∞ (cid:90) ∞−∞ x W ( x, θ ) dxdθ, (B8) < θ > = (cid:90) ∞−∞ (cid:90) ∞−∞ θ W ( x, θ ) dxdθ, (B9) < xθ > = (cid:90) ∞−∞ (cid:90) ∞−∞ xθW ( x, θ ) dxdθ. (B10) Appendix C: Transport of Wigner function undernon-linear maps
For propagation of the Wigner function, Bazarov(Property 5) uses the Hamiltonian H = ˆ p m + V ( x ) forthe Quantum Mechanics case of a particle in a potential.We have used a different Hamiltonian for the optics case.The problem of general transport of the Wigner functionremains. The evolution equation for the Wigner functionunder a general Hamiltonian is given as follows [13] ∂W ( x, θ x , y, θ y ; s ) ∂s = [ W, H ] ∗ , (C1)where the Moyal bracket is defined for arbitrary phasespace functions f and g as[ f, g ] ∗ = 1 iλ ( f ∗ g − g ∗ f ) (C2)and the Moyal star is given by ∗ = e iλ (cid:16) ←− ∂ x −→ ∂ θ −←− ∂ θ −→ ∂ x (cid:17) (C3)with the arrows representing action of the derivative, ei-ther to the left or right, depending on arrow orientation. Appendix D: Analytical Calculation of a GaussianWigner Distribution Through a Single Slit
The passage of a Gaussian Wigner distribution througha slit can be evaluated in terms of the integral I ( x, θ ) = (cid:90) ∞−∞ dτ sin τ qτ exp (cid:26) − ( θ − τ ) σ θ (cid:27) (D1)4It is convenient to normalize all the variables to σ θ , sothat ˆ θ = θ/σ θ , ˆ τ = τ /σ θ , and ˆ q = σ θ q , so that theintegral becomes: I ( x, θ ) = (cid:90) ∞−∞ d ˆ τ sin ˆ τ ˆ q ˆ τ exp (cid:26) −
12 (ˆ θ − ˆ τ ) (cid:27) (D2)This integral can be rewritten as the imaginary part ofan indefinite integral with respect to q , I ( x, θ ) = Im (cid:20) i (cid:90) d ˆ q (cid:90) ∞−∞ d ˆ τ e i ˆ τ ˆ q exp (cid:26) −
12 (ˆ θ − ˆ τ ) (cid:27)(cid:21) (D3)Expanding the Gaussian argument brings a Gaussian θ envelope out front: I ( x, θ ) = e − ˆ θ Im (cid:20) i (cid:90) d ˆ q (cid:90) ∞−∞ d ˆ τ exp (cid:26) −
12 (ˆ τ + 2 i ˆ τ (ˆ q + i ˆ θ )) (cid:27)(cid:21) . (D4) The argument can be simplified by completing thesquare, noting thatˆ τ + 2(ˆ θ + i ˆ q )ˆ τ = (ˆ τ + ˆ θ + i ˆ q ) − (ˆ θ + i ˆ q ) (D5)which then gives the integral as I ( x, θ ) = e − ˆ θ Im (cid:20) i (cid:90) d ˆ qe (ˆ θ + i ˆ q ) (cid:90) ∞−∞ d ˆ τ exp (cid:26) −
12 (ˆ τ + i ˆ q + ˆ θ ) (cid:27)(cid:124) (cid:123)(cid:122) (cid:125) √ π (cid:21) (D6)and the integral becomes I ( x, θ ) = √ πe − ˆ θ Im (cid:20) i (cid:90) d ˆ q e (ˆ θ + i ˆ q ) (cid:21) (D7)and which is then given by I ( x, θ ) = (cid:114) π (cid:20) i erf (cid:32) ˆ q + i ˆ θ √ (cid:33)(cid:21) (D8) [1] M. Sanchez del Rio, R. Celestre, M. Glass, G. Pirro, J. R.Herrera, R. Barrett, J. C. da Silva, P. Cloetens, X. Shi,and L. Rebuffi, J. Synchrotron Radiat. , 1887 (2019).[2] X. Shi, R. Reininger, R. Harder, and D. Haeffner, in Ad-vances in Computational Methods for X-Ray Optics IV ,Vol. 10388 (International Society for Optics and Photon-ics, 2017) p. 103880C.[3] M. S. del Rio, N. Canestrari, F. Jiang, and F. Cerrina,Journal of synchrotron radiation , 708 (2011).[4] O. Chubar and P. Elleaume, Particle accelerator. Pro-ceedings, 6th European conference, EPAC’98, Stock-holm, Sweden, June 22-26, 1998. Vol. 1-3 , Conf. Proc.
C980622 , 1177 (1998), [,1177(1998)].[5] J. W. Goodman,
Statistical Optics, 2nd Edition (Wiley,2015).[6] G. Geloni, E. Saldin, E. Schneidmiller, and M. Yurkov,Nucl. Instrum. Methods Phys. Res., Sect. A , 463(2008).[7] O. Chubar, L. Berman, Y. S. Chu, A. Fluerasu, S. Hul-bert, M. Idir, K. Kaznatcheev, D. Shapiro, Q. Shen, andJ. Baltser, Advances in Computational Methods for X-Ray Optics II , 814107 (2011).[8] M. Glass and M. S. del Rio, EPL , 34004 (2017).[9] R. Li and O. Chubar, in
Advances in ComputationalMethods for X-Ray Optics V , Vol. 11493 (InternationalSociety for Optics and Photonics, 2020) p. 114930G.[10] R. Co¨ısson and S. Marchesini, J. Synchrotron Radiat. ,263 (1997).[11] K.-J. Kim, in Insertion Devices for Synchrotron Sources ,Vol. 0582, edited by I. E. Lindau and R. O. Tatchyn,International Society for Optics and Photonics (SPIE,1986) pp. 2 – 9.[12] T. L. Curtright, D. B. Fairlie, and C. K. Zachos,
Quan-tum Mechanics in Phase Space: An Overview with Se-lected Papers (World Scientific) .[13] T. L. Curtright and C. K. Zachos, arXiv (2011),10.1142/S2251158X12000069, 1104.5269. [14] E. Wigner, Phys. Rev. , 749 (1932).[15] I. V. Bazarov, Phys. Rev. ST Accel. Beams , 050703(2012).[16] A. E. Siegman, Lasers (University Science Books, 1986).[17] C. Ferrero, D.-M. Smilgies, C. Riekel, G. Gatta, andP. Daly, Appl. Opt. , E116 (2008).[18] T. Tanaka, Phys. Rev. Spec. Top. Accel. Beams ,060702 (2014).[19] L. Mandel and E. Wolf, Optical Coherence and QuantumOptics (Cambridge University Press, 1995).[20] B. Nash, N. Goldring, J. Edelen, C. Federer, P. Moeller,and S. Webb, in
Advances in Computational Methods forX-Ray Optics V , Vol. 11493 (International Society forOptics and Photonics, 2020) p. 114930C.[21] C. J. Rom´an-Moreno, R. Ortega-Mart´ınez, andC. Flores-Arvizo, Revista Mexicana de F´ısica , 290(2003).[22] R. P. Walker, Phys. Rev. Accel. Beams , 050704(2019).[23] R. Khubbutdinov, A. P. Menushenkov, and I. A. Var-tanyants, J. Synchrotron Radiat. , 1851 (2019).[24] B. Nash, O. Chubar, N. Goldring, D. L. Bruhwiler,P. Moeller, R. Nagler, M. Rakitin, S. Gwo, D.-J. Huang,and D.-H. Wei, AIP Conf. Proc. , 060080 (2019).[25] B. Nash, O. Chubar, D. L. Bruhwiler, M. Rakitin,P. Moeller, R. Nagler, and N. Goldring, in Advancesin Laboratory-based X-Ray Sources, Optics, and Applica-tions VII , Vol. 11110 (International Society for Opticsand Photonics, 2019) p. 111100K.[26] B. Nash, J. Edelen, N. Goldring, and S. Webb, in (2019) p. TUPAG17.[27] B. Nash, J.-H. Wu, and A. W. Chao, Phys.Rev.ST Ac-cel.Beams (2006), 10.1103/PhysRevSTAB.9.032801.[28] M. Sands, Conf.Proc.C , 257 (1968).[29] M. Glass, Statistical optics for synchrotron emission : nu-merical calculation of coherent modes , Ph.D. thesis, Uni- versit´e Grenoble Alpes (2017).[30] O. Chubar, A. Fluerasu, L. Berman, K. Kaznatcheev,and L. Wiegart, J. Phys. Conf. Ser. , 162001 (2013).[31] A. J. Dragt and S. Habib, arXiv (1998), quant-ph/9806056.[32] J. Healy, M. Kutay, H. Ozaktas, and J. Sheridan, LinearCanonical Transforms: Theory and Applications , Vol.198 (2016).[33] “xrt (XRayTracer) — xrt 1.3.3 documentation,” (2020),[Online; accessed 23. Apr. 2020].[34] R. Chernikov and K. Klementiev, Advances in Compu-tational Methods for X-Ray Optics IV , 1038806 (2017).[35] J. Bahrdt and U. Flechsig, in
Gratings and GratingMonochromators for Synchrotron Radiation , Vol. 3150(International Society for Optics and Photonics, 1997)pp. 158–170.[36] O. Chubar, P. Elleaume, S. Kuznetsov, and A. A. Sni-girev, in
Optical Design and Analysis Software II , Vol.4769 (International Society for Optics and Photonics,2002) pp. 145–151.[37] J. W. Goodman,