Tracing optical force fields within graded-index media

The mechanical interaction between light and graded index media (both isotropic and anisotropic) is presented from the geometrical optics (GO) perspective. Utilizing Hamiltonian equations to determine ray trajectories combined with a description of the Lorentz force exerted on bound currents and charges, we provide a general method that we denote ‘force tracing’ for determining the direction and magnitude of the bulk and surface force density in arbitrarily anisotropic and inhomogeneous media. This technique provides the optical community with machinery which can give a good estimation of the force field distribution in different complex media, and with significantly faster computation speeds than full-wave methods allow. Comparison of force tracing against analytical solutions shows some unusual limitations of GO, which we also illustrate.


Introduction
Photon wave-particle duality, which states that photons, the quanta of light, exhibit both wavelike and particle-like characteristics, is responsible for almost all the subtleties attached to the behavior of light. Reflection, refraction and interference of light are all a consequence of wave-more illustrative than full-wave descriptions. Second, in the GO limit we need only to solve Hamilton's equations to calculate ray equations and consequently the force felt by each point of the medium under investigation. These equations are much easier to handle than a wave-based field description of the propagating and/or scattered fields, which requires construction of the energy-momentum density tensor and integrating it over the domain of interest. Third, similar to what was introduced in references [44] and [45], in the GO domain various tomographic methods can be applied to reconstruct proper electromagnetic profiles to potentially engineer mechanical interactions. As a matter of fact, the force-tracing method could potentially yield a method of forcing an optical force distribution within a medium to be aligned with our interests; in other words, it may lead to a sort of'transformation optics' for forces. Such a design methodology is likely possible only in the ray-optics limit, because the wave-optics process of constructing the energy-momentum density tensor and integrating it creates a problem much less amenable to inversion.
In the next section, we first introduce the concept of bulk and surface optical force densities and briefly explain the Abraham-Minkowski controversy. We then provide a concise review of GO and ray-tracing techniques. Then we reach the main part of the work, where a general formulation of force tracing in isotropic and anisotropic media is presented, followed by a few worked examples. is the bound magnetic current density. Equation (2) is actually a demonstration of conservation of momentum in electrodynamics that was first derived by Minkowski. If, in a dielectric with refractive index n, we write the electromagnetic energy density u in terms of photons as ω = ℏ u q V, where q is the average number of photons with angular frequency ω in a volume V , and relate it to the electric field intensity, we can immediately obtain the momentum of a single photon in the dielectric medium as np, where ω = ℏ p cand c is the velocity of light in vacuum.
However, in general the stress tensor (or the energy-momentum tensor) in Minkowski's formulation is not symmetric and it seems to be in contradiction with the conservation of angular momentum. This concern made Abraham introduce a different stress tensor and hence a different momentum density in the momentum conservation equation is the Abraham stress tensor and ⃗ f 1 is called the Abraham force density. In fact, ⃗ f 1 is the price to pay for having symmetry in the stress tensor in Abraham's formulation, and is exactly what is required to make the two formulations equivalent to each other; i.e., . Like the previous case, if we express the energy density in terms of photons and use ⃗ g A for the momentum density, we find the momentum of a photon in the dielectric medium as p n. Now we are left with a situation where two seemingly equivalent formulations give two different values for the momentum of a photon within a dielectric; this is the place where the Abraham-Minkowski dilemma appears. From a theoretical point of view, both formulations are accurate in some sense and are grounded in Maxwell's equations. From an experimental point of view, unfortunately the momentum of a photon in a medium cannot be measured directly and it must be interpreted through measurement of other entities. In experiments reported, the momenta measured were affected by the interpretation used, which again can make both interpretations seem correct [28,50]. Barnett [32,50] solved this dilemma by recalling the fact that two types of momentum can be attributed to photons: kinetic momentum as the product of the photon-associated mass and its velocity (owing to the particle nature of photons), and canonical momentum as Planck's constant divided by the photon de Broglie wavelength (owing to the wave nature of photons). The total momentum within a medium then is either the summation of the kinetic momentum of a photon and the kinetic momentum assigned to the matter or the summation of the canonical momentum of the photon and the canonical momentum assigned to the matter. According to Barnett where clearly Abraham's momentum is identified as the kinetic momentum of light and Minkowski's momentum is identified as the canonical momentum of light. In fact, all historical experiments seeming to favor either the Abraham or Minkowski momenta were missing a comprehensive inference of the photon-matter interaction within the medium under observation; this can be traced to equation (5). Last, from a macroscopic point of view, any formulation, e.g. the formulation given by Abraham or Einstein-Laub, for conservation of momentum derived directly from Maxwell's equations is equivalent to the Minkowski formulation provided it yields a correct force density (like ⃗ f 1 in Abraham's formulation) exerted within or on a medium. This general equivalence was proven in [28] and its references.

A review of GO
GO represents the asymptotic situation in optics where the wavelength is very small such that it can be approximated as λ → 0 and the energy transport can be described by the language of geometry [51]. This approximate theory is reliable only where the electromagnetic fields have rapid changes in their phases and gradual variations in their amplitudes. For the sake of accuracy in the limit of GO the media parameters (permittivity, permeability, conductivity, etc) should fluctuate slowly on the scale of at least one wavelength or, more precisely where ξ ⃡ stands for the electromagnetic constitutive tensors, σ for the electromagnetic loss or gain and ⃗ F for the fields and π λ = k 2 0 . The propagating fields can be expressed in terms of locally plane waves (quasi-plane waves), and with Fermat's principle and Lagrangian calculus they can be traced beautifully as pencils of rays [45,51].
In the domain of GO, we are allowed to write the electromagnetic fields as where ω is the angular frequency, ⃗ k is the wave vector and the magnitudes of both ⃗ E 0 and ⃗ H 0 are assumed to be approximately constant. In this domain, with the use of the corresponding in a specific problem, we arrive at Hamilton's set of differential equations, , which together form a path through space. Traditionally a ray of light is defined as a geometrical path which is optically the shortest possible path (Fermat's principle) and through which energy is transported. However, we believe that not only the energy but also the momentum of light can be tracked by rays, and this is something which has not been dealt with considerably in the rich literature of GO (though papers such as [53][54][55][56] have studied the optical force and mechanical interactions between light and matter partly from a GO point of view); this will be explored in the rest of this paper. We believe that the proposed machinery along with the equations derived in the next section can form a useful toolset in the study of mechanical interaction between light and media; this'forcetracing' involves simple equations-an a priori alternative to the usual complex calculation of optical force via full-wave analysis and integration of the stress tensor.

Force tracing
As said earlier, the total Lorentz force felt by a medium is the summation of the force acting on free sources ⃗ ( ) f s and the force acting on bound sources ⃗ ( ) f b . In the absence of free sources within a medium, ⃗ f s vanishes and ⃗ f b is the force that affects the medium. The time-averaged force density within a source-free medium therefore would be [46][47][48][49] where ⁎ stands for the complex conjugate, and ε 0 and μ 0 are the permittivity and permeability of free space, respectively. In lossless media we have 0 0 Let us consider isotropic media first. In the domain of GO we can write the electromagnetic fields as in equation (7). In isotropic media we know that ⃗ E 0 , ⃗ H 0 and ⃗ k make a right triplet [46], and if we know enough about the directions of two of them then we can calculate the direction of the third. In this sense, we can construct an arbitrary orthogonal coordinate system in which one of the basis vectors is along ⃗ k . If we call two other basis vectorŝ e z andê t , then we can decompose the electric field into its components, parallel (in-plane component) and perpendicular (out-of-plane component), to the plane defined by ⃗ k andê t . For these two polarizations, we have 0 0 0 0 However, in homogeneous isotropic media both ⃗ E 0 and ⃗ H 0 are vectors of constant amplitude and direction, so their divergences are zero and there exists a null bulk force density within such media. In homogeneous media the only thing which matters is the surface force originating from refraction and reflection of the fields at the media interfaces.
Then for the in-plane wave we have , the first term on the right-hand side of equation (15) does not have any contribution to the bulk force density, and the time-averaged force density for the in-plane polarization can be expressed with the help of equation (9) as Similarly, for the out-of-plane polarization we can obtain the following expression for the force density: where E 0 is the magnitude of the electric field of the unpolarized light. After some mathematical simplification, we obtain the following expression for the bulk force density within an isotropic medium: x y x y In equation (19), the partial derivatives can be replaced by simpler terms, since we have where α and β can take x or y. Hence we can simplify equation (19) as . It should be mentioned that in (21) we have normalized the force density by ε E 2 0 0 2 . Furthermore, the x and y components of the normalized bulk force density in a spherically symmetric isotropic medium [45] would be where b is the impact parameter of the incident ray. It should be recalled that in GO fields are locally plane waves and hence, as can be seen in (21)-(23), we are always able to decompose the force into two components along a plane of interest and the force does not have any component perpendicular to this plane. However, the orientation of the plane changes along a ray's trajectory, and thus so do the directions of the introduced basis vectors.
For the last word on the isotropic case we refer back to equation (19). With the help of (20), we can write the magnitude of the normalized bulk force density as With knowledge of the isotropic Hamiltonian and the set of ray differential equations, we are allowed to make replacements such as k x = ½ dx/dτ, k y = ½ dy/dτ, 2n ∂n/∂x = dk x /dτ = ½ d 2 x/ dτ 2 and 2n ∂n/∂y = dk y /dτ = ½ d 2 y/dτ 2 . Utilizing these replacements, we can simplify equation (24) as Interestingly, the fraction on the right-hand side of equation (25) is an important geometrical quantity-the curvature (κ) of rays inside an isotropic medium. In other words, In consequence, the bulk force density is directly related to the curvature of the optical rays and therefore we can state a rule of thumb: the greater the curvature of a light ray, the stronger the bulk force density inside an isotropic medium! Now let us consider the anisotropic case. We know that in anisotropic media ⃗ D, ⃗ B and ⃗ k make a right triplet [46]. In this part then we assign the in-plane and the out-of-plane polarizations with respect to the displacement field ⃗ D. In an anisotropic medium we have Accordingly, by inserting the electric and magnetic fields into equation (9) and having equation (10) in mind, we obtain the following expression for the bulk force density of the inplane polarization: In addition to the bulk force, there exists another force, the surface force. This force is due to discontinuities at media interfaces which reflect or refract incident light. At the interface of a medium, the force density equation (9) can be written as [46][47][48][49]  After setting up the theoretical foundation of the force-tracing process, in this part we work some examples. As a first example we consider the bulk force density distribution in an Eaton lens. An Eaton lens is a spherically symmetric graded-index lens of radius r 0 with a profile index = − ( ) n r r r ( ) 2 1 0 which acts like an isotropic reflector [45]. The ray trajectories in the Eaton lens are parts of ellipses which share a focal point at the center, and if the ambient medium is vacuum the incident rays do not feel any disturbance at the Eaton lens boundary (i.e., = n r ( ) 1 at = r r 0 ). Due to this transparency, there exists no surface force at the boundary. However, by gradual change of the index, the trajectories of the rays bend gradually and therefore an optical force proportional to the curvature of the rays would be felt inside the lens.
In order to explain what really happens to a photon traveling inside an Eaton lens, we borrow the idea of Balazs [10] about photon-dielectric interaction. In the Balazs model, a photon of frequency ω is propagating towards a transparent dielectric box of rest mass M and refractive index n. The box is assumed to be at rest on a frictionless surface. Before entering the box, the energy and the momentum of the photon are ω ℏ and ω ℏ c (the associated mass of the photon is ω ℏ c 2 ) and therefore the total energy and momentum of the system are ω = ℏ + E Mc 2 and ω = ℏ P c. (It should be noted that this P is different from what we used for the polarization vector in equation (9)). When the photon enters the box, its kinetic momentum would be ω ℏ nc and hence the box picks up a momentum of ω ℏ − ( )( ) c n 1 1 and starts to move with a velocity v, which can be calculated through conservation of energy. When the photon later exits the box, its momentum will be ω ℏ c and hence the box stops moving. We can consider a similar story for a photon's journey within an Eaton lens. Consider two adjacent points on a ray trajectory such as Q i and Q i + 1 in an Eaton lens, as shown in figure 1. At point Q i the photon has the kinetic momentum ω ℏ n c Qi , where n Qi is the value of the Eaton index profile at Q i . When the photon travels from this point to point Q i+1 , its kinetic momentum changes to ω ℏ where this momentum is proportional to the total force exerted on the lens. In the limiting case of → ∞ N , equation (35) can be written as ∫ ⃗ = ⃗ ∝ p dp optical force (36) AB C where C refers to the path of the photon. In order to find the optical force distribution inside an Eaton lens, we need to use equations (22) and (23). To do so, we can calculate the optical force density along an arbitrary ray, which can be parameterized with parameter τ. Shown in figures 2(a) and (b) are the bulk force density illustrated by black arrows along several ray trajectories and the bulk force distribution within the lens, respectively. One of the rays is shown in purple. The magnitude of the corresponding force density along this ray is depicted in figure 2(c). In addition to this, in order to show the validity of equation (26), the curvature (2κ) of the ray traced in purple is shown in figure 2(c), and as seen it completely overlaps the normalized magnitude curve.
The trajectories of light rays in isotropic media can be determined by the Euler-Lagrange equation [45] (which can be broken into the set of the Hamilton equations), and based on this equation the ray tracing parameter τ is given by dτ = dl/n, where dl is the infinitesimal increment of the geometrical path length and n is the refractive index. In order to find the total force along a ray's path we need to integrate the force density along the entire path, and according to the previous mentioned note the total force would be As seen in figures 2(a) and (b), the optical force increases at near the center of an Eaton lens where there is an index singularity. In other words, the optical force is proportional to the curvature of the rays, and in the vicinity of a singularity the optical force in this region blows up. This is depicted at the centers of figures 2(a) and (b). This can also be vividly interpreted in figure 2(c), where the curve shows a sharp peak near the middle of the ray curve and obviously the peak gets sharper and sharper as the impact parameter of the ray decreases. The other interesting thing about the Eaton lens is the symmetry of the force distribution with respect to the horizontal axis (x axis). As seen in figure 2(b), the rays with tiny impact parameters initially cause very small forces, and as they approach the singularity they face a sudden change in their trajectory and their mechanical interactions with the lens jumps up abruptly. However, this fluctuation in the force density becomes smoother as the impact parameter of the rays increases.
In fact, what determines the direction of the imparted force to the lens is the direction of the light's momentum change. Since in the graded-index medium the ray trajectory is curved, the momentum change is related to this change in the curvature of the ray trajectory. In order to validate our analysis, we performed a full-wave simulation with a finite element technique using commercial software COMSOL to numerically calculate the force density inside an Eaton lens. The simulation result is given in figure 2(d). It is seen that the force-tracing analysis is in a good agreement with the full-wave simulation; nevertheless, in the full-wave simulation, some force arrows along each ray have different heights. This is related to wave interference within the lens, and the locations of the regions of reduced force density are wavelength dependent (the wavelength simulated was slightly too large for an exact comparison to the geometric optics case-this was unfortunately necessary, however, because full-wave simulations with shorter wavelengths would require additional/significant computational resources).
As another example, we consider the Luneburg lens, in which incident parallel rays are focused at the back of the lens. The Luneburg lens is a spherically symmetric lens with index profile = − n r r ( ) 2 2 where ray trajectories are also (different) sections of ellipses [45]. Making use of equations (22) and (23), we can calculate the force distribution inside a Luneburg lens as depicted in figure 3. Like the previous case, in figures 3(a) and (b) the force density along individual ray trajectories and the bulk force distribution are depicted, respectively. In figure 3(a) one ray is chosen arbitrarily and traced in purple. The magnitude of the force density and the curvature (2κ) for the ray highlighted in purple are depicted in figure 3(c), which again confirms the validity of equation (26). The lens feels maximum pressure at initial and final points along the ray's trajectory and the minimum pressure exactly at the midpoint of the trajectory, for the curvature of the elliptic trajectory inside is smallest there. It is seen that the rays with higher impact parameters impose more force than those with lower impact parameters; the ray with zero impact parameter does not have any interaction with the lens. Again, this can be explained by the curvature of the trajectories. Shown in figure 3(d) is the full-wave calculation of force density inside the Luneburg lens in COMSOL, which is an indication of agreement between the proposed force-tracing method and the full-wave analysis.
As shown in figure 3(a), when many parallel rays enter a Luneburg lens, they propagate along trajectories to accumulate at a single point and exit at various angles corresponding to their impact parameters, and as a result the lens is pushed by the incoming rays. This is somehow vivid, for the projection of the outgoing rays' momenta along the horizontal axis is less than the corresponding momenta of the incoming rays and hence this reduction of their net momenta causes a pushing force. Conversely, if many rays were to enter a Luneburg lens from a single point, they would exit in parallel, and their momenta along one axis would increase during their trip within the lens and would cause a pulling force on the lens along that axis [57][58][59]. These two facts are illustrated schematically in figures 4(a) and (b).
In our final example we examine the optical force within an invisible cloak [45,60] designed through transformation optics; this is an anisotropic index profile. It is known that one can design an infinite number of cloaks to make an object invisible [61], but in this analysis we assume a simple case, where diag stands for a diagonal tensor, and R 2 and R 1 are the radii of the outer and inner boundaries of the cloak, respectively. It should be noted that, since ε μ ⃡ = ⃡ bulk force, in the case of this cloak a surface force also appears. Utilizing equations (31)-(34), we calculate the force distribution along the θ π = 2 plane for a cloak with dimensions = R 1 2 and = R 0.25 1 . Shown in figure 5 is the force distribution within the spherical cloak for in-plane polarized light.
As seen in figure 5(a), near the focal points of the ray trajectory (where the curvature changes its sign) we have a maximum in the force distribution, and due to the spherical symmetry there are two of these maximum force density points along each ray trajectory. In addition, as we go further from the cloaked region the rays feel less force; this makes sense, because the outer rays are not as curved as the inner rays. It is interesting to see that the surface force (the green arrows) is directed inward while the bulk force is pointing outward, and also the surface force is an order of magnitude less than the bulk force. The plot of the bulk force density is shown in figure 5(b). At points near the cloaked region the force is much larger than it is near the outer lens boundary, and this is related to the presence of a singularity at = r R 1 . It should be noted that our calculation of the force density at points very close to the singularity may not be accurate enough-at such points we are out of the domain of GO. However, as can be seen from figure 5(a), all the force density arrows cancel each other out and the total force on the perfect cloak is zero as it should be. To ensure the validity of our force-tracing calculations (as shown in figure 5(d)), we calculate the bulk and surface force densities from the analytical fullwave expressions given in [49]. It is seen that the agreement between our GO based analysis and the analytical analysis is quite acceptable. However, the force density distributions are clearly not the same. In the ray optics analysis that we have developed, the maximal bulk force density occurs when ray curvature is maximal and the maximal surface force density occurs where refraction is greatest; in the full-wave analysis these both occur along a vertical cutline through the cloak center. While the difference is not large between the two analyses, we find it to be a very interesting limitation of GO, since all of the conditions that would normally permit the use of GO have otherwise been met. Although we are unable to explain the specific reason for the differences in the two analyses in this example, we have found an even more remarkable example of this type of breakdown of GO, which we will now describe.
As is evident in figure 5(a), the force field is mirror symmetric about a vertical cutline through the cloak. When a ray enters the cloak from the left side, in its first half trip it causes a force with a negative horizontal component and in its second half trip it exerts a force with a positive horizontal component, and these two components cancel each other out. The surface force densities are at most one-tenth of the bulk force densities. Now consider cutting the cloak in half so that only the left side remains. With the right half gone, one might conclude that a half-cloak could be pulled by an incident plane wave! It would be remarkable if this were the case. However, due to the inability of GO to account for some phenomena such as diffraction, part of the truth has been hidden in such an interpretation (not to mention that it would violate conservation of momentum). In order to make it clear why such a half-cloak tractor beam is impossible even though force tracing says it is possible, we compare the ray analysis and full- wave analysis of a half-cloak in figures 6(a) and (b). In the ray analysis, a half-cloak produces a completely dark shadow, or in other words a perfect hollow beam, and the outgoing rays do not have any interference with one another. However, in reality the light exiting a half-cloak is diffracted, which makes the shadow formed behind the half-cloak incomplete. So the diffraction, not predicted by the ray analysis, is responsible for creating some pushing forces which aggregately overcome the negative horizontal components of the bulk force density. Therefore, the net force on a half-cloak by an incident plane wave is pushing. As it is beyond the scope of this paper, we will not calculate this pushing force, though it would be noteworthy if someone carried out the full computation of such a pushing force.

Conclusion
In conclusion, we have reviewed the conservation of light momentum in general and discussed the Abraham-Minkowski dilemma. With the aid of the GO toolset, we proposed a formulation for calculating the time-averaged bulk and surface Lorentz force densities within linear optical materials. Owing to the genuine simplicity of GO machinery, the formulation is much simpler than full-wave analysis and the standard process of integrating the stress tensor(s). In order to demonstrate the validity of our formulation, we considered two isotropic lenses (Eaton and Luneburg lenses) and one anisotropic device (perfect spherical cloak) and calculated their bulk and surface force densities. We hope that this presentation of light-media interactions from a simple perspective can become a useful tool in calculating the optical forces within media under the constraints of GO; in many practical situations, it gives significantly identical solutions to much more cumbersome full-wave methods.