The metasurface is a category of artificially designed structures, which have powerful abilities to simultaneously manipulate multiple parameters of light including amplitude, phase, polarization, orbit-angular-momentum, etc. [1–3]. The multi-dimensional wavefront modulation based on metasurfaces is usually unachievable by traditional optical elements. Such intriguing property has attracted enormous interest in related areas such as the achromatic metalens [4,5], holography [6–8], display [9,10], quantum optics [11,12], and polarization imaging [13,14].

In particular, the holographic images reconstructed by metasurfaces can have high resolution with good quality, large field of view, and freedom of the disturbance by the unwanted diffraction orders [15,16]. Meanwhile, metasurfaces have also promoted the investigations from scalar holography to vectorial holography [17]. Segmented metasurfaces have been demonstrated by dividing the metasurfaces into different segments according to the desired polarization states, which can reconstruct holographic images with various specific types of polarization states [18,19]. Single layer [20–23] or bilayered metasurfaces [24] are utilized to manipulate partial or all degrees of freedom of a $2\times 2$ Jones matrix. Hence, different independent holographic images can be obtained in different polarization channels. For making the holographic images exhibit more abundant polarization characteristics, the relations between the desired polarization distribution and the meta-atoms of the metasurfaces are established [25–27]. And the obtained vectorial holographic image can exhibit arbitrary polarization distribution based on the designed metasurface without spatial multiplexing. Meanwhile, the methods of generating multiple vectorial holographic images along the propagation direction are previously neglected, which is in favor of generating high-dimensional vectorial fields and leading various related applications including sensing, polarization imaging, and optical trapping [28–30].

Now, the metasurfaces usually consist of metallic or dielectric meta-atoms fabricated by electron beam lithography, focused ion beam lithography, and nanoprinting methods [31–33]. Meanwhile, these fabrication methods are limited in the fabrication of three-dimensional structures with freeform and controllable shapes. Femtosecond laser direct writing (Fs-LDW), as a flexible and high-precision 3D processing technology, has been widely used in the fabrication of micro-optical elements [34–36]. It can also be adopted to fabricate metasurfaces and exhibit some unique advantages. The femtosecond laser can achieve single point exposure of the photoresist through nonlinear two-/multi-photon polymerization, which guarantees sub-micrometer fabrication resolution. Arbitrary 3D structures can be processed based on Fs-LDW, which is in favor of providing more design freedoms of demonstrated meta-atoms to realize desired wavefront modulations [37,38]. It also brings the possibilities of fabricating the metasurfaces with freeform meta-atoms [39,40] or on the curved substrates [41]. Meanwhile, the metasurfaces fabricated by Fs-LDW are easier to integrate with other optical elements such as an optical fiber [42,43], photodetector [44], or complementary metal–oxide semiconductor (CMOS) [45,46] to form more complex optical devices with advanced functionalities. Furthermore, high aspect ratio structures can be fabricated by selecting the photoresist with larger Young modulus and optimizing the processing parameters of Fs-LDW.

Here, we demonstrated a novel method for realizing multi-plane vectorial holography based on the metasurface fabricated by Fs-LDW. We establish a Jones matrix framework and impose different amplitude and polarization restrictions at each image plane for empowering the desired polarization features. By encoding the generated unitary matrix hologram to the metasurface, each reconstructed image can exhibit inhomogeneous and customized polarization distribution which is verified by rotating a polarization analyzer in the experiment. Meanwhile, the demonstrated metasurface is composed of IP-L nanofins with different lengths, widths, heights, as well as orientation angles. The fabrication of the demonstrated metasurface represents the advantages of laser direct writing including large area fabrication as well as 3D structure processing. Furthermore, it is easy to fabricate structures that can be integrated with an optical fiber, photodetector, CMOS, or other devices, and promote the development of related areas including optical trapping, sensing, imaging, and so on.

The basic design principle of our demonstrated method is illustrated in Fig. 1. By integrating the multi-plane holographic algorithm with the Jones matrix method, a Jones matrix framework is established with the capability of tailoring the polarization distribution of each image plane. Under $x$-linearly polarized light illumination, various independent holographic images can be successfully obtained with customized polarization distributions by imposing the corresponding amplitude and polarization restrictions in the framework. Hence, the modulation of multiple parameters with longitudinally varied features can be realized. And the unitary part of the generated Jones matrix hologram is encoded to the fabricated single layer metasurface composed of IP-L nanofins with different lengths, widths, and heights by analyzing the corresponding eigenvalues and eigenvectors at each pixel.

Based on the matrix polar decomposition, any $2\times 2$ Jones matrix $\tilde{\mathit{J}}$ can be represented by the production of a Hermitian matrix $\tilde{\mathit{H}}$ (diattenuator item, polarizer behavior) and a unitary matrix $\tilde{\mathit{U}}$ (retarder item, waveplate behavior). The matrix $\tilde{\mathit{H}}$ and $\tilde{\mathit{U}}$ can be considered as the analogues or high-dimensional representations of the amplitude and phase of the optical field [47]. Hence, the polarization-dependent amplitude and phase modulation represented by the matrix $\tilde{\mathit{H}}$ and $\tilde{\mathit{U}}$ can be utilized to realize the desired vectorial optical field.

The flowchart of our demonstrated Jones matrix framework for generating multi-plane vectorial holographic images is shown in Fig. 2. At first, a matrix distribution ${\tilde{\mathit{U}}}_{\mathrm{meta}}$ composed of $640\times 640$ random unitary matrices is generated for the process of initialization. The propagation from the hologram plane to each image plane is obtained based on the Fresnel diffraction formula in the form of Fourier transform (FrT) [48]: $${\tilde{\mathit{D}}}_n(x_o,y_o)=\frac{\exp (jk{z}_n)}{j\lambda z_n}\exp (\frac{jk}{2z_n}(x_o^2+y_o^2))\times F({\tilde{\mathit{U}}}_{\mathrm{meta}}(x_h,y_h)\exp (\frac{jk}{2z_n}(x_h^2+y_h^2)\left)\right),$$(1)where ${\tilde{\mathit{D}}}_n$ ($n=1,2,\dots ,N$ represents different image planes) is the Jones matrix distribution at the $n$th plane, $\lambda$ is the wavelength of incident light, $z_n$ represents the predefined reconstructed distance of the $n$th plane, and ($x_h$, $y_h$) and ($x_o$, $y_o$) refer to the coordinates of the hologram plane and the image plane. At each image plane, the Jones matrix distribution ${\tilde{\mathit{D}}}_n$ is calculated parallelly and the corresponding unitary part ${\tilde{\mathit{U}}}_n$ is acquired based on matrix polar decomposition. We extract the upper left phase item of the unitary Jones matrix at each pixel of different image planes and constitute the random phase profile ${\phi }_{11\_n}$ which can slightly relieve the crosstalk between neighboring image planes. Meanwhile, the amplitude distribution $A_n$ and polarization response ${\tilde{\mathit{P}}}_n$ are imposed as essential restrictions according to the desired amplitude and polarization distributions of the target images. Then, the Jones matrix distribution ${\tilde{\mathit{J}}}_{\mathrm{meta}}$ at the hologram plane can be obtained after the backward propagation process based on the inverse Fresnel transform (IFrT) method expressed by Eq. (2): $${\tilde{\mathit{J}}}_{\mathrm{meta}}=\sum _{n=1}^N\mathrm{IFrT}\{A_n{e}^{i{\phi }_{11\_n}}{\tilde{\mathit{P}}}_n\}.$$(2)

At the hologram plane, the unitary part of ${\tilde{\mathit{J}}}_{\mathrm{meta}}$ is extracted and used as the updated unitary Jones matrix distribution for the next iteration after normalization. After 30 iterations, the desired unitary matrix hologram ${\tilde{\mathit{U}}}_{\mathrm{meta}}$ is generated by discarding the corresponding Hermitian part of ${\tilde{\mathit{J}}}_{\mathrm{meta}}$. Such a discarding process is similar to generating the phase-only hologram from the calculated complex amplitude optical field, which is in favor of realizing the desired polarization features along propagation direction within a single layer metasurface.

At each pixel of the hologram plane, the unitary matrix ${\tilde{\mathit{U}}}_{\mathrm{meta}}(x_h,y_h)$ can be realized by an elaborately designed nano-structure [49]: $${\tilde{\mathit{U}}}_{\mathrm{meta}}(x_h,y_h)=R(\theta (x_h,y_h))\left[\begin{array}{cc}{e}^{i{\phi }_x(x_h,y_h)}& 0\\ 0& e^{i{\phi }_y(x_h,y_h)}\end{array}\right]R(-\theta (x_h,y_h)),$$(3)where ${\phi }_x(x_h,y_h)$ and ${\phi }_y(x_h,y_h)$ represent the polarization-dependent phase shifts acquired by calculating the eigenvalues. And the rotation matrix $R(\theta (x_h,y_h))$ is determined according to the corresponding eigenvectors. Hence, the obtained unitary matrix hologram ${\tilde{\mathit{U}}}_{\mathrm{meta}}$ can be simply encoded to the fabricated metasurface composed of nano-structures with $C_2$ in-plane rotational symmetry.

As shown in Fig. 3(a), an IP-L nanofin on a glass substrate is chosen as the unit cell to compose the designed metasurface. Under $x$- and $y$-linearly polarized light illumination, the IP-L nanofin can implement the desired phase modulations ${\phi }_{xx}$ and ${\phi }_{yy}$ by adjusting the length ($L$), width ($W$), and height ($H$) of the nanofin. Such form-birefringent phenomenon is caused by tailoring the effective refraction index of the nanofin. The rigorous coupled wave analysis (RCWA) method is adopted to realize a 3D parameter optimization by sweeping the $L$, $W$, and $H$ of the nanofin. The length $L$ and width $W$ are both changed in the range of 350–800 nm, with an interval of 10 nm. And the height $H$ is changed in the range of 3–5 μm, with an interval of 0.2 μm. The period $P$ ($P=1.1\mathrm{\mu m}$) is determined by considering the difficulty of fabrication as well as the range of realized phase modulations. The realization of multi-height levels of the nanofin benefits from 3D arbitrary structure processing offered by the Fs-LDW method. Meanwhile, the multi-height levels of the nanofin also provide extra design freedom that can expand the range of desired phase modulations and relieve the deviations between the simulated phase shifts and desired phase shifts. The working wavelength is set as 800 nm in the simulation. And the corresponding refractive indices of the IP-L photoresist and substrate are set as $n_{\mathrm{IP}-\mathrm{L}}=1.508$ and $n_{\mathrm{sub}}=1.45$, respectively. The simulated amplitude [$\mathrm{abs}(t_{xx})$, $\mathrm{abs}(t_{yy})$] and phase (${\phi }_{xx}$ and ${\phi }_{yy}$) of the transmission coefficients $t_{xx}$ and $t_{yy}$ of IP-L nanofins with different dimensions are shown in Figs. 3(b)–3(e). At each pixel of the designed metasurface, the corresponding orientation angle and dimension of the nanofin are determined based on Eq. (3) according to the generated unitary matrix hologram ${\tilde{\mathit{U}}}_{\mathrm{meta}}$. The orientation angle of each nanofin is determined by the rotation matrix $R(\theta (x_h,y_h))$. Meanwhile, the dimension of nanofin is determined from the simulated results by guaranteeing the minimum of the error $\epsilon =\mathrm{abs}(t_{xx}-\exp (i{\phi }_x))+\mathrm{abs}(t_{yy}-\exp (i{\phi }_y))$.

For verifying the feasibility of our demonstrated method, we fabricated the designed height tunable metasurface on the indium tin oxide-coated soda lime glass substrate by utilizing a commercial femtosecond laser direct writing system (Photonic Professional GT2, Nanoscribe). The Fs-LDW is based on the principle of two-photon polymerization which converts small unsaturated molecules in a liquid state into solid macromolecules through the polymerization reaction of two-photon absorption. The process of two-photon absorption is always confined to the volume of the focal point. Hence, only in the photoresist domain is the laser focusing cured. And the unpolymerized resins can be easily removed with a developer, leaving the remaining 3D structure after femtosecond laser exposure. Meanwhile, the printing resolution of femtosecond laser direct writing can exceed the diffraction limit.

The Photonic Professional GT2 adopts the femtosecond laser with the wavelength of 780 nm, repetition frequency of 80 MHz, and 80–100 fs pulse width. A high numerical aperture oil immersed objective lens ($63\times \text{/NA}=1.4$) is used to realize high accuracy of fabrication. We choose the IP-L photoresist to constitute the demonstrated metasurface due to its high Young’s modulus, which is more suitable for fabricating nanofins with high aspect ratio. In the process of fabrication, the optimized power of the femtosecond laser was 49.5 mW and the scanning speed of the galvanometer was 4000 μm/s. In order to improve the mechanical strength of the IP-L nanofin, small processing parameters of hatching and slicing distances were used to fabricate the metasurface. (The movement steps of the laser in the transverse and longitudinal directions were set to 20 nm and 50 nm.) After the process of exposure, the sample was soaked in the propylene glycol monomethyl ether acetate for 20 min, and then placed in isopropyl alcohol for 3 min. Finally, the sample was dried by a supercritical dryer for 3 h in order to decrease the influence of capillary force and surface tension caused by liquid evaporation from nanofins.

The fabricated height tunable metasurface is composed of $640\times 640$ IP-L nanofins with different dimensions ($L$ and $W$, 350–800 nm, $H$, $3–4.6\mathrm{\mu m}$, and $P$, $1.1\mathrm{\mu m}$) as well as orientation angles ($\theta$, $-90°\mathrm{to}90°$). The realized aspect ratio of fabricated nanofins can reach up to 13. The area of the fabricated metasurface is $704\mathrm{\mu m}\times 704\mathrm{\mu m}$, and the corresponding scanning electron microscopy images from a top view and a side view are depicted in Figs. 4(a) and 4(b). We randomly select 10 nanofins in the same region and compare their size deviation between the fabricated nanofins and designed nanofins. The average deviations of the length and width of the fabricated nanofins are 3.58% and 23.17%, respectively. The larger deviation of width is due to the width of the fabricated nanofins being closer to the resolution limit of printing. Meanwhile, the overall fabrication deviation is less than 80 nm, which is acceptable in our current design. The experimental setup utilized to capture the reconstructed multi-plane vectorial holographic images is illustrated in Fig. 4(c). Light from a supercontinuum laser source (NKT Photonics Superk EVO) passing through a linear polarizer LP1 and a half-wave plate HWP illuminates the fabricated metasurface sample. LP1 and HWP are used together to manipulate the polarization orientation of incident linearly polarized light without changing its intensity. An objective lens ($20\times \text{/NA}=0.45$) as well as a CCD constitutes the imaging system. By adjusting the distance between the metasurface and the imaging system based on a manual displacement stage (Daheng Optics, GMM-T13M2C), the reconstructed images located at predefined distances can be captured. Another linear polarizer LP2 is used as a polarization analyzer to discriminate the inhomogeneous polarization distributions of the reconstructed multi-plane vectorial holographic images.

In the process of designing the metasurface, the capitalized alphabets “3D,” “PR,” and “INT” are chosen as the target images. The reconstructed images are located at different positions ($z_1=1.2\mathrm{mm}$, $z_2=1.45\mathrm{mm}$, and $z_3=1.7\mathrm{mm}$) along the longitudinal direction. In our algorithm, we impose the desired polarization restrictions to the three reconstructed images. Hence, the images “3D,” “PR,” and “INT” can exhibit the polarization behaviors of half-wave plates with continuously varied orientations of fast axis that changed from 0° to 45° in the horizontal, diagonal, and vertical direction, respectively. And the imposed polarization restriction of the pixel ($x_o$, $y_o$) at the $n$th image plane can be expressed as $${\tilde{\mathit{P}}}_n(x_o,y_o)=\left[\begin{array}{cc}i\cos (2{\theta }_n(x_o,y_o))& i\sin (2{\theta }_n(x_o,y_o))\\ i\sin (2{\theta }_n(x_o,y_o))& -i\cos (2{\theta }_n(x_o,y_o))\end{array}\right],$$(4)where ${\theta }_n$ denotes the orientation of the fast axis. Furthermore, the polarization distribution of the reconstructed image $|{\lambda }_{\mathrm{out}}^n(x_o,y_o)⟩$ ($|{\lambda }_{\mathrm{out}}^n(x_o,y_o)⟩={\tilde{\mathit{P}}}_n(x_o,y_o)|{\lambda }_{\mathrm{in}}$) at the $n$th image plane is determined by the imposed polarization restrictions $({\tilde{P}}_1,{\tilde{P}}_2,{\tilde{P}}_3)$ as well as the polarization of incident light $|{\lambda }_{\mathrm{in}}⟩$. Under the illumination of $x$-linearly polarized light, the desired inhomogeneous polarization distributions (continuously varied linear polarization states) of the reconstructed images are illustrated in Fig. 1. If the desired vectorial holographic images need to be reconstructed under other incident polarization states, our demonstrated scheme is still effective by adjusting the imposed polarization response ${\tilde{\mathit{P}}}_n(x_o,y_o)$ according to the desired $|{\lambda }_{\mathrm{out}}^n(x_o,y_o)⟩$.

The simulated and experimental results of the multi-plane vectorial holographic images under different configurations of LP1 and LP2 are shown in Fig. 5. By adjusting the distance between the fabricated metasurface and the imaging system, three different vectorial holographic images “3D,” “PR,” and “INT” located 1.16 mm, 1.42 mm, and 1.69 mm away from the metasurface are successfully captured. The reconstructed distances as well as the quality of reconstructed images are as expected. For verifying the inhomogeneous polarization distributions of the reconstructed vectorial holographic images in the experiment, the transmission axis of polarization analyzer LP2 is rotated to 0°, 30°, 60°, and 90°, successively. As shown in the right four columns of Fig. 5, specific parts of the reconstructed images become extinct when rotating LP2. And the experimental results agree well with the simulated results. Meanwhile, the correct extinct phenomena successfully prove the realization of the desired inhomogeneous polarization distributions of reconstructed multi-plane vectorial holographic images. We also note that the background noise of the left two columns of Fig. 5 is relatively obvious. And its polarization is consistent with the incident polarization. The observed background noise can be interpreted by the diffraction of the splicing deviations between neighbor spliced regions when fabricating the metasurface based on Fs-LDW. It can be relieved by compensating the splicing error in the process of fabrication or substituting the galvonometer mirror with a larger scanning field. For the case in which the same ($x_o$, $y_o$) coordinates in different planes require different polarization states, the obtained polarization distribution of the $n$th image plane will be slightly affected by the depth of field (DoF) of neighboring reconstructed images. Introducing random phase [50], using structured light [51] and adopting neural networks method [52] may improve the performance in such cases. In our current design, the imposed polarization response ${\tilde{\mathit{P}}}_n(x_o,y_o)$ is determined based on the desired polarization distribution of the reconstructed holographic image. And we fix the incident polarization as $x$-linearly polarized light. For the situations in which the metasurface is illuminated under other incident polarization states, the obtained polarization distribution can be calculated.

In the experiment, the measured transmission efficiency defined by the ratio of transmitted power to the incident power of the fabricated metasurface is 30.5%. Furthermore, our demonstrated method is not restricted to generating linear polarization states with different orientations. Circular or elliptical polarization states can also be obtained by imposing specific polarization restriction through the polarization transformation expressed by a $2\times 2$ symmetrical Jones matrix. And the number of reconstructed image planes is not limited to three. A three-dimensional vectorial holographic image can also be realized based on our demonstrated algorithm.

In summary, we have demonstrated a novel method for achieving multi-plane vectorial holographic images based on the metasurface fabricated by Fs-LDW. By integrating the multi-plane holographic algorithm with the Jones matrix method, a Jones matrix framework is established with the capability of tailoring the polarization distribution of each image plane by imposing desired amplitude and polarization restrictions. Various vectorial holographic images with customized polarization distributions can be reconstructed. Meanwhile, the calculated unitary matrix hologram can be easily implemented by the single layer metasurface composed of IP-L nanofins with different lengths, widths, heights, as well as orientation angles through analyzing the corresponding eigenvalues and eigenvectors. Our demonstrated method provides a feasible and convenient way to generate high-dimensional structured light with polarization, amplitude, and phase variations during propagation from the perspective of the holographic method. Such functionality is highly desired for generating a complex optical field based on the compact optical systems.

Furthermore, the Fs-LDW method empowers the possibilities of fabricating a 3D arbitrary structure with high aspect ratio, fabricating a large area metasurface, and fabricating a curved substrate. It is also in favor of integrating the designed metasurface with optical fibers and photodetectors, which may pave the way to achieve many meaningful applications including polarization imaging, sensing, optical trapping, and so on.