Propagation algorithms for Wigner functions

We propose an algorithm to remove a parabolic wavefront from a convergent or divergent beam in the Wigner function. Using this approach we numerically collimate the beam. This avoids a dense sampling in phase space to describe a convergent wavefront. Thereby we reduce the required computer memory, but maintain computational accuracy and physical effects. Furthermore, we compare two algorithms, shearing and Radon transform, to propagate the Wigner function in free space. We use the fast Fourier transform to accurately perform shearing. However, zero-padding is necessary to circumvent aliasing. We prove that the Radon transform is a more efficient approach for a long propagated distance.


Background
The Wigner function is a helpful tool to analyze optical signals in phase space [1,2].It includes information about ray optics and wave optics [3].In particular for partially coherent light, the Wigner function visualizes the coherence effects in a straightforward manner [4].However, the computation of the Wigner function has certain difficulties.We require a two-dimensional Wigner function to describe light field with one transverse dimension.For a field with two transverse dimensions, the Wigner function spans four dimensions.Consequently, the question of how to efficiently implement the Wigner function becomes an important issue in practice.The main goal of this work is to propagate light by using Wigner functions, while saving computer memory and keeping computational accuracy.
In optical systems, we often encounter strongly convergent or divergent beams.To represent such a beam in phase space requires a dense sampling grid.At first we introduce a method based on numerical collimation, to represent a convergent or divergent beam in phase space with as few sampling points as possible.This preserves all the diffraction effects.The corresponding details are introduced in Removing a parabolic wavefront in phase space.
The Wigner function can be propagated by using the ABCD matrix formalism.The propagation in free space corresponds to a shearing of signals in phase space.However a straightforward implementation of this can severely burden the computer memory.In Shearing and Radon transform we describe another algorithm based on the Radon transform with a better computational efficiency.
In this paper we assume that all the light is within a small numerical aperture (NA).The propagation operators are paraxial.The incident light fields are partially coherent in space with one transverse dimension in this paper.All the algorithms can be extended to fields with two transverse dimensions.

Method 1: Removing a parabolic wavefront in phase space
The phase space of a monochromatic partially-coherent beam is given by the Wigner function, where x and u denote the spatial and angular variables in phase space, Γ(x 1 , x 2 ) representing the correlation function between every two arbitrary transverse positions given by x 1 and x 2 , x = (x 1 + x 2 )/2, Δx = x 1x 2 .As the angle u and the spatial distance Δx are Fourier conjugated by their definition in Eq. 1, the samplings on the spatial and angular axes follow the following relation, where u max denotes the maximum value on the angular axis and λ is the wavelength of the light field, the symbols d(Δx) and d(x) defining the spatial distance between two adjacent sampling points on the Δx and x axes respectively, N x representing the total sampling points on the x axis.Here we use the term (N x -2) because we define one less sampling point on the nonnegative x values than on the negative x values.Based on Eq. 2, the sampling on the angular axis u is automatically defined once the sampling on the spatial axis x is chosen.If a beam is strongly convergent or divergent, a large angular range on the angular axis is required to fully describe the beam of a large angle cone.Thus, according to Eq. 2, we are forced to have a dense sampling on the spatial axis x in order to fulfill the range requirement on the angular axis.Eventually this leads to aliasing (Fig. 1a).However, if the optical component does not need so many sampling points to describe its spatial structure, a dense sampling on the spatial axis x is a waste of computer memory.Therefore we convert the convergent or divergent beam into a quasi-collimated beam by removing a parabolic wavefront.For a quasi-collimated beam, the angular axis does not need a big range any more.Then the sampling density on the spatial axis x can be reduced.
According to Siegman [5], a beam with a convergent wavefront with a radius of curvature R propagating for a distance of L generates the same diffraction effects as this beam without the convergent wavefront propagating for a transformed distance.The equation for the transformed distance is written as follows.
After the propagation we need to scale the collimated beam width by a factor of R/(R + L), so that it is comparable to the original beam width.
Figure 1b depicts a simple system as an example.In the original system, the kinoform lens with a focal length of 16 mm (i.e., R = −16 mm) focuses a quasi-collimated input beam into a convergent beam.The kinoform lens has a diameter of 2 mm and a uniform groove height Δz = 3.5λ, where λ denotes the wavelength.We define the refractive index of the lens as n = 2.0 with λ = 0.6328 μm.The phase difference generated by the groove height is 2π(n-1)Δz/λ = 7π, leading to destructive interference.Thus the outgoing convergent beam carries additional diffraction effects.The beam propagates in free space with a distance L = 10 mm after leaving the kinoform lens.In the transformed system, we place a paraxial negative lens right behind the kinoform lens to collimate the convergent beam.According to Eq. 3, the collimated beam propagates in free space with a transformed distance L trans = 26.7 mm.After this propagation we desire the same diffraction effects as in the original system.
To apply this transformation in the Wigner function, we insert a paraxial negative lens directly into Eq. 1.This lens only removes a convergent curvature from the wavefront, without changing any other phase effects from the kinoform lens.It is expressed in an equation as follows.
where W denotes the Wigner function of light leaving the paraxial negative lens, Γ source being the correlation function of the incident light onto the kinoform lens, t j representing the phase modulation function of a surface (i.e., either the kinoform lens or the paraxial negative lens), Δz j referring to the height of the corresponding surface.For the paraxial negative lens we define Δz par = exp[−iπx 2 /(λf par )], where f par is its focal length.The profile of the kinoform lens is shown in Fig. 2a.
Results and Discussion of Method 1 Thus the convergent beam is collimated.We propagate this collimated beam with the distance L trans = 26.7 mm and scale the transverse beam width after the propagation.Figure 2g shows the nearly identical transverse intensities, with a Pearson correlation coefficient [6] of 0.9996, between the original and the transformed systems.It indicates that the diffraction effects caused by the kinoform lens are preserved in the transformed system.In addition, the angular axis in Fig. 2e has half the range of the angular axis in Fig. 2c.An angular range of −0.04 ≤ u ≤ 0.04 (radian) is sufficient for the collimated beam, whereas a range of −0.08 ≤ u ≤ 0.08 (radian) is required for the convergent beam.Thus the grid in the Wigner function is decreased from 1024 × 1024 pixels to 512 × 512 pixels by using the transformed system.This method saves a factor of 4 of the computer memory.
Figure 3 depicts more intermediate results to explain the phase effects introduced by the paraxial negative lens.Based on Eq. 4 and 5, the product between t kino and t par includes a sum of the surface heights given by the kinoform lens and the paraxial negative lens, seen in Eq. 5.
The remaining height in Eq. 6 (red curve in Fig. 3a) is step jumps with a uniform height given by the kinoform Furthermore, our method offers an alternative to magnify the convergent beam near the focus region.As the collimated beam in the transformed system has a larger transverse width, we obtain a beam profile of more pixels filled with non-zero values.When we are interested in the exact focus of the convergent beam, we need to propagate the collimated beam to infinity (i.e., far field).A far-field propagation in phase space is performed by a 90°rotation, i.e., a Fourier transform.
It is worth noting that the Wigner function can describe non-paraxial light [7].However, the numerical collimation method presented in this section has a limitation.Since we use the thin element approximation to process surfaces and apply ABCD matrices to propagate light in free space, our propagation operators are paraxial.Therefore, the algorithms are only valid for systems within a small NA.If the condition of a small NA is fulfilled, the convergent beam is allowed to contain a wavefront deviating from a parabola, e.g., a non-perfect focusing beam.

Method 2: Shearing and Radon transform
The paraxial propagation of the Wigner function is associated with an ABCD matrix 1 z 0 1

!
, where z defines the propagated distance [8].The signals in phase space before and after the propagation share a shearing relation, i.e., where W and W' denote the Wigner functions before and after the propagation.To perform shearing numerically, each row in W(x, u) is shifted horizontally in phase space by a distance of x o = u o z, where x o marks the shifted distance and u o is the angular coordinate of a horizontal row.
A shift in space corresponds to multiplying a phase factor in the Fourier domain.Thus the spatial shift of each row can be represented by two Fourier transforms, where x and ν are Fourier conjugated.
The advantage of employing a Fourier transform in shearing is that it returns an accurate value at each shifted pixel.However, shearing with a Fourier transform also has a disadvantage.A large propagated distance z in free space leads to a large shifted distance of x o in phase space.Some signals in phase space are sheared outside the original region.The discrete Fourier transform mirrors these signals back into the original region and causes aliasing.To avoid this error, one must make sure that the computational region in phase space is wide enough to support shearing.Otherwise zero-padding should be used to increase the computational region.
As we are interested in light intensity along propagation, another method called Radon transform [1] reaches the same goal as shearing.A Radon transform projects an image along a radial line oriented at a specific angle to derive the intensity distribution.Figure 4a-c give a qualitative illustration.Figure 4a is the phase space at the propagated distance z generated by shearing the phase space at the propagated distance 0. The vertical lines in Fig. 4a indicate the integration direction of signals for the transverse intensity at the propagated distance z. Figure 4b results from an inverse shearing of Fig. 4a in the directions of the red arrows.This corresponds to a back propagation from the distance z to 0. Thus Fig. 4b returns the phase space at the propagated distance 0. The integration lines in Fig. 4a are sheared into a tilted orientation in Fig. 4b.They still imply the integration direction for the transverse intensity at the propagated distance z.Alternatively we may perform a rotation of the phase space of Fig. 4b into Fig.4c, and integrate the signals along the integration lines vertically in Fig. 4c to derive the transverse intensity at the propagated distance z.This is what we call as a Radon transform.
Figure 4d and e depict the schematic diagrams of the geometries for shearing and for the Radon transform.Assume there are two points A = (0, u o ) and B = (0, −u o ) in phase space at the propagated distance 0. For the propagated distance z, the shearing angle is defined as θ = arctan(z).Points A and B in phase space are sheared to coordinates of (u o z, u o ) and (−u o z, −u o ) respectively.The distance between these two points in the transverse intensity at the propagated distance z is denoted by Δx AB = 2u o z.The total width of the transverse intensity is Δx max = 2x m + u m z, where x m and u m define the maximum values on the spatial and angular axes in phase space at the propagated distance 0.
In the Radon transform (Fig. 4e), at first we divide the angular axis by a factor of ϒ = u m /x m , so that the two axes in phase space have the same unit and remain uniform during rotation.The projected intensity after a rotation of angle α yields a maximum width of Δx max ' = 2x m cosα + 2(u m /ϒ)sinα.Inside the projected intensity, the distance between A and B is Δx AB ' = 2(u o /ϒ)sinα.In order to implement the Radon transform, we have to find out the relation between the rotation angle α and the propagated distance z.This is achieved by keeping the ratio between Δx max and Δx AB , and the ratio between Δx max ' and Δx AB ' equal.
From Eq. 9 we derive the rotation angle α = arctan(ϒz).The ratio factor for the transverse intensity width between shearing and the Radon transform is ϒz/sinα.

Results and Discussion of Method 2
Figures 5 and 6 show an example of using the two algorithms to propagate the Wigner function.We use a quasicollimated Schell beam [9] with a nearly flat-top transverse intensity profile as the incident light source.The beam is propagated through an asymmetric lens with half aspheric and half kinoform profiles (Fig. 5b) to generate an asymmetric phase space (Fig. 5c), so that the phase space orientation after a rotation is more distinguishable.Figure 6a-c compare the two algorithms for a short propagated distance (z = 2 mm).In this case both algorithms produce nearly identical results.However, with a larger propagated distance (Fig. 6d-f, with z = 8 mm), the shearing stretches the signals in phase space outside and causes aliasing (highlighted by the black parallelogram in Fig. 6d).This kind of error can be avoided by using the rotation algorithm.A phase space with a width of 2 ffiffiffiffiffiffiffiffi 2x m p is large enough to support a rotation of any angle.For a propagated distance from zero to infinity, the rotation angle α corresponds to 0°up to 90°.According to Larkin [10] and Lohmann [11], an arbitrary rotation matrix can be expressed as a product of three shearing matrices.When the fast Fourier transform is used to implement these shearing operations, this way produces accurate results for the rotated image.
Figure 7 sketches the intensity along the propagated beam based on the results in Fig. 6.With the aliasing error in Fig. 6d, the transverse beam width in Fig. 7a is always confined to limited extent when z > 6 mm.This unphysical result is avoiding in Fig. 7b by using the Radon transform, as predicted by Fig. 6.
Both in shearing and the Radon transform, there are regions filled with zeros (shown as white arears in Fig. 4a  and c).These regions do not influence the results of the calculation but require computer memory.The larger the zero regions are, the less efficient the calculation becomes.Figure 8 compares the area of zero regions between the two algorithms based on the phase space given by Fig. 5c.For a short propagated distance where zero-padding is unnecessary for shearing, shearing yields less zero area than the Radon transform.Thus in our example in Figs. 5, 6 and 7, for a propagated distance of z < 3 mm, shearing has advantage of requiring less memory than the Radon transform.For a larger propagated distance, the area of zeros in the Radon transform does not increase after exceeding a rotation angle α of 45°.However, the area of zeros in shearing grows linearly with the propagated distance.It thus demands a linear rise of the computer memory, making this method impractical.

Conclusions
Our methods contribute to a fast implementation to paraxially propagate beams in phase space, particularly for partially coherent light.
We show how to apply a method known for working with coherent light to problems of partial coherence.A parabolic wavefront is subtracted from a convergent beam.The beam is converted into quasi-collimated form.The required sampling density to represent this beam in phase space is thus reduced.The diffraction effects in propagation are preserved after the parabolic wavefront is removed.Furthermore, this approach offers an alternative to magnify the convergent beam near the focus region by observing the collimated beam at the physically equivalent distance.
Besides, we compare the efficiency of two algorithms, shearing and Radon transform, for propagating the Wigner function in free space.For a large propagated distance, the shearing method requires zero-padding to avoid aliasing.Its demand on the computer memory grows linearly with the propagated distance.In contrast to this, the Radon transform keeps the computer memory within a finite limit.
The advantages of our proposed methods are even greater with the operation of a four dimensional phase space, even though this is not specifically shown in this paper.

Fig. 1
Fig. 1 Removing a parabolic wavefront from a convergent beam to overcome the under-sampling problem in phase space.a Under-sampling problem in phase space, b Collimating a convergent beam produced by a kinoform lens

Figure 2
Figure 2 compares the phase space of the original and the transformed systems given by Fig. 1b.In the original system, the focusing effect produced by the kinoform lens is expressed by a general tilt of all the signals in Fig. 2c.The additional diffraction effects are indicated by the oscillatory ripples in Fig. 2c.In the transformed system the paraxial negative lens introduces an extra defocusing effect to the convergent beam.It brings the signals back to the horizontal orientation in Fig. 2e.

Fig. 2
Fig. 2 Phase space and propagated beam paths in the original and transformed systems.The lens profile in a does not have perfectly uniform groove heights due to finite samplings.a Profile of the kinoform lens, b Phase space of the incident light, c Phase space of the convergent beam leaving the kinoform lens, d Intensity of the propagated beam path in the original system, e Phase space of the collimated beam leaving the paraxial negative lens, f Intensity of the propagated beam path in the transformed system, g transverse intensity after propagation

Fig. 3
Fig. 3 Correlation functions of light inside the transformed system.The incident light is a Gaussian Schell-model beam [12] with a beam width of 400 μm and a coherence length of 100 μm at the wavelength of 0.6328 μm. a Surface height, b Real values of the correlation function of the incident light, imaginary parts being zero, c Phase of the correlation function of light leaving the Kinoform lens, d Phase of the correlation function of light leaving the paraxial negative lens

Fig. 4
Fig. 4 Geometry of shearing and the rotation.a Phase space at the distance z, b Reverse shearing of (a), phase space at the distance 0, c Rotation of (b), d Details of the shearing geometry, e Details of the rotation geometry

Fig. 5
Fig. 5 Phase space results of a quasi-collimated Schell beam passing an asymmetric lens.a Phase space of incident light, b Profile of the lens, c Phase space of outgoing lights.(a) and (c) share the same colorbar

Fig. 7 Fig. 8
Fig. 7 Intensity along the propagated beam obtained by shearing (a) and the Radon transform (b)