Optical chiral sorting forces and their manifestation in evanescent waves and nanofibres

Optical fields can exert forces of chiral nature on molecules and nanoparticles, which would prove extremely valuable in the separation of enantiomers with pharmaceutical applications, yet it is inherently complex, and the varied frameworks used in the literature further complicate the theoretical understanding. This paper unifies existing approaches used to describe dipolar optical forces and introduces a new symmetry-based `force basis' consisting of twelve vector fields, each weighted by particle-specific coefficients, for a streamlined description of force patterns. The approach is rigorously applied to evanescent waves and dielectric nanofibres, yielding concise analytical expressions for optical forces. Through this, we identify optimal strategies for enantiomer separation, offering invaluable guidance for future experiments.


I. INTRODUCTION
The idea that light may exert a force on matter was proposed by Johannes Kepler in the 1600s to explain the observation that a comet's tail always points away from the sun.Using the phenomenon of optical forces for the levitation and control of small objects in a designed way had to wait many centuries.It was first explored by Arthur Ashkin in the 1970s, who proposed and developed the optical trap [1][2][3] (also called optical tweezers), earning him the 2019 Nobel Prize in Physics.These optical tweezers proved instrumental in the micromanipulation of molecules, cells, viruses, and atoms and have revolutionised research fields with the advent of laser cooling and quantum control of macroscopic objects.
In the optical manipulation of small particles, two main optical forces are considered: the optical gradient, which tends to attract particles towards the maxima in electromagnetic energy, and the optical pressure force, which tends to push particles in the direction of light propagation.Optical tweezers exploit the gradient force, usually stronger for small particles.By tightly focusing a laser beam, a particle is trapped in the high-intensity focus of the beam.The stronger the focusing, the higher the intensity gradient, and the stronger the trapping forceincreasing the stiffness of the trap.
Given the success in the manipulation of matter by optical forces, a question has arisen in recent years: can light be used to sort chiral particles and molecules, in other words, to separate two enantiomers (which, apart from being a mirror-reflection of one another are otherwise identical)?This is a challenging problem precisely because two enantiomers are identical in many ways, and because-our body's building blocks being chiral molecules-the opposite enantiomers of a given pharmaceutical molecule may have drastically different effects [4][5][6].The answer to the question is positive-if light contains some chirality itself.
The expressions for chiral optical forces have been given in many works [7][8][9][10][11][12][13][14][15][16], and different proposals for chiral separation of enantiomers using light have been put forth [11][12][13][14][15][16][17][18][19], with experimental success found only for chiral microparticles that are large compared to the wavelength [17].A new researcher in this field reading the literature is faced with a difficult task due to the many alternative ways of expressing the forces because of the different notation and symbols, different unit systems, different definitions, and most importantly, different ways of grouping terms that make the expressions look very different from one another.The expressions and theoretical analysis of these chiral optical forces are not as straightforward as the gradient and pressure forces for conventional optical manipulation.There are many terms, and great subtlety is needed.In this work, our aim is threefold: first, to be as clear as possible, going over the main different forms of the expressions used in the literature for chiral optical forces in dipolar particles, pointing out how and why they are equivalent, and making sure they are as general as possible.Second, in the process of clarifying the expressions, we uncover and put forward a novel way of expressing the force, which brings additional clarity, as it highlights the different symmetries involved in the phenomenon.We write the net force as a linear combination of different vector force fields, acting as a basis that describes all possible force patterns, each of them weighted by particle-dependent coefficients.Third, we exemplify how this simplified notation elegantly applies to two simple example cases: an evanescent wave and a cylindrical dielectric nanofibre.In the process, we present generally valid, but surprisingly concise, ex-pressions for diverse electromagnetic quantities in these modes, and based on this formulation, describe the best strategies for enantiomer separation available in these geometries.
This paper is intended to serve as a valuable aid to experimentalists trying to achieve chiral separation of small (dipolar) enantiomers.As such, and contrary to a large fraction of the existing literature, we shall adopt SI units (Système international d'unités) together with the convention for writing the speed of light in vacuum as c 0 , while the phase velocity of electromagnetic radiation in an arbitrary medium will be denoted c = ω/k = c 0 /n = 1/ √ µε.Throughout the paper, we further assume a dispersion-free, lossless linear background medium and monochromatic fields.We shall denote time-dependent real fields with scripted letters, e.g., A(r, t), while their time-independent phasor representation with regular Latin letters, A(r), such that A(r, t) = ℜ[A(r)e −iωt ].One of the advantages of this representation is that any time-averaged quantity quadratic in the fields can be obtained by simply ⟨AB⟩ = ℜ[A * B] rather than integrating over time.
We organise this work as follows.In section II we review the most general expression for a time-averaged optical force on a small particle.In section III we provide a clear overview of all alternative but equivalent ways of writing this expression for the case of a linear bi-isotropic dipolar particle, and we introduce a new classification of the forces based on the symmetries of the particles and the fields.In the same section, we also introduce the concept of force basis, which can be used to design the fields for separating particles which break certain symmetries and in sections IV and V we apply this concept to two simple yet practical examples of evanescent waves and dielectric nanofibre, respectively.

II. OPTICAL FORCES
The most infallible tool at our disposal to calculate optical forces is to apply the law of conservation of momentum.Any mechanical momentum that a particle acquires must come at the expense of the reduced electromagnetic momentum in the fields.Therefore, in a steady-state time-harmonic field, the net time-averaged force acting on any piece of matter can be computed as the total electromagnetic momentum entering, per unit of time, any closed volume that includes the object.This is achieved by integrating the momentum flux over a closed surface that surrounds the volume.Derived from first principles [20][21][22], the time-averaged flux density of electromagnetic momentum is given by the Maxwell stress tensor: where E and H are the total electric and magnetic fields, the operator ⊗ is the outer (or tensor) product of two vectors, I is the 3 × 3 identity matrix, an asterisk indicates complex conjugation, and ε and µ are the absolute permittivity and permeability of the medium, respectively.
Knowing ⟨T ⟩, one can compute the net time-averaged force acting on any material body via the flux integral: where ds = n ds is the surface form, n is the outward normal vector of a closed surface S enclosing the body and ds is the surface element.We will drop the angular bracket notation for time averages, and hereafter, each observable will be assumed to be time averaged.This general method works on any body, regardless of its size.In this work, we are interested in small particles that undergo Rayleigh scattering and hence behave as dipoles.The criterion for which technique applies on which domain is generally described as follows, where k = ω/c = 2π/λ is the medium wavenumber and a is the particle size: In the Rayleigh limit, any particle under illumination scatters light in an identical form to an electric p and/or a magnetic m dipole.Because the analytical expression of the electric and magnetic fields of a general combination of an electric and magnetic dipole is well known, one can substitute such fields into eq.( 2) and, analytically performing the integration [23][24][25], arrive at the expression for the force acting on a dipole in the Rayleigh limit where p and m are electric and magnetic dipole moment vectors, considered as phasors of time-harmonic dipoles, E and H represent some external applied electric and magnetic fields (excluding those produced by the dipoles) and η = µ/ε is the impedance of the surrounding environment.Equation ( 4) constitutes the most general dipole force equation.All the following equations for the force are derived from eq. ( 4) by assuming the linearity of the particle, such that the induced dipoles are proportional to the applied fields.The linear response of a particle to the applied external fields can be described using complex polarisabilities.These polarisabilities (α e , α m , α c , α t ) are, in general, second-rank tensors that can be represented by 3 × 3 square matrices.The linear response of dipole moments to external fields can be written as follows: where components of these electric-magnetic vectors are chosen so that they have identical units, e.g., [ The literature is not unified on the definitions of these polarisabilities.In particular, different authors choose different units for these polarisabilities by choosing whether they are direct responses to external fields E and H (used, for example, in [7,9,[13][14][15][16]18]) or E and B (used in [11,12,15]).In this work, following Mun et al. [10], we define all polarisabilities so that they all have dimensions of volume, 3 , making the force expressions more elegant.Polarisabilities are a characteristic of the given particle and typically are proportional to its volume.In this work, we keep the most general form of polarisability, allowing the electric-magnetic coupling terms to have both an antisymmetric part (the chiral polarisability α c ) and a symmetric part (the non-reciprocal polarisability α t ).The latter is usually neglected (Bliokh et al. [8] includes it but neglects the recoil force) because it only exists for non-reciprocal particles that break timereversal symmetry.We keep this response not only for the sake of generality, but also because it completes a clear framework relating the symmetries of the particle and the optical fields to the forces that we can expect.
Knowing the incident fields and the polarisability of the particle, the exact calculation of the optical force is straightforward in two steps: (i) calculate the induced dipoles using eq.( 5), and (ii) substitute them into eq.( 4).This substitution may be carried out analytically to arrive at a single-step equation, but the resulting expression is only manageable under certain simplifying assumptions, as discussed in the following section.It is also worth pointing out that if the particle is in an environment that reflects the dipolar fields back into itself, this must be taken into account in both steps.In the first step, the dipole moments must be calculated selfconsistently, as the incident fields appearing in eq. ( 5) have a new term, the self-reflected fields, which depend on the dipole moments themselves.This can also be addressed through the concept of an effective polarisability that depends on the environment [26,27].In the second step, the dipolar fields that are reflected back must be included in the fields in the interaction terms of eq. ( 4).This complex interaction allows for a plethora of different phenomena to arise, such as repulsion from nearby surfaces [28][29][30][31], pulling forces [26,32], self torques [33], trapping forces [34], and lateral force recoils due to excitation of nearby guided modes [13,[35][36][37].These cases greatly complicate the calculations, but ultimately follow the same physics as included in eqs.( 4) and (5).In the rest of this work, we neglect back-reflection, as if the particle was in a homogeneous background environment.

III. FORCE ON BI-ISOTROPIC DIPOLAR PARTICLE
In practice, the polarisability tensors will depend on how the particle is oriented in space relative to the incident fields.However, unless the particles in question are somehow aligned, for instance, by some external static field or by nanofabrication, we can assume that for an average randomly-oriented particle within the sample, the effective polarisability will be (bi-)isotropic, and therefore can be described by complex (pseudo)scalars (α e , α m , α c , α t ).In this case eq. ( 5) becomes Effective polarisabilities for the specific case of a spherical chiral bi-isotropic particle can be found in appendix A. Equation ( 6) can be substituted into eq.( 4), and after separating the real and imaginary parts of polarisabilities and further simplifications, one can write the force in terms of eight local observables that are quadratic in the fields.These quantities are: two real scalars which represent electric and magnetic energy densities, respectively; a pair of real pseudovectors representing electric and magnetic spin angular momentum densities; and a complex pseudoscalar and vector where ℜW c /ω = S is the time-averaged helicity density, and ℜΠ is the time-averaged Poynting vector representing the flow of active power.Note that quantities ℑW c and ℑΠ, while often overlooked, are non-zero in general and also lead to forces felt by particles that break certain symmetries.The vector ℑΠ represents the flow of reactive power [20], while the pseudoscalar ℑW c represents a static "chirality"1 due to ℜE ̸ ⊥ ℜH and ℑE ̸ ⊥ ℑH, in this sense, it can be thought of as a measure of instantaneous collinearity of the electric with the magnetic field [8].
As mentioned before, the resulting force can be written in more than one way.One way is to use the complex particle polarisabilities to define the real scalar quantities that are often for their dimensions called cross-sections: where σ ext is indeed an extinction cross-section.For chiral forces, we can also define pseudoscalar cross-sections: These can be used to write the force as in [10,13,15,38] where U is a scalar potential defined as The chiral terms of this force are those that are a product of pseudoscalar polarisabilities α c and cross sections γ with pseudovectors such as ∇W c , S e , S m and ∇ × Π.
If we reflect either the particle or the fields in a mirror, the pseudoscalars or their corresponding pseudovectors change sign, respectively.The rest of the force is achiral.Some of the literature that uses eq. ( 12) gives names to individual terms.The potential U is sometimes called the free energy and ∇U are collectively named the gradient forces.Forces proportional to ℜΠ are referred to as radiation pressure forces 2 while the term with ∇ × ℜΠ is known as the vortex force.Spin density forces are the terms directly proportional to S e and S m , while spin-curl forces are proportional to their curls.The last term to have a name is ℑΠ which is called the flow force because it is related to the alternating flow of "stored energy" [20].The term related to the curl of this vector is only present for nonreciprocal particles, it is not included in the literature works and hence does not have a name.While eq. ( 12) is an elegant expression relating kinetic properties of light [eqs.(7) to (9)] with observables related to the particle [eqs.(10) and (11)], it has one disadvantage: the expression no longer allows distinguishing the interaction and recoil parts of the dipole force that appeared in eq. ( 4).In order to recover the two terms, we introduce two more vectors: which are often called canonical (or sometimes orbital) linear momentum densities, on top of which we may define a 'complex chiral momentum density' pseudovector: 2 even though they include the self-recoil of the dipole, where ℜp c is the chiral and ℑp c is the magnetoelectric momentum densities introduced in Bliokh et al. [8].These canonical momentum densities are directly related to the previously defined (pseudo)vector observables Π, S e , S m and their curls as follows [39,40]: which, when substituted into eq.( 12), allows us to rewrite the force in a different way.Confusion may arise due to the fact that a fraction of the literature uses eq. ( 12) while another uses eq. ( 17) below.The two alternative expressions give the same value for the force, but the terms are grouped differently (and some terms in one expression are divided into two terms in the other).While eq. ( 12) is simpler in terms of the electromagnetic quantities involved (so we may call it the EM-centric formulation), eq. ( 17) is simpler in terms of the particle polarisabilities involved (with no cross sections needed, so we may call it the particle-centric formulation).This alternative formulation is used by [9,14,16,18,23] but missing the terms involving α t .To our knowledge, no previous publication lists both non-reciprocal interaction terms and the recoil force.The full force reads: where we can see three distinct terms.The first two are generalisations of the conservative gradient and radiation pressure forces, while the third term is purely the recoil force from eq. ( 4).The advantage of this form is that it is easier to interpret.Notice that the gradient forces will always point towards (or away from) the maximal value of energy densities, and they depend only on real polarisabilities, which means that a particle will feel it even if there is no absorption of photons, it is a conservative force.The radiation pressure points along the canonical momenta of the wave and depends on imaginary polarisabilities, hence extinction cross sections, and it will require photons to be absorbed or scattered.The recoil force points along the flow of active power ℜΠ, reactive power ℑΠ and the electric S e and magnetic S m spin angular momentum densities and is always proportional to squares of polarisabilities, while the other two terms were linearly proportional to the polarisabilities.
One more formulation is worth mentioning for its simplicity and relationship with symmetries, which was introduced in Bliokh et al. [8] (although neglecting the recoil force).The key observation is that in free space, Maxwell's equations are symmetric under the parity inversion (P): (x, y, z) ⊺ → (−x, −y, −z) ⊺ , time reversal (T): t → −t, and the the duality transformation (D): which is a symmetry leading to the conservation of optical helicity [41].One can then redefine all the quantities involved in the force such that they are either symmetric or antisymmetric under these transformations.We shall label these quantities with an index A ∈ {0, 1, 2, 3}, such that they have the same behaviour under these transformations as Stokes parameters, i.e., A = 0 will be symmetric under all these, A = 1 will be antisymmetric under (D), A = 2 will be antisymmetric under all and A = 3 will be antisymmetric under parity.The energy densities appearing in the gradient forces can be written as: the canonical momenta featured in the pressure force: and the spin-like quantities in the recoil force: Notably, the quantities in eqs.(19) to (21) are not independent.In particular, the energy densities eq. ( 19) are time components of momentum densities eq. ( 20) in the relativistic sense, while the spin-like quantities eq. ( 21) can be thought of as fluxes of the aforementioned energy densities (for more details, see appendix B).One can show (see appendix C) that the chiral energy, momentum, and spin densities are simply differences between these quantities carried by the positive/negative helicity components of the electromagnetic field, i.e., This also justifies why they have dimensions of energy, momentum, and spin densities, respectively.The same procedure can be followed for the polarisabilities, which leads to the following expressions: We can also define square polarisabilities β A , which are equivalent to the recoil cross sections eqs.( 10) and (11): The minus signs in front of quantities with the label A = 2 are chosen to match the signs of Stokes parameters later; however, in eq. ( 24), they will cancel out.Expressed in this basis, the optical force takes a surprisingly simple form where we have included the recoil terms, which were not previously calculated.While eqs. ( 12), ( 17) and ( 24) are only true in the dipole approximation an advantage of the symmetry-based approach is that we can always divide the force into components, F = F 0 + F 1 + F 2 + F 3 , such that F 0 is present for every particle and the remaining components F A appear only for particles breaking corresponding symmetries.We remark that eqs.( 12), ( 17) and ( 24) are all alternative formulations of the same force.Ultimately, one may notice that the force on a dipolar particle can be written as a sum of (pseudo)scalar coefficients multiplied by a (pseudo)vector basis.In other words, one can write where eqs.( 12), ( 17) and ( 24) are simply different choices of the basis vectors and coefficients.For the representation in eq. ( 24) we list all the coefficients and the basis in table I. Note, however, that this is not formally a vector basis unless we neglect recoil terms, since the coefficients in eq. ( 23) are products of the other coefficients.Despite that, the concept illustrated by eq. ( 25) is quite profound: it tells us that we have coefficients λ i that carry all the dependence on the particle (and have units of volume), which multiplied by the basis vectors V i that depend only on the incident electromagnetic fields (and whose dimensions are force per volume or force density) gives us the force that a given particle feels.This can be used in designing optical forces to separate particles with different symmetries.For example, using this idea, one can plot different basis vectors V i (r) for a given electromagnetic mode (as we will do later) to get an idea of the types of forces (the directions and amplitudes) that are possible in a particular waveguide-independently of the particle.Looking at table I, one can see that both the force densities V i and associated coefficients λ i have the same Coefficients λi (dimensions of volume) Basis Vi (dimensions of force/volume) symmetries broken by the particle scalar pseudoscalar vector pseudovector

Interaction force
Gradient Showing the (pseudo)scalars and (pseudo)vectors involved in the dipolar force and illustrating the idea of a vector basis for the force.The last column lists what symmetry (D for duality, P for parity and T for time reversal) has to be broken by both the particle and the wave for the force component to be non-zero.When the recoil force is negligible, we can say that the force forms an eight-dimensional vector space spanned by basis Vi.Formally, the addition of the recoil terms does not promote the force into a twelve-dimensional vector space because the coefficients of the recoil terms are fully determined by the products of the coefficients of the interaction force.
behaviour under symmetries of the electromagnetic field, i.e., if one breaks the symmetry the other does as well such that their product (i.e. the force term) is a dual symmetric vector (hence parity odd) which is time reversal even/odd for interaction/recoil force respectively.This gives the formulation in eq. ( 24) a very useful interpretation.The terms associated with the total energy W 0 , the total canonical momentum p 0 , and the active power ℜΠ will be symmetric under D, P, and T, and will be felt by all particles.Other terms will be felt only by particles that break relevant symmetries (see the last column of table I).For example, a non-magnetic dielectric particle breaks the D symmetry, so it will also feel D asymmetric forces associated with W 1 = W e − W m , p 1 = p e − p m and the reactive power ℑΠ.If we had a particle that breaks duality in the opposite sense (a dominant magnetic response), then it would feel these forces in the opposite direction, and one could design an optical field to separate electric and magnetic particles.In our case, we are interested in chiral forces that break the P symmetry, i.e., forces associated with ℜW c = W + − W − , ℜp c = p + − p − and spin recoil forces S e and S m (note that a dual symmetric chiral particle would only feel S 0 = S e + S m ).Therefore, these forces are at our disposal if we want to separate enantiomers.As before, each enantiomer will feel a chiral force of the same magnitude but in the opposite direction.As a curiosity, a non-reciprocal particle, i.e., a particle that breaks the T symmetry (as well as D and P), would additionally feel two new force densities ℑW c and ℑp c .

IV. EVANESCENT WAVE
We may now apply the formalism to simple though practical examples.Consider a particle illuminated by an evanescent wave in a linear lossless medium with a complex wave vector k = k ′ + ik ′′ = k p k′ + iκ k′′ , where κ > 0 is a real parameter that describes how evanescent the wave is. Figure 1 illustrates the orthonormal directions of k′ (phase advance) and k′′ (decay direction) in an evanescent wave.Because k The electric field satisfying the transversality condition ∇ • E = 0 can be written as: where A s and A p are arbitrary complex amplitudes of two allowed polarisations, while the magnetic field can be acquired from the Maxwell-Faraday equation ∇ × E = iωµH.The duality transformation of eq. ( 26) can be written as follows: All observables are quadratic in the fields.Hence, the dependence on polarisation can be written in terms of Stokes parameters, which can be defined in terms of the polarisation amplitudes Since A p is a complex scalar and A s a pseudoscalar, we can say that S 0 and S 1 are scalars, while S 2 and S 3 are pseudoscalars.Additionally, one can confirm using eq.( 27) that S 1 and S 2 are dual asymmetric while the rest are symmetric.The energy densities for this field can be calculated using eqs.( 7) and ( 19), leading to an elegant form, simply proportional to the Stokes parameters: where f (r) = exp −2k ′′ • r /2k 2 .Notice that if a wave is not evanescent κ = 0, then W 1 = W 2 = 0, which is directly related to the fact that a plane wave-unlike an evanescent wave-does not break dual symmetry nor time-reversal symmetry.It is also apparent that the parity can be broken even by a circularly polarised plane wave since W 3 is proportional to S 3 .Knowing the energy densities, remarkably simple equations lead to the energy density gradients ∇W A = −2k ′′ W A , and to the canonical momenta p A = k ′ W A /ω, as defined in eq. ( 20).Substituting these into eq.( 24), we obtain a very simple expression for the interaction force (gradient and pressure terms): which depends only on the energy densities (Stokes's parameters) eq. ( 28) and relevant polarisabilities.The interaction force has no component in the lateral k ′ × k ′′ direction; forces in this direction are coming purely from the recoil part of the force, which depends on the spin quantities eq.(21).It is common knowledge that evanescent waves possess transverse spin [42], which is independent of helicity (S = W 3 /ω).The total spin can be calculated from eqs. ( 8) and ( 21), and has two components, longitudinal and transverse: This transverse spin will lead to a chiral lateral recoil force, as used in Hayat et al. [14].It is also known that the real Poynting vector acquires transverse components in the case of a circularly polarised evanescent wave (see, e.g.Wei and Rodríguez-Fortuño [43]) explicitly given by Therefore, even an achiral particle will feel an achiral lateral force if irradiated by a circularly polarised evanescent wave, as confirmed experimentally by Antognozzi et al. [44].In addition to transverse spin (contained in S 0 ) and transverse Poynting vector (contained in S 3 ), the two remaining spin-like quantities from eq. ( 21) (which are the reactive power proportinal to S 2 and the difference in electric and magnetic spins S 1 ) will similarly have transverse contributions: notice that this time, the other contributions point in the direction k′′ rather than k′ .The full recoil force (third term in eq. ( 24)) is then The total force F = F int + F rec exerted by an evanescent wave on a particle can be decomposed into the three directions given by the evanescent wave propagation, decay, and lateral directions.This force depends only on the polarisation of the wave (Stokes' parameters), which controls the energy densities eq. ( 28), and on the evanescence factor κ, since k p = √ k 2 + κ 2 , and is explicitly written as follows: The chiral force will only consist of the components that have pseudoscalar coefficients, which (assuming a reciprocal particle, i.e., α 2 = 0) are Therefore, the chiral force of a reciprocal particle in an evanescent wave is: From this, one can see that in order to maximise the lateral 3 chiral force, one needs an evanescent wave with a linear polarisation S 1 = ±S 0 , which also forces nonchiral lateral forces to zero and is precisely the configuration proposed by Hayat et al. [14] and depicted in fig. 1, 3 Pointing in direction k′ × k′′ .
particle (a).Note that in that case, the transverse spinrecoil force performs the sorting.Hence, this method will be more successful for larger particles with non-negligible recoil force.One cannot perform lateral sorting for particles that are too small to feel recoil forces.If one wants to sort smaller particles that only feel gradient and pressure forces (F int from eq. ( 29)), one would have to choose a circular polarisation S 3 = σS 0 (where σ = ±1) and sort the enantiomers in the plane defined by the evanescent wave.This type of sorting of small enantiomers (suitable for chiral molecules) in an evanescent field is a novel proposal, to the best of our knowledge.The circularly polarised evanescent wave has non-zero helicity W 3 ̸ = 0, but unlike a plane wave, this helicity is not uniform but decays exponentially like the energy.This non-uniformity creates a helicity gradient force that will attract the enantiomer with σℜ(α c ) > 0 to the region with higher helicity and repel the other.This is shown in fig. 1, particle (b).The magnitudes of the total interaction force acting on the positive/negative enantiomer (the positive enantiomer is the one with ℜα c > 0) in terms of the chiral polarisability of the positive enantiomer α c will be For this method, typically the achiral force is stronger than the chiral forces, and so the particles within the evanescent wave will stick to the surface due to gradient forces.But one of the enantiomers feels this force stronger than the other, due to the additive or subtractive chiral force which reinforces or weakens the achiral force.This allows enantiomer separation proposals based on the difference in net gradient forces.If the chiral force is strong enough to overcome the achiral force, then one enantiomer would feel attraction, and the other repulsion, from the surface.This is illustrated in fig. 1, particle (b).Particles with non-negligible ℑα c will also feel a chiral pressure force in the direction of wave propagation [see fig. 1, particle (c)], which will introduce an angle θ between the two forces given by V. CYLINDRICAL NANOFIBRE Intuitively, the modal fields outside a nanofibre are equivalent to 'radial evanescent waves' and we can expect similar behaviour to the previous section, but with the added richness of having many modes.The eigenmodes of cylindrical waveguides are well studied, since they can be calculated analytically: a particularly concise formulation is provided in Picardi et al. [45].A dielectric nanofibre aligned with the z-axis can be modelled as a nonmagnetic medium of radius r 0 , characterized by its permittivity and permeability ε = ε 1 ε 0 , for ρ < r 0 ε 2 ε 0 , for ρ > r 0 and µ = µ 0 .
The solution to Maxwell's equations yields electric and magnetic fields that can exist in and around the dielectric fibre.As Picardi et al. [45] showed, it is convenient to choose the spin basis ê±1 = (x ± iŷ)/ √ 2 and ê0 = ẑ to represent the solution.On this basis, any vector can be expressed as F = F +1 ê+1 + F 0 ê0 + F −1 ê−1 , where the components F s are spin-weighted functions. 4The advantage of this representation is that rotation around the z-axis by an arbitrary angle θ involves just a phase shift in the spin basis by sθ, ês → e isθ ês , and a phase shift by −sθ in the vector components, F s → e −isθ F s [46].Picardi et al. [45] also showed that the total longitudinal angular momentum of the guided modes in the cylindrical fibre is quantised, that is, J z = ℓW 0 /ω, where ℓ is an integer.The electric field, E, of a single mode with a definite azimuthal number ℓ can be represented by where k is the wave number in the medium, k z > k 0 is the mode propagation constant (k 0 is the wave number in vacuum), κ = k 2 − k 2 z is the radial wave number.Function Z α (x) is defined piecewise as the Bessel/Hankel function of the first kind inside/outside of the fibre, respectively.The corresponding magnetic solution, H, can be obtained via the duality transformation: 4 spin-weighted functions have wide use in mathematical physics, particle physics and cosmology.In particular, vector spherical harmonics and Bessel functions can be generalised and simplified by replacing them with their spin-weighted counterparts [46].
The dispersion relationship for the propagation constant k z given parameters (ω, r 0 , ε 1 , ε 2 , ℓ) is found from the transcendental characteristic equation (see fig. 2), whereas the complex constants A, B both inside and outside the fibre are determined from the boundary conditions at ρ = r 0 (see appendix D).Notice that for a mode of a given ℓ (denoted by each colour in fig.2), there is Normalised radius (or frequency) Numerically calculated dispersion relationships for a dielectric fibre made of silicon nitride, ε1 = 4.3, immersed in water, ε2 = 1.7 (both values are at 20 • C valid for λ0 from 750 nm to 1750 nm [47]).We plot it for the angular momentum eigenmodes ℓ ∈ {0, ±1, ±2} (in blue, orange, and green).Each ℓ has an additional mode number n.We plot values n ∈ {0, 1, 2} (using solid, dashed, and dot-dashed lines).For ℓ ̸ = 0 the mode number n is related to the spin number of the dominant mode σ = (−1) n and the number of radial nodes in the energy density nr = n/2 + (σ − 1)/4.The orbital angular momentum of the dominant mode is then ℓ − σ.The dashed black line represents the optimal radius (see eq. ( 56)), in this case r0 = 0.18 λ0, this value is also used in figs.3 and 5.
more than one dispersion relationship (solid, dashed, and dash-dotted lines).We will label these additional modes with a non-negative integer n, which labels the modes from left to right.The solution for ℓ = 1, n = 0 corresponds to a fundamental circularly polarised eigenmode (CP), whose fields can be seen in the first row of fig. 3. It turns out that the spin-weighted components F s (either √ εE s or √ µH s eq. ( 35)), although not being individually physical, are eigenfunctions of definite longitudinal momenta operators with eigenvalues p z = ℏk z , J z = ℏℓ, S z = ℏs and L z = J z − S z = ℏ(ℓ − s).The momentum operator acts on these eigenfunctions as where prime represents the derivative of the function with respect to its argument.Expression in parentheses is a complex wave vector for each component where κZ ′ ℓ−s /Z ℓ−s is real both inside and outside the fibre.Notice that the simple longitudinal and azimuthal components of the momentum come from the fact that F s are eigenfunctions of p z and L z .In particular, the longitudinal orbital angular momentum operator gives which means that ẑ • (r × k s ) = ℓ − s, which explains the azimuthal component of eq. ( 38). 5  The orthonormality of the spin basis ês ensures that the scalar product of any two vectors F and G is which means that an eigenmode F (either given by n and ℓ can be thought of as a linear combination of independent eigenfunctions F s ês of definite spin Equation ( 40) ensures that the energy densities of any eigenmode can be written as a sum where W s A are energy densities calculated for the individual vector components F s ês if all other components were zero.To calculate these, one can define two real and one complex coefficient quadratic in scalar A and pseudoscalar B such that they have the same symmetries as W 0 , W 3 and W c : One can then write the energy densities [eq.( 19)] associated with each spin-weighted component.The transverse spin (s = 0) component has 5 Recall that in cylindrical coordinates r = z ẑ + ρ ρ.
while the remaining components will be mixed where k z changes sign under a parity transformation.For a circular polarisation, from the symmetries, one would expect that w − and ℑw c will vanish, but that is not the case for a dielectric fibre due to the fact that there is electromagnetic asymmetry A ̸ = ±iB.This asymmetry comes from the refractive index being due only to the permittivity ε 1 .If instead we have a fibre made of material that has ε 1 = µ 1 = n 2 , then w − and ℑw c will vanish.
By decomposing the mode into spin-weighted functions, we may write the energy density gradients identically to a sum of three evanescent waves with ∇W s A = −2k ′′ s W s A , each having wavevector given by eq. ( 38) and energy densities given by eqs.( 43) and ( 44): This is valid for all energy density quantities A = 0, 1, 2, 3 defined in eq.(19).Remarkably, the same holds for canonical momenta, Notice that if we introduced factors of ℏ in both the numerator and denominator in the above sums, then the canonical momenta (gradients) would look like a weighted average of real (imaginary) momentum carried by each spin-weighted function where weights n s A = W s A /(ℏω) for A = 0 can be intuitively thought of as the number density of photons with the momentum given by eq.(37).For A = 3, it would be the number density of these photons that break parity, duality for A = 1, or both and time reversal for A = 2. Notice that for A ̸ = 0, these number densities can be positive or negative, representing whether the photons are more right/left-handed, electric/magnetic, etc.Using the same logic, eq. ( 42) is a weighted average of the energies ℏω carried by each spin component.
The interaction force for a pure fibre eigenmode (later we study mode combinations) will then be a sum of forces identical to three evanescent waves eq. ( 29) with wavevector given by eq. ( 38) for each of the spin components Unfortunately, this is where the similarities with the evanescent waves end.If a particle is big enough (with respect to the wavelength), it will also feel the recoil force, which depends on the fluxes of power and helicity.Those are related to the spin operator, but we can only repeat the trick we used in eq. ( 37) for the z component.Therefore, the spin densities have to be obtained from the definition eq. ( 21).The total spin density will be where we can see that the longitudinal component is again just the weighted average of the eigenvalues of each spin component.Notice that there is no longitudinal contribution from the eigenfunction with s = 0. Instead, this function will lead to a transverse spin, which depends only on W 0 A and k 0 .Notice that when we talk about longitudinal and transverse spin in this scenario, we mean with respect to the z-direction along the fibre.The spin quantity related to the real Poynting vector ωS 3 = ℜΠ/c will be very similar: and so will be the spin asymmetry vector One can see that all of these quantities have only longitudinal and azimuthal components; the only exception to this will be the spin quantity related to the imaginary Poynting vector which also has a radial contribution.The analytical form of the total force F = F int + F rec on a particle near the  fibre is then, in cylindrical basis: Note that the chiral force (assuming reciprocal particles α 2 = 0) are only the components which feature ℜ(α 3 ), ℑ(α 3 ), β 0 and β 1 , i.e., For a visual representation of the different force terms, the total force can still be expressed using the concept of a basis F = i λ i V i as introduced in eq. ( 25).It is then very instructive to plot the twelve different basis vector fields V i as given in table I.These are the twelve plots shown in the bottom rows of fig.3, and they are entirely determined by the fields of the circularly polarised mode (ℓ = 1, n = 0) while being independent of the particle.To supplement this figure, we also include a Python script [48] with interactive Jupyter notebooks capable of reproducing these twelve plots for any material parameters and for any mode.
Examining these basis vector fields V i , one can deduce that small particles (with negligible recoil, and therefore ignoring the last row) near the fibre will experience a total energy gradient force ∇W 0 attracting all particles towards the surface of the fibre, added to a helicity gradient force ∇W 3 which attracts one enantiomer towards the fibre, while repelling the other one away from it, exactly as was the case of an evanescent wave described earlier.This is due to the circular polarisation of the mode having a large helicity density in the fibre, which radially decays away, creating a gradient in the helicity.Added to these gradient forces are possible longitudinal (z-directed) pressure forces, as shown in p 0 and p 3 in fig. 3 which are relevant for particles with ℑ(α e +α m ) ̸ = 0 and ℑα c ̸ = 0, respectively.This can be easily deduced by looking at the coefficients λ i that multiply the relevant basis functions in table I.These gradient and pressure forces are completely equivalent to the in-plane separation by the evanescent wave (exactly as in fig. 1 particles   b and c), with the only difference being the geometry of the surface to which the particles will stick.For example, a particle with radius 300 nm, relative permittivity ε p = 2.5, permeability µ p = 1 and chiral parameter |κ| = 0.5 (considered in Li et al. [49]) would feel a chiral force of 205 fN mW −1 near a silicon nitride fibre that guides the fundamental circular mode with wavelength of 1310 nm and radius 240 nm.In this case, the chiral forces would be strong enough to separate enantiomers in order of seconds using just mW of power (taking the same assumptions as in Martinez-Romeu et al. [50]).Therefore, this would be a viable configuration for separating small enantiomers.The provided code [48] can calculate and plot the force fields for these and any other particle parameters.For bigger particles, the recoil forces become stronger.Looking at the last row in fig.3, corresponding to the recoil forces, we can see there is a strong azimuthal component in both the achiral recoil force kωS 3 = kℜΠ/c and the chiral spin recoil force kωS 0 = kω(S e + S m ) which would suggest a possible enantiomer sorting mechanism in this circular mode.In practice, however, the particle-dependent polarizabilities (λ i ) multiplying the achiral basis are much stronger than the chiral one, and so the achiral response dominates, making this circular mode not well suited for recoil force sorting.To achieve strong chiral recoil forces that are not overtaken by achiral forces, one must rely on mode combinations to synthesise linearly polarised modes, as shown later.
The cylindrical fibre can support linear combinations of eigenmodes E ℓ,n eq. ( 35) with complex coefficients a ℓ,n and similarly for the magnetic field.Modes with different ℓ and n, unlike the mode components with different s, are not linear in the quadratic observables.Therefore, energy densities will not be simply the sum of individual energy densities.
Instead, there will also be interference cross terms.For example, take the electric energy density where the summand is guaranteed to be real only if ℓ ′ = ℓ and n ′ = n, and in that case it would be the electric energy density of that particular mode W ℓ,n e .Other terms will be generally complex, and we only get a real energy density once we perform the sum.Let us, however, label these, in general complex, interference terms W ℓ,ℓ ′ ,n,n ′ ,s A , so that, for instance, if A = e (electric energy) we have for the fundamental circular polarisation modes with ℓ = 1.Solid lines represent chiral forces, while dashed lines represent achiral forces.Note that these curves are invariant in these units under a change of wavelength or power.The black dashed line represents the optimal radius of the fibre which maximises the chiral gradient force (for CP) and the spin recoil (for LP) (see eq. ( 56)).
where s again labels the component in the spin basis.Similarly, we can define them for A = 0, 1, 2, 3. Canonical momenta and gradients of energy densities will be, in this most general case, as follows: where the complex weights in these weighted averages can be thought of as modified number densities of photons with complex momentum ℏk ℓ,n,s of the species A, but taking into account the interference with all the other eigenmodes.This quantity is, in general, complex and its phase will redistribute (rotate) between how much of k ′ ℓ,n,s and k ′′ ℓ,n,s contributes to the gradient and pressure forces.In theory, one can use these interferences to design gradient and pressure forces of arbitrary directions and amplitudes by controlling coefficients a ℓ,n in eq. ( 52).However, in practice, we don't have such fine experimental control over these coefficients.
An example of an experimentally easy-to-obtain mode superposition is that of a linear polarisation.Two orthogonal linear polarisation modes are usually called the quasi-transverse electric (TE) and quasi-transverse magnetic (TM) modes.These modes are obtained in general as the superpositions The fields of the linearly polarised TE mode with ℓ = 1, n = 0 can be seen in the top row of fig. 5.In this case, the two modes TE and TM are degenerate in the sense that their fields are identical, just rotated by an angle π/2, and the same is true for all the force densities.
To visualise the possible forces on this TE mode, we once again plot the twelve different basis vector fields V i from table I, shown in fig. 5.In this linear polarisation, the gradient and pressure forces are clearly dominated by the non-chiral total energy gradient ∇W 0 and momentum p 0 , hence this mode is not well suited to smaller molecules in which the gradient and pressure forces dominate.Instead, for larger particles on which the recoil forces dominate, this linearly polarised mode shows more promise.Indeed, looking at the recoil terms (last row in fig.5) one can see an achiral Poynting recoil kωS 3 = kℜΠ/c that is mostly longitudinal (z-directed), while there is a strong in-plane azimuthal spin angular momentum recoil force kωS 0 = kω(S e + S m ) pointing around the fibre.This means that opposite enantiomers will be pushed in opposite directions around the fibre, resulting in a potential lateral sorting mechanism.This is caused by the transverse spin, and is exactly equivalent to the lateral sorting in linearly polarised evanescent waves from the previous section.Interestingly, this force competes with an achiral force proportional to the flux of reactive power −kωS 2 = kℑΠ/c which depends on the coefficient β 2 ∼ ℑ(α * e α m ) and which is present for any particle that is breaking dual symmetry (i.e. a dominant electric or magnetic response).For the same particle as before (a = 300 nm, ε p = 2.5, µ p = 1, and |κ| = 0.5) and linear polarisation with 1310 nm in a silicon nitride fibre with r 0 = 240 nm the maximal chiral force on the particle would be 218 fN mW −1 , which would again make it possible to separate these particles in order of seconds.Finally, it is interesting to study the effect that the fibre radius has on the different force terms.To this end, fig. 4 plots the maximum magnitude of each of the twelve force term basis V i in both the circularly polarised and linearly polarised fibre modes as a function of normalised radius.One can see that there is clearly an optimal fibre radius-to-wavelength ratio.This optimum is the sweet spot between the radius being too large (hence mode not interacting with the particles outside) and too small (hence mode fields very spread outside, with low gradients).In the limit where the fibre is in a medium for which n 2 = 0 this radius coincides with the point at which the mode with ℓ = 0 starts.The condition for a single-mode fibre operation is k 0 r 0 n 2 1 − n 2 2 < j 0,1 [51], where j 0,1 is the first root of a Bessel function J 0 with the approximate value of 2.40483.While the radius at which the fibre can support additional modes depends on n 2 , the optimal radius does not, leading to a simple expression for the optimal radius to wavelength ratio for any value of n 2 : which only depends on the wavelength and the refractive index of the material of the fibre.

VI. CONCLUSIONS
The total optical force acting on a general dipolar particle is a relatively large analytical expression that can be written in many alternative but equivalent ways.We provided a clear view of the different alternatives used in the literature and introduced our own alternative that classifies forces based on the symmetries broken by the particle and the fields.The chiral forces that separate enantiomers rely on particles and fields that break parity-but other forces can separate, for instance, electric from magnetic particles, which break duality symmetry.The concept of a force basis, relying on twelve force fields, each with its own particle-dependent coefficient, was introduced-and their manifestation in evanescent waves and nanofibre modes was exemplified.We also developed very concise analytical expressions for the optical forces in such modes-including the case of arbitrary nanofibre modes and their combination.These analytical expressions suggest that the separation of bigger enantiomers, where the recoil force dominates, should rely on transverse-spin-based lateral forces requiring linearly polarised modes, while separation of smaller (moleculesized) enantiomers should rely on gradient and pressure forces-so the use of circularly polarised modes provides a helicity gradient to attract or repel opposite enantiomers towards or away from the nanofibre.An insight that seems to be also true in the case of a rectangular waveguide [50].An optimal radius was found for maximising forces in nanofibres, which depends only on their material and wavelength.Finally, we provide an interactive Python script [48] that can reproduce all figs.3 to 5 for an arbitrary choice of parameters and modes.We hope that our theoretical work provides guidance and clarity for the design of future experimental attempts at optical separation of chiral enantiomers near waveguides, which is of enormous practical importance in the pharmaceutical domain.

(A3)
Appendix B: Fluxes of energy and helicity densities Interestingly, the four spin-like quantities, S A , in eq. ( 21) can be seen as the fluxes of the energy densities W A in eq. ( 19).Starting with the real Poynting vector ℜΠ = ωcS 3 , the Poynting theorem (continuity equation for active power) relates the time derivative of the total active energy density W 0 to ℜΠ which can be interpreted as the flow of active power [20].Somewhat less known is the recently proposed complex Poynting theorem [20,52] whose real part gives the usual time-averaged Poynting theorem for a time-harmonic field, where the term with W 0 vanishes since the average energy density does not change with time.The imaginary part relates the reactive or stored energy density W 1 and an alternating flow of reactive power ℑΠ = −ωcS 2 , which in our notation reads: where the term on the right-hand side is the reactive power.Similarly, it is a well-known fact that the conservation and flow of integrated optical helicity (the continuity equation) relates helicity density S = W 3 /ω to the spin (helicity flux) density S 0 [53,54].For time-averaged quantities, the W 3 will be again missing like in the case of the real Poynting theorem because the time-averaged helicity will be conserved and we will get notice the appearance of pseudoscalar quantity in units of power on the right-hand side, which seems to be a chiral equivalent of active power.A similar relationship to eq. (B2) then exists between the remaining two quantities, suggesting that S 1 is the flow of W 2 One can show that chiral energy density W 3 , chiral momentum p 3 and chiral spin angular momentum den-sity S 3 are just differences between energies, momenta and spins carried by the right-and left-handed fields.In order to do that one can write the angular spectrum decomposition of the electric field to separate it into positive and negative helicity components

(C2)
Using the definition of W 3 one can show that where W ± is the energy density carried by the right-or left-handed component of the field.The same can be done for the canonical momentum density and for the spin angular momentum density that is related to the Poynting vector ℜα e W e +ℜα m W m +ℜα c ℜW c +ℜα t ℑW c ) gradient force +2ω(ℑα e p e +ℑα m p m +ℑα c ℜp c +ℑα t ℑp c ) radiation pressure force −(σ rec ℜΠ +σ im ℑΠ)/c − ω(γ e rec S e +γ m rec S m ) dipole recoil force .

FIG. 1 .
FIG. 1. Separation of enantiomers using an evanescent wave.An illustrative figure, forces not to scale.For the larger particle (a) we only show the lateral recoil force, while for the lossless particle (b) and the lossy particle (c) we show the interaction forces.

Pressure forces 2ωp 0 FIG. 3 .
FIG.3.Circular force density basis for optical force around a dielectric fibre guiding total angular momentum eigenmode with ℓ = 1.Arrow maps represent transverse components, while colour maps represent longitudinal components of vector quantities plotted in a cross-sectional view of the fibre.This force basis plot is independent of wavelength.It depends on the material of the fibre SiN (ε1 = 4.3), the surrounding medium water (ε2 = 1.7) and the radius chosen to be r0 = 0.18 λ0.

FundamentalFieldsFIG. 5 .
FIG.5.Linear force density basis for optical force around a dielectric fibre.Arrow maps represent transverse components, while colour maps represent longitudinal components of vector quantities plotted in a cross-sectional view of the fibre.This force basis plot is independent of wavelength.It depends on the material of the fibre SiN (ε1 = 4.3), the surrounding medium water (ε2 = 1.7) and the radius chosen to be r0 = 0.18 λ0.