Lineshape study of optical force spectra on resonant structures

: Understanding the frequency spectrum of the optical force is important for controlling and manipulating micro-and nano-scale objects using light. Spectral resonances of these objects can significantly influence the optical force spectrum. In this paper, we develop a theoretical formalism based on the temporal coupled-mode theory that analytically describes the lineshapes of force spectra and their dependencies on resonant scatterers for arbitrary incident wavefronts. We obtain closed-form formulae and discuss the conditions for achieving symmetric as well as asymmetric lineshapes, pertaining, respectively, to a Lorentzian and Fano resonance. The relevance of formalism as a design tool is exemplified for a conceptual scheme of the size-sorting mechanism of small particles, which plays a role in biomedical diagnosis.


Introduction
Light field can exert an optical force on an object.Optical force is widely used in atomic physics [1] for manipulating atoms.In the area of nanoscience, optical force can be used for trapping and positioning of nanoparticles [2,3], with important implications in biomedical applications especially in drug delivery [4] and cell sorting [5].
The frequency dependence (spectrum) of the optical force exerted on a particle is determined by the presence of resonances that the particle can support.Near a resonance frequency, the optical force can be strongly enhanced [6][7][8][9][10], and the force spectrum may exhibit an asymmetric lineshape [11].In the case where the sign of the force varies from positive to negative as the frequency is swept across the resonance, one can choose the appropriate frequency of the incident light to repulse, attract, or trap the particle [12][13][14][15].Therefore, a theoretical understanding of the spectrum of the optical force is essential.
There has been considerable interest in the theoretical study of optical forces on small particles.For particles with dimensions much smaller than the wavelength of the incident light, the electromagnetic response of the particle can be treated as that of a point dipole [16], and a simple formula for the force then naturally follows.For particles of arbitrary size and shapes, the optical force spectra need to be computed by numerical methods that solve the Maxwell equations, such as the finite-difference time-domain [17] or the finite-element methods [18].For particles that are either cylindrically or spherically symmetric, the force spectra can be computed analytically or semi-analytically using the Mie theory [19].However, none of these methods directly yield an insightful closed-form formula for the optical forces that applies to general structures.
In this paper, we provide an analytic study that captures the general characteristics of optical forces that emerge due to the presence of resonance.For our study, we use the temporal coupled-mode theory [20,21], which has been successfully employed to study resonant behaviors of photonic structures.In particular, this formalism was applied to study the spectra of absorption and scattering cross-sections for nanoparticles [22][23][24].Here we extend this formalism and apply it to study the optical force spectrum.We show that the temporal coupled-mode theory enables us to derive simple closed-form formulae that describe general aspects of the optical force and expose unique features of the optical force spectrum.
Our work here complements existing theoretical and computational works on the optical force.Unlike direct numerical simulations [10,25], where the starting point is the permittivity distribution of a given optical structure, in temporal coupled-mode theory the starting point is prior knowledge of the existence of optical resonances that dominate the optical response [20,21].The temporal coupled-mode theory then links optical properties of the structure, such as the force responses and cross-section, to the properties of these resonances such as their resonant frequencies and linewidths.For the computation of the optical force on a specific single structure, in many cases there is no advantage of using temporal coupled-mode theory, since the frequencies and the linewidths of the resonances would still need to be determined through numerical simulations, and the computational cost of such simulations often is comparable to what is required for determining the spectra in the first place [26].The importance of the temporal coupled-mode theory, rather, is that it provides a way to develop an understanding of the general behavior of the optical structures, in terms of a few parameters that have direct physical meanings [21].As such the temporal coupled-mode theory not only provides a better insight into the governing mechanisms but also has high merit as a design tool that has been very widely used phenomenologically to guide the design of complex optical structures [26][27][28][29][30][31].For example, one may use coupled-mode theory as a guide to obtaining a specific line shape of the force by controlling the positions and the linewidth of the resonances.Moreover, for a given structure the coupled-mode theory can be used to efficiently design an optimal wavefront for specific applications.We foresee that our efforts here in adopting temporal coupled-mode theory for the treatment of optical force will find similar use in the design of optical structures for the optomechanical control.
Our paper is structured as follows: In Section 2, we first derive the temporal coupled-mode equations governing the optical force on a general resonant scatterer under arbitrary coherent wavefronts.In Section 3, we apply the force formalism for elongated small particles.We exemplify our formalism for an azimuthally symmetric cylindrical scatterer under one (Section 3.1) and two (Section 3.2) incident plane waves, in which we demonstrated a case where the scattering cross-section deviates from the force spectrum.We finally show a conceptual application of sorting elongated particles in Section 4. Our findings have potential opportunities in precise optical switches, sensors, and tweezers.

Force operator and the coupled-mode theory framework
The numerical calculations of the force for specific cases do not necessarily provide insights into the general underlying physics.In contrast, we develop here a model that elucidates the general characteristics of the force spectrum and their relations to optical resonances of the particle.The results do not depend on the origins of the resonance (structural, material) and thus are broadly applicable to different physical systems supporting resonances.
The scheme for our formulation is as follows.The total electromagnetic field in the embedded medium of the particle is decomposed into a sum of input and output channels (channel can be any degree of freedom or eigenmode of the field).Using this representation, the force acting on the particle is extracted from the fields by the Maxwell's stress tensor and converted to a "force operator" acting between the aforementioned input and output channels.Then, to extract the amplitude of the output channels, the temporal coupled-mode theory is employed which solves the engagement between the input to output channels as mediated by the resonant internal modes of the scatterer.Although our paper is theoretical, the connection between main theory variables to experimental data is straightforward.The coupled-mode parameters such as resonant frequency, decay rates, coupling coefficients, can be obtained by fitting the experimental scattering/absorption cross-section spectrum.The complex amplitude of the output channels can be extracted experimentally by interfering with the output wavefront with a plane-wave reference on the face of a photodetector array [32].These establish a direct connection between our theoretical formalism and experimental measurements.
As shown in Fig. 1 (a), a scatterer is placed near the origin and is surrounded by air.In the presence of a monochromatic incident field at angular frequency ω, the total phasor of the magnetic field ⃗ H outside the scatterer can be decomposed as where H 0 is a normalization factor, and |c in,n | 2 and |c out,n | 2 represent the power density carried by the incoming and outgoing waves in the n-th channel, respectively.The n-th incoming (outgoing) channel is described by the wave basis functions ⃗ V in(out),n (⃗ r), and ⃗ r is the position vector.As an example, the explicit forms of the incoming and outgoing waves in the cylindrical basis are given by Eqs. ( 32) and ( 33) in Section 6.Throughout the paper, n is an integer indexing the channel and ranges from −∞ to +∞.The actual number of channels that are involved in the scattering process is related to the ratio between the size of the scatterer and the wavelength, and the permittivity of the scatterer.T over any surrounding surface A that encloses the scatterer, where n is the local normal vector of the surface.With the decomposition in Eq. ( 1), the optical force density can be expressed by a force operator ⃗ P in a compact form (see derivation in Appendix), where c in,out is the vector of input/output channel amplitude.Equation (3) can be used to solve for the forces exerted by an arbitrary input light field on any scatterer in either 2D or 3D geometries.
For 3D structures, one can use for example the vector spherical harmonic wave basis [33], and the resulting force operator can be derived from the 3D Maxwell's stress tensor.Typically, for an input wave described by c in , the output wave c out can be solved with a suitable solver as described in the introduction.Here, to extract the general aspects of the forces exerted on a resonant scatterer with geometries of different dimensions, we employ the temporal coupled-mode theory formalism [21], which applies to both 2D and 3D systems, to connect c in and c out .
In the framework of the temporal coupled-mode theory formalism, we start from the knowledge that the scatterer supports M resonant modes.The interaction of the incident light with these resonant modes is assumed to dominate the optical response of the scatterer in the frequency range of interest around the resonant frequencies.The amplitude of the resonant modes is represented as an M-dimensional vector a, where |a i | 2 represents the energy density stored in the i-th resonant mode.We note that each resonant mode can couple to several different incoming/outgoing channels.Then, for this general system, the temporal coupled-mode theory equations can be written as: is the resonant frequency of the i-th resonant mode, and i is the mode index.Throughout the paper we use the e jωt convention.In the subsequent part of the paper, when we consider scatterers with mirror symmetry, we will use negative (−|i|) and positive (|i|) indices to label the even and odd modes, respectively.Γ describes the radiative decay rate of these resonant modes due to their coupling to the outgoing channels [34].
) is half of the intrinsic energy loss rate of the i-th resonant mode due to, for example, material loss.The background scattering matrix B describes the scattering process in the absence of the resonant scatterer.This process involves a direct scattering from the incoming to outgoing channels without coupling to the particle resonant modes.The matrix elements of B vary slowly as a function of frequency across the resonant linewidth.The optical force exerted on the scatterer can be evaluated by its signature on the momentum transfer between the N incoming and outgoing channels.These incoming channels are assumed to couple significantly to M resonant modes, described by the N-by-M input coupling matrix K. Similarly, the resonances couple with N outgoing channels via the N-by-M output coupling matrix D, where D is equal to K by time-reversal symmetry [21].Using these coupled-mode theory parameters, c out can be connected to c in via the scattering matrix S, where T(ω) = D[jωI − jΩ + Γ + Γ 0 ] −1 D T describes the resonant scattering process.We also introduce a matrix B 0 which is the scattering matrix of the system in the absence of the scatterer [23].For the general case where the scattering matrix S deviates from B 0 (detailed in Eq. ( 39)), we can define the scattered wavefront vector by subtracting the directly transmitted part B 0 c in from the outgoing wavefront vector: The coupled-mode theory formalism introduced in this section is generally applicable for the optical force spectrum of any scatterer and arbitrary incident wavefront and is independent of the physical source of the resonance.However, by the nature of the temporal coupled-mode theory, this theoretical study of optical force lineshape is applicable only when the total loss rates, represented by Γ 0 + Γ, are significantly smaller than the resonant frequencies.
By using the modal decomposition of Eq. ( 7) in Eq. ( 3), the optical force can be computed as, The first term in the right-hand side of Eq. ( 8) vanishes by P + B 0 P T B 0 = 0. We note that our decomposition of incoming and scattered field contribution is consistent with the result in [13].Substituting the coupled-mode theory result for c scat in Eq. ( 7), the optical force can be divided into three terms, ⃗ Equations ( 9a)-(9c) describe three distinct contributions to the optical force.⃗ F bg describes the force due to background scattering, which is non-resonant and vanishes in the absence of background scattering when B = B 0 for a scatterer size that is much smaller than the wavelength.It has only a weak frequency dependence and forms a smooth background in the force spectrum.⃗ F res is attributed to the resonances and exhibits a Lorentzian lineshape around each resonant frequency ω i .⃗ F inter results from the interference between the background scattering and the resonance wavefront.Since ⃗ F inter includes spectral components that are odd functions of the frequency difference (ωI − Ω) shown in the matrix T(ω), it exhibits an asymmetric lineshape around each resonant frequency.Equation ( 9) is particularly useful for analyzing the force spectrum for a resonant scatterer.

Optical forces on elongated particles
To show the use of our formalism in extracting the physics behind the optical forces, we derive analytical results for the case of elongated particles that can be modeled as infinitely long cylinders with an arbitrary cross-section.This scenario requires simpler mathematics compared to the full-fledged 3D case and thus more directly illustrates the main general aspects of our formalism.The scenario is depicted in Fig. 1 (b).The infinite cylinder axis is along z-axis.We deal with optical forces due to incident waves propagating in the xy plane with magnetic field vector along the z-axis ("Transverse magnetic (TM) waves").A convenient basis to describe both the channels and the internal modes of the scatterer is the cylindrical harmonic basis (see Eqs. ( 32) and (33).In this configuration, the force and scattering/extinction cross-sections are defined as per unit z-length.Therefore, we calculate ⃗ F by integration of the Maxwell's stress tensor over a circle in the xy plane.The matrix element of the force operator of Eq. ( 3) is indexed according to the channels m and n, where c is the speed of light in vacuum, and δ m,n is the respective Kronecker delta.The detailed derivation is in the Section 6.1 in the Appendix.This form of the force operator in Eq. (10) where the contributions to the force are nonzero only between channels with indices different by ±1, is consistent with the result in [35].As shown in Fig. 1 (b) where a general incident wavefront exert a force on a two-dimensional scatterer, from Eq. ( 3), the optical force is given by, where we obtain ⃗ P mn from Eq. ( 10) and combine it with the incoming and outgoing wavefront vectors c in,n and c out,n .To apply Eq. ( 11) for an arbitrary incident wavefront, one can obtain c in by decomposing the incoming wave as r).Before introducing the numerical examples below, we briefly comment on the resonance mode in the system.For simplicity, we assume that these resonant modes form an orthogonal basis.In general, the normal modes in an open resonator system may not be orthogonal.However, it is possible that an open resonator system can support orthogonal modes if the system has high symmetry, and the modes have different symmetries, which will be the case for the numerical examples shown in the paper.A general discussion on the theoretical conditions for the open resonator to support orthogonal modes can be found in [21].

Single incident plane wave
We first consider the case where the incident wave consists of a linearly polarized plane wave , where ⃗ k is the wave vector in vacuum with a magnitude of k = ω/c, the total magnetic field outside the scatterer is, Using the temporal coupled-mode analysis, we can obtain the total scattered and extinction powers as P scat = c † scat c scat and , respectively.Following the definition of the scattering and extinction cross-section [19], σ scat = P scat /I inc and σ ext = P ext /I inc , respectively, we obtain their expressions, As a first example shown in Fig. 1(c), we consider the incidence of a plane wave that propagates along +x direction onto a sub-wavelength cylindrical scatterer with azimuthal symmetry.This azimuthal symmetry imposes that B should take the form and the phase Φ −n = Φ +n ∈ [−π, π) by mirror symmetry.The resonant modes are standing waves.
The radial parts can be expressed as Bessel functions, and the angular parts are of the form cos(mϕ) or sin(mϕ).Taking into account the time-reversal symmetry and energy conservation constraint [21], we obtain the form of D matrix where D ni represents the coupling between the i-th sinusoidal resonant mode and the n-th outgoing exponential channel: Therefore, the n-th component of the outgoing wavefront vectors, as given by Eq. ( 6), is For a sub-wavelength scatterer, the resonant scattering from modes with |n|>1 as well as 0-th order mode is negligible.Therefore, the resonant scattering wavefront, where the channel vector component corresponds to the second term in Eq. ( 17), is composed of only two outgoing waves in the ±1 resonant modes, i.e., for otherwise, c out,n = (−1) n+1 c in,−n e jΦ n .Such scattered waves carry no net momentum along x-or y-axis, namely, ⃗ F res = 0 in Eq. (9b).This is also evident from the force operator ⃗ P in Eq. ( 10) that couples only the channels of neighboring orders (|n − m| = 1), but not between the +1 and −1 orders.Therefore only ⃗ F inter and ⃗ F bg , the respective force resulting from the background scattering and its interference with the resonant scattering, contribute to the resonant features of the force spectrum.
By considering the background scattering phases only from ±1 channels, as well as Φ −1 = Φ 1 , the optical force ⃗ F, and the extinction cross-section σ ext , are therefore: F y = 0 due to the mirror symmetry about the xz plane.We are ready now to illustrate the general lineshapes of optical force and extinction crosssection spectra as described by Eqs. ( 18) and (19), in both lossless and low-loss cases.For the lossless (γ p = 0) case, σ ext and σ scat are equal, while in the lossy case, σ ext includes contributions from both the scattered and absorbed power.Fig. 2 plots F x spectra with different background scattering phases Φ 1 for a lossless scatterer γ 0 1 = 0 in panel (a) and a lossy scatterer γ 0 1 = γ 1 in panel (b).From Eq. ( 19), σ ext can be obtained from F x spectra by multiplying with 1/c.The background scattering phase Φ 1 has a strong influence on the spectra of F x and σ ext as follows.
First, in the absence of background scattering, i.e., Φ 1 = 0, F x and σ ext spectra have a Lorentzian lineshape.At the resonant frequency ω 1 , F x and σ ext maximize to 2λ/πc and 2λ/π, respectively, which are the maximum possible values when the scattering process involves only the ±1 channels.In the presence of loss, the peaks of the optical force and the extinction cross-section are reduced.
Second, for a background scattering phase Φ 1 = π, the background scattering contribution, F bg , maximizes at all frequencies across the resonance.In the absence of loss, the Fano interference between the background scattering wavefront vector Bc in and the resonant wavefront vector Tc in causes a dip with null values in both F x and σ ext spectra at ω 1 .In the presence of loss, it is no Fig. 2. Optical force x component spectrum exerted on a small scatterer by a single plane wave given by Eq. ( 18).(a) presents a lossless case and (b) presents a lossy case where the internal and external loss rates are equal.The background scattering phase is chosen as small, intermediate, and dominant in the optical force spectra.The extinction cross-section spectra can be obtained from F x spectra by multiplying with 1/c.longer possible to completely eliminate the net optical force at a single frequency through the use of a Fano interference.Thus, the optical force and the extinction cross-section cannot reach zero.
Third, for all cases where Φ 1 is neither 0 nor π, e.g., Φ 1 = π/2, F x and σ ext spectra exhibit an asymmetric Fano-like lineshape, which varies asymmetrically in the vicinity of ω 1 .The value of Φ 1 determines the shift of both the minimum and maximum in the spectra of F x and σ ext with respect to ω 1 .
To validate the results discussed above, which were obtained just by assuming the prescribed symmetry and the existence of the resonances, we compare our results for two specific examples, with exact calculations of optical force and extinction cross-section spectra based on the Mie solution [19].For both examples, we consider a scatterer with the geometry shown in Fig. 1 (c), which consists of two concentric layers of dielectric and plasmonic materials.The permittivity of the dielectric core is ϵ d = 12.96, and the permittivity ϵ m of the metal [36] is described by a Drude model, ϵ m = 1 − ω 2 p /(ω 2 − jγ p ω) where ω p is the plasma angular frequency and γ p characterizes the internal loss rate, which is equal to zero for a lossless scatterer.The influence of background scattering is studied by varying the size of the scatterer.
The first case consists of a small scatterer with the geometrical parameters of r 1 = 0.2λ p , r 2 = 0.135λ p , where λ p = 2πc/ω p .These geometrical parameters are chosen so that the background scattering is weak.As we will see in the discussion below, this system exhibits the resonance at ±1 modes.Throughout the paper, we refer to this structure as the "small scatterer".The optical force component along x-axis is plotted in Fig. 3 (a).For such a small scatterer, as described by Eq. ( 18), the spectrum peaks near 0.1552ω p with a Lorentzian lineshape.From Eq. ( 19), the extinction spectrum is proportional to the force spectrum in the frequency range under study.
To quantitatively compare the exact solutions with the coupled-mode results, we extract the coupling parameters by fitting the spectrum of c out with the parametric response from the coupled-mode theory in Eq. ( 17).This fitting is done by applying the least-square method to minimize the difference between the numerical and analytical c out spectrum for all frequencies.The coupled-mode parameters are ω 1 = 0.1552ω p , γ 1 = 2 × 10 −4 ω p , Φ 1 = 0.0516.The fact that Φ 1 ≪ π verifies that the background scattering is weak.Using these coupled-mode fitting parameters, the analytically fitted optical force spectrum agrees excellently with the exact Mie solutions as shown in Fig. 3(a).
The second case is the same lossless structure with a new set of geometrical parameters: r 1 = 0.76λ p , r 2 = 0.27λ p .As we will see in the discussion below, this system exhibits the resonance at ±1 modes and a strong background scattering from multiple channels.Throughout the paper, we refer to this structure as the "large scatterer".The spectrum of the optical force is plotted in Fig. 3 (b), as obtained from exact Mie solutions of the scatterer.We observe that in the frequency range under study, the optical force spectrum exhibits asymmetric lineshapes due to the Fano interference from the background scattering.
As before, we obtain the coupled-mode parameters from a least-square fitting procedure.The parameters for the ±1-th order resonances are ω 1 = 0.1568ω p and γ 1 = 3 × 10 −5 ω p .The resonance frequencies of all other orders are sufficiently far from ω 1 , thus they do not contribute significantly to the lineshape of the spectrum which is considered here.Thus, we set all other parameters to 0. For the background scattering, the fitted parameters are Φ 0 = −0.4390,Φ 1 = 0.7773, Φ 2 = 0.065.From these scattering phases, the matrix B in Eq. ( 6) can be determined.With these coupled-mode fitting parameters, the analytically fitted spectra agree well with the exact Mie solution.In comparison with the small scatterer, the large values of Φ 0 and Φ 1 cause an asymmetric Fano-like lineshape around ω 1 .
The theoretical analysis allows us to link the properties of the scatterer to a few other aspects of the scattering and force behaviors.When only a single pair of channels excite the resonant mode with the same angular mode index, the force and the extinction spectra are strictly proportional, as we see in the small scatterer case.For the large scatterer, these two spectra are no longer proportional since multiple channels with different angular indexes contribute.For both small and large scatterers, the force is positive across all frequencies, i.e., the direction of the force is parallel to the Poynting vector as expected for the case of a single incident plane wave.The strong background scattering, which results in the Fano-like lineshape, does not change the sign of the force -thus only repulsive (or zero) force is exhibited here.

Two incident plane waves
The behavior of the optical force is significantly richer when the incident wave is more complex than a single plane wave [37,38].As an illustration, we proceed one step in complexity to study the optical force spectrum with two incident plane waves.
The two waves have the same frequency and amplitude and propagate along angles ±θ measured from the x-axis, as schematically shown in Fig. 1 (c).Additionally, the two waves can have a relative phase ϕ.The incoming wavefront is represented in the Hankel basis by c in,n = (−1) n+1 j n 2 (︁ e −jnθ + e j(φ+nθ) )︁ . For such an incident plane wave, the Poynting vector is always along the x-direction.For the small scatterer where the ±1-th order modes dominate, the scattered wavefront vector c scat is given by Eq. ( 7).
Due to the gradient force contribution [39] from the standing-wave pattern [40], the optical force is in general not parallel to the incident Poynting vector.By applying the optical force operator in Eq. ( 8) and normalizing by the intensity I inc of a single plane wave, the optical force is given by: ]︄ sin ϕ sin θ cos 2θ.
(21) Therefore, the lineshapes of the x and y components of the force include an odd function of (ω − ω 1 ).In this case, the F x is always positive, for all incident wavefront θ, ϕ and background scattering phase Φ 1 and it is not possible for the force spectrum to change the sign for the small scatterer in this region.In order to observe negative force, the background scattering phase needs to be nonzero in multiple channels.
As the force spectra lineshape depends on the angle θ and phase difference ϕ, it is possible to control the force spectrum by changing the incident wave angles and relative phases.We consider the force on the small and large scatterers that were discussed in the previous section.We study the optical force when the incident angle θ = π/3 and the relative phase of the two incident plane waves at the origin is ϕ = π/2.The results for the small scatterer are depicted in Fig. 4(a, b).In the frequency range under study, the force component along the Poynting vector, i.e., F x , exhibits a nearly Lorentzian spectral lineshape.This is consistent with Eq. ( 20) with a weak background scattering.We recall from the previous section that Φ 1 = 0.0516 ≪ π for this case.On the other hand, in panel (b), the force component orthogonal to the Poynting vector, i.e., F y exhibits an asymmetric lineshape, where the force is very close to an odd function of (ω − ω 1 ).The sign of this gradient force can change from positive to negative, as the frequency varies across the resonant frequency.Using the same coupled-mode fitting parameters extracted earlier, we again see an excellent agreement with the numerical results.
This force behavior for a lossless structure persists in the presence of a modest amount of loss.Figure 4 (c, d) shows the optical force spectra when a material loss characterized by γ p = 2 × 10 −3 ω p is introduced to approximate the property of gold [41].Since γ p ≪ ω 1 , the previous coupled-mode fitting parameters still apply, including the resonant frequency ω 1 , the background scattering phase Φ 1 , and the external loss rate γ 1 .Then, an internal loss rate is determined as γ 0 1 = 9 × 10 −4 ω p by repeating the same parameter fitting process as outlined above for the lossless case.The main features of F x and F y spectra are still present with this intrinsic loss.However, compared with the lossless case, we observe that the peak value is diminished and the bandwidth is broadened as can be expected.
The interplay between the background and resonance scattering can enable diverse lineshapes of the force spectrum.To this end, we again consider the large scatterer in the absence of material loss.The phase and angle of incidence of waves 1 and 2 are chosen to be the same as for the small lossless scatterer described above.As shown in Fig. 4(e, f), the dots are fitting results using Eqs.( 20) and ( 21), with the same parameters as above when considering the single incident plane wave case.The solid lines are exact numerical calculations via the Mie solution.We again see an excellent agreement between these results.
For the large scatterer, the optical force spectra along x direction exhibit Fano-like lineshapes, shown in Fig. 4(e, f) in the frequency range under study.We recall that the large background scattering phases of Φ 0 = −0.4390and Φ 1 = 0.7773 result in a Fano-like lineshape around the resonant frequency ω 1 .For F y spectra, due to the stronger background scattering, the asymmetric lineshape no longer has an odd symmetry around ω 1 .This example also verifies the applicability of the temporal coupled-mode model for large scatterers under complex illumination.
The next example is important because it pertains to the case where the force spectrum is substantially different from the extinction spectrum.For this, we use two incident plane waves that have different amplitudes (the resulting force equation is discussed in the Section 6.3 in the Appendix).We present in Fig. 5 an example of two such counter-propagating plane waves that achieve significant asymmetric force line shapes, while the extinction cross-section spectrum is still symmetric.The incident wavefront is composed of two waves as in Fig. 1 (c) where wave 1 is propagating along the +x direction, while wave 2 is propagating along −x direction.The magnitude of wave 2 is 0.8 of that of wave 1 and wave 2 has a phase ϕ = π/2 delay relative to wave 1. Fig. 5. Extinction cross-section and optical force under two counter-propagating plane waves for the small scatterer with r 1 = 0.2λ p , r 2 = 0.135λ p without material loss.Wave 1 is propagating along the +x direction; while wave 2 is along −x direction.The magnitude of wave 2 is 0.8 of wave 1 and wave 2 has phase ϕ = π/2 delay compared with wave 1.The optical force along x-axis (a), and extinction cross-section (b) are calculated by the Mie solution and the temporal coupled-mode theory fitting.The blue solid lines are the Mie theory solution, and the red dotted line is from the temporal coupled-mode theory fitting.
As shown by the spectra, in the frequency range under study, the extinction cross-section shows a Lorentzian line shape.In contrast, the force component along the Poynting vector, i.e., F x spectrum, shows an asymmetric line shape and varies asymmetrically near the resonant frequency ω 1 as shown in panel (b).We interpret the lineshape based on our analysis in this section, which takes into account the total force including both the gradient and radiation forces.With the two counter-propagating incident waves of similar magnitudes, the radiation force is reduced significantly.Furthermore, due to the small size of the scatterers where the background scattering is zero, the only contribution for the asymmetry of the line shape is the gradient force.This is the first demonstration in this paper where the gradient force is present along the direction of the incident Poynting vector.

Size-sorting mechanism
To further explore the implications of the phase control over the force spectrum, we fix the incident angle to π/3 and place the scatterer center at an arbitrary point (x, y) to obtain different relative phases (measured at the center of the scatterer).The two incident plane waves are selected to have no phase difference at the origin, thus at any observation point (x, y) on the z = 0 plane, the local phase difference of the two incident plane waves is ϕ(x, y) = 2ky sin θ.We calculate the force exerted on the scatterer at different locations based on the temporal coupled-mode theory Eqs. ( 20) and (21).We find null values of F y located at points along y according to the spectral resonances.These null points can be either stable or unstable along the y direction.From Eq. ( 21), one can see that the y component of the force is zero where the phase 2ky sin θ is an integer multiple of π.Similar force behavior has been seen in optical lattices [42].
In Fig. 6, we present the quiver plot of the force density in the xy plane exerted on the previously described small scatterer in Fig. 6 (a) and large scatterer in Fig. 6 (b), under two incident plane waves.As an example, the frequency is chosen as ω = 0.1565ω p , near the resonance frequencies of the force spectrum.Given this frequency, the y locations that fulfill the geometric condition for null F y correspond to integer multiples of y 0 = 1.8446λ p .Material loss γ p = 1 × 10 −3 ω p is introduced to both scatterers to approximate the property of realistic materials.We assume that the scatterers are diluted to the desired concentration, such that each time we only have one scatterer in the field of view.Then, the optical force exerted on each scatterer can be considered individually.We observe for the small scatterer in panel (a), F y points towards the even multiples of y 0 , which are flow channels along the x direction that collect the small scatterers.On the contrary, for the large scatterer in panel (b), the odd multiples of y 0 are flow channels along the x direction that collect the large scatterers.The results depicted in Fig. 6 can be used conceptually for sorting scatterers (such as viruses combined with nanoparticles [43]) based on their sizes or more precisely on their resonant frequencies.To realize precise size sorting, the force spectrum should change the sign across its resonant frequency, and the resonant frequencies of scatterers with different sizes should be well separated.In our case, the small and large particles exhibit resonant frequencies around 0.155 ω p and 0.157 ω p respectively.Both scatterers exhibit dramatic sign changes in F y as indicated by Fig. 4 (b, f).Similarly, the signs of F y also vary as a function of y (i.e., phase difference) for the two scatterers on the spatial plane.Thus, different scatterers get collected at different flow channels.
This sorting mechanism is different from previous works on sorting co-streaming scatterers [44] with a five-beam interference pattern, sorting on-chip scatterers with whispering gallery mode [45] and sorting nanoparticles with near-field potential on-chip [46].The specific particle sorting we present here is restricted to highly elongated particles, with z direction stabilization (e.g. by an electric field) to prevent rotation.However, the main point is the fact we can provide the designer the generalized force maps as a function of very few parameters for specific applications.

Conclusion
We have developed a theoretical framework for analyzing optical force resonant lineshapes based on the temporal coupled-mode theory.This theory links optical properties of the structure to the properties of optical force resonances.With a single plane wave, we show the lineshape in the force and extinction cross-section spectra as determined by the resonances.With two plane waves, we demonstrate tunable asymmetric lineshapes of the force spectrum along and orthogonal to the incident Poynting vectors.We verify our theoretical analysis with numerical examples.Our theory highlights the engineering of the optical force [47] by utilizing the resonance from the scattering process, which applies to both lossless and low-loss systems.It would be of interest to apply the theory in three dimensions to treat the resonant forces on nanoparticles.And for this purpose, one should expand the fields in the vector spherical waves, but the basic formalism will remain the same.We have also proposed a potential method of sorting metallic cylindrical scatterers based on our theory.This theory may also be important for the application of precise nano-manipulation with optical tweezers and size sorting with biological sensors.

Force matrix derivation
In this section we provide a derivation of Eq. ( 3) and Eq.(10).We derive Eq. ( 3) first.Assume a magnetic field distribution ψ.(Similar derivation can be done using electric field distribution.)The computation of the force ⃗ F for such a field distribution can be written in a bilinear form as The bilinear form ⃗ O has the following properties: We now consider a magnetic field distribution ψ = −ψ * , which corresponds to the time reversal of the field distribution ψ.Since the force is even with respect to time reversal, the magnetic field distribution ψ should generate the same force as ψ.Therefore, Thus, the bilinear form satisfies Let us now consider the field distribution, The corresponding real-valued force is then Let us denote the (m, n)-th element of the matrix ⃗ P as Suppose we choose the outgoing wave basis such that ⃗ V out,m = − ⃗ V * in,m , i.e., the incoming and outgoing wave bases form a time-reversal pair.Further, as the force operator is Hermitian, we ).Then in Eq. ( 28), we apply Eq. ( 26) and obtain, On the other hand, momentum conservation requires that the integration of the Maxwell's stress tensor over the integration surface surrounding the scatterer, i.e., ⃗ F, be independent of r, i.e., the distance from the center of the scatterer.As an illustration, consider the cylindrical coordinate systems.For an integration surface far from the center of the scatter, the incoming and outgoing fields in the m-th channel can be approximated as 1 √ r e ±j ⃗ k m •⃗ r in its radial dependency, where ⃗ k m is the wavevector along the radial direction.With these fields, one can see that ⃗ F in and ⃗ F out has no explicit r dependency whereas ⃗ F int has an explicit r dependency.Therefore, ⃗ F int needs to vanish to maintain momentum conservation.Therefore, which proves Eq. ( 3).
In the following, we derive the explicit form of ⃗ P, i.e., Eq. ( 10), for the quasi-2D particles.One possible basis function in polar coordinates (r, φ) centered at the origin is ⃗ V in,m (r, φ) = −H (1)  m (kr)e −jmϕ ẑ, where H (1)  n and H (2)  n are the n-th order Hankel functions of the first and second kind, respectively.For the force calculation, Eq. ( 29) can be explicitly written as where the integration is taken on a circle centered on the particle.ϵ 0 and µ 0 denote the vacuum permittivity and permeability, respectively.⃗ E in,m and ⃗ H in,m are the electric and magnetic fields for the incoming wave in the m-th channel.⃗ H in,m = H 0 ⃗ V in,m with ⃗ V in,m from Eq. ( 32) and H 0 = √︂ ωϵ 0 2 as the normalization factor.Then ⃗ E in,m = 1 jωϵ 0 ∇ × ⃗ H in,m .By momentum conservation, we can take the radius R of the circle to be sufficiently large, such that ⃗ H in,m can be approximated as: Keeping only the terms with R − 1 2 order, the corresponding electric field is then Therefore, Eq. ( 34) becomes which gives the form of the force operator matrix in Eq. ( 10) as,

Temporal coupled-mode theory
In this Section, we provide some of the details for the temporal coupled-mode theory.The incoming waves can be transmitted completely to the outgoing channels.For the TM-polarized field, B 0 is expressed as, The term δ m,−n is attributed to the conservation of angular momentum between the incoming wave in the n-th channel and the outgoing wave in the −n-th channel.The (−1) n+1 factor can be seen from the expansion of a plane wave with the wavevector ⃗ k in cylindrical harmonics (Eqs.( 32) and ( 33)): where J n is the Bessel function of order n.Using the decomposition from Eq. ( 40) we obtain c in,n = − 1 2 j −n , c out,n = 1 2 j −n + c scat,n .In Eqs. ( 4) and ( 5), Ω, Γ, Γ 0 are all Hermitian matrices.The N-by-N matrix B in Eq. ( 5) describes the background scattering, i.e., the scattering of the incoming wave without the excitation of the resonant modes.In this paper, we assume that the background scattering is lossless and hence B † B = I.The matrices B, D, K, and Γ are related by: as imposed by time-reversal symmetry and Lorentz reciprocity arguments [21].We note that Eq. ( 41) is strictly valid only for a lossless scattering system but can approximate a system with low loss [23,24].Using the same temporal coupled-mode analysis as in Eqs. ( 18) and (19), we obtain the scattering cross-section spectrum described by γ 2 1 (1 + cos Φ 1 ) + 2γ 1 (ω − ω 1 ) sin Φ 1 (ω − ω 1 ) 2 + (γ 1 + γ 0 1 ) 2 . (44)

Arbitrary superposition of two plane waves
In the main text of the paper, we considered a cylindrical scatterer illuminated by two incident plane waves of equal amplitude.Now we generalize to an arbitrary (along θ direction from the x axis) incident TM-polarized wavefront.To compute the force from such a wavefront, we introduce a same-frequency plane-wave basis which is complete and orthogonal, c(θ), where c n (θ) = (−1) n+1 j n 2 e jnθ .For the force component along an arbitrary direction with an angle α from the x-axis (below as the α-direction), the corresponding force operator is denoted as P α .From Eq. ( 10), the matrix element of P α , indexed according to the channels m and n is, P α,mn = 1 ]︁ .Considering scatterers having ±1-th resonant modes, with the resonant frequency ω 1 , and external and internal loss rates, γ 1 and γ 0 1 , the exerted force along the α-direction is obtained by replacing ⃗ P in Eq. ( 3) with P α .To explore the symmetry of the lineshapes in the optical force spectrum, we decompose the force into the symmetric component F s and the asymmetric component F a , which is an even and an odd function of frequency difference ∆ω = ω − ω 1 , respectively: where Q α = −B −1 ⃗ P T α B. For convenience, we denote a new complete and orthogonal basis v(θ) = B † c * (θ).Using this new basis, in Eq. ( 45), an arbitrary incident TM-polarized wavefront c in can be decomposed as c in = ∫ dθq(θ)v(θ) where q(θ) = 2 π v † (θ)c in .We then can determine the direction α s at which the force has a purely even-symmetric spectrum, as well as the direction α a at which the odd-symmetric force component maximizes.By setting the odd-symmetric component F a of the force is zero, we obtain α s as,
For any direction other than that given by Eq. ( 46), the force spectrum exhibits some degree of asymmetry.By setting the derivative of F a with respect to α as zero, we obtain α a where the odd-symmetric component of the force maximizes, From α a and α s , for a small scatterer where the ±1-th modes dominate, there is always a direction α s at which the force exhibits an even-symmetric Lorentzian lineshape, as well as a direction α a at which the force exhibits a maximally asymmetric lineshape around the resonant frequency ω 1 .
As an application of α a and α s , we examine the case of two incident plane waves with the same amplitude where the propagating directions of the two waves from the x axis are denoted as −θ 1 and θ 2 , respectively.Two basis functions v(−θ 1 ) and v(θ 2 ) are enough to describe the incoming wavefront c in .Then, α a and α s are simplified to, α s = (θ 2 − θ 1 )/2 and α a = π/2 + (θ 2 − θ 1 )/2.This is consistent with the numerical observations in the previous sections.α s and α a do not depend on the detail of the scatterer and are controlled only by the directions of incident waves, provided that the ±1-th resonant modes dominate the optical response.

Fig. 1 .
Fig. 1.Optical systems consisting of a cylindrical scatterer under (a) a general incident wavefront for an ensemble of scatterers, (b) a general incident wavefront for a quasi-2D scatterer, and (c) one plane wave as the incident wavefront.The spherical surface (a) and dashed lines (b,c) are the integration surface used to evaluate Maxwell's stress tensors.(d) Schematic of an optical resonator system coupled with multiple incoming and outgoing ports.The field plots represent the incoming and outgoing channels V in(out) from Eqs. (32) and (33).The time-averaged optical force ⃗ F on the scatterer, as shown in Fig. 1 (a), is calculated via an integral of the time-averaged Maxwell's stress tensor

Fig. 3 .
Fig. 3. Optical force and extinction cross-section spectra of the lossless scatterer under single plane wave incidence.The composition of the concentric scatterer is shown in Fig. 1(b), with the permittivity ϵ d = 12.96 in the core, and ϵ m = 1 − ω 2 p /(ω 2 − jγ p ω) with γ p = 0 in the shell.(a) The optical force for a small scatterer (r 1 = 0.2λ p , r 2 = 0.135λ p ) along x direction.(b) The optical force for a large scatterer (r 1 = 0.76λ p , r 2 = 0.27λ p ) along x direction.The blue solid lines are the Mie theory solution, and the red dotted line is from the temporal coupled-mode theory fitting.

Fig. 4 .
Fig. 4. Optical force with two incident plane waves at a relative phase ϕ = π/2 and the incident angle difference θ = π/3.The composition of the concentric scatterer is shown in Fig. 1(c), with the permittivity ϵ d = 12.96 in the core, and ϵ m = 1 − ω 2 p /(ω 2 − jγ p ω) in the shell.(a, b) The optical force for a small and lossles (γ p = 0) scatterer with r 1 = 0.2λ p , r 2 = 0.135λ p .(c, d) The optical force for a small and lossy (γ p = 2 × 10 −3 ω p ) scatterer.(e, f) The optical force for a large and lossles (γ p = 0) scatterer with r 1 = 0.76λ p , r 2 = 0.27λ p .The blue solid lines are the Mie theory solution, and the red dotted line is from the temporal coupled-mode theory fitting.

Fig. 6 .
Fig.6.Quiver plots of the horizontal force experienced by a cylindrical scatterer under two incident plane waves at frequency ω = 0.1565ω p .The geometry is as Fig.1(b) with θ = π/3 and waves 1 and 2 have the same phase at the origin.(a) The small scatterer with r 1 = 0.2λ p , r 2 = 0.135λ p is collected at y = 0 and y = ±2y 0 if the scatterer is in a flow channel, where y 0 = 1.8446λ p , as shown by the small red dots (drawn to scale) on the left.(b) The larger cylindrical scatterer with r 1 = 0.76λ p , r 2 = 0.27λ p is collected at y = ±y 0 if the scatterer is in a flow channel, as shown by the large red dots (drawn to scale) on the left.Both scatterers have a material loss γ p = 0.001ω p .