## Abstract

It is well known that incident photons carrying momentum ℏ**k** exert a positive photon pressure. But if light is impinging from a negative refractive medium in which ℏ**k** is directed toward the source of radiation, should light exert a photon “tension” instead of a photon pressure? Using an ab initio method that takes the underlying microstructure of a material into account, we find that when an electromagnetic wave propagates from one material into another, the electromagnetic stress at the boundary is, in fact, indeterminate if only the macroscopic parameters are specified. Light can either pull or push the boundary, depending not only on the macroscopic parameters but also on the microscopic lattice structure of the polarizable units that constitute the medium. Within the context of an effective-medium approach, the lattice effect is attributed to electrostriction and magnetostriction, which can be accounted for by the Helmholtz stress tensor if we use the macroscopic fields to calculate the boundary optical stress.

- Photon pressure
- electromagnetic stress
- optical pulling force
- metamaterial
- negative-refractive-index material
- effective medium

## INTRODUCTION

Electromagnetic (EM) wave or light-induced forces have been extensively studied for decades, generating fruitful applications for manipulating microscopic particles, biological materials, and macroscopic mechanical components (*1*–*12*). It is well known that a plane wave in vacuum always exerts a photon pressure on an object in the direction of momentum flux. With the advent of metamaterials (*13*–*15*), which can have negative refractive indices (NRIs), can light pressure become light tension when light propagates in such a material (*13*)?

For example, consider the situation shown in Fig. 1A, which depicts a boundary separating two materials, each specified by their own macroscopic permittivity (ε) and permeability (μ). How would the optical stress (for example, pressure or tension) at the boundary depend on the ε and μ of the media on either side of the boundary? The answer to such a question is much more difficult to find than it might seem. Even if we assume that the media are truly isotropic and homogeneous, we run into difficulties such as the correct interpretation of EM momentum density (the Abraham/Minkowski controversy) and the related issue of what type of EM stress tensors should be used in calculations (*16*–*25*). Another major problem is that metamaterials are always inhomogeneous. All metamaterials with unusual optical properties have an underlying microstructure, usually in the form of arrays of subwavelength resonators. Here, we will adopt an “ab initio” approach and we calculate the boundary stress with a medium treated as an array of coupled scatterers embedded in vacuum with a lattice constant much less than the wavelength. This approach kills two birds with one stone. First, it takes the effect of the microstructure fully into account. Second, it has the advantage that the time-averaged force for a time harmonic field acting on each scatterer (and hence the force density everywhere) can be determined unambiguously using the Maxwell tensor approach as each scattering unit is embedded in air. The boundary stress can thus be determined. We will see that the boundary optical stress is, in fact, indeterminate if we only specify the macroscopic parameters as in Fig. 1A. Light can push or pull on the boundary, depending on the details of the microstructure, although the total force is always positive. We will see that if we want to consider the optical stress problem within the context of an effective-medium approach, we can use the Helmholtz stress tensor, which carries the information of microstructure in the electrostriction and magnetostriction terms. These terms do not appear in standard effective-medium theories but can be calculated with some additional effort (*26*).

## RESULTS

### EM stress at the boundary of the microscopic model system

Let us start with the simplest possible configuration as shown in Fig. 1A, where a boundary separates two semi-infinite materials of different optical properties as specified by their respective permittivity and permeability. An EM plane wave (with electric field amplitude of 1 V/m) propagates from the left-hand side (LHS) to the right-hand side (RHS) and induces a local optical stress on the boundary, and we would like to know whether the light exerts a pressure or a tension if the media have unusual (say NRI) optical properties. Because metamaterials with strange values of refractive indices derive their optical properties from an underlying lattice of scatterers or resonators, a more realistic representation of the boundary region would be the schematic picture shown in Fig. 1B, in which the discreteness of the microscopic sublattice is explicitly shown. In our calculations, we shall assume that the lattice consists of subwavelength polarizable units schematically shown in Fig. 1B, where points of different colors represent units of different materials. For simplicity, we consider a two-dimensional (2D) configuration in which each polarizable unit is characterized by electric and magnetic polarizabilities α^{e} and α^{m}. In the case of *E*-*z* polarization (electric field along *z* direction), each unit then represents an out-plane electric monopole **p** = α^{e}**E** and an in-plane magnetic dipole **m** = α^{m}**H**. The polarizable units form a lattice (for example, square), with the size of each unit cell being *a*. In the limit of λ ≫ *a*, the microscopic lattice can be treated as an effective medium with relative permittivity ε_{eff} = 1 + α^{e}/(*A*ε_{0}) and relative permeability μ_{eff} = (2*A* + α^{m})/(2*A* − α^{m}), where *A* denotes the area occupied by a single unit cell (*27*). The EM field within the microscopic lattice system can be determined using the multiple scattering theory detailed in Materials and Methods. The force exerted on each polarizable unit can be unambiguously calculated by integrating the time-averaged Maxwell stress tensor (*28*) over a closed loop as shown in Fig. 1B. We note that the loop resides in vacuum, and hence, the Maxwell stress tensor approach is “exact.”

In the usual sense of effective-medium or homogenization theories as accepted by the optics and metamaterial communities, the microscopic model system depicted in Fig. 1B is macroscopically equal to Fig. 1A as long as λ ≫ *a*. However, some subtleties still arise even in that ideal limit. The boundary in Fig. 1A is a geometric line (the black line), but when the microstructure is taken into consideration as in Fig. 1B, where exactly should we draw the boundary? The homogeneous system depicted in Fig. 1A is translational invariant along the direction of the boundary, whereas the discrete system is not and the calculated result will depend on the relative alignment of the left and right lattice structures. We will use *d* to denote the relative shift distance between the two sublattices, with *d* = 0 corresponding to perfect alignment. It is clear in Fig. 1B that the boundary encloses a finite region of polarizable units from both sides. To define the boundary stress, this region has to be determined, and it can be done in the following way. We will focus on discrete systems with impedance-matched (μ_{eff}/ε_{eff} the same for the left and right sublattices) configurations and leave the impedance-mismatched cases for later discussion. We first consider an example and choose the values of α^{e} and α^{m} so that the effective permittivity and permeability are ε_{eff,1} = μ_{eff,1} = −ε_{eff,2} = −μ_{eff,2} = −2 and the wavelength is fixed to be λ = 100*a* unless otherwise specified. The results obtained with the multiple scattering theory are shown in Fig. 1D, where each open circle denotes the calculated force acting on the corresponding polarizable units, with *x* = 0 marking the middle line separating the two lattices. Four configurations with different relative lattice shift *d* are calculated. We see that for all the configurations, the forces are slowly vanishing except near the interface region marked by blue color. This is because for impedance-matched configurations, momentum transfer happens at the interface region and it is then natural for us to define the domain within the blue region as the “boundary.” The boundary stress is then defined as the total force per unit length acting on this region. We note that although the force distribution within the interface region does vary with *d*, the total force in the blue domain (that is, boundary stress) adds up to be the same constant as shown in the inserted panel in Fig. 1D. This means that the “boundary stress” is well defined to the extent that it does localize near the geometric boundary and it does not depend on how the microlattices are aligned relative to each other.

### Indeterminate stress

We note that the relationship between a macroscopic effective medium (Fig. 1A) and the underlying microstructure (Fig. 1B) is not a one-to-one mapping. For example, the underlying lattice does not need to be a square lattice as shown in Fig. 1B. We can find a triangular lattice that gives exactly the same effective ε_{eff} and μ_{eff}. Here, we calculate the boundary stress for two different types of microscopic lattice structures (*C*_{4v} and *C*_{6v}) that correspond to the same effective macroscopic system in Fig. 1A. The results are summarized in Fig. 2 (A to D), in which the dots show the calculated boundary stress for different combinations of lattice types as illustrated in the insets. Figure 2A shows the results for the boundary formed by two square sublattices, where we fix *n*_{eff,1} (*n*_{eff,1} = ε_{eff,1} = μ_{eff,1}) and vary *n*_{eff,2} (*n*_{eff,2} = ε_{eff,2} = μ_{eff,2}) to see how the stress changes. We see that when *n*_{eff,1} = 1 (corresponding to air), the boundary stress is always negative, meaning that a plane wave propagating in air always pulls the surface of a medium whose microscopic polarizable units form a square lattice. However, for *n*_{eff,1} = −1, the boundary stress can be either positive or negative. Figure 2B shows that the results are very different when the two sublattices are triangular. In contrast to the square-lattice case, a plane wave in air (*n*_{eff,1} = 1) always pushes the boundary. The σ_{x} versus *n*_{eff,2} relationship is completely different for the two kinds of lattice, even though they are purposely designed to have exactly the same effective ε_{eff} and μ_{eff}. We can also combine a square lattice with a triangular lattice. The results are shown in Fig. 2C, with the square lattice on the left and the triangular lattice on the right, and Fig. 2D shows the results when the left/right positions are exchanged. We see that the boundary stress changes sign when the positions of the two sublattices are exchanged. In addition, we see that the hybrid systems give different boundary stress compared with the previous two cases shown in Fig. 2 (A and B). The four cases considered in Fig. 2 (A to D) have very different boundary stresses, even though they have, by design, the same effective macroscopic permittivity and permeability. This can be illustrated in Fig. 2E, where the case with *n*_{eff,1} = −1, *n*_{eff,2} = 4 is taken as an example. The four types of lattice correspond to the same macroscopic effective system as depicted in the center panel. And yet, the stress exerted on the boundary is different in sign (indicated by the arrow) and magnitude in the four cases. These results show that the light-induced boundary stress for a metamaterial system depicted in Fig. 1A is indeterminate even if the macroscopic permittivities and permeabilities are completely specified. The results depend on microscopic details. We will see in the second part of the paper that this rather counterintuitive phenomenon is due to the existence of electrostriction and magnetostriction effects.

If the sublattices’ impedances do not match each other, the wave reflection at the boundary will result in a field gradient that induces a gradient force on each polarizable unit. In that case, the boundary stress is not necessarily well defined because the optical force is not localized near the boundary and there is no unique way to attribute part of the force to the boundary. However, if the LHS sublattice is air (see the insets of Fig. 3), the boundary stress is still well defined even if the surface impedance of the RHS half space is not matched with air. In this case, the reflected field on the LHS does not affect the boundary stress. Here, we consider *n*_{eff,2} < 0 with ε_{eff,2} ≠ μ_{eff,2} so that its impedance does not match that of air on the LHS. The dots in Fig. 3 (A and B) show the boundary stress for the square lattice and triangle lattice, respectively, where we fix ε_{eff,2} = −1. The stress is seen to pull in the square lattice but push in the triangle lattice. If we fix μ_{eff,2} = −1 and change ε_{eff,2}, the results are similar (Fig. 3, C and D). Hence, even in the general case of an impedance-mismatched boundary, whether the light exerts a pressure or a tension on the boundary depends on microscopic details, even though the macroscopic effective-medium parameters are the same.

We note that if the microscopic system is a slab of finite thickness, the total force induced on the system by an incident plane wave is the same for different lattice types, as long as their effective permittivity and permeability have the same value. Take for example the configuration depicted in Fig. 1B, and suppose that the medium on both sides of the lattice (yellow- and green-colored regions) is air. This system corresponds to a finite-sized slab. Again, we calculate the force acting on each polarizable unit in this case. The dots in Fig. 4 (A and B) show the internal force distributions for a square-lattice and a triangle-lattice system, respectively. We have set ε_{eff,1} = μ_{eff,1} = −1 and ε_{eff,2} = μ_{eff,2} = 4. We see that only the polarizable units at the boundary feel a force due to impedance-matched configuration, and the boundary force is very different for the two kinds of lattice. In Fig. 4C, we show that the cumulative force summed starting from the left-most layer of the slab and the final value (ending on the RHS and marked by a blue dot) gives the total force acting on the slab. The total forces (blue dot) are equal for both cases of Fig. 4 (A and B), whereas the forces on individual boundaries are different. Figure 4 (D and E) shows an impedance-mismatched case with ε_{eff,1} = μ_{eff,1} = −1 and ε_{eff,2} = 1, μ_{eff,2} = 4, where the optical forces extend into the bulk due to the field gradient arising from Fabry-Perot reflections. We see that the force distributions in the bulk and at the interface regions are very different for different lattices. In this case, the boundary stress is hard to define, but we see from Fig. 4F that the total force (end value of the lines in Fig. 4F marked by the blue dot) again adds up to the same value for the two lattice configurations. These results show that the total forces acting on an object are the same as long as the macroscopic parameters are the same, but the internal force can be very different in sign and magnitude for different underlying lattices. We further consider the total force exerting on a semi-infinite reflecting metallic material when the plane wave is incident from an NRI material. The metallic material on the right of the boundary is modeled by a microscopic lattice structure, which gives a negative effective permittivity and μ_{eff} = 1 in the effective-medium limit so that the field decays exponentially into the bulk of the metallic material. The results for a square lattice and a triangle lattice are shown in Fig. 4 (G and H), respectively. We set ε_{eff,1} = μ_{eff,1} = ± 1 (that is, air and NRI) and ε_{eff,2} = −4, μ_{eff,2} = 1. We note that no matter the reflecting metal is forming a boundary with air or NRI, the total forces (sum of the forces acting on all points with *x* > 0) exerting on the metallic material are always positive for both the square lattice and the triangle lattice. These results show that light always pushes on a metallic material, independent of whether it is incident from a positive or negative index medium.

We now consider the situation in which the underlying lattice is amorphous as shown in Fig. 5A, where both sublattices consist of randomly positioned polarizable units. In this case, we use a numerical solver COMSOL (*29*) to calculate the fields and hence the boundary stress. In the simulations, each polarizable unit is defined by a cylinder with radius *r*, relative permittivity ε_{c}, and relative permeability μ_{c}. In the long wavelength limit and *E*-*z* polarization, the system is essentially a collection of in-plane magnetic dipoles and out-of-plane electric monopoles with effective parameters given by ε_{eff} = 1 + ρ(ε_{c} − 1), μ_{eff} = −1 + 2/(1 − ρ*M*), where *M* = (μ_{c} − 1)/(μ_{c} + 1) and ρ is the filling ratio. We limit our discussion to the impedance-matched cases, and other cases can be considered similarly. Because of the fluctuations induced by the randomness, the force acting on each cylinder is generally nonzero and shows fluctuation even for those located outside the boundary region, but the boundary stress can still be meaningfully determined if we average out the fluctuation. We calculate the stress by averaging the total force per unit length over domains of different sizes (each containing the same number of cylinders from the two sublattices) centered at the geometric interface (see fig. S1). This process is carried out for 10 sample realizations for each combination of *n*_{eff,1} and *n*_{eff,2}, and the boundary stress is taken to be an ensemble average. For each sample configuration, each sublattice consists of 300 cylinders randomly distributed in an area of *l* × *w* = 30*a* × 10*a*. The radii of the two kinds of cylinders are set to be *r*_{1} = 0.1*a* and *r*_{2} = 0.16*a*. Periodic boundary condition is applied on the upper and lower boundaries of the “supercell.” We fix *n*_{eff,1} = 2 and vary *n*_{eff,2}. The yellow circles in Fig. 5B denote the boundary stress calculated by COMSOL for the 10 realizations, and the dotted yellow line is the ensemble average. Figure 5B also shows for comparison the boundary stress of the square lattice, denoted as blue circles. A clear difference between the amorphous lattice and the square lattice is noted. To take a further step, we studied the transition of the boundary stress from a square lattice to an amorphous lattice. Let δ be the random amplitude that characterizes the deviation of a cylinder from the unit-cell center of the square lattice (δ = 0 corresponds to the unperturbed square lattice). For all cases, the distance between two arbitrary cylinders meets the condition *d* ≥ 0.8*a* so that the cylinders will not come too close together. The results are shown in Fig. 5B. Each circle denotes the boundary stress of a sample, and dotted lines are the ensemble averages, which show that the boundary stress gradually approaches that of the amorphous configuration as δ is increased. It is clear from the above results that the amorphous-lattice system gives a different boundary stress compared to that of a square-lattice system.

### EM stress at the boundary from an effective-medium perspective

The results in Figs. 2 to 5 show that specifying the standard effective-medium parameters (ε_{eff},μ_{eff}) does not provide sufficient information to determine the boundary optical stress. We iterate that there is nothing “wrong” per se in the effective-medium description. The lattice system in Fig. 1B, with λ = 100*a*, can be represented well as an effective-medium system (Fig. 1A). This can be justified by the comparison of the *E*_{z} field in the two kinds of system shown in Fig. 1C, where we considered square lattice with ε_{eff,1} = μ_{eff,1} = −ε_{eff,2} = −μ_{eff,2} = −2 (the dots in Fig. 1C mark the position of the polarizable units). The field shows very similar phase and magnitude distribution in the two systems. It can be further improved if we go to a smaller *a*/λ ratio, but the lattice structure information is still required to describe the boundary stress. If we want to stick with an effective-medium approach and avoid the burden of doing a full-wave ab initio calculation taking the lattice structure explicitly into account (either multiple scattering or with numerical solvers), we need an EM stress tensor formulation that takes into account the lattice symmetry. Such a formulation exists. The Helmholtz stress tensor (*30*–*34*) incorporates the underlying symmetry information through the electrostriction and magnetostriction effects. For a lattice system, it takes the expression (*26*, *30*)(1)where *u*_{ij} is the strain tensor and *E* and *H* are the fields in the corresponding effective-medium system. The electrostriction and magnetostriction terms (∂ε_{eff}/∂*u*_{ij}, ∂μ_{eff}/∂*u*_{ij}) are explicitly lattice symmetry–dependent. They can be analytically derived by multiple scattering theory or numerically retrieved from the band diagram of the lattice (*26*). For a 2D system and *E*-*z* polarization, one can show that(2)where φ is the angle between the direction of the local wave vector and *x* axis (for the normal incidence φ = 0). The coefficients γ and ν are symmetry-dependent parameters, and for a square lattice of *C*_{4v} symmetry, γ_{square} = 1.297, ν_{square} = − 0.596, whereas for a triangular lattice of *C*_{6v} symmetry, γ_{triangle} = 0.5, ν_{triangle} = 0 (*26*). With this stress tensor, we can analytically determine the boundary stress as(3)with γ = γ_{square} for the square-lattice system and γ = γ_{triangle} for the triangle-lattice system. Here, *H* is the magnetic field at the boundary. For a hybrid system formed by a square sublattice on the left and a triangle sublattice on the right, the boundary stress is expressed as(4)

For the triangle-square case, the positions of γ_{square} and γ_{triangle} are interchanged in the above equation. Equations 3 and 4 show that, on one hand, the boundary stress can be determined by the macroscopic field obtained using standard effective parameters, but the coefficient γ arising from the underlying lattice structure must be known in order for the boundary stress to be determinate. The solid lines in Figs. 2 and 3 mark the boundary stress calculated analytically with Eqs. 3 and 4, which agree very well with that of the microscopic model system (dots).

The Helmholtz stress tensor for the amorphous system is (*30*)(5)where the electrostriction and magnetostriction terms depend on the filling ratio ρ and . Appling the above stress tensor, one can analytically obtain the boundary stress with the same expression as Eq. 3 (with γ = 0). The result is shown in Fig. 5B (red lines) and agreement with the result of the microscopic amorphous system is again noted.

We showed in the above discussions that it is possible to reproduce the boundary stress of the microscopic system within the context of effective-medium theory. To do this, one needs to know the microscopic detail of the lattice and derive the corresponding electrostriction and magnetostriction corrections incorporated in the Helmholtz stress tensor for evaluation of the boundary stress. Expressions of Eqs. 2 to 5 for the *H*-*z* polarization can also be obtained by invoking the duality relationship and substitute μ → ε and *H* → −*E*, and an example is given in fig. S2. We note that although the Helmholtz stress tensor can predict the total boundary stress, it obviously cannot reproduce the finer details at the boundary such as the dependence of the local force density distribution on the relative shift of the lattice on the left and right as shown in Fig. 1D. There is simply no information within the context of effective-medium theory that would allow the determination of such details.

### EM stress at the boundary of a metamaterial system

We have so far considered idealized metamaterials in which the underlying microstructures are point dipoles arranged in a lattice with *a* ≪ λ. We do this on purpose so that the results we obtained cannot be attributed to the “coarse-grainedness” of the lattice or the complexity of the internal structure of the resonators, but rather to something more intrinsic. In experimentally realizable metamaterials, the internal dimension of the resonators is usually not small compared with *a*, and *a* is not that small compared with wavelength. For completeness, we also consider such configurations here by evaluating a model 2D metamaterial system composing a periodic array of core-shell cylinders (Fig. 6A) with a perfect electric conductor core of radius *r*_{inner} = 0.144*a* and a high dielectric outer shell of *r*_{outer} = 0.272*a* and ε_{shell} = 50, μ_{shell} = 1. Such a system can be used to realize an isotropic NRI material under the *H*-*z* polarization (*35*, *36*). The band diagram for a square lattice of such core-shells is shown in Fig. 6B, which shows that the second band has a negative group velocity. Using standard retrieval procedures (*37*, *38*), we obtained the effective parameters corresponding to this band, which is shown in Fig. 6C. ε_{eff} and μ_{eff} are both negative so that the system behaves as an NRI material at that frequency range. For force calculations, we use 10 layers of core-shell cylinders on the LHS to represent the NRI material and equal area of positive refractive index material on the RHS represented by a fine-grained array (lattice constant = *a*/10) of dielectric cylinders of radius *r*_{c} = 0.02*a*, as shown in Fig. 6A. Periodic condition in the *y* direction is assumed. We pick the normalized frequency ω = 0.327 at which ε_{eff} = μ_{eff} = −0.135 for the core-shell sublattice and calculate the boundary stress while varying ε_{c}, μ_{c} of the cylinders on the RHS. We note that at this frequency and with *r*_{outer} = 0.272*a*, we have a “coarse-grained” configuration and the size of core-shell cylinders is not small compared with *a*. The results of the boundary stress show that it is positive for square lattice (Fig. 6D). We then repeat the calculation for a triangular-lattice arrangement (for both LHS/RHS), with the lattice constant adjusted so that ε_{eff} = μ_{eff} = −0.135 at the same frequency. As shown in Fig. 6D, the stress now becomes negative. The dots in the figure represent the numerical results of the microscopic systems, whereas the red lines denote the analytical results predicted by using the *H*-*z* version of Eq. 3. We note agreement between the two. Figure 6E shows the *H*_{z} field distribution for the square-lattice case with *n*_{eff,2} = 1 (air), where, on the basis of the phase variation, we confirm that the Bloch wave vector does satisfy the condition of *k*_{bloch}*a* ≪ 1, which is a necessary condition for the effective-medium theory to work.

## DISCUSSION

In summary, we studied the EM stress induced at the boundary formed by two kinds of materials using a generic microscopic model. We gave explicit examples to show that square, triangle, and amorphous lattices, which are designed to have the same macroscopic ε_{eff} and μ_{eff}, give very different boundary optical stresses, both in sign and in magnitude. Hence, one cannot tell whether the optical stress at the boundary region is compressive or tensile even if macroscopic ε_{eff} and μ_{eff} are completely specified. The underlying lattice structure strongly affects the boundary stress, and one can determine whether the boundary is pushed or pulled only if the microscopic lattice symmetry is specified in addition to ε_{eff} and μ_{eff}.

We need to emphasize again here that the indeterminate value of the boundary stress is not due to the inaccuracy of the effective permittivity and permeability. In the numerical results we presented, we considered λ = 100*a*, and in this small *a*/λ regime, ε_{eff} and μ_{eff} are very accurate in the usual sense of effective-medium or homogenization theories. The lattice structure dependence will persist even if *a*/λ tends to an arbitrarily small number. Such a dependence can be understood from the perspective of the virtual work principle. Because the optical stress can, in principle, be determined by calculating the change of free energy under a virtual deformation of the lattice (*30*), the strain tensor and its underlying symmetry should come into the expression of the stress tensor, which manifests in the form of electrostriction and/or magnetostriction correction within the context of an effective-medium description. When treating such a system with effective-medium theory, the Helmholtz stress tensor can be used to predict correctly the boundary stress as the symmetry information is taken into account as electrostriction and magnetostriction terms in the Helmholtz tensor formulation. Other forms of stress tensors such as Abraham, Minkowski, or Einstein-Laub cannot predict the correct results for boundary optical stresses because of the absence of lattice structure–dependent terms (see Materials and Methods). We note that our conclusion applies to all metamaterials that have been made in practice, because all such materials are artificial composites consisting of an underlying structure.

Last but not least, we note that the macroscopic effective parameters do provide sufficient information to determine the total optical force acting on a finite-sized object, and the total force does not depend on the details of the microstructure as long as the effective parameters are accurate, which is usually the case when *a*/λ is small. It is the stress on the boundary that is indeterminate if only the macroscopic parameters are specified. If an object is immersed in an NRI medium, the total optical force is always positive (that is, light pushes the object), whereas the boundary stress can be either a pressure or a tension depending on the details of the microstructure that makes up the NRI medium. In some special configurations, a beam of light can attract an object (*6*, *39*–*44*) but that is due to the special properties of the light beam rather than the refractive index of the medium in which the beam travels. In addition to metamaterials, the current study may also find applications in the field of biophysics (for example, how the EM stress will deform a biological cell). The relationship between microscopic lattice symmetry and boundary stress in three dimensions also warrants further study in the future.

## MATERIALS AND METHODS

### Multiple scattering formulation of the microscopic model

The multiple scattering theory was formulated here for the *E*-*z* polarization, and the *H*-*z* polarization could be obtained easily by invoking duality relations. We first assumed that the semi-infinite homogeneous media on both sides of the lattice in Fig. 1B were air, in which case the system would correspond to a slab of metamaterial. We denoted the field acting on the polarizable unit locating at *i*th column and *j*th row as **E**_{ij} and **H**_{ij}, *i* = 1, 2, …, *N*_{1} + *N*_{2}, with *N*_{1}, *N*_{2} being the number of column layers of the two sublattices. The induced dipole/monopole moments were , and the following self-consistent equation can be formulated(6)
where δ_{ii′} denotes the Kronecker delta and are the position vectors of the target and source units, respectively. is the dyadic Green’s function, with and being the zeroth-order cylindrical Hankel function of the first kind. Substituting the expression of Green’s function into the above equation and using the periodic condition in the *y* direction, we obtained(7)where *n* = *i* − *i*′ and *d* denotes the transverse shift of the target layer with respect to the source layer along the *y* direction. *F*_{i}(*n*, *d*) are lattice sums that take the following expressions (*45*)(8)where *K*_{m} = 2π*m*/*a*, *K*_{x} = *q*_{0}, and , with *k*_{y} being the Bloch wave vector along *y* direction. *m* is the diffraction order, and the summation over *m* is truncated to certain value in implementation. Sign(*n*) is the sign function. Substituting the above equation into Eq. 7 makes it possible to solve the resulted linear equations for the fields at the positions of the polarizable units.

If homogeneous media instead of air were attached to the lattice on both sides as in Fig. 1B, corresponding to the case that the boundary is formed by two semi-infinite materials, both the incident and the scattered fields from the dipoles would be multireflected by the boundaries of the medium. This becomes a plane wave transmission problem because each channel of the scattered fields of the dipoles/monopoles (that is, the lattice sums in Eq. 8) takes a simple plane wave form, and the results can be found in many textbooks such as by Chew (*46*). The final form of the self-constant equations is similar to Eq. 7, except that the lattice sums now have extra contributions from the multiple reflections because of the two boundaries. Note that the boundary stress was not affected by the width of lattice as long as enough numbers of column layers were considered.

Once the field acting on the polarizable units was obtained, the induced dipole/monopole moments **p**_{ij}, **m**_{ij} could be calculated straightforwardly and the field at arbitrary position could also be easily obtained using the Green’s function. The optical force acting on each polarizable unit was then evaluated by integrating the Maxwell stress tensor over a closed loop as shown in Fig. 1B.

### Numerical simulation

For the amorphous-lattice system, we did not have a semi-analytical method that could be used to calculate the boundary stress. So, we performed the full-wave EM simulations with a commercial finite element package COMSOL Multiphysics (*29*). The positions of the cylinders were first generated by a random number generator and then imported into COMSOL to form the model geometry. Periodic boundary condition was applied in the *y* direction for the supercell, whereas scattering boundary condition was applied in the *x* direction, with a background plane wave excited from the LHS. The EM force acting on each point was then calculated by the Maxwell stress tensor method as a postprocessing step.

### Boundary stress given by Abraham, Minkowski, and Einstein-Laub tensors

In homogeneous and isotropic media, the time-averaged Abraham and Minkowski tensors use the same equation (*16*–*18*, *24*)(9)whereas the Einstein-Laub tensor can be denoted using this equation(10)For the *E*-*z* polarization of the effective-medium system, the former contributes to a boundary stress of(11)which does not have a dependence on the lattice symmetry and predicts the same result for different microscopic lattice systems. The latter simply vanishes at the boundary because of the continuity of the parallel field components. Hence, none of the Abraham tensor, Minkowski tensor, or Einstein-Laub tensor can explain the boundary stress phenomena in the main text.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/2/3/e1501485/DC1

Fig. S1. Boundary stress evaluation by Maxwell stress tensor method in the amorphous-lattice system.

Fig. S2. Boundary stress in square-lattice system under *H*-*z* polarization.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank Z. Q. Zhang and W. J. Sun for their valuable comments and suggestions.

**Funding:**This work was supported by the Hong Kong Research Grants Council (AoE/P-02/12).

**Author contributions:**S.W. did the calculations. J.N. helped on the theory of Helmholtz stress tensor. M.X. helped on the Green’s function method. C.T.C. conceived the idea and supervised the project. S.W. and C.T.C. co-wrote the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are presented in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2016, The Authors