Epsilon-near-zero (ENZ)-based optomechanics

Optomechanics deals with the control and applications of mechanical effects of light that stems from the redistribution of photon momenta in light scattering. As an example, light-induced levitation of an in ﬁ nitesimally small dipolar particle is expected in front of epsilon-near-zero (ENZ) metamaterials. However, a theoretical understanding of these effects on single-material and multi-material larger particles is still lacking. Here, we investigate, analytically and numerically, optical forces on polarizable particles with size ranging from 20 nm to a 1 μ m in proximity of ENZ metamaterials. We look at the general features of the repulsive-attractive optomechanics from the nano to the microscale exploiting different theoretical methods (dipole approximation, ﬁ nite elements calculations, transition (T)matrix). We discuss the role of realistic layered materials, as our ENZ substrate, on optical forces and analyze the in ﬂ uence of composition and shape by studying a range of complex particles (dielectric, core-shell, plasmonic ellipsoids). Physical insights into the results are discussed and future research directions are forecasted. Our results provide possibilities in exploiting engineered materials and surfaces for the manipulation and tailoring of light-induced forces in optomechanics.


Introduction
Recent developments in the field of metamaterials and metasurfaces have provided useful platforms for manipulating and tailoring light-matter interaction with numerous applications ranging from cloaking [1,2], enhanced spontaneous emission [3], sensing [4], signal processing and information handling [5,6], and nonreciprocity [7], just to name a few.Among various classes of metamaterials, the epsilon-near-zero (ENZ) and near-zeroindex (NZI) structures have attracted increasing attention due to their unique features in light-matter interaction [8,9,10,11].In such structures, relative permittivity and/or relative permeability attain values near zero, thus making the effective refractive index of the structure near zero.Consequently at the operating frequency, the wavelength in these media is "stretched", making the phase of the signal approximately uniform across this structure [12].As a result, the waves exhibit "static-like" spatial distributions, while temporally they are dynamic.This has led to numerous exciting wave phenomena, with several potential applications [8,9,10,11,13].One such feature is the possibility of levitation of electrically-polarized nanoparticles in the vicinity of ENZ substrates [14].In our earlier work, we theoretically showed that an infinitesimally small nanoparticle, when electrically polarized at a given frequency, could be levitated when placed near an ENZ substrate.This phenomenon, which was inspired as a classical analogue to the Meissner effect (levitated magnets in proximity of superconductors), can provide a new approach in optomechanics when manipulation of electrically polarizable particles is desired in the presence of optical fields.
Careful manipulation of particles with light, which has a long history dating back to the pioneering work of Ashkin in 1970s [15,16], has played important roles in various areas, from biology [17] to nanoscience and nanotechnology [18].At the nanoscale, various methodologies have been used for such optical manipulation, including trapping [19,20], pushing [21,22], and binding [23,24], with different materials such as dielectrics, semiconductors, plasmonic, and biological [25].The surrounding media can be vacuum, air, or liquid.Optical tweezing [26,25] is usually achieved using optical beam shaping to generate desired potential traps [27].Recently, new approaches to optical manipulation of objects without a beam-shaping were proposed.Soljacic and coworkers [28] proposed that the motion of a Janus particle with spatially asymmetric absorption can be controlled by changing the incident wavelength.Ilic and Atwater [29] proposed self-stabilizing optical manipulation of macroscopic objects by controlling the anisotropy of the scattered light from the structured object's surface.Both approaches, however, rely on structuring the object in lieu of the incident light.
In the present work, we merge the two fields of ENZ metamaterials and of optical trapping, providing a new platform, which we name ENZ-based optomechanics, for manipulating and controlling mechanical motion of particles in vicinity of ENZ structures.We explore, numerically and analytically, how various parameters, such as the size, shape and composition of the particle and its distance to the ENZ substrate affect the optomechanical forces on such particles.We consider both homogeneous and layered structures as our ENZ substrates.In recent years, researchers have been able to tailor the effective permeability and permittivity of composite media by engineering the electric and magnetic resonances of nanostructures.Together with related developments in nanophotonics, metamaterials provide unprecedented freedom to define and sculpt electromagnetic modes.Metamaterials allow to alter the topology of photonic isofrequency surfaces -which govern the momentum and energy of optical modes inside a medium -contrarily to conventional bounded spherical and ellipsoidal isofrequency surfaces in natural dielectrics [30].Among many other extreme optical features, unbounded iso-frequency surfaces in hyperbolic dispersion metamaterials [31] and point-like vanishing surfaces in epsilonnear-zero (ENZ) media [32] constitutes two examples of advanced modal engineering.In particular, epsilonnear-zero metamaterials provide extended modes with uniform phase over micrometer length scales inducing profound effects on nanoscopic light-matter interactions.These deeply subwavelength structured surfaces support unique electromagnetic modes that can be used in sub-diffraction imaging, [33] and waveguiding [8], spontaneous emission engineering [34] and biosensing [35].
In the following, we introduce the geometry of the problem, and discuss the electromagnetic modeling for the structure, along with the dipole approximation.We also present extensive numerical results, based on the finite-element method (using the commercial software COMSOL Multiphysics ® ) and on the T-matrix methods.We also present a series of results for various parameters involved in the problem.Physical insights into the results are presented and future directions are discussed.
Geometry of the problem.Figure 1a presents the geometry of our problem.A polarizable particle, made of a single nonmagnetic material (or multilayered materials), surrounded by an external medium (e.g., water) of We consider a generic particle, in principle of any shape and composition, in front of a metamaterial surface at an edge-to-edge distance, h , immersed in an external medium (e.g., water) of refractive index n m .The origin of the coordinate system is placed on the surface so that the z axis is positive in the semi-infinite space where the particle resides.A monochromatic optical wave illuminates the particle and surface at normal incidence so that the total field is the superposition of incident (E I ), reflected (E R ) and scattered fields (E S , E SR ).The resulting optical force can be either attractive (negative) or repulsive (positive).b) Near-field force (R , φ) map.We explore the range of reflectivity, R , and phase, φ, related to the reflection of the incident wave on a generic surface.The ideal ENZ surface (R = 1, φ = 0)) is in the top right corner of the map and shows a repulsive force.Here the near-field force component is calculated in the dipole approximation for a polystyrene (dielectric constant p =2.543 at 560 nm)) particle of radius a = 20 nm at a fixed edge-to-edge distance of h = 10 nm in water.c) Curves of (R , φ) for different substrates consisting of alternating metal and dielectric layers, in the thin layer limit where effective medium theory is valid.The color of each curve indicates the metal filling fraction.The zero force lines from panels d and e are superimposed as dashed lines.d-f) Total force (R , φ) maps as calculated from T-matrix methods for a dielectric particle size of a = 20 nm (d), a = 220 nm (e), and a = 1000 nm (f), respectively, and at a fixed h = 10 nm.The total force maps have a structure that is strongly dependent on particle size.This is due to the increase of the scattering force component that, for large particles, overcomes any other gradient-like force component that is dominant for nanoparticles.refractive index n m , is located at an edge-to-edge distance h above a metamaterial substrate.The particle can be spherical (or other shapes as will be discussed later in the manuscript), and it is made of a (dielectric or metallic) material with a relative permittivity p .The substrate can be considered as a homogenized nonmagnetic medium with relative permittivity near zero at the frequency of operation, or a layered structure engineered to function as ENZ.
A monochromatic optical wave is illuminating this structure at normal incidence.The goal is to evaluate the optical force on the particle and to investigate how various parameters, radius a , edge-to-edge distance h , particle's permittivity, and signal frequency affect the optical force's magnitude and direction, i.e. whether it is a repulsive (positive) or an attractive (negative) force.In the next section, in order to gain some physical insight we start by assuming the polarizable particle to be represented by an infinitesimally small electric dipole, and discuss the analytical approach for evaluating the force acting on this particle.In the subsequent sections, we will expand our approach to include the full-wave numerical simulations of the problem, allowing to consider realistic sizes and shapes for this particle.
Dipole approximation.We first consider a particle size much smaller than the light wavelength (a << λ) so that optical forces can be calculated analytically within the dipole approximation (DA) [36,37,25].Due to its simplicity, the dipole approximation can provide useful results that can be compared with more complex light scattering approaches (T-matrix, finite elements methods) at the nanoscale [38].
We start our analysis from the near-field force component.It has been shown [14] that in front of an ENZ surface an emitting point dipole source is subjected to a near-field repulsive force, reminiscent of the Meissner effect in superconductors [14].This portion of the force, which we refer to as the "near-field" force, is due to the interaction of the emitting dipole with the substrate (excluding the force due to the presence of the incident and reflected waves).When all forces are considered (including the forces caused by the incident and reflected waves), the forces are called "total force".We can extend the result of Ref. [14] to a finite-sized polarizable particle illuminated by an incident field by considering the radiated power upon scattering, P rad = σ scat I (z ), in terms of the scattering cross section, σ scat , and light intensity, I (z ).Thus, the near-field force component is (see Supp.Info.): where z is the axial coordinate (z = h +a , a is the radius of the particle, h is the edge-to-edge distance of the particle from the surface), c is the vacuum speed of light, m = 0 n 2 m is the permittivity of the surrounding medium, n m is the refractive index of the medium, and s is the complex dielectric permittivity of the ENZ surface.
In Fig. 6, a panel summarizing the results of the calculation of the near-field force (d-f) on a 20 nm dielectric bead in water is shown.Three different surfaces are considered: lossless (Im( s )=0), with medium loss (Im( s )=0.5), and with high loss (Im( s )=0.8).The comparison with the results obtained for a point dipole in vacuum [14] shows that in this work the presence of a medium (water) broadens the repulsive near field force region from −1 < s < 1 to −1.77 < s < 1.77; moreover, as already observed [14], even in surfaces with high loss there is still a repulsive near-field force.
In order to explore how the ENZ surface can influence the near-field and total forces on the particle, we evaluate such effects in terms of the amplitude ρ and phase φ of the complex reflection coefficient of an incident wave from this surface.In Fig. 1b, the near-field force on a a =20 nm radius dielectric bead at h =10 nm from the surface has been calculated as a function of the surface reflectivity R = ρe i φ 2 = ρ 2 and phase angle φ which are connected to the surface complex refractive index ñ = n s + i k s by [39]: Here we use e −i ωt as our time harmonic convention.The calculated near-field force can reach a fraction of femtonewton for an incident intensity of approximately 5.6•10 8 W/m 2 (corresponding to a typical experimental configuration, see Sect.1.2 of the Suppl.Info.) and changes character from attractive to repulsive when the reflection phase angle changes from φ=-π to φ = 0. Metals such as Au or Ag, having a certain amount of absorption  d,e) Size dependence of the total optical force for fixed edge-to-edge distance, h =10 nm.For small particles the gradient force component of the partially reflected plane wave dominates, resulting in a dependence with particle size, while for large particles radiation pressure has a major contribution resulting in a negative force pushing the particle towards the surface.
(k s in Eq. 2), are in the attractive region of the near-field force (compare Figs. 1b and c).On the contrary, in front of an ideal ENZ surface, having R =1 and φ=0, the near-field force is repulsive.Substrates of alternating metal and dielectric layers can span a broader range of φ and R values.In the limit of layers much smaller than the incident wavelength, where effective medium theory (EMT) is valid, we show the (R , φ) results for four different metal / dielectric mixtures in Fig. 1c (see Sect. 4 of the Supp.Info.for more details on the EMT calculation).
The metal filling fraction is indicated by the color of the curve.Depending on the fraction, we can switch the sign of the force from attractive to repulsive and vice versa.If we go beyond EMT and take into account the finite thickness of layers in real structures, as described in the discussion of Fig. 3 below, we can achieve an even wider range of φ and R values.
We now consider the total optical force from an incident field on a nanoparticle calculated in DA.This is the sum of two main components: a gradient force, F grad , and a scattering force, F scat [25].For plane wave illumination (for Gaussian beams see Sect.1.2 of Supp.Info.) impinging normally to the ENZ surface, the force components are influenced by incident and reflected fields.Thus, considering only the axial direction z , they are written as (see Supp.Info.): where α is the particle complex polarizability [40], α 0 is the Clausius-Mossotti polarizability and σ ext = k m Im(α) is the extinction cross-section, related to the particle absorption and scattering [37,25], with k = 2πn m /λ the wave number and λ the wavelength.
The gradient force, F grad , drives the particle towards the maximum (minimum) of the modulated light intensity profile for positive (negative) real part of the polarizability.On the other hand, the scattering force, F scat , is constant with respect z , and it always pushes the particle along the beam propagation direction.
In Fig. 1d, the (R , φ) map of calculated total axial force on a polystyrene a =20 nm bead (dielectric constant p =2.543 at 560 nm) at h =10 nm from the surface is shown.The total force is one order of magnitude larger than the near-field force (Fig. 1b) and shows a change in the repulsive-attractive character when the phase angle changes from −π to 0 , respectively.This is due to the gradient force (see also Fig. 2a, short dots) that dominates the optomechanical response and drives the particle towards the high field intensity regions.The change of the phase of the reflection coefficient shifts the intensity modulation resulting from the interference between the incident and reflected wave.Thus, the high intensity points shift accordingly and the sign of the force changes around φ ∼ −π/4.
We now calculate the total optical force on the dielectric bead in front of glass, Ag and ENZ surfaces as a function of distance, h .The ENZ material is chosen so that n s ≈ 0.476 and k s ≈ 0.511, in order to obtain a real part of complex permittivity close to zero and an imaginary part close to 0.5 to include unavoidable losses of realistic systems.This choice leads to values of R and φ similar to those in the experimentally fabricated layered substrates described below, corresponding to the point marked with a red star in Fig. 3.The strong modulation resulting from the standing wave is clearly visible.The points with zero force and negative slope are trapping points that correspond to equilibrium positions for the particle dynamics (arrows in Fig. 2a).For the case of the ENZ surface the equilibrium point closest to the surface occurs at h ∼10 nm, while for the glass and Ag surface they occur at h ∼ 85 nm and h ∼ 60 nm, respectively.By linearizing the force at the equilibrium points, F (z ) ≈ −κz , a trap spring constant κ can be calculated.The trap spring constants κ ENZ and κ Ag calculated in front of ENZ and Ag surfaces can be compared to the spring constant κ S calculated, in a standard single-beam optical tweezers setup, with the same particle and at the same light intensity (see Supp.Info.).The spring constants are κ ENZ = 15 fN/µm and κ Ag = 27 fN/µm, while the trap spring constant κ S in a standard optical tweezers setup is two order of magnitude lower, κ S = 0.24 fN/µm.The beneficial effect of the ENZ and Ag reflective surfaces on the trapping is evident.The increasing size of the particles corresponds to larger optical forces and different trapping points (see Fig. 5 in Supp.Info.for DA calculations on larger size nanoparticles at 50 and 100 nm).
Full-wave simulations.In order to calculate optical forces on larger particles, we use two different full-wave modeling approaches based on the transition (T-)matrix formalism [41,42] and on finite-elements methods using the commercial software COMSOL Multiphysics ® , respectively.In particular, electromagnetic scattering from particles near to or deposited on a plane surface that separates two homogeneous media of different optical properties in the T-matrix formalism [42,43,44,45] can give account on the role of the different multipoles in the particle-surface interaction (see Supp.Info.for details).Indeed, the presence of the surface can have a striking effect on the scattering pattern from the particles, because the field that illuminates the particle is partly or totally reflected by the surface and the reflected fields contribute both to the exciting and to the observed field.Moreover, the field scattered by the particle is reflected by the interface and thus contributes to the exciting field.In other words, there are multiple scattering processes between the particles and the interface.As a result, the field in the accessible half-space includes the incident field E I , the reflected field from the interface E R (as we would have if no particle were present), the scattered field from the particle E S and, finally, the field that after scattering by the particle is reflected by the surface, E SR , related to E S by the reflection condition (see Fig. 1a).Thus, the observed field, superposition of E S and E SR , includes all the scattered and scattered-reflected multipole contributions (see Supp.Info.for more details).
It is possible to define the T-matrix for particles in the presence of the interface that is the starting point to calculate optical forces and torques either by direct integration of the Maxwell stress tensor (MST) over a closed surface containing the particle [25] or by exploiting the general expressions of optical force and torque in terms of multiple expansion [46,47,48].The optical force is obtained in COMSOL by direct integration of the MST that is calculated based on the total electric field and magnetic field which include the incident fields, the scattered fields by the particle, and all the reflected fields by the surface (see Supplementary Information for details on full wave methods).
The results obtained in DA from the different surfaces are compared in Fig. 2a with those obtained by using full electromagnetic calculations based on COMSOL (circles), and T-matrix methods (continuous lines).A very good agreement is clearly observed.In all approaches, the total optical force on small particles is modulated by the sinusoidal term of the gradient force.Its magnitude is larger (in the fN range) on more reflective surfaces and its phase changes sign going from a Ag to an ENZ substrate, leading to the formation of optical trapping points at different distances (arrows in Fig. 2a).In brief, the gradient force dominates the ENZ-optomechanics for small particles even in proximity of the surface.
T-matrix and COMSOL allow the calculation of optical forces for larger particles than in DA.In Fig. 1e and f, the (R , φ) map of total axial force calculated with the T-matrix approach on a a =220 nm bead (Fig. 1e) and a =1 µm bead (Fig. 1f), at h =10 nm from the surface, are shown.The comparison with Fig. 1d highlights the strong dependence of the total optical force on the bead size.The repulsive-attractive behaviour is driven by the competition between gradient force and scattering force, which may give repulsive behaviour, for intermediate size beads, in surfaces having large reflectivity (see Fig. 1e); however, at large bead size (Fig. 1f), scattering force overcomes gradient force, and the total optical force is attractive in front of every type of surface.
In Fig. 2b and c the T-matrix calculations of the total optical force on larger particles are shown as a function of the edge-to-edge distance from ENZ, Ag and glass surfaces.The larger size of these particles with respect to the nanosized bead in Fig. 2a highlights the increased contribution of the scattering force on the gradient force.The scattering force is detrimental towards stable equilibrium positions in front of glass surface for the 220 nm radius bead and in front of both ENZ and glass surfaces for 1 µm radius bead.The lower reflectivity of these surfaces as compared to the reflection from the Ag surface does not allow an efficient balance between scattering force from incoming and reflected beams, increasing the scattering force contribution with respect gradient force and hindering the trapping.
In Fig. 2d and 2e the results are reported for increasing bead size at fixed distance, h =10 nm, from the ENZ, Ag or glass surfaces.It is shown that at small bead size (below approximately 300 nm radius), the gradient force modulates the total force.At increasing bead size, the particle extinction cross section increases, consequently the scattering force is predominant on the gradient force, inhibiting equilibrium points and inducing an effective attractive force directed towards the surfaces.
Epsilon-Near-Zero Metamaterials.Regarding layered ENZ materials, we have demonstrated experimentally that it is possible to control the optical topology and to induce the ENZ behavior by designing and fabricating subwavelength layered lattice structures as a result of interlocking noble metals and dielectric thin films [49].Upon selecting metal-dielectric bilayers, the thickness of each layer, the filling fraction and the number of bilay-  Curve labeled EMT is the effective medium approximation to the system, which corresponds to n → ∞, d → 0 with n d = 500 nm.In all the above cases the superstrate is water and the substrate is glass.For comparison we show points indicating the R and φ values for a simple interface between a water superstrate and a pure material substrate (Ag, Au, Ge, TiO 2 , Al 2 O 3 , and an ideal ENZ).We also show experimental results (green stars, details in the Supp.Info.) involving a water superstrate and 5 trilayers (Al 2 O 3 /Ag/Ge from top to bottom, where Ge is present as a thin wetting layer to ensure fabrication quality).The dotted green trend line corresponds to keeping the Ag and Ge layer thicknesses fixed at 15 nm and 2.5 nm respectively, while varying Al 2 O 3 thickness from 80 nm to 20 nm (left to right).In order to compare EMT calculation with the full-wave analysis, COMSOL is used (black diamonds) to calculate R for different layered structures with various metal's filling fraction (0.4 and 0.6), layer's thickness (50 nm and 100 nm) while keeping the total thickness of the layered structure unchanged (500 nm) and also different order of material in the stack (metal on top, blue labels, and dielectric on top, red labels).
ers, the frequency of the optical topological transition in the iso-frequency surface leading to the epsilon-nearzero behavior can be tailored.The lattice structure is fabricated as a five tri-layer system using Al 2 O 3 , Ag, and Ge from top to bottom.The Ag layer thicknesses were in the range of 10-25 nm, with a thin Ge layer (1-3 nm) underneath to ensure surface wetting.The Al 2 O 3 layer thicknesses were systematically varied between roughly 20 nm and 80 nm across different material systems (Fig. 3), subsequently tuning the frequency of the topological transition.In previous studies we used effective medium theory to calculate the dielectric permittivity of the entire structure, as opposed to more recent inverse design approaches to account for a wider material parameters space.We perform spectroscopic ellipsometry measurements to evaluate the dielectric tensor components and the dispersive behavior of the layered structure.By fitting the measured angular reflectance and the ellipsometry parameters ψ and ∆, we can directly obtain the effective optical constants of the multilayer slab.Using the transfer matrix method, we can then predict the magnitude and phase of reflection at normal incidence with a water superstrate.The green stars in Fig. 3 represent these predicted values from 6 samples consisting of a 5 bilayer Al 2 O 3 /Ag thin-film stack with a Ge seed layer to ensure the uniformity of the Ag films.By varying the thickness of the Al 2 O 3 layers, we covered a phase range of ∆Φ ≈ 180 • and reflectance range of ∆R ≈ 0.5.The full range of accessible R and φ values is even larger if we expand the design space of the substrate to include different numbers of bilayers and metal filling fractions.The thick curves in Fig. 3 show transfer matrix calculations of (R , φ) for Al 2 O 3 /Ag stacks with different structural parameters indicated by the labels.In all cases the total thickness of the stack was kept fixed at 500 nm.The color at each point along the curves corresponds to the metal filling fraction.In the limit of many thin bilayers we approach the EMT result of Fig. 1c, which is also reproduced here for comparison.Note that actual layered materials can achieve positive values of φ, while homogeneous materials (for example those described by EMT) are confined to the φ < 0 subspace.
Complex particles (core-shell, ellipsoids, ENZ) In addition to spherical beads, we evaluate the optical forces on different types of particles in front of dielectric, metallic or ENZ surfaces.We consider spherical core-shell particles based on SiO 2 and Ag, an Ag prolate spheroid and a spherical particle made by an ENZ material.
We first used COMSOL simulation to calculate the forces on core-shell structures in front of layered ENZ material at 560 nm illumination.The total particle radius a tot is fixed at 20 nm.The particles had alternatively SiO 2 or Ag as the core, with the other material as the shell.In Fig. 4a the total force on the core-shell particles as a function of the distance from the ENZ surface is shown.It is clearly observed that the presence of Ag in the outer shell enhances the total force with respect of the inverse structure having SiO 2 as the shell, but also with respect to the pure Ag sphere.The highest value of the force is found (red curve in Fig. 4a) for a SiO 2 -Ag core-shell structure having a core radius of a 1 = 16.1 nm and an Ag shell 3.9 nm thick which, as shown in Fig. 4c, is at the resonance condition at the ENZ wavelength.
As shown in Figures 4d, 8d and 8e, the particle resonance at 560 nm enhances the optical force to the piconewton range (Fig. 4d) but only at very short distances from the surfaces, being repulsive in the ENZ case (Fig. S4d) and attractive in the Ag case (Fig. 8e).Otherwise, the total optical force is at the fN range.
Specifically, at the resonance F enz is in the piconewton range close to the ENZ surface (from h =0 nm to roughly 10 nm).The gradient force, F grad , has an oscillating character, but its amplitude is smaller (≈ 1 fN) than F enz , due to the small real part of the polarizability at resonance (Re(α) = 0.04 • 10 −32 Fm 2 ).On the contrary, F scatt is large (tens of femtonewton), because of the large extinction cross section at resonance.Thus, at 560 nm (black curve in Fig. 9a), the total force is repulsive and in the piconewton range close to the surface, but becomes attractive and approximately constant as the F enz contribution fades off with distance.
The behaviour of the forces on the core-shell particle can also be studied for wavelengths smaller and larger than the particle plasmon resonance (Suppl.Info.).The calculation has been made for 552 nm, on the blue side of the plasmon resonance, and at 566 nm, on its red side.At these wavelengths, the scattering force is slightly lower than at resonance, while F grad increases by at least one order of magnitude.For this reason, its oscillating character shows up in the total force (Fig. 9a, blue and red curves).Moreover, as the polarizability changes sign from one side to the other of the resonance, also the gradient force inverts its phase from the blue to the red side of the resonance.Similar discussions hold for the optical forces in front of Ag surface (Figure 9b); however, in this case, the F enz is attractive close to the surface.
We now consider an Ag prolate spheroid as a prototypical non-spherical particle.This is chosen with a long axis a 1 = 56.8nm and short axes a 2 = a 3 = 20 nm.As shown in Figure 4e, the particle has, in water, a long axis resonance at 560 nm and a short axis resonance at 360 nm.For the calculation of the total optical forces we considered the case in which the spheroid has the long axis aligned with the wave polarization, and the short semiaxis as the size parameter in Eq. 1.We obtain a further enhancement of the total optical force (tens of piconewton, Fig. 4f) which, as in the core-shell structure, is repulsive in front of ENZ surface and attractive in front of Ag surface.In Figure 4f a contour plot of the total optical force, calculated as a function of the surface reflectivity R and phase shift φ, namely, in front of all possible surfaces, is shown.We clearly see that the repulsive force can be close to 200 pN in front of an "ideal" ENZ surface, having the maximum reflectivity and a vanishing phase shift.
In the case of ENZ particles, we used the same n and k values used for the ENZ surface.We calculated optical forces in front of glass, Ag or ENZ surfaces.The calculation has been made for ENZ beads having radii a =20, 50 and 100 nm.As shown in Fig. 7, the forces are about five times larger than the ones observed in dielectric bead counterparts.The larger scattering force of ENZ particle hinders its trapping in front of glass surface, for all radii.Moreover, the 100 nm radius ENZ particle cannot be trapped also in front of ENZ surface (Fig. 7c).Results are shown in Supplementary Information.
Finally, we have studied the total optical force in case of a focused (NA=1.3)Gaussian beam, typical of optical tweezers experiments (Section S1.2).The calculations, made for a 20 nm radius polystyrene bead, show that, both in front of ENZ (Fig. 11a) and Ag (Fig. 11b) surfaces, the beam focusing induces a fading of the total force with the distance h (see Fig. 11).The extension of the calculations for beads with larger radius (contour plots of the total optical force in front of ENZ, Fig. 11c, and Ag, Fig. 11d, surfaces) shows that the total force increases at increasing bead radius, reaching the range of tens of femtonewton in front of ENZ and hundreds of femtonewton in front of Ag surface.The modulation induced by the gradient force is clearly visible.It is worth noting that when Gaussian beams are used, for a direct comparison, the beam power is reduced with respect to the plane wave case in order to maintain the intensity at the beam focus similar to the plane wave intensity.

Conclusions.
In conclusion, ENZ-based optomechanics represents a novel way to manipulate and tailor mechanical effects of light exploiting flat surfaces.We focused our study on the repulsive-attractive optomechanics for particles in front of an ENZ surface in realistic conditions for a wide range of parameters (particle size and shape, ENZ surface structure, etc.) in the axial direction.Combining the unique optical properties of ENZ metamaterials with patterning capabilities will also enable further manipulation and control in the transverse direction towards a full dynamical engineering of ENZ-based optical forces.Various potential applications for future study include particle sorting due to the strong dependence of ENZ-based optical forces on the size and material composition of particles, biomolecular trapping and sensing, wavelength multiplexing of optical forces, and chiral optical sorting, just to name a few.

S1 Optical forces in the dipole approximation in front of epsilon-near-zero Materials
Dipole approximation (DA) is an easy and quick method to calculate optical forces on a nanoparticle.It is valid when the size of the particle is very small compared to the wavelength of the field [36,37,25], and due to its simplicity it provides useful results that can be compared with more complex light scattering approaches (Tmatrix, DDA) in the limit of small particles [38].
In DA the total optical force on a particle is usually split in a gradient force F grad and a scattering force F scat [25]: Here, k is the wave propagation direction that for an axially directed plane wave coincides with the axial coordinate ẑ , r is the radial coordinate, c is the speed of light, m = 0 n 2 m is the medium permittivity, 0 is the vacuum permittivity, n m is the refractive index of the medium, I (r, z ) is the wave intensity, α is the particle complex polarizability, where α 0 is the polarizability in the static field limit (Clausius-Mossotti), and σ ext is the extinction cross-section, related to the particle absorption and scattering [37,25]: with k = 2πn m λ the wave number and λ the wavelength.F grad drives the particles towards the maximum of light intensity if they have positive polarizability; otherwise, the particles are repelled from it.On the contrary, F scat always pushes the particles along field propagation direction, k .Another contribution to the total force may come from the spin-curl force [25], but only when beams having spatial polarization gradients are used [18], which is not the case in this work.
Recently, it has been proposed [14] that in front of an -near-zero (ENZ) surface a point dipole source is subjected to a near-field repulsive force, reminding the Meissner effect in superconductors [14].In the quasistatic approximation the near-field force is [14]: where σ is a prefactor accounting for the orientation of the dipole (σ=1, horizontal dipole; σ=2, vertical dipole), s is the complex dielectric permittivity of the surface, z is the height of the dipole above the surface and P rad is the radiated power of the dipole in free space.

S1.1 Plane wave illumination
Here, we calculate the total optical force on a finite-size particle in front of an arbitrary reflective surface in the dipole approximation.In this case, the exciting field E E is the superposition of the incident E I and reflected E R electromagnetic waves which produces a standing wave that, in the simplest case of plane waves travelling in the z direction, can be written as: with I 0 = n m 0 c E 2 0 /2.Note that z is taken positive in the direction of the reflected beam and ρ and φ are the amplitude and phase, respectively, of the complex reflection coefficient of the surface r m = ρe i φ , which is connected to the surface complex refraction index ñ = n s + i k s by [39] Thus, the gradient force along the axial direction on a finite-size particle in front of a reflective surface can be written as: where z = h + a is the axial coordinate, h is the edge-to-edge distance of the particle from the surface and a is the particle radius.The scattering force is the sum of the opposite contributions due to the incident and reflected plane waves [50,39,51]: where ρ is related to the surface reflection coefficient Finally, the near-field force on the particle is: where the radiated power P rad =σ ext I (z ) is related to the light scattering process.
We can now add Eqs.S13-S15 to calculate the total optical axial force on different types of particles in front of dielectric, metallic or ENZ surfaces.As the dipole is induced by a linearly polarized wave travelling ortogonally to the surface, an horizontal dipole (σ=1 in Eq.S15) is used.We used an incident light intensity of about 5.6•10 8 W/m 2 , corresponding to a beam power of 10 mW and a beam waist of approx.3.5 µm, both of which can be realized in a typical experimental configuration in our laboratories.Stable equilibrium points for the particle dynamics are found at z values in which the total optical force vanishes with a negative slope.For small displacements from these points, the particles are subjected to a restoring force that can be linearized as F z ≈ −κ z z , with κ z the trap spring constant.We consider four types of model particles: a homogeneous dielectric (polystyrene) spherical bead, a spherical particle with parameters equivalent to an ENZ material, a spherical core-shell particle (SiO 2 core, Ag shell), and an Ag prolate spheroid.The different surfaces have been considered in the calculation by means of their complex refractive index values at 560 nm; in the ENZ case, we have chosen n s ≈ 0.476 and k s ≈ 0.511 in order to obtain a real part of complex permittivity close to zero and an imaginary part close to 0.5.The same values have been used for the complex permittivity of the ENZ particle.Dielectric bead.We calculated the optical forces under λ=560 nm in water (n m =1.33) on spherical polystyrene (relative permittivity 2.54) beads having radii a =20, 50 and 100 nm.In this case, the Clausius-Mossotti polarizability is [52]: In Fig. 2a of the main text, the results (short dots) obtained in DA approximation for the 20 nm dielectric bead as a function of its distance from the different surfaces are compared with those obtained by using more sophisticated approaches (COMSOL, cicles, and T-matrix, continuous lines).A very good agreement is clearly observed.In all approaches, the total optical force on small particles is modulated by the sinusoidal term in the gradient force.It is larger (in the fN range) on more reflective surfaces and, going from Ag to ENZ, it changes phase, leading to stable traps at different distances (arrows in Fig. 2).and c) 100 nm radius dielectric beads.The beads are in front of ENZ (blue curve), Ag (magenta curve) or glass (red curve) surfaces.Note how for smaller particles the dominant contribution to the optical force comes from the gradient force, while for larger particles the greater scattering force shifts downwards the force modulation resulting from the interference between incident and reflected field.In d) and e), 3D plots of the total force as a function of the particle radius and of the distance from ENZ (d) and Ag (e) surfaces.
The axial trap spring constants κ ENZ and κ Ag in front of ENZ and Ag surfaces have been calculated by a linear fit of the total force at the equilibrium points.They are κ ENZ = 15 fN/µm and κ Ag = 27 fN/µm, that can be compared to the trap spring constant in the axial direction obtained in a standard optical tweezers setup, based on a single Gaussian beam.In this case, where w 0 is the beam waist, w (z ) = w 0 1 is the Rayleigh range, z 0 is the position of the beam waist and I 0 = 2P /πw 2 0 is the on-axis intensity at the waist of a beam having total power P .To evaluate the beam waist, we used the Abbe criterion, w 0 = 0.5λ N A , with NA=1.3 the numerical aperture, as in typical single beam optical tweezers.Eqs.S6 and S7 can be used to calculate the axial component of the total force, and the corresponding k S at the equilibrium point is obtained by a linear fit.We consider as before the particle in water and illuminated at λ=560 nm.We find, at the same light intensity used in front of ENZ and Ag surfaces, a two order of magnitude lower κ S = 0.24 fN/µm axial spring constant.
In Fig. 5 the optical forces in DA on dielectric beads at increasing bead radius (20, 50 and 100 nm) are shown.The increasing size of the particles corresponds to larger optical forces and different trapping points.However, Figure 6: Contour plots of the total optical force (a-c) and of the near-field force (d-f) in front of (a,d) lossless, (b,e) medium loss, Im( s )=0.5 and (c,f) very high loss, Im( s )=0.8 surfaces on a 20 nm dielectric bead as a function of of the real part of the surface permittivity and of the particle normalized height h /λ above the surface (λ=560 nm).The maximum optical force is in the fN range.
the 100 nm bead is not trapped in front of glass surface, whereas it is trapped in front of ENZ and Ag surfaces, whose higher reflectivity with respect glass surface better counteracts the scattering force due to the incoming beam.
In Fig. 6, a panel summarizing the results of the calculation of the total optical force (a-c) and of the near-field force (d-f ) on a 20 nm dielectric bead in water is shown.Three different surfaces are considered: lossless, with medium loss (Im( s )=0.5) and with high loss (Im( s )=0.8).The comparison with the results obtained for a point dipole in vacuum [14] shows that in this work the presence of a medium (water) broadens the repulsive near field force region from −1 < s < 1 to −1.77 < s < 1.77; moreover, as already observed [14], even in surfaces with high loss there is still a repulsive near field force.However, the calculation of the total optical force gives values not higher than 1 fN, which is found only in front of lossless surfaces.
ENZ particle.Optical forces on spherical beads made by ENZ material have been calculated in front of glass, Ag or ENZ surfaces.ENZ beads having radii a =20, 50 and 100 nm have been considered.As shown in Fig. 7, the forces are always larger than the ones observed in dielectric bead counterparts.The larger scattering force of ENZ particle hinders its trapping in front of glass surface, for all radii.Moreover, the 100 nm radius ENZ particle cannot be trapped also in front of ENZ surface (Figure 7c).

Core-shell particle.
To enhance the optical force, we used a SiO 2 -Ag core-shell particle designed to be resonant at approximately 560 nm and having a total radius a t o t = 20 nm.The calculation of the extinction cross-section shows (Fig. 8a) that the resonance condition is fulfilled if the core-shell structure has a core radius a 1 = 16.1 nm, and the Ag shell thickness is 3.9 nm.The particle polarizability is [52]: In this equation, a t o t is the core-shell total radius, 1 and 2 are the core and shell complex permittivity, respectively, and f = a 1 a t o t is the ratio between the core radius a 1 and the total particle radius a t o t .As shown in Figure 8, the resonance at 560 nm enhances the optical force to the pN range but only at very short distances from the surfaces, being repulsive in the ENZ case and attractive in the Ag case.Otherwise, the total optical force is at the fN range.
More specifically, at the resonance F enz is in the pN range close to the ENZ surface (from h =0 nm to roughly 10 nm).F grad has an oscillating character, but its amplitude is smaller (≈ 1 fN) than F enz , due to the small real part of the polarizability at resonance Re(α) = 0.04 • 10 −32 Fm 2 .On the contrary, F scatt is large (tens of fN), because of the great extinction coefficient at resonance.Thus, at 560 nm (black curve in Fig. 9 a), the total force is repulsive and in the pN range close to the surface, but becomes attractive and approximately constant as the F enz contribution fades off.
The behaviour of the forces on the core-shell particle can also be studied for wavelengths smaller and larger than the particle plasmon resonance.The calculation has been made for 552 nm, on the blue side of the plasmon resonance, and at 566 nm, on its red side.At these wavelengths, the scattering force is slightly lower than at resonance, while F grad increase by at least one order of magnitude.For this reason, its oscillating character now can be better noticed in the total force (Fig. 9 a, blue and red curves).Moreover, as the polarizability changes sign from one side to the other of the resonance, also the gradient force is "out of phase" going from the blue to the red side of the resonance.Similar discussions hold also for the calculation of forces in front of Ag surface (Figure 9 b); however, in this case, the F enz is attractive close to the surface.

Ag prolate spheroid.
We choose an Ag prolate spheroid having long axis a 1 = 56.8nm and short axes a 2 = a 3 = 20 nm.The particle polarizability is In this equation, p is the particle permittivity and L i is a geometric factor relative to the spheroid axis a i .In case of a prolate spheroid, L 1 is and ).As shown in Figure 10a, the particle has, in water, a long axis resonance at 560 nm and a short axis resonance at 360 nm.For the calculation of the total optical force we considered the case in which the spheroid has the long axis aligned with the wave polarization, so to use α 1 for the polarizability in Eqs.S13 and S9, and the short semiaxis as the size parameter in Eq.S15.We obtain a further enhancement of the total optical force (tens of pN) which, as in the core-shell structure, is repulsive in front of ENZ surface and attractive in front of Ag surface.In Figure 10b a contour plot of the total optical force, calculated as a function of the surface reflectivity R and phase shift φ, namely, in front of all possible surfaces, is shown.We easily see that the repulsive force can be close to 200 pN in front of an "ideal" ENZ surface, having the maximum reflectivity and a vanishing phase shift.

S1.2 Gaussian beams
In optical tweezers, light beams are tightly focused in order to increase F grad with respect to F scat .We can introduce this condition in our calculations by using Gaussian beams [50] instead of plane waves: Figure 9: Total force on a SiO 2 -Ag core-shell particle at three different wavelengths: at resonance (560 nm, black curve), at 552 nm (blue-shifted with respect resonance, blue curve) and at 566 nm (red-shifted with respect resonance, red curve).The total force is in the pN range close to the surface, due to the F enz contribution.The sinusoidal behaviour of the gradient force is visible only out of resonance (blue and red curves), while it is negligible at resonance, where, far from the surface, only the scattering force drives the total force.Close to the surface, the total force is repulsive in front of ENZ and attractive in front of Ag.Here, w 0 is the beam waist, z R = n m πw 2 0 λ is the Rayleigh range, z 0 is the position of the beam waist, R i and R r are the wave curvature radii of the incident and reflected wave, respectively, and w i (z ) and w r (z ) are the beam widths at z distance: For the sake of simplicity, we restrict ourselves to the calculation of the optical force along the beam propagation axis.The light intensity distribution I (z ) is [50]: Here, I 0 = 2P /πw 2 0 is the on-axis intensity at the waist of a beam having total power P and is a factor due to the phase shift of the beam on reflection from the surface.Thus, optical forces are calculated from Eqs. S6, S7 and S10.As above, our calculations consider a particle in water (n m =1.33) and under illumination at λ=560 nm; moreover, to evaluate the beam waist, we use the Abbe criterion w 0 = 0.5 λ N A , where the numerical aperture (NA) of the beam is NA=1.3, as in typical optical trapping experiments.The comparison between the results obtained with both plane wave and Gaussian beam on a small dielectric bead (radius 20 nm) are shown in Fig. 11.It is worth noting that when Gaussian beams are used, the beam power is reduced with respect to the plane wave case in order to maintain fixed the intensity at the beam focus.

S2 Electromagnetic scattering theory and T-matrix formalism in front of epsilonnear-zero materials
We use two different modeling approaches based on the T-matrix formalism and on finite elements methods (COMSOL), respectively.In particular, electromagnetic scattering from particles near to or deposited on a plane surface that separates two homogeneous media of different optical properties in the T-matrix formalism [42,43,44,45] can give account on the role of the different multipoles in the particle-surface interaction.Indeed, the presence of the surface can have a striking effect on the scattering pattern from the particles since the exciting field does not coincide with the incident plane wave and the observed field does not coincide with the field scattered by the particle.The field that illuminates the particles is partly or totally reflected by the surface and the reflected field contributes both to the exciting and to the observed field.Moreover, the field scattered by the particles is reflected by the interface and thus contributes to the exciting field.In other words there are multiple scattering processes between the particles and the interface.As a result, the field in the accessible half-space includes the incident field E I , the reflected field E R , the scattered field E S and, finally, the field E SR that after scattering by the particles is reflected by the surface.
The mathematical difficulties that are met in calculating the scattering pattern are due to the need that the field in the accessible half-space satisfy the boundary conditions both across the (closed) surface of the particles and across the (infinite) interface.In other words, even by assuming that we are able to impose the boundary conditions across the surface of the particle, the problem still remains of imposing the boundary conditions across the interface [43,44,45].It is possible to define the transition matrix for particles in the presence of the interface that is the starting point to calculate optical forces and torques either by direct integration of the Maxwell stress tensor or by exploiting the general expressions of optical force and torque in terms of multiple expansion [46,47,48].The focusing induces a fading of the total force with h .Note that the Gaussian beam power is reduced with respect to the plane wave case in order to maintain fixed the intensity at the beam focus.(c,d) Contour plots of the total optical force on a dielectric bead under focused Gaussian beam illumination in front of ENZ (c) and Ag (d) surfaces, as a function of the distance h from the surface and of the bead radius.The modulation due to the sinusoidal term in the gradient force is clearly visible.The total force increases at increasing bead radius, reaching the range of tens of fN in front of ENZ and hundreds of fN in front of Ag surface.
Incident and Reflected Fields.The reflection of a plane wave on a plane surface can be dealt with in general terms, i.e., without specifying whether the medium that fills the not accessible half-space is a dielectric or a metal.This information can, indeed, be supplied at the end of the algebraic manipulations.Let us thus assume that the interface is the plane z = 0 of a Cartesian frame of reference and that the half-space z > 0, which we take as the accessible half-space, is filled by a homogeneous medium of (real) refractive index n m .The half-space z < 0 is assumed to be filled by a homogeneous medium with (possibly complex) refractive index ñ .Figure 12 shows the adopted geometry.The plane wave field which propagates within the accessible half-space, is reflected by the interface into the plane wave where k I = k kI and k R = k kR are the propagation vectors of the incident and of the reflected wave, respectively, k = n m k and êI and êR are the respective unit polarization vectors.The polarization is analyzed with respect to the two pairs of unit vectors ûIη and ûRη that are parallel (η = 1) and perpendicular (η = 2) to the plane of incidence that, as usual, is defined as the plane that contains k I , k R and the z axis.Our choice of the orientation is defined by the equations with ûI2 ≡ ûR2 .In terms of the projections on the polarization basis, the incident and the reflected field can be written and In the preceding equations the incident field E I and the reflected field E R are decomposed into their components parallel and orthogonal to the plane of incidence and can be referred to each other by means of the Fresnel coefficients F η for the reflection of a plane wave with polarization along ûη .
Requiring the continuity of the normal and tangential components of the fields, the reflection condition [53] yields the equation where the Fresnel coefficients are defined as in which ϑ I is the angle between kI and the z axis, n = ñ /n m .The reflected wave can be rewritten as The incident and the reflected field, solutions of Helmholtz equation in accessible free space, can be expanded in terms of a series of spherical vector multipole fields centered on a suitable common origin, O .To ensure the regularity of the fields at the origin, we choose J-multipole fields defined in terms of spherical radial Bessel functions j l (k r ) [53,42].The result is where the incident and reflected amplitudes are respectively: Because of the reflection condition due to the presence of the surface, the incident and reflected amplitudes are not mutually independent.Infact, as the polar angles of ûR1 and ûR2 are we get In this way the amplitudes of the reflected field never need to be explicitly considered, and conveniently we can define the exciting field as the superposition of incident and reflected fields As a consequence the multipole expansion of E E can be written in a more compact form as Eηl m (39) with Scattering from a Sphere on a Plane Surface.We assume that a spherical scatterer lies entirely within the accessible half-space and is illuminated by a plane wave.Outside the scatterer the total field is where E E = E I + E R is the same as we would have if no particle were present.E S , the field scattered by the sphere, and E SR , the field that after scattering by the particle is reflected by the surface, are related to each other by the reflection condition.Their superposition represent the observed scattered field that we indicate with E Obs .The field that is scattered by a sphere that lies entirely in the accessible half-space can be expanded in a series of vector H-multipole fields that satisfy the radiation condition at infinity.The multipole fields are defined in terms of spherical radial Hankel functions [42] h l (k r ) Choosing for the scattered field the origin O within the particle, we obtain where the unknown amplitudes can be determined by applying the boundary conditions at the particle's surface.The asymptotic expression of E Sη can be written easily as follows These are the multipole fields that enter in the definition of scattering amplitude of the system.The scattered field E Sη impinges on the plane surface and, by reflection yields a reflected-scattered field in the vicinity of the surface of the particle.Thanks to the reflection rule of H-vector multipole fields [54], that proves the fields are given by a superposition of J-multipole vector fields with origin at O , we get: The quantities l l ;m can be understood as the elements of a diagonal matrix F that effects the reflection of the H-multipole fields on the plane interface giving the formal solution to the problem.Assuming that the scattering Figure 12: Geometry adopted for electromagnetic scattering from a sphere in the vicinity of a surface.
particle is a homogeneous sphere with (possibly complex) refractive index n p and radius a , also the field regular at O within the sphere can be expanded in the form The boundary conditions at the surface of the sphere between the external total field, E E +E S +E SR , and the field within the scatterer, E T , can be applied provided that the exciting field E E is referred to the center of the sphere, O .This can be done resorting to the appropriate phase factors: exp(i k I • R ) and exp(i k R • R ).For each p , l , and m, we obtain four equations among which the amplitudes of the internal field C can be easily eliminated.As a result, we get, for each m, a system of linear nonhomogeneous equations for the amplitudes where (p ) The quantities R (1) l coincide with the Mie coefficients b l and a l , respectively, for a homogeneous sphere of refractive index n p embedded into a homogeneous medium of refractive index n m .We remark that our theory can easily deal also with sphere sustaining longitudinal waves (plasmonic particles) or with radially nonhomogeneous spheres [42].
Once the amplitudes (p ) ηl m of E Sη have been calculated by solving (46), the reflected-scattered field E SRη is also determined by (50).A brief comment on the expression of the reflected-scattered field is in order.E SRη is valid only in the vicinity of the surface of the sphere as it includes multipole fields that do not satisfy the radiation condition at infinity, for this reason to get the reflected-scattered field that would be observed by an optical instrument in the far zone it is necessary to cast E SRη in its asymptotic form.At any point of the accessible halfspace, E FSRη is given by the equation [55] where for a sphere on or near the surface, this H-vector multipole fields with the origin at O can be considered as the the mirror image of the source of the original H fields. From the superposition of scattered and reflectedscattered fields, all referred to a common origin, eqs.42-46 and 49-50, we get the field with Eqs. (51)(52) lead us to the definition and derivation of the transition matrix for a scatterer in the presence of a plane interface [43,44,45].The advantages yielded by the use of the transition matrix is evident if we had to deal with the problem of a random dispersion of non spherical particles deposited on a plane surface.Moreover, the amplitudes of the observed field are the key quantities for calculating the radiation force of which we will discuss later.

S2.1 Optical force in front of a substrate
In this section we briefly recall our approach to determine the radiation force exerted by a plane waves, with a definite polarization, on a scatterer (of any shape and composition) placed in a homogeneous medium of (real) refractive index n m .We refer to the geometry sketched in Fig. 12 in which Σ is the customary laboratory frame and Σ is a frame of reference whose axes are parallel to the axes of Σ and whose origin O lies within the particle.
The vector position of O with respect to Σ is R O .The conservation laws applied to the electromagnetic scattering problem leads to the optical force acting on the particle [53,42,25]: where the integration is over the full solid angle, r is the radius of a sphere with center at R O surrounding the particle, and 〈T M 〉, the averaged Maxwell stress tensor (MST), describes the mechanical interaction of light with matter.The general expression of the MST in a medium in the Minkowski form [53,42,25] is where E is the electric field, D is the electric displacement, H is the magnetic field, B is the magnetic induction, evaluated in the frame Σ as indicated by the prime, ⊗ represents the dyadic product, and I is the dyadic unit.
We assume that all the fields are harmonic, propagating in a homogeneous, linear, and non-dispersive medium, and depend on time through the factor e −i ωt that is omitted.So, we can simplify the expression for the MST by using the complex amplitudes of the fields, E = E (r) and B = B (r), as [56,46,25]: where the fields are the superposition of the incident and of the scattered field.In presence of a plane surface that separates two homogeneous media with different refractive indexes, the role of the incident field is played by the exciting field E E = E I + E R while the superposition of E S and E SR acts like the observed field due to the presence of particle.It is possible to simplify [42] equation (53) since the dyadic products in the expression of 〈T M 〉 give a vanishing contribution to the radiative force [56,46].For these reason, the component of the radiation force along the direction characterized by the unit vector vζ turns out to be where E Obs and B Obs are the superposition of the fields scattered by the particle and the reflected-scattered fields.Obviously, since the exciting field is a plane wave, the integral (56) gets no contribution from the terms E E •E * E , and B E • B * E that, accordingly, have been omitted.At this stage, using the orthogonality properties of vector spherical harmonics through which we develop the fields, see eqs. ( 39)-( 40) and eqs.( 42)-( 50), we obtain the Borghese equations for the optical force components [48]: where where the matrix elements can be analytically valuated.We notice that F I l m .This dependence is analogous to that of the scattering cross section and of the extinction cross section, respectively, so that F (Sca) Rad ζ can be somewhat related to scattering properties of the particle, whereas F (Ext) Rad ζ can be related to its extinction.Similar considerations hold true also for the radiation torque [47,57].

S3 Finite elements methods
To evaluate the accuracy of the analytical results, we have computed the optical force on the macro particle using numerical simulations (Fig. 13).We used software package COMSOL Multiphysics 5.4 which uses finite element method (FEM) to solve Maxwell's equation and calculate the optical force.To increase the accuracy of simulation we used periodic boundary condition; however, the micro-particle radius, a , is significantly smaller than the unit cell size, L, to prevent the mutual coupling between adjacent cells (a /L =0.02) .In order to reduce the computation time the simulation is run for silicon particle with radius a = 20 nm while the wavelength of incoming plane wave is λ=560 nm.The surrounding medium is water.We used different substrates as reflecting surfaces and compared the computed force for all substrates with the analytical results.We used silver, glass, layered structure (silver and aluminum dioxide) and ENZ surfaces.The optical properties used for all surfaces are measured values at λ=560 nm.The thickness of all substrates are considerably larger than λ to mimic the semi-infinite medium.The simulation region should be meshed finely specially in three regions: i) the plasmonic layers (Ag) in the layered structure ii) in the near-field of the substrate (0 − λ/10) to capture the near-field effects on the calculated force and iii) the region surrounding the particle that the force is calculated.
To calculate the force, we used Maxwell's stress tensor.It is known that the total time-averaged force acting on any material objects can be found by calculating the integral of Maxwell's stress tensor on any surface that defines a volume containing the objects where T M is the Maxwell's stress tensor calculated based on the total electric and magnetic fields, S is the surface surrounding the volume containing the object and n is the unit vector perpendicular to the surface S .In the simulation, we chose a cylinder as a surrounding volume (Fig. 13b).In our simulation, we are interested to calculate the force in z direction consequently 〈T M,z z (r, t )〉 on top and bottom of the cylinder and 〈T z r (r, t )〉 on Figure 13: Numerical computation of the optical force: (a) Simulation region and the geometry of the problem.A polystyrene particle is placed above a substrate, the distance h is from the bottom of the particle to the substrate.The particle is surrounded with water and the incoming plane-wave is illuminated from the top and propagates in z direction.(b) Calculation of the optical force by integration of Maxwell's stress tensor on a surface of cylindrical volume surrounding the particle.The Maxwell's stress tensor consists of information about the total electric and magnetic field which includes incident field (E I ), reflected field (E R ), field scattered by the particle due to the excitation filed (E S ), and field scattered by the particle and reflected by the surface (E RS ).As we are interested on the optical force in z direction it is enough to integrate T M,z z on top and bottom surface and T M,r z on the circumferential surface.(c) [Fig.2a of the main text] Numerical computation (circles) versus analytical T-matrix calculation (solid line) and dipole approximation calculation (dots) of optical force for different distances h above different substrates (Silver, ENZ, glass).its circumference should be calculated (Fig. 13b).As the normal vector on the top and bottom of the cylinder has opposite direction, the total time-averaged force is: As integrals in Eq.S61 are calculated over the surface of the cylinder, the cylinder is meshed densely (λ/200) to avoid the numerical error.Since the field scattered by the particle causes no singularity, the height and radius of the cylinder can be as close as possible to the diameter and radius of the sphere respectively, however the smaller is the cylinder the denser should be the mesh to capture the intensity of the electric and magnetic field.To avoid computational error within reasonable computation time, the height of the cylinder is chosen to be 10 nm bigger than the sphere's diameter and the radius of the cylinder to be 5 nm bigger than the radius of the sphere.The numerical computation of total force is done for various spacing (h) of the sphere from different reflecting surfaces (Fig. 13a).h is the distance from the bottom of the sphere to the surface as shown in Fig. 13a.The numerical calculation is done for several scenarios; various spacing h , various radii a , different substrates, and particles with different polarizabilities and geometries.The comparison between analytical solution and numerical analysis for various spacing above various substrates is shown in Fig. 2a.There is a very good agreement between analytical calculation and numerical simulation.
Experimental comparison.The experimental (φ, R ) points shown as green stars in Fig. 3 are based on six different fabricated trilayer thin film stack systems (5 trilayers of Al 2 O 3 / Ag / Ge from top to bottom).The Ag layer thicknesses were in the range 10-25 nm, with a thin Ge layer (1-3 nm) underneath to ensure surface wetting.The Al 2 O 3 thicknesses were systematically varied between roughly 20 nm and 80 nm across the different systems.These stacks were deposited on a glass substrate (Corning Inc.) using electron-beam evaporation for Ge (0.5 Å/s) and Al 2 O 3 (0.3 Å/s), and thermal evaporation for Ag (0.5 Å/s).All materials were purchased from Kurt J. Lesker.The sample's ellipsometric properties (amplitude Ψ and phase difference ∆) were measured in an air superstrate using a Variable-angle, high-resolution spectroscopic ellipsometer (J. A. Woollam Co., Inc, V-VASE) for incident angles 45 • , 50 • , 55 • and wavelength range 300 − 1000 nm.Individual fits to the ellipsometric data for each system yielded best-fit results for the thicknesses and optical constants, which were then used to estimate the values of R and φ at normal incidence with λ = 560 nm shown in Fig. 3.
Effective medium theory.The effective medium theory (EMT) used to calculate the (φ, R ) curves in Fig. 1c and the curve labeled EMT in Fig. 3 takes the following form.We consider the interface between a water superstrate with refractive index n 0 and an underlying metamaterial which is a mixture of dielectric with index n d and a metal with complex index ñM = n M + i k M .If f is the filling fraction of the metal versus the dielectric, the approximate EMT permittivity of the metamaterial is given by ε EMT = (1 − f )n 2 d + f ñ 2 M .This allows us to calculate the effective refractive index n EMT and extinction coefficient k EMT as: The corresponding complex Fresnel reflection coefficient is given by: The associated reflectance R = |r EMT | 2 and phase φ = − arg r EMT .

Figure 1 :
Figure1: Geometry and optical force parameter space maps.a) Sketch of the geometry.We consider a generic particle, in principle of any shape and composition, in front of a metamaterial surface at an edge-to-edge distance, h , immersed in an external medium (e.g., water) of refractive index n m .The origin of the coordinate system is placed on the surface so that the z axis is positive in the semi-infinite space where the particle resides.A monochromatic optical wave illuminates the particle and surface at normal incidence so that the total field is the superposition of incident (E I ), reflected (E R ) and scattered fields (E S , E SR ).The resulting optical force can be either attractive (negative) or repulsive (positive).b) Near-field force (R , φ) map.We explore the range of reflectivity, R , and phase, φ, related to the reflection of the incident wave on a generic surface.The ideal ENZ surface (R = 1, φ = 0)) is in the top right corner of the map and shows a repulsive force.Here the near-field force component is calculated in the dipole approximation for a polystyrene (dielectric constant p =2.543 at 560 nm)) particle of radius a = 20 nm at a fixed edge-to-edge distance of h = 10 nm in water.c) Curves of (R , φ) for different substrates consisting of alternating metal and dielectric layers, in the thin layer limit where effective medium theory is valid.The color of each curve indicates the metal filling fraction.The zero force lines from panels d and e are superimposed as dashed lines.d-f) Total force (R , φ) maps as calculated from T-matrix methods for a dielectric particle size of a = 20 nm (d), a = 220 nm (e), and a = 1000 nm (f), respectively, and at a fixed h = 10 nm.The total force maps have a structure that is strongly dependent on particle size.This is due to the increase of the scattering force component that, for large particles, overcomes any other gradient-like force component that is dominant for nanoparticles.

Figure 2 :
Figure 2: Total optical force as a function of edge-to-edge distance, h , for polystyrene particles of different size: a) a = 20 nm), b) a = 220 nm, and c) a = 1 µm.Different approaches for the calculation of the force are compared: dipole approximation (short dots), COMSOL (circles) and T-matrix (continuous lines).Different surfaces , glass (red), silver (magenta), ENZ (blue), yield very different optomechanical interactions in terms of force amplitude, modulation with respect to h , and phase-shifts.Arrows indicate self-binding points, where particles are stably trapped in front of the surface.(d,e)Size dependence of the total optical force for fixed edge-to-edge distance, h =10 nm.For small particles the gradient force component of the partially reflected plane wave dominates, resulting in a dependence with particle size, while for large particles radiation pressure has a major contribution resulting in a negative force pushing the particle towards the surface.

#
bilayers × bilayer thickness: red labels: diel.on top blue labels: metal on top 0

Figure 3 :
Figure 3: Accessing the full range of reflectance (R ) and reflected phase (φ) via layered metamaterials.R versus φ is illustrated for light at normal incidence with wavelength λ = 560 nm reflected from the surface of a thin film stack.The thick curves of varying color are transfer matrix numerical calculations labeled n ×d , where n refers to the number of bilayers in the stack, and d the thickness of each bilayer.The total thickness nd = 500 nm is kept constant.The bilayers consist of individual Ag and Al 2 O 3 layers, with the fraction of Ag in the bilayer indicated by the metal filling fraction color.The red n × d labels correspond to systems where the dielectric is the upper layer in each bilayer (the one closest to the surface), while the blue labels are the ones where the metal is on top.Curve labeled EMT is the effective medium approximation to the system, which corresponds to n → ∞, d → 0 with n d = 500 nm.In all the above cases the superstrate is water and the substrate is glass.For comparison we show points indicating the R and φ values for a simple interface between a water superstrate and a pure material substrate (Ag, Au, Ge, TiO 2 , Al 2 O 3 , and an ideal ENZ).We also show experimental results (green stars, details in the Supp.Info.) involving a water superstrate and 5 trilayers (Al 2 O 3 /Ag/Ge from top to bottom, where Ge is present as a thin wetting layer to ensure fabrication quality).The dotted green trend line corresponds to keeping the Ag and Ge layer thicknesses fixed at 15 nm and 2.5 nm respectively, while varying Al 2 O 3 thickness from 80 nm to 20 nm (left to right).In order to compare EMT calculation with the full-wave analysis, COMSOL is used (black diamonds) to calculate R for different layered structures with various metal's filling fraction (0.4 and 0.6), layer's thickness (50 nm and 100 nm) while keeping the total thickness of the layered structure unchanged (500 nm) and also different order of material in the stack (metal on top, blue labels, and dielectric on top, red labels).

Figure 4 :
Figure4: Role of polarizability and shape on optical forces.a) COMSOL calculation of total optical force on coreshell particles based on SiO 2 and Ag as a function of the distance from the metamaterial surface.Different ratios between core radius a 1 and total particle radius a tot =20 nm have been considered.Moreover, both materials have been considered as a core.b) Maximum force found in each core-shell structure considered, as pointed out by the dashed blue line in a).c) Extinction spectrum of the SiO 2 -Ag core-shell particle (total radius a tot =20 nm and core radius a 1 =16.1 nm) in water.d) (R , φ) contour plot of the total optical force for the core-shell particle at h =10 nm from the surface.The ENZ and Ag surfaces used for the calculation of the optical forces in DA approximation are shown as circles.e) Extinction spectra of Ag prolate ellipsoid in water oriented with the long axis parallel to the field (black solid line) and oriented with the short axis parallel to the field (red dashed line).The resonances relative to the long and short axes are indicated.f) (R , φ) contour plot of the total optical force on the Ag ellipsoid at h =10 nm distance from the surface.In the calculation, the spheroid is aligned with the long axis in the direction of the wave polarization.The ENZ and Ag surfaces used for the calculation of the optical forces in DA approximation are shown as circles.The optical force is in the order of tens of pN in front of ENZ (repulsive) and Ag (attractive).The force can be close to 200 pN if the spheroid is in front of an ideal ENZ surface, having R=1 and φ =0.

Figure 5 :
Figure 5: Optical forces on a) 20 nm, b) 50 nm and c) 100 nm radius dielectric beads.The beads are in front of ENZ (blue curve), Ag (magenta curve) or glass (red curve) surfaces.Note how for smaller particles the dominant contribution to the optical force comes from the gradient force, while for larger particles the greater scattering force shifts downwards the force modulation resulting from the interference between incident and reflected field.In d) and e), 3D plots of the total force as a function of the particle radius and of the distance from ENZ (d) and Ag (e) surfaces.

Figure 7 :
Figure 7: (a-c) Total force on spherical beads made by ENZ material.The forces are shown for 20 nm radius (a), 50 nm radius (b) and 100 nm radius (c).ENZ (blue curve), Ag (magenta) and glass (red) surfaces are considered for the calculation.(d,e) 3D plots of the total force as a function of the ENZ particle radius and of the distance from the ENZ (d) surface and Ag (e) surface.

Figure 8 :
Figure 8: (a) Extinction spectrum of the SiO 2 -Ag core-shell particle (total radius a t o t =20 nm and core radius a 1 =16.1 nm) in water.(b,c) Total optical force of the core-shell particle at fixed distance h =10 nm from ENZ (b) and Ag (c) surfaces as a function of the a 1 to a t o t ratio.(d,e)Contour plots of the optical force with respect to the a 1 to a t o t ratio and the distance h from the surface.The force on the core-shell particle is in the pN range only at short distances from the surfaces and repulsive in front of ENZ (d) while attractive (e) in front of Ag.

Figure 10 :
Figure 10: (a) Extinction spectra of Ag prolate ellipsoid in water oriented with the long axis parallel to the field (black solid line) and oriented with the short axis parallel to the field (red dashed line).The resonances relative to the long and short axes are indicated.(b) Contour plot of the total optical force on the Ag ellipsoid at h =10 nm distance from a surface as a function of the surface reflectivity R and phase shift φ.In the calculation, the spheroid is aligned with the long axis in the direction of the wave polarization.The ENZ and Ag surfaces used for the calculation of the optical forces in DA approximation are shown.The optical force is in the order of tens of pN in front of ENZ (repulsive) and Ag (attractive).The force can be close to 200 pN if the spheroid is in front of an ideal ENZ surface, having R=1 and φ=0.

Figure 11 :
Figure 11: (a,b) Total optical force on a 20nm dielectric bead under plane wave (blue curves) and focused Gaussian beam (red curves), in front of ENZ (a) and Ag (b) surfaces, as a function of the distance h from the surface.The focusing induces a fading of the total force with h .Note that the Gaussian beam power is reduced with respect to the plane wave case in order to maintain fixed the intensity at the beam focus.(c,d) Contour plots of the total optical force on a dielectric bead under focused Gaussian beam illumination in front of ENZ (c) and Ag (d) surfaces, as a function of the distance h from the surface and of the bead radius.The modulation due to the sinusoidal term in the gradient force is clearly visible.The total force increases at increasing bead radius, reaching the range of tens of fN in front of ENZ and hundreds of fN in front of Ag surface.
Rad ζ depends on the amplitudes A (p ) l m of the scattered field only, whereas F (Ext) Rad ζ depends jointly on the amplitudes of the scattered field A (p ) l m and on those of the incident field W (p )

S
T o p 〈T M,z z (r, t )〉d S − S B o t t o m 〈T M,z z (r, t )〉d S + S C i r c u m〈T M,r z (r, t )〉d S(61)