## Abstract

This tutorial gives an overview of the use of the Wigner function as a tool for
modeling optical field propagation. Particular emphasis is placed on the spatial
propagation of stationary fields, as well as on the propagation of pulses through
dispersive media. In the first case, the Wigner function gives a representation of
the field that is similar to a radiance or weight distribution for all the rays in
the system, since its arguments are both position and direction. In cases in which
the field is paraxial and where the system is described by a simple linear relation
in the ray regime, the Wigner function is constant under propagation along rays. An
equivalent property holds for optical pulse propagation in dispersive media under
analogous assumptions. Several properties and applications of the Wigner function in
these contexts are discussed, as is its connection with other common phase-space
distributions like the ambiguity function, the spectrogram, and the Husimi,
*P*, *Q*, and Kirkwood–Rihaczek functions. Also
discussed are modifications to the definition of the Wigner function that allow
extending the property of conservation along paths to a wider range of problems,
including nonparaxial field propagation and pulse propagation within general
transparent dispersive media.

©2011 Optical Society of America

## 1. Introduction

One of the most used mathematical tools in physics and engineering is the Fourier
transformation, which converts a function of a given variable into a function of
another, *Fourier conjugate* variable. Particularly in situations where
the mathematical description of the physical process is independent of the choice of
origin for the initial variable, the Fourier conjugate variable has clear physical
significance. Consider, for example, a function of time representing the local air
pressure for a sound or the electric field for optical radiation. Fourier transformation
allows us to express this function in terms not of time but of frequency, which our
senses perceive as pitch or color. When considering the spatial distribution of a
monochromatic beam, on the other hand, Fourier transformation replaces the spatial
variable with one associated with direction of propagation. In quantum mechanics,
Fourier transformation over the position coordinates of a particle’s wave function leads
to a distribution in momentum, while Fourier transformation over time gives an energy
distribution.

In any of these specific contexts, the well-known mathematical properties of the Fourier transformation acquire a physical interpretation. In particular, a function and its Fourier transform cannot be simultaneously arbitrarily narrow, since the product of the widths of their squared moduli (defined as the standard deviations of their variables) must be equal to or greater than a given constant. This purely mathematical property, sometimes referred to as the space–bandwidth product theorem, has profound physical implications in the different cases mentioned earlier. Within the context of sound, it tells us that a musical note cannot be arbitrarily short and yet have a well-defined pitch. When applied to monochromatic wave beams over a transverse coordinate, it implies that a beam cannot be simultaneously arbitrarily thin and well collimated. In quantum mechanics (according to the Copenhagen interpretation), this relation gives rise to the famous Heisenberg uncertainty principle that limits the simultaneous knowledge about the position and momentum of a particle. Due to this latter interpretation, many authors refer to the general space–bandwidth relation of Fourier transform pairs as the uncertainty relation.

In physical applications, the Fourier transformation provides a change in the
representation of a system, from a description in terms of one physically meaningful
variable (e.g., time or position) into one in terms of another (e.g., frequency,
direction, or momentum). Nevertheless, in many of these contexts, our intuition is more
inclined towards a representation in terms of *both* Fourier conjugate
variables. For example, we think of a musical piece as a sound distribution in terms of
both time and frequency, represented as a musical score. Similarly, when designing and
modeling macroscopic imaging or illumination optical systems, one thinks of light as a
distribution over a bundle of rays, each of these specified by both its position and
direction. The same is true in mechanics, where the classical model, closer than the
quantum one to our everyday perception of the world, requires the specification of both
the position and momentum of a particle. Notice that these three simplistic models
(musical notation, ray optics, and classical mechanics) disregard the limitations
imposed by the space–bandwidth product theorem or uncertainty relation.

There is a family of mathematical transformations that allow the description of a
function in terms of both a variable and its Fourier conjugate. When applied to the
physical situations mentioned earlier, these transformations mimic the so-called
*phase-space distributions* used in the simplified physical models
(e.g., a musical score, a radiance distribution giving the weight of all rays in an
optical system, or a distribution of classical particles in terms of both their position
and velocities). In this context, the term *phase space* refers to the
extended space composed of both sets of variables (time–frequency, position–direction,
position–momentum). For this reason, these mathematical transformations are sometimes
referred to as *phase-space quasi-probability distributions* or
*mock-phase-space distributions* [1], although most authors drop the cautionary qualifiers and refer to them
simply as *phase-space distributions*, as we will do here.

One of the most widely used distributions of this type is the Wigner function, proposed
by Wigner [2] to represent a quantum state as a
distribution in terms of both position and momentum, reminiscent of a
statistical–mechanical distribution of classical particles. The same mathematical
transformation was independently proposed in the signal analysis community by
Ville [3] for studying time signals, and by
both Dolin [4] and Walther [5] to represent optical wave fields as ray distributions analogous to
what is known in the theory of radiometry as the *radiance* or
*specific intensity* [6] and in
image science as the *light field* [7]. Excellent reviews on the Wigner function and other phase-space
distributions are given in the books by Cohen [8]
and Mecklenbräuker and Hlawatsch [9], and the
review articles by Lee [10], Agarwal and
Wolf [11–13], Balazs and Jennings [1], and
Hillery *et al.* [14], to name a
few. For applications in optics, the reader can also consult, for example, the book
edited by Testorf *et al.* [15],
the book by Torre [16], and the review articles
by Dragoman [17] and Sheridan
*et al.* [18].

This tutorial gives a review of the Wigner function and its use within several contexts in classical optics. While many other phase-space representations have also been proposed (as discussed in Sections 3 and 6), the Wigner function presents two properties that make it unique in many applications: (i) its so-called marginal projections yield observable physical quantities, and (ii) under suitable situations, it satisfies the same conservation rules under spatial propagation or temporal evolution as the simpler physical model it mimics. For example, in the case of paraxial fields propagating through simple optical systems, the Wigner function plays the role of a distribution of weights for the rays, because its value for each ray remains constant under propagation across the system, and because its integral over all rays passing through a point leads to the local optical intensity. The main emphasis of this work is precisely on the property of conservation along paths under propagation or evolution, its usefulness in practical applications, and range of validity. Also presented here are generalizations to the standard definition of the Wigner function that allow extending this property to more general situations, including nonparaxial propagation and pulse propagation in arbitrary transparent dispersive media.

The outline of this article is as follows. Some unitary transformations (including the Fourier transformation) that are useful in the description of certain simple optical systems are described in Section 2. A very simple phase-space representation, the spectrogram, is introduced in Section 3. The definition and some properties of the Wigner function are presented in Section 4 for the case of deterministic functions, and the more general definition of this representation for the case of correlations is given in Section 5. An overview of several phase-space distributions is given in Section 6. The Wigner function’s marginal and conservation properties are discussed in Section 7 for the case of classical paraxial optical fields propagating through simple systems and in Section 9 for the case of optical pulses propagating in media presenting linear dispersion. In both these contexts (as well as in quantum optics), a method can be used to recover the coherence properties of the field based on the measurement of the marginal projections of the Wigner function at different stages. This scheme, known as phase-space tomography, is discussed briefly in Section 10. The conservation property of the Wigner function under classical propagation, however, is valid only under certain limiting cases, as discussed in Section 8. In Section 11 we discuss a procedure for defining generalizations to the Wigner function that satisfy the desired propagation and marginal properties for other systems, including nonparaxial propagation through homogeneous media and pulse propagation through more general dispersive media. Finally, some brief concluding remarks are given in Section 12.

## 2. The Fourier Transform and Other Unitary Linear Transformations

Before introducing the Wigner function and other phase-space distributions, let us first review some simple integral transformations that appear naturally in the approximate description of certain optical systems. All these transformations depend linearly on the function they are applied to.

#### 2.1. The Fourier Transform

Let us start by introducing the generic definition of Fourier transformation that is
used in this article. Consider a scalar, square-integrable function $f\left(\mathbf{q}\right)$, where $\mathbf{q}=({q}_{1},{q}_{2},\dots ,{q}_{N})$ is an *N*-dimensional variable. The Fourier
transform $\stackrel{\u0303}{f}=\stackrel{\u02c6}{\mathcal{F}}f$ of this function is defined here as

*K*is a real constant [with units of ${(\mathbf{q}\cdot \mathbf{p})}^{-1}$] dictated by the physical context. Throughout, integrals run from minus to plus infinity unless otherwise specified. The inverse Fourier transformation has the form

- (a) time
*t*and frequency*ω*, for which $N=1$, $K=-1$, $q=t$ and $p=\omega $; - (b) the transverse position and transverse direction of a beam at a fixed frequency
*ω*, where $N=2$, $K=k=\omega /c$ (the free-space wavenumber, given by the ratio of the frequency and the speed of light*c*), $\mathbf{q}=\mathbf{x}=(x,y)$ (the components of the position vector perpendicular to the optical axis), and**p**is the transverse part of a vector indicating a direction of propagation; - (c) two quantized field quadratures, i.e., the “real” and “imaginary” parts of the complex representation of a Cartesian component of the quantized electric field, with $N=K=1$ and
*q*and*p*being the two quadratures.

**q**and

**p**are the position and momentum variables, respectively, $K={\u0127}^{-1}$ (the inverse of the reduced Planck constant), and

*N*is the spatial dimensionality of the problem.

There are several well-known properties of the Fourier transformation (see, e.g., Ref. [19] for a comprehensive list). Let us mention here only a couple:

(i) *Parseval–Plancherel theorem*. This
theorem states that the inner product of two functions *f* and
*g* is the same in either representation, due to the unitarity of
the Fourier transformation, i.e.,

*Space–bandwidth product theorem or uncertainty relation*. The standard deviations of $\left|f\right(\mathbf{q}){|}^{2}$ and $\left|\stackrel{\u0303}{f}\right(\mathbf{p}){|}^{2}$ for each Cartesian component of their variables cannot simultaneously be arbitrarily small, since their product must satisfy for $i=1,2,\dots ,N$, where the standard deviations are defined according to

#### 2.2. The Fractional Fourier Transform

A generalization of the Fourier transformation that gives access to a continuum of
functions joining the original function and its Fourier transform was proposed by
Condon [21–24]. This transformation, referred to as the fractional Fourier
transformation of degree *θ*, is defined here as

*σ*(assumed to be positive) with units of $\sqrt{q/p}$ was introduced. In applications where

**q**and

**s**have the same units, we can set $\sigma =1$. The fractional Fourier transformation is unitary, so an analog of the Parseval–Plancherel theorem also exists for it:

Notice that, for $\theta =0$ and $\theta =\pi /2$, the fractional Fourier transform reduces, respectively, to scaled versions of the original function and its Fourier transform, i.e.,

(The verification of the first relation requires the use of the method of stationary phase [25] under the limit $\theta \to 0$.) The fractional Fourier transform presents an additivity property, which takes a particularly simple form in cases where*σ*can be set to unity: two consecutive transformations with degrees ${\theta}_{1}$ and ${\theta}_{2}$, respectively, are equivalent to a single transformation with degree ${\theta}_{1}+{\theta}_{2}$. As will be discussed later in this article, the spatial field propagation through a quadratic gradient index waveguide or the temporal evolution of the quadratures of a quantized field are described mathematically by fractional Fourier transformation. In these contexts, the additivity property has a straightforward physical interpretation, since

*θ*is proportional to the propagation distance or time. Due to both Eqs. (9a) and (9b) and to this additivity property, many authors use, instead of the degree

*θ*, the order $\gamma =2\theta /\pi $, so that $\gamma =1$ corresponds to the standard Fourier transformation [24]. This way,

*γ*can be thought of as the power to which the Fourier transformation operator $\stackrel{\u02c6}{\mathcal{F}}$ is raised before acting on a function.

It must be noted that a more general version (sometimes called anamorphic) of the
fractional Fourier transformation can be written in the multidimensional case, where
a different *θ* is used for each Cartesian component of the
transformation variables [26]. This more
general definition is in turn a special case of an even more general family of
transformations, described in what follows.

#### 2.3. Linear Canonical Transformations

A more general class of unitary transformations that includes the regular and
fractional Fourier transformtions is now described. These transformations, known as
*linear canonical transformations*, are used in paraxial optics and
nonrelativistic quantum mechanics to describe simple systems [27–30]. Consider a linear
transformation whose kernel has constant amplitude, like those of the standard and
fractional Fourier transforms:

*a*is a constant that ensures that the transformation is unitary, and $v({\mathbf{q}}^{\prime},\mathbf{q})$ is a real function. Let us assume that this function is well approximated by its second-order Taylor expansion in both variables, i.e.,

**v**, and the matrices ${\mathbb{V}}^{\u2033}$, ${\mathbb{V}}^{\prime}$, and $\mathbb{V}$, are constant expansion coefficients. Clearly, ${\mathbb{V}}^{\u2033}$ and $\mathbb{V}$ must be symmetric matrices. Since the only effect of the constant component ${v}_{0}$ is an unimportant global phase factor of the transform in Eq. (10), we will set this constant to zero. Further, it can be shown that the analog of Eq. (3) is given by

Let us now assume that *v* satisfies the symmetry property $v(-{\mathbf{q}}^{\prime},-\mathbf{q})=v({\mathbf{q}}^{\prime},\mathbf{q})$, so that ${\mathbf{v}}^{\prime}=\mathbf{v}=\mathbf{0}$ (where **0** denotes a zero vector). It turns out to be
convenient to rewrite Eq. (11) as

Note that in Eq. (17), the transformation $\stackrel{\u02c6}{\mathcal{C}}$ depends on $\mathbb{S}$, which is a $2N\times 2N$ matrix in which $\mathbb{A}$, $\mathbb{B}$, and $\mathbb{D}$ are embedded:

*symplectic*, i.e., that it satisfies the property

At a first glance, characterizing the linear canonical transformation in terms of the $2N\times 2N$ matrix $\mathbb{S}$ instead of the three $N\times N$ matrices ${\mathbb{V}}^{\u2033}$, ${\mathbb{V}}^{\prime}$, and $\mathbb{V}$ might seem unnecessarily convoluted. However, this parametrization is convenient, since it confers on the linear canonical transformation a group property similar to (and more general than) that of the fractional Fourier transformation: two consecutive linear canonical transformations are equivalent to a single one whose parameter is the product of the parameters of the individual transformations, i.e.,

Let us now discuss four special cases that appear often in the study of optical systems, starting with the fractional Fourier transform.

### 2.3a. Fractional Fourier Transformation

The particular case of the fractional Fourier transformation corresponds to $\left\{{\stackrel{\u02c6}{\mathcal{F}}}_{\theta}f\right\}\left(\mathbf{s}\right)=exp(\mathrm{i}\theta /2)\left\{\stackrel{\u02c6}{\mathcal{C}}\right[{\mathbb{S}}_{\mathrm{FrFT}}(\theta ,\sigma )\left]f\right\}\left(\sigma \mathbf{s}\right)$, where the symplectic matrix is given by

### 2.3b. Scalings

Norm-preserving scalings can be defined according to

*b*is a real nonzero constant. This is a linear canonical transformation corresponding to the matrix

It is straightforward to see that ${\mathbb{S}}_{\mathrm{Sc}}\left({b}_{2}\right){\mathbb{S}}_{\mathrm{Sc}}\left({b}_{1}\right)={\mathbb{S}}_{\mathrm{Sc}}\left({b}_{1}{b}_{2}\right)$. More general, anamorphic versions of this transformation can be written where a different scaling applies in each direction. If these directions are aligned with the coordinate axes, this anamorphic scaling corresponds to a linear canonical transformation where $\mathbb{S}$ is diagonal.

### 2.3c. Fresnel Transformation

The Fresnel transformation, used to model paraxial free propagation in optics or nonrelativistic free evolution in quantum mechanics, is given by

### 2.3d. Chirping

Chirping corresponds to a simple multiplication by a quadratic phase factor:

#### 2.4. Phase-Space Shifts

A different type of linear unitary transformation corresponds to combinations of shifts in both sets of Fourier conjugate variables. Let us define an operator that causes a shift in the origin of a function as

*phase-space shift*. Note, however, that these two types of shift do not commute:

The action of many simple systems is expressible as a combination of linear canonical
transformations and phase-space shifts. This includes the general linear
transformation in Eqs. (10) and
(11) in which the linear
coefficients ${\mathbf{v}}^{\prime}$ and **v** do not vanish.

## 3. Simple Phase-Space Distributions: the Windowed Fourier Transform and the Spectrogram

In this and the next section, we discuss different mathematical definitions of
phase-space representations, i.e., representations of a function in terms of both sets
of Fourier conjugate variables. Perhaps the most intuitive representation of this type
is the *windowed Fourier transform*. Here, the part of a function
*f* around a prescribed location **q** is isolated through
multiplication by a window function that differs significantly from zero only within a
small neighborhood around **q**. This product is then Fourier-transformed, in
order to give the “local frequency content” of the function at that region. The windowed
Fourier transform is then a function of **q**, the central position of the
window, and **p**, the Fourier variable. Mathematically, the windowed Fourier
transform can be written as

*f*and

*w*has a similar form:

*f*with a phase-space-shifted version of the window function by using the phase-space-shift operation defined in Eq. (35):

The windowed Fourier transform is linear in *f* and complex. For many
practical applications, it is more useful to use the squared modulus of this
transformation, referred to in the signal analysis community as the
*spectrogram*:

In order to visualize the definition of the windowed Fourier transform and the spectrogram, consider a one-dimensional ($N=1$) function given by the sum of two identical Gaussians, one shifted by the action of $\stackrel{\u02c6}{\mathcal{T}}({q}_{0},{p}_{0})$ and multiplied by a constant phase factor $exp\left(\mathrm{i}\Phi \right)$:

Consider the case ${q}_{0}=8$, ${p}_{0}=4$ and $\Phi =0$. Also, let us use $K=1$ for simplicity, and employ a Gaussian window, $w\left(q\right)={\pi}^{-1/4}exp(-{q}^{2}/2)$. The first row of the movie frame (Media 1) in Fig. 1 shows the real (blue) and imaginary (red) parts of $f\left({q}^{\prime}\right)$ for ${q}^{\prime}\in [-16,16]$. The second row shows the window function $w({q}^{\prime}-q)$ over the same range, for varying*q*. The third row shows the products of each of the curves in the first row with that in the second row, i.e., the real and imaginary parts of a local detail of

*f*isolated by the window function. The spectrogram for the corresponding value of

*q*, given by the squared modulus of the Fourier transform of the function in the third row, is shown in the fourth row. It is clear that the spectrogram has significant contributions only when the window overlaps the regions where

*f*is significant.

The movie frames (Media 2 and
Media 3) in Figs. 2 and 3 present the spectrogram of this
function over phase space for several values of ${q}_{0}$, ${p}_{0}$ and *Φ*, for a rectangular window $w\left(q\right)=\mathrm{rect}\left(2q\right)/\sqrt{2}$ (which equals $1/\sqrt{2}$ for $\left|q\right|\le 1$ and vanishes otherwise) and for a Gaussian window $w\left(q\right)={\pi}^{-1/4}exp(-{q}^{2}/2)$, respectively. The blue dot indicates the coordinates $({q}_{0},{p}_{0})$. Unless both ${q}_{0}$ and ${p}_{0}$ are of the order of unity or smaller, the spectrogram is composed of
two separate and similar distributions, one centered at the origin (associated with ${f}_{1}$) and one at $({q}_{0},{p}_{0})$ (associated with ${f}_{2}$). This example shows that the action of $\stackrel{\u02c6}{\mathcal{T}}$ is indeed to shift a distribution in phase space by the amounts
specified by its two arguments. The movies also show that the effect of the phase
*Φ* between the two separate contributions on the spectrogram is
negligible except when the two contributions are so close that they overlap
significantly.

One drawback of the windowed-Fourier transform and the spectrogram as representations of local frequency content is that they do not have a unique definition, since they depend on the choice of a window function. In some cases, the window function is dictated by the physics of the problem, as we will discuss later. In general, though, the choice of window function limits the range of frequencies that can be described. For this reason, it can be useful to have a phase-space distribution that does not rely on the choice of an ancillary function. One such distribution is the Wigner function, described in the next section.

## 4. The Wigner Function

We now introduce the Wigner function as a way to represent the function
*f* jointly in terms of both **q** and **p**,
without the need for an ancillary function:

*f*in a very similar form:

The definition of the Wigner function in Eq. (41) is illustrated for the case $N=1$ in Media 4 (Fig. 4), for the same function as in Fig. 1, namely, the function in Eqs. (40) with $({q}_{0},{p}_{0})=(8,4)$ and $\Phi =0$. The first and second rows show plots, respectively, of $f(q+{q}^{\prime}/2)$ and ${f}^{\ast}(q-{q}^{\prime}/2)$ as functions of ${q}^{\prime}\in [-16,16]$. The value of *q* is varied from −12 to 12 in
Media 4. In these plots, the real and imaginary
parts are shown as blue and red curves, respectively. Note that, since they are plotted
as functions of ${q}^{\prime}$, these functions are stretched by a factor of 2. The third row shows
the product of these two functions, and the fourth row shows the Fourier transform of
this product, corresponding to the Wigner function as a function of *p*
for the specified value of *q*. A comparison of the graphic
representation in Fig. 4 with that for the
spectrogram in Fig. 1 suggests that, in the case of
the Wigner function, ${f}^{\ast}(q-{q}^{\prime}/2)$ is somewhat analogous to a window function acting on $f(q+{q}^{\prime}/2)$. Note that the product $f(q+{q}^{\prime}/2){f}^{\ast}(q-{q}^{\prime}/2)$ has the same local frequency (i.e., the same rate of change of the
phase in terms of ${q}^{\prime}$) for small ${q}^{\prime}$ as $f(q+{q}^{\prime})$. Note also that, for this choice of *f*, the Wigner
function is significant in three regions: when the contributions due to ${f}_{1}$ in the first two rows overlap (i.e., for $q\approx 0$), when the contributions due to ${f}_{2}$ in the first two rows overlap (i.e., for $q\approx {q}_{0}$), and when the contributions due to ${f}_{1}$ in one row overlap with those due to ${f}_{2}$ in the other row (for $q\approx {q}_{0}/2$). This third contribution (visible in the static frame of Fig. 4) will be discussed in more detail later.

Comprehensive summaries of the properties of the Wigner function are given in Refs. [1,8–14]. Here we concentrate on a few of them that are of particular interest for the applications mentioned in what follows.

#### 4.1. Reality and Possible Negativity

The Wigner function is always real, because it is the Fourier transform over ${\mathbf{q}}^{\prime}$ of the Hermitian function $f(\mathbf{q}+{\mathbf{q}}^{\prime}/2){f}^{\ast}(\mathbf{q}-{\mathbf{q}}^{\prime}/2)$. That is, changing the sign of the variable of integration ${\mathbf{q}}^{\prime}$ in the integrand in Eq. (41) is equivalent to complex conjugation, so the Wigner function equals its complex conjugate. It must be noted, however, that the Wigner function is not necessarily non-negative; for example, the Wigner function plotted in the fourth row of Fig. 4 presents positive and negative oscillations when $q\approx {q}_{0}/2=4$. The possible negativity of the Wigner function for some values of its arguments is sometimes seen as a conceptual problem when trying to ascribe to it a literal physical meaning as a power or probability density. However, as will be discussed later, these negative values are needed for the description of coherent effects. Note also that multiplying $f\left(\mathbf{q}\right)$ by a global constant phase factor does not change the resulting Wigner function, given its bilinear definition.

#### 4.2. Marginals

Among the main properties of the Wigner function are the so-called marginal relations: the integral over one set of variables of the Wigner function gives the square modulus of the function in the representation associated with the remaining variable:

**q**and

**p**have different units, we introduce a positive scaling constant

*σ*with units of $\sqrt{q/p}$. Then, the marginal projection over the rotated direction in each $({\sigma}^{-1}{q}_{i},\sigma {p}_{i})$ subspace at an angle

*θ*from the $\sigma {p}_{i}$ axis gives the squared modulus of the fractional Fourier transform defined in Eq. (7):

From any of the marginal relations given earlier, it is easy to see that the integral
of ${W}_{f}$ over all phase space gives the squared norm of *f*:

#### 4.3. Linear Mapping under Linear Canonical Transformations and Phase-Space Shifts

The Wigner representation for the result of a linear canonical transformation on a given function is just a simple mapped version of the Wigner representation of the original function, namely,

### 4.3a. Fractional Fourier Transformation

The substitution of Eq. (24) into Eqs. (47) gives

*θ*of the Wigner function in phase space. This is illustrated Media 5 (Fig. 5) for the Wigner function in Eq. (49) with $\sigma =1$.

### 4.3b. Scalings

Similarly, a scaling transformation causes a phase-space volume-preserving scaling of the Wigner function:

**q**and

**p**is divided and multiplied, respectively, by the corresponding scaling factor.

### 4.3c. Fresnel Transformation

Fresnel transformation corresponds to a horizontal shearing of the Wigner function:

### 4.3d. Chirping

Similarly, chirping corresponds to a vertical phase-space shearing of the Wigner function:

### 4.3e. Phase-Space Shifts

A similar simple phase-space rearrangement occurs in the case of phase-space shifts: the Wigner function of $\left\{\stackrel{\u02c6}{\mathcal{T}}\right({\mathbf{q}}_{0},{\mathbf{p}}_{0}\left)f\right\}\left(\mathbf{q}\right)$ is given by

#### 4.4. Inner Product

Suppose that we have two arbitrary square-integrable functions $f\left(\mathbf{q}\right)$ and $g\left(\mathbf{q}\right)$. The integral over all phase space of the product of their corresponding Wigner functions is always a non-negative real number, given by

*f*is the wave function for the state of a system and

*g*is that for a measuring instrument, then the expression on the right-hand side of this equation is related to the result of a measurement of the wave function with the corresponding instrument [40,41]. A similar interpretation exists for the classical measurement of a monochromatic wave field [42–45]. In this context, the left-hand side of Eq. (55) has an intuitive interpretation, as will be discussed in Subsection 7.5.

Notice that for $f=g$, we get

#### 4.5. Superposition

Given the bilinear character of this transformation, the Wigner function of the sum of two functions, $f={f}_{1}+{f}_{2}$, is not, in general, the sum of the corresponding two Wigner functions. Instead we get

*Φ*) between the two contributions. The origin of this interference contribution can be appreciated from Media 4 (Fig. 4): the main contributions to the Wigner function occur not only for values of

*q*at which a significant part of

*f*overlaps with itself within the two top rows of that figure, but also when one significant part of

*f*in the first row overlaps with a different significant part in the second row. At these values of

*q*, the product shown in the third row of the figure is composed of two identical (up to an inversion and complex conjugation) contributions, roughly separated by the same distance as that between the two contributions of

*f*. The Fourier transform of this product (i.e., the Wigner function) is then oscillatory, with frequency proportional to this distance.

The dependence of the Wigner function on the location of the second contribution and
the phase difference between the two contributions is illustrated
in Media 10 (Fig. 10) for the same variation in ${q}_{0}$, ${p}_{0}$, and *Φ* as in Figs. 2 and 3. It is seen there that the
oscillation fringes are parallel to the straight line joining the two contributions.
The fringe structure indicates the phase between the two contributions, regardless of
their separation: the central fringe is positive if the two contributions are in
phase, and negative if they are completely out of phase. The role of these fringes is
clear from considering the marginal properties in Eq. (43a) and (43b),
shown in Fig. 11(a) for
the case $f\left(q\right)={\pi}^{1/4}exp[-{(q+4)}^{2}/2]+{\pi}^{1/4}exp[-{(q-4)}^{2}/2]$: the interference term cancels almost entirely following
integration in *p*, leading to a marginal distribution with two maxima
at $\pm {q}_{0}$; on the other hand, the interference term is responsible for the
oscillations of the marginal in *q*. That is, the negative values of
the Wigner function are needed to account for coherent effects like destructive
interference. Figure 11(b) shows the Wigner function of a superposition of three Gaussians, one
of them multiplied by a linear phase factor in order to displace it from the
*q* axis. Notice that there is an interference term between any two
contributions, and that the frequency of oscillations increases with the distance
between them.

It was discussed in Subsection 4.4 that the effect of a linear canonical transformation and/or phase-space shift is a simple linear mapping of the Wigner function over phase space. These mappings are the only ones that not only preserve phase-space area, but also map any straight line in phase space onto another straight line. In other words, these mappings are the only ones that guarantee that the half-point interference contributions of the Wigner function remain as half-points following the transformation. Therefore, transformations that are not combinations of linear canonical transformations and phase-space shifts cannot correspond to mappings of the Wigner function in general. Notice that following a linear canonical transformation and/or phase-space shift, the fringes of an interference contribution remain aligned with the straight line joining the two contributions giving rise to it, and an increase/decrease in the distance between two contributions brings with it the corresponding decrease/increase in the transverse spacing of the fringes, given the constraint of phase-space area conservation.

#### 4.6. Oscillatory Character of the Wigner Function

The interference terms described in Subsection 4.5 can cause the Wigner function to be highly oscillatory within some
regions of phase space. From a computational point of view, this sometimes
constitutes a major shortcoming of this representation: the Wigner function
represents a function of *N* variables in terms of a potentially
complicated, highly oscillatory function of $2N$ variables. Therefore, except for a limited number of cases that
admit of an analytic description, the numerical modeling of a system by using the
Wigner function requires both doubling the number of variables over which to sample
and potentially increasing significantly the sampling rate in each direction in order
to be able to resolve the fast oscillations.

The problem with oscillations is particularly critical for functions whose Wigner
functions extend over phase space regions much larger than the
minimum-uncertainty-limit cell size, $|2K{|}^{-N}$. One such case of particular relevance in optics and quantum
mechanics is that of functions of the form $f\left(\mathbf{q}\right)=A\left(\mathbf{q}\right)\hspace{0.17em}exp\left[\mathrm{i}K\Omega \right(\mathbf{q}\left)\right]$, where both *A* and *Ω* are real, and
*A* varies slowly compared to the oscillations of $exp\left[\mathrm{i}K\Omega \right(\mathbf{q}\left)\right]$. In this case, one would expect the Wigner function to be localized
in phase space around the *N*-dimensional manifold described by the
equation

As described in detail by Berry and collaborators [46–48], ${W}_{f}$ does indeed present a positive ridge structure that mimics the
manifold defined by Eq. (60), but due
to the interference contributions discussed in the previous section, it also displays
an intricate pattern of fast oscillations involving large values, both positive and
negative, at many other regions of phase space. These regions are composed of all
phase space locations that are midway between any two points along the manifold. Let
us consider the $N=1$ case for simplicity, and use the construction proposed in
Ref. [49] to visualize the location of
these contributions. Consider a plot of the curve $p={\Omega}^{\prime}\left(q\right)$, shown as a black curve in Fig. 12(a). Now superpose on this plot a series of
one-half-scale replicas of the same curve, where, for each replica, the anchor point
for the scaling is a point along the curve, as illustrated in this figure and the
accompanying movie. Let the density of these replicas be proportional to $\left|A\right(q){|}^{2}$, where *q* is the location of the anchor point. The
region occupied by these replicas indicates the region where the Wigner function is
significant, as can be seen from the comparison of this picture with the Wigner
function in Fig. 12(b).
As explained in Ref. [49], the amplitude,
phase, and orientation of these oscillations can be inferred from this geometrical
construction: the amplitude is proportional to the density of replica intersections.
The phase of the oscillations, on the other hand, is proportional to
*K* times the area enclosed by the phase space curve and a straight
line joining the anchor points of the two replicas crossing at the corresponding
point, and the local direction of the fringes is approximately parallel to this line,
as shown in Fig. 12. The complexity of the
pattern depends on the geometry of the curve; for example, more complicated patterns
result for curves including inflection points (like the one shown in Fig. 12), for which more than one pair of points along
the curve can have their half-point at the same phase-space location.

The same type of construction can be used for functions associated with phase-space
manifolds that are not single-valued in **q**. Consider the simplest such
case, given by a one-dimensional Hermite–Gaussian function of the form

*n*corresponds to the number of zeros in the radial direction. As shown in Fig. 13, this result is consistent with the replica-based construction described above. In particular, the amplitude of the oscillations is related to the density of replica intersections, with a maximum at the origin where all replicas cross.

## 5. Partial Coherence

#### 5.1. Correlations and Coherent Mode Expansions

An advantage of bilinear representations like the Wigner function and the spectrogram over linear ones like the windowed Fourier transform is that they can be used to represent not only deterministic functions $f\left(\mathbf{q}\right)$, but also the bilinear autocorrelation of stochastic functions at two values of their argument, i.e.,

*F*is Hermitian, i.e., $F({\mathbf{q}}_{1},{\mathbf{q}}_{2})={F}^{\ast}({\mathbf{q}}_{2},{\mathbf{q}}_{1})$. Correlations of this type are used, for example, in quantum mechanics and quantum optics to describe mixed states, or in classical optics to describe partially coherent fields.

A useful representation of *F* is given in terms of its eigenvalues ${\Lambda}_{n}$ and eigenfunctions ${f}_{n}\left(\mathbf{q}\right)$, resulting from considering the integral Fredholm equation

*F*is Hermitian, its eigenvalues ${\Lambda}_{n}$ must be real, and its eigenfunctions ${f}_{n}\left(\mathbf{q}\right)$ can be made to be orthonormal and form a complete set:

*n*, and using Eq. (67),

*F*can be written as a diagonal sum known as a

*Mercer*or

*coherent mode expansion*[50]:

**q**, regardless of the functional form of the eigenfunctions ${f}_{n}$.

#### 5.2. Wigner Function and Spectrogram for Mixed States

The Wigner function for a partially coherent field or mixed state described by
*F* is given by

#### 5.3. Overall Degree of Purity

A measure of the overall level of coherence or purity of a field is given by [51]

*m*equally weighted coherent modes, this measure would be $\nu ={m}^{-1/2}$. Therefore, $1/{\nu}^{2}$ can be thought of as a measure of the effective number of significant coherent modes in

*F*[52]; $\nu =1$ implies that there is only one coherent mode, that is, that

*F*can be factored out as the product ${f}^{\ast}\left({\mathbf{q}}_{1}\right)f\left({\mathbf{q}}_{2}\right)$. The analog of Eq. (56) for mixed states is given by

#### 5.4. Quasi-homogeneous Fields

As will be discussed later, there is a type of optical field known as
*quasi-homogeneous*, for which the following factorization exists
for the correlation:

*μ*is dimensionless and must satisfy the properties

*μ*is significantly different from zero only over a range of values of its arguments of size ${q}_{\mathrm{c}}$ (the coherence width), then ${F}_{0}$ must have insignificant variation over changes of its arguments that are of the order of ${q}_{\mathrm{c}}$. Note that, for this type of field, the overall degree of purity can be written as

*μ*. Note that $\stackrel{\u0303}{\mu}$ is real, given the second relation in Eq. (77). Since

*μ*is much narrower than ${F}_{0}$, $\stackrel{\u0303}{\mu}$ is much wider than ${\stackrel{\u0303}{F}}_{0}$, so ${W}_{{F}_{\mathrm{QH}}}$ is significantly different from zero over a region of phase space much larger than the minimum uncertainty cell size, $|2K{|}^{-N}$.

#### 5.5. Smoothness of the Wigner Function for Highly Mixed States

As discussed in Subsection 4.6, the Wigner function of a deterministic function can be highly oscillatory and include many negative regions. However, for a significantly mixed state involving a very large number of coherent modes (i.e., for $\nu \ll 1$), the Wigner function tends to be smooth and non-negative. While not a rigorous proof, we now give an argument of why this is the case, based on the following two observations: first, the Wigner function of each coherent mode of a mixed state can oscillate between positive and negative values, although overall it is “more positive than negative” given the combination of Eqs. (46) and (66):

To illustrate this argument, consider a Gaussian correlation of the form

The overall degree of purity can be found either from the substitution of Eq. (84) into Eq. (73) and the use of the geometric series formula, or from the substitution of Eq. (83) into Eq. (75). The result is simply

so ${\sigma}_{2}={\sigma}_{1}$ corresponds to a fully coherent Gaussian beam. This result can also be found by substituting Eq. (83) into Eq. (75). Notice that ${\nu}^{-2}={\sigma}_{1}/{\sigma}_{2}$ equals the ratio between the cross-sectional phase-space area of the Wigner function for*F*, given in Eq. (83), to that for a fully coherent Gaussian function. That its, overall, the average relative variation of ${W}_{F}$ in any phase-space direction is smaller by a factor of

*ν*than that for the Wigner function of a coherent Gaussian. Media 12 (Fig. 14) shows (a) the Wigner function for the $n\mathrm{th}$ mode of this expansion given in Eq. (86), and (b) the partial sum up to

*n*of the superposition in Eq. (71). One can appreciate how the distribution becomes smoother as more modes are included.

It was argued in Subsection 4.6 that for
“pure” or “coherent” systems described by a function $f\left(\mathbf{q}\right)$, the use of the Wigner function can be problematic from a numerical
point of view, due to the fast oscillations and the doubling of the number of
variables. However, for “mixed” states this is not the case. First, the correlation
*F* depends on $2N$ variables, as does the Wigner function. Second, as discussed
earlier, the oscillations of the Wigner function tend to wash out for significantly
mixed states.

## 6. Other Phase-Space Representations

The Wigner function and the spectrogram are not the only representations of a function in terms of both its original variables and the Fourier conjugate ones. Extensive reviews of the many alternatives that have been proposed are given in some of the references mentioned earlier [1,8,10–14,54]. Here, I just give a brief summary of some of the standard representations and their properties.

#### 6.1. Distributions Related to the Spectrogram

The similarity between Eqs. (55) and (39), together with the shift property in Eq. (54), means that the spectrogram is related to the Wigner function according to

where the asterisk denotes convolution in both**q**and

**p**. That is, the spectrogram is a “blurred” version of the Wigner function, where the “blur” is the Wigner function of the window. This relation also holds for the spectrogram and Wigner function of a mixed function

*F*, since it holds for each of its coherent modes: One effect of this blurring is to smooth out the phase-space distribution, removing the fast oscillations and negative values sometimes presented by the Wigner function. In particular, it was mentioned earlier that for functions of the form $f\left(\mathbf{q}\right)=A\left(\mathbf{q}\right)exp\left[\mathrm{i}K\Omega \right(\mathbf{q}\left)\right]$, where

*A*varies slowly within the scale of the oscillations of the exponential, the Wigner function has a ridge structure around the manifold $\mathbf{p}=\partial \Omega /\partial \mathbf{q}$, but also presents high oscillations within the concavities of this ridge. For these functions, such oscillations are not present in the spectrogram, which reduces to a ridge localized around the manifold if an appropriate window function is used, as will be shown later.

As mentioned earlier, the window function used in the definition of the spectrogram
is sometimes dictated by the physics of the problem. This is the case, for example,
of measurement processes that are described by an overlap integral like that in
Eq. (55). Note that the definition
of the spectrogram in Eqs. (39) and
(38) has this same form, where the
measuring device is characterized by a wave function in the ${\mathbf{q}}^{\prime}$ space corresponding to $g\left({\mathbf{q}}^{\prime}\right)=\left\{\stackrel{\u02c6}{\mathcal{T}}\right(\mathbf{q},\mathbf{p}\left)w\right\}\left({\mathbf{q}}^{\prime}\right)\propto w({\mathbf{q}}^{\prime}-\mathbf{q})exp(\mathrm{i}K{\mathbf{q}}^{\prime}\cdot \mathbf{p})$. This “device” then makes a measurement of the process described by
*f*, focusing around the prescribed values **q** and
**p**, but respecting the limits established by the uncertainty relation.
For example, consider the measurement of a function of time through absorption by
resonant vibration modes of the measuring device (e.g., our hearing); in the simplest
case where the mode can be modeled as a damped harmonic oscillator, the window
function has the form $w({t}^{\prime}-t)\propto \Theta ({t}^{\prime}-t)exp\left[\gamma \right({t}^{\prime}-t\left)\right]$, where *Θ* is the Heaviside step function and
*γ* is the damping coefficient. For a detailed treatment of the use
of the spectrogram (and other distributions) to describe realistic measurements of
stochastic time signals, see Ref. [55].

In cases where the window function is not dictated by the physical situation, it can
be chosen for mathematical convenience. For example, it is well known that the
function that reaches the lower bound set by the uncertainty relation is a Gaussian,
whose Fourier transform is also a Gaussian. Therefore, a Gaussian window minimizes
the size of the blur in phase space. When the window is chosen as a normalized
Gaussian of width *a*, i.e.,

*Husimi function*[56]:

*Q function*[57–59]: This distribution corresponds to the overlap of a quantum density matrix associated with

*F*and what is known as a coherent state. An associated distribution also used in quantum optics is the

*Glauber–Sudarshan P function*[60–62], defined mathematically such that its convolution with the Wigner function of a Gaussian window of width ${a}_{0}$ gives the Wigner function of

*F*: It must be noted that, for pure states

*f*and even for many mixed states

*F*, the

*P*function is highly singular. It is therefore not often used for graphical representation over phase space of a function, but instead as a theoretical tool for derivations involving the coherent mode representation of quantized fields.

#### 6.2. The Kirkwood–Rihaczek and Margenau–Hill Distributions

Many other phase-space representations have been defined, including the Kirkwood–Rihaczek distribution [63,64], which for pure states is given by

*Q*functions) and the

*P*function, the Kirkwood–Rihaczek distribution (as well as its real part) satisfies marginal relations analogous to those satisfied by the Wigner function and given in Eqs. (43a) and (43b):

#### 6.3. The Ambiguity Function

Another commonly used representation that is defined in a
phase space is the ambiguity function [66–68]. However, this function is
quite different from those described earlier, since the phase space over which it is
defined is the Fourier conjugate of that for the Wigner function and other
representations. The ambiguity function is the *characteristic
function* for the Wigner function, i.e., its Fourier (or inverse Fourier)
transform over all variables:

The ambiguity function satisfies a phase-space mapping under linear canonical transformations similar to that of the Wigner function. To see this, consider the substitution of Eqs. (47) into the first expression in Eq. (97):

#### 6.4. The Cohen Class of Distributions

Except for the ambiguity function, all of the representations described so far are members of the so-called Cohen class of phase-space distributions [54,69,70], which in terms of the ambiguity function can be defined as

*ϕ*determines the representation. This kernel corresponds to the ratio of the characteristic function for the distribution in question to that for the Wigner function (i.e., the ambiguity function):

Since all Cohen-class representations are bilinear, they all present interference
contributions like those of the Wigner function. However, these effects can be
suppressed, amplified, or relocated depending on the kernel, as can be seen from
Eq. (101). Since the interference
contributions to the Wigner function are oscillatory, they contribute to the
ambiguity function at regions away from the origin in the $({\mathbf{q}}^{\prime},{\mathbf{p}}^{\prime})$ space. If the kernel $\varphi ({\mathbf{q}}^{\prime},{\mathbf{p}}^{\prime})$ decays away from the origin, as in the case of the spectrogram and
its special cases, the Husimi and *Q* functions, then these
interference terms are largely suppressed. For the *P* function, on
the other hand, for which the kernel grows extremely fast away from the origin, the
interference terms are greatly enhanced, causing the distribution to be singular
except in cases where ${A}_{F}$ decays sufficiently fast as to overcome the kernel’s divergence. In
the case of the Kirkwood–Rihaczek distribution, the kernel is a phase factor with
unit magnitude, so the interference effects are of similar significance as those for
the Wigner function, but shifted to a different region of phase space, as will be
discussed later.

Notice that, by using the convolution theorem, Eq. (101) states that all members of the Cohen class are given by the convolution over phase space of the Wigner function with a blur function, i.e.,

where the “blur”*Φ*is proportional to the Fourier transform of the kernel:

*P*function,

*Φ*is not well defined. The blurs for the different Cohen-class distributions are also presented in Table 1.

The inner-product property of the Wigner function in Eq. (55) is a special case of a more general relation for Cohen-class phase-space distributions, given by

*P*function is the dual of the

*Q*function, and the Kirkwood–Rihaczek distribution is the dual of its complex conjugate. Similarly, the moments of these different representations correspond to average values of the operators associated with the variables according to different ordering schemes, as discussed in detail in the reviews in Refs. [1,10,14].

To illustrate some of the distinctive features of these different representations,
let us consider the function $f\left(q\right)=exp(\mathrm{i}K{q}^{3}/12)$. As mentioned in Subsection 4.6, it would be expected that the significant values of a phase-space
distribution should be concentrated around the phase-space curve $p={q}^{2}/4$. Figure 15 shows (a) the
Wigner function, (b), (c) the Husimi function for two values of the width
*a*, (d) the magnitude, and (e) the real part of the
Kirkwood–Rihaczek distribution, all of which can be calculated analytically in terms
of Airy functions. The Wigner function presents the characteristic interference
contributions within the concavities of the curve. For the Husimi function, on the
other hand, these contributions are largely suppressed, except for some small effects
near the vertex of the parabola in 15(c),
because the width of the blur in *q* is not sufficiently narrow to
resolve this vertex properly. Figures 15(d) and
(e) illustrate the fact that it is the real part of the Kirkwood–Rihaczek
distribution (i.e., the Margenau–Hill distribution) and not its magnitude that holds
some resemblance to the curve. However, even the real part of this distribution
contains significant oscillatory interference contributions over large regions of
phase space. These regions correspond to the intersection of all values of
**q** occupied by $f\left(\mathbf{q}\right)$ and all values of **p** occupied by $\stackrel{\u0303}{f}\left(\mathbf{p}\right)$. A second example is given in Media
14–16 (Fig. 16), which show the transformation of several of these distributions under
fractional Fourier transformation for a sum of two Gaussians of the form $f\left(q\right)={\pi}^{1/4}\{exp[-{(q-4)}^{2}/2]+exp[-{(q+4)}^{2}/2\left]\right\}$. Note that, of all the representations shown in the figure, only
the Wigner and the *Q* function rotate rigidly in phase space with
*θ*. While the *P* function would also transform in
this way, it is highly singular for this case and therefore not plotted. Notice also
the counterrotating interference contributions in the case of the real part of the
Kirkwood–Rihaczek distribution.

Table 2 shows whether or not several of the
Cohen-class distributions transform as a simple mapping under different linear
unitary transformations. In the last column, the transformation $\stackrel{\u02c6}{\mathcal{G}}$ represents a general unitary transformation that is
*not* a combination of phase-space shifts and linear canonical
transformations.

Figures 15 and 16 illustrate the potentially very different
appearances of all the bilinear phase-space representations described so far for the
case of a pure function $f\left(\mathbf{q}\right)$. However, it is important to note that these differences tend to
wash away in the case of mixed functions $F({\mathbf{q}}_{1},{\mathbf{q}}_{2})$: the more incoherent or mixed the function is, the more similar all
these representations become. This is particularly easy to see in the
quasi-homogeneous limit, where *F* satisfies the factorization in
Eq. (76). For such fields, the
ambiguity function takes the form

*μ*and ${\stackrel{\u0303}{F}}_{0}$ are sufficiently sharply peaked, the approximation ${\stackrel{\u0303}{F}}_{0}\left({\mathbf{p}}^{\prime}\right)\mu \left({\mathbf{q}}^{\prime}\right)\varphi ({\mathbf{p}}^{\prime},{\mathbf{q}}^{\prime})\approx {\stackrel{\u0303}{F}}_{0}\left({\mathbf{p}}^{\prime}\right)\mu \left({\mathbf{q}}^{\prime}\right)\varphi (\mathbf{0},\mathbf{0})={\stackrel{\u0303}{F}}_{0}\left({\mathbf{p}}^{\prime}\right)\mu \left({\mathbf{q}}^{\prime}\right)$ is valid. That is, the dependence on the kernel

*ϕ*disappears, and the expression reduces to the result in Eq. (79):

*μ*are sharply peaked functions, ${F}_{0}$ and $\stackrel{\u0303}{\mu}$ are wide, slowly varying functions. Therefore, ${C}_{{F}_{\mathrm{QH}}}(\mathbf{q},\mathbf{p})$ varies slowly in all phase-space variables, occupying a region of phase space much greater than the minimum uncertainty cell. To illustrate this, consider the correlation function in Eq. (82). The Wigner function [already calculated in Eq. (83)], as well as the

*Q*,

*P*, and Kirkwood–Rihaczek functions, can be found to be given by where in the last expression $\nu =\sqrt{{\sigma}_{2}/{\sigma}_{1}}$ is the overall degree of purity. Notice that the

*P*function is only well defined if ${a}_{0}$ falls within the range ${\sigma}_{1}>{a}_{0}>{\sigma}_{2}$. In the limit ${\sigma}_{1,2}\to {a}_{0}$ (corresponding to a fully coherent Gaussian function of width ${a}_{0}$) the

*P*function tends to a delta function. On the other hand, in the highly incoherent case where ${\sigma}_{1}\gg {\sigma}_{2}$, the Kirkwood–Rihaczek function converges towards the Wigner function, as do also the

*P*and

*Q*functions if ${a}_{0}$ is chosen such that ${\sigma}_{1}\gg {a}_{0}\gg {\sigma}_{2}$. The convergence of all the Cohen-class representations occurs not only for functions satisfying the factorization in Eq. (76), but for any function $F({\mathbf{q}}_{1},{\mathbf{q}}_{2})$ that differs significantly from zero only for values of $|{\mathbf{q}}_{2}-{\mathbf{q}}_{1}|$ much smaller than the scale of variation of $F(\mathbf{q},\mathbf{q})$.

As mentioned in Section 1, the main topic of
this tutorial is the conservation of the Wigner function along paths associated with
a simple physical model. As discussed in Subsection 4.3, there are two general classes of linear unitary transformation on a
function *f* whose effect in phase space is a simple mapping
(i.e., rearrangement of arguments) of the Wigner function: linear canonical
transformations $\stackrel{\u02c6}{\mathcal{C}}\left(\mathbb{S}\right)$ and phase-space shifts $\stackrel{\u02c6}{\mathcal{T}}({\mathbf{q}}_{0},{\mathbf{p}}_{0})$. The most general class of unitary transformation that corresponds
to a mapping of the Wigner function is an arbitrary combination of these two types of
transformations. Therefore, the conservation property of the standard Wigner function
is restricted to systems where propagation or evolution is described by a combination
of phase-space shifts and linear canonical transformations. However, as will be seen
in Section 11, it is possible to retain the
conservation property for other physical systems by modifying the definition of the
Wigner function.

## 7. Wigner-Based Modeling of Beam Propagation

In optics and quantum mechanics, one of the main properties of the Wigner function is that, for a certain class of simple systems, its spatial propagation or temporal evolution follows the same rules as the simpler physical formalisms it mimics (i.e., geometrical optics or classical mechanics). In the case of paraxial classical optical propagation, this class of systems includes, but is not limited to, the so-called ABCD systems [71], which are characterized in both the ray and the wave models by a 4 × 4 symplectic matrix. In this section, we first give a brief review of the ray-optical and radiometric description of this type of system, and then show that, by using the Wigner function, the propagation of wave fields through these systems follows the same rules, namely, conservation along rays.

#### 7.1. Paraxial ABCD Systems in the Ray Domain

At a plane of constant *z*, a ray is fully characterized by its
intersection with this plane and by its direction. The intersection’s position is
specified by its transverse coordinate, $\mathbf{x}=(x,y)$, while direction is specified by the optical momentum two-vector, $\mathbf{p}=({p}_{x},{p}_{y})$, corresponding to the product of the local refractive index $n(\mathbf{x},z)$ and the transverse part of a unit vector in the ray’s direction
(see Fig. 17). The ray is then completely
characterized by the *state vector*
$(\mathbf{x},\mathbf{p})$, which gives the coordinates of the ray in phase space. The
paraxial approximation implies that $\left|\mathbf{p}\right|=nsin\theta \approx ntan\theta \ll n$, where *θ* is the angle between the ray’s tangent
vector and the *z* axis. For a simple class of optical systems, the
propagation of a paraxial ray is described mathematically by a simple relation
involving linear combinations of **x** and **p**. These systems,
known as linear systems or ABCD systems [71],
are fully characterized by a 4 × 4 symplectic matrix $\mathbb{S}$ composed of the 2 × 2 matrices $\mathbb{A}$, $\mathbb{B}$, $\mathbb{C}$, and $\mathbb{D}$, so that the initial state vector of the ray, $({\mathbf{x}}^{\prime},{\mathbf{p}}^{\prime})$, corresponding to a plane intersecting the optical axis at ${z}^{\prime}$, transforms into the final state vector at $z>{z}^{\prime}$, $(\mathbf{x},\mathbf{p})$, through multiplication by this matrix:

The fact that, in this context, $\mathbb{S}$ must be symplectic follows from a simple derivation of Eq. (111), based on the Hamiltonian
formalism of ray optics [72–74]. The part of an optical system contained
between the planes at ${z}^{\prime}$ and *z* can be described by Hamilton’s point
characteristic $V({\mathbf{x}}^{\prime},{z}^{\prime};\mathbf{x},z)$, which equals the optical length of an extremal path (i.e., a ray)
joining the points $({\mathbf{x}}^{\prime},{z}^{\prime})$ and $(\mathbf{x},z)$. (This description is complete provided that there is one and only
one extremal path joining these points; otherwise, a different type of characteristic
function must be used [74].) In the paraxial
regime and for a simple class of systems, sufficiently accurate results are obtained
from the quadratic expansion of the point characteristic. Further simplification
results from assuming that the *z* axis is an axis of symmetry of the
system, i.e., that the system has mirror symmetry around two perpendicular planes
intersecting at this axis. (This includes, as a special case, systems with rotational
symmetry around the *z* axis.) In this case, the quadratic expansion
does not include linear terms, and can be written, without loss of generality, in a
form analogous to Eq. (14):

**x**and

**p**: the first gives, straightforwardly, $\mathbf{x}=\mathbb{A}{\mathbf{x}}^{\prime}+\mathbb{B}{\mathbf{p}}^{\prime}$, as expected, while the second can be written, after eliminating

**x**from it by using the first relation, as $\mathbf{p}=\mathbb{C}{\mathbf{x}}^{\prime}+\mathbb{D}{\mathbf{p}}^{\prime}$, where $\mathbb{C}$ is defined precisely by Eq. (21), meaning that $\mathbb{S}$ is indeed symplectic.

To illustrate the type of system that can be modeled through a linear relation involving a symplectic ABCD matrix, consider the following three simple examples:

- (a) Propagation through a homogeneous transparent medium of index
*n*leaves the ray’s slope unchanged but changes the ray’s transverse position in an amount proportional to both the slope (approximately equal to $\mathbf{p}/n$ in the paraxial approximation) and the propagation distance $\Delta z$, so that the matrix is given by$${\mathbb{S}}_{\mathrm{HM}}\left(\Delta z\right)=\left(\begin{array}{cc}\hfill \mathbb{I}\hfill & \hfill \Delta z\mathbb{I}/n\hfill \\ \hfill \mathbb{O}\hfill & \hfill \mathbb{I}\hfill \end{array}\right).$$Notice that this matrix is similar to the one in Eq. (28). - (b) Propagation through a perfect thin lens of focal distance
*f*leaves the transverse position unchanged, but redirects the rays in an amount proportional to**x**, so that the matrix takes the form in Eq. (29): - (c) Propagation by a distance $\Delta z$ through a gradient index (GRIN) medium with refractive index $n\left(\mathbf{x}\right)={n}_{0}\sqrt{1-{n}_{2}^{2}|\mathbf{x}{|}^{2}}$ is given by the matrix$${\mathbb{S}}_{\mathrm{GRIN}}\left(\Delta z\right)=\left(\begin{array}{cc}\hfill \mathbb{I}cos\left({n}_{2}\Delta z\right)\hfill & \hfill \mathbb{I}sin\left({n}_{2}\Delta z\right)/{n}_{2}{n}_{0}\hfill \\ \hfill -\mathbb{I}{n}_{2}{n}_{0}sin\left({n}_{2}\Delta z\right)\hfill & \hfill \mathbb{I}cos\left({n}_{2}\Delta z\right)\hfill \end{array}\right).$$Like the matrix in Eq. (24), this matrix is, up to a scaling factor, a rotation matrix by an angle $\theta =-{n}_{2}\Delta z$.

#### 7.2. Radiometry and the Radiance in the Paraxial Regime

Ray optics gives a simple approximate description of the paths followed by optical
radiation, but not of the amount of power they carry. For this purpose, the theory of
radiometry supplements the rays with a weight, through the definition of the
*radiance* or *specific intensity*, which
corresponds to the amount of power per unit transverse area per unit solid
angle [75]. A closely related and perhaps
more fundamental quantity, known as the *basic radiance* [76], is equal to the radiance divided by the
square of the refractive index, i.e., to the amount of power per unit volume of phase
space:

*Φ*is the optical power, ${\partial}^{2}x=\partial x\partial y$ and ${\partial}^{2}p=\partial {p}_{x}\partial {p}_{y}$. When integrated over all the rays through a point, each specified by its direction

**p**, this quantity gives the

*optical intensity*or

*irradiance*at that point:

When losses at interfaces are negligible, the basic radiance is conserved under propagation along rays [75], and can therefore be regarded as a measure of the power density carried by the rays in the system. For an ABCD system, this conservation can be written as

Consider, for example, propagation through free space from an initial plane at ${z}^{\prime}$, for which $\mathbb{S}$ is as given in Eq. (114) with $n=1$ and $\Delta z=z-{z}^{\prime}$. In this case, this propagation law reduces to

#### 7.3. Paraxial ABCD Systems in the Wave Domain

In the wave domain, the propagation of a monochromatic paraxial scalar field from an
initial plane at ${z}^{\prime}$, traveling exclusively toward larger values of *z*,
can be described by a linear expression of the form

*k*is the free-space wavenumber, $V({\mathbf{x}}^{\prime},{z}^{\prime};\mathbf{x},z)$ is Hamilton’s point characteristic given by the optical path length of a ray joining the points specified by the arguments, ${\mathbb{M}}_{V}$ is the stability matrix

*z*, due to the changes of sign of $\mathrm{Det}\left({\mathbb{M}}_{V}\right)$ that occur at caustics [i.e., at places where the rays emanating from $({\mathbf{x}}^{\prime},{z}^{\prime})$ intersect].

The expression for the propagator in Eq. (127) was first proposed in optics by Walther [77,78]. An analogous
propagator was proposed earlier by Van Vleck [80] to describe the evolution of solutions to Schrödinger’s equation in
terms of classical particle trajectories, where instead of *V* one
uses the classical action. This result can be considered as a form of the
Huygens–Fresnel principle: the field at the initial plane can be decomposed into
secondary sources, each emitting a field proportional to $\mathcal{K}({\mathbf{x}}^{\prime},{z}^{\prime};\mathbf{x},z)$ whose phase delay at the observation point is proportional to the
optical distance from the secondary source, and whose optical intensity is
proportional to the density of rays around the observation point arriving from the
secondary source point [i.e., $\mathrm{Det}\left({\mathbb{M}}_{V}\right)$]. Note that the original expressions by Walther and Van Vleck did
not include the Maslov index correction, which was added by Gutzwiller [81] within the quantum context.

For situations where the ABCD formalism is valid in the ray-optical regime, sufficiently accurate results are obtained from using in Eq. (127) the quadratic expansion of the point characteristic in Eq. (112). It is straightforward to see that, in this approximation, the propagated field is proportional to a linear canonical transformation as defined in Eqs. (17) and (14) with $K=k$, applied to the initial field:

#### 7.4. Wave Propagation Through ABCD Systems in Terms of the Wigner Function

The Wigner function of the field $U(\mathbf{x},z)$ over the transverse variables gives a representation in terms of
both transverse position **x** and transverse momentum **p**:

To illustrate these ideas, let us consider the simplest physical situation,
corresponding to propagation through a region of free space surrounding the plane $z=0$. The substitution in Eq. (124) of ${z}^{\prime}=0$ and with *L* replaced with ${W}_{U}$ in Eq. (132) gives
the expression

*y*, for which the phase space $(x,p)$ is two-dimensional. In this figure, the initial field ${U}_{0}\left(x\right)=U(x,0)$ is chosen to be constant for $\left|x\right|\le a/2$ and zero otherwise. Within the paraxial regime, this field corresponds, for example, to a plane wave traveling in the

*z*direction that is diffracted a slit of width

*a*. The Wigner function, which can be calculated analytically, is shown in Fig. 18(b). Its marginal in

*p*does return the desired rectangular distribution. As shown in Media 17, propagation in

*z*is described by a horizontal linear shearing of the Wigner function, whose marginal in

*p*gives the intensity profile at the corresponding propagation distance. This example shows the importance of the negative values of the Wigner function in order to account for destructive interference. It also gives a simple visualization of the fact that, at the propagation distance ${z}_{\mathrm{peak}}=0.0545k{a}^{2}$, a peak in the axial intensity is achieved, associated with the vertical alignment of some of the main positive regions of the Wigner function. This is precisely the distance at which, for example, a pinhole camera with a square pinhole of side

*a*would achieve its best resolution. (For a circular pinhole of diameter

*a*, the numerical factor is slightly smaller, 0.044, as found by Rayleigh.)

An alternative presentation of the same example is given in Fig. 19. Figure 19(a) shows the intensity distribution over the (scaled) physical plane. Figure 19(b) shows a superposition of lines representing rays, which are either black or white depending on whether the Wigner function for the corresponding ray is negative or positive, and whose opacity is proportional to the magnitude of the Wigner function, such that rays for which the Wigner function is close to zero are nearly invisible. These lines are superposed over a grey background (representing zero), so that the final shade of grey at each point depends on the number and weight of the positive and negative rays going through it. Media 18 shows partial superpositions of these rays for the sake of illustration. It is seen there that, once a sufficient number of rays is included, the ray superposition in Fig. 19(b) reproduces the features of the wave pattern in Fig. 19(a).

It must be noted that ABCD systems are not the most general type of system for which the rules of paraxial propagation for the Wigner function exactly mimic the radiometric model. As discussed in Section 4, the most general type of transformation that corresponds to a simple mapping of the Wigner function in phase space is a combination of linear canonical transformations and phase-space shifts. This corresponds to systems where Hamilton’s characteristic function can be approximated by a quadratic expansion like that in Eq. (112), but where this expansion can include linear terms that account for asymmetries of the system. These linear terms lead to phase-space shifts. Lateral spatial shifts can be introduced, for example, by a tilted étalon, or by a deliberate lateral shift of the optical axis in order to model a misaligned optical element. Directional shifts can be introduced by, say, a transparent wedge, or by a deliberate small rotation of the optical axis in order to model a tilted optical element.

#### 7.5. Phase-Space Interpretation of Measurement

As discussed in [45], an ideal measuring
device can be thought of as a source that emits a characteristic field whose
amplitude and phase are such that it cancels as much of the incident field as
possible. (For example, a point detector that measures the local intensity can be
thought of as a point source emitting a spherical wave.) Let us consider the paraxial
case. In general, let the field specific to some detector be given, at a plane
*z* (after the detector), by the function $g\left(\mathbf{x}\right)$, times an undetermined constant complex amplitude, ${A}_{0}$. The power absorbed by the detector is then the difference between
the power of the incident field *U* and that of the field past the
detector, $U-{A}_{0}g$. This absorbed power is given by

Of course, this interpretation should not be taken literally, since rays are not real, independent physical objects, but rather they are just a mathematical construct. Further, the Wigner function of the field can be negative for some thin bundles of rays. However, due to the uncertainty relation, the Wigner function of the detector cannot be arbitrarily narrow, so that it cannot be focused on arbitrarily narrow bundles of rays of the field. In particular, it cannot focus on narrow bundles of rays over which the field’s Wigner function is purely negative; a sufficient amount of “positive rays” is always also accepted into the detector, so that the result of the measurement is always non-negative. Something analogous happens with possible negative regions of the Wigner function for the measuring device.

#### 7.6. Applications of the Wigner and Ambiguity Functions in the Paraxial Regime

In the past three decades, a large number of articles have been published whose topic is the application of the Wigner and ambiguity functions to the study and visualization of optical systems within the paraxial regime. It is impossible to do justice to all of this body of work within the limited space provided by this article, and I apologize in advance to the many authors whose interesting work is not mentioned or cited here. Since the main emphasis of this article is the use of the simple conservation properties of the Wigner function, the applications singled out in this section are ones that exploit these properties.

### 7.6a. Imaging Systems and the Optical Transfer Function

The expression for the intensity of a field propagating through free space in terms of the Wigner function is given in Eq. (135). The corresponding expression in terms of the ambiguity function is given by

*et al.*[88], for the field generated by a point source in object space, the ambiguity function at “radial” planar sections defined by ${\mathbf{x}}^{\prime}=-z{\mathbf{p}}^{\prime}$ corresponds to the OTF of the system for different amounts of defocus

*z*, as can be seen from Eq. (140). Based on this realization, both the Wigner and ambiguity functions have been employed in the design of optical systems, e.g., to study axial resolution [89] or to image transparent objects with partially coherent illumination [90].

### 7.6b. Propagation-Invariant Beams and Extended Depth of Field

There is a class of solutions of the paraxial wave equation whose intensity profile is exactly invariant under propagation. These solutions do not correspond to realistic fields, since they require an infinite amount of optical power. However, they can be achieved approximately, leading to fields whose transverse intensity profile is nearly invariant under propagation within some range of distances. Perhaps the best known such solution is that of Bessel beams [91], given by the expression

**x**and the

*x*axis, ${J}_{n}$ is an

*n*th order Bessel function of the first kind,

*n*is the vortex charge of the beam, and ${p}_{0}$ is a constant that sets the numerical aperture of the beam. In general, all propagation-invariant beams have a transverse Fourier transform (with $K=k$) that is restricted to a ring of radius ${p}_{0}$, i.e.,

*A*is an arbitrary function of

**p**. The special case of Bessel beams corresponds to $A\left(\mathbf{p}\right)=exp\left(\mathrm{i}n{\varphi}_{p}\right)$, where ${\varphi}_{p}$ is the angle between

**p**and the ${p}_{x}$ axis. By using the definition of the Wigner function in terms of the Fourier transform, one can easily find

*Θ*is the Heaviside function and $\stackrel{\u0304}{\mathbf{p}}\left(\mathbf{p}\right)$ is a vector perpendicular to

**p**, defined as

**x**is in a dot product with $\stackrel{\u0304}{\mathbf{p}}$, which is perpendicular to

**p**; so replacing

**x**with $\mathbf{x}-z\mathbf{p}$ has no effect. This invariance under horizontal shearings, though, implies that the Wigner function extends over all values of

**x**, meaning that its integral is infinite. This is consistent with the fact that propagation-invariant fields have infinite power.

There is an interesting second class of propagation-invariant field, whose
intensity profile shape is invariant under paraxial propagation, although it does
not travel along a straight line but along a parabolic path. The archetypal beam
of this type is the Airy beam [93]. Let us
consider first the version of this beam for a space with only one transverse
variable, *x*, in addition to the propagation distance,
*z*. The formula for an Airy beam is then given by

*a*is a parameter that regulates both the width of the lobes of the beam’s intensity as well as the path these lobes follow. The intensity pattern of this beam, shown in Fig. 20(a), is given by

*z*up to a lateral displacement proportional to ${z}^{2}/\left({k}^{2}{a}^{3}\right)$. As in the case of the Bessel beams, the Wigner function is easier to calculate from the field’s transverse Fourier transform at $z=0$, given by a simple cubic phase:

*z*is equivalent to a phase-space shift. By substituting $x\to x-zp$ in Eq. (147), it is easy to show that this phase-space shift is given by $\stackrel{\u02c6}{\mathcal{T}}[{z}^{2}/(4{k}^{2}{a}^{4}),z/(2{k}^{2}{a}^{3}\left)\right]$; the spatial part of the shift is responsible for the curvature of the beam’s path, while the directional part leads to the phase factor in Eq. (144). By using the projection properties of Airy functions [94], one can show that the marginal in

*p*of this Wigner function equals the intensity distribution in Eq. (145). Note that the intriguing feature of Airy beams of following a curved trajectory becomes less mysterious under this phase-space picture: all rays composing this beam travel along straight lines, but the rays that give rise to a maximum at a given value of

*z*are not the same ones as those that give rise to the maximum at a different propagation distance. An analogous interpretation can be given in terms of caustics [93]. A treatment that is similar to the one presented here was published recently [95].

A three-dimensional version of the Airy beam can be composed as the product of the
field in Eq. (144) times a
similar solution in terms of *y* [divided by ${u}_{0}exp\left(\mathrm{i}kz\right)$], and the corresponding Wigner function is proportional to the
products of the Wigner functions in Eq. (147) with arguments $(x,{p}_{x})$ and $(y,{p}_{y})$, respectively. Approximations to these beams can be
generated [96,97] by placing at the pupil of an optical system a transparent
mask whose phase has a cubic dependence [as in Eq. (146)] in both Cartesian directions. Coincidentally, this
cubic phase mask was proposed by Dowski and Cathey [98] as a way to extend the depth of field of imaging systems.
With such a cubic phase, the PSF corresponds approximately to the intensity of an
Airy beam, which remains nearly invariant for a range of defocus distances.
Further, given this invariance and known profile, this PSF can be deconvolved
computationally, roughly independently of defocus. Of course, there is a small
deviation of the defocused PSFs with respect to the ideal geometrical location due
to the curvature of the path followed by the maximum, but this deviation is
negligibly small in imaging applications. Approximations to Bessel beams,
particularly those of zeroth order, have also been used to extend depth of
focus [99–101]. These are achieved, for example, by the inclusion of
either annular pupils or axicons within the imaging system. Other studies of
extension of the depth of field of an imaging system based on phase-space
distributions are given in [102,103].

### 7.6c. The Talbot Effect

The ray-like propagation of the Wigner function also provides a simple, visual
description of the Talbot effect. In this effect, fields that are periodic over
the initial transverse plane form images of themselves at certain distances of
propagation in free space, without the use of lenses or mirrors. Consider again,
for simplicity, a two-dimensional field (independent of *y*). Let
the field at the plane $z=0$, ${U}_{0}\left(x\right)$, be periodic with period *Λ*, i.e., ${U}_{0}(x+\Lambda )={U}_{0}\left(x\right)$ for any *x*. The Fourier transform of ${U}_{0}$ is then a discrete sum of weighted Dirac deltas:

*p*, differing from zero only when

*p*is an integer multiple of $\pi /k\Lambda $, as shown in Fig. 21. (Again, for a realistic field, these delta functions must be replaced by narrow finite distributions.) Notice that this spacing in

*p*is half of that between the delta functions in ${\stackrel{\u0303}{U}}_{0}\left(p\right)$ in Eq. (148), due to the presence of half-point interference contributions in the Wigner function. At the discrete set of values of

*p*where it differs from zero, the Wigner function has different periodicity in

*x*depending on whether

*p*is an even or odd multiple of $\pi /k\Lambda $; it can be seen from Eq. (149) that

*n*. This can be appreciated in Fig. 21, where the values of the Wigner function at even rows has periodicity $\Lambda /2$, while for odd rows it oscillates between positive and negative values with period

*Λ*.

Given its periodicity in *x* and discreteness in
*p*, one can see that the Wigner funcion is also periodic under
propagation in *z*, due to the linear shear relation

*Λ*. Since the Wigner function vanishes unless

*p*is an integer multiple of $\pi /k\Lambda $, then the Wigner function returns to its original configuration when $z={z}_{n}=nk{\Lambda}^{2}/\pi $ (the Talbot distances) where

*n*is an integer. This is illustrated by Media 20 (Fig. 21). As discussed in Refs. [104,106], this phase-space picture also helps visualize the so-called fractional Talbot effect. Consider, for example, the case when $z={z}_{1}/2=k{\Lambda}^{2}/2\pi $. Given the phase-space shearing, all even rows of the Wigner function, corresponding to $p=2\pi n/k\Lambda $, shift by an integer number of periods, returning to their state at $z=0$. On the other hand, the odd rows, corresponding to $p=(2n+1)\pi /k\Lambda $, shift by a half-integer number of periods. Given the periodicity properties in Eq. (150), the resulting Wigner function is identical to the initial one except for a global shift in

*x*of $\Lambda /2$. This explains why, at this propagation distance, the intensity pattern is identical to that at $z=0$ except for a half-period shift. Extensions of self-imaging to polychromatic fields have also been studied based on the Wigner function [107].

Note that this phase-space picture of the Talbot effect has several features in common with that of the propagation-invariant beams discussed in Subsection 7.6b. In both cases, the Wigner function after horizontal shearing is identical to the initial one (or a shifted version of it). The only difference is that, in this case, this reproduction of the initial state happens only at a discrete set of propagation distances. Nevertheless, this similarity helps explain other common features, like the one referred to as ”self-healing.” Consider the case where a section of the initial field (a Bessel or Airy beam, or a periodic field) is perturbed somehow, e.g., through blocking by an obstacle. After some propagation, this perturbation washes away, so that the beam returns to its original unperturbed shape. As in the case of Bessel and Airy beams, local features of the field after propagation should not be thought of as being the same physical entity as the corresponding features of the initial field, even if they are identical in shape. Such misinterpretation can lead to unphysical conclusions, like the violation of special-relativistic causality.

### 7.6d. Sampling and Holography

The Wigner function also gives a useful visualization of why the well-known
sampling theorem can be relaxed for chirped signals [108,109], e.g., for
paraxial beams that have propagated away from their waist plane. Shannon’s
sampling theorem states that the maximum sampling spacing over the transverse
variable *x* that allows a complete reconstruction of the field is
inversely proportional to the support of the Fourier transform of the field at the
initial plane. Paraxial free propagation does not change this support, given the
relation $\stackrel{\u0303}{U}(\mathbf{p},z)={\stackrel{\u0303}{U}}_{0}\left(\mathbf{p}\right)exp\left[\mathrm{i}kz\right(1-|\mathbf{p}{|}^{2}/2)]$. One could therefore think that the maximum sampling spacing
remains constant under propagation. Gori [110] showed that, for a field that is strictly limited to an interval,
say, at $z=0$, the maximum transverse sampling spacing needed to fully recover
the field increases linearly with *z*. Similarly, for fields whose
Fourier transform ${\stackrel{\u0303}{U}}_{0}\left(\mathbf{p}\right)$ is (at least approximately) restricted to a region $\left|\mathbf{p}\right|<{p}_{\mathrm{M}}$, the sampling frequency needed to approximately reproduce the
field can be relaxed under propagation. A simple interpretation of this relaxation
in terms of the Wigner function follows a derivation that is mathematically
similar to that mentioned earlier for the Talbot effect, but with
*x* and *p* exchanged. Consider sampling a
two-dimensional field at a discrete number of equidistant points over the initial
line ($z=0$), i.e., ${U}_{\mathrm{s}}\left(x\right)={\sum}_{n}U(x,z)\delta (x-n\Lambda )$, where *Λ* is the sampling distance. The Wigner
function of this sampled field can be found to be given by [108,109]

*p*direction at a separation of $\pi /k\Lambda $, which is half the separation between the corresponding replicas of the Fourier transform. This is because the replicas corresponding to odd ${n}^{\prime}$ are due to half-point interference effects. Given the alternating sign in front of them, these intermediate replicas can be cancelled through a suitable projection, so they are not considered in the schematic analysis that follows. Figure 22(a) shows schematically the replicas corresponding to even ${n}^{\prime}$, separated vertically by $2\pi /k\Lambda $. To be able to reconstruct the original field, this separation must be sufficient so that the different replicas do not overlap. If the field is not chirped, this condition reduces to $2\pi /k\Lambda \ge 2{p}_{\mathrm{M}}$, where $2{p}_{\mathrm{M}}$ is the support of the angular spectrum. This relation can be written as $\Lambda \le \pi /k{p}_{\mathrm{M}}$, which is the standard sampling theorem. If the field is chirped, e.g., due to propagation away from its waist plane, its Wigner function is sheared, and the replicas can be stacked more densely without overlapping significantly, relaxing the minimum value of

*Λ*.

The revised sampling theorem described earlier has been applied to digital holography [111,112]. Similarly, the Wigner function has been used to explain the principles of regular [113–116] and volume holography [117]. Other recent applications of the Wigner function include descriptions of the moiré effect [118] and the detection of subsurface targets closely located beneath a randomly rough surface [119].

#### 7.7. Partially Coherent Fields

The fact that paraxial propagation through ABCD systems corresponds to a simple remapping of the Wigner function makes this formalism a useful numerical tool in the case of partially coherent fields. In the frequency domain, a partially coherent field is described by its cross-spectral density, defined as the correlation of the field at two points:

*W*, but we denote it here by

*J*instead to prevent confusion with the Wigner function.) At any plane ${z}_{1}={z}_{2}=z$, this function depends on four variables (${\mathbf{x}}_{1}$ and ${\mathbf{x}}_{2}$). At the same plane, the Wigner function also depends on four variables (

**x**and

**p**):

*ν*), the smoother its Wigner function tends to become; so its sampling over phase space can be relaxed. Given this, if, say, the intensity of a field needs to be computed at a large number of points $\mathcal{N}$, traditional calculations based on propagation integrals would require the computation of a quadruple integral of an oscillatory integrand for each point. The computation time would then be proportional to ${D}^{4}\mathcal{N}$, where

*D*is the number of samples required in each variable for the integration. On the other hand, if the Wigner function is used, there is an initial computational investment in calculating this function at $m\times m\times m\times m$ sample points over the four-dimensional phase space, each requiring the evaluation of a double oscillatory integral that involves ${M}^{2}$ operations. Once this calculation is completed and its results stored in memory, propagation corresponds to a simple remapping of the Wigner function, and the evaluation of the intensity at any point involves only integration over the directional variables, i.e., ${m}^{2}$ operations. Therefore, the total computation time is proportional to ${m}^{4}{M}^{2}+{m}^{2}\mathcal{N}$, versus ${D}^{4}\mathcal{N}$ for the traditional approach. For fairly incoherent fields,

*M*is comparable with

*D*, and both are significantly larger than

*m*, so the Wigner-based computation becomes advantageous when $\mathcal{N}$ is sufficiently large [120]. In both the traditional and Wigner approaches, computational savings can be achieved by using fast Fourier transform algorithms or asymptotic integral estimations, whenever appropriate. For partially coherent fields involving only a reduced number of coherent modes, however, it might be more efficient to use a modal decomposition and propagate each mode independently.

While the Wigner function gives direct access to one-point field properties such as the optical intensity, its ray-like propagation properties have also been exploited for the propagation of two-point quantities like the cross-spectral density, through the inversion of Eq. (154) [121,122].

#### 7.8. Wigner Functions that Account for Polarization in the Paraxial Regime

Wigner functions have also been proposed to model polarization effects in the paraxial regime. These functions take the form of a 2 × 2 matrix [123,124], corresponding to the substitution of the 2 × 2 cross-spectral density matrix into Eq. (154):

*x*or

*y*. Alternatively, Luis [125,126] defined a Stokes–Wigner function that takes the form of a 4-vector whose elements are a generalization of the Stokes parameters, calculated as the trace of the product of ${\mathbb{W}}_{\mathbb{J}}$ and each of the four Pauli matrices

**p**of this Wigner function gives the four Stokes parameters at the corresponding point.

## 8. The Use of the Wigner Function for More Complicated Systems and/or Nonparaxial Fields

The simple propagation properties of the Wigner function are limited to paraxial systems in which rays evolve according to linear relations. For more complicated systems involving optical elements introducing phase delays whose dependence on the transverse position are more complicated than quadratic (e.g., lenses with aberrations), propagation does not correspond to a simple rearrangement of the arguments of the Wigner function; it must be described instead either through an integral propagation kernel [70] involving all phase-space variables, or by a series of differential operators [10,14,41] acting on the Wigner function. This is also the case for diffraction by planar obstacles, where one can use an integral relation to model diffraction [127,128]. In the case of fairly spatially incoherent fields, a radiometric-like model can be used, where the primary effect of the obstacle is to block the rays incident on it (i.e., the value of the Wigner function for these rays is set to zero), while those missing the obstacle continue to propagate (i.e., the value of the Wigner function for these rays remains unchanged). Additionally, the edges of the obstacle behave as secondary sources emitting rays [129], whose Wigner function is related to that of the rays incident on the edges. Approximate paraxial propagation models based on the Wigner function have also been proposed for the case of turbulent media [130–133], as well as for nonlinear Kerr media [134–138].

Even for propagation in free space (or a homogeneous transparent medium), the Wigner function is not conserved along rectilinear paths for fields that are not paraxial. To show this, let us use the angular spectrum representation of a field traveling in free space:

**x**of ${U}_{0}\left(\mathbf{x}\right)=U(\mathbf{x},0)$, and ${p}_{z}\left(\mathbf{p}\right)=\sqrt{1-|\mathbf{p}{|}^{2}}$. Notice that, even though the field is nonparaxial, we still choose to treat

*z*as a propagation parameter, separate from the other two spatial coordinates. Let us assume that ${\stackrel{\u0303}{U}}_{0}\left(\mathbf{p}\right)=0$ for $\left|\mathbf{p}\right|>1$, so that the field contains no evanescent components. The Wigner function for this field at any plane of constant

*z*can be written in terms of either

*U*or ${\stackrel{\u0303}{U}}_{0}$:

#### 8.1. Generalized Radiometry for Nonparaxial Fields

The first use of the Wigner function in classical optics was within the context of nonparaxial partially coherent fields. The objective of such work was to give a wave-optical foundation to the heuristic theory of radiative transfer [140,141] and in particular to the radiometric formalism [6]. These efforts concentrated on deriving a wave-based radiance. This was first done, independently, by Dolin in 1964 [4] and by Walther in 1968 [5]. In particular, Walther proposed the use of the standard Wigner function over the transverse variables, and introduced an extra obliquity factor ${p}_{z}\left(\mathbf{p}\right)=\sqrt{1-|\mathbf{p}{|}^{2}}$ to account for nonparaxiality of the field:

*generalized radiances*, as these wave-based definitions of the radiance are referred to in the literature, have been proposed. Many articles on this topic are collected in the volume edited by Friberg [142], and useful summaries on the subject are given in the articles by Wolf [143] and Friberg [144]. Many of these generalized radiances result from the use of other members of the Cohen class [145]. For example, Walther proposed a second generalized radiance that turns out to be the Kirkwood–Rihaczek distribution times an obliquity factor [146].

There are two objections to the wave-based definition of the radiance in Eq. (159): it can be negative, and it is not rigorously conserved along rays beyond the paraxial regime. However, it is easy to show that, in the limiting case of significantly incoherent fields of the quasi-homogeneous type, these two problems disappear [145,147–148]. Consider a field whose cross-spectral density at a plane $z=0$ has the form

*z*axis. This relation holds not only for the factorizable quasi-homogeneous fields in Eq. (160) but also for any field where the directional correlation ${\stackrel{\u0303}{\stackrel{\u0303}{J}}}_{0}({\mathbf{p}}_{1},{\mathbf{p}}_{2})$ differs significantly from zero only for separations of its two arguments that are much smaller than unity.

As discussed in Subsection 6.4, any definition of the generalized radiance based on a standard Cohen-class distribution would lead to the same result in this limit (including Walther’s second generalized radiance definition [146]). Light generated by thermal sources typically presents coherent widths that are small compared with the scale of variation of the intensity, so any Cohen-class generalized radiance can be used to justify the radiometric model for the case of thermal light.

## 9. Short Pulses and Propagation in Linear Dispersion Media

Another area of optics where phase-space distributions have been applied is the
characterization of pulses, and the description of their propagation through dispersive
media [149–153]. Most work on this topic has been limited to a one-dimensional
treatment, useful for studying highly collimated fields traveling through homogeneous
media or propagation modes in waveguides, particularly optical fibers, where the
approximately invariant transverse structure of the field can be factorized from the
more interesting longitudinal/temporal behavior. The vast majority of this work is based
on the use of the phase space composed of time and frequency [149,151,154,155],
often referred to as the *chronocyclic phase space*. However, some
authors consider instead the phase space composed of the propagation distance and
wavenumber [150,152,153].

As discussed in this section, within some approximations, the propagation of pulses
through certain optical elements turns out to be mathematically analogous to paraxial
beam propagation through ABCD systems [156,157]. Useful reviews of this analogy and its
applications to the characterization of ultrashort pulses is given in the book chapter
and review article by Dorrer and Walmsley [154,155]. For simplicity, most of the
treatment below considers single pulses, described by a component of the electric field, $E(z,t)$, where *z* is the spatial propagation variable.
However, for the characterization of pulsed lasers used for ultrafast applications, it
is convenient to use a statistical formalism analogous to that of the theory of partial
coherence. Here, each pulse in the train is treated as a member of an ensemble. Since
the different pulses are generally not identical, the pulse train is characterized by
the correlation function [158,159]

Within the linear domain and neglecting the effects of absorption, a dispersive medium
or system accepts monochromatic solutions of the form $exp\left[\mathrm{i}\right(kz-\omega t\left)\right]$. The temporal frequency *ω* and wavenumber (or spatial
frequency) *k* are related through the dispersion relation

*c*is the speed of light in vacuum. The refractive index

*n*is assumed to be purely real, although in reality this is not possible (except in free space) due to the constraints imposed by causality and expressed mathematically through the Kramers–Kronig relations. However, for many cases of practical interest, absorption is negligible over the range of frequencies occupied by the field, so its effects can be neglected. The plot of

*ω*versus

*k*is known as the dispersion curve.

A pulse traveling in a dispersive medium of this type can be expressed as a superposition of monochromatic components:

*E*at $z=0$:

*ω*versus

*k*plane, the propagation/evolution of the pulse is mathematically analogous to the free-space propagation of a paraxial beam. These two parabolic approximations lead to similar but nonequivalent treatments, as we now discuss.

#### 9.1. Wigner Function of Time versus Frequency

Consider the Taylor expansion of *k* up to second order around the
central frequency of the pulse, ${\omega}_{0}$, i.e.,

*Γ*is the correlation defined in Eq. (164).] From the substitution of Eq. (166) into Eq. (170), it is easy to show that, if the approximation in Eq. (168) is valid, the Wigner function propagates as where, given Eq. (168), ${k}^{\prime}\left(\omega \right)={k}_{2}\omega +{k}_{1}-{k}_{2}{\omega}_{0}$ is a linear function of frequency. That is, propagation in

*z*is described within this approximation by a simple linear shearing of the Wigner function over the phase space $(t,\omega )$. One way to understand this follows from writing Eq. (166) in the following form:

*z*retains its role as a propagation parameter and

*τ*is analogous to the transverse position variable. As the last step in this equation states, propagation in

*z*corresponds to a linear canonical transformation in

*τ*with a matrix analogous to that in Eq. (114), i.e.,

As discussed in Refs. [154,155], devices such as quadratic temporal phase modulators introduce chirp factors and hence play a role analogous to lenses; i.e., they correspond to linear canonical transformations with matrices analogous to that in Eq. (115). The combination of dispersive media and quadratic temporal phase modulators then allows the composition of general “temporal ABCD systems” that can perform transformations such as imaging, magnification, and Fourier and fractional Fourier transformation.

#### 9.2. Wigner Function of Position versus Wavenumber

Consider instead the parametrization of *ω* in terms of
*k*, and assume that, over the range of spatial frequencies
occupied by the pulse, a quadratic Taylor expansion of the frequency around ${k}_{0}$ (the central wavenumber) is sufficient:

*z*of the field:

*t*then corresponds to a linear shearing in the phase space $(z,k)$, since ${\omega}^{\prime}\left(k\right)={\omega}_{2}k+{\omega}_{1}-{\omega}_{2}{k}_{0}$ is linear in

*k*. Again, this is because the field can be written as a solution to a Schrödinger-like equation, although on different variables. Let us again factor out the carrier phase and express the remaining function in terms of regular and retarded times,

*t*then corresponds to a Fresnel-type linear canonical transformation in

*τ*of the form

## 10. Phase-Space Tomography

Phase-space distributions like the Wigner and ambiguity functions have been used for the recovery of the coherence properties of mixed states from measurements of the marginal projections through a set of techniques known as phase-space tomography, first proposed by Bertrand and Bertrand [35]. This approach works particularly well for systems where phase space is two-dimensional and for which propagation or evolution corresponds either to a Fresnel-type integral (i.e., to a linear shearing of phase space) or to a fractional Fourier transform (i.e., to a rotation of phase space). The idea is that, if sufficient marginal projections are obtained through measurements at different stages, the Wigner function can be estimated through standard tomographic methods like the inverse Radon transform. Once the Wigner function is known, the correlation function $F({q}_{1},{q}_{2})$ can be recovered. Much of the pioneering work on this type of technique came from the same research group [36–38], who applied it to three different physical contexts, as discussed in what follows.

#### 10.1. Phase-Space Tomography in Quantum Optics

The first experimental implementation of phase-space tomography was by Smithey
*et al.* [38] (following
work by Vogel and Risken [160]), who used
this technique to recover of the quantum state of a light field. To let us discuss
this technique, I give a very brief description of the use of the Wigner function in
quantum optics. To make the presentation as similar as possible to that in the rest
of this article, the standard quantum notation of operators and density matrices is
avoided. For more standard treatments of the topic, the reader can consult the books
by Mandel and Wolf [50], Gerry and
Knight [161], Leonhardt [39], and Schleich [162].

In classical optics, the electric field is the “wave function” of interest, and its variables are position and/or time. Consider a Cartesian component of the electric field at a fixed point for a monochromatic mode within a cavity. In the classical picture, this component can be written as

where*ω*is the frequency of oscillation. Note that we are using the complex (or analytic signal) representation of the field. Let us define the normalized quadratures

*α*is a constant (inversely proportional to the square root of $\u0127\omega $) that guarantees that ${q}^{2}+{p}^{2}$ gives, instead of the intensity, the total energy of the mode divided by $\u0127\omega $, i.e., the number of photons in the mode. Classically,

*q*and

*p*can be fully determined simultaneously, so they can be represented by a single point in the classical phase space $(q,p)$, as shown in Fig. 23(a). The angular coordinate of this point in this plane corresponds to the phase of the field, while the radial one is proportional to the amplitude, and therefore related to the intensity or number of photons. As time advances, this point rotates around the origin along a circular path, with frequency

*ω*. (Note that this phase-space motion is mathematically analogous, although entirely different physically, to that for a ray propagating paraxially in a quadratic GRIN medium, where

*q*is the transverse position,

*p*is the transverse optical momentum, and ${n}_{2}z$ replaces $\omega t$.)

In the quantum-optical description of the field, on the other hand, the quadrature
*q* is not a simple function of time, but the variable of a wave
function $\Psi (q;t)$. Further, *q* and *p* cannot be
prescribed simultaneously, since they are Fourier conjugate variables (with $K=N=1$). However, by using a phase-space representation like the Wigner
function, we can represent the state of the field at any given time as a function of
both *q* and *p*, shown in
Fig. 23(b). That
is, at a fixed time, instead of having a delta-like localized distribution, the field
is represented by a distribution ${W}_{\Psi}(q,p;t)$ of “possible quadrature pairs” or “possible amplitudes and phases”
whose extent in phase space is limited by the uncertainty relation: the narrower the
distribution is in *q*, the wider it is in *p*, and
vice versa. Since we have, instead of a single point, an extended distribution, the
field does not have well-defined amplitude and phase: the smaller the uncertainty in
the phase (related to the angular spread of the distribution with respect to the
origin) of the field is, the larger the uncertainty in the number of photons (related
to the radial spread), and vice versa.

As mentioned earlier, temporal evolution in the classical picture corresponds to
uniform rigid rotation around the origin. The same is true in the quantum picture
when using a phase-space representation like the Wigner, *Q*, or
*P* functions, since temporal evolution corresponds to fractional
Fourier transformation in *q* of $\Psi (q;0)$, with scaling $\sigma =1$ and degree $\theta =\omega t$. Recall from Subsection 5.2
that fractional Fourier transformation corresponds to a rigid clockwise rotation for
these three representations.

Let us consider the Wigner function. Since Gaussian wave functions achieve the
maximum possible joint localization in *q* and *p*, the
Wigner representation of a quantum field that is closest to the phase-space
representation of a classical field (a single point in phase space) corresponds to a
Gaussian wave function, shifted to the desired phase-space point. This Wigner
function is itself Gaussian in both *q* and *p*. If
this Gaussian has circular cross-section, then the state is called a *coherent
state*, meaning that its width in *q* and
*p* does not change under temporal evolution. If the phase-space
Gaussian is instead elongated in one direction, then the state is referred to as a
*squeezed state*, and its marginal distributions will not only
oscillate harmonically but also “breathe,” expanding and contracting at a frequency $2\omega $ as time progresses. If the elongation is in the radial phase-space
direction [as shown in Fig. 23(c)], then the phase of the field (related to the angular spread
subtended by the distribution from the origin) is better determined than for a
coherent state of the same centroid, at the cost of the number of photons being more
uncertain.

By using a series of balanced homodyne measurements, one can estimate the marginal in
*p* of the Wigner function of the quantum state corresponding to
different stages of its evolution. These projections at different times are
equivalent to projections in different directions of the Wigner function at a fixed
time, from which the Wigner function is inferred by use of the inverse-Radon
transformation. An extensive review of this method is given in the book by
Leonhardt [39].

#### 10.2. Phase-Space Tomography of Classical Fields

Phase-space tomography has also been applied for the spatial or temporal
characterization of classical fields. In the temporal domain, it has been used for
characterizing trains of ultrashort pulses by estimating the chronocyclic Wigner
function [37,163]. Within the context of the spatial characterization of
stationary fields, the first proposition of using phase-space tomography to recover
the cross-spectral density of a paraxial optical field is due to Nugent [164], and the first experimental implementation
of this approach was by McAlister *et al.* [36]. This approach has been used for the characterization of
x-ray sources [165,166] whose main variation is in one transverse coordinate. Tu
and Tamura [167,168] realized that the same method can be cast in much simpler
terms if, instead of the Wigner function, one employs the ambiguity function. The key
to this approach is evident from the analogue of Eq. (140) for fields that depend only on one transverse variable,
i.e.,

*z*. Therefore, by sampling the intensity at many planes, both before and after the focal plane ($z=0$), one can estimate the ambiguity function at sufficient radial slices to allow accurate reconstruction through interpolation. The cross-spectral density is then calculated as

*et al.*[169,170].

Hazak [171] and Gori
*et al.* [172] realized
that, in the more general case where the transverse cross-spectral density depends on
both *x* and *y*, the information contained in the
field’s intensity over a volume (a function of three variables) is not sufficient for
recovering the initial cross-spectral density (a function of four variables). Within
the approach of Tu and Tamura, this problem can be appreciated from the fact that
Eq. (140) gives access to the
ambiguity function only within a three-dimensional “slice” in which ${\mathbf{x}}^{\prime}$ is parallel to ${\mathbf{p}}^{\prime}$. To address this problem, Raymer *et al.* [173] proposed introducing a fourth independent
parameter in the form of the varying focal distance of a cylindrical lens inserted
before the measurement of the intensity over a volume. The effect of such lenses is
to give access to marginal projections of the Wigner function that correspond to the
square modulus of different anamorphic fractional Fourier transforms. A system of
this type has been used for beams where the field is separable in the two transverse
coordinates [174]. Recently, Rojas
*et al.* [175] performed
simultaneous tomographic measurements over the phase spaces of position–direction and
time–frequency.

## 11. Generalized Wigner Functions: Achieving Exact Conservation

As discussed in the previous sections, the Wigner function allows a description of the propagation or evolution of optical wave fields in terms of a distribution that is conserved along the paths associated with simpler physical models. The condition for this description is that wave propagation in the system corresponds to a combination of linear canonical transformations and phase-space shifts. However, even in a situation as simple as free-space propagation beyond the paraxial regime, this conservation law is no longer satisfied exactly. In this section, we discuss modifications to the definition of the Wigner function for which the property of conservation along paths under propagation is preserved. The two cases that we concentrate on are nonparaxial field propagation and pulse propagation through general dispersive transparent media.

#### 11.1. General Procedure for Constructing Conserved Wigner Functions

We now give a simple prescription for the definition of generalizations of the Wigner function that present the property of describing propagation/evolution exactly as a rearrangement of the arguments [176,177], which is valid whenever the solutions to the wave equations in question can be expressed as superpositions of exponentials of imaginary linear combinations of the spatial or spatiotemporal variables. We concentrate on two situations:

*i. Nonparaxial propagation in free-space of a field
without evanescent components.* For the sake of illustration, let the
field be monochromatic, fully coherent, scalar, and limited to two dimensions, so
that it satisfies the Helmholtz equation:

*A*is the complex amplitude of the plane-wave component traveling in a direction at an angle

*θ*with respect to the

*z*axis.

*ii. One-dimensional propagation of pulses through
dispersive transparent media.* According to Eq. (166), pulses traveling through
dispersive media can be written as

Notice that, in both these situations, the field of interest is of the form

**g**is restricted to a curve (or manifold, for higher dimensions). For the case of nonparaxial propagation, this curve is the unit circle (or the unit sphere in three dimensions) of plane-wave directions, while for pulse propagation it is the dispersion curve

*k*versus

*ω*. Another physical situation, not discussed here, whose solutions can be written in the form in Eq. (191) is that of free-space relativistic quantum mechanics, where the manifold described by

**g**is a hyperboloid in the energy–momentum 4D space [177].

The goal is to define a generalized Wigner function ${\mathcal{W}}_{\mathcal{E}}(\mathbf{r},\xi )$ for a field of the form in Eq. (191) that is conserved under propagation/evolution and whose
marginal projection gives a desired bilinear property, e.g., the optical intensity.
Here, *ξ* is a suitable variable, e.g., the angle *θ*
with respect to the *z* axis of a *ray* passing through
**r** in the case of nonparaxial fields, or the uniform velocity
*v* of a *light particle* at position
*z* and time *t* in the case of dispersive
propagation. Mathematically, the conservation condition can be written as a
first-order differential equation of the form

In order to find ${\mathcal{W}}_{\mathcal{E}}$, the first step is to substitute Eq. (191) into the definition of the intensity (or the bilinear quantity of interest), leading to

**r**is within the exponent, it is easy to see that this definition of the generalized Wigner function satisfies the conservation requirement in Eq. (192) as long as the following relation holds:

*η*. Then, in order to evaluate ${\mathcal{W}}_{\mathcal{E}}$ at a given value of

*ξ*, the integration in

*α*must parametrize all pairs of points $\mathbf{g}\left({\eta}_{1}\right)$ and $\mathbf{g}\left({\eta}_{2}\right)$ that are joined by a line segment normal to $\mathbf{u}\left(\xi \right)$. That is, ${\eta}_{1}$ and ${\eta}_{2}$ must be chosen as corresponding to all pairs of points along the curve spanned by

**g**that are intersections with straight lines whose slopes are normal to $\mathbf{u}\left(\xi \right)$. In other words,

*ξ*regulates the slope of the line containing $\mathbf{g}\left({\eta}_{1}\right)$ and $\mathbf{g}\left({\eta}_{2}\right)$ while, for each

*ξ*, the displacement of this line is regulated by

*α*.

In the case of dispersion, the curve can present inflection points, in which a case
there might be more than two intersections with a straight line. If this is the case,
all possible pairs must be accounted for by the parametrization. On the other hand,
for higher-dimensional problems where **g** spans not a curve but a surface
or some higher-order manifold, the condition in Eq. (196) is not sufficient to find a unique change of variables,
i.e., a unique definition of the generalized Wigner function with the desired
properties. As discussed below, some other physically motivated criterion must be
imposed in this case to remove this ambiguity.

#### 11.2. Nonparaxial Field Propagation

For two-dimensional nonparaxial free propagation, $\mathbf{g}\left(\theta \right)=(sin\theta ,cos\theta )$ spans a unit circle. Note that, in this case $\mathbf{u}\left(\theta \right)=\mathbf{g}\left(\theta \right)$, so Eq. (196) takes the form

*θ*. It has been shown [45] that it also satisfies an inner-product overlap property analogous to that in Eq. (55).

Analogous generalized Wigner functions can be defined in the three-dimensional
case [182,183], where the manifold spanned by **g** is the surface of a
unit sphere. Let this sphere be parametrized directly by the unit vector
**u** [that is, $\mathbf{g}\left(\mathbf{u}\right)=\mathbf{u}$]. The change of variables is slightly more complicated in this
case, since it is necessary to include one more integration variable in the
definition of the Wigner function. The condition in Eq. (196) now becomes

*ϕ*is the extra integration variable needed for the change of variables. That is,

**u**must be constrained to the plane normal to ${\mathbf{u}}_{2}-{\mathbf{u}}_{1}$. For the sake of symmetry and in order to remove the ambiguity mentioned earlier, let us also choose

**u**to be coplanar with ${\mathbf{u}}_{1}$ and ${\mathbf{u}}_{2}$. This implies that

**u**must bisect ${\mathbf{u}}_{1}$ and ${\mathbf{u}}_{2}$, as shown in Fig. 25(b). As in the 2D case, we choose

*α*as the angle between ${\mathbf{u}}_{1}$ and ${\mathbf{u}}_{2}$. The resulting generalized Wigner function then has the form

**u**and to each other.

Wigner functions of this type have been used to describe the propagation not only of
the intensity of scalar fields but also of the electric and magnetic energy densities
and Poynting vector [184,185], polarization [185,186], and
cross-spectral density [187] of nonparaxial
electromagnetic fields. In the case of polarization, the Wigner functions are 3 × 3
matrices [analogous to the 2 × 2 matrix Wigner function in Eq. (155)] whose elements correspond to
correlations of different Cartesian components of the field. Further, upon
propagation across or reflection from interfaces between homogeneous media, these
generalized Wigner functions transform following approximately a radiometric model:
to leading order, they follow ray-optical laws of refraction and reflection and are
multiplied by a Fresnel transmissivity or reflectivity, with corrections taking the
form of derivatives of the leading term, which are less significant for more smoothly
varying, less coherent fields [188,189]. Generalized Wigner functions have also
been derived to model propagation through homogeneous anisotropic (e.g., birefringent
and/or chiral) transparent media [190]. For
these media, **g** spans not one but two manifolds, one associated with each
eigenpolarization of the medium. In general, these manifolds are ellipsoids, given
the dependence of the refractive index on the propagation direction.

The generalized Wigner functions in Eqs. (198) and (200) are valid
for fields that include components propagating in any direction. However, in the case
of nonparaxial fields composed exclusively of plane waves traveling towards larger
values of *z* (i.e., within the forward hemisphere of directions),
these Wigner functions can be written in a form that makes the description of
propagation formally equivalent to that for paraxial fields [191]. To see this, consider the generalized Wigner function for
two dimensions in Eq. (198). The
field is assumed to be composed of only forward-propagating plane waves; i.e., the
integral in Eq. (189) extends only
from $-\pi /2$ to $\pi /2$. Therefore, the Wigner function in Eq. (198) vanishes for $cos\theta <0$. Now perform a change in the directional parameter, so that instead
of the angle *θ*, we specify direction by the slope $\tau =tan\theta $. After also including a Jacobian factor $\left|\partial \right(\theta )/\partial (\tau \left)\right|$, we find the expression

*τ*. This is illustrated in Fig. 26 for a truncated cylindrical wave, corresponding to $A\left(\theta \right)$ being constant for $\left|\theta \right|\le \pi /4$ and zero otherwise. The alternative representation in Eq. (202) allows, among other things, the definition of a generalized ambiguity function for nonparaxial fields, given by the Fourier transform over

*x*and the inverse Fourier transform over

*τ*of ${\stackrel{\u0304}{\mathcal{W}}}_{U}(x,\tau ;z)$. Given the full analogy with the paraxial regime, the standard phase-space tomography techniques used in the paraxial 2D case can be extended to the nonparaxial regime [191].

The generalized Wigner functions in Eqs. (198) and (200) are
formally equivalent to wave-based definitions of the radiance proposed by Ovchinnikov
and Tatarskii [192] and Pedersen *et
al.* [193–195]. The derivation followed by those authors, however, was
quite different: first, the Wigner function over all three spatial variables is
calculated, yielding a distribution over a six-dimensional phase space. This
distribution turns out to be confined to the interior of a sphere in the Fourier
variables, and to diverge at the surface of this sphere. Then, a radial projection of
this distribution over the 3D Fourier variable is performed, with an appropriate
weight factor, to eliminate the redundant variable corresponding to the magnitude of
the direction vector. The resulting distribution presents the desired properties of
conservation along rectilinear rays, as well as a directional marginal related to the
field’s flux. However, to compute this distribution, one must first know the field in
all space. A very different derivation for what is essentially the same distribution
(up to an obliquity factor) was also proposed by Littlejohn and Winston [196]. Wolf *et al.* [178] also proposed a Fourier-subspace radial
projection for the case of two spatial dimensions, but by expressing the field as a
weighted superposition of plane waves, they found the expression in Eq. (198). This result was then extended to
fields in 3D space and with any state of coherence, leading to a family of
generalized Wigner functions [179,182,183]
whose marginals yield different local properties like the intensity, flux, and energy
density. Sheppard and Larkin [180,181] found a more direct derivation in the
two-dimensional case through a change of the angular variables to centroid and
difference, consistent with the construction discussed here. The generalized Wigner
functions in both two and three dimensions can also be calculated in terms of a
series of derivatives acting on the standard Wigner function [197].

#### 11.3. Pulse Propagation Through Transparent Dispersive Media

In the case of propagation through dispersive media, the generalized Wigner function
behaves as a statistical–mechanical distribution of independent “light particles”
traveling at different velocities *v*, since it is defined to obey the
propagation and marginal properties [176,177]

### 11.3a. Lorentz Model Dispersion

Consider first a nonabsorbing approximation to the Lorentz model of dispersion for a medium presenting a single resonance at ${\omega}_{\mathrm{R}}$, given by

*β*is a positive constant related to the plasma frequency of the medium. As shown in Ref. [176], the change of variables that results from imposing Eq. (207) can be written as

*k*given in Eq. (208). The corresponding Wigner function then takes the form

*v*, equal to the spatial intensity profile, is also shown.

### 11.3b. Waveguide Dispersion

As a second example, consider the case of a perfect metallic hollow waveguide, whose dispersion is given by

for $\omega <{\omega}_{\mathrm{co}}$, where ${\omega}_{\mathrm{co}}$ is the cutoff frequency. A convenient change of variables that satisfies Eq. (207) for this dispersion relation is given by

In the examples we just presented, we plotted the Wigner function over the phase
space of *z* versus *v* for different propagation
times *t*. However, there are other ways to present these Wigner
functions. For example, if the temporal intensity profiles are of more interest
than the spatial ones, it might be more useful to use a variant of these Wigner
functions in terms of *t* and a “slowness” variable $s=1/v$, as was done in Ref. [198]. This way, propagation in *z* corresponds to a
horizontal linear shearing in the phase space $(t,s)$. Alternatively, if it were preferable to express these Wigner
functions in the chronocyclic phase space, one can substitute $v\left(\omega \right)=1/{k}^{\prime}\left(\omega \right)$ and include a Jacobian factor for normalization, so that the
conserved chronocyclic Wigner function is given by

#### 11.4. Reduction to the Standard Wigner Function for Parabolic Manifolds

Notice that, for the simple case in which **g** spans a parabola with any
location and orientation, the generalized Wigner function reduces to a form of the
standard Wigner function. For a general parabolic curve, this vector can be
parametrized as

**a**,

**b**, and

**c**are three constant vectors. After making the change of variables ${\eta}_{1,2}(\xi ,\alpha )=\stackrel{\u0304}{\eta}\left(\xi \right)\mp \alpha /2$, which has a Jacobian equal to $\left|{\stackrel{\u0304}{\eta}}^{\prime}\right(\xi \left)\right|$, Eq. (196) can be written as

The reduction of the generalized Wigner function to the standard one for parabolic
curves or manifolds explains the usefulness of the standard Wigner function in the
case of paraxial free propagation and narrow-bandwidth pulses. In these cases, only a
small segment of the manifold spanned by **g** is of interest, which can
then be accurately approximated by a parabola: for paraxial optics we can approximate $\mathbf{g}=(sin\theta ,cos\theta )\approx (\theta ,1-{\theta}^{2}/2)$, while for narrow-band pulse propagation we can use either of the
quadratic approximations in Eq. (168) or (176). The same is
true in the case of quantum mechanics where, following the Klein–Gordon or Dirac
equations, the manifold for free evolution of a particle is a (two-sheeted)
hyperboloid of “radius” equal to the particle’s mass in the energy–momentum space. In
the nonrelativistic limit, only the region near the (positive energy) vertex of this
hyperboloid is of importance, so it can be approximated by a paraboloid. The standard
Wigner function is therefore useful in the nonrelativistic limit.

## 12. Concluding Remarks

An overview of the use of the Wigner function in optics was presented, with an emphasis on its application to model field propagation. As discussed briefly in Section 6, there are many other phase-space distributions that are commonly used in physics and engineering, each of them with a unique set of properties that make it well suited for some purposes. However, there are two properties that, together, confer the Wigner function a particular connection to physical intuition. The first is that the Wigner function for an optical field after spatial propagation or temporal evolution corresponds to a simple mapping of arguments of the Wigner function for the initial field. This property is valid in cases where propagation or evolution correspond mathematically to an arbitrary combination of linear canonical transformations and phase-space shifts. Such cases include propagation of paraxial stationary fields through optical systems that are linear in the ray domain, propagation of time-dependent fields through transparent media with linear dispersion, and the temporal evolution of quantum states of light. This mapping of the Wigner function can be interpreted as conservation along paths associated with a simpler physical model, like ray optics, free particle dynamics, or classical field evolution. The second key property is that, at any point of propagation over such systems, the marginal projection of the Wigner function over the Fourier variables is associated with an observable, physically meaningful local property of the field. Given these two properties, the Wigner function’s behavior mimics that of the phase-space distributions used in the simpler physical models for these situations (e.g., the basic radiance from the theory of radiometry). However, the Wigner function incorporates the fundamental limitations imposed by the space–bandwidth product theorem, and this causes it to be potentially negative within some regions of phase space, and to have strict restrictions in its functional form.

By modifying the definition of the Wigner function according to a simple geometrically
motivated prescription, the exact satisfaction of its two key properties can be extended
at least for certain physical situations that are not described by linear canonical
transformations and phase-space shifts. These situations include all cases where the
wave equation accepts a complete set of solutions corresponding to purely imaginary
exponentials of linear combinations of the variables. This is the case of nonparaxial
field propagation (if evanescent waves are excluded) and of pulse propagation through
homogeneous media with arbitrary dispersion (if absorption is neglected). This extension
brings to mind the following question: what is the most general class of system and
field for which a generalized Wigner function can be defined that satisfies exactly both
the marginal relation to the intensity and the property of conservation along paths?
This general question was not addressed in this article. Let us consider the case of
propagation in general transparent inhomogeneous linear media. There are inhomogeneous
refractive index distributions that allow a rigorous description of the propagation of
intensity as a conserved weight distribution for rays, even beyond the paraxial regime.
One trivial example of such case is that of inhomogeneous and anisotropic media that are
mappings of homogeneous media through what is known as transformation optics [202], since the Wigner function would then be the
result of the corresponding mapping of the Wigner function for a homogeneous medium. On
the other hand, there are other distributions for which this type of representation is
certainly impossible. Take the case of a gradient index waveguide that is symmetric with
respect to the *z* axis but that is anharmonic, i.e., that present
nonlinear mode dispersion. In the ray picture, propagation in *z* would
correspond to a swirling of phase space, since different rays oscillate around the axis
with different periods. In the wave domain, on the other hand, fields that are
concentrated spatially and directionally far from the axis of the waveguide are known to
return to their initial state after some propagation distance, due to the phenomenon
known as a wave packet revival [203]. Therefore,
there seems to be an irreconcilable difference in the behavior of the rays and the
waves.

In many cases, phase-space distributions like the Wigner function are valuable more for conceptual reasons than for computational ones. They provide an intuitive framework that leads to useful ideas and interpretations, even if the methods resulting from them can then be expressed in a more direct way that bypasses the phase-space picture. Computationally, the increased number of variables and complicated functional form of the Wigner function or its generalizations discussed here can be important obstacles towards numerical implementation, particularly in the case of coherent fields or pure states. For partially coherent fields or mixed states, on the other hand, the value of Wigner functions can be not only conceptual but also computational, because using standard propagation formulas can become computationally prohibitive in this case. If used properly, the Wigner function can allow significant computational savings in these cases, as discussed in Ref. [120].

This article is dedicated to the memory of Luis Edgar Vicent (1972–2011), friend, colleague and coauthor of Ref. [120].

## References and Notes

**1. **N. L. Balazs and B. K. Jennings, “Wigner’s function and other distribution functions
in mock phase spaces,” Phys. Rep. **104**, 347‒ 391 (1984).
[CrossRef]

**2. **E. P. Wigner, “On the quantum correction for thermodynamic
equilibrium,” Phys. Rev. **40**, 749‒ 759 (1932).
[CrossRef]

**3. **J. Ville, “Thèorie et applications de la notion de signal
analytique,” Cables Transm. **2A**, 6174 (1948).

**4. **L. S. Dolin, “Beam description of weakly inhomogeneous wave
fields,” Izv. Vyssh. Uchebn. Zaved. Radiofiz. **7**, 559‒ 563 (1964).

**5. **A. Walther, “Radiometry and coherence,”
J. Opt. Soc. Am. **58**, 1256‒ 1259 (1968).
[CrossRef]

**6. **R. W. Boyd, *Radiometry and the Detection of Optical Radiation*
(Wiley, 1983), pp.
13‒ 27).

**7. **P. Moon and D. E. Spencer, *The Photic Field* (MIT
Press, 1981).

**8. **L. Cohen, *Time–Frequency Analysis* (Prentice
Hall, 1995).

**9. **W. Mecklenbräuker and F. Hlawatsch, *The Wigner Distribution: Theory and Applications in Signal
Processing* (Elsevier,
1997).

**10. **H. W. Lee, “Theory and application of the quantum phase-space
distribution functions,” Phys. Rep. **259**, 147‒ 211 (1995).
[CrossRef]

**11. **G. S. Agarwal and E. Wolf, “Calculus for functions of noncommuting operators and
general phase-space methods in quantum mechanics. I. Mapping theorems and ordering
of functions of noncommuting operators,” Phys. Rev.
D **2**, 2161‒ 2186 (1970).
[CrossRef]

**12. **G. S. Agarwal and E. Wolf, “Calculus for functions of noncommuting operators and
general phase-space methods in quantum mechanics. II. Quantum mechanics in phase
space,” Phys. Rev. D **2**, 2187‒ 2205 (1970).
[CrossRef]

**13. **G. S. Agarwal and E. Wolf, “Calculus for functions of noncommuting operators and
general phase-space methods in quantum mechanics. III. A generalized Wick theorem
and multitime mapping,” Phys. Rev. D **2**, 2206‒ 2225 (1970).
[CrossRef]

**14. **M. Hillery, R. F. O’Connell, M. O. Scully, and E. P. Wigner, “Distribution functions in physics:
Fundamentals,” 106 121‒ 167 (1984).

**15. **M. E. Testorf, B. M. Hennelly, and J. Ojeda-Castañeda, *Phase-Space Optics: Fundamentals and Applications*
(McGraw-Hill,
2009).

**16. **A. Torre, *Linear Ray and Wave Optics in Phase Space: Bridging Ray and
Wave Optics via the Wigner Phase-Space Picture*
(Elsevier, 2005).

**17. **D. Dragoman, “The Wigner distribution function in optics and
optoelectronics,” in *Progress in Optics XXXVII*,
E. Wolf (Elsevier, 1997), pp.
1‒ 56.

**18. **J. T. Sheridan, W. T. Rhodes, and B. M. Hennelly, “Wigner Optics,” Proc.
SPIE **5827**, 627‒ 638 (2005).

**19. **A. Papoulis, *The Fourier Integral and its Applications*
(McGraw-Hill,
1962).

**20. **G. Folland and A. Sitaram, “The uncertainty principle: a mathematical
survey,” J. Fourier Anal. Appl. **3**, 207‒ 238 (1997).
[CrossRef]

**21. **E. U. Condon, “Immersion of the Fourier transform in a continuous
group of functional transformations,” Proc. Natl. Acad.
Sci. **23**, 158‒ 163 (1937).
[CrossRef]

**22. **V. Namias, “The fractional order Fourier transform and its
application to quantum mechanics,” IMA J. Appl.
Math. **25**, 241‒ 265 (1980).
[CrossRef]

**23. **A. W. Lohmann, D. Mendlovic, and Z. Zalevsky, E. Wolf, “Fractional transformations in
optics,” in *Progress in Optics XXXVIII*1998), pp. 263‒ 242.

**24. **H. M. Ozaktas, Z. Zalevsky, and M. A. Kutay, *The Fractional Fourier Transform with Applications in Optics
and Signal Processing* (John Wiley &
Sons, 2001).

**25. **A. Erdélyi, *Asymptotic Expansions*
(Dover, 1956).

**26. **D. Mendlovic, Y. Bitran, R. G. Dorsch, C. Ferreira, J. Garcia, and H. M. Ozaktaz, “Anamorphic fractional Fourier transform: optical
implementation and applications,” Appl. Opt. **34**, 7451‒ 7456 (1995).
[CrossRef] [PubMed]

**27. **S. A. Collins, “Lens-system diffraction integral written in terms of
matrix optics,” J. Opt. Soc. Am. **60**, 1168‒ 1177 (1970).
[CrossRef]

**28. **M. Moshinsky and C. Quesne, “Linear canonical transformations and their unitary
representations,” J. Math. Phys. **12**(8), 1772‒ 1783
(1971). [CrossRef]

**29. **K. B. Wolf, *Integral Transforms in Science and Engineering*
(Plenum Press, 1979), Ch.
9, 10.

**30. **T. Alieva and M. J. Bastiaans, “Properties of the linear canonical integral
transformation,” J. Opt. Soc. Am. A **24**, 3658‒ 3665 (2007).
[CrossRef]

**31. **D. Mustard, “The fractional Fourier transform and the Wigner
distribution,” J. Aust. Math. Soc. B-Appl. Math. **38**, 209‒ 219 (1996)
Published earlier as Applied Mathematics Preprint AM89/6 School of
Mathematics, UNSW, Sydney, Australia (1989). [CrossRef]

**32. **A. W. Lohmann, “Image rotation, Wigner rotation, and the fractional
Fourier transform,” J. Opt. Soc. Am. A **10**, 2181‒ 2186 (1993).
[CrossRef]

**33. **A. W. Lohmann and B. H. Soffer, “Relationships between the Radon–Wigner and
fractional Fourier transforms,” J. Opt. Soc. Am.
A **11**, 1798‒ 1801 (1994).
[CrossRef]

**34. **See Chapter 4 by G. Saavedra and W. Furlan in Ref. [15], pp. 107–164.

**35. **J. Bertrand and P. Bertrand, “A tomographic approach to Wigner’s
function,” Found. Phys. **17**, 397‒ 405 (1987).
[CrossRef]

**36. **D. F. McAlister, M. Beck, L. Clarke, A. Mayer, and M. G. Raymer, “Optical phase retrieval by phase space tomography
and fractional-order Fourier transforms,” Opt.
Lett. **20**, 1181‒ 1183 (1995).
[CrossRef] [PubMed]

**37. **M. Beck, M. G. Raymer, I. A. Walmsley, and V. Wong, “Chronocyclic tomography for measuring the amplitude
and phase structure of optical pulses,” Opt.
Lett. **18**, 2041‒ 2043 (1993).
[CrossRef] [PubMed]

**38. **D. T. Smithey, M. Beck, M. G. Raymer, and A. Faridani, “Measurement of the Wigner distribution and the
density matrix of a light mode using optical homodyne tomography: application to
squeezed states and the vacuum,” Phys. Rev. Lett. **70**, 1244‒ 1247 (1993).
[CrossRef] [PubMed]

**39. **U. Leonhardt, *Measuring the Quantum State of Light*
(Cambridge U. Press,
1997).

**40. **C. Cohen-Tannoudji, B. Diu, and F. Laloe, *Quantum Mechanics*, Vol. 1
(Wiley, 1977), pp.
214‒ 227.

**41. **J. E. Moyal, “Quantum mechanics as a statistical
theory,” Math. Proc. Camb. Phil. Soc. **45**, 99‒ 124 (1949).
[CrossRef]

**42. **R. G. Littlejohn and R. Winston, “Generalized radiance and
measurement,” J. Opt. Soc. Am. A **12**, 2736‒ 2743 (1995).
[CrossRef]

**43. **D. Dragoman, “Wigner-distribution-function representation of the
coupling coefficient,” Appl. Opt. **34**, 6758‒ 6763 (1995).
[CrossRef] [PubMed]

**44. **A. Wax and J. E. Thomas, “Optical heterodyne imaging and Wigner phase space
distributions,” Opt. Lett. **21**, 1427‒ 1429 (1996).
[CrossRef] [PubMed]

**45. **M. A. Alonso, “Measurement of Helmholtz wave
fields,” J. Opt. Soc. Am. A **17**, 1256‒ 1264 (2000).
[CrossRef]

**46. **M. V. Berry, “Semi-classical mechanics in phase space: a study of
Wigners function,” Philos. Trans. R. Soc. Lond. **287**, 237‒ 271 (1977).
[CrossRef]

**47. **M. V. Berry and N. L. Balazs, “Evolution of semiclassical quantum states in phase
space,” J. Phys. A Math. Phys. **12**, 625‒ 642 (1979).
[CrossRef]

**48. **M. V. Berry, “Quantum scars of classical closed orbits in phase
space,” Proc. R. Soc. Lond. Ser. A **423**, 219‒ 231 (1989).
[CrossRef]

**49. **M. A. Alonso and G. W. Forbes, “Phase-space distributions for high-frequency
fields,” J. Opt. Soc. Am. A **17**, 2288‒ 2300 (2000).
[CrossRef]

**50. **L. Mandel and E. Wolf, *Optical Coherence and Quantum Optics*
(Cambridge University Press, 1995),
Sec. 4.7.

**51. **M. J. Bastiaans, “Uncertainty principle for partially coherent
light,” J. Opt. Soc. Am. **73**, 251‒ 255 (1983).
[CrossRef]

**52. **A. Starikov, “Effective number of degrees of freedom of partially
coherent sources,” J. Opt. Soc. Am. **72**, 1538‒ 1544 (1982).
[CrossRef]

**53. **See Ref. [50], p.
261.

**54. **L. Cohen, “Time–frequency distributions—a
review,” Proc. IEEE **77**, 941‒ 981 (1989).
[CrossRef]

**55. **J. H. Eberly and K. Wódkiewicz, “The time-dependent physical spectrum of
light,” J. Opt. Soc. Am. **67**, 1252‒ 1261 (1977).
[CrossRef]

**56. **K. Husimi, “Some formal properties of the density
matrix,” Proc. Phys. Math. Soc. Jpn. **22**, 264‒ 314 (1940).

**57. **Y. Kano, “A new phase-space distribution function in the
statistical theory of the electromagnetic field,” J. Math.
Phys. **6**, 1913‒ 1915 (1965).
[CrossRef]

**58. **C. L. Mehta and E. C. G. Sudarshan, “Relation between Quantum and Semiclassical
Description of Optical Coherence,” Phys. Rev. **138**, B274‒ B280 (1965).
[CrossRef]

**59. **R. J. Glauber, “Optical coherence and photon
statistics,” in *Quantum Optics and Electronics*,
C. Dewitt, A. Blandin, and C. Cohen-Tannoudji (Gordon and Breach,
1965), p. 65).

**60. **R. J. Glauber, “Coherent and incoherent states of the radiation
field,” Phys. Rev. **131**, 2766‒ 2788 (1963).
[CrossRef]

**61. **J. R. Klauder, “Continuous representation theory. I. Postulates of
continuous representation theory,” J. Math. Phys. **4**, 1055‒ 1058 (1963).
[CrossRef]

**62. **E. C. G. Sudarshan, “Equivalence of semiclassical and quantum mechanical
descriptions of statistical light beams,” Phys. Rev.
Lett. **10**, 277‒ 279 (1963).
[CrossRef]

**63. **J. G. Kirkwood, “Quantum statistics of almost classical
ensembles,” Phys. Rev. **44**, 31‒ 37 (1933).
[CrossRef]

**64. **A. W. Rihaczek, “Signal energy distribution in time and
frequency,” IEEE Trans. Info. Theory **14**, 369‒ 374 (1968).
[CrossRef]

**65. **H. Margenau and R. N. Hill, “Correlation between measurements in quantum
theory,” Prog. Theoret. Phys. **26**, 722‒ 738 (1961).
[CrossRef]

**66. **P. M. Woodward, *Probability and Information Theory with Applications to
Radar* (Pergamon,
1953).

**67. **A. Papoulis, “Ambiguity function in Fourier
optics,” J. Opt. Soc. Am. **64**, 779‒ 788 (1974).
[CrossRef]

**68. **L.-P. Guigay, *Ambiguity Function in Optical Imaging*,
Chapter 2 of Ref. [15], pp.
45–62.

**69. **L. Cohen, “Generalized phase-space
distributions,” J. Math. Phys. **7**, 781‒ 786 (1966).
[CrossRef]

**70. **M. J. Bastiaans, *Wigner Distribution in Optics*, Chapter 1
of Ref. [15], pp.
1–44.

**71. **R. K. Luneburg, *Mathematical Theory of Optics*
(University of California Press,
1966), pp. 246‒ 257.

**72. **See Ref. [71], pp.
103–110.

**73. **M. Born and E. Wolf, *Principles of Optics*, 7th Ed.
(Cambridge University Press, 1999),
pp. 142‒ 144.

**74. **H. A. Buchdahl, *Hamiltonian Optics*
(Dover, 1993), pp. 7‒
12.

**75. **See Ref. [6], pp.
13–25.

**76. **See Ref. [6], pp.
75–80.

**77. **A. Walther, “Lenses, wave optics and eikonal
functions,” J. Opt. Soc. Am. **59**, 1325‒ 1333 (1969).
[CrossRef]

**78. **A. Walther, *The Ray and Wave Theory of Lenses*
(Cambridge University Press, 1995),
pp. 169‒ 187.

**79. **M. A. Alonso and G. W. Forbes, “Semigeometrical estimation of Green’s functions and
wave propagators in optics,” J. Opt. Soc. Am. A **14**, 1076‒ 1086 (1997).
[CrossRef]

**80. **J. H. Van Vleck, “The correspondence principle in the statistical
interpretation of quantum mechanics,” Proc. Natl. Acad.
Sci. USA **14**, 178‒ 188 (1928).
[CrossRef]

**81. **M. C. Gutzwiller, “Periodic orbits and classical quantization
conditions,” J. Math. Phys. **12**, 343‒ 358 (1971).
[CrossRef]

**82. **D. Mendlovic and H. M. Ozaktas, “Fractional Fourier transforms and their optical
implementation: I,” J. Opt. Soc. Am. A **10**, 1875‒ 1881 (1993).
[CrossRef]

**83. **M. J. Bastiaans, “Wigner distribution function applied to optical
signals and systems,” Opt. Commun. **25**, 26‒ 30 (1978).
[CrossRef]

**84. **M. J. Bastiaans, “Transport equations for the Wigner distribution
function,” J. Mod. Opt. **26**, 1265‒ 1272 (1979).

**85. **M. J. Bastiaans, “The Wigner distribution function and Hamilton’s
characteristics of a geometric–optical system,” Opt.
Commun. **30**, 321‒ 326 (1979).
[CrossRef]

**86. **M. J. Bastiaans, “Wigner distribution function and its application to
first-order optics,” J. Opt. Soc. Am. **69**, 1710‒ 1716 (1979).
[CrossRef]

**87. **J. W. Goodman, *Introduction to Fourier Optics*
(McGraw-Hill, 1988), pp.
101‒ 136.

**88. **K.-H. Brenner, A. W. Lohmann, and J. Ojeda-Castañeda, “The ambiguity function as a polar display of the
OTF,” Opt. Commun. **44**, 323‒ 326 (1983).
[CrossRef]

**89. **S. B. Oh and G. Barbastathis, “Axial imaging necessitates loss of lateral shift
invariance: proof with the Wigner analysis,” Appl.
Opt. **48**, 5881‒ 5888 (2009).
[CrossRef] [PubMed]

**90. **S. B. Mehta and C. J. R. Sheppard, “Using the phase-space imager to analyze partially
coherent imaging systems: bright-field, phase contrast, differential interference
contrast, differential phase contrast, and spiral phase contrast,”
J. Mod. Opt. **57**, 718‒ 739 (2010).
[CrossRef]

**91. **J. E. Durnin, J. J. Miceli, and J. H. Eberly, “Diffraction-free beams,”
Phys. Rev. Lett. **58**, 1499‒ 1501 (1987).
[CrossRef] [PubMed]

**92. **W. P. Schleich, J. P. Dahl, and S. Varró, “Wigner function for a free particle in two
dimensions: a tale of interference,” Opt. Commun. **283**, 786‒ 789 (2010).
[CrossRef]

**93. **M. V. Berry and N. L. Balazs, “Nonspreading wave packets,”
Am. J. Phys. **47**, 264‒ 267 (1979).
[CrossRef]

**94. **M. V. Berry and F. J. Wright, “Phase-space projection identities for diffraction
catastrophes,” J. Phys. A. Math. Gen. **13**, 149‒ 160 (1980).
[CrossRef]

**95. **R.-P. Chen, H.-P. Zheng, and C.-Q. Dai, “Wigner distribution function of an Airy
beam,” J. Opt. Soc. Am. A **28**, 1307‒ 1311 (2011).
[CrossRef]

**96. **G. Siviloglou and D. Christodoulides, “Accelerating finite energy Airy
beams,” Opt. Lett. **32**, 979‒ 981 (2007).
[CrossRef] [PubMed]

**97. **G. Siviloglou, J. Broky, A. Dogariu, and D. N. Christodoulides, “Observation of accelerating Airy
beams,” Phys. Rev. Lett. **99**, 213901-1-213901-4 (2007). [CrossRef]

**98. **E. R. Dowski and W. T. Cathey, “Extended depth of field through wave-front
coding,” Appl. Opt. **34**, 1859‒ 1866 (1995).
[CrossRef] [PubMed]

**99. **W. T. Welford, “Use of annular apertures to increase focal
depth,” J. Opt. Soc. Am. **50**, 749 (1960). [CrossRef]

**100. **C. J. R. Sheppard, D. K. Hamilton, and I. J. Cox, “Optical microscopy with extended depth of
field,” Proc. R. Soc. Lond. A **387**, 171‒ 186 (1983).
[CrossRef]

**101. **J. Ojeda-Castañeda, L. R. Berriel-Valdos, and E. Montes, “Bessel annular apodizers: imaging
characteristics,” Appl. Opt. **26**, 2770‒ 2772 (1987).
[CrossRef] [PubMed]

**102. **J. Ojeda-Castañeda, L. R. Berriel-Valdos, and E. Montes, “Ambiguity function as a design tool for high focal
depth,” Appl. Opt. **27**, 790‒ 795 (1988).
[CrossRef] [PubMed]

**103. **Q. Yang, L. Liu, J. Sun, Y. Zhu, and W. Lu, “Analysis of optical systems with extended depth of
field using the Wigner distribution function,” Appl.
Opt. **45**, 8586‒ 8595 (2006).
[CrossRef] [PubMed]

**104. **M. Testorf and J. Ojeda-Castañeda, “Fractional Talbot effect: analysis in phase
space,” J. Opt. Soc. Am. A **13**, 119‒ 125 (1996).
[CrossRef]

**105. **M. Testorf, “Designing Talbot array illuminators with phase-space
optics,” J. Opt. Soc. Am. A **23**, 187‒ 192 (2006).
[CrossRef]

**106. **M. E. Testorf, *Self-imaging in Phase Space*Chapter 9 of Ref. [15], pp.
279–307.

**107. **J. Lancis, E. E. Sicre, A. Pons, and G. Saavedra, “Achromatic white-light self-imaging phenomenon: an
approach using the Wigner distribution function,” J. Mod.
Opt. **42**, 425‒ 434 (1995).
[CrossRef]

**108. **A. Stern and B. Javidi, “Sampling in the light of Wigner
distribution,” J. Opt. Soc. Am. A **21**, 360‒ 366 (2004)
errata, **21**, 1602–1612 (2004). [CrossRef]

**109. **B. M. Hennelly, J. J. Healy, and J. T. Sheridan, *Sampling and Phase Space*, Chapter 10 of
Ref. [15], pp.
309–336.

**110. **F. Gori, “Fresnel transform and sampling
theorem,” Opt. Commun. **39**, 293‒ 297 (1981).
[CrossRef]

**111. **A. Stern and B. Javidi, “Improved-resolution digital holography using the
generalized sampling theorem for locally band-limited fields,”
J. Opt. Soc. Am. A **23**, 1227‒ 1235 (2006).
[CrossRef]

**112. **A. Stern and B. Javidi, “Space-bandwidth conditions for efficient
phase-shifting digital holographic microscopy,” J. Opt.
Soc. Am. A **25**, 736‒ 741 (2008).
[CrossRef]

**113. **K. B. Wolf and A. L. Rivera, “Holographic information in the Wigner
function,” Opt. Commun. **144**, 36‒ 42 (1997).
[CrossRef]

**114. **A. W. Lohmann, M. E. Testorf, and J. Ojeda-Castañeda, “Holography and the Wigner function,”
in *The Art and Science of Holography: A Tribute to Emmett Leith and Yuri
Denisyuk*, H. J. Caulfield (SPIE Press, 2004), pp.
127‒ 144.

**115. **A. W. Lohmann, M. E. Testorf, J. Ojeda-Castañeda, and A. W. Lohmann, “The space–bandwidth product, applied to spatial
filtering and to holography,” in *Selected Papers on
Phase-Space Optics* (SPIE Press,
2006), pp. 11‒ 32.

**116. **M. Testorf and A. W. Lohmann, “Holography in phase space,”
Appl. Opt. **47**, A70‒ A77 (2008).
[CrossRef] [PubMed]

**117. **S. B. Oh and G. Barbastathis, “Wigner distribution function of volume
holograms,” Opt. Lett. **34**, 2584‒ 2586 (2009).
[CrossRef] [PubMed]

**118. **M. Testorf, “Analysis of the moiré effect by use of the Wigner
distribution function,” J. Opt. Soc. Am. A **17**, 2536‒ 2542 (2000).
[CrossRef]

**119. **N. Morelle, M. E. Testorf, N. Thirion, and M. Saillard, “Electromagnetic probing for target detection:
rejection of surface clutter based on the Wigner distribution,”
J. Opt. Soc. Am. A **26**, 1178‒ 1186 (2009).
[CrossRef]

**120. **L. E. Vicent and M. A. Alonso, “Generalized radiometry as a tool for the propagation
of partially coherent fields,” Opt. Commun. **207**, 101‒ 112 (2002).
[CrossRef]

**121. **E. Wolf, “Radiometric model for propagation of
coherence,” Opt. Lett. **19**, 2024‒ 2026 (1994).
[CrossRef] [PubMed]

**122. **A. T. Friberg and S. Yu. Popov, “Radiometric description of intensity and coherence
in generalized holographic axicon images,” Appl.
Opt. **35**, 3039‒ 3046 (1996).
[CrossRef] [PubMed]

**123. **K. Duan and B. Lü, “Wigner-distribution-function matrix and its
application to partially coherent vectorial nonparaxial beams,”
J. Opt. Soc. Am. B **22**, 1585‒ 1593 (2005).
[CrossRef]

**124. **R. Castañeda and J. Garcia-Sucerquia, “Electromagnetic spatial coherence
wavelets,” J. Opt. Soc. Am. A **23**, 81‒ 90 (2006).
[CrossRef]

**125. **A. Luis, “Spatial-angular Mueller matrices,”
Opt. Commun. **263**, 141‒ 146 (2006).
[CrossRef]

**126. **A. Luis, “Ray picture of polarization and coherence in a Young
interferometer,” J. Opt. Soc. Am. A **23**, 2855‒ 2860 (2006).
[CrossRef]

**127. **R. Castañeda, J. Carrasquilla, and J. Herrera, “Radiometric analysis of diffraction of
quasi-homogeneous optical fields,” Opt. Commun. **273**, 8‒ 20 (2007).
[CrossRef]

**128. **R. Castañeda and J. Carrasquilla, “Spatial coherence wavelets and phase-space
representation of diffraction,” Appl. Opt. **47**, E76‒ E87 (2008).
[CrossRef] [PubMed]

**129. **M. A. Alonso, “Diffraction of paraxial partially coherent fields by
planar obstacles in the Wigner representation,” J. Opt.
Soc. Am. A **26**, 1588‒ 1597 (2009).
[CrossRef]

**130. **A. Wax and J. E. Thomas, “Measurement of smoothed Wigner phase-space
distributions for small-angle scattering in a turbid medium,”
J. Opt. Soc. Am. A **15**, 1896‒ 1908 (1998).
[CrossRef]

**131. **H. T. Yura, L. Thrane, and P. E. Andersen, “Closed-form solution for the Wigner phase-space
distribution function for diffuse reflection and small-angle scattering in a
random medium,” J. Opt. Soc. Am. A **17**, 2464‒ 2474 (2000).
[CrossRef]

**132. **C.-C. Cheng and M. G. Raymer, “Propagation of transverse optical coherence in
random multiple-scattering media,” Phys. Rev. A **62**, 023811 (2000). [CrossRef]

**133. **A. C. Fannjiang, “White-noise and geometrical optics limits of
Wigner–Moyal equation for wave beams in turbulent media,”
Commun. Math. Phys. **254**, 289‒ 322 (2005).
[CrossRef]

**134. **D. Dragoman, “Wigner distribution function in nonlinear
optics,” Appl. Opt. **35**, 4142‒ 4146 (1996).
[CrossRef] [PubMed]

**135. **D. Dragoman, J. P. Meunier, and M. Dragoman, “Beam-propagation method based on the Wigner
transform: a new formulation,” Opt. Lett. **22**, 1050‒ 1052 (1997).
[CrossRef] [PubMed]

**136. **B. Hall, M. Lisak, D. Anderson, R. Fedele, and V. E. Semenov, “Statistical theory for incoherent light propagation
in nonlinear media,” Phys. Rev. E **65**, 035602R (2002). [CrossRef]

**137. **M. Lisak, L. Helczynski, and D. Anderson, “Relation between different formalisms describing
partially incoherent wave propagation in nonlinear optical media,”
Opt. Commun **220**, 321‒ 323 (2003).
[CrossRef]

**138. **H. Gao, L. Tian, B. Zhang, and G. Barbastathis, “Iterative nonlinear beam propagation using
Hamiltonian ray tracing and Wigner distribution function,”
Opt. Lett. **35**, 4148‒ 4150 (2010).
[CrossRef] [PubMed]

**139. **S. Abe and J. T. Sheridan, “Wigner optics in the metaxial
regime,” Optik **114**, 139‒ 141 (). [CrossRef]

**140. **Yu. A. Kravtsov and L. A. Apresyan, “Radiative transfer: new aspects of the old
theory,” in *Progress in Optics*Vol.
XXXVI, E. Wolf (North Holland, 1996),
pp. 179‒ 244.

**141. **L. A. Apresyan and Yu. A. Kravtsov, *Radiation Transfer: Statistical and Wave Aspects*
(Gordon and Breach,
1996).

**142. **A. T. Friberg, *Selected Papers on Coherence and Radiometry*,
Milestone Series Vol. MS69 (SPIE
Optical Engineering Press, 1993).

**143. **E. Wolf, “Coherence and radiometry,”
J. Opt. Soc. Am. **68**, 6‒ 17 (1978). [CrossRef]

**144. **A. T. Friberg, “Effects of coherence in radiometry,”
Opt. Eng. **21**, 927‒ 936 (1982).

**145. **G. S. Agarwal, J. T. Foley, and E. Wolf, “The radiance and phase-space representations of the
cross-spectral density operator,” Opt. Commun. **62**, 67‒ 72 (1987).
[CrossRef]

**146. **A. Walther, “Radiometry and coherence,”
J. Opt. Soc. Am. **63**, 1622‒ 1623 (1973).
[CrossRef]

**147. **A. T. Friberg, “On the generalized radiance associated with
radiation from a quasihomogeneous planar source,” Opt.
Acta **28**, 261‒ 277 (1981).
[CrossRef]

**148. **J. T. Foley and E. Wolf, “Radiometry as a short-wavelength limit of
statistical wave theory with globally incoherent sources,”
Opt. Commun. **55**, 236‒ 241 (1985).
[CrossRef]

**149. **J. Azaña, “Time–frequency (Wigner) analysis of linear and
nonlinear pulse propagation in optical fibers,” EURASIP J.
Appl. Signal Process. **2005**, 1554‒ 1565 (2005).
[CrossRef]

**150. **P. Loughlin and L. Cohen, “A Wigner approximation method for wave
propagation,” J. Acoust. Soc. Am. **118**, 1268‒ 1271 (2005).
[CrossRef]

**151. **J. Ojeda-Castañeda, J. Lancis, C. M. Gómez-Sarabia, V. Torres-Company, and P. Andrés, “Ambiguity function analysis of pulse train
propagation: applications to temporal Lau filtering,” J.
Opt. Soc. Am. A **24**, 2268‒ 2273 (2007).
[CrossRef]

**152. **P. Loughlin and L. Cohen, “Approximate wave function from approximate
non-representable Wigner distributions,” J. Mod.
Opt. **55**, 3379‒ 3387 (2008).
[CrossRef]

**153. **L. Cohen, P. Loughlin, and G. Okopal, “Exact and approximate moments of a propagating
pulse,” J. Mod. Opt. **55**, 3349‒ 3358 (2008).
[CrossRef]

**154. **C. Dorrer and I. A. Walmsley, *Phase space in ultrafast optics*, Chapter
11 of Ref. [15], pp.
337–383.

**155. **I. A. Walmsley and C. Dorrer, “Characterization of ultrashort electromagnetic
pulses,” Adv. Opt. Phot. **1**, 308‒ 437 (2009).
[CrossRef]

**156. **B. H. Kolner, “Space–time duality and the theory of temporal
imaging,” IEEE J. Quantum Electron. **30**, 1951‒ 1963 (1994).
[CrossRef]

**157. **B. E. A. Saleh and M. C. Teich, *Fundamentals of Photonics*, 2nd
ed. (Wiley, 2007), p.
188.

**158. **H. Lajunen, J. Tervo, J. Turunen, P. Vahimaa, and F. Wyrowski, “Spectral coherence properties of temporally
modulated stationary light sources,” Opt. Express **11**, 1894‒ 1899 (2003).
[CrossRef] [PubMed]

**159. **B. J. Davis, “Observable coherence theory for statistically
periodic fields,” Phys. Rev. A **76**, 043843 (2007). [CrossRef]

**160. **K. Vogel and H. Risken, “Determination of quasiprobability distributions in
terms of probability distributions for the rotated quadrature
phase,” Phys. Rev. A **40**, 2847‒ 2849 (1989).
[CrossRef] [PubMed]

**161. **C. C. Gerry and P. L. Knight, *Introductory Quantum Optics*
(Cambridge University Press, 2005),
pp. 56‒ 71.

**162. **W. P. Schleich, *Quantum Optics in Phase Space*
(Wiley-VCH, 2001).

**163. **C. Dorrer and I. Kang, “Complete temporal characterization of short optical
pulses by simplified chronocyclic tomography,” Opt.
Lett. **28**, 1481‒ 1483 (2003).
[CrossRef] [PubMed]

**164. **K. A. Nugent, “Wave field determination using three-dimensional
intensity information,” Phys. Rev. Lett. **68**, 2261‒ 2264 (1992).
[CrossRef] [PubMed]

**165. **C. Q. Tran, A. G. Peele, A. Roberts, K. A. Nugent, D. Paterson, and I. McNulty, “X-ray imaging: a generalized approach using
phase-space tomography,” J. Opt. Soc. Am. A **22**, 1691‒ 1700 (2005).
[CrossRef]

**166. **C. Q. Tran, A. P. Mancuso, B. B. Dhal, K. A. Nugent, A. G. Peele, Z. Cai, and D. Paterson, “phase space reconstruction of focused x-ray
fields,” J. Opt. Soc. Am.
A **23**, 1779‒ 1786 (2006).
[CrossRef]

**167. **J. Tu and S. Tamura, “Wave field determination using tomography of the
ambiguity function,” Phys. Rev. E **55**, 1946‒ 1949 (1997).
[CrossRef]

**168. **J. Tu and S. Tamura, “Analytic relation for recovering the mutual
intensity by means of intensity information,” J. Opt. Soc.
Am. A **15**, 202‒ 206 (1998).
[CrossRef]

**169. **R. Horstmeyer, S. B. Oh, and R. Raskar, “Iterative aperture mask design in phase space using
a rank constraint,” Opt. Express **18**, 22545‒ 22555 (2010).
[CrossRef] [PubMed]

**170. **“Partially coherent ambiguity
functions for depth-variant point spread function design,”
presentation during PIERS, Marrakesh, March,
2011.

**171. **G. Hazak, “Comment on ‘Wave field determination using
three-dimensional intensity information’,” Phys. Rev.
Lett. **69**, 2874‒ 2874 (1992).
[CrossRef] [PubMed]

**172. **F. Gori, M. Santarsiero, and G. Guattari, “Coherence and the spatial distribution of
intensity,” J. Opt. Soc. Am. A **10**, 673‒ 679 (1993).
[CrossRef]

**173. **M. G. Raymer, M. Beck, and D. F. McAlister, “Complex wave-field reconstruction using phase-space
tomography,” Phys. Rev. Lett. **72**, 1137‒ 1140 (1994).
[CrossRef] [PubMed]

**174. **A. Cámara, T. Alieva, J. A. Rodrigo, and M. L. Calvo, “Phase space tomography reconstruction of the Wigner
distribution for optical beams separable in Cartesian
coordinates,” J. Opt. Soc. Am. A **26**, 1301‒ 1306 (2009).
[CrossRef]

**175. **P. Rojas, R. Blaser, Y. M. Sua, and K. F. Lee, “Optical phase-space-time-frequency
tomography,” Opt. Express **19**, 7480‒ 7490 (2011).
[CrossRef] [PubMed]

**176. **J. C. Petruccelli and M. A. Alonso, “Phase space distributions tailored for dispersive
media,” J. Opt. Soc. Am. A **27**, 1194‒ 1201 (2010).
[CrossRef]

**177. **Generalized Wigner Functions, Ph.D. Thesis, University of Rochester
(Rochester, NY, 2010).

**178. **K. B. Wolf, M. A. Alonso, and G. W. Forbes, “Wigner functions for Helmholtz wave
fields,” J. Opt. Soc. Am. A **16**, 2476‒ 2487 (1999).
[CrossRef]

**179. **M. A. Alonso, “Radiometry and wide-angle wave fields. I. Coherent
fields in two dimensions,” J. Opt. Soc. Am. A **18**, 902‒ 909 (2001).
[CrossRef]

**180. **C. J. R. Sheppard and K. G. Larkin, “Wigner function for nonparaxial wave
fields,” J. Opt. Soc. Am. A **18**, 2486‒ 2490 (2001).
[CrossRef]

**181. **C. J. R. Sheppard and K. G. Larkin, “Wigner function for highly convergent
three-dimensional wave fields,” Opt. Lett. **26**, 968‒ 970 (2001).
[CrossRef] [PubMed]

**182. **M. A. Alonso, “Radiometry and wide-angle wave fields. II. Coherent
fields in three dimensions,” J. Opt. Soc. Am. A **18**, 910‒ 918 (2001).
[CrossRef]

**183. **M. A. Alonso, “Radiometry and wide-angle wave fields. III. Partial
coherence,” J. Opt. Soc. Am. A. **18**, 2502‒ 2511 (2001).
[CrossRef]

**184. **M. A. Alonso, “Exact description of free electromagnetic wave
fields in terms of rays,” Opt. Express **11**, 3128‒ 3135 (2003).
[CrossRef] [PubMed]

**185. **M. A. Alonso, “Wigner functions for nonparaxial, arbitrarily
polarized electromagnetic wave fields in free-space,” J.
Opt. Soc. Am. A. **21**, 2233‒ 2243 (2004).
[CrossRef]

**186. **J. C. Petruccelli, N, J. Moore, and M. A. Alonso, “Two methods for modeling the propagation of the
coherence and polarization properties of nonparaxial fields,”
Opt. Commun. **283**, 4457‒ 4466 (2010).
[CrossRef]

**187. **J. C. Petruccelli and M. A. Alonso, “Ray-based propagation of the cross-spectral
density,” J. Opt. Soc. Am. A **25**, 1395‒ 1405 (2008).
[CrossRef]

**188. **J. C. Petruccelli and M. A. Alonso, “Propagation of partially coherent fields through
planar dielectric boundaries using angle-impact Wigner functions I. Two
dimensions,” J. Opt. Soc. Am. A **24**, 2590‒ 2603 (2007).
[CrossRef]

**189. **J. C. Petruccelli and M. A. Alonso, “Propagation of nonparaxial partially coherent fields
across interfaces using generalized radiometry,” J. Opt.
Soc. Am. A **26**, 2012‒ 2022 (2009).
[CrossRef]

**190. **J. C. Petruccelli and M. A. Alonso, “Generalized radiometry model for the propagation of
light within anisotropic and chiral media,” J. Opt. Soc.
Am. A **28**, 791‒ 800 (2011).
[CrossRef]

**191. **S. Cho and M. A. Alonso, “Ambiguity function and phase-space tomography for
nonparaxial fields,” J. Opt. Soc. Am. A **28**, 897‒ 902 (2011).
[CrossRef]

**192. **G. I. Ovchinnikov and V. I. Tatarskii, “On the problem of the relationship between coherence
theory and the radiation-transfer equation,” Radiophys.
Quantum Electron. **15**, 1087‒ 1089 (1972).
[CrossRef]

**193. **H. M. Pedersen, J. H. Eberly, L. Mandel, and E. Wolf, “Exact geometrical description of free space
radiative energy transfer for scalar wavefields,” in
*Coherence and Quantum Optics VI*
(Plenum, 1990), pp.
883‒ 887.

**194. **H. M. Pedersen, “Exact theory of free-space radiative energy
transfer,” J. Opt. Soc. Am. A **8**, 176‒ 185 (1991)
errata, **8**, 1518 (1991). [CrossRef]

**195. **H. M. Pedersen, “Geometric theory of fields radiated from
three-dimensional, quasi-homogeneous sources,” J. Opt.
Soc. Am. A **9**, 1626‒ 1632 (1992).
[CrossRef]

**196. **R. G. Littlejohn and R. Winston, “Corrections to classical
radiometry,” J. Opt. Soc. Am. A **10**, 2024‒ 2037 (1993).
[CrossRef]

**197. **S. Cho, J. C. Petruccelli, and M. A. Alonso, “Wigner functions for paraxial and nonparaxial
fields,” J. Mod. Opt. **56**, 1843‒ 1852 (2009).
[CrossRef]

**198. **M. A. Alonso, T. Setälä, and A. T. Friberg, “Optimal pulses for arbitrary dispersive
media,” J. Eur. Opt. Soc. R.P. **6**, 1100 (2011).

**199. **A. Papandreou-Suppappola, F. Hlawatsch, and G. Boudreaux-Bartels, “Quadratic time–frequency representations with scale
covariance and generalized time-shift covariance: A unified framework for the
affine, hyperbolic, and power classes,” Digital Signal
Processing **8**, 3‒ 48 (1998). [CrossRef]

**200. **F. Hlawatsch, A. Papandreou-Suppappola, and G. Boudreaux-Bartels, “The power classes-quadratic time–frequency
representations with scale covariance and dispersive time-shift
covariance,” IEEE Trans. Signal Process. **47**, 3067‒ 3083 (1999).
[CrossRef]

**201. **A. Papandreou-Suppappola, R. L. Murray, B.-G. Iem, and G. F. Boudreaux-Bartels, “Group delay shift covariant quadratic time–frequency
representations,” IEEE Trans. Signal Process. **49**, 2549‒ 2564 (2001).
[CrossRef]

**202. **J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,”
Science **312**, 1780‒ 1782 (2006).
[CrossRef] [PubMed]

**203. **R. W. Robinett, “Quantum wave packet revivals,”
Phys. Rep. **392**, 1‒ 119 (2004).
[CrossRef]

**Miguel A. Alonso** received the degree of Engineer in Physics in 1990 from
the Universidad Autonoma Metropolitana in Mexico City and a Ph.D. in Optics in 1996 from
The Institute of Optics, University of Rochester. After
postdoctoral work at Macquarie University in Sydney, Australia, and a faculty position
at the National Autonomous University of Mexico (UNAM) in Cuernavaca, Mexico, he joined
the faculty of The Institute of Optics in 2003, where he is currently an Associate
Professor. His main research interest is in finding new mathematical models to describe
the propagation of waves through a variety of systems. He is a Deputy Editor of
*Optics Express*, Chair of Spotlight on Optics, and a fellow of the
Optical Society of America.