Optical forces in nanorod metamaterial

Optomechanical manipulation of micro and nano-scale objects with laser beams finds use in a large span of multidisciplinary applications. Auxiliary nanostructuring could substantially improve performances of classical optical tweezers by means of spatial localization of objects and intensity required for trapping. Here we investigate a three-dimensional nanorod metamaterial platform, serving as an auxiliary tool for the optical manipulation, able to support and control near-field interactions and generate both steep and flat optical potential profiles. It was shown that the ‘topological transition’ from the elliptic to hyperbolic dispersion regime of the metamaterial, usually having a significant impact on various light-matter interaction processes, does not strongly affect the distribution of optical forces in the metamaterial. This effect is explained by the predominant near-fields contributions of the nanostructure to optomechanical interactions. Semi-analytical model, approximating the finite size nanoparticle by a point dipole and neglecting the mutual re-scattering between the particle and nanorod array, was found to be in a good agreement with full-wave numerical simulation. In-plane (perpendicular to the rods) trapping regime, saddle equilibrium points and optical puling forces (directed along the rods towards the light source), acting on a particle situated inside or at the nearby the metamaterial, were found.

limitations on the flexibility of optical manipulation by reducing potential degrees of freedom, available for optomechanical control. On the other hand, three-dimensional artificially created nanostructures or metamaterials could provide additional benefits and flexibility by configuring near-field interactions in large volumes 15 .
Hyperbolic metamaterials 16 are one class of artificially created electromagnetic structures, capable to enhance efficiencies of various light-matter interaction processes owning to the unusual dispersion regimes of eigenmodes, supported by the structure -namely its hyperbolic dispersion. Among various possible designs of this type of metamaterials it is worth mentioning composites made of vertically aligned nanorods 17,18 , periodic metal dielectric layers 19 and semiconductor quantum structures 20,21 . While the far-field interactions of waves with hyperbolic composites were proven to be well characterized in terms of the effective medium approximation, this description could be questioned if near-field mediated processes are involved 22 .
The general criterion for estimating impact of near-field contributions to an interaction is based on comparison of k-vector spectra with reciprocal vector of the metamaterial lattice. For example, scattering from objects within hyperbolic metamaterials involves consideration of the near-field effects 23 .
Analysis of near-and far-field contributions to optical force, acting on objects embedded in the nanorod metamaterial is the central topic of the manuscript. In particular, optical forces, acting on nano-sized spherical particle embedded inside the metamaterial assembly, are investigated both numerically and by using a semi-analytical approach, considering the finite size nanoparticle as a point dipole and neglecting re-scattering between the particle and nanorod array. The impact of the finite structure of the metamaterial unit cell and the relative arrangement of the particle in respect to it was analyzed as a function of the system's geometry and frequency of incident illumination. A semi-analytical approach based on dipole near-field interaction is developed and shown to be in a good agreement with results of the full-wave numerical analysis. The interplay between near-and far-field effects in the context of effective medium approximation is discussed.

Nanorod Metamaterial: far-field characteristics
The geometry under investigation is schematically represented in Fig. 1a -it shows an array of vertically aligned gold nanorods and a gold nanoparticle placed inside it. Material parameters of the constitutive elements were taken from widely used sources 24 . The parameters of the structure are indicated in the figure caption. While the nanorods in this model are situated in vacuum, substrate effects and host material filling the space between the rods could be taken into account straightforwardly. Similar structures have already found use in various multidisciplinary applications, among them bio-sensing 25 , enhancement of nonlinearities 26 , acoustic waves detection 27 , thin optical elements 18 and others. The key properties of this auxiliary nanostructure leading to enhanced performance are large surface area and unusual collective optical response of the system, enabling control over both far-and near-field interactions. Hence, investigation of optical forces, mediated by nanorod metamaterials, has a profound potential interest.
Far-field interactions between electromagnetic waves and metamaterials, under certain circumstances, can be described within the effective medium approximation. The main idea of this homogenization procedure is to average the electromagnetic field over a unit cell of a structure. Therefore, the field inside the structure is assumed to be uniform. In the context of optical forces, as it will appear in the next section, the non-uniformities play a major role and, in fact, predefine the spatial structure of optical potentials. Recently, a phenomenological approach taking into account the finite size of the metamaterial unit cell was proposed 28 . Inclusion of a depolarization volume around optically manipulated particles enabled investigations of far-field contributions to optical forces. However, near-field interactions, being strongly dependent on a specific metamaterial design, were not included explicitly.
One of the key properties of hyperbolic metamaterials, making them attractive for electromagnetic applications, is their unique ability to support an unusual regime of dispersion caused by having permittivity tensor components of opposite signs ε ε ( < ) ⊥ 0 . An immediate implication of this hyperbolic dispersion regime is the high density of photonic states, available for both emission and scattering 29,30 .
The effective permittivity tensor of nanorod metamaterial is given by: Here ε ⊥ and ε are effective permittivities perpendicular and along the wires, respectively. Dispersion of the tensor components for the structure under consideration was calculated with the approach developed in 31 . The transition between elliptic and hyperbolic dispersion regimes occurs at the wavelength around 523 nm (see Fig. 1b). The transition point is called epsilon-near-zero (ENZ) regime, as the real part of the permittivity along the rods is vanishing, if the spatial dispersion effects are ignored. The high density of photonic states as well as the strong scattering emerges in the hyperbolic and ENZ regimes. The wavelength of the external illumination, exploited for optical manipulation in the subsequent investigations, is chosen around this ENZ point in order to distinguish between various dispersion regimes and their impact on optical forces.

Optical force distribution
Numerical model. The distribution of the optical forces, acting on the gold nanoparticle placed inside the wire medium is analyzed hereafter. The hyperbolic, ENZ, and elliptic dispersion regimes of the bulk metamaterial and their impacts on optical forces are compared and discussed.
In the first case, normal incidence scenario is considered -the illumination is chosen to be linearly polarized along the y-axis and it propagates along the z-axis (ϕ = 0°, see Fig. 1a). Oblique incidence with ϕ = 45° will be investigated hereafter. Full 3D numerical analysis, based on finite elements method 32 , is performed in order to calculate self-consistent electromagnetic fields in the system. Consequently, optical forces acting on the nanoparticle are calculated by integrating the Maxwell's stress tensor components over an imaginary spherical surface surrounding the nanoparticle.
The presence of a single nanoparticle breaks the inherent translation symmetry of the initial metamaterial geometry. In order to overcome the computation complexity of large systems modeling, Floquet periodical boundary conditions were imposed on finite size geometries. This type of model corresponds to a periodic system with variable unit cell, which consists of a square array of nanorods and the nanoparticle. If the electromagnetic coupling between the particles in adjacent cells is minor, this type of analysis recovers the behavior of the infinite system with a single particle.
The numerical procedure is as follows: the number of rods in the unit cell is increased gradually and the convergence of a certain quantity (optical forces in our case) is checked. Recently, a similar approach was applied in studies of the Purcell effect in nanorod 33 and wire 34 metamaterials. Square unit cells containing 4, 9, and 16 nanorods were considered and the convergence of optical force values at different points of the metamaterial volume was checked. A unit cell of 4 nanorods (the smallest one) was shown to predict the behavior of an infinite array within the accuracy of several percent. All the subsequent results were obtained for this size of the unit cell. The direct consequence of this calculation is that (i) only nearest neighbor rods define the value of optical force and (ii) nanoparticles in different unit cells almostly do not interact with each other.
It should be noted, however, that the collective macroscopic behavior of the array is taken into account by imposing periodical Floquet boundary conditions. Lateral force component. All the subsequent calculations were done for a particle of 10 nm in diameter. The optical force F, in the most general case, has three non-zero components (F x , F y , F z ). The lateral force F ⊥ = (F x , F y , 0) will be analyzed first. Values of optical forces are normalized to the intensity of the incident wave and volume of the particle in order to perform direct comparisons with other optical manipulation schemes.
The resulting normalized forces at the cut-plane z = 10 nm, calculated for wavelengths λ = 450 and 600 nm, corresponding to the elliptic and hyperbolic dispersion regimes respectively, are shown in Fig. 2(a). Numerical simulation shows that spatial distribution of both electric field and optical force for different wavelengths of excitation (different dispersion regimes) are qualitatively similar and the quantitative difference is shown with different color bars on top panel in Fig. 2. The qualitative explanation of such a behavior is two fold: (i) small size of the sphere and nanorods (in comparison with wavelength) enables considering those structures as point dipoles; (ii) the same material of the sphere and nanorods providing similar frequency behaviour of their polarizability.
Force maps at various cut-planes (with different z-coordinate) show qualitatively similar behavior too. The reasoning is as follows: at normal incidence, lateral optical force is mainly determined by first (gradient) term in Eq. (2). Within a good approximation |E| 2 can be factorized as f(z)g(x, y), where f(z) is Fabry-Perot envelope function and g(x, y) is the in-plane field distribution. Therefore, the lateral force distribution is similar for different cut-planes while its absolute value is modulated by Fabry-Perot envelope function.
It should be noted, that the values and directions of forces are attributed to the geometrical center of the particle, hence certain regions (dark blue shells around nanowires with thickness equal to the radius of nanoparticle) on Fig. 2(a,c) are blank, as this center cannot approach the boundaries of the rods.
It can be seen from Fig. 2(a) that the force distribution has saddle points at the center of the unit cell and at its edges between the rods. These places correspond to the the saddle points of electromagnetic field magnitude distribution [ Fig. 2(d)], and, consequently, to the unstable equilibrium positions of the nanoparticle. Some peculiarities in the optical force distribution appear on cut-planes near the edges of the nanorods (z = 350 nm and z = 0 nm), but they are attributed to the longitudinal (z-component) force component and will be discussed further. where the perturbation of the field by the particle is neglected. Panel (d) shows the numerical simulation of the electric field magnitude distribution without the particle. Panel (a) should be compared with (c), while (b) with (d). Dark blue shells around nanorods on panels (a), (c) have the width equal to particle's radius. Optical forces are not calculated at those areas, as the nanoparticle's center cannot approach the nanorods that close. Upper and lower scales of the color bars correspond to the hyperbolic (λ = 600 nm) and elliptic (λ = 450 nm) dispersion regimes of the metamaterial, respectively. Electric field amplitude of the incident wave is 1 V/m. Electric field E of the incident wave is parallel to the y-axis. Wavevector of the incident wave k = (0, 0, k 0 ). The similarity of the spatial distribution of the forces at hyperbolic and elliptic dispersion regimes results from the dominating near-field coupling between the nanorods and the particle. Figure 2(b) shows the magnitude distribution of total electric field , while the arrows show its direction. One can see that the field map is formed by electrical dipoles induced on the rods and the particle by the incident wave. Orientations of the dipoles nearly coincide with the polarization of the incident wave. Minor deviations from the above description are related to the higher multipole contribution and the interaction between the particle and the nanorods.
Lateral force distribution analysis can be provided with the following semi-analytical approach. First, the total electric field distribution in the nanorod array under external incident wave without the particle is calculated numerically with the periodic boundary conditions applied. Results of the simulation are shown in Fig. 2(d). The knowledge of the spatial distribution of the electric field magnitude enables calculation of optical forces with two assumptions: (i) the nanoparticle is represented by a structureless point electric dipole with a moment μ (ii) the dipole is assumed to act as a small perturbation to the fields of the standalone metamaterial. This means, that only collective scattering properties of the nanorod array were taken into account, while the mutual re-scattering of the field between the particle and nanorods was neglected. Comparison between Fig. 2(b,d) verifies this approximation -both the structure and values of the field magnitude are similar.
The time averaged optical force acting on the point dipole is given by 7 : 2 0 0 where α = α′ + iα′ ′ is the complex particle's polarizability. The polarizability of the spherical particle is given by 7 : The resulting optical force map, calculated using the dipolar approximation [Eq.
(2)] appears on Fig. 2(c). The arrows show the direction of the force at corresponding points. The color pattern corresponds to the absolute value of the force. The remarkable similarities between Maxwell's stress tensor calculations [ Fig. 2(a)] and the approximate analytical model [ Fig. 2(c)] suggest the validity of the dipolar model and highlights the impact of near-fields on the optical force. It should be noted, however, that overall values of optical forces, calculated within those approaches, have about 20% difference, which is related to the finite size of the particle and the mutual field re-scattering between the particle and nanorods.
Distribution of optical force in the case of oblique incidence (ϕ = 45°, see Fig. 1a) obtained with full numerical simulations and semi-analytical approach is shown in Fig. 3. One can see that in contrast to the case of normal incidence electric field distribution around the nanorods is asymmetric [ Fig. 3(b,d)]. It results in asymmetry of optical force distribution. The asymmetry can be explained by excitation of non-dipole modes in the nanorods at oblique illumination 35 . Vertical force component. The distribution of the lateral optical force component (perpendicular to the nanorods) was analyzed in the previous section. Longitudinal force component F z (parallel to the nanorods) is analyzed here.
As it was already mentioned, the homogenization procedure averages the near-fields over the unit cell, hence, it is inapplicable for estimation of gradient optical force in the lateral plane. Nevertheless, the field distribution along the z-axis can be roughly estimated considering the slab of the nanorod metamaterial as a Fabry-Perot resonator in z-direction 33 . Therefore, it is reasonable to expect a standing wave in the slab (along the z-axis) and maxima of the electric field resulting in in-plain trapping of the particle. The results of the numerical calculation suggest the validity of this hypothesis. Estimation of the field maxima position deeply inside the slab can be provided by the effective medium approximation 36,37 but a more detailed analysis, that takes into account boundary effects, demands numerical simulation.
The profiles of the total electric field magnitude along the line parallel to the rods and passing through the point x = 30 nm and y = 5 nm [see inset in Fig. 4(a)] calculated without nanoparticle for the elliptic (λ = 450 nm), ENZ (λ = 523 nm), and hyperbolic (λ = 600 nm) regimes are shown in Fig. 4(a). The insets show distribution of the total electric field magnitude in xz-plane passing through y = 5 nm. One can see that electric field distribution along the z-axis strongly depends on the wavelength of the incident wave. In the hyperbolic regime (λ = 600 nm), three distinct field maxima are observed -one inside the slab and two in the vicinity of its boundaries. In the elliptic (λ = 450 nm) and ENZ (λ = 523 nm) regimes, field decays inside the metamaterial and weak oscillations do not possess sharp field maxima. Additional contribution to those differences (apart from the interplay of dispersion regime and geometry, namely Fabry-Perot conditions) comes from a strong wavelength dependence of losses in gold 24 : Optical losses cause the reduction in quality factors of the modes, smearing out the sharp peaks, as could be seen in the case of elliptic dispersion.
Distributions of the z-component of the optical force along the nanorod for elliptic (λ = 450 nm), ENZ (λ = 523 nm) and hyperbolic (λ = 600 nm) regimes are shown in Fig. 4(b). Positions of the stable trapping in transverse planes are marked with arrows (note, that the force derivative should be negative in order to obtain a stable equilibrium). Shaded areas on the figure show the regions within the metamaterial, where the optical force component F z is directed towards the light source (for λ = 600 nm). Optical pulling forces or optical attraction gained considerable attention over the last decade, as it provides additional flexible degree of freedom in optical manipulation 38 . The hyperbolic regime supports three regions of optical attraction, while the elliptic and ENZ have only one, as could be seen in Fig. 4(b). This occurrence could be understood as follows: in the elliptic regime both weak gradient of the electric field magnitude [see Fig. 4(a)] and high material losses of the particle result in the domination of radiation pressure [second term in Eq. (2)] over the gradient force. The radiation pressure is co-directional with the Poynting vector of the incident radiation, so the optical attraction cannot be obtained in this case. Nevertheless, the first term of Eq. (2) overcomes the second one in the vicinity of the nanorod's edge where strong gradient of the electric field intensity is observed [see Fig. 4(a)]. In the hyperbolic regime, on the other hand, there are several regions where the optical force component is directed to the light source -that's the result of high quality factor Fabry-Perot modes and dominating real part of the particle's polarizability.
For the hyperbolic regime (λ = 600 nm), optical potential in the vicinity of nanorod's edge is stronger than in elliptical and ENZ regimes. For example, optical traping potential of 26 meV for 5nm radius particle can be achieved with electric field intensity ~10 7 V/m that corresponds to focusing of 1 W beam to 3 μm spot in diameter.
As a separate case, the particle situated over the metamaterial slab will be considered next. This scenario describes the case where the metamaterial is used as a substrate for advanced optical manipulation. Results of numerical studies appear in Fig. 5, showing the distribution of the vertical optical force acting on the nanoparticle in the lateral plane of z = 10 nm above the nanorods. It could be seen, that  the maximal attraction force on the particle emerges in the vicinity of nanorods edges (the sample is illuminated from below -see Fig. 5). Arrows indicate the direction of the optical force. Color pattern shows the distribution of the optical force component parallel to the nanorods (F z ). Solid white lines correspond to F z = 0. Remarkable behaviour of forces above the metamaterial substrate could suggest the later as an auxiliary nanostructure or metasurface, providing additional flexibility in optical manipulation.

Conclusion
In this work, comprehensive analysis of the optical forces acting on a metal nanoparticle placed inside or in the vicinity of three-dimensional nanorod metamaterial slab was performed. Numerical simulations of finite size square unit cells with periodical Floquet boundary conditions enable to take into account all collective effects in the metamaterial and estimate optical forces on small particles. Unit cells containing 4, 9, and 16 nanorods were analyzed and the convergence of the optical forces for different positions of the particle was checked. It was shown that the smallest unit cell already reproduces the effect of optical forces on a particle, situated within the infinite metamaterial. Therefore, only four neighboring nanorods nearest to the particle make the dominant contribution to the optical forces. This statement has been confirmed with the developed semi-analytical model which neglects the particle's interior and the re-scattering effects between the particle and nanorods. Furthermore, it was shown that the 'topological transition' from the elliptic to hyperbolic dispersion regime of the metamaterial, usually having an impact on various light-matter interaction processes, is less important for optical forces.
In-plane optical trapping and optical pulling forces were observed. The comprehensive numerical modeling enables estimation of optical forces values, normalized to incident power and particle's volume. Values as high as 2.3 × 10 3 pN/W/nm for both lateral and optical pulling forces were predicted. Those results overcome other reported values 39,40 .
The remarkable structure of predicted optomechanical interactions (in particular pulling forces), mediated by the metamaterial, makes the later to be a promising platform for large span of multidisciplinary applications, involving demands for precise nanoscale mechanical manipulation, including trapping sorting, mixing and more.