The Fourier U(2) Group and Separation of Discrete Variables

The linear canonical transformations of geometric optics on two-dimensional screens form the group $Sp(4,R)$, whose maximal compact subgroup is the Fourier group $U(2)_F$; this includes isotropic and anisotropic Fourier transforms, screen rotations and gyrations in the phase space of ray positions and optical momenta. Deforming classical optics into a Hamiltonian system whose positions and momenta range over a finite set of values, leads us to the finite oscillator model, which is ruled by the Lie algebra $so(4)$. Two distinct subalgebra chains are used to model arrays of $N^2$ points placed along Cartesian or polar (radius and angle) coordinates, thus realizing one case of separation in two discrete coordinates. The $N^2$-vectors in this space are digital (pixellated) images on either of these two grids, related by a unitary transformation. Here we examine the unitary action of the analogue Fourier group on such images, whose rotations are particularly visible.


Introduction
The real symplectic group Sp(4, ) of linear canonical transformations is widely used in paraxial geometric optics [1, Part 3] on two-dimensional screens, and in classical mechanics with quadratic potentials. Its maximal compact subgroup is the Fourier group U(2) F = U(1) F ⊗ SU(2) F [2,3]. These are classical Hamiltonian systems with two coordinates of position (q x , q y ) ∈ 2 , and two coordinates of momentum (p x , p y ) ∈ 2 , which form the phase space 4 , and which are subject to linear canonical transformations that preserve the volume element. The central subgroup U(1) F contains the fractional isotropic Fourier transforms, which rotate the planes (q x , p x ) and (q y , p y ) by a common angle. The complement group SU(2) F contains anisotropic Fourier transforms that rotate the planes (q x , p x ) and (q y , p y ) by opposite angles; it also contains joint rotations of the (q x , q y ) and (p x , p y ) planes; and thirdly, gyrations, which 'cross-rotate' the planes (q x , p y ) and (q y , p x ). The rest of the Sp(4, ) transformations shear or squeeze phase space, as free flights and lenses, or harmonic oscillator potential jolts.
This classical system can be quantizedà la Schrödinger into paraxial wave or quantum models. Indeed, the group of canonical integral transforms was investigated early by Moshinsky and Quesne in quantum mechanics [4,5], [6,Part 4] and by Collins in optics [7], and yielded the integral transform kernel that unitarily represents the (two-fold cover of the) group Sp(4, ) Figure 1. Action of the Fourier group on the sphere (from [21]).
In two dimensions (q x , q y ) ∈ 2 and with the corresponding momenta (p x , p y ) ∈ 2 , one can build ten independent quadratic functions that will close into the real symplectic Lie algebra sp (4, ). Four linear combinations of them generate transformations that include the fractional Fourier transform (FT), which are of interest in optical image processing, isotropic FT F 0 := 1 4 p 2 x + p 2 y + q 2 x + q 2 y =: 1 2 (h x + h y ), (2.3) anisotropic FT F 1 := 1 4 p 2 x − p 2 y + q 2 x − q 2 y =: 1 2 (h x − h y ), (2.4) gyration F 2 := 1 2 (p x p y + q x q y ), (2.5) rotation F 3 := 1 2 (q x p y − q y p x ). (2.6) Their Poisson brackets close within the set, {F 1 , F 2 } P = F 3 , and cyclically, {F 0 , F i } P = 0, i = 1, 2, 3, and characterize the Fourier Lie algebra u(2) F = u(1) F ⊕ su(2) F , with F 0 central. This is the maximal compact subalgebra of (2.2), u(2) F ⊂ sp(4, ) [3]. In Fig. 1 we show the transformations generated by (2.4)-(2.6) on the linear space of the su(2) algebra, which leaves invariant the spheres F 2 1 + F 2 2 + F 2 3 = F 2 0 . The Lie algebra u(2) and group U(2) will appear in several guises, so it is important to distinguish their 'physical' meaning. The classical Fourier group U(2) F of linear optics will be matched by a 'Fourier-Kravchuk' group U(2) K acting on the position space of finite 'sensor' arrays in the finite oscillator model. In one dimension, this model is based on the algebra u(2), whose generators are position, momentum, and energy; while in two, the algebra is su(2) x ⊕ su(2) y = so (4). In all, the well-known properties of su(2) = so(3) will of use [22].

Finite quantization in one dimension
The usual Schrödinger quantization (q →Q = q·, p →P = −iηd/dq) of the phase space coordinates leads to the paraxial model of wave optics (for the reduced wavelength η = λ/2π) or to oscillator quantum mechanics (for η = ). In one dimension, (2.2) yields the three generators of the (double cover of the) group Sp(2, ) of canonical integral transforms [4] acting on functionscontinuous infinite signals -in the Hilbert space L 2 ( ). In two dimensions, the Schrödinger quantization of (2.2) yields the generators of Sp(4, ), and its Fourier subgroup is generated by the quadratic Schrödinger operatorsF 0 ,F 1 ,F 2 ,F 3 corresponding to (2.3)-(2.6), that we indicate with bars, and are up-to-second order differential operators acting on the Hilbert space. Finite quantization on the other hand, asigns self-adjoint N × N matrices q → Q, p → P to the coordinates of phase space.
Since the oscillator algebra is noncompact, it cannot have a faithful representation by finite self-adjoint matrices [12]; we must thus deform the four-generator osc into a compact algebra, of which there is a single choice: u(2). We should keep the geometry and dynamics of the harmonic oscillator contained in the two Hamilton equations in (2.1), with commutators i[ · , · ] in place of the Poisson brackets {·, ·} P , and in place of the oscillator Hamiltonian h, a matrix K added with some constant times the unit matrix. We call K the pseudo-Hamiltonian. The third -and fundamental -commutator [Q, P ] will then give back K. In [13] we proposed the assignments of self-adjoint N × N matrices given by the well-known irreducible representations of the Lie algebra su(2) of spin j, thus of dimension N = 2j + 1, where the matrix representing position is chosen diagonal. Using the generic notation {J 1 , J 2 , J 3 } for generators of su(2), the matrices are The commutation relations of these matrices are i.e., [J 1 , J 2 ] = −iJ 3 , and cyclically. (3.4) The central generator 1 of osc is corresponded with the unit N × N matrix that generates a U(1) central group of multiplication by overall phases. This completes the algebra u(2) realized by matrices acting on the linear space C N of complex N -vectors. The spectra Σ of the matrices Q of position {q}, P of momentum { }, and K of pseudoenergy {κ} are thus discrete and finite, with j integer or half-integer. We define the mode number as n := κ + j ∈ {0, 1, . . . , 2j}, and energy is n + 1 2 . The Lie group generated by {1; Q, P, K} is the U(2) group of linear unitary transformations of C N . Finally, the sum of squares of these matrices on C N is the Casimir invariant C := Q 2 + P 2 + K 2 = j(j + 1)1.
In Dirac notation, C N is a spin-j u(2) representation, where the eigenstates of position q and of mode n = κ + j satisfy The position eigenvectors |q 1 form a Kronecker basis for C N where the components of |n 3 will be provided, from (3.3), by a three-term relation, i.e., a difference equation that is satisfied by the Wigner little-d functions [22] for κ = n − j. (The reader may be disconcerted for having J 1 diagonal, whereas almost universally J 3 is declared the diagonal matrix; formulas obtained by permuting 1 → 2 → 3 → 1 are unchanged.) The overlaps between the two bases in (3.5) are the finite oscillator wavefunctions, This overlap is expressed here [23] in terms of the square root of a binomial coefficient in q, which is a discrete version of the Gaussian, and a symmetric Kravchuk polynomial of degree n, K n (q) ≡ K n (q; 1 2 , N − 1) = K q (n − j) [24]. We have called the Ψ n (q)'s Kravchuk functions [13]; they form a real, orthonormal and complete basis for C N . The Lie exponential of the self-adjoint matrix K generates the unitary U(1) group of fractional Fourier-Kravchuk transforms, whose realization we shall indicate by U(1) K .
For future use we give the general expression of the Wigner little-d functions for angles β ∈ [0, π] as a trigonometric polynomial [12,22] where due to the denominator factorials, the summation extends over the integer range of max(0, m − m ) ≤ k ≤ min(j − m , j + m); for m − m > 0 the reflection formulas (3.9) apply.

Two-dimensional systems
The two-dimensional classical oscillator algebra osc 2 is sometimes taken to consist of the Poisson brackets between 1, q x , q y , p x , p y , and h x , h y ; and sometimes an angular momentum m = q 1 p 2 − q 2 p 1 is included, Poisson-commuting with the total Hamiltonian h = h x + h y . The finite quantization of two-dimensional systems deforms osc 2 to the Lie algebra su(2) ⊕ su(2) = so(4) in both cases. (The subalgebra of 1 need not concern us here.) To work comfortably with the six generators of so(4), we consider the customary realization of J i,i = −J i ,i by self-adjoint operators that generate rotations in the (i, i ) planes, which obey the commutation relations that we summarize with the following pattern: A generator J i,i has non-zero commutator with all those in its row i and its column i (reflected across the i = i line); all its other commutators are zero. The asignments of position, momentum and pseudo-energy with the su

The Cartesian coordinate system
Passing from one to two dimensions can be achieved prima facie by building the direct sum algebra su(2) x ⊕ su(2) y , which is accidentally equal to so(4). This isomorphism is shown in patterns by where the square super-index of all generators, X 0 i , indicates the identification between the so(4) generators J i,i with the Cartesian coordinates and observables.
Since the x-generators commute with the y-generators in (4.2), using (3.5) a Cartesian basis of positions and also a basis of modes can be simply defined as direct products |q x , q y 0 1 := |q x 1x |q y 1y and |n x , n y 0 where the pseudo-energy eigenvalues are n i − j = κ i | j −j . The two-dimensional finite oscillator wavefunctions in Cartesian coordinates are thus given as a product of two Ψ n (q)'s from (3.6), (3.7), Ψ 0 nx,ny (q x , q y ) = 0 1 q x , q y |n x , n y 0 on positions q x , q y | j −j and of mode numbers n x , n y | 2j 0 . The positions can be accomodated in the square pattern of Fig. 2a, and the modes in the rhombus pattern of Fig. 2b. The Cartesian eigenstates of the finite oscillator are shown in Fig. 3.
Since K 0 x and K 0 y generate independent rotations in the (Q x , P x ) and in the (Q y , P y ) planes respectively, their sum K := K 0 x + K 0 y = J 1,2 transforms phase space isotropically. Thus we identify K ∈ so(4) as the generator of a group U(1) K of isotropic fractional Fourier-Kravchuk transforms in C N × C N = C N 2 , and corresponding with the operator 2F 0 ∈ U(2) F in (2.3). That sum commutes with A := K 0 x − K 0 y = J 3,4 ∈ so(4), which generates skew-symmetric Fourier rotations by opposite angles in the x-and y-phase planes; so A corresponds with 2F 1 ∈ u(2) F in (2.4). However, a counterpart of the 2F 3 generator of rotations in the (q x , q y ) and (p x , p y ) planes cannot be found within so(4), and neither can the gyrations in (2.5). These will be imported in the next section.

The polar coordinate system
The six generators of the Lie algebra so(4) can be identified with positions and modes following an asignment different from the Cartesian direct sum of the previous section. We indicate the  new generators with a circle super-index, X • i ; as in the classical case, a generator of rotations M between the x-and y-axes should satisfy the commutation relations while the isotropic Fourier generator K = J 1,2 should rotate between position and momentum operators, The commutator (4.10) asserts that K and M can be used to define a basis for C N 2 with the quantum numbers of mode and angular momentum, that will be given below. Expressed in the pattern (4.1), a new so(4) generator assignment that fulfills these requirements is This assignment satisfies the conditions (4.6)-(4.10), but implies the further commutators Of these, (4.12) echoes the su(2) nonstandard commutator in (3.4), while the commutator (4.13) is also nonstandard, and indicates that Q • x and Q • y cannot be simultaneously diagonalized. First, we find operators with quantum numbers corresponding to radius and angle; we use the subalgebra chain When both su(2)'s in (4.2) have the same Casimir eigenvalue j(j + 1), the principal Casimir invariant is i<i J 2 i,i = 2j(j + 1)1. In so(4) there is a second invariant, is identically zero in these 'square' cases [22]. Let us now consider the Casimir operator of the subalgebra so(3) ⊂ so(4) in (4.14), which is The Gel'fand-Tsetlin branching rules [25] determine that C N 2 then decomposes into subspaces C ρ that are irreducible under this so (3), where (4.15) exhibits the eigenvalues [22] ρ(ρ + 1), Although R is not an element of the algebra so (4), we shall identify it as the radius operator. The final link in the reduction (4.14) is M , whose eigenvalues in the ρ-representation of so (3)  (2ρ + 1) = (2j + 1) 2 = N 2 , the same as for the N × N square grid, and shown in Fig. 4a. Thus we define the radius-angular momentum (ra) eigenvectors and we adopt ρ as the radial position coordinate.
In the last step we use the discrete Fourier transform matrix, F ρ in each (2ρ + 1)-dimensional subspace, to pass between angular momenta and angles, and thus build the Kronecker basis of states |ρ, φ localized at a definite radius ρ and the 2ρ + 1 equidistant angles φ k as |ρ, φ k : where the ψ ρ 's are fixed but arbitrary phases. In Fig. 4b we show the resulting arrangement of N 2 points (ρ, φ k ) thus defined. Note that the Fourier transformation in (4.17) is linear and unitary but is not an element of the group SO (4): it has been imported to act on each of the 2ρ + 1-dimensional so(3) irreducible representation subspaces ρ ∈ {0, 1, . . . , 2j}.
Having the position Kronecker eigenstate basis |ρ, φ k for polar coordinates, we now define the basis of mode and angular momentum ma, eigenbasis of the commuting operators K = K 0 x + K 0 y with total mode number n| 4j 0 (pseudo-energy κ = n − 2j), and of M = K 0 x + K 0 y with angular momentum m, placed in the pattern (4.11). We build these eigenstates |n, m • MA asking for We observe that whereas the ra states in (4.16) are classified by the Gel'fand-Tsetlin chain of subalgebras so(4) ⊃ so(3) ⊃ so (2), the ma states in (4.18) follow the chain so(4) = su(2)⊕su(2) ⊃ u(1) ⊕ u(1), which in ordinary quantum theory entails the coupling of two spin-j representations to total spin ρ. The overlaps between the ma and ra states should therefore be Clebsch-Gordan coefficients ∼ C j , m 1 , j , m 2 , j m , with m 1 = 1 2 (κ + m) and m 2 = − 1 2 (κ − m), adding to the total m [22]. In the present construction though, the subalgebra chain (4.14) reduces along the 'lower' subalgebras, and this differs from the original Gel'fand-Tsetlin reduction that reduces along the upper ones; also, generally one counts su(2) multiplet states from the top down, and we have ordered them from the bottom up. The overlap of ra and ma states entails an extra phase that must be computed carefully. It is [18,19] • RA ρ, m|κ + j, m • MA = ϕ(j, ρ, κ, m) C j ,  In Fig. 5 we show the C N 2 basis of ma states. Note that due to the Clebsch-Gordan selection rules, states of a given angular momentum m are nonzero only at radii ρ ≥ |m|. The generator K in (4.11) generates rotations between both position and momentum operators, and corresponds to twice the isotropic Fourier transform generator 2F 0 in (2.3). The action of K(ω) := exp(−iωK) on the ma basis of C N 2 is thus Under these rotations, images f (ρ, φ k ) = ρ, φ k | f on the polar-pixellated screen will thus transform into where for each radius ρ| 2j 0 there is a (2ρ + 1) × (2ρ + 1) matrix representing the same rotation . (4.23) These are circulating matrices, functions of φ k − φ k = 2π(k − k )/(2ρ + 1) modulo 2π, and periodic in k, k modulo 2ρ + 1. For each radius, the 'unit' rotation angle is θ = 2π/(2ρ + 1), and for multiples l thereof, the matrix (4.23) is nonzero at the diagonal k = k + l. In Fig. 6 we give an example of such rotation. Isotropic Fourier transformations and rotations on the polar screen are produced by generators within the so(4) algebra in the pattern (4.11). However, the pattern also informs us that with linear combinations of K and M , we can not gyrate the planes (Q • x , P • y ) and (Q • y , P • x ) jointly as with 2F 2 in (2.5); also missing is the anisotropic Fourier transform generated by 2F 1 in (2.4), which was natural in the Cartesian basis. These transformations will be imported on C N 2 in Section 7.

Importation of rotations on the Cartesian screen
We noted above that in the subalgebra chain (4.14), the generators of isotropic Fourier transformations and rotations, K ↔ 2F 0 and M ↔ 2F 3 are domestic to so(4), while anisotropic Fourier transformations and gyrations, corresponding to 2F 2 , are foreign. Now, in the Cartesian subalgebra decomposition of so(4) in (4.2), the two independent Fourier transform generators (2.3), (2.4), K 0 x ↔F 0 +F 1 and K 0 y ↔F 0 −F 1 are domestic to so(4), while those of gyrations and rotations, 2F 2 and 2F 3 , are foreign. Since we cannot complete a fully domestic U(2) K Fourier-Kravchuk group in correspondence with the Fourier group U(2) F ⊂ Sp(4, ), we must import the missing transformations. Such importation was used in (4.17) with the (2ρ + 1) × (2ρ + 1) discrete Fourier transform matrix.
The domestic and imported transformations properly mesh, and we see that the Fourier-Kravchuk group U(2) K is indeed represented unitarily and faithfully by (6.4) on C N 2 . Its action on images f (q x , q y ) pixellated on the Cartesian screen can be found from here, writing Ω := (φ, θ, ψ), through a real similarity transformation by the matrix formed with the C N 2 basis (4.5) of Cartesian mode Kravchuk functions, f ω,Ω (q x , q y ) := D(ω, Ω) f (q x , q y ) = q x ,q y D 0 (q x , q y ; q x , q y ; ω, Ω) f (q x , q y ), (6.5) where, with n = n x + n y = n x + n y and n| 4j 0 the kernel is D 0 (q x , q y ; q x , q y ; ω, Ω) = nx,ny n x n y Ψ 0 nx,ny (q x , q y ) e −i(n−2j)ω D n/2 (nx−ny)/2,(n x −n y )/2 (Ω)Ψ 0 n x ,n y (q x , q y ), (6.6) and these N 2 × N 2 matrices are unitary representations of U(2) on the Cartesian grid of points in Fig. 2b. With this result we now turn to the polar screen.

The Fourier group on polar screens
In the subalgebra chain of so(4) that produces the polar screen (4.11), the commuting rotations and isotropic Fourier-Kravchuk transforms are domestic. To complete the Fourier-Kravchuk group U(2) K on the polar screen, we must import either gyrations or the anisotropic transform. Since these transformations have been realized already in the Cartesian screen basis, we should find a unitary map between both screens. This transformation has been studied in [18,19], and consists in identifying the eigenstates of mode and angular momentum in both bases; in the polar basis these are |ρ, m • RA in (4.16) shown in Fig. 5. In the Cartesian basis such states will be constructed now as eigenvectors of R(θ) = exp(−iθM ) with eigenvalues e −iθm , or equivalently ofF 3 with eigenvalues 1 2 m, corresponding to the eigenvalues of 1 2 M in (2.6). These Cartesian states will be linear combinations, respecting total mode number n = n x + n y , of all states From (5.1), the coefficients obey the difference equation n y (n x + 1)C n,m nx+1,ny−1 − 1 2 m C n,m nx,ny + n x (n y + 1)C n,m nx−1,ny+1 = 0, which is the ubiquitous su(2) three-term recursion relation of the Wigner little-d functions for angle 1 2 π (around the 1-axis [22]), now with j ↔ 1 2 n = 1 2 (n x + n y ) and 1 2 m ↔ 1 2 (n x − n y ). We thus define ma states having angular momentum m on the square grid as  Fig. 4]. The basis of states (7.2) is shown in Fig. 8, which can now be identified with the basis -also of mode and angular momentum -in Fig. 5. We can now import the equivalence between the ma bases |n, m 0 A pixellated image on the Cartesian screen, f 0 (q x , q y ) ∈ C N 2 , can be thus unitarily transformed into an image on the polar screen, f • (ρ, φ k ) ∈ C N 2 through the transformation where the transform kernel, using (7.3), is The transformation of images inverse to (7.4), from the polar to the Cartesian screen, is We note that these transformations map real functions onto real functions. In Fig. 9 we show an example of the map (7.4) on the 0-and-1 image of the letter R. The action of the Fourier-Kravchuk group U(2) K on images on the polar screen can be now derived from its action (6.5), (6.6) on the Cartesian screen: where the kernel is a representation of U(2) K on the polar screen We have thus far not found a more compact analytic form for the representation of the Fourier group on the discrete coordinates of radius and angle.
The N 2 × N 2 matrix U = (V ) † that is needed above, is arduous to calculate but its numerical values need be computed only once for each size N of the screen, and will serve to transform any image from the Cartesian to the polar screen and back. For N = 32 there are ≈ 10 6 values to be stored, and the transformation of an image as that in Fig. 9 between the screens involves this number of sums and products. Commercial image-rotating algorithms store the original image and interpolate the few Cartesian pixels around each geometrically determined point. Certainly, the Fourier-Kravchuk group of transformations of images presented here does not provide a fast algorithm, but it is unitary, and hence does not loose information, as interpolation invariably does; it is exact.

Concluding remarks
Willard Miller Jr. is recognized as having initiated the modern study of systems whose governing equation allows for separation in various coordinate systems using Lie algebraic reduction methods [26]. These equations are differential because space has been considered continuous. Finite spaces with discrete coordinates, used as foundation for finite Hamiltonian systems characterized by governing difference equations, have not been considered -to our knowledge -in the same context. The present paper reviews our work in the particular case of a two-dimensional Hamiltonian system whose two governing Hamilton equations identify as a finite harmonic oscillator, mothered by the compact algebra so (4). And since it happens that algebra has two distinct subalgebra reductions, we can relate them with Cartesian an polar coordinates.
The difference equations whose solutions are the finite oscillator wavefunctions and the Clebsch-Gordan coefficients are contained in so(3) and so(4), and entail special Kravchuk and Hahn polynomials. Other finite polynomials, such as those of Meixner and Pollaczek, also appear in related harmonic oscillator models [27], and thus could be interpreted in terms of two-dimensional physical and optical models. Further, finite three-dimensional coordinate systems beyond Cartesian and circular cylindric ones may be of interest in applications. But we should end with a word of caution regarding applications: the two-dimensional arrays that purportedly measure the quality of Laguerre-Gauss (polar) laser beams generally use sampled Laguerre-Gauss functions; the question of whether these, or the mode-angular momentum finite functions, are the most efficient to find the mode and angular momenta of the actual beams, seems to favor the sampled functions [28,29], although they do not form orthonormal bases for the space of pixellated images. But only the discrete function bases synthesize all images exactly.
The continuous two-dimensional harmonic oscillator is a superintegrable system, which is known to separate generally in elliptic coordinates, of which the Cartesian and polar are limits. We have been unable so far to find a corresponding finite array that would allow for pixellation following ellipses and hyperbolas, although there are three enticing leads: diagonalization of the so(3) operator J 2 + αJ 3 [26]; a line of functions with a parameter which joins continuously Kronecker deltas with Clebsch-Gordan coefficients [30]; and the production by holographic means of laser beams whose nodes follow elliptic coordinates [31]. There are also 'geometrical' reasons to doubt that this is possible, however: whereas in Cartesian coordinates all pixel sizes are equal, and in polar ones they are almost so (except very near to the center), trying to visualize these two pixellations as limits of a finite elliptic pixellation presents some problems that we have not yet been able to overcome [32].