2.1 Design principle for generation of THz achromatic Airy beams

Rather than using a fixed spatial wavefront to generate a chromatic Airy beam, here we exploit a wavelength-dependent spatial phase based on a silicon metasurface to generate an achromatic Airy beam with the same trajectory. A schematic of the achromatic Airy beam generator is shown in Figure 1(a). The incident THz wave impinges vertically on the metasurface device and the transmitted THz wave follows the trajectory along the white dash line f(z). To realize the achromatic Airy beam, here we study the truncated Airy beam to acquire the wavelength-dependent spatial phase of the achromatic Airy beam. The Airy beam is generated by an initial field distribution in the metasurface plane at z = 0, given by ψ(x,z=0)=Ai(x) . Here, we employ the semiclassical approximation to analyze the Airy function. The standard asymptotic forms show that ψ(x,z=0) is exponentially small for x > 0, and hence negligible, and is oscillatory for x < 0, the precise expression being ψ(x,z=0)≈(−π2x)−1/4sin[23(−x)3/2+π4]=A+eikϕ++A−eikϕ−), where ϕ±=∓23(−x)3/2+π/4 and A±=±(i/2)(−π2x)−1/4 are the initial phases and amplitude functions. The initial phases ϕ+ and ϕ− give rise to rays that propagate along the +x and opposite −x directions. Specifically, the ray species ϕ+ emerges sideways, converging to a caustic along the trajectory f(z)=(az2)/4, where the variable a is equal to 1/(k2x03), and k and x₀ are the wavenumber and an arbitrary transverse scale. The details of the parabolic trajectory are derived in Note S1 in the Supplementary material. The geometric construction deriving the phase profile is illustrated in Figure 1(b). The parabolic trajectory is indicated by the black dash curve f (z) as shown, and we obtain the corresponding spatial phase function φ(x) at the metasurface plane with z = 0 that will generate the curve as a caustic. The caustic is constructed as the envelope to a family of rays ϕ+ such that each point x at the plane z = 0 can be functionally related to a point on the caustic via a tangent of slope θ, where tanθ=f′(z)=df(z)/dz. Since the tangent can be parameterized in terms of x using f′(z)=[f(z)−x]/z, we can determine the desired phase distributions by integrating the phase derivative condition:

Figure 1:

Schematic and phase profile of an achromatic Airy beam.

(a) Schematic of the experimental generation of terahertz (THz) achromatic Airy beams by a silicon metasurface. The transmitted left circular polarization (LCP) THz intensities in the x–z plane fit the parabolic trajectory f(z) for the right circular polarization (RCP) incident light. (b) To determine the phase distribution φ(x,f) of the metasurface under design, we draw an auxiliary line (black solid line) tangent to the parabolic asymptotical trajectory f(z) of the desired Airy beam (black dash line) at an arbitrary point, which intersects with x axis exhibiting an angle θ. With θ(x) function obtained by repeating such a process, we then obtain φ(x,f) based on Eq. (1). (c) Phase profile of an achromatic Airy beam within the frequency range f∈[fmin,fmax], where fmin=0.4 THz, fmax=0.8 THz, φshift=820∘ at 0.8 THz, and the arbitrary transverse scale x₀ = 650 µm in our design.

When the analysis in the preceding discussion is simplified in the paraxial approximation, we obtain the following wavelength-dependent phase profile:

where the value of φshift(f) represents an integration constant at the frequence of f. To realize an achromatic Airy beam within f∈[fmin=0.4 THz,fmax=0.8 THz] in Figure 1(c), we need to design a metasurface exhibiting the phase profile φ(x,f), as shown in Eq. 2, for every frequency within the band [fmin,fmax].

We now illustrate how to design such a metasurface exhibiting the desired phase profile. With φshift(f) term dropped for the moment, we decompose the remaining part of phase profile φ(x,f) into two parts, namely a basic phase profile φ1(x,f) representing the phase profile required by the device at f_min and φ2(x,f) the phase difference at other frequencies. Straightforward calculations yield that:

and

Now our task is to design a series of meta-atoms, which not only yield the required phases required by Eq. (3) at the frequency f_min, but also exhibit the different frequency dispersions as dictated by Eq. (4). However, φ2(x,f) exhibits negative slopes against frequency at every point x, which can not be realized by any resonating structures exhibiting normal frequency dispersions (i.e., transmission phases are increasing functions of frequency). Fortunately, there is an additional term φshift(f) in Eq. (2), which can be freely chosen to solve this issue. Specifically, set φshift(f)=α⋅(f−fmin), we find that:

Choosing an appropriate α value, we can make the term [−4πa1/23cx3/2+α] taking positive values at every position x, so that such frequency-dependent phase profiles can be realizable using resonating meta-atoms.

We proceed to search a series of meta-atoms exhibiting the desired phases (Eq. (3)) at frequency f_min and frequency-variation slopes as required by Eq. (5). The building structures of the silicon metasurface are presented in Figure 2(a) and (b), consisting of solid and inverse rectangular structures on a silicon substrate. To fulfill the complex requirements on the phase responses, we design the metastructures based on a combination of two mechanisms, namely, the resonance mechanism yielding a frequency-dependent transmission phase and the Pancharatnam-Berry (PB) one yielding a frequency-independent phase. Since the PB mechanism is adopted, we assume the incident THz wave to exhibit RCP, and employ finite-difference time-domain (FDTD) simulations to determine the width w and length l of our metastructures (with fixed lattice period p and etching depth t) based on two criterions: 1) transmission phases exhibit linear dependences on frequency with desired slopes; and 2) they possess relative high conversion efficiencies to LCP within the frequency band of interest. The section of method provides more simulation details. Finally, we carefully select 32 different structural parameters containing solid and inverse rectangular structures in Table S1 and present the individual phase responses and conversion efficiencies in Figure S1.To directly present the coversion efficiencies, we provide an efficiency map (as shown in Figure S2) of the selected meta-atoms with respect to different frequencies. Here as an example, Figure 2(c) presents the simulated phase differences of the solid structure (sequence number 15, 580°) and the inverse structure (sequence number 27, 725°) with a same etching depth t = 350 μm, showing the linear phase response as a function of f (or 1/λ). Note that the inverse structures in Table S1 have a wider phase response than the solid structures in Table S1 and can compensate for a larger phase difference, so that the positions of the inverse structures are close to x = 0 and those of the solid structures are close to the left edge of the sample, as shown in Figure 1(a).

Figure 2:

Phase distributions and conversion efficiencies of the solid and inverse structure.

Side views of (a) a solid and (b) an inverse silicon structures with length l, width w, etching depth t, and lattice constant p on a silicon substrate. (c) Simulated phase responses (blue curves) and RCP-to-LCP polarization conversion efficiencies (red curves).

Before closing this section, we emphasize that only the LCP component of the transmitted wave can acquire the PB phases under RCP incidence [49], [, 50], so that only this wave component can generate the desired achromatic Airy beam. Moreover, the metadevice thus constructed does not work for linear-polarization incidence under which the meta-atoms do not generate any PB phase.

2.2 Numerical calculation of THz achromatic Airy beams

To confirm the desired phase profile φ (x,f) in Eq. (5) derived by the geometrical construction method and reveal the underlying mechanism of the achromatic Airy beam, we performed numerical calculations based on Fraunhofer diffraction integration. In our numerical results, the structural region ranged along x axis from x = −7 mm to x = +4 mm. Figure 3(a) shows the phase profiles at frequencies f_min = 0.4 THz, f = 0.6 THz, and f_max = 0.8 THz. Note that the phase profile at 0.8 THz is smaller than that at 0.4 THz, which is in conflict with the relation between transmission phase response and frequency. This impasse is broken by introducing an additional phase shift φshift(f) that is larger than the phase difference Δφ=φ0.4THz−φ0.8THz=1350∘ at the position x = −7 mm. To compensate for the large phase difference, the building structures are no longer subwavelength. In this case, we resort only to the geometric phases of three different angles of rotation with respect to the laboratory coordinate system to match the phase profiles of Airy beams with different frequencies at 0.4, 0.6, and 0.8 THz. Table S2 in the Supplementary material shows the structural parameters and the rotated angle distributions. Subsequently, we performed a numerical calculation based on the initial field ℱ{eikφ(x)}, in which all metasurfaces in the array with period 70 μm were considered as subsources radiating cylindrical surface waves with designed initial phase profiles. Figure 3(b–d) shows the numerically calculated Airy beams and Figure S3(a–c) in the Supplementary material show the FDTD simulated beams at frequencies 0.4, 0.6, and 0.8 THz. We find that the trajectories of the Airy beams are in good agreement with the parabolic curve f(z), although they are both imperfect owing to the limited number of metasurface elements and the nonideal phase profiles.

Figure 3:

Fraunhofer diffraction integration of the phase profiles in Eq.(2) realized by PB phase.

(a) Phase profile of an achromatic Airy beam with additional phase shift φshift(f). (b), (c) and (d) Numerically calculated Airy beams at frequencies 0.4, 0.6 and 0.8 THz. The white dash lines are the parabolic trajectory curve f(z).

2.3 Experimental setup and fabricated samples

To validate the above theoretical analysis, we designed a silicon metasurface device for generating a broadband THz achromatic Airy beam. However, the theoretical design of achromatic Airy beams requires a phase compensation of 1350°, which cannot be compensated owing to the large aspect ratio in experimental fabrication. Therefore, we reduced the sample size, which ranges along the x-axis from −3.2 to 4 mm. In our design, we chose the additional phase shift φshift(f) to be 820° at 0.8 THz. To encode the phase profiles of the achromatic Airy beam within the frequency band, we placed the 32 selected structures at x coordinates with approximately equal phase differences, as shown in Figure S4, in the Supplementary material. The starting points of each vertical line is the initial phase φ₀ at 0.4 THz, and the height of each vertical line is the phase coverage (or phase difference). The unresolved problem is to realize the phase profile at 0.4 THz. Here, the phase profile is taken as the basic phase φ₁(x,f), determining the parabolic trajectory of the achromatic Airy beam within the frequency range. The basic phase φ₁(x,f) is related solely to f_min and is independent of the working frequency f. The PB phase, which is independent of the transmission phase, is used to acquire the phase profile at 0.4 THz. To obtain the correct basic phase φ₁(x, f), the spatial rotation angle θ of each structural unit is taken as θ=(φ1−φ0)/2 (for further details, see Table S3 in the Supplementary material)

To verify the reliability of our designed metasurface devices, we fabricated samples to generate achromatic Airy beams and explore the corresponding characteristics using THz near-field scanning microscopy (NFSM), as shown in Figure 4(a). The polarizer controls the polarization of the light output by the fiber laser, which is focused by the lens and reflected by the mirror, and finally illuminates the light on the probe. In the real system shown in Figure S5 in the Supplementary material, the collimated THz waves radiating from a 100 fs (λ=780 nm) laser-pulse-pumped photoconductive antenna emitter are modulated with an appropriately polarized state. A commercial THz near-field probe is positioned 2 mm away from the sample to detect the E_x of the transmitted LCP light, with the RCP light illuminating the metasurface device. The sample is fabricated on a silicon wafer by conventional lithography together with deep reactive ion etching. Further details of sample fabrication are given in methods. Figure 4(b) shows an optical microscope image of the whole of the fabricated sample. Figure 4(c) and (d) shows partially enlarged scanning electron microscope (SEM) images of the inverse and solid structures, respectively. Note that the inverse structures with a larger phase compensation are located toward the right of the metasurface device, while the solid structures with a smaller phase compensation are located toward the left side.

Figure 4:

Experimental setup of the terahertz (THz) near-field scanning microscopy (NFSM) and sample images.

(a) Schematic of THz NFSM. (b) Optical image of the whole sample. (c) Partial enlargement on the right side (red desh box) of the inverse structures. (d) Partial enlargement on the left side (red box) of the solid structures.

2.4 Experimental detection of THz achromatic Airy beams

Figure 5(a)–(c) shows the normalized abs(E_x) distributions of the transmitted LCP light obtained from the THz NFSM at the three frequencies. Note that there are no data for the first 2 mm in the experimental results because of a safety distance imposed during the experiment to prevent collisions. In these figures, the normalized abs(E_x) distributions show a good match with the parabolic trajectory (the white dash line) at 0.4 THz. To verify the experimental results, Figure 5(d)–(f) presents the normalized intensity distributions E_x simulated by FDTD at frequencies 0.4, 0.6, and 0.8 THz, respectively. Furthermore, in Figure S6 in the Supplementary material, we present the results of Fraunhofer diffraction integration at more frequencies, the results of simulation, and experimental data, all of which verify the achromatism of the Airy beams. For a more intuitive comparison between the simulation results and experimental data, we obtain the lateral offsets x of the maximum normalized intensity (with error bars corresponding to 95% of the maximum normalized intensity) at different propagation distances z and compare them with the numerical deflection. Figure 5(g)–(i) shows the results of this comparison at the three frequencies, and it can be seen that there is good agreement. In the simulation, we calculated the working efficiency of the device, as shown in Figure S7 in the Supplementary material. The working efficiency ranges from 20 to 60% and is higher at frequencies near 0.6 THz. Furthermore, we present the simulation results for the self-healing properties of the achromatic Airy beam in the case of scattering by silicon particles in Figure S8 in the Supplementary material.

Figure 5:

Achromatic Airy beam experiments.

(a)–(c) Experimental results on near-field detections at the three frequencies, showing that the intensity distributions of the left circular polarization (LCP) field transmitted through the metasurface devices match well with the white dash line indicating the parabolic trajectory. (d)–(f) Simulation results at the three frequencies. (g)–(i) Comparison of experiment, simulation, and theoretical deflections, where the error bars correspond to 95% of the maximum normalized intensity.

2.5 THz achromatic Airy-beam-based metalens

Airy-beam-based metalenses have many advantages for detection, such as large depth of focus (DOF) and self-healing properties, and so here we design a metalens based on autofocusing achromatic Airy beams. To realize this, the achromatic Airy-beam-based metalens, a metasurface device composed of two symmetrical structures is designed to generate two counter-propagating achromatic Airy beams. Figure S9 in the Supplementary material shows optical and SEM images of the metasurface device. To highlight the robustness of the Airy-beam-based metalens, an aluminum foil cover is placed over a 1 mm region (ranging along the x-axis from x = −0.5 mm to x = +0.5 mm) around the center of the experimental sample. Figure 6(d)–(f) shows the normalized abs (E_x) distributions of the metalens at frequencies 0.4, 0.6 and 0.8 THz, respectively. Experimental results are presented on the left and FDTD results on the right. Experimental results at other frequencies are shown in Figure S10 in the Supplementary material. In particular, it should be noted that there is less energy in the region covered by aluminum foil, which is quite different from what occurs with a traditional metalens. The intensity distributions in the vertical direction at the center of the two cases can be found in Figure S11 in the Supplementary material. A comparison shows that the Airy-beam-based metalens has the advantage of a larger DOF and detection with a high signal-to-noise ratio. Interestingly, the starting points of the autofocusing effect fit well with the intersection of the parabolic trajectory. When two counter-propagating achromatic Airy beams begin to intersect, constructive interference occurs, and the trajectory curves coincide with the lower edge of the normalized intensity distributions of the focal spots. Furthermore, because the full width at half maximum (FWHM) of the focal spots is the main indicator to evaluate the performance of Airy-beam-based metalenses, we compared the experimental and theoretical distributions at a focal length f = 10 mm and measured the FWHM at the three frequencies, as shown in Figure 6(a)–(c). The FWHM is 0.84 mm at 0.4 THz (λ = 0.75 mm), 0.52 mm at 0.6 THz (λ = 0.5 mm), and 0.38 mm at 0.8 THz (λ = 0.375 mm). The experimental results agree well with those of the simulations, and the FWHMs are approximately equal to the wavelength of incidence of the THz wave. The self-healing property is typically verified by placing an obstacle on the travel path of the parabolic trajectory. The trajectory recovers quickly if two silicon obstacles of diameter 400 μm are located at (x, z) = (−0.8 mm, 6 mm) and (0.8 mm, 6 mm), indicating a rather robust self-healing property of the generated Airy-beam-based metalens. The self-bending, diffraction-free, and self-healing achromatic properties are illustrated in Figure S12 in the Supplementary material.

Figure 6:

Achromatic Airy beam-based metalens.

(a)–(c) Comparison of simulated and experimental results for the intensity distributions at the focal length. The full width at half maximum (FWHMs) of the intensity distribution are 0.84 mm at 0.4 THz, 0.52 mm at 0.6 THz, and 0.38 mm at 0.8 THz. (d)–(f) Demonstration of the generation of two counter-propagating achromatic Airy beams in experiments and simulations. The focal length f = 10 mm.

4.1 Simulation details

FDTD was used to find and simulate the structural units and the transmission E_x field of the metasurface devices. In the simulation of the unit, the mesh grid size was set to 10 nm. In the simulation of the whole metasurface device, the grid size was set to 1 μm

4.2 Sample fabrication

All the silicon metasurface devices were fabricated using conventional lithography together with deep reactive ion etching. The refractive index n of silicon is 3.45 and its resistance R is higher than 10⁴ Ω cm. First, a 2-μm-thick silica layer was grown on a 500-μm-thick double-side-polished high-resistivity silicon wafer. Next, by conventional photolithography, the left photoresist and the silica became a double layer of protection above the silicon posts. Conventional deep reactive ion etching was then employed to make silicon posts of thickness t = 350 μm. Finally, the remaining photoresist and the remaining silica were cleaned off separately.

4.3 Experimental characterization

THz near-field scanning microscopy was employed in the experimental tests owing to its high scanning speed and high resolution. Here, only the x-polarized electric field E_x component was measured with RCP illumination. The electric field was detected at 0.2 mm intervals from −3.2 to 4 mm in both the x and z directions. Note that we reserved 2 mm to avoid collisions between the probe and the metasurface devices.