**Dipolar SP excitation under CP incidences.** In general, a dipole source consists two monopole sources of equal strength but opposite phase and separated by a distance far smaller than wavelength. Suppose an SP dipole source located at *xy*-plane, as depicted in Fig. 1a, rotated by an angle of *θ* with respect to the *x*-axis. This dipole can only be excited by incidence component polarized along its orientation direction and then produce SPs ripple around. In the case of a CP normally incident from the back side with \(\sigma \in \left\{ {+, - } \right\}\) standing for the spin direction LCP or RCP, the SP field at an arbitrary point *P* can be calculated as [42]:

$${E_P}=\eta \frac{{\sqrt 2 }}{2}\frac{{\exp \left( {{\text{i}}{k_{{\text{SP}}}}\left| {\mathbf{r}} \right|} \right)}}{{{\text{i}}\sqrt {{\lambda _{{\text{SP}}}}\left| {\mathbf{r}} \right|} }}\left( {\cos \theta +{\text{i}}\sigma \sin \theta } \right)\cos \left( {\theta - \xi } \right)$$

1

Here, *η* is the complex coefficient that describes the conversion efficiency of SP excitation from Ein and the phase shift acquired from the dipole resonance; *λ*SP is the SP wavelength; *k*SP = 2π/*λ*SP is the SP wave number; r is the vector from the dipole position to the point *P*; *ξ* is the azimuth angle of point *P*, i.e., the angle of r with respect to the *x*-axis. By performing trigonometric transformation, Eq. (1) turns to:

$${E_P}=\eta \frac{{\sqrt 2 }}{2}\frac{{\exp \left( {{\text{i}}{k_{{\text{SP}}}}\left| {\mathbf{r}} \right|} \right)}}{{{\text{i}}\sqrt {{\lambda _{{\text{SP}}}}\left| {\mathbf{r}} \right|} }}\left[ {\frac{1}{2}{e^{{\text{i}}\sigma \left( {2\theta - \xi } \right)}}+\frac{1}{2}{e^{{\text{i}}\sigma \xi }}} \right]$$

2

Besides the propagation evolution, the initial excitation can be broken down into two parts: exp[i*σ*(2*θ-ξ*)]/2 and exp(i*σξ*)/2, named as term 1 and term 2. Intriguingly, the excitation phases of these two terms are all related to the azimuth angle *ξ* with the spin dependence of *σ*. Figure 1b illustrates the calculated SP field distributions emitted by a dipole source with orientation angle of *θ = +*45°, by using Eq. (1). Correspondingly, Figs. 1c and d illustrate the excitation component results from term 1 and term 2, respectively, which are two spiral-like sources with opposite swirl directions. In addition, term 1 is also related to *θ* with spin dependence of 2*σ*, which means that the initial excitation phase of term 1 for arbitrary target direction can be freely controlled by changing *θ*. In contrast, the SP excitation of term 2 is independent with *θ*. For instance, Figs. 1e-g illustrate the calculated SP field distributions in the case of *θ* = − 45°. Clearly, the SP excitation in Fig. 1f has a phase difference of π comparing with that in Fig. 1c, while the SP excitation in Fig. 1g is the same as that in Fig. 1d. The opposite spin will also reveal this phenomenon due to the *σ* only changing the swirl direction in both term1 and term 2 (see Fig. S1, Supporting Information).

For a set of dipole sources located at given positions, the overall SP excitation results from term 1 can be joint steered to form target focusing beams by manipulating their orientation angles, while the SP excitation results from term 2 are fixed and uncontrollable. Intriguingly, the latter can be significantly suppressed by arranging the dipole sources as a two-dimensional periodic array with a deep subwavelength period, where the SP excitations from different dipole sources eliminate each other due to the destructive interference. For reference, we compared the SP excitation results from term 1 and term 2 of a single focus metalens, where the latter is at least 2 orders of magnitude weaker than the former (see Section B and Fig. S2, Supporting Information). As thus, we adopt subwavelength lattice arrangement and ignore the influence of term 2 in the design of following multiplexed metalenses.

**Holographic approach for multiplexed metalens design.** Computer-generated holography, widely employed for tailoring the free-space waves with desirable field distributions, has become a handy tool to realize complex functionalities in many areas [43–49]. Here, we utilize a two-dimensional holographic approach to design the multi-channel SP metalenses. In simplicity, the basic principle of holography is the reversibility of the light path. As schematically shown in Fig. 1h, suppose there are *m* dipole sources, inside a circular region of radius *R*, composing a hologram to record the scattering information from all imaginary foci. The dipole sources are arranged in a square lattice with the subwavelength period of Λ along both *x* and *y* directions. In general, the imaginary foci can be arbitrarily placed. Suppose there are *n* imaginary foci respectively working for the LCP incidences of *n* wavelengths, as depicted by the yellow points in Fig. 1h, where each imaginary focus can be treated as a point SP source of the corresponding SP wavelength *λ**n*. The superpositions at the *m*th dipole source can be calculated as:

$$H_{m}^{+}=\sum\limits_{n} {\frac{{S_{n}^{+}\exp \left[ { - {\text{i}}{k_n}\left| {{{\mathbf{r}}_{nm}}} \right|+{\text{i}}\arg \left( {{{\mathbf{r}}_{nm}}} \right)} \right]}}{{{\text{i}}\sqrt {{\lambda _n}\left| {{{\mathbf{r}}_{nm}}} \right|} }}}$$

3

here is the complex coefficient of the *n*th imaginary point source; *k**n* = 2π/*λ**n* is the SP wave number; r*nm* is the vector from the *n*th imaginary point source to the *m*th dipole source. Similarly, suppose there are *n’* imaginary foci respectively working for the RCP incidences of *n’* wavelengths, as depicted by the green points in Fig. 1h, the superpositions at the *m*th dipole source can be calculated as:

$$H_{m}^{ - }=\sum\limits_{{n'}} {\frac{{S_{{n'}}^{ - }\exp \left[ { - {\text{i}}{k_{n'}}\left| {{{\mathbf{r}}_{n'm}}} \right| - {\text{i}}\arg \left( {{{\mathbf{r}}_{n'm}}} \right)} \right]}}{{{\text{i}}\sqrt {{\lambda _{n'}}\left| {{{\mathbf{r}}_{n'm}}} \right|} }}}$$

4

here and *k**n’* is the complex coefficient and SP wave number for the *n’*th imaginary point source, and the r*n’m* is the vector from the *n’*th imaginary point source to the *m*th dipole source. Once the superpositions from all the imaginary foci are obtained, the orientation angle of the *m*th dipole source can be directly calculated as:

$${\theta _m}={{\arg \left[ {H_{m}^{+}+{{\left( {H_{m}^{ - }} \right)}^*}} \right]} \mathord{\left/ {\vphantom {{\arg \left[ {H_{m}^{+}+{{\left( {H_{m}^{ - }} \right)}^*}} \right]} 2}} \right. \kern-0pt} 2}$$

5

here the superscript * stands for conjugate operation. Once all the dipole sources are rotated by *θ**m*, as schematically shown in Fig. 1i, the multiplexed metalens design is finished. The superposed SP field at an arbitrary point *P**l* (see Fig. 1i) from all the dipole sources, under a specific CP incidence (wavelength of *λ*SP and spin direction of *σ*) can be calculated as:

$${E_l}=\eta \frac{{\sqrt 2 }}{2}\sum\limits_{m} {\frac{{\left( {\cos {\theta _m}+{\text{i}}\sigma \sin {\theta _m}} \right)\exp \left( {{\text{i}}{k_{{\text{SP}}}}\left| {{{\mathbf{r}}_{ml}}} \right|} \right)\cos \left[ {{\theta _m} - \arg \left( {{{\mathbf{r}}_{ml}}} \right)} \right]}}{{{\text{i}}\sqrt {{\lambda _{{\text{SP}}}}\left| {{{\mathbf{r}}_{ml}}} \right|} }}}$$

6

here r*ml* is the vector from the center of the *m*th dipole source to the point *l*. By calculating the superpositions over a region of interest, the SP excitation and focusing performance for different incidences can thus be obtained. Notable, the complex coefficient *η* is wavelength-dependent, which is determined by the realistic resonator design.

**Realistic resonator design.** In order to implement the above design into realistic metalenses, it needs to seek a proper resonator that has the broadband dipole response in SP excitation (see Fig. S3, Supporting Information). The connected-double-ring slit resonator, as depicted by Fig. 2a, was chosen as the realistic resonator. The resonator is made from 200 nm thick Aluminum film patterned on a 2-mm-thick high-resistance silicon wafer. Figure 2b illustrates the simulated electric field distributions of a single connected-double-ring slit resonator under the LCP and RCP incidences at wavelength of 375 µm, 450 µm, 525 µm, and 600 µm, respectively. Clearly, the emitted SP fields under these different CP incidences all exhibit dipolar responses, that is, of equal strength but opposite phase. As thus, we adopt this resonator as the element unit in the following metalens designs.

**Metalens I.** Figure 3a illustrates the schematic view of an eight-channel metalens design, named as Metalens I. The focusing positions are designated at a circle with the radius of 6 mm, where the LCP incidences at working wavelengths of 600 µm, 525 µm, 450 µm, and 375 µm can be coupled into focusing beams along azimuth angles of 5π/8, 7π/8, − 7π/8, and − 5π/8, respectively. The focusing positions for RCP incidences are opposite to that of LCP cases. The connected-double-ring slit resonators are arranged in a square lattice with period of 100 µm along both *x* and *y* directions, inside a circular region with the radius of 3 mm. The center of this circular region is the origin of the in-plane *xy* coordinate system. For simplicity, each imaginary point source *S* + *n* or *S– n'* is assigned to 1. In this condition, the orientation angles of all the resonators can be calculated by using Eqs. (3)-(5), as shown in Fig. 3b. We fabricated Metalens I by using conventional photolithography and metallization processes (see Methods for the details). The microscope image in Fig. 3c shows the excitation part of the fabricated Metalens I. By using Eq. (6), the focusing performances for different incidences were calculated and illustrated in the Fig. 3d. Clearly, eight high-quality focusing beams can be respectively formed along the target directions, numerically verifying our proposed holographic approach. Furthermore, the fabricated Metalens I was characterized by a near-field scanning terahertz microscopy (see Methods for the details) [22], and corresponding measurements are illustrated in Fig. 3e. Obviously, the SP excitations and propagation features agree very well with the calculation results, which confirms the broadband in-plane dipole functionality of our selected resonator and the feasibility of our proposed design strategy. Further extracted from the experimental SP profiles, the intensities round on the radius of 6 mm in Fig. 3e can be plotted to show the crosstalk quantitative analysis at each expected channel in the polar coordinates (see Fig. S6a, Supporting Information), where each of the focusing beams exhibits good quality and has very low crosstalk in terms of both directionality and wavelength isolation.

**Metalens II.** In addition to steering the different incidences into different on-chip directions, in some application scenarios, it needs to selectively couple more than one kind of incidence into SP focusing beams along the same direction, performing demultiplexing and multiplexing operations simultaneously. Figure 4a illustrates the schematic view of another eight-channel metalens design, named as Metalens II. In this design, the focusing positions are also designated at a circle with the radius of 6 mm. Each two of the LCP incidences can be coupled into focusing beams along the same direction: 600 µm and 450 µm toward the azimuth angle of 3π/4, 525 µm and 375 µm toward the azimuth angle of − 3π/4. The focusing positions for RCP incidences are opposite to that of LCP cases. Figures 4b and c show the calculated resonator orientation angle distributions and the microscopic image of the fabricated Metalens II. The calculated and measured results of Metalens II are shown in Figs. 4d and e, which agree well with each other that eight high quality SP focusing beams can be formed as expected. Similarly extracted from the corresponding experimental SP profiles round on the radius of 6 mm in the Fig. 4e, Figure S6b indicates that each two selected incidences can be coupled into SP focusing beams along the same direction, exhibiting the superior performance of our proposed design scheme.

To corroborate the generality of our proposed design scheme, two additional eight-channel metalenses were fabricated and characterized: Metalens III and Metalens IV (see Fig. S4 and Fig. S5, Supporting Information). Metalens III can couple the eight incidences into SP focusing beams along two opposite directions according to the polarization of LCP and RCP, which functions as an on-chip polarization spitter and wavelength multiplexer simultaneously (see Fig. S4a). Metalens IV can couple all the eight kinds of incidences into SP focusing beams along the same direction, which functions as an eight-channel on-chip multiplexer (see Fig. S5a). Figures S4b-e and Figures S5b-e illustrate the corresponding calculation and measurement results of Metalens III and Metalens IV, respectively, which agree well with each other that can satisfactorily perform predesignated functionalities. The corresponding quantitative analyses are likewise extracted in the Figs. S6c and d, which also perform expected directionality for different kinds of outgoing beams with less crosstalk. Despite we only demonstrate eight-channel metalenses in the experiment, the selected resonator and the proposed design scheme also hold promise for even more multiplexing channels. For instance, we carried out several simulation studies with the identical metasurface configuration to show the capability of achieving ten-channel focusing of SPs (see Fig. S7, Supporting Information), while the constant capacity of multiplexing energy only allows each channel energy to be divided with the multiplexing channels increasing on numerical results (see Fig. S8, Supporting Information). Moreover, the focusing positions for different multiplexing channels can be freely arranged, and the focusing directions for different spins at the same wavelength are not limited to be opposite (see Fig. S9, Supporting Information), where high-quality SP focusing beams are obtained in each calculation result for the metalenses of different focusing position configurations as expected. These results further verify the robustness and versatility of our proposed design strategy.