### Principle

The PM-based DIM scheme is depicted in Fig. 1. It consists of four parts, namely, multi-channel bit streams, a microcontroller unit (MCU), a PM-based transmitter, and user equipment (UE) located in different directions. The scheme directly transmits the digital information, and the mechanism is presented as follows. Firstly, the phase coding library of the PM is determined by an efficient discrete optimization algorithm, and the aim is to generate the symbols of the specified modulation scheme in the desired directions. Then, the bit streams are translated into symbols, which are further mapped into the phase coding of PM by indexing the library. The phase coding is then imposed on the tunable components of PM via MCU. Finally, UEs extract the information bits by comparing the received signals with the corresponding constellation symbols. The desired UEs receive the regular constellation diagrams, which have a phase shift due to the EM wave propagation and can be compensated by adding an identical phase bias for each received signal. On the contrary, the undesired UEs will receive distorted constellation diagrams. Since the structure of the received signals is disturbed, the correct information acquisition is nearly impossible, regardless of the receiver sensitivity.

We briefly present the element of the PM utilized in this work. The geometric size is 13 × 13 × 2.5 mm (0.48 × 0.48 × 0.09 λ at 11 GHz), exhibiting the low-profile feature which has not been reported in the previous DIM schemes based on the transmissive or reflective metasurfaces. According to the principle of the DIM scheme, the hardware can work in the transmitting or receiving mode (see Supplementary Note 3 for more details). The main difference between the transmitting and receiving modes is the coding pattern. The former is updated in real-time with the input constellation symbol, while the latter is fixed to ensure a high signal-to-noise ratio. It is worth noting that the receiving capability of the design in the current work was established by two important modifications to the PM. Firstly, we improved the feed network and used a Wilkinson power divider with high isolation. In the previous design, the feed network employed a T-type power divider, which made the received signals distorted due to serious interference between different channels. Secondly, we replaced the varactor with the PIN diodes at the feed network. This is necessary due to the limitations of the varactor, which requires high tunable voltages, complex control circuits, and large power supplies. Each element comprises three parts, including the radiation layer, the direct current bias layer, and the feed network layer, as depicted in Fig. 2a. The three parts are laminated together with two pieces of prepreg substrates. In the top radiation layer, dual PIN diodes bridge two rectangular metal patches and a small square patch. The direct current bias layer is located between the two grounding plates, creating an EM-shielded space and thus reducing the influence of the two bias lines. One bias line with the scalloped branch connects to the radiation layer, and the other links to the feed network layer. The feed network layer consists of a two-port reflection-type phase shifter that integrates two PIN diodes. The 2-bit phases (S1: − 90°, S2: 0°, S3: 90°, S4: 180°) are achieved by adjusting the voltages imposed on the PIN diodes. It is worth noting that the implementation of 2-bit phases does not need external phase shifters, thereby reducing the complexity and cost.

Figure 2b demonstrates the magnitude and phase responses of the meta-element. The magnitude responses of the four states are close to one in the wide frequency range of 10.7 ~ 11.7 GHz, and the phase difference of the four states remains almost 90° over the whole frequency band, exhibiting the broadband performance of the meta-element. More geometric parameters and operating principles of the element are detailed in our previous work^{31}.

### Algorithm

As the above-mentioned model, we consider the situation where a PM-based transmitter, located at the xoy plane and loaded with *N* tunable elements, simultaneously serves *K* single-antenna users by transmitting constellation symbols to each of them. The reference constellation sets \({{{\rm{O}}}}\), including 8PSK, 16QAM, and 64QAM, are encoded with Gray code to reduce the error bit rates (see Supplementary Note 4). The received signal at the *k*th user is

$${y}_{k}={{{{\bf{h}}}}}_{k}^{H}{{{\bf{x}}}}+n,\, k=1,\,{{\mathrm{..}}}.,\,K,$$

(1)

where \({{{\bf{x}}}}\in {{{{\rm{C}}}}}^{N},\, n \sim {{{\rm{C}}}}{{{\rm{N}}}}(0,\, {\sigma }^{2})\), and *H* are the coding sequences, the complex additive white Gaussian noise, and the conjugate transpose operator, respectively; and \({{{{\bf{h}}}}}_{k}^{H}\in {{{{\rm{C}}}}}^{N}\) is the channel vector

$${{{{\bf{h}}}}}_{{{{\bf{k}}}}}={e}^{-j\beta {d}_{k}}\sqrt{G({\theta }_{k},\, {\varphi }_{k})}{\left({e}\, ^{j{{{{\bf{v}}}}}_{1,1}^{H}{{{\bf{u}}}}},{{\mathrm{..}}}.,{e}\, ^{j{{{{\bf{v}}}}}_{p,q}^{H}{{{\bf{u}}}}},{{\mathrm{..}}}.,{e}^{j{{{{\bf{v}}}}}_{{N}_{x},{N}_{y}}^{H}{{{\bf{u}}}}}\right)}^{T},$$

(2)

where *d*_{k}, *θ*_{k}, *φ*_{k}, and \(G({\theta }_{k},\, {\varphi }_{k})\) are the distance, the elevation angle, the azimuth angle, and the directional gain of the element with respect to the *k*th user, respectively. In addition,\({{{{\bf{v}}}}}_{p,q}={(p{d}_{x},\, q{d}_{y})}^{T}\) and \({{{\bf{u}}}}=2\pi /\lambda {(\sin {\theta }_{k}\cos {\varphi }_{k},\sin {\theta }_{k}\sin {\varphi }_{k})}^{T}\) are the auxiliary vectors to simplify Eq. (2). We clarify that although each meta-element can only provide finite possible discrete phases of \({x}_{i}\), their propagation distances to the user are different. As the scale of array increases, there will be many \({h}_{i}\) with different phases, and the role of \({x}_{i}\) is to combine those \({h}_{i}\) to provide the desired symbol at the *k*th user. We stress that this principle has been analyzed and demonstrated in ref. ^{32}.

Our goal is to make the received signals \({{{\bf{y}}}}={({y}_{1},{{\mathrm{..}}}.,{y}_{K})}^{T}\) as close to the reference constellation symbols \({{{\bf{s}}}}={({s}_{1},{{\mathrm{..}}}.,{s}_{K})}^{T}\in {{{\rm{O}}}}\) as possible. To this end, the coding sequences **x** of the PM are elaborately optimized to minimize signal variation. In this context, we formulate the optimization problem as

$$\begin{array}{cc}{\min }_{{{{\bf{x}}}},\beta } & {\Vert {{{\bf{s}}}}-\beta {{{\bf{Hx}}}}\Vert }_{2}^{2}+K{\beta }^{2}{\sigma }^{2}\\ {{{\rm{s}}}}.\,{{{\rm{t}}}}. & x(i)\in {{X}},i=1,\,{{\mathrm{..}}}.,\,N\end{array},$$

(3)

where \({{{\bf{H}}}}=({{{{\bf{h}}}}}_{1}^{H},{{{{\bf{h}}}}}_{1}^{H},\,{{\mathrm{..}}}.,\,{{{{\bf{h}}}}}_{K}^{H})\in {{{C}}}^{K\times N}\) and \({{X}}=\{1/\sqrt{N}{e}^{j{w}_{i}}|{w}_{i}=2\pi \cdot {{{\rm{i}}}}/{2}^{P},i=1,2,\,{{\mathrm{..}}}.,\,{2}^{P}\}\) are the channel matrix and the *P*-bit phase quantization set, respectively, and \(\beta \in {{\mathbb{R}}}^{+}\) is the precoding factor to rescale the magnitude of the received signal *y*_{k} to the same level as the constellation symbol *s*_{k}. In this work, we assume that the *β* factor transmitted by the pilot technology is perfectly known to all users.

Problem (3) involves the widely recognized issue of finite resolution precoding in massive multi-input multi-output (MIMO) systems^{33,34}. The problem is non-convex due to the discrete constraint. The problem can be firstly converted into a mixed integer linear program by using the special ordered set technique and then solved by using the branch-and-bound algorithm^{35}. However, the computational complexity is unaffordable for large-scale arrays. In this study, we employ the ADMM framework to address the problem. The framework adopts the idea of divide and conquer that effectively decomposes a non-convex problem into several convex problems and exhibits good performance in terms of fast convergence and stability^{36,37}. By introducing an auxiliary variable, we rewrite the problem as

$$\begin{array}{cc}{\min }_{{{{\bf{x}}}},\tilde{{{{\bf{x}}}}},\beta \,} & {\Vert {{{\bf{s}}}}-\beta {{{\bf{H}}}}{{{\bf{x}}}}\Vert }_{2}^{2}+K{\beta }^{2}{\sigma }^{2}\\ {{{\rm{s}}}}.\,{{{\rm{t}}}}. & \begin{array}{c}\tilde{{{{\bf{x}}}}}={{{\bf{x}}}},\\ \tilde{{x}}(i)\in {{X}},i=1,\,{{\mathrm{..}}}.,\,N\end{array}\end{array}.$$

(4)

The augmented Lagrangian of the problem is

$${{{\rm{ L}}}} ({{{\bf{x}}}},\, \beta,\, \tilde{{{{\bf{x}}}}},\, {{{\bf{u}}}})={\Vert {{{\bf{s}}}}-\beta {{{\bf{H}}}}{{{\bf{x}}}}\Vert }_{2}^{2}+K{\beta }^{2}{\sigma }^{2}+\Re \{{{{{\bf{u}}}}}^{H}(\tilde{{{{\bf{x}}}}}-{{{\bf{x}}}})\}+\frac{1}{2}\rho {\Vert \tilde{{{{\bf{x}}}}}-{{{\bf{x}}}}\Vert }_{2}^{2},$$

(5)

where \(\Re \{\cdot \}\), \({{{\bf{u}}}}\in {{{{\rm{C}}}}}^{N}\), and *ρ* are the real-part operator, the dual vector, and the penalty parameter, respectively. By using the alternating optimization strategy, we decompose the problem as

$${\tilde{{{{\bf{x}}}}}}^{t+1}={{\mbox{arg}}\,\min }_{\tilde{{{{\bf{x}}}}}} {{{\rm{L}}}} ({{{{\bf{x}}}}}^{t},\, {\beta }^{t},\, \tilde{{{{\bf{x}}}}},\, {{{{\bf{u}}}}}^{t}),$$

(6)

$${{{{\bf{x}}}}}^{t+1}={{\mbox{arg}}\,\min }_{{{{\bf{x}}}}} {{{\rm{L}}}} ({{{\bf{x}}}},\, {\beta }^{t},\, {\tilde{{{{\bf{x}}}}}}^{t+1},\, {{{{\bf{u}}}}}^{t}),$$

(7)

$${\beta }^{t+1}={{\mbox{arg}}\,\min }_{\beta } {{{\rm{L}}}} ({{{{\bf{x}}}}}^{t+1},\, \beta,\, {\tilde{{{{\bf{x}}}}}}^{t+1},\, {{{{\bf{u}}}}}^{t}),$$

(8)

$${{{{\bf{u}}}}}^{t+1}={{{{\bf{u}}}}}^{t}+\rho ({\tilde{{{{\bf{x}}}}}}^{t+1}-{{{{\bf{x}}}}}^{t+1}).$$

(9)

The problems (6), (7), and (8) have closed-form solutions, namely,

$${\tilde{{{{\bf{x}}}}}}^{t+1}=\frac{1}{\sqrt{N}}{\left({e}\, ^{j{\tilde{w}}_{1}^{t+1}},\, {e}\, ^{j{\tilde{w}}_{2}^{t+1}},{{{\mathrm{..}}}}.,{e}\, ^{j{\tilde{w}}_{N}^{t+1}} \right)}^{T},$$

(10)

$${{{{\bf{x}}}}}^{t+1}=\frac{1}{2}{({{{{\bf{R}}}}}_{H}^{t})}^{-1}{\phi }^{t+1},$$

(11)

$$\beta=\left|\frac{\Re \{{{{{\bf{s}}}}}^{H}{{{\bf{H}}}}{{{{\bf{x}}}}}^{t+1}\}}{{\Vert {{{\bf{H}}}}{{{{\bf{x}}}}}^{t+1}\Vert }_{2}^{2}+K{\sigma }^{2}}\right|,$$

(12)

where \({{{{\bf{R}}}}}_{H}^{t}={({\beta }^{t})}^{2}{{{{\bf{H}}}}}^{H}{{{\bf{H}}}}+\frac{\rho }{2}{{{{\bf{I}}}}}_{N}\), \({\phi }^{t+1}=2{{\beta }}^{t}{{{{\bf{H}}}}}^{H}{{{\bf{s}}}}+{{{{\bf{u}}}}}^{t}+{\rho }{\tilde{{{{\bf{x}}}}}}^{t+1}\), and \(|\cdot |\) are the absolute value operators. Moreover, \({\tilde{w}}_{i}^{t+1}\) is the optimized discrete phase of the *i*th element in the *t* + 1 iteration. The detail is given in Supplementary Note 5. The optimization process is summarized in Table 1, and the detailed analysis of the computational complexity is also given in Supplementary Note 5.

It is worth noting that the discrete optimization algorithm remains available to larger-scale arrays owing to its excellent properties, including low complexity and fast convergence (see Supplementary Note 6 for more details). We also provide the compare-and-contrast study with other algorithms reported in the field of DIM, including the definition of the performance metrics and the discussion on the scalability to demonstrate the advances of our method, as shown in Supplementary Note 6.

### Experimental results

To demonstrate the PM-based DIM scheme, we conducted measurements on the fields radiated by a fabricated PM and compared the measured results with reference constellation diagrams. The experiments were conducted in a far-field anechoic chamber shown in Fig. 3a, b. The PM was excited by a sinusoidal carrier signal at the SMA port, and the coding sequences were imposed on PM via MCU (ST STM32F103C8T6). A horn antenna (HZ-90HA20NZJL) operating in 9.8 ~ 15.0 GHz was connected to the vector network analyzer (VNA, model Agilent N9010A) for signal detections. The fabricated PM is illustrated in Fig. 3c, which is composed of 8 × 8 elements with an overall size of 144 × 144 × 2.5 mm^{3}. Moreover, it provides the far-field control capability, while keeping low cost and power consumption. All received signals were rescaled using the power normalization and the factor *β*, and a uniform phase was added to compensate for the propagation effect (see Supplementary Note 7 for more details).

To clearly demonstrate the working mechanism and measurements of our scheme, three representative experiments are presented in the main text, including the single-channel 8PSK, single-channel 64QAM, and dual-channel 16QAM modes. However, the capability of our scheme is not limited to dual channels, and the measurements for the four-channel DIM are demonstrated in Supplementary Note 8 (even more channels can be achieved if the scale of the metasurface increases). In the experiments, the PM was firstly configured to work in a single-channel mode, giving rise to 8PSK symbols to be transmitted in the desired direction \((\theta,\, \varphi )=(-19^{\circ},\, 0^{\circ} \, \, )\). Figure 3d, e shows the phases and magnitudes of the measured signals, which are in excellent agreement with the reference symbols. Then, the PM was configured to give rise to 64QAM symbols to be transmitted in the desired direction \((\theta,\, \varphi )=(21.5^{\circ},\,0^{\circ} \, \, )\). Such a configuration is more difficult than the previous one since the modulation has a higher order, and the distances between adjacent symbols are closer. Figure 3f, g demonstrates the phases and magnitudes of the measured signals. The measured results exhibit good agreement with the reference symbols.

As shown in Fig. 4, we utilize the error vector magnitude (EVM) to quantify the performance of the directional security. The metric of EVM is more intuitive and commonly used in wireless communications to assess the quality of constellation diagrams. The EVM is defined at a given time as the vector difference between the measured and reference symbols: \(EVM=\sqrt{{\sum }_{i=1}^{M}{|{{{{\bf{s}}}}}_{ref}^{i}-{{{{\bf{s}}}}}_{mea}^{i}|}^{2}/{\sum }_{i=1}^{M}{|{{{{\bf{s}}}}}_{ref}^{i}|}^{2}}\), where **s**_{ref}, **s**_{mea}, and *M* are the reference, measured signals, and the number of symbols, respectively. The lower the EVM value is, the closer the measured symbols are to the reference. Figure 4a, c present the EVM distribution and signal structures of the 8PSK case, while Fig. 4b, d show the measured results of the 64QAM case. Clearly, the up-and-right subfigure of Fig. 4c and the left-and-down subfigure of Fig. 4d indicate that the measured signals in the desired directions are correct, while the other subfigures demonstrate that the signal structures in the unintended directions are wrong, which validates the directional security of the proposed method. We further propose the concept of a secure zone, defined by the value of EVM, to separate the desired users and eavesdroppers. The in-depth analyses are conducted on the minimum angular width and the maximum number of supported channels (see Supplementary Note 13 for more details).

To further validate the 2D and dual-channel DIM capabilities of the proposed scheme, we then conducted experiments where the PM independently and simultaneously radiated the 16QAM symbols towards two users located respectively in the directions of \(({\theta }_{1},\, {\varphi }_{1})=(28.5^{\circ},\, 0^{\circ} \, )\) and \(({\theta }_{2},\, {\varphi }_{2})=(28.5^{\circ},\, 90^{\circ} \, \, )\). Figure 5a, b show the phases and magnitudes of the measured signals, which match well with the reference symbols. Figure 5c, d present the radiation powers and phases of the measured fields in the *φ* = 0° plane associated with the 16QAM, respectively. The radiation powers form three clusters near the desired directions \(\theta=28.5^{\circ}\), and the average powers of each cluster are presented on the right. From the above results, we clarify that the discrete optimization algorithm can effectively engineer the main-lobe beams to align with the desired directions. The main lobe powers have exhibited strong correlations with the magnitudes of the transmitted constellation symbols, addressing the effectiveness of energy allocation. The radiation phases of the measured fields are regularly arranged in the vicinity of the target direction. Figure 5e, f presents the results in the *φ* = 90° plane, where we can observe similar features in the target direction. The above results demonstrate that the simultaneous dual-channel 16QAM is realized. Finally, the optimized digital coding sequences of the 8 × 8 PM that implements the dual-channel 16QAM are shown in Fig. 5g.

For a multi-user communication system, the evaluation of the crosstalk of different users is crucial because it determines whether the channels are independent. To obtain the complete cross-talk matrix (see “**Methods**” for the definition and calculation), we have conducted two experiments. Specifically, In the first group of experiments, the sixteen different symbols are sent to user 1, and the same symbols coded as “0000” are sent to user 2. In the second group of experiments, the sixteen different symbols are sent to user 2, and the same symbols coded as “0011” are sent to user 2. The magnitudes and phases of the measured fields are demonstrated in Supplementary Note 9, which is well consistent with the transmitted symbols. The measured constellation diagrams of each user are presented in Fig. 6a, b. Taking Fig. 6a as an illustrative example, the symbols received by user1 are in good agreement with the 16QAM, while the signals of user2 are clustered near the corresponding reference point. We then calculate the cross-talk matrix **C** according to Eq. (8), and the results are shown in Fig. 6c. The cross-talk values between different users are relatively low, and the maximum is 0.12 (−18.4 dB), implying acceptable performance under the dual-channel transmission. We remark that the transmission efficiency of different channels can be further improved by three methods: (1) increase the phase quantization level of the metasurface; (2) increase the scale of the metasurface to reduce the multi-user interference; and (3) employ an improved channel model to take the mutual coupling of the metasurface into account^{38}.

Figure 7a, c demonstrate the EVM distribution and signal structures in the *φ* = 0° plane, while Fig. 7b, d show the measured results in the *φ* = 90° plane. The EVM values in the vicinity of the target direction \(\theta=28.5^{\circ}\) are relatively lower than in other directions. The left-and-down subfigures in Fig. 7c, d indicate that the measured signals in the desired directions are correct, while the other subfigures demonstrate that the signals in the undesired directions are distorted, which validates the directional security of the proposed method.

DIM has the advanced characteristics of providing a secure zone wherein the received signals exhibit a similar structure to the constellation diagrams, which improves the robustness of directional communication. To assess the secure zone, we measured the symbols in the vicinity of the desired directions. Figure 8a–d present the signal structures of the 8PSK case, the 64QAM case, the 16QAM case in the *φ* = 0° plane, and the 16QAM case in the *φ* = 90° plane, respectively. The received signals in different directions gather near the reference constellation symbols, which verifies the secure-zone feature of the proposed DIM scheme. Furthermore, we have also presented the distorted constellation diagrams in other directions as a reference to prove the directional security (see Supplementary Note 10).

To further verify the availability of the DIM at different frequencies, we also measured the symbols in the vicinity of the operating frequency (namely, 11 GHz). Figure 8e–h presents the signal structures of the 8PSK case, the 64QAM case, the 16QAM case in the *φ* = 0° plane, and the 16QAM case in the *φ* = 90° plane, respectively. The received signals are clustered near the reference constellation symbols in the broadband frequency range, indicating that the proposed DIM scheme supports next-generation high-throughput wireless communication.