Mathematical Fundamentals of Modern Linear Optics

Abstract

All known quantum-mechanical approaches to wave and statistical optics are united into a single theory, using Feynman's path integral as a fundamental principle. In short-wave approximations, this principle, the Fourier transformations, and concepts of the theory reproduce Fermat's principle, the Legendre transformations, and concepts of Hamilton's optics and radiometry in a one-to-one fashion.

1. Introduction

There is a fruitful analogy between mechanics and optics. For example, geometrical optics and classical mechanics can be expressed in the same mathematical language of Hamilton’s theory [1]. Moreover, this analogy played a substantial role in the development of wave and statistical optics as a theory that has the same relations to geometric optics and radiometry as quantum mechanics has to classical mechanics. The advent of quantum mechanics radically changed the classical way of physical thinking [2]. Physics that incorporates elements of quantum mechanics is said to be modern physics [3]. Thus, classical optics based on mathematical fundamentals of quantum and statistical mechanics can be called “modern linear optics.” The adjective “linear” emphasizes its difference from modern optics that embraces quantum optical phenomena.

Although the traditional wave optics is based on Maxwell’s electrodynamics [4], in the literature there are four approaches to wave optics which are based on the analogy with quantum mechanics.(1)Generalizing the Huygens-Fresnel principle [5], Feynman has created a new mathematical tool: the path integral [6] which can be used as a fundamental principle of wave and statistical optics [7]. In the stationary-phase approximation Feynman’s integral over all paths transforms into an integral over one path which is satisfied by Fermat’s principle [7], that is, the fundamental principle of Hamilton’s optics.(2)It is known that the Legendre transformations are a framework of Hamilton’s optics [1, 4, 8], and the Fourier transformations are a framework of wave optics [9–14]. Using the stationary phase approximation Walther proved [12–14] that these transformations which before were considered independent, are connected to each other: if the complex wave functions are connected by the Fourier transformations, then the phases of these functions are connected by the Legendre transformation.(3)It is known that geometrical optics is only a short-wave approximation of wave optics and cannot describe optical phenomena as completely as wave optics does, because all equations of geometrical optics can be derived from the wave equation by neglecting certain terms. But Gloge and Marcuse [15, 16] proved that the wave equation can indeed be recovered from the Hamilton equation of geometrical optics by quantization, that is, by formal replacing of all of the variables by corresponding operators.(4)In 1968 Walther proved [17] that the two divisions of optics, which were considered before independently—the theory of the partial spatial coherence (statistical optics) and radiometry of non-Lambertian sources—have in fact the general subject for study, namely, the radiation energy transfer. The basic laws of radiometry can be derived from the theory of partial coherence by using the Wigner definition function (WDF) [17, 18].

Note that according to Bohr’s correspondence principle [19] modern optics in short-wave approximation must reproduce Hamilton’s optics and radiometry.

The purpose of the present paper is to unite these four approaches to wave and statistical optics in a single theory, modern optics. In the short-wave approximation the Huygens-Fresnel principle (in form of Feynman’s path integral), the Fourier transformations, and concepts (e.g., WDF) of this modern optics must reproduce Fermat’s principle, Legendre transformations, and concepts (e.g., radiance) of Hamilton’s optics [8] and radiometry [20] in a one-to-one fashion.

2. Feynman’s Path Integral Reproduces Fermat’s Principle

It is known that all results of geometrical optics can be derived from only one principle—the Fermat principle [1, 4, 8], and all results of wave optics (in the scalar quasi-monochromatic approximation) can be derived from the another principle—the Huygens-Fresnel principle, which can be mathematically expressed in the form of Feynman’s path integral [5]. Let us show that the Fermat principle can be deduced from the Feynman path integral by the method of stationary phase [5, 7].

In a scalar quasi-monochromatic approach, that is, neglecting polarisation and spectral structure of radiation, the distribution of a wave field in the plane 𝑧 can be described by the complex amplitude 𝑈(𝑥,𝑦;𝑧). Consider the propagation of a scalar quasi-monochromatic field with wavelength 𝜆_𝑜 in vacuum along the 𝑧-direction through a layer of an inhomogeneous optical medium with the refractive index distribution 𝑛(𝑥,𝑦,𝑧) restricted by the input {(𝑥,𝑦)}={(𝑥,𝑦,𝑧)∶𝑧=0} and the output {(𝑥^′,𝑦^′)}={(𝑥,𝑦,𝑧)∶𝑧=𝑍} planes (see Figure 1). This layer can be described as a linear system which transforms the input complex amplitude 𝑈(𝑥,𝑦) into the output complex amplitude 𝑈^′(𝑥^′,𝑦^′) by the superposition integral𝑈^′𝑥^′,𝑦^′=_𝑅2𝑈𝑥,𝑦⋅𝐾_𝑆𝑥,𝑦;𝑥^′,𝑦^′𝑑𝑥𝑑𝑦. (1) The kernel 𝐾_𝑆(𝑥,𝑦;𝑥^′,𝑦^′) of this input-output relation is called an optical propagator. Following Feynman [6, 7], the propagator can be written as the integral over all possible virtual paths connecting points (𝑥,𝑦) and (𝑥^′,𝑦^′):𝐾_𝑆𝑥,𝑦;𝑥^′,𝑦^′=_allpath𝐴⋅exp𝑖𝑘_𝑜𝑆_path𝑥,𝑦;𝑥^′,𝑦^′𝐷path, (2) where the symbol 𝐷(path) is the “path differential measure” signifying integration over all paths, 𝑆_path(𝑥,𝑦;𝑥^′,𝑦^′) is called the optical length of the virtual path connecting points (𝑥,𝑦) and (𝑥^′,𝑦^′), 𝑘_𝑜≡2𝜋/𝜆_𝑜is the wave number of the scalar quasi-monochromatic wave in vacuum, and 𝐴 is a normalization constant.

Figure 1 

Illustration of the path integral applied to a layer of an inhomogeneous medium. The path integral contains three elements. (1) A collection of virtual paths connecting an input point (𝑥,𝑦) with an output point (𝑥′,𝑦′) can be labelled with a parameter 𝑢∶𝛾∶(𝑧,𝑢)→(𝑥,𝑦;𝑧)∈𝑅3, that is, 𝛾𝑢(0)=(𝑥,𝑦;0)=(𝑥,𝑦) and 𝛾𝑢(𝑍)=(𝑥,𝑦;𝑍)=(𝑥′,𝑦′). (2) All virtual paths going through the optical inhomogeneous medium 𝑛(𝑥,𝑦,𝑧) from the input point (𝑥,𝑦) to the output point (𝑥′,𝑦′) have the optical length 𝑆∶𝑢→∫𝛾𝑢𝑛(𝑥,𝑦,𝑧)1+𝑥̇2+𝑦2̇𝑑𝑧. (3) The path integral 𝐾𝑆(𝑥,𝑦;𝑥′,𝑦′)=∫∞−∞𝐴⋅exp[𝑖𝑘𝑜𝑆𝑢(𝑥,𝑦;𝑥′,𝑦′)]𝑑𝑢 is the optical propagator. The phase function 𝑆𝑢(𝑥,𝑦;𝑥′,𝑦′) in the neighbourhood of the stationary point  𝑢, where 𝜕S𝑢/𝜕𝑢|𝑢=𝑢=0 can be approximated by a parabola with the vertex at point (𝑢,𝑆𝑢). In this case the optical propagator can be approximated by the function 𝐾𝑆(𝑥,𝑦;𝑥′,𝑦′)=𝐶⋅|𝐾𝑆(𝑥,𝑦;𝑥′,𝑦′)|⋅exp[𝑖𝑘𝑜𝑆𝑢(𝑥,𝑦;𝑥′,𝑦′)], where 𝑆𝑢(𝑥′,𝑦′;𝑥,𝑦) is a point eikonal, that is, the optical length of the path which satisfies Fermat’s principle.

Let us consider the simplest case. The collection of all virtual paths connecting points (𝑥,𝑦) and (𝑥^′,𝑦^′) can be labelled with a parameter 𝑢 [8]. All these paths can be projected onto the 𝑧 axis so that they can be described by functions 𝑥(𝑧) and 𝑦(𝑧). Then the entire assembly of paths running from the input point (𝑥,𝑦) to the output point (𝑥^′,𝑦^′) is given by two functions of two parameters: 𝛾_𝑢(𝑧)={𝑥(𝑧,𝑢),𝑦(𝑧,𝑢)} (see Figure 1). Following Luneburg [1], we use the 𝑧-axis as an independent variable, similar to the time axis in analytical mechanics. So the direction in space can be characterized by “velocities”:𝑥̇=𝑑𝑥𝑑𝑧=tan𝜃⋅sin𝜑,𝑦̇=𝑑𝑦𝑑𝑧=tan𝜃⋅cos𝜑, (3) in terms of the angles (𝜃,𝜑) of the spherical coordinate system (see Figure 2). In this case the optical length of the path takes the following form 𝑆_𝑢(𝑥,𝑦;𝑥^′,𝑦^′):𝑆_𝑢𝑥,𝑦;𝑥^′,𝑦^′=_𝛾𝑢𝐿𝑥,𝑦;𝑥̇,𝑦;̇𝑧𝑑𝑧, (4) where 𝐿𝑥,𝑦;𝑥̇,𝑦;̇𝑧=𝑛𝑥,𝑦,𝑧1+𝑥̇²+𝑦²̇ (5) is the optical Lagrangian [1, 8]. In this case the optical propagator (2) can be rewritten as𝐾_𝑆𝑥,𝑦;𝑥^′,𝑦^′=^∞_−∞𝐴⋅exp𝑖𝑘_𝑜𝑆_𝑢𝑥,𝑦;𝑥^′,𝑦^′𝑑𝑢. (6)

Figure 2 

Direction of the tangent to a path 𝛾𝑢(𝑧) at a given point 𝑃∗(𝛾𝑢(𝜁),𝜁) can be characterized by angles (𝜃,𝜑) of the spherical coordinate system (𝜑∈(0,2𝜋),𝜃∈(0,𝜋/2)), or by “velocities” (𝑥̇,𝑦)̇.

Note that for sufficiently large 𝑘_𝑜 (the short-wavelength approximation), the exponential functions exp[𝑖𝑘_𝑜𝑆_𝑢(𝑥,𝑦;𝑥^′,𝑦^′)] fluctuate so rapidly that they compensate each other, and the final amplitude becomes close to zero. Compensation does not occur in the neighbourhood of the stationary point 𝑢 (see Figure 1):𝜕𝑆_𝑢𝜕𝑢_𝑢=𝑢=0, (7) where the phase function 𝑆_𝑢(𝑥,𝑦;𝑥^′,𝑦^′) varies relatively slowly, the exponential function constructively interferes, and the final amplitude, the optical propagator equation (6), is rather large.

Since outside the neighbourhood of the stationary point 𝑢 the resulting amplitude nearly equals zero the function 𝑆_𝑢(𝑥,𝑦;𝑥^′,𝑦^′) can be approximated by a parabola with the vertex at point (𝑢,𝑆_𝑢) (see Figure 1):𝑆_𝑢𝑥,𝑦;𝑥^′,𝑦^′≈𝑆_𝑢𝑥,𝑦;𝑥^′,𝑦^′+12𝛼(𝑢−𝑢)², (8) where 𝜕𝛼≡²𝑆_𝑢𝑥,𝑦;𝑥^′,𝑦^′𝜕𝑢²_𝑢=𝑢. (9) Thus, in (6) we can substitute the phase function 𝑆_𝑢(𝑥,𝑦;𝑥^′,𝑦^′) by its approximation (8):𝐾_𝑆𝑥,𝑦;𝑥^′,𝑦^′=𝐴⋅exp𝑖𝑘_𝑜𝑆_𝑢𝑥,𝑦;𝑥^′,𝑦^′⋅^∞_−∞exp𝑖𝑘_𝑜𝛼2(𝑢−𝑢)²𝑑𝑢. (10) As the integral has the closed form [10, 21, 22]^∞_−∞exp𝑖𝑘_𝑜𝛼2(𝑢−𝑢)²𝑑𝑢=2𝜋𝑘_𝑜𝛼𝑖𝜋exp4. (11) Equation (8) can be rewritten as𝐾_𝑆𝑥,𝑦;𝑥^′,𝑦^′=𝐾_𝑆𝑥,𝑦;𝑥^′,𝑦^′⋅exp𝑖𝑘_𝑜𝑆_𝑢𝑥,𝑦;𝑥^′,𝑦^′. (12) Equation (12) is the stationaryphase approximation of the propagator (6) [5].

Note that (12) creates the impression that the wave runs along the path 𝛾_𝑢(𝑧) whose optical length is stationary with respect to small perturbations of 𝑢, Equation (6), that is, whose optical length satisfies Fermat’s principle and the phase function 𝑆_𝑢(𝑥^′,𝑦^′;𝑥,𝑦) is a point eikonal [5]. Thus, using the method of stationary phase, we have shown that Feynman’s path integral reproduces Fermat’s principle and the point eikonal (see sections 2 and  4 in [8]).

3. 𝑘_𝑜-Fourier Transformations Reproduce the Legendre Transformations

For true understanding of problems, it is necessary to consider them from various points of view. This possibility is given by one-to-one mappings such as the Legendre transformations in geometrical optics and the Fourier transformations in wave optics. Note that to understand optics as a whole, one needs to know not only these transformations alone, but also the one-to-one mapping between them which is given by the method of stationary phase [5, 12, 13, 21, 22].

The direct 𝑘_𝑜-Fourier transformation of the complex-valued function 𝑈(𝑥,𝑦) of real variables (𝑥,𝑦) is a new complex-valued function 𝑈(𝑝,𝑞) of new real variables (𝑝,𝑞), which is defined as𝑈=𝑝,𝑞_𝑘𝑜𝐹_{𝑥→𝑝𝑦→𝑞}𝑈≡𝑘𝑥,𝑦_𝑜2𝜋𝑖_𝑅²𝑈𝑥,𝑦⋅exp−𝑖𝑘_𝑜𝑥⋅𝑝+𝑦⋅𝑞𝑑𝑥𝑑𝑦. (13a)

The inverse 𝑘_𝑜-Fourier transformation of the complex-valued function 𝑈(𝑝,𝑞) of real variables (𝑝,𝑞) is a new complex-valued function 𝑈(𝑥,𝑦) of new real variables (𝑥,𝑦), which is defined as𝑈=𝑥,𝑦_𝑘𝑜𝐹⁻¹_{𝑝→𝑥𝑞→𝑦}𝑈≡−𝑘𝑝,𝑞_𝑜2𝜋𝑖_𝑅²𝑈𝑝,𝑞⋅exp𝑖𝑘_𝑜𝑥⋅𝑝+𝑦⋅𝑞𝑑𝑝𝑑𝑞. (13b)

The composition of the direct 𝑘_𝑜-Fourier transformation and the inverse 𝑘_𝑜-Fourier transformation (or the composition of the inverse 𝑘_𝑜-Fourier transformation and the direct 𝑘_𝑜-Fourier transformation) is the identity transformation 𝐼:^𝑘𝑜𝐹_{𝑥→𝑝𝑦→𝑞}∘^𝑘𝑜𝐹⁻¹_{𝑝→𝑥𝑞→𝑦}=𝐼_𝑘𝑜𝐹⁻¹_{𝑝→𝑥𝑞→𝑦}∘^𝑘𝑜𝐹_{𝑥→𝑝𝑦→𝑞}=𝐼. (14)

If the complex functions are connected by the 𝑘_𝑜-Fourier transformations, the phases of these functions are connected by the Legendre transformations [21, 22].

Note that if the complex functions 𝑈(𝑥,𝑦) and 𝑈(𝑝,𝑞) are factorized into the amplitude and phase factors:𝑈𝑥,𝑦=𝑈𝑥,𝑦⋅exp𝑖𝑘_𝑜Ψ,𝑈=𝑈𝑥,𝑦𝑝,𝑞𝑝,𝑞⋅exp𝑖𝑘_𝑜Ψ𝑝,𝑞, (15) then their direct and inverse 𝑘_𝑜-Fourier transformations (13a), and (13b) take the forms𝑈=𝑘𝑝,𝑞^𝑜2𝜋𝑖_𝑅²𝑈𝑥,𝑦⋅exp𝑖𝑘_𝑜Ψ,𝑈=−𝑘𝑥,𝑦−𝑥⋅𝑝−𝑦⋅𝑞𝑑𝑥𝑑𝑦𝑥,𝑦^𝑜2𝜋𝑖_𝑅²𝑈𝑝,𝑞⋅exp𝑖𝑘_𝑜Ψ𝑝,𝑞+𝑥⋅𝑝+𝑦⋅𝑞𝑑𝑝𝑑𝑞. (16)

The phase functions Ψ(𝑥,𝑦)−𝑥𝑝−𝑦𝑞 and Ψ(𝑝,𝑞)+𝑥𝑝+𝑦𝑞 have stationary points (𝑝,𝑞) and (𝑥,𝑦), respectively, that is, 𝜕Ψ𝑥,𝑦−𝑥⋅𝑝−𝑦⋅𝑞𝜕Ψ𝜕𝑥=0,𝑥,𝑦−𝑥⋅𝑝−𝑦⋅𝑞⟹=𝜕𝑦=0𝑝,𝑞𝜕Ψ𝑥,𝑦,𝜕𝑥𝜕Ψ𝑥,𝑦𝜕Ψ𝜕𝑦𝑝,𝑞+𝑥⋅𝑝+𝑦⋅𝑞𝜕Ψ𝜕𝑝=0,𝑝,𝑞+𝑥⋅𝑝+𝑦⋅𝑞⟹=−𝜕Ψ𝜕𝑦=0𝑥,𝑦𝑝,𝑞𝜕Ψ𝜕𝑝,−𝑝,𝑞.𝜕𝑞 (17)

These stationary points are nondegenerate: det^{Ψ′′𝑥𝑥Ψ′′𝑥𝑦}_{Ψ′′𝑦𝑥Ψ′′𝑦𝑦}≠0 and Ψdet^{′′𝑝𝑝}Ψ^′′_𝑝𝑞Ψ^′′_𝑞𝑝Ψ^′′_𝑞𝑞≠0.

Using the stationaryphase approximation (see Appendix A), the Fourier transformations can be estimated by the equations^𝑘𝑜𝐹_{𝑥→𝑝𝑦→𝑞}𝑈𝑥,𝑦exp𝑖𝑘_𝑜Ψ≈𝑥,𝑦_{𝑘𝑜→∞}Ψdet^{′′𝑝𝑝}Ψ^′′_𝑝𝑞Ψ^′′_𝑞𝑝Ψ^′′_𝑞𝑞^1/2⋅𝑈−𝜕Ψ𝜕Ψ𝜕𝑝,−𝜕𝑞⋅exp𝑖𝑘_𝑜Ψ𝑝,𝑞_𝑘𝑜,𝐹⁻¹_{𝑝→𝑥𝑞→𝑦}𝑈𝑝,𝑞exp𝑖𝑘_𝑜Ψ≈𝑝,𝑞_{𝑘𝑜→∞}Ψdet^{′′𝑥𝑥}Ψ^′′_𝑥𝑦Ψ^′′_𝑦𝑥Ψ^′′_𝑦𝑦^1/2⋅𝑈𝜕Ψ,𝜕𝑥𝜕Ψ𝜕𝑦⋅exp𝑖𝑘_𝑜Ψ𝑥,𝑦. (18a)(18b)

The phase function Ψ(𝑝,𝑞) is the direct Legendre transformation of the continuously differentiable real-valued phase function Ψ(𝑥,𝑦), that is, a new real-valued phase function Ψ(𝑝,𝑞) of new real variables (𝑝,𝑞), which is defined asΨ𝑝,𝑞=𝐿_{𝑥→𝑝𝑦→𝑞}Ψ𝑥,𝑦≡Ψ𝑥,𝑦−𝑥𝑝−𝑦𝑞,where𝑝≡𝜕Ψ𝜕𝑥,𝑞≡𝜕Ψ.𝜕𝑦 (19a) The phase function Ψ(𝑥,𝑦) is the direct Legendre transformation of the continuously differentiable real-valued phase function, that is, a new real-valued phase function Ψ(𝑝,𝑞) of new real variables (𝑥,𝑦), which is defined asΨ𝑥,𝑦=𝐿⁻¹_{𝑝→𝑥𝑞→𝑦}Ψ≡Ψ𝜕Ψ𝑝,𝑞𝑝,𝑞+𝑥𝑝+𝑦𝑞,where𝑥≡𝜕Ψ𝜕𝑝,𝑦≡.𝜕𝑞 (19b)

If the function Ψ(𝑥,𝑦) is positive definite, det^{Ψ′′𝑥𝑥Ψ′′𝑥𝑦}_{Ψ′′𝑦𝑥Ψ′′𝑦𝑦}>0, (or negative definite, det^{Ψ′′𝑥𝑥Ψ′′𝑥𝑦}_{Ψ′′𝑦𝑥Ψ′′𝑦𝑦}<0), then the function Ψ(𝑝,𝑞) is negative definite, Ψdet^{′′𝑝𝑝}Ψ^′′_𝑝𝑞Ψ^′′_𝑞𝑝Ψ^′′_𝑞𝑞<0 (or positive definite, Ψdet^{′′𝑝𝑝}Ψ^′′_𝑝𝑞Ψ^′′_𝑞𝑝Ψ^′′_𝑞𝑞>0). In this caseΨdet^{′′𝑥𝑥}Ψ^′′_𝑥𝑦Ψ^′′_𝑦𝑥Ψ^′′_𝑦𝑦Ψ⋅det^{′′𝑝𝑝}Ψ^′′_𝑝𝑞Ψ^′′_𝑞𝑝Ψ^′′_𝑞𝑞=−1. (20) If the functions Ψ(𝑥,𝑦) and Ψ(𝑝,𝑞) are positive or negative definite functions, then the composition of the direct Legendre transformation and the inverse Legendre transformation (or the composition of the inverse Legendre transformation and the direct Legendre transformation) is the identity transformation 𝐼:𝐿_{𝑥→𝑝𝑦→𝑞}∘𝐿⁻¹_{𝑝→𝑥𝑞→𝑦}=𝐼,…𝐿⁻¹_{𝑝→𝑥𝑞→𝑦}∘𝐿_{𝑥→𝑝𝑦→𝑞}=𝐼. (21)

4. Dual Meaning of Spatial Frequencies

In Fourier optics, the spatial frequency has a dual meaning. On the one hand it is a sinusoidal grating, and on the other hand it is the momentum of the monochromatic plane wave which goes through this grating. To analyze this in detail, let us consider the simplest optical system, that is, a layer of an optically homogeneous medium with refractive index 𝑛 and thickness 𝑍 bounded by parallel input (𝑥,𝑦)=(𝑥,𝑦,0) and output planes (𝑥_𝑧,𝑦_𝑧)=(𝑥,𝑦,𝑍) (see Figure 3).

Figure 3 

A layer of an optically homogeneous medium as a linear system.

It is known that a scalar quasi-monochromatic plane wave travelling through this layer in a direction specified by the vector (𝑝,𝑞,𝑛²−𝑝²−𝑞²) (see Figure 4) can be described by the equationexp𝑖𝑘_𝑜𝑥𝑝+𝑦𝑞+𝑧𝑛²−𝑝²−𝑞². (22)

Figure 4 

Direction of the “ray vector” can be characterized by angles (𝜃,𝜑) of the spherical coordinate system (𝜑∈(0,2𝜋), 𝜃∈(0,𝜋/2)), or by momenta (𝑝,𝑞).

Using the direction cosines 𝑝 and 𝑞 as two independent variables, we can describe the complex amplitude on a plane 𝑍, 𝑈(𝑥,𝑦;𝑍), as a linear superposition of plane waves with weight function 𝑈_in(𝑝,𝑞) [9–13]:𝑈=−𝑘𝑥,𝑦;𝑍_𝑜2𝜋𝑖_𝑅²𝑈_in𝑝,𝑞exp𝑖𝑘_𝑜𝑥𝑝+𝑦𝑞+𝑍𝑛²−𝑝²−𝑞²𝑑𝑝𝑑𝑞. (23)

Let us consider the layer of a homogeneous optical medium as a linear system (see Figure 3) [9–13] which transforms the input complex amplitude 𝑈_in(𝑥,𝑦)≡𝑈(𝑥,𝑦;0) to the output complex amplitude 𝑈_out(𝑥_𝑧,𝑦_𝑧)≡𝑈(𝑥,𝑦;𝑍). Then using the inverse 𝑘_𝑜-Fourier transformation (14), (24a) and (24b) can be rewritten twice𝑈_in=𝑥,𝑦_𝑘𝑜𝐹⁻¹_{𝑝→𝑥𝑞→𝑦}𝑈_in𝑈𝑝,𝑞,if𝑧=0,_out𝑥_𝑧,𝑦_𝑧=^𝑘𝑜𝐹⁻¹_{𝑝→𝑥𝑧𝑞→𝑦𝑧}𝑈_in𝑝,𝑞⋅exp−𝑖𝑘_𝑜𝑍⋅𝐻𝑝,𝑞,if𝑧=𝑍, (24a)(24b)

where the function𝐻𝑝,𝑞=−𝑛²−𝑝²−𝑞² (25) is known as the optical Hamiltonian [1, 5, 8].

Thus, the weight function 𝑈_in(𝑝,𝑞)=^𝑘𝑜𝐹_{𝑥→𝑝𝑦→𝑞}{𝑈_in(𝑥,𝑦)} can be interpreted as an input frequency spectrum, and the variables (𝑝,𝑞)∈𝑅²—as spatial frequencies. According to ((24b), the input frequency spectrum and the output frequency spectrum 𝑈_out(𝑝,𝑞)=^𝑘𝑜𝐹_{𝑥𝑧𝑦→𝑝𝑧→𝑞}{𝑈_out(𝑥_𝑧,𝑦_𝑧)} are related by the equation𝑈_out𝑝,𝑞=𝑈_in𝑝,𝑞⋅𝐾_𝑍𝑝,𝑞, (26) where 𝐾_𝑍(𝑝,𝑞)=exp[−𝑖𝑘_𝑜𝑍⋅𝐻(𝑝,𝑞)] is the transfer function of the layer of a homogeneous optical medium with refractive index 𝑛 of thickness 𝑍.

Thus, in region 𝑝²+𝑞²≤𝑛² variables (𝑝,𝑞) can be interpreted not only as spatial frequencies but also as the direction cosines of the wave vector. This manifests the so called dual meaning of spatial frequencies [11].

According to the convolution theorem [9, 10], Equation (26) can be rewritten as𝑈_out𝑥_𝑧,𝑦_𝑧=_𝑅2𝑈_in𝑥_𝑧−𝑍𝑥̇,𝑦_𝑧=−𝑍𝑦̇⋅𝐺𝑥̇,𝑦̇𝑑𝑥̇𝑑𝑦_𝑅2̇𝑈_in𝑥,𝑦⋅𝐺_𝑧𝑥_𝑧−𝑥,𝑦_𝑧−𝑦𝑑𝑥𝑑𝑦, (27) where𝐺_𝑧𝑥_𝑧−𝑥,𝑦_𝑧=1−𝑦𝑍²𝐺𝑥_𝑧−𝑥𝑍,𝑦_𝑧−𝑦𝑍 (28a) is the impulse response function of the layer of a homogeneous optical medium with refractive index 𝑛 of thickness 𝑍, and𝐺≡𝑘𝑥̇,𝑦̇^𝑜𝑍⋅2𝜋𝑖^𝑘𝑜𝑍𝐹⁻¹_{𝑝→̇𝑥𝑞→̇𝑦}exp−𝑖𝑘_𝑜𝑍⋅𝐻𝑝,𝑞 (28b)

is the impulse response function of the layer of a homogeneous optical medium with refractive index 𝑛 of unit thickness. Using the stationaryphase approximation of the inverse Fourier transformation (18b), the impulse response function (28b) can be rewritten as:𝐺≡𝑘𝑥̇,𝑦̇^𝑜𝑍⋅2𝜋𝑖^𝑘𝑜𝑍𝐹⁻¹_{𝑝→̇𝑥𝑞→̇𝑦}exp−𝑖𝑘_𝑜𝑍≈⋅𝐻𝑝,𝑞_{𝑘𝑜𝑍→∞}𝑘_𝑜𝑍⋅𝐿2𝜋𝑖det^{′′̇𝑥̇𝑥}𝐿^′′_̇𝑥̇𝑦𝐿^′′_̇𝑦̇𝑥𝐿^′′_̇𝑦̇𝑦^1/2⋅exp𝑖𝑘_𝑜𝑍⋅𝐿𝑥̇,𝑦̇. (29) Here𝐿𝑥̇,𝑦̇=−𝐿⁻¹_{𝑝→̇𝑥𝑞→̇𝑦}𝐻𝑝,𝑞=𝑥̇⋅𝑝+𝑦̇⋅𝑞−𝐻𝑝,𝑞, (30) where𝑥̇=𝜕𝐻𝑝,𝑞𝜕𝑝,𝑦̇=𝜕𝐻𝑝,𝑞.𝜕𝑞 (31) Note that (31) can be interpreted as Hamilton’s equations which describe propagation of a light ray through an inhomogeneous optical medium with the refractive index distribution 𝑛(𝑥,𝑦,𝑧).

In the optically homogeneous medium (𝑛=const) the Legendre transformation (30) of the Hamiltonian (25) is the optical Lagrangian 𝐿(𝑥̇,𝑦)̇=𝑛1+𝑥̇²+𝑦²̇ [1, 8, 15, 16]. Since𝐿det^{′′̇𝑥̇𝑥}𝐿^′′_̇𝑥̇𝑦𝐿^′′_̇𝑦̇𝑥𝐿^′′_̇𝑦̇𝑦𝑛=det1+𝑦̇1+𝑥̇²+𝑦̇^23/2𝑛−𝑥̇𝑦̇1+𝑥̇²+𝑦̇^23/2𝑛−𝑥̇𝑦̇1+𝑥̇²+𝑦̇^23/2𝑛1+𝑥̇1+𝑥̇²+𝑦̇^23/2=𝑛1+𝑥̇²+𝑦̇²², (32) the impulse response function (29) takes the form𝐺𝑥̇,𝑦̇≈𝑛𝑘_𝑜𝑍2𝜋𝑖exp𝑖𝑛𝑘_𝑜𝑍1+𝑥̇²+𝑦²̇1+𝑥̇²+𝑦²̇, (33) or, according to (28a),𝐺_𝑧𝑥_𝑧−𝑥,𝑦_𝑧≈−𝑦𝑛𝑘_𝑜𝑍2𝜋𝑖exp𝑖𝑛𝑘_𝑜𝑍²+𝑥_𝑧−𝑥²+𝑦_𝑧−𝑦²𝑍²+𝑥_𝑧−𝑥²+𝑦_𝑧−𝑦². (34) Let us compare the approximation (34) to the well-known exact solution of (28a), and (28b), the so-called Rayleigh-Sommerfeld impulse response function given by [9–13]𝐺_𝑧𝑥_𝑧−𝑥,𝑦_𝑧1−𝑦=−𝑑2𝜋𝑑𝑍exp𝑖𝑛𝑘_𝑜𝑍²+𝑥_𝑧−𝑥²+𝑦_𝑧−𝑦²𝑍²+𝑥_𝑧−𝑥²+𝑦_𝑧−𝑦²=11−𝑖𝑛𝑘_𝑜𝑍²+𝑥_𝑧−𝑥²+𝑦_𝑧−𝑦²⋅𝑛𝑘_𝑜𝑍×2𝜋𝑖exp𝑖𝑛𝑘_𝑜𝑍²+𝑥_𝑧−𝑥²+𝑦_𝑧−𝑦²𝑍²+𝑥_𝑧−𝑥²+𝑦_𝑧−𝑦². (35) For large values of 𝑍, the second term between the parentheses is smaller than the first term, and in this most important case approximation (34) is very good.

Using the method of stationary phase we have shown that the Fourier transformations between the transfer function and the impulse response function of the layer of a homogeneous optical medium with refractive index 𝑛 of the unit thickness reproduce the Legendre transformations between the Hamiltonian and the Lagrangian (see Section  3 in [8]).

5. Input-Output Transformations of Wave Optics Reproduce Eikonals Theory

From the possibility of both coordinate and frequency descriptions of the wave field at the input and output of the optical system follows the existence of four equivalent ways of wave description of this system in terms of “input-output transformations.” Using the method of stationary phase, Walter showed that each of these transformations of the wave field reproduces the eikonal which describes the transformation in terms of geometrical optics.

Consider a linear optical system bounded by parallel input (𝑥,𝑦) and output (𝑥^′,𝑦^′) planes lying in layers of optically homogeneous media with refractive indices 𝑛 and 𝑛^′, respectively, and a layer of an optically inhomogeneous medium between them (see Figure 5). Variables with primes refer to the output plane.

Figure 5 

Linear optical system.

As the radiation field in the input and output planes of the linear optical system can be specified either in the coordinate representations 𝑈(𝑥,𝑦) and 𝑈^′(𝑥^′,𝑦^′) or in the frequency representations 𝑈(𝑝,𝑞) and 𝑈^′(𝑝^′,𝑞^′), according to the superposition principle, the linear optical system can be described by any of the four integral equations [5, 12–14]:𝑈^′𝑥^′,𝑦^′=_𝑅2𝐾_𝑆𝑥,𝑦;𝑥^′𝑦^′𝑈⋅𝑈𝑥,𝑦𝑑𝑥𝑑𝑦,_′𝑝^′,𝑞^′=_𝑅2𝐾_𝑉𝑥,𝑦;𝑝^′𝑞^′𝑈⋅𝑈𝑥,𝑦𝑑𝑥𝑑𝑦,_′𝑥^′,𝑦^′=_𝑅2𝐾_𝑉′𝑝,𝑞;𝑥^′𝑦^′⋅𝑈𝑈𝑝,𝑞𝑑𝑝𝑑𝑞,_′𝑝^′,𝑞^′=_𝑅2𝐾_𝑇𝑝,𝑞;𝑝^′𝑞^′⋅𝑈𝑝,𝑞𝑑𝑝𝑑𝑞. (36a)(36b)(36c)(36d)

The kernels of these linear input-output transformations can be called the point 𝐾_𝑆(𝑥,𝑦;𝑥^′,𝑦^′), point-angle 𝐾_𝑉(𝑥,𝑦;𝑝^′,𝑞^′), angle-point 𝐾_𝑉′(𝑝,𝑞;𝑥^′,𝑦^′), and angle 𝐾_𝑇(𝑝,𝑞;𝑝^′,𝑞^′) propagators.

The complex amplitudes, 𝑈(𝑥,𝑦) and 𝑈^′(𝑥^′,𝑦^′), and frequency spectra, 𝑈(𝑝,𝑞) and 𝑈^′(𝑝^′,𝑞^′), are connected by 𝑘_𝑜-Fourier transformations. As a result, the propagators of the optical system are connected by the Fourier transformations too [5, 12–14]:𝐾_𝑉𝑥,𝑦;𝑝^′,𝑞^′=^𝑘𝑜𝐹_{𝑥′→𝑝′𝑦′→𝑞′}𝐾_𝑆𝑥,𝑦;𝑥^′𝑦^′,𝐾_𝑉′𝑝,𝑞;𝑥^′,𝑦^′=^𝑘𝑜𝐹⁻¹_{𝑥→𝑝𝑦→𝑞}𝐾_𝑆𝑥,𝑦;𝑥^′𝑦^′,𝐾_𝑇𝑝,𝑞;𝑝^′,𝑝^′=^𝑘𝑜𝐹^𝑘−1_{𝑥→𝑝𝑦→𝑞}^𝑜𝐹_{𝑥′→𝑝′𝑦′→𝑞′}𝐾_𝑆𝑥,𝑦;𝑥^′𝑦^′. (37)

In these expressions, the frequency representations are transformed into the coordinate representations by the direct Fourier transformations in the input plane and by the inverse Fourier transformations in the output plane.

Note that propagators can be factorized into amplitude and phase factors:𝐾_𝑠𝑥,𝑦;𝑥^′,𝑦^′=𝐾_𝑠𝑥,𝑦;𝑥^′,𝑦^′⋅exp𝑖𝑘_𝑜𝑆𝑥,𝑦;𝑥^′,𝑦^′, (38a)𝐾_𝑉𝑥,𝑦;𝑝^′,𝑞^′=𝐾_𝑉𝑥,𝑦;𝑝^′,𝑞^′⋅exp𝑖𝑘_𝑜𝑉𝑥,𝑦;𝑝^′,𝑞^′, (38b)𝐾_𝑉′𝑝,𝑞;𝑥^′,𝑦^′=𝐾_𝑉′𝑝,𝑞;𝑥^′,𝑦^′⋅exp𝑖𝑘_𝑜𝑉^′𝑝,𝑞;𝑥^′,𝑦^′, (38c)𝐾_𝑇𝑝,𝑞;𝑝^′,𝑞^′=𝐾_𝑇𝑝,𝑞;𝑝^′,𝑞^′⋅exp𝑖𝑘_𝑜𝑇𝑝,𝑞;𝑝^′,𝑞^′. (38d)

Since the propagators of the optical system are connected by the Fourier transformations (37), then according to (18a), and (18b), their phase functions are connected by the Legendre transformations:𝑉𝑥,𝑦;𝑝^′,𝑞^′=𝐿_{𝑥′⟶𝑝′𝑦′⟶𝑞′}𝑆𝑥,𝑦;𝑥^′𝑦^′,𝑉^′𝑝,𝑞;𝑥^′,𝑦^′=𝐿⁻¹_{𝑥⟶𝑝𝑦⟶𝑞}𝑆𝑥,𝑦;𝑥^′𝑦^′,𝑇𝑝,𝑞;𝑝^′,𝑞^′=𝐿⁻¹_{𝑥⟶𝑝𝑦⟶𝑞}𝐿_{𝑥′⟶𝑝′𝑦′⟶𝑞′}𝑆𝑥,𝑦;𝑥^′𝑦^′. (39)

Because the phase function 𝑆(𝑥,𝑦;𝑥^′,𝑦^′) is a point eikonal 𝑆_𝑢(𝑥,𝑦;𝑥^′,𝑦^′) (4), the phase functions 𝑉(𝑥,𝑦;𝑝^′,𝑞^′), 𝑉^′(𝑝,𝑞;𝑥^′,𝑦^′), and 𝑇(𝑝,𝑞;𝑝^′,𝑞^′) can be interpreted as the point-angle, angle-point and angle eikonals [8, 12–14].

It is known that the point-angle eikonal of the perfect optical system has the form𝑉^perfect𝑥,𝑦;𝑝^′,𝑞^′=−𝑚𝑥𝑝^′+𝑦𝑞^′, (40) and the point-angle eikonal of the isoplanar (or shift invariant) optical system has the form𝑉^isoplanar𝑥,𝑦;𝑝^′,𝑞^′=−𝑚𝑥𝑝^′+𝑦𝑞^′+𝑊𝑝^′,𝑞^′, (41) where 𝑊(𝑝^′,𝑞^′) is the wave aberration [23, 24]. Note that in an optical system with telecentric stop in object space (see Figure 6) the wave aberration is determined in the region of moments limited by the telecentric stop with the radius 𝑅: 𝑝^′2+𝑞^′2≤𝑅/𝑓, where 𝑓 is the back focal length (distance between front focal 𝐹 and front nodal points 𝑁 of the optical system). In the isoplanar optical system, the point eikonal 𝑆^isoplanar(𝑥,𝑦;𝑥^′,𝑦^′) and the point propagator 𝐾^isoplanar_𝑠(𝑥,𝑦;𝑥^′,𝑦^′) depend only on the difference variables:𝑆^isoplanar𝑥,𝑦;𝑥^′,𝑦^′=𝐿_𝑝⁻¹_{′→𝑥′𝑞′→𝑦′}𝑉^isoplanar𝑥,𝑦;𝑝^′𝑞^′=−𝑚𝑥𝑝^′+𝑦𝑞^′+𝑊𝑝^′,𝑞^′+𝑥^′𝑝^′+𝑦^′𝑞^′=𝑊𝑝^′,𝑞^′+𝑥^′−𝑚𝑥𝑝^′+𝑦^′−𝑚𝑦𝑞=𝐿_𝑝⁻¹_{′→(𝑥′𝑞−𝑚𝑥)′→(𝑦′−𝑚𝑦)}𝑊𝑝^′𝑞^′=𝑆𝑥^′−𝑚𝑥,𝑦^′,𝐾−𝑚𝑦_{isoplanar𝑆}𝑥,𝑦;𝑥^′,𝑦^′=^𝑘𝑜𝐹_𝑝⁻¹_{′→𝑥′𝑞′→𝑦′}𝐾_𝑉𝑥,𝑦;𝑝^′𝑞^′=^𝑘𝑜𝐹_𝑝⁻¹_{′→𝑥′𝑞′→𝑦′}𝐾𝑝^′𝑞^′⋅exp𝑖𝑘_𝑜𝑉^isoplanar𝑥,𝑦;𝑝^′,𝑞^′=−𝑘_𝑜2𝜋𝑖_𝑅²𝐾𝑝^′𝑞^′⋅exp𝑖𝑘_𝑜−𝑚𝑥𝑝^′+𝑦𝑞^′+𝑊𝑝^′,𝑞^′+𝑥^′⋅𝑝^′+𝑦^′⋅𝑞^′𝑑𝑝^′𝑑𝑞^′=−𝑘_𝑜2𝜋𝑖_𝑅²𝐾𝑝^′𝑞^′⋅exp𝑖𝑘_𝑜𝑊𝑝^′,𝑞^′+𝑥^′−𝑚𝑥⋅𝑝^′+𝑦^′−𝑚𝑥⋅𝑞^′𝑑𝑝^′𝑑𝑞^′=^𝑘𝑜𝐹_𝑝⁻¹_{′→(𝑥′𝑞−𝑚𝑥)′→(𝑦′−𝑚𝑦)}𝐾𝑝^′𝑞^′⋅exp𝑖𝑘_𝑜𝑊𝑝^′𝑞^′≡𝐾𝑥^′−𝑚𝑥,𝑦^′−𝑚𝑦. (42)(43)

Figure 6 

An isoplanar optical system with telecentric stop in object space.

In this case the functions 𝐾(𝑝^′𝑞^′) and 𝐾(𝑥^′−𝑚𝑥,𝑦^′−𝑚𝑦) are called the transfer function and the impulse response function of the isoplanar optical system, respectively.

Using the method of stationary phase, we have shown that the Fourier transformations between propagators of optical system (37) reproduce the Legendre transformations between its eikonals (39) (see Section  4 in [8]). In the same way, Fourier transformations between the transfer function and the impulse response function of an isoplanar optical system (43) reproduce the Legendre transformations between its wave aberration function and its point eikonal (42).

6. Quantization: Hamiltonian Reproduces Helmholtz Equation

According to Gloge and Marcuse a formal quantization of the Hamiltonian of geometrical optics leads to the Helmholtz equation of wave optics [15, 16]. Let as show that this formal procedure can be substantiated by using the relations between the Fourier and the Legendre transformations, (18a), and (18b).

Note that the Fourier transformations between input-output relations (37) and the Legendre transformations between eikonals (39) are valid for any position of the input and output planes. Let the input and output planes be separated by a layer of a medium of small thickness 𝜀. Then the differential equation of wave optics (the Helmholtz equation) can be derived from the integral superposition (2).

Consider the particular case where the input plane 𝜁 and the output plane 𝜁+𝜀 are separated by a plane-parallel layer of medium with smooth homogeneities (see Figure 7). Let the thickness 𝜀 be less than the characteristic size of the optical homogeneities, but substantially greater than the radiation wavelength 𝜆_𝑜/𝑛 [5]. In this case the medium in a neighbourhood of a point (𝑥,𝑦) can be considered as quasi-homogeneous and the integral superposition equation (36a) with the propagator (38a) can be approximated by the equation [5, 6]𝑈𝑥^′,𝑦^′𝑘≡𝑈𝑥,𝑦;𝜁+𝜀=const⋅^𝑜2𝜋𝑖_𝑅²exp𝑖⋅𝜀𝑘_𝑜𝐿𝑥̇,𝑦,̇𝑥,𝑦;𝜁×𝑈𝑥−𝜀𝑥̇,𝑦−𝜀𝑦;̇𝜁𝑑𝑥̇𝑑𝑦,̇ (44) because 𝑆(𝑥,𝑦,𝑥^′,𝑦^′)=𝑆(𝑥^′−𝑥,𝑦^′−𝑦)=𝜀𝐿((𝑥^′−𝑥)/𝜀,(𝑦^′−𝑦)/𝜀,𝑥,𝑦).

Figure 7 

Thin layer of an optically quasi-homogeneous medium.

Using the operator representation of the Taylor series (Appendix B), the function 𝑈(𝑥−𝜀𝑥̇,𝑦−𝜀𝑦;̇𝜁) can be given as𝑈𝜕𝑥−𝜀𝑥̇,𝑦−𝜀𝑦;̇𝜁=exp−𝜀𝑥̇𝜕𝜕𝑥−𝜀𝑦̇𝑈.𝜕𝑦𝑥,𝑦;𝜁 (45) Substituting this expression into the left-hand side of (44), we get𝑈𝑘𝑥,𝑦;𝜁+𝜀=const⋅^𝑜.2𝜋𝑖_𝑅²exp𝑖𝜀𝑘_𝑜𝐿−𝑖𝑥̇,𝑦,̇𝑥,𝑦;𝜁−𝑥̇𝑘_𝑜𝜕−𝑖𝜕𝑥−𝑦̇𝑘_𝑜𝜕𝜕𝑦×𝑈𝑥,𝑦;𝜁𝑑𝑥̇𝑑𝑦̇ (46) or, using the definition of the Fourier transformation (13a),𝑈(𝑥,𝑧;𝜁+𝜀)=const⋅^𝑘𝑜𝐹_{̇𝑥→(−𝑖/𝑘𝑜𝜕/𝜕𝑥)̇𝑦→(−𝑖/𝑘𝑜𝜕/𝜕𝑦)}exp𝑖𝜀𝑘_𝑜𝐿.𝑥̇,𝑦,̇𝑥,𝑦;𝜁×𝑈𝑥,𝑦;𝜁 (47)

If 𝐿(𝑥̇,𝑦,̇𝑥,𝑦;𝜁) is a smoothly varying function, then, according to the asymptotic (18a)), the phases of the functions connected by the Fourier transformation are connected by the Legendre transformation:^𝑘𝑜𝐹_{̇𝑥→(−𝑖/𝑘𝑜𝜕/𝜕𝑥)̇𝑦→(−𝑖/𝑘𝑜𝜕/𝜕𝑦)}exp𝑖𝜀𝑘_𝑜𝐿𝑥̇,𝑦,̇𝑥,𝑦;𝜁≅const^′⋅exp𝑖𝜀𝑘𝐻_𝑜^⌢−𝑖𝑘_𝑜𝜕𝑖𝜕𝑥,−𝑘_𝑜𝜕,𝜕𝑦,𝑥,𝑦;𝜁 (48) where𝐇^⌢−𝑖𝑘_𝑜𝜕𝑖𝜕𝑥,−𝑘_𝑜𝜕𝜕𝑦,𝑥,𝑦;𝜁=𝐿_{̇𝑥→(−𝑖/𝑘𝑜𝜕/𝜕𝑥)̇𝑦→(−𝑖/𝑘𝑜𝜕/𝜕𝑦)}𝐿−𝑖𝑥̇,𝑦,̇𝑥,𝑦;𝜁=𝑥̇𝑘_𝑜𝜕−𝑖𝜕𝑥+𝑦̇𝑘_𝑜𝜕.𝜕𝑦−𝐿𝑥̇,𝑦,̇𝑥,𝑦;𝜁 (49)

The approximation of (48) makes it possible to simplify the right-hand side of (47) and to write it in the differential form𝜕𝜕𝜁+𝑖𝑘𝐻_𝑜^⌢−𝑖𝑘_𝑜𝜕𝑖𝜕𝑥,−𝑘_𝑜𝜕𝑈𝜕𝑦,𝑥,𝑦;𝜁𝑥,𝑦;𝜁=0. (50)

Comparing (49) and the optical Hamiltonian (25) [1, 8, 15, 16], we can see that the operator 𝐻^⌢((−𝑖/𝑘_𝑜)(𝜕/𝜕𝑥),(−𝑖/𝑘_𝑜)(𝜕/𝜕𝑦),𝑥,𝑦;𝜁) is obtained from the optical Hamiltonian by replacing coordinates 𝑥, 𝑦 and momenta 𝑝 and 𝑞 with the corresponding operators in the position representation [5, 25]𝑥^⌢=𝑥,𝑦^⌢=𝑦,^⌢𝑖𝑝=−𝑘_𝑜𝜕,𝜕𝑥^⌢𝑖𝑞=−𝑘_𝑜𝜕,𝜕𝑦 (51) that is, 𝐻^⌢−𝑖𝑘_𝑜𝜕𝑖𝜕𝑥,−𝑘_𝑜𝜕𝜕𝑦,𝑥,𝑦;𝜁=−𝑛²+1𝑥,𝑦;𝜁𝑘_𝑜²𝜕²𝜕𝑥²+1𝑘_𝑜²𝜕²𝜕𝑦². (52) This makes it possible to interpret (50) as the symbolic representation of the Helmholtz equation [5, 15, 16]𝜕²𝑈𝜕𝑥²+𝜕²𝑈𝜕𝑦²+𝜕²𝑈𝜕𝜁²+𝑘²_𝑜𝑛²𝑥,𝑦;𝜁𝑈=0. (53)

So we substantiate the “quantum theory of light rays” (see Appendix C) [15, 16]: the operator representation of the Helmholtz equation can be obtained from the optical Hamiltonian by replacing the momenta 𝑝 and 𝑞 with the corresponding momentum operators in the position representation.

Note that the symbolic representation of the Helmholtz equation (50) has the complex conjugation analog𝜕𝜕𝜁+𝑖𝑘𝐻_𝑜^⌢∗−𝑖𝑘_𝑜𝜕𝑖𝜕𝑥,−𝑘_𝑜𝜕𝑈𝜕𝑦,𝑥,𝑦;𝜁^∗𝑥,𝑦;𝜁=0, (54) where the asterisk ∗ denotes complex conjugation.

It is shown that formal quantization, that is, replacing coordinates 𝑥, 𝑦 and momenta 𝑝 and 𝑞 with the corresponding operators, (52), can be substantiated by using the relation between the Fourier and the Legendre transformations.

7. Position Representation and Momentum Representation of Operators

It is known that a wave field of radiation can be described in the position and the momentum representations which are connected to each other by Fourier transformations [25].

The fundamental symmetry of the operator theory is the duality principle: the 𝑘_𝑜-Fourier transformation (13a) connects the complex amplitude 𝑈(𝑥,𝑦;𝑧) of the wave field with its spectrum 𝑈(𝑝,𝑞;𝑧), and also replaces the operators in position representation (51) with the operators in the momentum representation [25]:=𝑖𝑥̆𝑘_𝑜𝜕𝑖𝜕𝑝,𝑦̆=𝑘_𝑜𝜕,̆𝜕𝑞𝑝=𝑝,𝑞̆=𝑞. (55)

More precisely,^𝑘𝑜𝐹_{𝑥→𝑝𝑦→𝑞}=𝑖𝑥𝑈𝑥,𝑦;𝜁𝑘_𝑜𝜕𝜕𝑝^𝑘𝑜𝐹_{𝑥→𝑝𝑦→𝑞}𝑈𝑥,𝑦;𝜁_𝑘𝑜,𝐹_{𝑥→𝑝𝑦→𝑞}=𝑖𝑦𝑈𝑥,𝑦;𝜁𝑘_𝑜𝜕𝜕𝑞^𝑘𝑜𝐹_{𝑥→𝑝𝑦→𝑞}𝑈𝑥,𝑦;𝜁_𝑘𝑜,𝐹_{𝑥→𝑝𝑦→𝑞}−𝑖𝑘_𝑜𝜕𝑈𝜕𝑥𝑥,𝑦;𝜁=𝑝⋅^𝑘𝑜𝐹_{𝑥→𝑝𝑦→𝑞}𝑈𝑥,𝑦;𝜁_𝑘𝑜,𝐹_{𝑥→𝑝𝑦→𝑞}−𝑖𝑘_𝑜𝜕𝑈𝜕𝑦𝑥,𝑦;𝜁=𝑞⋅^𝑘𝑜𝐹_{𝑥→𝑝𝑦→𝑞}𝑈𝑥,𝑦;𝜁. (56)

Therefore the 𝑘_𝑜-Fourier transformation (13a) replaces the Hamilton operator in the position representation (52) with the Hamilton operator in the momentum representation:𝐹𝐻^𝑘𝑜_{𝑥→𝑝𝑦→𝑞}^⌢−𝑖𝑘_𝑜𝜕𝑖𝜕𝑥,−𝑘_𝑜𝜕𝑈=𝐇̆𝑖𝜕𝑦;𝑥,𝑦;𝜁𝑥,𝑦;𝜁𝑝,𝑞;𝑘_𝑜𝜕,𝑖𝜕𝑝𝑘_𝑜𝜕𝜕𝑞;𝜁^𝑘𝑜𝐹_{𝑥→𝑝𝑦→𝑞}𝑈𝑥,𝑦;𝜁. (57a) That is true for the complex conjugation Hamilton operator as well:𝐹𝐇^𝑘𝑜_{𝑥→𝑝𝑦→𝑞}^⌢∗−𝑖𝑘_𝑜𝜕𝑖𝜕𝑥,−𝑘_𝑜𝜕𝑈=𝐇̆𝜕𝑦;𝑥,𝑦;𝜁𝑥,𝑦;𝜁_∗𝑖−𝑝,−𝑞;−𝑘_𝑜𝜕𝑖𝜕𝑝,−𝑘_𝑜𝜕×𝜕𝑞;𝜁_𝑘𝑜𝐹_{𝑥→𝑝𝑦→𝑞}𝑈𝑥,𝑦;𝜁. (57b)

Thus, the Hamilton operator in the momentum representation𝐻̆𝑝,𝑞;𝑖𝑘_𝑜𝜕,𝑖𝜕𝑝𝑘_𝑜𝜕𝜕𝑞;𝜁=−𝑛²𝑖𝑘_𝑜𝜕,𝑖𝜕𝑝𝑘_𝑜𝜕𝜕𝑞;𝜁−𝑝²−𝑞² (58) can be derived from the optical Hamiltonian (25) by replacing coordinates 𝑥, 𝑦 and momenta 𝑝 and 𝑞 with the corresponding operators in the momentum representation (55). If the Helmholtz equation in the position representation has form (50), then in the momentum representation it takes the following form:𝜕𝜕𝜁+𝑖𝑘_𝑜𝐇̆𝑖𝑝,𝑞;𝑘_𝑜𝜕,𝑖𝜕𝑝𝑘_𝑜𝜕𝑈𝜕𝑞;𝜁𝑝,𝑞;𝜁=0. (59)

8. Transport Equation for the Wigner Distribution Function

It is known that in a plane 𝑧=𝜁 the spatial coherence of a scalar quasi-monochromatic wave is characterized by two-point correlation functions—the mutual intensityΓ𝑥₁,𝑦₁,𝑥₂,𝑦₂;𝜁≡𝑈𝑥₁,𝑦₁;𝜁𝑈^∗𝑥₂,𝑦₂;𝜁 (60a) or the mutual coherence spectrumΓ𝑝₁,𝑞₁,𝑝₂,𝑞₂;𝜁≡𝑈𝑝₁,𝑞₁;𝜁𝑈^∗𝑝₂,𝑞₂;𝜁. (60b)

Here, the angle brackets denote averaging over an ensemble.

By introducing mean variables 𝑥=(𝑥₁+𝑥₂)/2, 𝑦=(𝑦₁+𝑦₂)/2, 𝑝=(𝑝₁+𝑝₂)/2, 𝑞=(𝑞₁+𝑞₂)/2 and difference variables Δ𝑥=𝑥₂−𝑥₁, Δ𝑦=𝑦₂−𝑦₁, Δ𝑝=𝑝₂−𝑝₁, Δ𝑞=𝑞₂−𝑞₁ and using the Fourier transform from difference variables, it is easy to obtain from the correlation functions (60a) and (60b) the mathematically equivalent but more convenient expressions of the spatial coherence of radiation—the WDF [17, 26–34]:𝑊=𝑥,𝑦;𝑝,𝑞;𝜁_𝑘𝑜𝐹_{Δ𝑥→𝑝Δ𝑦→𝑞}𝑈𝑥+Δ𝑥2,𝑦+Δ𝑦2;𝜁×𝑈^∗𝑥−Δ𝑥2,𝑦−Δ𝑦2;𝜁 (61a) or𝑊=𝑥,𝑦;𝑝,𝑞;𝜁_𝑘𝑜𝐹⁻¹_{Δ𝑝→𝑥Δ𝑞→𝑦}𝑈𝑝+Δ𝑝2,𝑞+Δ𝑞2×𝑈;𝜁_∗𝑝−Δ𝑝2,𝑞−Δ𝑞2;𝜁. (61b)

Thus, the WDF describes the partially coherent radiation field in the position and in the momentum representations simultaneously.

The fundamental equations of optics are the equations describing propagation of light through an optically inhomogeneous medium. They are Hamilton’s equations in geometrical optics, (31), and the Helmholtz equation in wave optics, (53). Propagation of the WDF through an optically inhomogeneous (and random) media is described by the so-called transport equation [32, 33]. Let us show that the duality between the position and momentum representations of the Helmholtz equation, (50) and (59), allows us to prove it.

If the WDF (61a) is differentiated with respect to the longitudinal variable 𝜁, then we have𝜕𝑊𝑥,𝑦;𝑝,𝑞;𝜁=𝜕𝜁_𝑘𝑜𝐹_{Δ𝑥→𝑝Δ𝑦→𝑞}𝜕𝑈𝑥+Δ𝑥/2,𝑦+Δ𝑦/2;𝜁𝜕𝜁×𝑈^∗𝑥−Δ𝑥2,𝑦−Δ𝑦2;𝜁+𝑈𝑥+Δ𝑥2,𝑦+Δ𝑦2×;𝜁𝜕𝑈^∗𝑥−Δ𝑥/2,𝑦−Δ𝑦/2;𝜁.𝜕𝜁 (62) Using the Helmholtz equation in the position representation (50) and its complex conjugate analogue (54) we can rewrite (62) in the following form:𝜕𝑊𝑝,𝑞;𝑥,𝑦;𝜁𝜕𝜁=−𝑖^𝑘𝑜𝐹_{Δ𝑥→𝑝Δ𝑦→𝑞}𝑘𝐇_𝑜^⌢−𝑖𝑘_𝑜𝜕𝜕(),−𝑖𝑥+Δ𝑥/2𝑘_𝑜𝜕𝜕𝑦+Δ𝑦/2;𝑥+Δ𝑥2,𝑦+Δ𝑦2−;𝜁_⌢𝐇^∗−𝑖𝑘_𝑜𝜕𝜕(),−𝑖𝑥−Δ𝑥/2𝑘_𝑜𝜕𝜕;𝑦−Δ𝑦/2𝑥−Δ𝑥2,𝑦−Δ𝑦2×𝑈;𝜁𝑥+Δ𝑥2,𝑦+Δ𝑦2;𝑧𝑈^∗×𝑥−Δ𝑥2,𝑦−Δ𝑦2;𝜁. (63)

According to Appendix D𝜕𝜕()=1𝑥±Δ𝑥/22𝜕±𝜕𝜕𝑥,𝜕𝜕Δ𝑥𝜕=1𝑦±Δ𝑦/22𝜕±𝜕𝜕𝑦,𝜕Δ𝑦 (64) so (65) takes the following form:𝜕𝑊𝑝,𝑞;𝑥,𝑦;𝜁𝜕𝜁=−𝑖^𝑘𝑜𝐹_{Δ𝑥→𝑝Δ𝑦→𝑞}𝑘𝐇_𝑜^⌢−𝑖𝑘_𝑜12𝜕+𝜕𝜕𝑥,−𝑖𝜕Δ𝑥𝑘_𝑜12𝜕+𝜕𝜕𝑦;𝜕Δ𝑦𝑥+Δ𝑥2,𝑦+Δ𝑦2−;𝜁_⌢𝐇^∗−𝑖𝑘_𝑜12𝜕−𝜕𝜕𝑥,−𝑖𝜕Δ𝑥𝑘_𝑜12𝜕−𝜕𝜕𝑦;𝜕Δ𝑦𝑥−Δ𝑥2,𝑦−Δ𝑦2𝑈;𝜁𝑥+Δ𝑥2,𝑦+Δ𝑦2;𝜁𝑈^∗×𝑥−Δ𝑥2,𝑦−Δ𝑦2. (65)

Using (57a), (57b) and (65) can be written as𝜕𝑊𝑝,𝑞;𝑥,𝑦;𝑧𝜕𝜁=−𝑖𝑘_𝑜𝐇̆𝑖𝑝−2𝑘_𝑜𝜕𝑖𝜕𝑥,𝑞−2𝑘_𝑜𝜕𝑖𝜕𝑦;𝑥+2𝑘_𝑜𝜕,𝑖𝜕𝑝𝑦+2𝑘_𝑜𝜕−𝐇̆𝜕𝑞;𝜁_∗𝑖𝑝−2𝑘_𝑜𝜕𝑖𝜕𝑥,𝑞−2𝑘_𝑜𝜕;𝑖𝜕𝑦𝑥+2𝑘_𝑜𝜕𝑖𝜕𝑝,𝑦+2𝑘_𝑜𝜕.𝜕𝑞;𝜁×𝑊𝑝,𝑞;𝑥,𝑦;𝜁 (66)

Using the identity−𝑖𝐇̆−𝐇̆^∗=2Im𝐇̆, (67) where the symbol Im𝐇̆ denotes the imaginary part of the Hamilton operator 𝐇̆, (66) can be rewritten in the symmetric form [32, 33]:𝜕𝑊𝑝,𝑞;𝑥,𝑦;𝑧𝜕𝜁=2𝑘_𝑜𝐇̆𝑖Im𝑝−2𝑘_𝑜𝜕𝑖𝜕𝑥,𝑞−2𝑘_𝑜𝜕𝑖𝜕𝑦;𝑥+2𝑘_𝑜𝜕,𝑖𝜕𝑝𝑦+2𝑘_o𝜕𝑊.𝜕𝑞;𝜁𝑝,𝑞;𝑥,𝑦;𝜁 (68)

We can approximate the Hamiltonian in (68) via a Taylor series expansion by the Appendix B𝐇̆𝑝−𝑖2𝑘_𝑜𝜕𝑖𝜕𝑥,𝑞−2𝑘_𝑜𝜕𝑖𝜕𝑦;𝑥+2𝑘_𝑜𝜕𝑖𝜕𝑝,𝑦+2𝑘_𝑜𝜕𝑖𝜕𝑞;𝜁=𝐻𝑝,𝑞;𝑥,𝑦;𝜁⋅exp2𝜕𝐻𝜕𝜕𝑥+𝜕𝑝𝜕𝐻𝜕𝜕𝑦−𝜕𝑞𝜕𝐻𝜕𝜕𝑝−𝜕𝑥𝜕𝐻𝜕𝜕𝑞.𝜕𝑦 (69) In that case the transport equation (68) takes the form𝜕𝑊𝑝,𝑞;𝑥,𝑦;𝜁𝐻𝑖𝜕𝜁=−2Im𝑝,𝑞;𝑥,𝑦;𝜁⋅exp2𝜕𝐻𝜕𝜕𝑥+𝜕𝑝𝜕𝐻𝜕𝜕𝑦−𝜕𝑞𝜕𝐻𝜕𝜕𝑝−𝜕𝑥𝜕𝐻𝜕𝜕𝑞𝜕𝑦×𝑊𝑝,𝑞;𝑥,𝑦;𝜁 (70) or, in the symbolic notation of Besieris and Tappert [34],𝜕𝑊𝑝,𝑞;𝑥,𝑦;𝜁𝐻𝑖𝜕𝜁=−2Im𝑝,𝑞;𝑥,𝑦;𝜁⋅exp2^⟵𝜕𝜕⃗𝜕𝑥+𝜕𝑝^⟵𝜕𝜕⃗𝜕𝑦−𝜕𝑞^⟵𝜕𝜕⃗𝜕𝑝−𝜕𝑥^⟵𝜕𝜕⃗𝜕𝑞.𝜕𝑦×𝑊𝑝,𝑞;𝑥,𝑦;𝜁 (71)

In this expression the arrows over the partial derivatives indicate the direction in which the derivatives act, that is, those with a left arrow act on 𝐻(𝑝,𝑞;𝑥,𝑦;𝜁), and those with a right arrow act on 𝑊(𝑝,𝑞;𝑥,𝑦;𝜁).

In a weakly inhomogeneous medium the transport equation (71) reduces to𝜕𝑊𝑝,𝑞;𝑥,𝑦;𝜁𝐻×𝑖𝜕𝜁=−2Im𝑝,𝑞;𝑥,𝑦;𝜁1+2^⟵𝜕𝜕⃗𝜕𝑥+𝜕𝑝^⟵𝜕𝜕⃗𝜕𝑦−𝜕𝑞^⟵𝜕𝜕⃗𝜕𝑝−𝜕𝑥^⟵𝜕𝜕⃗𝜕𝑞.𝜕𝑦×𝑊𝑝,𝑞;𝑥,𝑦;𝜁 (72) Thus, in the usual notation, this equation takes the form of the Liouville equation [27–31]𝜕𝑊()[],𝜕𝜁=−2Im𝐻⋅𝑊−Re𝐻,𝑊 (73) where Re denotes the real part and[Re𝐻,𝑊]=𝜕Re𝐻𝜕𝑥𝜕𝑊+𝜕𝑝𝜕Re𝐻𝜕𝑦𝜕𝑊−𝜕𝑞𝜕Re𝐻𝜕𝑝𝜕𝑊−𝜕𝑥𝜕Re𝐻𝜕𝑞𝜕𝑊𝜕𝑦 (74) is the “Poisson bracket”. The Liouville equation (73) reduces to the ordinary differential equation [27–31]𝑑𝑊𝑑𝜁−2Im𝐻^⌢⋅𝑊=0. (75)

In the special case that 𝐻(𝑝,𝑞;𝑥,𝑦;𝜁) is a real function of 𝑥,𝑦,𝑝, and 𝑞, Equation (75) can be formulated in geometric optical terms as follows: along a light ray the WDF has a constant value:𝑑𝑊𝑑𝜁=0. (76)

It is known that the WDF with some natural assumptions can be interpreted as a wave analogue of the radiance [1]. So equation (76) is the theorem about invariance of radiance along the light ray, that is, the most fundamental theorem of radiometry. Note that this theorem in the form of the Liouville equation (73) can be derived from the point of view of Hamilton’s optics, too [20].

9. Conclusion

We have unite all known quantum-mechanical approaches to wave and statistical optics in a single theory using Feynman’s path integral as a fundamental principle. In particular, we have substantiated the quantization which allows one to translate the Hamiltonian of geometrical optics into the Helmholtz equation of wave optics.

Appendices

A. The Stationary-Phase Approximation

Wave optics often takes account of oscillatory integrals of the type𝐼=_𝑅2𝑔𝑥,𝑦exp𝑖𝑘_𝑜Ψ𝑥,𝑦𝑑𝑥𝑑𝑦 (A.1) in which 𝑔(𝑥,𝑦) is a slowly varying complex-valued function and Ψ(𝑥,𝑦) is a slowly varying real-valued function of the real variables (𝑥,𝑦). We assume that the phase function Ψ(𝑥,𝑦) has one stationary point (𝑥_𝑜,𝑦_𝑜) in the interior of the region of integration:𝜕Ψ𝑥,𝑦𝜕𝑥^𝑥=𝑥_{𝑜𝑦=𝑦𝑜}=0,𝜕Ψ𝑥,𝑦𝜕𝑦^𝑥=𝑥_{𝑜𝑦=𝑦𝑜}=0. (A.2) The stationary point (𝑥_𝑜,𝑦_𝑜) is nondegenerateΨdet^{′′𝑥𝑥}𝑥_𝑜,𝑦_𝑜Ψ^′′_𝑥𝑦𝑥_𝑜,𝑦_𝑜Ψ^′′_𝑦𝑥𝑥_𝑜,𝑦_𝑜Ψ^′′_𝑦𝑦𝑥_𝑜,𝑦_𝑜≠0. (A.3) Then, according to the stationary phase method [5, 12–14, 21, 22], the integral (A.1) can be approximated by the equation_𝑅2𝑔𝑥,𝑦exp𝑖𝑘_𝑜Ψ≈𝑥,𝑦𝑑𝑥𝑑𝑦_{𝑘𝑜→∞}2𝜋𝑖𝑘_𝑜⋅Ψdet^{′′𝑥𝑥}𝑥_𝑜,𝑦_𝑜Ψ^′′_𝑥𝑦𝑥_𝑜,𝑦_𝑜Ψ^′′_𝑦𝑥𝑥_𝑜,𝑦_𝑜Ψ^′′_𝑦𝑦𝑥_𝑜,𝑦_𝑜^−1/2⋅𝑔𝑥_𝑜,𝑥_𝑜exp𝑖𝑘_𝑜Ψ𝑥_𝑜,𝑦_𝑜. (A.4)

B. The Taylor Series in Operator Form

Using the Taylor series, a function 𝑈(𝑥−Δ𝑥,𝑦−Δ𝑦) can be expanded in a power series up to the linear term:𝑈𝑥−Δ𝑥,𝑦−Δ𝑦;𝜁≅𝑈𝑥,𝑦;𝜁−𝜕𝑈𝜕𝑥Δ𝑥−𝜕𝑈𝜕𝑦Δ𝑦. (B.1) We can rewrite (B.1) in the operator form [5, 6]𝑈𝑥−Δ𝑥,𝑦−Δ𝑦;𝜁≅𝑈𝑥,𝑦;𝜁−𝜕𝑈𝜕𝑥Δ𝑥−𝜕𝑈=𝜕𝜕𝑦Δ𝑦1−𝜕𝜕𝑥Δ𝑥−𝑈.𝜕𝑦Δ𝑦𝑥,𝑦;𝜁 (B.2)

Using the Taylor series, we can expand the exponent exp(−(𝜕𝑈/𝜕𝑥)Δ𝑥,−(𝜕𝑈/𝜕𝑦)Δ𝑦) in a power series up to the linear term:−𝜕exp𝜕𝜕𝑥Δ𝑥−𝜕𝜕𝑦Δ𝑦≅1−𝜕𝜕𝑥Δ𝑥−𝜕𝑦Δ𝑦. (B.3) Comparing (B.2) with (B.3), we see that the right-hand sides of these equations are equal, so we can write the approximate operator equation [5, 35]:𝑈−𝜕𝑥−Δ𝑥,𝑦−Δ𝑦;𝜁=exp𝜕𝜕𝑥Δ𝑥−𝑈.𝜕𝑦Δ𝑦𝑥,𝑦;𝜁 (B.4)

C. Formal Quantum Theory of Light Rays

In the quantum theory of light rays, the coordinates 𝑥, 𝑦 retain their meaning as numbers, but the canonically conjugate variables 𝑝, 𝑞 (the momenta), (52), and the Hamiltonian 𝐻, (51), become differential operators𝜆𝑝=−𝑖^𝑜𝜕2𝜋,𝜆𝜕𝑥𝑞=−𝑖^o𝜕2𝜋𝜆𝜕𝑦,𝐻=𝑖^o𝜕2𝜋,𝜕𝑧 (C.1a)(C.1b)(C.1c)

where 𝜆_𝑜 is the free-space wavelength. Thus, the operator equation of the quantum theory of light rays can be obtained from the optical Hamiltonian (25) and operator equation (C.1c) by replacing momenta 𝑝 and 𝑞 with the momentum operators (C.1a) and (C.1b)−𝑛²+𝜆𝑥,𝑦;𝑧^𝑜2𝜋²𝜕²𝜕𝑥²+𝜆_𝑜2𝜋²𝜕²𝜕𝑦²𝜆=𝑖^𝑜𝜕2𝜋𝜕𝑧. (C.2)

A wave function 𝑈(𝑥,𝑦;𝑧) of the quantum theory of light rays is the usual distribution of the scalar wave field over the plane 𝑧. Applying the operator equation (C.2) to the wave function 𝑈(𝑥,𝑦;𝑧) results in the equation𝑖𝜕𝑈𝜕𝑧=−2𝜋𝜆_𝑜²𝑛²+𝜕𝑥,𝑦;𝑧²𝜕𝑥²+𝜕²𝜕𝑦²⋅𝑈 (C.3) that takes into account the symmetry of the optical medium relative to the 0𝑧 axis. Equation (C.3) can be interpreted as the symbolic form of the Helmholtz equation:𝜕²𝑈𝜕𝑥²+𝜕²𝑈𝜕𝑦²+𝜕²𝑈𝜕𝑧²+2𝜋𝜆_𝑜²𝑛²𝑥,𝑦;𝑧𝑈=0. (C.4)

D.

One has(a)𝜕𝑓𝜕(𝑥+Δ𝑥/2)=lim_𝜀→0𝑓(𝑥+Δ𝑥/2+𝜀)−𝑓(𝑥+Δ𝑥/2)𝜀=lim_𝜀→0𝑓𝜀𝑥+2+(Δ𝑥+𝜀)/2−𝑓(𝑥+Δ𝑥/2)𝜀=lim_𝜀→0𝑓(𝑥+𝜀/2+(Δ𝑥+𝜀)/2)−𝑓(𝑥+(Δ𝑥+𝜀)/2)𝜀+lim_𝜀→0𝑓(𝑥+(Δ𝑥+𝜀)/2)−(𝑥+Δ𝑥/2)𝜀=12𝜕𝑓+𝜕𝑥𝜕𝑓.𝜕Δ𝑥 (D.1)

One has the following:(b) 𝜕𝑓𝜕(𝑥−Δ𝑥/2)=lim_𝜀→0𝑓(𝑥−Δ𝑥/2+𝜀)−𝑓(𝑥−Δ𝑥/2)𝜀=lim_𝜀→0𝑓(𝑥+𝜀/2−(Δ𝑥−𝜀)/2)−𝑓(𝑥−Δ𝑥/2)𝜀=lim_𝜀→0𝑓(𝑥+𝜀/2−(Δ𝑥−𝜀)/2)−𝑓(𝑥−(Δ𝑥−𝜀)/2)𝜀+lim_𝜀→0𝑓(𝑥−(Δ𝑥−𝜀)/2)−(𝑥−Δ𝑥/2)𝜀=12𝜕𝑓−𝜕𝑥𝜕𝑓.𝜕Δ𝑥 (D.2)

Acknowledgments

It is pleasure to thank Dr. Wilhelm Becker and Dr. Mikhail Kalashnikov for helpful discussions and conversations.