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.

Keywords

( A ) Boundary formed by two different homogeneous media described by effective parameters. ( B ) For metamaterials, the effective medium has an underlying microscopic structure, taken to be a square lattice in this example. ( C ) The effective-medium (A) is a good representation of the microscopic system (B) if the fields calculated for them are nearly the same as shown in this panel, where we compare the calculated electric fields (E z ). ( D ) Force distribution in the model system with ε eff,1 = μ eff,1 = − ε eff,2 = −μ eff,2 = −2 and λ = 100a, with a being the size of the unit cell. d is the relative shift between two sublattices along the y direction. The inserted panel shows the total force (boundary stress) acting on the boundary, which is the sum of the nonvanishing force values in the blue shadowed region.

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 = αeE and an in-plane magnetic dipole m = αmH. 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 = (2A + αm)/(2A − α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 λ = 100a 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. Fig. 2 Boundary optical stress in impedance-matched lattice systems. (A to D) Boundary stress in the square-square (A), triangle-triangle (B), square-triangle (C), and triangle-square (D) systems. We fix the refractive index on the left n eff,1 = −1 or n eff,1 = 1 (corresponding to air) and vary the refractive index n eff,2 on the RHS. The dots correspond to the boundary stresses of the microscopic systems calculated with multiple scattering theory, whereas the lines denote the analytical values obtained using the Helmholtz stress tensor. The insets show the configurations under consideration. (E) Four types of lattice systems correspond to the same effective-medium system (n eff,1 = −1 on the left, n eff,2 = 4 on the right) but have entirely different boundary stresses. The direction of the stresses is shown in the panels by an arrow. 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. Fig. 3 Boundary optical stress in impedance-mismatched lattice systems. (A and C) Boundary stress as a function of μ eff,2 , ε eff,2 in the “air/square-lattice” configuration, respectively. (B and D) Boundary stress as a function of μ eff,2 , ε eff,2 in the “air/triangular-lattice” configuration, respectively. In (A) and (B), we fix ε eff,2 = −1, whereas in (C) and (D), we fix μ eff,2 = −1. 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. Fig. 4 Optical force distribution in finite slab and total reflecting systems. Optical force distribution in finite square-lattice (A and D) and triangle-lattice (B and E) slab systems. For (A) and (B), the impedances of the lattice systems match that of air, whereas for (D) and (E), they do not. The solid and dashed lines in (C) and (F) denote the accumulated force counting from the first layer, and the final value (as marked by a blue dot) corresponds to the total force exerting on the slab. We note that whereas the forces on the boundary depend strongly on the microscopic details, the total force acting on the slab is exactly the same independent of the microscopic details as long as the macroscopic ε eff and μ eff are the same. Optical force distribution for a semi-infinite square (G) [triangle (H)] lattice with ε eff,2 = −4, μ eff,2 = 1 attached to a semi-infinite square (triangle) lattice with ε eff,1 = μ eff,1 = −1 or air with ε eff,1 = μ eff,1 = 1. Here, x denotes the position of the unit, with a being the column layer distance. 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 = 30a × 10a. The radii of the two kinds of cylinders are set to be r 1 = 0.1a and r 2 = 0.16a. 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.8a 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. Fig. 5 Boundary optical stress in amorphous-lattice system. (A) Part of an amorphous-lattice sample drawn to scale. (B) Comparison between the boundary stress of a square lattice and an amorphous lattice. Circles denote the stress calculated numerically for the microscopic lattice systems. Solid lines denote the stress calculated analytically using the Helmholtz stress tensor in the effective-medium system. Dashed lines are the average of the values represented by circles. δ is the random amplitude defined in the main text.

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 λ = 100a, 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.