## Abstract

In this paper, the focal shift issues of squared metasurface lenses are investigated. Axial intensity distribution formula of the squared planar lens model is obtained by utilizing traditional diffraction theories combining with the Fresnel approximation method. Fresnel integral and Cornu spiral are adopted to obtain the relations between the relative focal shift and Fresnel number, which can be used to predict the shift for different Fresnel numbers. In the far infrared region (at 10.6μm), a group of C-shaped nanoantennas are designed to cover the phase shift from 0 to 2π and simulations also performed by using finite-difference time-domain (FDTD) method. Several lenses are arranged by those resonators, and simulation results of focusing performance are in good agreement with the theoretical prediction. It’s expected that this work will provide a better solution for the design of lens in the infrared integrated optical systems.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Metasurfaces is a kind of novel two-dimensional artificial thin layer structure composed of subwavelength units. Owing to the flexible modulation capability on electromagnetic waves, it has attracted extensive attentions. Yu et al. deduced generalized Snell’s law based on Fermat’s principle in 2011 [1]. A group of V-shaped nanoantennas were designed, anomalous reflection and refraction phenomena have been observed and demonstrated experimentally. On this basis, in recent years, this flexible method in controlling electromagnetic wave is widely applied to designing a series of new type of optical components such as vortex plates [1,2], high resolution holograms [3–5], plasmonic coupled devices [6,7], and focusing lenses [8–16]. Moreover, Cong et al. [17] also realized the functionality of any desired output polarization state with the dispersion-free operation by using grating-based metadevice in a broadband range of the far infrared region. As an important optical component, micro focusing lens has been used in the field of integrated optical circuit, optical communications, and optical measurement etc. Since the size of each antenna unit is usually at the subwavelength scale, thus the metasurface lens composed by these antennas is of wavelength magnitude, which can satisfy the requirements of miniaturization and integration in modern optics very well. Besides, it’s much easier to fabricate compared with conventional microlens. However, due to the fact that when the size of lens is at scale of wavelength, actual focusing performance will be affected by the diffraction effect of lens edge, leading to a deviation of focus position from the designing position, which has already been demonstrated in [18–20]. Nowadays, amplitude-phase retrieval algorithms are commonly adopted to design diffractive optical components, the metasurface lenses are also included [11,14]. The focal length in the theoretical designing have also been validated in good agreement with the simulation and experiment results. However, this kind of iterative technique is of heavy workload and time-consuming, which will waste lots of energy of researchers. It’s easy and convenient to design metasurface lens by diffraction method, while affected by edge diffraction effect, the focal shift is still a problem in real application [9,13,15].

Some research have been conducted on the focal shift issues caused by diffraction effect. In [18], a detailed analysis was carried out on the diffraction of converging spherical wave under a circular aperture based on the Huygens-Fresnel principle by Li and Wolf. The relationship curve between Fresnel number and relative focal shift was also obtained. After that, both Fresnel number and *f* number were considered to analyze the focal shift behavior by Li [19], and the results showed that the former having a qualitative influence on the shift than the later. Soon afterwards Li and Platzer [20] investigated the field distribution of focusing system in the low-Fresnel-number situations, which demonstrated the existence of the focal shift in the experiment. Later, Li and Wolf [21] also discussed the general trend of the changes in the structure of the focal region as the Fresnel number decreases. Furthermore, the Fresnel number was also introduced by Yu and Zappe [22] in their research on focal shift of plasmonic lens which based on metallic nanoslit arrays, some suggestions were given to design the lenses. Besides, they also compared three methods of predicting the focal length provided by previous papers in low-Fresnel number systems, and an active plasmonic lens was also designed to observe its focusing performance [23]. Meanwhile, based on the research conducted by Yu and Zappe, Gao et al. [24] calculated the axial field intensity distributions of lens composed of metallic nanoslit arrays based on scalar diffraction theory. The relation curves between Fresnel number and relative focal shift were also plotted, which can be used to predict the real focal length of nanoslit lens. Although the factors influencing focal shift of lens have been investigated by previous studies and some theoretical analysis were given, to the best of author’s knowledge, there are few research on squared metasurface lens. As a kind of novel artificial structure, there is a bright future prospect for metasufaces. Therefore, it’s meaningful to carry out this study.

In the present paper, focal shift issues based on squared metasurface lens designed by diffracted method are discussed with the traditional diffraction theories. Fresnel diffraction integral is adopted to obtain the axial field distribution of squared planar lens, and the relations between Fresnel number and relative focal shift is provided. A group of C-shaped resonators are designed in the far infrared region, which can realize the phase modulation on transmitted wavefront from 0 to 2π with constant amplitude transmittance. At the specific wavelength of 10.6μm, several squared metasurface lenses are also designed, and simulated by FDTD method. Simulation results are compared with theoretical analysis, finding that the predicted focal lengths by theoretical all possess a high precision for different Fresnel numbers. The results illustrated in this paper may provide an instructive guide for the design of squared metasurface lens.

## 2. Squared planar lenses

#### 2.1 Theoretical analysis

In order to analyze the focusing performance of squared planar lens, assume the lens is located in the *xoy*-plane, where point *o* is the geometrical center of lens and optical axis is *z*-axis. An arbitrary point is chosen as a point source at *z* > 0, and Fresnel diffraction method can be employed to obtain complex amplitude distribution of squared lens. Corresponding phase distribution can be used to design the lens, and theoretically the focus would be observed at the preset location when the lens is stimulated by a plane wave. Different from conventional optical lens to shape the wavefront by gradual phase accumulation, realization of function of metasurface lens is based on the discontinuous phase modulation of separated antenna units, therefore, the model of planar lens is considered here, as shown in Fig. 1.

As is well known to all, when a monochromatic, uniform, converging spherical wave is diffracted at a circular aperture, as the condition (*a* / *f*)^{2} << 1 is satisfied simultaneously, the Fresnel number, which is commonly defined as the number of Fresnel zones that fill the aperture when the aperture is viewed from the geometrical focus, can be expressed as:

Where *a* is the radius of the circular aperture, *f* is the distance from the geometrical focus to the plane of aperture. Besides, when the *a* >> *λ* is also satisfied, owing to the influence of edge diffraction, as described in the previous papers [18], the real position of the maximum intensity point along the axis for the circular aperture can be predicted by introducing the definition of Fresnel number in Eq. (1). However, for the Fresnel number of the squared planar lens here, it is necessary to give it a realistic definition.

To obtain the definition of Fresnel number of squared aperture, the Fresnel number for the part of maximum inscribed circle of the square and the rest area are calculated respectively. Assume that *a* is the half length of the aperture. The Fresnel number for the area of maximum inscribed circle *S*_{1} can be obtained through Eq. (1). The area of boundary region *S*_{2} is acquired by subtracting the *S*_{1} from the square, which is about 0.8584*a*^{2}. For a circular aperture, since the area of each Fresnel waveband is identical, the effect of boundary area *S*_{2} can be treated as the equivalent number of wavebands of the circular aperture. It can be calculated the ratio of *S*_{2} to *S*_{1} is about 0.2732, which suggests that the equivalent number of wavebands of *S*_{2} is about 0.2732 times that of the actual Fresnel number of *S*_{1}. Therefore, the equivalent formula of Fresnel number for squared lens is given:

Here *a* is the half length of the squared lens, *λ* is the wavelength in vacuum, and *f* is the preset focal length. By comparing the pattern of Fresnel wavebands on circular and squared aperture, it is noted that the contribution of the wavebands in *S*_{2} to the focus along the axis is very complicated in calculation, since the amplitude distribution at the focus caused by these incomplete annular wavebands depends not only on the number of the wavebands, but also on the proportion of each waveband. Besides, the area of *S*_{2} is relatively small in the squared aperture. Therefore, when there is a significant focal shift with condition of *a << f*, it is reasonable to ignore the influence of *S*_{2} on the actual focusing performance. Accordingly, when analyzing focal shift issues of squared planar lens, the definition of Fresnel number of circular aperture will be adopted. It could be proved in the following discussion that, though this approximate calculation method is taken, the prediction results from theoretical analysis are still have good agreement with the simulations.

When the focus lies on the optical axis, take Fresnel approximation into consideration, for a squared lens, the complex amplitude distribution of arbitrary point *z*(0,0,*z*) along the optical axis at *z* > 0 can be described as:

*U*

_{0}is the complex amplitude distribution of the lens at

*z*= 0, when the plane wave incident, field distribution on the transmitted side of lens should be:where

*A*represents amplitude constant, phase distribution of the lens can be obtained by geometrical optics, and it’s also determined by the location of preset focus along the axis:

When the condition *a* << *f* is satisfied, Eq. (5) can be approximated and simplified by using binomial expansion:

Substitute Eq. (6) into Eq. (4), then Eq. (3) can be transformed into:

It is clear from Eq. (7) that a new phase factor appears in the exponential term, which can be defined as the relative position of arbitrary point *z* along the optical axis:

When *z* value is considered as the theoretical position of focus under the influence of diffraction, due to the focal shift is always moves toward the lens [18], then the values of *p* are always positive. According to derivation method of Fresnel diffraction integral in [25], together with Eq. (7), light intensity distribution of arbitrary point along the optical axis is expressed as:

Besides, *C*(*w*) and *S*(*w*) are Fresnel cosine and sine integral, respectively:

*w*is upper limit of the two integrals, through calculation, which can be presented as:

Substitute Eq. (1) and Eq. (8) into Eq. (12), *w* is transformed into:

#### 2.2 Focal shift

It can be seen from above, axial intensity distribution is just a function of the relative position *p* of observation point. Since the Fresnel number is used to characterize the number of wavebands on the diffraction screen, here it can be utilized to describe the corresponding planar lenses with different parameters. Besides, the ratio of *I*/*I*_{0} in Eq. (9) equals to the biquadrate of the distance from arbitrary point on the Cornu spiral to the coordinate origin. To acquire the theoretical position of focus of squared planar lens, the point of maximum intensity *I*_{max} along the optical axis is taken. For intensity formula *I*, the *I*_{0} is a function of *p*, so if the maximum of *I*/*I*_{0} is expected, then it is required to find a point on the Cornu spiral, to ensure the longest distance to the coordinate origin. As described in [26], when the value of upper limit *w* takes 1.2 in Eq. (11), the coordinate of [*C*(*w*),*S*(*w*)] is (0.7154,0.6234), and is (0.6386,0.6863) when *w* takes 1.3. Moreover, complex Simpson formula is adopted to plot the Cornu spiral with the posteriori error of 10^{−12}, which is shown in Fig. 2(a). The spiral is zoomed in partially, as shown in Fig. 2(b). The point on the curve which satisfies maximum distance to the origin located in (0.7069,0.6332). Meanwhile, due to *w* represents the arc length from origin to arbitrary point on the curve, by using curve integration, *w* value of this point is ensured as 1.21.

With Eq. (13), relation between *N* and *p* is obtained:

From the Eq. (14), the relative focal shift of focus only has to do with Fresnel number *N*, thus the function curve between *p* and *N* is able to be plotted, which is shown in Fig. 3.

It is worth noting that in the previous literatures, in order to obtain the relation curves between the relative focal shift *p* and Fresnel number *N*, methods by taking the derivative of axial intensity distribution [18–20] and by using series expansion [24] were commonly adopted. However, with the help of Fresnel integrals and Cornu spiral, function of axial maximum intensity is ensured from Eq. (9), certain relations between *N* and *p* can also be deduced. The method is simple and feasible, which can be used for predicting the real focal length of squared planar lens.

In Fig. 3, it can be seen that relative focal shift decreases gradually as the Fresnel number increases, nonlinear inverse relationship is satisfied as well. More detailed, when incident wavelength is fixed, the shift has to do with both lens size and the preset focal length, from Eq. (9). Meanwhile due to the fact that *p* values are always positive, so *N-p* curve lies above the horizontal axis. For the incident plane wave with a certain wavelength, *N* value of a planar lens can be calculated with the definite *a* and *f* by Eq. (1), and its relative focal shift can be obtained from Eq. (14). Substitute *p* value into Eq. (8), and theoretical prediction of focal length can be calculated. The relative focal shift will vary with *N* values. As a result, when the parameters of a lens change, the corresponding diffraction effect will exert different influence on their focusing performance. Although the shift can be decreased by increasing *N* values, constrained by the size of lens, boundary diffraction coming from *x* and *y* directions always exists, there is still 4.9% of shift when *N* = 14, greater than circular aperture [18] and nanoslit [24] situations, which also mean the squared planar lens yields a stronger diffraction effect for its sharp boundary on the four sides.

## 3. Design of metasurface lenses

In order to realize focusing function by metasurfaces, a group of suitable phase gradient units are needed. Here the C-shaped resonator is employed as a basic unit structure, which shown in Fig. 4(a).

The resonator is made from gold with the conductivity of 4.561 × 10^{7} S·m^{−1}, thickness along the *z*-axis is 0.1μm, which is placed on a silicon substrate (*ε*_{si} = 11.9) with thickness of 10μm, both of them composing a unit block. Parameter *α* is the half of placket angle, *β* is the angle between the *x*-axis and symmetrical line of resonator, which is 45°, *r* is the inner radius, *b* is the width of the ring, *P* is the period of each unit, here taking *b* = 0.1μm, *P* = 2μm. Aiming at the wavelength of 10.6μm, with the formula of resonance, equivalent radius of C-shaped resonator is estimated to be 0.8μm. Full wave simulations were performed for the chosen resonator by using commercial software CST Microwave Studio. Periodic boundary conditions were used in both *x* and *y* directions whereas perfectly matched layers boundary condition was applied for *z* direction. In the simulation, finite-difference time-domain (FDTD) method was selected and the *x*-polarized plane wave was used to excite a single resonator from the silicon side at normal incidence, besides, the resulting *y*-polarized wave under different geometric parameters were recorded. Then these eight resonators with different parameters of (*r*,2*α*) are chosen from calculated results at 28.3THz, which can modulate the phase of transmitted *y*-polarized wave at the gradient of π/8 from 0 to π, with amplitude transmittance remains constant. The geometrical parameters of resonators are *r* = 0.78, 0.75, 0.73, 0.74, 0.73, 0.75, 0.78, and 0.75μm corresponding 2*α* = 26°, 40°, 66°, 86°, 102°, 132°, 158°, and 176°, respectively. Each of these resonators is rotated by 90° clockwise around its circular centre, thus additional phase shift of π can be acquired respectively [27], with transmission amplitude nearly unchanged. The phase shifts and amplitudes of transmitted wave are shown in Fig. 4(b).

Based on the result in Fig. 4(b), the author is also try to testify the anomalous refraction characteristics of this set of resonators. When these sixteen resonators are excited under a normal incidence of linear polarized plane wave, the cross-polarized wave of transmitted will emerge anomalous refraction phenomenon in the theory. This set of resonators are conducted as a synthetic unit with the period of *P _{x}* = 32μm and

*P*= 2μm. The FDTD solutions was chosen for full wave simulations, periodic boundary conditions were used in both

_{y}*x*and

*y*directions whereas perfectly matched layers boundary condition was for

*z*direction. Figure 5(a)-5(c) show the scattered electric field distribution polarized in the

*y*direction which are excited by

*x*-polarized plane wave at normal incidence from the silicon substrate at the frequencies 27.8, 28.3 and 28.8THz, respectively. It can be seen that these three transmitted

*y*-polarized plane waves are well deflected, showing a broadband property. The deflected angles are 19.88°, 19.21°, and 19.12°, respectively, which all agree well with the theoretical results of 19.71°, 19.35°, and 19.01° from generalized Snell’s law [1], further demonstrating the broadband property of the metasurface. Meanwhile, the same results could also be obtained for the

*x*-polarized wave detection under a normal incidence of

*y*-polarized plane wave.

Suppose metasurface lens located in *xoy*-plane, as presented in Fig. 1, geometrical focus lies on the optical axis, the distance of which to the lens is acted as the preset focal length. Fresnel diffraction method is adopted to obtain the complex amplitude distribution on the metasurface lens. After discretization, the size of each pixel of lens is equal to the period of single resonator. The metasurface lens will be arranged by these designed resonators according to the calculated phase distribution, then the lens with certain theoretical focal length can be obtained.

To verify the theoretical analysis for focal shift above, at the wavelength of 10.6μm, several metasurface lenses with different parameters, corresponding to different Fresnel numbers, were arranged by the sixteen C-shaped resonators artificially, full wave simulations were also performed to observe their focusing and focal shift characteristics. When the side length 2*a* is fixed at 120μm, three lenses with different preset focal lengths of 100, 200, and 300μm were designed, respectively. For these *N* = 3.396, 1.698, and 1.132, with the values of *p* = 0.1773, 0.3012, and 0.3927 from Eq. (14), which all maintained 4-bit significant digits. Corresponding predicted focal lengths are 82.27, 139.76, and 182.19μm, respectively, from Eq. (9). The well-arranged lenses structure were simulated with FDTD method, open (and space) boundary conditions were applied in the *x* and *y* directions while perfectly matched layer boundary condition was for the *z* direction. The *x*-polarized plane wave was used to excited the lens structure normally from silicon side, electric field distribution of transmitted *y*-polarized wave were recorded, frequency of field monitor was set to 28.3THz. Simulation results of electric field distribution in the *yz*-plane are shown in Figs. 6-8.

Simulation results of the focal length for these three lenses are *f* ′ = 80, 136, and 179μm, which all agree well with theoretical results. Similarly, for the case that when preset focal length *f* is fixed at 100μm, lenses with the side length 2*a* = 100 and 140μm, were also designed, for *N* = 2.358 and 4.623, corresponding *p* = 0.2369 and 0.1367, predicted focal lengths are 76.31 and 86.33μm from Eq. (9), results of simulation are *f* ′ = 75 and 82.5μm, as shown in Figs. 9-10, which are also in good accordance with predicted values.

In order to revalidate the previous analysis, method of calculating focal shift provided above could equally be applied to designing the plasmonic lens with specific focal length. With the incident wavelength *λ* = 10.6μm, side length 2*a* = 120μm, assuming that a lens with the real focal length of 100μm is desired, thus Eqs. (1), (8), and (14) are being used to calculate its preset focal length and the result is *f* = 127.48μm. Then the lens was arranged and simulation result of field intensity distribution in *yz*-plane is shown in Fig. 11.

As can be seen from Fig. 11, simulation result presents a strong focusing performance with the focal length *f* ′ = 98.5μm, which is almost coincidence with our expectation, verified the validity of early analysis again.

From above, it is worth noting that there exist a little deviations between the theoretical prediction and simulation results, while the values of simulated are less than the predicted slightly. However, since the amplitude transmission remains identical for all C-shaped resonators on the metasurface lens, so the overall transmission of the lens won’t influence its axial field distribution. Furthermore, the deviations may be caused by the approximation, which is introduced from Eq. (5) to Eq. (6), for the expression of spherical wave is replaced by the quadric surface wave, leading to a longer focal length in the process of theoretical calculation. Besides, the relative deviations will decrease when the preset focal length is increased or the side length is reduced, as shown in Figs. 6-8 along with Figs. 9 and 10, due to they are better satisfied the condition of *a* << *f*.

## 4. Discussion

Several lenses designed above are all simulated at the wavelength of 10.6μm, considering the application in CO_{2} laser. As discussed in the Section 3, a group of C-shaped resonators were designed at the central frequency of 28.3THz, and they have also exhibited similar response at 27.8-28.8THz from the simulation results, showing a broadband property. As a result, it can be stated that these lenses are frequency scalable. If an researcher wants to investigate their focal shift performance, he must calculate its Fresnel number according to Eq. (1), and the shift will be further obtained by substitute the *N* value into Eq. (14).

Next, another important parameter is diffraction efficiency. In general, the efficiency of a diffractive optical element is defined as the ratio of power on the focal plane to the total power of the incident beam. For the case of single point focusing lenses designed in this paper, the author first integrate the intensity of the electric field on the focal plane over the round area centered in the peak of the focus, and then the ratio of its value to the total incident power can be calculated. Here diameter of focal spot is defined by the full-width at half-maximum (FWHM). The simulation results of intensity distribution of the several lenses on each focal plane are shown in Fig. 12(a)-12(f), both of them are normalized to their focal spot. The corresponding calculation results of diffraction efficiency are all listed in Table 1.

From the table it can be found that the diffractive efficiency of these lenses is about 2% at the frequency of 28.3THz. The low efficiency is mainly caused by the lower conversion efficiency from one polarization state to its cross-polarized (*x*-polarized to *y*-polarized). The theoretical limit of the conversion efficiency under the incidence of linearly polarized wave is only 25% for the transmissive type of metasurfaces. Besides, the loss of silicon substrate also contributes to decreasing the diffractive efficiency. In the simulation, a silicon substrate with the thickness of 10μm was used, and if the effect of the substrate is not considered, the theoretical diffraction efficiency of the single-layer metasurface lenses designed in this paper can reach to 5%. At the same time, it can also be concluded from the table that, the focal shift caused by the diffraction effect of the lens edge will result in some of power losses when focusing. For example, when the side length of lens is fixed at 120μm, with the increase of preset focal length, the approximate condition of *a << f* will be better satisfied. Therefore, the focal shift caused by diffraction is more obvious. The efficiency is gradually reduced from 2.79% to 1.91%, this also suggests that when the shift increases, the short focal length lens would be more efficient than the long. The low efficiency of the long focal length lens may be caused by the diffraction of the lens edge, leading to that a part of incident power is dispersed into the other free spaces. Similarly, when the preset focal length is fixed at 100μm, the area of the lens increases with the side length, the focal shift is getting smaller, and the efficiency also increases slightly, which also prove that the efficiency will decrease as the shift increase. And this circumstance is also consistent with the traditional lens, as the lens size becomes larger or the focal length becomes smaller, it is expressed as the greater numerical aperture (NA), then the better focusing performance of the lens is.

It is generally known that the focal shift phenomenon can be observed obviously only when the diffraction effect of lens edge is significant enough to influence the focusing performance. As illustrated in Eq. (9), the focal shift of lens only depends on the Fresnel number. However, the Fresnel number is determined by the side length, preset focal length of lens, and the wavelength of incident wave. In fact, the focal shift caused by the diffraction always exists for the lens. It can be seen from Eq. (14) and Fig. 3, when the incident wavelength is fixed, with the increase of Fresnel number, relative focal shift decreases rapidly, and the real focal length will be closer to preset value, thus the limit of shift tends to be zero. However, the prediction error will increase, since it is unable to meet the approximate condition of *a* << *f* very well. When the Fresnel number decreases, letting preset focal length *f* reach to infinity in Eq. (6), then the axial field distribution of light intensity of Eq. (9) is converted to:

where the upper limit of the two Fresnel integrals *w* is:

Under this circumstance, preset focal length of the lens is infinity, the Fresnel number *N =* 0, phase on the planar lens are all equal, the intensity distribution along the axis completely depends on the diffraction of lens edge toward the incident light. For example, when the side length of lens 2*a* = 120μm, incident wavelength *λ* = 10.6μm, and the *w* = 1.21, from Eq. (16), the position of maximum intensity along the axis is at a distant about 463.93µm. This also means that when the side length and wavelength are taken as the above values respectively, the maximal focal length of the metasurface lens designed by traditional diffraction method will not go over this maximum distance. Because the preset focal length of planar lens be seen as infinity at this moment, therefore it can also be considered that the focal shift of the lens is infinity as well. Thus, it can be concluded that the limits on the variation of focal shift of the designed lens is from zero to infinity.

In the process of sample fabrication on designed metasurface lens, the C-shaped gold nanostructures for the near-field focusing experiments can be fabricated on a 10-μm-thick silicon substrate by using electron beam lithography into a PMMA (polymethylmethacrylate) resist film. And a layer of 0.1μm gold film will be deposited by electron-beam evaporation in a highly vacuum chamber, then the lens sample is obtained with a following standard lift-off procedure to remove the rest of resist. Incident infrared light is coming from a grating tunable CO_{2} laser which can generate the light wave with a pure single frequency in a broadband range from 27.8 to 28.8THz, and collimated by the parabolic mirror to a plane wave [Fig. 13]. Besides, a linear polarizer should be placed in front of the lens sample to ensure only the *x*-polarized wave excite the structure from the silicon side while another polarizer also placed behind the sample to allow *y*-polarized light output. The HgCgTe infrared camera (FLIR SC7900VL) can be employed to characterize the focal spot of the lens sample in the focal plane. Moreover, a microscope objective can also be placed behind the polarizer to collected transmitted light for clear imaging, and one can adjust the position of camera along the optical axis to ensure the focal spot be imaged by camera clearly.

## 5. Conclusion

In summary, the focal shift issues of squared planar lenses based on metasurfaces are discussed theoretically and demonstrated by simulation. By adopting traditional diffraction theories, axial light intensity distribution is deduced with Fresnel approximation. Compared with previous research by taking the derivative of axial field distribution or by series expansion method, the Fresnel integral and Cornu spiral here are utilized to obtain an exact relationship between focal shift and Fresnel number of the lens, which provides a more convenient and intuitive way to analyze the focal shift issues caused by the diffraction effect. Meanwhile, reducing preset focal length or increasing the size of lens all contributes to decreasing the focal shift. A group of C-shaped phase gradient resonators are designed with phase shifts ranged from 0 to 2π in the far infrared region, which can be arranged to realize a metasurface lens, meanwhile, the broadband property is also demonstrated. Several lenses with different Fresnel numbers are also designed at the wavelength of 10.6μm and the simulation results all accord well with theoretical expectation. Small deviations between them may caused by the approximation which introduced from the theoretical analysis. Besides, as demonstrated above, the predicting formula owns a relatively high precision under the condition of (*a* / *f*)^{2} < 0.5. The diffraction efficiency of designed lenses is around 2% at 28.3THz from simulation results. Moreover, efficiency will decrease as the shift increase, which also indicate that the short focal length lens is more efficient than the long. Proposed method to predict focal shift in this paper may provide a reference for designing plasmonic lens in the future, for example, multi-focus lens which focusing in the focal plane or along the axis, or to compensate the chromatic aberration.

## Funding

National Natural Science Foundation of China (NSFC) (Grant Nos. 11274389, 61331005).

## Acknowledgments

The author would like to acknowledge the editor and two anonymous reviewers for their valuable advice to improve the quality of this manuscript. The author acknowledges Dr. Guangyao Shi at the University of Maryland for his help in the polish of the language. The funding is supported by the Metamaterials Laboratory, Department of Basic Education at the Air Force Engineering University, China.

## References and links

**1. **N. Yu, P. Genevet, M. A. Kats, F. Aieta, J. P. Tetienne, F. Capasso, and Z. Gaburro, “Light propagation with phase discontinuities: Generalized laws of reflection and refraction,” Science **334**(6054), 333–337 (2011). [CrossRef] [PubMed]

**2. **P. Genevet, N. Yu, F. Aieta, J. Lin, M. A. Kats, R. Blanchard, M. O. Scully, Z. Gaburro, and F. Capasso, “Ultra-thin plasmonic optical vortex plate based on phase discontinuities,” Appl. Phys. Lett. **100**, 013101 (2012). [CrossRef]

**3. **X. Ni, A. V. Kildishev, and V. M. Shalaev, “Metasurface holograms for visible light,” Nat. Commun. **4**, 2807 (2013). [CrossRef]

**4. **L. Huang, X. Chen, H. Mühlenbernd, H. Zhang, S. Chen, B. Bai, Q. Tan, G. Jin, K. W. Cheah, C. W. Qiu, J. Li, T. Zentgraf, and S. Zhang, “Three-dimensional optical holography using a plasmonic metasurface,” Nat. Commun. **4**, 2808 (2013). [CrossRef]

**5. **Q. Wang, X. Zhang, Y. Xu, J. Gu, Y. Li, Z. Tian, R. Singh, S. Zhang, J. Han, and W. Zhang, “Broadband metasurface holograms: toward complete phase and amplitude engineering,” Sci. Rep. **6**, 32867 (2016). [CrossRef] [PubMed]

**6. **J. Lin, J. P. B. Mueller, Q. Wang, G. Yuan, N. Antoniou, X. C. Yuan, and F. Capasso, “Polarization-controlled tunable directional coupling of surface plasmon polaritons,” Science **340**(6130), 331–334 (2013). [CrossRef] [PubMed]

**7. **L. Huang, X. Chen, B. Bai, Q. Tan, G. Jin, T. Zentgraf, and S. Zhang, “Helicity dependent directional surface plasmon polariton excitation using a metasurface with interfacial phase discontinuity,” Light Sci. Appl. **2**, e70 (2013). [CrossRef]

**8. **F. Aieta, P. Genevet, M. A. Kats, N. Yu, R. Blanchard, Z. Gaburro, and F. Capasso, “Aberration-free ultrathin flat lenses and axicons at telecom wavelengths based on plasmonic metasurfaces,” Nano Lett. **12**(9), 4932–4936 (2012). [CrossRef] [PubMed]

**9. **X. Li, S. Xiao, B. Cai, Q. He, T. J. Cui, and L. Zhou, “Flat metasurfaces to focus electromagnetic waves in reflection geometry,” Opt. Lett. **37**(23), 4940–4942 (2012). [CrossRef] [PubMed]

**10. **D. Hu, X. Wang, S. Feng, J. Ye, W. Sun, Q. Kan, P. J. Klar, and Y. Zhang, “Ultrathin terahertz planar elements,” Adv. Opt. Mater. **1**(2), 186–191 (2013). [CrossRef]

**11. **X. Y. Jiang, J. S. Ye, J. W. He, X. K. Wang, D. Hu, S. F. Feng, Q. Kan, and Y. Zhang, “An ultrathin terahertz lens with axial long focal depth based on metasurfaces,” Opt. Express **21**(24), 30030–30038 (2013). [CrossRef] [PubMed]

**12. **X. Chen, Y. Zhang, L. Huang, and S. Zhang, “Ultrathin metasurface laser beam shaper,” Adv. Opt. Mater. **2**(10), 978–982 (2014). [CrossRef]

**13. **Q. Wang, X. Zhang, Y. Xu, Z. Tian, J. Gu, W. Yue, S. Zhang, J. Han, and W. Zhang, “A broadband metasurface-based terahertz flat-lens array,” Adv. Opt. Mater. **3**(6), 779–785 (2015). [CrossRef]

**14. **J. He, J. Ye, X. Wang, Q. Kan, and Y. Zhang, “A broadband terahertz ultrathin multi-focus lens,” Sci. Rep. **6**, 28800 (2016). [CrossRef] [PubMed]

**15. **D. Jia, Y. Tian, W. Ma, X. Gong, J. Yu, G. Zhao, and X. Yu, “Transmissive terahertz metalens with full phase control based on a dielectric metasurface,” Opt. Lett. **42**(21), 4494–4497 (2017). [CrossRef] [PubMed]

**16. **Q. Yang, J. Gu, D. Wang, X. Zhang, Z. Tian, C. Ouyang, R. Singh, J. Han, and W. Zhang, “Efficient flat metasurface lens for terahertz imaging,” Opt. Express **22**(21), 25931–25939 (2014). [CrossRef] [PubMed]

**17. **L. Cong, N. Xu, J. Han, W. Zhang, and R. Singh, “A Tunable Dispersion-Free Terahertz Metadevice with Pancharatnam-Berry-Phase-Enabled Modulation and Polarization Control,” Adv. Mater. **27**(42), 6630–6636 (2015). [CrossRef] [PubMed]

**18. **Y. Li and E. Wolf, “Focal shifts in diffracted converging spherical waves,” Opt. Commun. **39**(4), 211–215 (1981). [CrossRef]

**19. **Y. Li, “Dependence of the focal shift on Fresnel number and *f* number,” J. Opt. Soc. Am. **72**(6), 770–774 (1982). [CrossRef]

**20. **Y. Li and H. Platzer, “An experimental investigation of diffraction patterns in low-Fresnel-number focusing systems,” Opt. Acta (Lond.) **30**(11), 1621–1643 (1983). [CrossRef]

**21. **Y. Li and E. Wolf, “Three-dimensional intensity distribution near the focus in systems of different Fresnel numbers,” J. Opt. Soc. Am. A **1**(8), 801–808 (1984). [CrossRef]

**22. **Y. Yu and H. Zappe, “Effect of lens size on the focusing performance of plasmonic lenses and suggestions for the design,” Opt. Express **19**(10), 9434–9444 (2011). [CrossRef] [PubMed]

**23. **Y. Yu and H. Zappe, “Theory and implementation of focal shift of plasmonic lenses,” Opt. Lett. **37**(9), 1592–1594 (2012). [CrossRef] [PubMed]

**24. **Y. Gao, J. Liu, X. Zhang, Y. Wang, Y. Song, S. Liu, and Y. Zhang, “Analysis of focal-shift effect in planar metallic nanoslit lenses,” Opt. Express **20**(2), 1320–1329 (2012). [CrossRef] [PubMed]

**25. **M. Born and E. Wolf, *Principles of Optics*, 7th (expanded) ed. (Cambridge, 1999).

**26. **E. Jahnke and F. Emde, *Tables of Functions with Formulae and Curves*, 4th ed. (Dover, 1945).

**27. **L. Liu, X. Zhang, M. Kenney, X. Su, N. Xu, C. Ouyang, Y. Shi, J. Han, W. Zhang, and S. Zhang, “Broadband metasurfaces with simultaneous control of phase and amplitude,” Adv. Mater. **26**(29), 5031–5036 (2014). [CrossRef] [PubMed]