Traditional projection and displays are no longer sufficient to meet the increasingly diverse demands, such as holographic head-up display (HUD) [1] and augmented reality (AR) [2]. Advanced imaging and display technologies, represented by Fourier optics with diffractive optics [3] as a prime example, provide an alternative approach. For instance, the advent of the computer-generated hologram (CGH) has provided various manifestations for holography, including three-dimensional (3D) reconstruction [4–8], high-resolution capabilities [5,7,9,10], and implementations in augmented reality (AR) and virtual reality (VR) [7,11,12]. In recent years, driven by increasingly mature design and process capabilities, higher precision and flexibility in optimizing and fabricating devices have been gained. This has unlocked boundless possibilities in manufacturing, imaging, and display, such as the utilization of freeform surfaces [13–15], diffractive optical elements [13,16,17], optical waveguides [2,11], and liquid crystal [18,19]. In particular, the advent of metasurfaces has opened versatile modalities for coherent imaging [20–23]. Beyond the ultra-thin merits inherent to traditional refractive devices, metasurfaces offer a substantial advantage in polarization control [24]. Metasurface holography, with its high information density in a planar format, holds tremendous potential for cutting-edge displays [23–28], paving the way for the realization of high-fidelity holography.

Despite continuous advancements in the design and fabrication methods of optical elements, without exception, these designs involving scalar diffractive optics invariably employ conventional Fourier methods. The angular spectrum method (ASM), single Fresnel transform (SFT), and their variants represent the prevailing ones. Despite their widespread utilization and computational efficiency with fast Fourier transform (FFT), they suffer inherent limitations. Fixed sampling intervals, dictated by the conservation of the space-bandwidth product (SBP), restrict their universality. ASM requires equal-sized object and image planes, while SFT correlates image plane size linearly with propagation distance under the operation of optical Fourier transformation [3,29,30]. In a considerable proportion of cases, ASM based or SFT based diffraction iteration processes may result in insufficient casting or resource waste. Challenges arise when intricate scale transformation operations are required such as AR and holographic HUD, or when different depths of 3D projection demand varying levels of detail. Moreover, in color imaging and display, conventional Fourier methods such as the Fraunhofer method (FM) may exhibit chromatic aberrations, necessitating complex correction and demanding additional computational resources. So, the pixel pre-scaling for different colors during the computational process will compromise fidelity. In addition, the widely acknowledged accurate Rayleigh–Sommerfeld integral (RSI) is computationally slow and thus less suitable for inverse design applications. Conventional Fourier methods offer a rapid and efficient tool for diffraction calculation, and associated sampling and calculation strategies have been continually proposed and refined [4,31–37]. But in general they did not break inherent scaling limitations and merely represented a slight extension of conventional ASM and SFT. At times, new strategies were accompanied by decreased image quality and narrowed application scope. For example, diffraction calculation algorithms can be based on non-uniform FFT [38]. However, they may not be well-suited for inverse design applications, so they are rarely used. Overall, the modulation of the light field still has a relatively low degree-of-freedom, even in relatively simple cases of coherent imaging and display, such as CGH imaging and automatic HUD. The fundamental reason for this lies in the fact that these methods are entrapped by the two local solutions of the Helmholtz equation in wave optics. New solutions are difficult to achieve through fast Fourier algorithms.

Finding concise solutions to address these longstanding issues directly from the Helmholtz equation may be a futile endeavor at present. Therefore, we seek breakthroughs from a higher-dimensional perspective to derive new diffraction computation methods suitable for fast Fourier algorithms. Phase space analysis, initially introduced in quantum mechanics through the behavior of the Wigner distribution function (WDF) [39], has found extensive applications in describing optical systems [40,41]. It has been instrumental in understanding various phenomena, including the relationship between coherence and radiation measurement [42], the connection between ambiguity function and holography [43], and the measurement of partial coherence [44]. The completeness of the WDF for signal description fully reveals the spatial distribution of light in phase space. In higher dimensions, the WDF based phase space analysis itself constitutes a form of Fourier analysis. Phase space analysis directly analyzes the modulation of the light field from the perspective of distribution and transmission. Conventional Fourier optics can be viewed as a degenerate form of phase space analysis, and they are not mutually exclusive. In recent research [45–47], the relationship between conventional Fourier methods and WDF propagation was discussed. However, there are few reports on FFT based Fourier algorithms designed for coherent light field diffraction calculations using phase space analysis.

Within the higher-dimensional phase space perspective, this paper explores a new course for coherent imaging computations. First, it explores a universal framework based on matrix decomposition that empowers the researchers to design diffraction calculation methods according to different application scenarios. Subsequently, a variable-scale model is delivered that allows the maximum spatial frequency of the image plane to be freely selected within certain constraints. This model still leverages Fourier analysis for fast computations and effectively addresses the size limitation issue inherent in conventional Fourier methods. The comparison and advantages of our model against other canon methods are illustrated in Table 1. Experimental validations, including advanced holographic and tomographic displays, demonstrate the model’s robustness in variable-scale capability and chromatic aberration correction. Additionally, our model is applied to the design of full-color, near-zero crosstalk holographic metasurfaces, showcasing their applicability at nanoscale. Therefore, the longstanding algorithmic chromatic aberration issue in full-color holography would be effectively addressed. Rooted in phase space analysis, our approach offers a new tool for coherent imaging and display, providing an effective diffraction computation scheme for the diffractive optics community.Table 1.

Comparison of Our Model against the Canon Methods

Under coherent illumination, a two-dimensional (2D) complex field distribution Uo(r) possesses a WDF characterized in phase space. The WDF of Uo(r) on the object or input plane is a four-dimensional (4D) function, encompassing information from both spatial and Fourier domains simultaneously, providing a comprehensive description of the light field signal: Wo(ro,ko)=∬−∞∞Uo(ro+ro′/2)Uo*(ro−ro′/2)⁢exp(−ikot·ro′)dro′,(1)where * denotes the conjugation, ro′=(Δxo,Δyo)t is the spatial offset relative to r in spatial domain, and ko=(kox,koy)t is the k-vector in the Fourier domain with t indicating the transpose operation. The point spread function (PSF) in free space characterizes the system’s response to a point source, while the transfer function analyzed in the Fourier domain represents the response to plane waves. In phase space, since the WDF is the complete signal representation, its propagation can respond to both the point source and a single frequency, simultaneously and mathematically. It not only seamlessly combines the spatial domain and Fourier domain but also establishes a strong connection with the concept of “light rays” in geometrical optics. The corresponding ray spread function, which manifests as the response of the product of two Dirac delta functions Wo(r,k)=δ(r−ro,k−ko), forms a double Wigner distribution.

In scenarios where the Hamiltonian is limited to quadratic terms at most, the analysis of wave propagation can be elegantly achieved through straightforward matrix transformations. Linear canonical transform (LCT) emerges as a robust instrument for characterizing the behavior of light during its propagation within first-order optical systems [48]. By employing the paraxial approximation, the trajectories of light rays within such optical systems can be efficiently described using a fundamental ABCD matrix, which is a specific form of the LCTs. The ray transformation matrix T4D=[A,B;C,D] in phase space is applied as [ri;ki]=T4D[ro;ko]=[A,B;C,D][ro;ko], and the propagation relationship between Wi on the image plane and Wo on the object plane is expressed as [40,49]Wi(Ar+Bk,Cr+Dk)=Wo(r,k).(2)

In addition, the WDF transport equation of coherent light in free space is (k/k)∂W∂r+((k2−|k|2)/k)∂W∂z=0,(3)where z is the propagation distance, λ is the wavelength, and k=2π/λ denotes the wave number. Equation (3) has the following solution: Wi(r,k,z)=Wo(r−kz/(k2−|k|2),k,0).(4)

Considering the low-band-limited characteristics of general objects and the applicability of sampling theorem, adapting |k|≪k and the transformation of f=k/2π, T4D has the following form: T4D≅[I2πzk−1I0I],(5)where I is the identity matrix. T4D contains the propagation characteristics of Wo to Wi, having the ability to analyze the complete propagation pattern. For convenience, we adopt its 2D form, which can be effortlessly extended to four dimensions: T2D≅[1λz01].(6)T2D is a concise Fresnel transform matrix, which is a two-order LCT. This straightforward matrix encapsulates extensive Fourier analysis features [50]. In Ref. [51], it has been proven that certain optical transformations, such as the Fourier transform, Fresnel transform, quadratic phase modulation represented by single-lens modulation, and fractional Fourier transform, all have their corresponding matrix forms in phase space. Moreover, these matrices can be decomposed into LCT matrix cascades based on the conservation of SBP for numerical tracking. The LCTs corresponding to various transformations or modulations applied to Uo(r) in phase space are detailed in Appendix A. In phase space analysis, to evaluate the accuracy of a diffraction calculation method in a single numerical diffraction calculation or tasks such as phase retrieval in inverse design, it is essential to provide an appropriate matrix cascade while ensuring the conservation of the SBP to correspond to the respective computational operations. Meanwhile, the LCT matrix cascades, which are employed for free-space diffraction, can be designed, and they presumably have the capability to degenerate into T4D under specific conditions. In essence, this leads to a universal framework where we can extend an infinite number of methods to calculate light propagation, catering to diverse needs, especially in free space diffraction within first-order optical systems. That is to say, in the phase space context, researchers or engineers can design diffraction computation methods independently and provide corresponding matrix cascades which must closely align with T4D or T2D. Alternatively, T4D or T2D can be directly decomposed into a few operations like Fourier transforms.

To validate our perspective and break the size limitations introduced by ASM and SFT, a decomposition approach with variable-scale characteristics would be proposed for the calculation of coherent imaging and display. To facilitate the operations of FFT in either spatial domain or Fourier domain, T2D can be decomposed as follows: T2D=Q[m−1mλz]M[m]F−π/2[1]Q[−λzm]Fπ/2[1]Q[1−mλz],(7)where Q[h]≅limϵ→0[1,ϵ;h,1+ϵh]=[1,0;h,1] and h is a non-zero constant. M[m]=[m,0;0,1/m] is the scaling matrix, with m being the scaling factor, i.e., a magnifier. Fπ/2[h]=[0,h;−1/h,0] and F−π/2[h]=[0,−h;1/h,0] are the Fourier transform matrix and the inverse Fourier transform matrix, respectively. The phase space diagram (PSD) is defined as the region within phase space where WDF is non-negative. In Fig. 1(a), a representative PSD is depicted, characterized by a spatial extent Lo and a bandwidth Bo. Although it can be approximated through methodologies such as signal superposition [46], finding its physical entity is challenging. In this context, our focus is solely on illustrating the transformation trajectory of it. Figure 1(a) illustrates the transformation process of the PSD. It clearly shows the process of clockwise (CW) rotation, counterclockwise (CCW) rotation, and magnifying and shearing of the PSD under the applications of each LCT according to Eq. (7). The same processes of ASM and SFT are easy to be delivered (see Appendix A). The decomposition in Eq. (7) transforms a simple single phase space transformation into a series of chirp modulation and Fourier transformation operations, which is efficient and computationally tractable for implementation on a computer. Therefore, in spatial domain, the diffraction process in free space can be described as Ui(xi,m)=CQp3(xi,m)M{[Uo(xo)Qp1(xo,m)]⊗Qp2(xo,m)},(8)where xo and xi correspond to the coordinate systems on the object plane and the image plane, respectively. C is a complex constant, and ⊗ denotes linear convolution. The chirp modulator Qp1(xo,m)≅exp[iφp1(xo,m)]=exp[iπ(1−m)xo2/λz] is equivalent to Q[(1−m)/λz] in phase space, and similarly, Qp3(xi,m) is equivalent to Q[(m−1)/mλz]. In addition, Qp2(xo,m) is the expression of Q[−λz/m] in the spatial domain. The transformation operator M{·} corresponds to M[m]. Although M{·} does not result in any computational operation during numerical computation, in physical terms it directly maps the field distribution of xo system to xi system. This results in the calculated field distribution possessing variable-scale characteristics. In the case of m=1, Eq. (8) collapses into a typical Fresnel form. Moreover, it shares similarities with Schmidt’s description [32]. In the context of a fixed object plane with fixed sampling procedures, we retain the freedom to control the spatial extent of the image plane, i.e., the tunability of spatial frequency, while conserving the SBP. We refer to this method, allowing precise governance over the spatial frequency, with the spatial frequency tunable method (SFTM). It is a tool with potential capabilities rooted in phase space analysis, applicable for solving problems in specific realms such as beam shaping, holography, metasurface design, and pixel mismatch in end-to-end optimization [52]. It provides a powerful and easily implementable computational tool for the diffractive optics community, enabling them to break free from the confines of traditional methods. Figure 1(b) demonstrates the variable-scale capability and the automatic aberration correction capability of the SFTM.

The relationships and constraints (CSTs) of m, z, and several other quantities are illustrated in Table 2 (the sampling criteria analysis for complex amplitude and intensity is detailed in Appendix B). CST 4 offers two distinct sampling approaches, sampling in the Fourier domain and sampling in the spatial domain, to accommodate various scale transformation requirements. The partitioning of the sampling region, introduced by CST 4, is referred to as the spatial frequency sampling region (SFSR) and spatial sampling region (SSR), respectively. Figure 2(a) illustrates the permissible region of allowed m and z values for λ=0.532μm, N0=N/2=1000, and δxo=8μm scenarios. Meanwhile, the restrictions of both ASM and SFT have been delineated, where they represent only a segment of a straight line in the m–z space. The threshold z0=N0(δxo)2/λ is the boundary between ASM and SFT. It is evident that, with proper design, values for m can be quite flexible, and it goes far beyond the conventional Fourier analysis represented by ASM, SFT, and their derivatives.Table 2.

Sampling CSTs for the SFTM

Figures 2(b) and 2(c) show the comparison of amplitude between the SFTM based numerical results and the analytical solutions for a 4mm square aperture illuminated by a converging spherical wave when δxo=8μm. The focus of the spherical wave is at the same distance from the square aperture as the diffraction distance z. When z is 200 mm and the scaling factor m is 0.8, the peak signal-to-noise ratio (PSNR) of the numerical result is 57.23 dB. When z is 800 mm and m is 5.0, the PSNR is 63.67 dB. Additionally, when δxo=2μm and the aperture width is 2 mm, the PSNR of the numerical result is 42.34 dB at 100 mm (m=0.6) and is 47.11 dB at 800 mm (m=2.5). The PSNR values are all greater than 40 dB, proving that the SFTM is quite robust.

We employ a reflective 2D phase-only spatial light modulator (SLM) with a pixel size of 8 μm to verify the accuracy and practicality of the SFTM. The inverse diffraction method has been employed to design 1000×1000-pixel CGHs. A slight perturbation from the correct calculation may result in the “butterfly effect” in inverse diffraction after interactions [53]. Although the singularity of inverse diffraction kernels can be disregarded in the design of CGHs for display under homogeneous light illumination, it is nearly impossible to accurately simulate and experiment with the correct image if the diffraction calculation model deviates from T4D. Traditional terminology refers to this situation as non-convergence. We utilize a three-stage iterative Fourier transform algorithm (IFTA) based on adaptive constraints in the Fourier domain (the algorithm flowchart is shown in Appendix C) instead of the traditional monotonous Gerchberg–Saxton (GS) algorithm to design CGHs with a signal window (SW) of 750×750 pixels. This choice is motivated partly by the high sensitivity of the GS algorithm to local optima. On the other hand, practically, the phase modulated by the SLM is quantized into 256 steps, and achieving continuous control of phase is often impractical in the manufacturing process of diffractive optical elements. Importantly, the three-stage IFTA demonstrates superior reconstruction accuracy compared to the traditional GS algorithm, exhibiting higher values of the structural similarity index measure (SSIM) and lower values of root-mean-square error (RMSE) after hundreds of iterations [54].

Figure 12 in Appendix F illustrates the experimental setup for the monochromatic display of the CGHs. Figure 3(a) shows a comparison of the images projected by CGHs designed with ASM and SFTM within a propagation distance shorter than the threshold z0=120.3mm. The actual size of the image is directly proportional to the scaling factor m. When m=1, the size of the image within the SW is 6 mm. Since m was set to 1 and can be seen within the oversampled region of transfer function for ASM at z=40mm, 50 mm, 60 mm, 80 mm, and 100 mm [specific parameter selections have been marked with solid points in Fig. 2(a)], the rabbits designed with ASM and SFTM exhibit almost identical effects, without apparent aliasing or twin images. Similarly, as indicated by the solid points in the SFT regime in Fig. 2(a), we compared SFT and SFTM at the propagation distances of z>z0; see Fig. 3(b). Letting SFTM have the same linear scaling characteristics as SFT, we captured the enlarged images at z=160mm, 200mm, 240 mm, 280 mm, and 320 mm, respectively. In Fig. 3(b), the remarkable congruence among the cat images, acquired through three-stage IFTA using SFT and SFTM, is evident. This substantiates that the inverse algorithm, governed by the SFTM, not only upholds m=1 for z<z0, where SFTM serves as a high-order alternative to ASM, but also showcases identical linear scaling characteristics to SFT for z>z0. Consequently, the SFTM can seamlessly substitute SFT in this regime. Additionally, holographic image projection experiments were conducted across different magnification ranges within the permissible regions of the m–z space. As depicted in Fig. 3(c), a 0.8× chick and a 1.5× chick were captured at z=90mm and z=110mm, respectively. Likewise, 0.8×, 0.6×, and 2.0× peacocks were projected at z=180mm, 210 mm, and 260 mm.

In reality, most constraints of the sampling criteria originate from the Nyquist sampling theorem [3], which can be slightly relaxed in practical operations. On the other hand, T4D allows the size of the support on the image plane extending to Lo+λzBo. By employing suitable filtering techniques, we can try to slightly break through the limitations of the allowed m–z space for the SFTM, leading to further breakaway from the limitations of conventional Fourier methods. Similarly, in Fig. 3(c), a 0.6× chick out of the allowed m–z space is projected at z=70mm, and at another location of z=300mm, a 3.5× peacock, just beyond the boundary of the constraints, is projected. They both yield excellent imaging results. Furthermore, the SFTM algorithm’s capability for long-distance projection is assessed by enlarging an image of a shark to 20 times its original size at z=1800mm, as depicted in Fig. 3(d). In conventional Fresnel or Fraunhofer models, the achievable magnification factor is typically limited to around 15.

Conventionally, the far-field projection of holographic images relies on SFT based IFTA or FM based IFTA. However, chromatic aberration exists, wherein the size of the image plane is directly proportional to the wavelength as depicted in Fig. 4(a). Holographic image distortion within a large field of view (FOV) was usually addressed through two methodologies: image pre-processing and hologram correction. Nevertheless, both are computationally intensive. They also do not fundamentally alter the intrinsic linear-scale characteristics. Importantly, fidelity will be compromised. As previously highlighted, our model is applicable for long-distance holographic image projection. By employing the SFTM based three-stage IFTA for holographic design, not only can the crosstalk be further minimized, but also the linear-scale mismatch introduced by Fresnel-FFT algorithms can be rectified. This approach eliminates the need for intricate and constrained pre-processing procedures. Moreover, the scale of the image plane is adjustable, offering substantial flexibility for the realization of full-color holographic displays. We decoupled a color windmill image based on the three primary colors and let them iterate within their respective color channels. In Fig. 4(d), at a propagation distance of z=400mm, a 3.5× magnification of a full-color holographic image projection is obtained with a PSNR of 21.71 dB. Its optical reconstruction result is shown in Fig. 4(e). Notably, this process was conducted without pre-processing or chromatic aberration correction. The SFT based numerical and optical reconstruction results are presented in Figs. 4(b) and 4(c), respectively. The numerical reconstruction exhibits a lower PSNR, and the quality of the optical reconstruction is inferior to the results obtained with the SFTM.

During the general 3D reconstruction process, different depths of a 3D object may require varying detail levels. Typically, SFT is used for reconstructing large 3D objects, obtaining 3D Fresnel holograms, as the image plane size of SFT linearly increases with the propagation distance. However, the size of an image plane at a particular depth remains fixed, while the 2D image at that depth may vary depending on the object and its position. Pre-scaling of the 2D target images at each depth is generally required for 3D reconstruction. This inevitably leads to fidelity reduction in the entire 3D scene. Recently, holograms of large objects with up to hundreds of planes were reported [4], offering an effective solution for future dynamic large-field 3D holography. However, this remarkable work is still based on the Fresnel method, requiring pre-processing of the 2D target image at each depth, which not only compromises fidelity but also consumes excessive computational resources. If our model can be utilized in tomography involving two or more image planes, significant improvements on addressing these issues can be expected.

In the previous section, the scaling capabilities and chromatic aberration correction abilities of the SFTM were validated. Therefore, we conducted experiments to evaluate the multi-plane imaging capabilities of the SFTM algorithm. Figure 5(a) depicts a schematic illustration of the SFTM based tomographic reconstruction. In this representation, a fixed-size pixel on the hologram plane can be mapped to different image planes, each with pixels of varying sizes. This emphasizes the variable-scale capability of the SFTM in tomographic reconstruction. Utilizing the tomography algorithm (see Appendix C for details), we imaged three pentagrams of different internal structures with magnification of 3×, 6×, and 10× at depths z=200mm, 400 mm, and 700 mm, respectively. As shown in Fig. 5(b), a comparison of three identical duck objects reveals significant variations in magnification among the distinct pentagrams. For 3D Fresnel holograms, the depth of focus (DOF) of each image plane is a crucial factor influencing 3D reconstruction [4]. It is directly affected by the relation DOFi∝λ(zi/NSWδxo)2, where i is the serial number of the image plane. The requirement zi+1−zi=γ(DOFi+DOFi+1) for low crosstalk on image planes imposes a significant constraint on the 3D reconstruction of 3D Fresnel holograms, where γ is an empirical parameter. To mitigate this limitation, that is, to decrease DOFs, one relatively straightforward approach, without considering additional optimization or rearrangement of Fresnel zone plate phases, is to increase the resolution of the hologram. However, it is often impractical and can lead to resource consumption and waste. We note that, under the same projection depth, the increase in magnification leads to an elongation of the DOF, as is commonly known in the case of Fresnel holograms. Figures 5(c) and 5(d) show the results of two-layer full-color tomography. The rainbow flower is positioned at z=150mm with a magnification of 1.5, and the cube is at z=300mm, enlarged by m=2.0. It can be seen that when z is appropriate and the magnification is not excessive, image crosstalk from the other plane is not significant. But in Fig. 5(b), we can see that the 6× pentagram on the middle image plane is impacted by crosstalk from the 10× pentagram on the farthest image plane, whereas the 3× pentagram on the nearest image plane experiences minimal interference. The settings for these magnifications of the pentagrams serve two purposes: first, to verify whether the SFTM possesses sufficient diffraction modulation capability, and second, for the convenience of capture. In numerous scenarios, achieving such extreme magnification may not be necessary. Therefore, it becomes viable to project images with varying spatial frequencies onto distant and multiple image planes. This approach can effectively overcome the constraints of 3D Fresnel holograms, paving ways for the design of high-quality, super-multi-plane 3D variable-scale holograms.

The superiority of metasurfaces over traditional refractive and diffractive devices is evident not only in their ultra-thin features but also in their powerful capability for controlling multiple wavelengths and facilitating efficient polarization reuse. In our previous work [24], we implemented a TiO2 metasurface design characterized by low crosstalk and featuring a tri-polarization-channel configuration for holographic display. Like most designs of full-color holographic metasurfaces, the previous design relied on the conventional Fraunhofer model, necessitating intricate pixel correction operations for different color channels. It will inevitably result in reduced resolution and fidelity. If the variable-scale model can be employed in the context of metasurface holography to validate its applicability at nanoscale, this problem will be effectively solved. Furthermore, researchers would have access to a more powerful tool for diffraction calculations in metasurface design, freeing themselves from the constraints of ASM and SFT.

Figure 6(a) illustrates the schematic of the near-zero crosstalk metasurface holography using the SFTM based algorithm. The letters “H,” “N,” and “U” corresponding to different colors should all have the same size. Figure 6(b) illustrates the TiO2 meta-atom design of the metasurface. The nanopillars, which can produce unique phases on mutually orthogonal linear polarizations, are fabricated on a square-shaped SiO2 substrate with a period of P=400nm. They possess three tunable freedoms, i.e., parameters D1, D2 and rotation angle θ, with a fixed height of H=800nm. The degree-of-freedom of the Jones matrix can be dynamically controlled by the three parameters of the meta-atom [24], achieving excellent broadband transmission and conversion efficiency characteristics of the metasurface. We decouple the light of the three colors, red (λ=0.633μm), green (λ=0.532μm), and blue (λ=0.450μm), into different polarization channels, minimizing crosstalk between various polarization states. Then, on the image plane, they will be superimposed based on the principles of the three primary colors. Leveraging the broadband transmission characteristics of the TiO2 meta-atom, the polarization states of incident green and red lights are set to x-polarized and blue light is set to y-polarized. The output red light is set to x-polarized, while green and blue lights are set to y-polarized. We decouple the target color image into red, green, and blue channels. Each channel component is allowed to undergo the three-stage IFTA process. The SFTM model is used for the reconstruction of the phase-only holograms (see Appendix D). The size of the metasurface is 400μm×400μm, encompassing N02=1000×1000 units of meta-atoms. For the incident light of red, green, and blue, the allowed m–z spaces are identified in Fig. 8. Significantly, when designing the metasurface, special consideration must be directed towards the allowed m–z space for blue light incidence. This is particularly crucial since blue light exhibits the weakest diffraction capability, as constrained by the size of Lo+λbzBo. Leveraging the lookup algorithm, we systematically identify and intricately arrange the appropriate meta-atoms with fixed parameters D1, D2, and θ. This meticulous arrangement allows precise control over the Jones matrix, ensuring compliance with the specifications of three channels, diverse polarization states, and distinct phase distributions. As a result, this systematic procedure leads to the assembly of the metasurface.

The finite-difference time-domain (FDTD) method is employed to calculate the electromagnetic field for meta-atoms. Subsequently, the metasurface design is translated into sample through the electron beam lithography (EBL) followed by the reactive ion etching (RIE) during the sample fabrication procedure (see Appendix E for details). We designed a holographic metasurface with m=1.2 and z=426.6μm. On the image plane, the same-size letters “H,” “N,” and “U” were projected, corresponding to the red, green, and blue colors, respectively. Figures 6(c) and 6(d) show the SEM images of the designed metasurface in oblique-view and top-view. As shown in Figs. 6(e) and 6(f), the simulation and experimental results closely match, with negligible crosstalk between the color channels (the characterization setup for the metasurface sample can be found in Appendix F). The fabricated metasurface samples may exhibit minor deviation from the design. These subtle differences lead to an actual wavefront passing through a metasurface not entirely the same as the designed one. Therefore, the experimental results shown in Fig. 6(f) exhibit a slight blurring phenomenon.

The observable low crosstalk superposition capability among the three primary colors highlights the effectiveness of the meta-holography system. Our demonstration underscores the applicability of the SFTM in nanoscale holography, effectively resolving the chromatic aberration issue caused by diffraction algorithms, thus paving the way for diverse applications in advanced coherent imaging and display.

The matrix cascade designed in our work provides powerful computational applications and extremely high-degree-of-freedom for diffraction calculations. It is evident that both the liquid crystal array in the SLM and the meta-atom in the holography metasurface provide phase modulation capabilities spanning nearly 0 to 2π. Consequently, the theoretical bandwidth of the WDF in phase space can be broadened. Conventional IFTA is band-limited; that is why there is not a “sufficiently large” region in the allowed m–z space. Breaking this limitation is not infeasible but may result in a reduction in image quality. If we need to consider ultra-wideband scenarios, the WDFs of certain transfer functions or PSFs may not be a simple Dirac delta function but could exhibit distorted tails on the space or frequency axis. In such cases, the highest order of the Hamiltonian might break the quadratic limitation, making the decomposition by LCTs exceedingly challenging. However, this is a separate topic and holds significant importance in the shaping and optimization of the PSF (or transfer function).

The design of the variable-scale model offers a novel paradigm for coherent imaging and display. Compared to ASM and SFT, it is not overly complex and can be implemented by FFT based algorithms, providing a powerful and efficient method for various diffraction calculations. However, SFTM is not always the optimal choice in every situation. An SFTM based single calculation requires two or three FFTs, whereas SFT only needs one. In far-field holography, if the image quality and magnification are not important but the generation time consumption is critical, the SFTM might not be the best option. In near-field holography, if the magnification is close to 1 and precise reconstruction of the diffraction field is required, ASM might be a better choice. This is because, in ultra-wideband scenarios, the form of the propagation matrix T might not be [1,λz;0,1].

In the design of full-color holography, each CGH corresponds to a specific color. If we want to design holographic videos, a control module for the synchronous switching of the three-color CGHs and the lasers is required. This reduces the holographic frame rate to one-third of the SLM refresh frame rate. Our future focus will be laid on the design of a single full-color hologram by SFTM. By simultaneously illuminating a single hologram with lasers of three colors, a full-color holographic image can be directly displayed, thereby maintaining the display frame rate of holographic videos. Complemented by GPU acceleration, the SFTM exhibits the potential to realize dynamic, high degree-of-freedom, and high frame-rate full-color displays in future applications.

In summary, from the perspective of WDF transportation characteristics, a framework was proposed under the principle of matrix decomposition. Using this framework, a Fourier analysis method named the SFTM with variable-scale diffraction calculation characteristics is derived, which offers much higher degrees of modulation freedom than traditional methods. The effectiveness of SFTM based algorithms have been demonstrated, highlighting their powerful scaling control and automatic chromatic aberration correction capability. Additionally, the SFTM is a superior method for 3D Fresnel holography compared to SFT or FM. A meta-atom design strategy is implemented that maps polarizations respectively to color channels. Near-zero crosstalk, variable image scaling, and full-color metasurface holography are achieved. This effectively resolves the longstanding issues of resolution and fidelity reduction in the metasurface community.

Finally, our variable-scale model is robust enough to provide a significant degree-of-freedom for most coherent imaging and displays applications. This makes it widely applicable in beam shaping, diffractive optical element design, CGH, 3D displays, HUDs, and meta-holography.