Time-dependent optical force theory for optomechanics of dispersive 3D photonic materials and devices

We present a position- and time-dependent optical force theory for optomechanics of dispersive 3D photonic materials and devices. The theory applies to media including material interfaces, waveguides, and general photonic crystal structures. The theory enables calculation of the dynamical state of the coupled field-material system and the interference of this state with other excitations of the material, such as surface acoustic waves or phonons. As an example, we present computer simulations of energy and momentum flows through a silicon crystal with anti-reflective structured interfaces. Using commercially available simulation tools, the theory can be applied to analyze optical forces in complex photonic materials and devices.


I. INTRODUCTION
The electric and magnetic fields of light in arbitrary 3D photonic materials and devices can be unambiguously solved from Maxwell's equations. Considering the interest to develop optomechanic materials and devices, it is astonishing that no generally applicable time-dependent optical force theory has been presented for inhomogeneous dispersive materials [1,2]. The optical force density is ultimately related to the conservation law of momentum. Optical forces having the same physical origin as those discussed in the present work have already led to important scientific discoveries and photonic technologies in optical trapping and laser cooling. Optical forces are also related to the centenary Abraham-Minkowski controversy of the difference of the momenta of light in a material and in vacuum [3][4][5][6][7]. In previous literature, most works assume the Maxwell or Helmholtz stress tensor and further take the time average of the resulting optical force density [1,2,[8][9][10][11][12]. Very few works, if any, are studying the exact position-and time-dependent optical force density in inhomogeneous dispersive materials. However, the knowledge of the exact position-and timedependence of the optical force density is of fundamental importance for quantitative analysis of the interaction of light with elastic waves or other mechanical eigenmodes of the material.
The recently introduced mass-polariton (MP) theory of light [13][14][15][16][17][18] considers a light pulse propagating in a material as a coupled state of the electromagnetic (EM) field and the material and splits the total momentum of light into the field and material components. The material component of the momentum is carried by an atomic mass density wave (MDW) driven forward by optical forces. In this work, we generalize the MP theory for structured interfaces and give a general expression for the optical force density in an inhomogeneous dispersive dielectric. The optical force density is derived from the continuity equations of the energy and momentum at interfaces and from the covariance principle of the special theory of relativity. The optical force density is tested by computer simulations of energy and momentum flows through a silicon crystal with anti-reflective structured interfaces [19].
The exact time-dependent force density enables accurate calculation of the time-and position-dependent dynamics of the material at structured interfaces and in the bulk, which can also be a photonic crystal. This provides an interesting possibility to discover eventual coupling of the field-driven material dynamics in photonic materials and devices with acoustic waves and a possibility to fine tune this coupling for improved or new operational characteristics of photonic devices. This paper is organized as follows. In Sec. II, we describe the physical foundations for the separation of the energy density, momentum density, and the stress tensor of the field-material system unambiguously into the field and material parts. In Sec. III, the differential forms of the conservation laws of energy and momentum are briefly presented for the field and the material. In Sec. IV, we solve for the quantities appearing in the conservation laws in the special case of the laboratory frame, i.e., in the rest frame of the material. In Sec. V, we present an example simulation of a Gaussian plane wave light pulse propagating through a silicon crystal block with anti-reflective structured interfaces. In Sec. VI, we present the relativistically invariant expressions of the stress-energy-momentum (SEM) tensors of the field and material parts of the system. Finally, conclusions are drawn in Sec. VII.

II. SEPARATION OF THE SYSTEM INTO FIELD AND MATERIAL PARTS
A. Separation of the energy density, momentum density, and the stress tensor To enable unambiguous separation of the energy density, momentum density, and the stress tensor of the field-arXiv:2112.11128v4 [physics.optics] 20 Jul 2022 material system into the field and material parts, we define the energy density of the material to include only the classical rest energy density and kinetic energy density of atoms. All other forms of the energy density, such as the energy density of the polarization field, are considered to be field energy density. It is well known that under the influence of classical EM field, the atomic trajectories follow very accurately classical Newtonian dynamics, and therefore, their rest energy and kinetic energy are strictly unambiguous and equal to their classical values. Therefore, within our definition of the energy density of the material, there cannot be any uncertainties in the splitting of the total energy density or momentum density into the field and material parts. This means that the energy and momentum densities and the stress tensor of the material in a general inertial frame are uniquely given by the classical formulas W mat = ρ a c 2 , G mat = ρ a v a , T mat = ρ a v a ⊗ v a , respectively, where ρ a is the mass density, v a is the atomic velocity, c is the speed of light in vacuum, and ⊗ denotes the outer product of vectors [20,21]. The mass density is given by ρ a = γ va m 0 n a , where m 0 is the rest mass of an atom, n a is the number density of atoms, and γ va = 1/ 1 − |v a | 2 /c 2 is the Lorentz factor. The relations above preserve Newton's equation of motion of the material, in all inertial frames, in the form n a d dt p a = f opt , where p a = γ va m 0 v a is the momentum of an atom and f opt is the optical force density [13,22]. The number density of atoms satisfies the continuity equation ∂ ∂t n a + ∇ · (n a v a ) = 0 [13,22]. Note that the elastic force density and the electro-and magnetostrictive force densities between atoms, discussed in Sec. II B, can, however, modify the energy and momentum densities and the stress tensor of the material above as briefly discussed in Sec. VI. The energy density, momentum density, and the stress tensor of the field in a general inertial frame are denoted by W EM , G EM , and T EM , respectively.
Note that all quantities in the general inertial frame are unambiguously determined by their values in the laboratory frame by their respective well-known Lorentz transformations. The Lorentz transformations needed are collected in Table V of Ref. [13], which also presents a complete record of the covariant theory of light in dispersive materials. In spite of using a different expression of the optical force and the SEM tensors, in this work, the covariance properties presented in Ref. [13] are as is strictly applicable and fulfilled for the theory presented. Since, in our theory, the energy density, momentum density, and the stress tensor of the material are uniquely determined by the classical formulas as explained above, the separations of the total energy density, momentum density, and the stress tensor of the field-material system into the field and material parts are also unique. In previous literature, some works claim that the separation of the momentum density of light into the field and material parts would be arbitrary [23], while some other works base their separation of the system into the field and material parts on arguments that are different from those of us [24,25].
The previous challenges in the unambiguous splitting of the system into the field and material parts may be related to adding to the energy density of the material, for instance, the energy density of the polarization field, which has its origin in the field energy.
B. Separation of the total field-induced force density into gradient forces and forces that carry the wave momentum of light In addition to the optical force density f opt acting between the propagating EM field and the material atoms along the wave vector of light, there are optically induced intra-material force densities, called electro-and magnetostriction, which appear in the form of a gradient force f st acting along optical energy density gradients [26][27][28][29]. Very recently, the electrostrictive force density has also been verified experimentally for an optical field [30]. The splitting of the total force density into gradient forces and forces associated with the wave propagation has also been discussed in a recent preprint [31]. The optical electroand magnetostrictive force densities f st rise from the dependence of the optical electric and magnetic energy densities on the atomic density through the permittivity and permeability of the material [26]. The integral of the gradient force densities over the volume of the material is zero at any time, in particular when a light pulse is entering the material. Thus, these force densities do not lead to the movement of the center of mass of the material even though they give rise to local accelerations of the material atoms. Therefore, f st differs from f opt , which is a topic of the present work. In the present work, we have preserved the term optical force to mean only forces that carry the wave momentum of light. Consequently, we focus on the investigation of the dynamical effects of f opt and leave the detailed time-dependent study of f st in dispersive materials as a topic of a future work.

III. CONSERVATION LAWS OF ENERGY AND MOMENTUM
We solve the optical force density in an inhomogeneous dispersive material starting from the conservation laws of energy and momentum for the EM field, written for a general inertial frame in a differential form as [13,22,26,32,33] 1 Here φ opt = f opt · v a is the optical power-conversion density of a lossless material. In the rest frame of lossless materials, i.e., the laboratory frame, the optical power-conversion density can be approximated to zero as φ (L) opt = 0, since it only converts an exceedingly small amount of EM energy to the kinetic energy of atoms [13]. The conservation laws for the energy and momentum of the material are, accordingly, given for a general inertial frame in a differential form by [13] 1 Here it is essential to note that the source terms on the right hand side are opposite to those of the EM field in the conservation laws in Eqs. (1) and (2). Thus, the total quantities of the field-material system, i.e., sums of the field and material parts, satisfy the conservation laws in the form, where the source terms are zero. When the conservation laws in Eqs. (1)-(4) are satisfied, it is automatically guaranteed that the energy and momentum are conserved at material interfaces and that the law of constant center of energy velocity is satisfied for an isolated system.

A. Energy and momentum densities of the field
The EM contribution to the total momentum density of light is continuous over material interfaces and given in the laboratory frame by [26,32] Here E (L) is the electric field and H (L) is the magnetic field solved in the laboratory frame from Maxwell's equations with position-and frequency-dependent permittivity and permeability. Using the EM momentum density in Eq. (5), the exact time-dependent energy density of the EM field can be solved from the conservation law of energy in Eq. (1) with φ (L) opt = 0. For a narrow frequency band in a locally isotropic lossless dispersive material, where the permittivity and permeability are real-valued scalars, the EM energy density becomes [13] W (L) Here ω 0 is the central angular frequency of the narrow frequency band, and ε = ε(ω 0 ) and µ = µ(ω 0 ) are the permittivity and permeability of the material in the laboratory frame at ω 0 . For a narrow frequency band, the permittivity can be approximated by the first two terms of its Taylor series ε(ω) ≈ ε(ω 0 ) + ∂ε(ω0) ∂ω0 (ω − ω 0 ). A similar expansion is valid for the permeability. The angle brackets in Eq. (6) . The brackets here and below are also equal to the time average of the pertinent quantities over the harmonic cycle, e.g., /ω 0 is the length of the harmonic cycle. The EM energy density in Eq. (6) is equal to that in Ref. [13], but it is expressed here in terms of the permittivity and permeability of the material and their derivatives. In previous literature, there are some works that add to the EM momentum density of the field in Eq. (5) additional dispersion related terms [34,35]. In contrast, in the present work, we rely on the EM momentum density in Eq. (5) and show that it leads to a consistent theory. In our theory, dispersion dependence appears in the momentum density of the EM field only in the general inertial frame as presented in Sec. VI.

B. Optical force density and stress tensor of the field
Having the energy density of the EM field solved from the conservation law of energy in Eq. (1), we are left with a problem of determining the EM stress tensor and the optical force density from the conservation law of momentum, given in Eq. (2). Appriori, if we are only given the conservation law of momentum in Eq. (2) with the EM momentum density in Eq. (5), both the EM stress tensor and the optical force density are unknown. Therefore, more information is needed to determine these quantities and to guarantee their uniqueness. In Ref. [13], we postulated for the optical force density, in the rest frame of a homogeneous material, a generalized expression of the Abraham force, given by f opt,A to determine the corresponding EM stress tensor [13]. In Ref. [13], it was proven that the force density f (L) opt,A leads to a consistent covariant theory. However, it was also pointed out that the covariance condition and the continuity equation are not sufficient to quarantee the uniqueness of the optical force density. In this work, we use a different approach by introducing an additional constraint to the available EM energy density and the total momentum density of light. This constraint is shown to lead to unambiguous expressions for the optical force density and the EM stress tensor.
To obtain the necessary constraint, we argue that, for a light pulse propagating in a homogeneous material, for every differential volume element, the total momentum density of light must have the same constant proportionality to the EM energy density. Without this local property, if an arbitrarily small amount of energy is absorbed from the EM field, the absorbed momentum would depend on the space-time point, where the absorption takes place. This would imply that the energy-momentum ratio of the non-absorbed part of the light pulse also changes, which is not physically meaningful. The constant of proportionality is required to be n p /c, which agrees with the high-precision experimental results for the radiation pressure in dispersive liquids by Jones et al. [36]. Thus, the total momentum of light G wherê v g is a unit vector in the direction of the group velocity of light. Using Newton's equation of motion, we can write G opt dt , which relates the given condition unambiguously to the expression of the optical force density. After some algebra, the optical force density is given by Here the angle bracket expression is given by . This time-and position-dependent quantity is also equal to the time average of the Poynting vector over the harmonic cycle. In Eq. (7), we have added the last two well-known gradient terms [26] to conform with the conservation law of momentum in Eq. (2) at material interfaces. The EM stress tensor corresponding to the optical force density in Eq. (7) is found to be given by Here I is the 3 × 3 unit matrix. Equation (8) is the symmetrized form of the conventional Abraham stress tensor [13,23,37]. The symmetrization is essential in birefringent materials, where the electric and magnetic flux densities D (L) and B (L) are not parallel to the electric and magnetic fields E (L) and H (L) , respectively [26]. The optical force density in Eq. (7) and the EM stress tensor in Eq. (8) lead to a theory that fulfills the same covariance conditions as the theory presented in Ref. [13]. Accordingly, the integrated value of the momentum density of a homogeneous material over the volume of the light pulse following from f (L) opt is equal to that obtained using f (L) opt,A .

C. Atomic mass density wave
We next consider the implications of the optical force density to the dynamical state of the material. We define the energy and momentum densities of the atomic MDW as differences between the energy densities and momentum densities of the material under the influence of the EM field and in the absence of it as W MDW = W mat − W mat,0 = ρ a c 2 − ρ a0 c 2 and G MDW = G mat − G mat,0 = ρ a v a − ρ a0 v a0 . Here ρ a0 is the atomic mass density in the absence of light and v a0 is correspondingly the atomic velocity in the absence of light, i.e., v (L) a0 = 0 in the laboratory frame. The momentum density of the material atoms, G opt dt , resulting from the optical force density, is obtained as a solution of Newton's equation of motion. Using the expression of the optical force density in Eq. (7) then gives the momentum density of the atomic MDW in the laboratory frame as The excess energy density of the material, W Other dynamical quantities of the material, such as the strain, are also obtained from the solution of Newton's equation of motion.

V. SIMULATION OF FORCE AND MOMENTUM DENSITIES
To investigate the physical implications of the unified optical force theory presented above, we have simulated light in various geometries including several material interfaces. All of these studies support the expression of the optical force density in Eq. (7). As an example, we present numerical simulations of the propagation of a Gaussian plane wave light pulse through a silicon crystal block with anti-reflective structured interfaces illustrated in Fig. 1. Minimizing reflection with structured interfaces has been discussed more extensively in Ref. [19].
In the simulations, we assume an incident Gaussian plane wave light pulse with a central vacuum wavelength λ 0 = 1310 nm. For this wavelength, the phase refractive index profile of silicon, based on a linear fit to the experimental data from Ref. [38], is n p = 3.5039 and the group refractive index is n g = 3.6840. The electric field of the Gaussian plane wave light pulse propagating in vacuum in the direction of the positive x axis and polarized along the y axis is given by Here k 0 = ω 0 /c is the wave number in vacuum, ∆k 0 is the standard deviation of the wave number in vacuum, and E 0 is the normalization factor corresponding to the We have tested our unified optical force theory by simulating the fulfillment of the energy and momentum conservation laws when a light pulse propagates through a silicon crystal block with anti-reflective structured interfaces.
peak value of the electric field of the incident Gaussian light pulse. In the simulations, we use ∆k 0 = 0.025k 0 and E 0 = 30 kV/m. The electric flux density, magnetic field, and magnetic flux density follow unambiguously from the electric field in Eq. (11) through Maxwell's equations. The simulations are carried out using Comsol Multiphysics simulation tool [39]. Periodic boundary conditions are applied in the direction of the y axis in Fig. 1. The periodicity of λ 0 /4 is determined by the periodicity of the structured vacuum-silicon interface. Figure 2 presents the optical volume and interface force densities due to the Gaussian light pulse (see Visualization 1). The snapshot is taken at the instance of time when the center of the light pulse is at the position of the entry interface of the silicon crystal. The volume force density in Fig. 2(a) and the interface force density in Fig. 2(b) have both x and y components that vary as a function of the position. The force density maxima are around positions where the field is focused by the local lensing effect of the vacuum-silicon interface. However, after the interface, the focusing fades out and the force density follows the plane wave form of the field. Figures 3(a)-(d) show the position dependencies of the EM momentum density, MDW momentum density, interface momentum density, and the MDW mass density, respectively, at the instance of time when the center of the light pulse is at the position of the entry interface of the silicon crystal corresponding to the force densities in Fig. 2 (see Visualization 2). The EM momentum density in Fig. 3(a) consists of the incident and transmitted contributions. The reflection is close to zero due to the anti-reflection property of the structured interface. The atomic MDW momentum density in Fig. 3(b) reminds the EM momentum density inside silicon on the right and there is no MDW in vacuum on the left. The interface momentum density at all positions in Fig. 3(c) is directed toward vacuum. Figure 3(d) presents the mass density of the atomic MDW, which is driven forward by the optical volume force density. This transfer of mass and the related rest energy is crucial for the fulfillment of the covariance principle of the special theory of relativity [14,15].
The propagation of the Gaussian light pulse through the silicon crystal is illustrated in Fig. 4, where the EM momentum density is plotted at three instances of time: first, before the crystal (t = 197 fs), second, inside the crystal (t = 393 fs), and third, after the crystal (t = 590 fs). The very small reflections from the first and the second interface are not visible because of their smallness. The EM energy reflected from the first interface is 0.14% of the total incident EM energy, and the EM energy reflected from the second interface is also 0.14% of the total incident EM energy. Thus, the total EM energy lost from the incident pulse in the full transmission process is 0.28%. One can conclude that the antireflection property of the interfaces is relatively well optimized. Inside the crystal the pulse in Fig. 4 is spatially much narrower than in vacuum. This is due to the lower group velocity of light in the crystal compared to that in vacuum. We can also see that the largest magnitudes of the EM momentum density in the crystal and in vacuum are equal. This follows from the continuity of the Poynting vector.
The total EM momentum of the pulse is the volume integral of the EM momentum density. Neglecting the small reflection, the ratio of these integrals in the crystal and in vacuum is equal to the ratio of the corresponding group velocities. Thus, the EM momentum in the crystal is p EM = p 0 /n g , where p 0 is the momentum of the incident EM field. The momentum of the atomic MDW is equal to p MDW = (n p − 1/n g )p 0 . Thus, the total momentum of the coupled MP state of light in the crystal becomes p MP = n p p 0 . This phase index proportionality of the total momentum of light is in agreement with the high-precision measurements of radiation pressure in liquids by Jones et al. [36] and the measurements of the recoil momentum of atoms in a dilute Bose-Einstein condensate gas by Campbell et al. [40].

VI. SEM TENSORS AND THEIR RELATIVISTIC INVARIANCE
The SEM tensor of a physical system collects the energy and momentum densities and the stress tensor in a single physical quantity. The general contravariant form of an arbitrary SEM tensor in the Minkowski space-time is defined by T = T αβ e α ⊗ e β , where the Einstein summation convention is used, and e α and e β are unit vectors of the space-time. The Greek indices range over the four components of the space-time, i.e., (ct, x, y, z). The corresponding matrix representation of T is given by [21,32,41] Here the superscript T denotes the transpose. Since the SEM tensor of the material, made from the quantities discussed in Sec. II, is symmetric by definition, so is also the EM SEM tensor since the total SEM tensor of the system must be symmetric. Therefore, all of our SEM tensors are symmetric and strictly based on the classical definition in Eq. (12). The SEM tensor T mat of the material can be formed by substituting the energy density, momentum density, and the stress tensor of the material from Sec. II into Eq. (12). Even more compactly, T mat can be written using the four-velocity of the material, given by U a = γ va (c, v a ). The resulting compact form of T mat is given by [20,21] Consequently, the SEM tensor of the field is also unambiguously determined in all inertial frames. Note that accounting for the elastic force density and the electroand magnetostriction will modify the SEM tensor of the material in Eq. (13). The resulting modified SEM tensor of the material is expected to be of the form presented in Ref. [42]. Its detailed dynamical study under the influence of the optical field is left as a topic of a future work.
Following the approach presented in Sec. IVB of Ref. [13] to use the Lorentz transformation to determine the SEM tensor of the field in an arbitrary inertial frame once the SEM tensor of the field in the laboratory frame is known, in the present case, we obtain Here U a0 = γ va0 (c, v a0 ) is the four-velocity of the rest frame of the material, Tr(x) is the trace of a matrix, and g = g αβ e α ⊗ e β , with g 00 = 1, g 11 = g 22 = g 33 = −1, is the diagonal matrix representation of the Minkowski metric tensor. The EM field tensor F and the EM displacement tensor D in Eq. (14) are given in contravariant forms F = F αβ e α ⊗ e β and D = D αβ e α ⊗ e β as The term ρ EM,disp appearing in Eq. (14) is an effective EM mass density term associated with dispersion. In the special case of the laboratory frame, it is defined through the last two terms of the EM energy density in Eq. This quantity transforms between inertial frames as a mass density and it is very convenient in expressing the relativistically invariant SEM tensor of the EM field in Eq. (14). The SEM tensor T mat,0 of the material in the absence of light is given by a similar form as the instantaneous SEM tensor T mat of the material under the influence of the EM field in Eq. (13). Thus, T mat,0 is given by The SEM tensor of the atomic MDW is defined as the difference T MDW = T mat − T mat,0 . Substituting the expressions of the tensors T mat and T mat,0 into this difference then gives Being given in terms of the EM field and displacement tensors, four-velocities, and Lorentz scalar factors, all SEM tensors presented in this section are transparently covariant, i.e., they have the same expressions in all inertial frames, thus, satisfying the conditions of relativistic invariance.
Note that the sum of the last three terms of T EM in Eq. (14) can be viewed as a relativistically consistent generalization of the Abraham SEM tensor, where the Abraham SEM tensor is allowed to depend on U a0 . This has been discussed in previous literature, see, e.g., Eq. (42) of Ref. [43], Eq. (2.10) of Ref. [44], and Eq. (41) of Ref. [13]. Consequently, for this generalization of the Abraham SEM tensor, the previous claim [24,37,[45][46][47] that the Abraham SEM tensor would be relativistically invalid is no longer true.

VII. CONCLUSIONS
In conclusion, we have presented an unambiguous position-and time-dependent optical force theory applicable to simulations of the propagation of light pulses in inhomogeneous dispersive materials. With existing simulation tools, the theory enables detailed modeling of the optomechanics of 3D photonic materials and devices. For example, the present theory has been applied to negative-index metamaterials in Ref. [48]. For absorbing materials, the well-known Lorentz force density should be added to the force density of the present work to describe also forces related to absorption. In the nonlinear optics regime, other phenomena, such as the Kerr effect, contribute and, accordingly, one must separately consider effects like electrostriction. The present classical field theory must be extended to the quantum domain to describe also the local torque originating from the interaction of the spin of light with the atoms.