Supplementary Figure S2 Results from the design procedure (MT_{2} and MT_{4}).a, b Numericallysynthesized phase and magnitude distributions, respectively, of the transmission coefficient pertaining to the single nanoholes in the supercell (shown on top) of the MT_{2} design (with parameters as given in Table 1), for the copolarized (blue square markers) and crosspolarized (red circle markers) components, assuming an infinite periodic array of period , under normallyincident polarized planewave illumination at . Element #1 is chosen as phase reference. Continuous curves are guides to the eye only. c, d Same as above, but for MT_{4} design.
We verified that the lookup maps in Supplementary Figure S1, calculated for a specific interelement spacing (), could also be utilized for designs featuring slightly smaller values of (e.g., MT_{2} and MT_{3}), yielding only slight distortions. Supplementary Figure S2 shows the phase and magnitude profiles pertaining to the MT_{2}and MT_{4} designs (not shown in the main text for brevity). We highlight that, for the MT_{5} design, any value of the nanoholesidelengths would produce the same phase gradient. However, the chosen design was found to provide a good tradeoff between small size (so as to reduce ) and high crosspolar transmission coefficient.
The farfield intensity profiles in Figure 4 (as well as Supplementary Figure S7below) are computed in two steps. First, we compute the transmitted field pertaining to a MT of finite width ~ along the direction, and assumed as infinitelyperiodic in the direction, illuminated by a 1D Gaussian beam with waist size of impinging from the silica region (so as to partially mimic the modal field of the SMF28 optical fiber).The structure is terminated by Blochtype periodic boundary conditions along the direction, and by ad hocperfectly matched layers (PML) in the remaining directions. The above assumptions, which neglect the nonuniform polarization distribution (along the direction) of a realistic 2D beam fiberoptic illumination, are instrumental to maintain the computational burden within affordable limits, while still providing a good estimate of the illuminationtapering effects.
The nearfield distribution is computed by means of COMSOL Multiphysics, with the simulated region extending up to a distance of in air and in the fiber region, in order to guarantee the computational affordability. Once again, an adaptive meshing is employed, with maximum element size of in the uniform dielectric regions, 50nm in the air regions of the array, and a minimum number two element per skindepth in the gold layer,resulting in about 4.5 million degrees of freedom. The PARDISO solver is utilized, with default parameters.As a second step, the computed nearfield is propagated to the farfield region using wellknown formulae for the radiation from planar apertures (e.g., Sec. 4.1 in Ref. 2). In order to account for the lack of polarization control in the measurements shown in Figure 4 (as well as Supplementary Figure S7 below), the farfield intensity profiles are averaged over the two limiting cases of and polarized incidence. Results are normalized with respect to the maximum values.
From the same numerical simulations, we also estimate the efficiency of the designed MT prototypes, i.e., the fraction of incident power that gets transferred to the anomalous beam. To this aim, assuming anpolarized incidence, we compute the flux of the Pointing vector (realpart) associated with the crosspolarized transmitted field through a planar surface at a distance from the metasurface. We then obtain the efficiency by normalizing this quantity by the same flux pertaining to the incident field only, calculated in the absence of the metasurface, and by truncating the fiber region with a PML.
The field maps shown in Figures79are also computed via COMSOL Multiphysics by assuming a periodic 2D array under normallyincidentpolarized planewave illumination. Once again, Blochtype periodic boundary conditions are assumed along the  and directions, with ad hoc PML terminations along the direction. In this case, the simulated region extends up to a distance of in air and in the fiber region, with maximum meshelement size of in the uniform dielectric regions and of 20nm in the air regions of the array, and a minimum number two element per skindepth in the gold layer (resulting in about 2 million degrees of freedom). Also in this case, the MUMPS solver is utilized. Results are normalized with respect to the incidentfield amplitude.
For better computational affordability, the reflectivity spectra for these configurations are computed by means of a publicdomain numerical code (sourceforge.net/projects/rcwa2d/files/) that implements a 2D rigorous coupled wave analysis (RCWA).^{3} Once again, a normallyincidentpolarized planewave illumination is assumed. Convergence is achieved by using modes. To realistically simulate the experimental reflection setup, the reflectivity pertaining to the zeroth diffraction order is considered, corresponding to light traveling back along the fiber axis. Moreover, in order to roughly mimic a realistic deposition process, the SiO_{x} overlay is assumed as conformal to the MT surface (i.e., filling the nanoholes and covering the unpatterned gold regions).
