Using the Belinfante momentum to retrieve the polarization state of light inside waveguides

Current day high speed optical communication systems employ photonic circuits using platforms such as silicon photonics. In these systems, the polarization state of light drifts due to effects such as polarization mode dispersion and nonlinear phenomena generated by photonic circuit building blocks. As the complexity, the number, and the variety of these building blocks grows, the demand increases for an in-situ polarization determination strategy. Here, we show that the transfer of the Belinfante momentum to particles in the evanescent field of waveguides depends in a non-trivial way on the polarization state of light within that waveguide. Surprisingly, we find that the maxima and minima of the lateral force are not produced with circularly polarized light, corresponding to the north and south poles of the Poincaré sphere. Instead, the maxima are shifted along the great circle of the sphere due to the phase differences between the scattered TE and TM components of light. This effect allows for an unambiguous reconstruction of the local polarization state of light inside a waveguide. Importantly, this technique depends on interaction with only the evanescent tails of the fields, allowing for a minimally invasive method to probe the polarization within a photonic chip.


Optical forces on Rayleigh particles in an evanescent field
The force acting on a Rayleigh particle in an evanescent field can be calculated analytically. The time-avaraged Lorentz force acting on such a spherical, dipolar particle with radius a consists of three terms 1 : Here, E 0 is the amplitude of the electric field with wave vector k, and α is the polarizability of the particle. For small particles (ka 1) this polarizability can be approximated by α = α 0 1 + 2/3ik 3 |α 0 | 2 with α 0 = a 3 (ε r − 1)/(ε r + 2), where ε r = ε p /ε b is the relative permittivity, ε p the permittivity of the particle, and ε b the permittivity of the background medium. The first term in Eq. (S1) corresponds to the gradient force, whereas the second and the third term are scattering forces. The second term corresponds to the radiation pressure whereas the last term is proportional to the spin density of the field. We consider an evanescent field generated by an internally reflected TM-polarized plane wave (see Fig. 1 in the main paper). In this case, the first term pulls the particle towards the surface, the second one pushes the particle along the surface in the direction of propagation, and the final term depends on the curl of the helicity of light but remains negligible in the dipole approximation. In the case of linearly polarized light, the final term disappears altogether. The ratio of the forces perpendicular to (z direction) and longitudinal along the surface (x direction) is thus given by: where . |F z /F x | is thus entirely determined by the ratio of the real and the imaginary parts of the particle's polarizability α and the evanescent field's wave vector: 1 F x F z x z y F y Figure S1. A visualization of the forces acting on a particle in the evanescent field of a total interally reflected beam. In the Rayleigh approximation, the gradient force pulls the particle in the z direction and the scattering force pushes the particle along the x direction. In the Mie regime, this picture become more complex and scattering forces appear in the three different directions. In a typical measurement of these forces, the particle is held at a fixed position using an optical trap (focused Gaussian beam from the objective) The values of q and k z are fixed by the angle of incidence. Using Draine's radiation correction for the polarizability of a small particle: α = α 0 1 − 2 3 ik 3 α 0 −1 ≈ α 0 1 + 2 3 ik 3 α * 0 , one finds that, for lossless particles, the ratio of the real and the imaginary parts of the polarizability equals The ratio of the force in the z and the x direction grows to infinity for decreasing particles radius a, as the static polarizability α 0 is proportional to a 3 .

Optical forces on Mie particles in an evanescent field
When the size of a dielectric particle is comparable to the wavelength of the incident light, the Rayleigh approximation is not valid anymore. In this case, one needs to solve the full electromagnetic problem, calculating the scattered fields imposing Maxwell's boundary equations at the surface of the particle. For homogeneous, isotropic, spherical particles, this problem was solved by Gustaaf Mie in 1908 and is now known as the generalized Lorenz-Mie theory 2 . The formalism starts by decomposing the incident fields-with subscript (i)-as a sum of spherical harmonics: Here, a is the radius of the particle, n b is the refractive index of the background, k 0 is the free space wave vector (k 0 = ω/c), E 0 and H 0 are the amplitudes of the electric and magnetic fields, Y l,m (θ , φ ) are the spherical harmonics, and ψ l (r) = r j l (r) are the Riccati-Bessel functions of the first kind. These functions are defined as ψ l (r) = r j l (r), with j l (x) the spherical Bessel function of the first kind. It is sufficient to consider only the radial parts of the electromagnetic field to carry out this decomposition, as the other components are directly related to the radial ones through Maxwell's equations. Using the orthogonality of the spherical harmonics, the coefficients A l,m and B l,m can be retrieved by evaluation of the following integrals at the surface of the particle:
Next, one considers a similar decomposition of the scattered fields-with subscript (s). However, in order to impose the correct boundary conditions at infinity, one traditionally uses modified Riccati-Bessel functions ξ (r) = rh are the spherical Hankel functions of the first kind: In general, one would have to write a similar decomposition for the fields inside the particle and subsequently impose the continuity of the tangential parts of the fields at the boundary of the particle. Following this procedure, Gustav Mie derived a closed form expression for the scattered fields as a function of the incident fields, or equivalently, the scattered coefficients a l,m and b l,m as a function of the incident coefficients A l,m and B l,m : where k b = n b k 0 and k p = n p k 0 .
At this point, the full electromagnetic problem is solved and the time-avaraged force acting on the particle can, in principle, be calculated by integrating the flux of the Maxwell stress tensor over a surface that surrounds the particle: where Unfortunately, this procedure is rather cumbersome as the integral converges slowly for larger particles and is very sensitive to discretization error. In order to avoid this last integration step, Barton et al. derived a formula that allows for the calculation of the force using an algebraic combination of the Mie scattering coefficients 3 . Using this algorithm, we calculated the time-averaged force as (S15) The infinite series in the Mie scattering algorithm was truncated at a value l for which the relative magnitude of the different force contributions dropped below 10 −7 .

The location of the maximum of the lateral force
It is interesting to take a closer look at the location on the Poincaré sphere where the maximal/minimal lateral force appears. In Fig. 2F of the main paper, we plot the angle φ of this polarization as a function of particle size. This angle uniquely defines the location because the polarization that generates the largest lateral force always lies on the great circle defined by ψ = π/4. Physically, this can be understood because of the fact that the lateral force is proportional to the product of the amplitudes of the TE-and the TM-polarizations. Recalling the definition of ψ and φ , we find that the force is proportional to sin(2ψ), which is maximal at ψ = π/4. It is also intriguing to note that the different transmission amplitudes for TE-and TM-polarizations at the interface only affect the φ component of the polarization maximum/minimum. Indeed, the electric field above the interface can be written as: where t TM and t TE are the transmission coefficients at the interface. It is interesting to note here that the amplitudes of the E T E and the E T M modes above the surface are determined by the parameter ψ, on the one hand, and the Fresnel transmission coefficients, on the other hand. However, the product of the TE and TM components is still maximized at ψ = π/4. The phase difference between the TE and TM transmission coefficients shifts the location of the maximum on the Poincaré sphere in the φ direction. This is visualized in Fig. S2, where we compare the lateral forces that act on a silicon particle (a = 100nm) between (A) a setup where the transmission coefficients are identical and (B) a setup where the transmission coefficients are given by the normal Fresnel values. The location of the maximal lateral force (red dot in both plots) is given by (π = π/4, φ = 1.52π) in (A) and by (π = π/4, φ = 1.30π) in (B). The shift of this location in the φ direction is in agreement wit the phase difference between t TE and t TM at θ = π/4, where arg(t TE ) − arg(t TE ) = 0.22π.

The evanescent field in Cartesian coordinates
In order to calculate the force components acting on a dielectric particle in the Mie regime, it is convenient to write the evanescent field in Cartesian coordinates. We can rewrite Eq. (S19) by expressing E TE and E TE as a function of 1 x , 1 y , and 1 z , using the expressions The full expression of the electric field above the interface is thus given by

5/8
In this section, we discuss the equilibrium position of a particle above a (birefringent) waveguide in detail. The geometry of the waveguide, shown in Fig. 4 of the main paper, consists of a silicon core (400 nm wide, 300 nm high) and a silica cladding. The lateral equilibrium position (y) of a particle above the waveguide is essentially determined by two different forces: (1), a gradient force, attracting the particle to the position of largest intensity, and (2), the helicity dependent lateral force (caused by the Belinfante momentum). Let us first consider the gradient force. This force essentially appears because of the "snaking" intensity profile of the field above the waveguide. This is a result of the different propagation constants of the different modes inside the waveguide, in combination with the fact that the field intensity above the waveguide shifts when you add the different modes with with different weights (ψ) and different relative phase (φ ). To obtain this profile, we calculate the field amplitude at the interface of the waveguide using full wave numerical simulations (COMSOL Multiphysics). In Fig. S3(A-C), we show a few examples of these profiles for different values of ψ and φ . It is clear that the maximum of the field profile shifts as a function of the polarization state inside the waveguide. This equilibrium location (y 0 ) was extracted for each different polarization state of light inside the waveguide and is visualized in Fig. S3(G).
The second contribution is the helicity-dependent lateral force. In contrast to the gradient force, this force is only non negligible for larger particles, with radii in the Mie regime. This extra lateral force acting on 300 nm radius silicon particles was calculated using the Mie scattering algorithm, where the incident wave is now a superposition of the two evanescent waves of the different modes inside the waveguide. This result is shown in Fig. S4(A). It is interesting to note that there is an asymmetry in this density plot, displaying larger positive forces than negative forces. This can be understood when we look the intensity of the mode profiles above the waveguide as a function of the polarization state, as shown in Fig. S3(H). The asymmetry of the lateral force is thus entirely determined by the differences in intensity of the light above the waveguide. In order to find the overal displacemnt of the 300 nm radius particle, we still need to know how large the spring constant is of the gradient force attracting the particle towards the largest intensity. This spring constant is also a function of the polarization of light. Using the gradients the field cross sections (Fig S3(D-F)), we calculated this spring constant as shown in Fig. S3(I). Not surprisingy, the polarization regions with largest field amplitudes in Fig. S3(H) are also the polarization regions with the largest spring constants in Fig. S3(I). As a result, we notice that the extra displacement ∆y = F lat /k y , shown in Fig. S4(B) is again more symmetric. Finally, in Fig. S4(C), we plot the overlap of the contour lines of y 0 and ∆ y , demonstrating that the points of intersection allow for the retrieval of the polarization state inside the waveguide. In general, the two roughly elliptical contour lines will have two intersections, from which the correct intersection can be chosen using the asymmetry of Fig. S3(H) in combination with a calibration of the scattering intensity of light from the particle.
We  helicity-dependent lateral force acting on a 300 nm radius silicon particle on top of a birefringent waveguide. (B) The displacement of the particle, found by balancing the lateral force with the spring force of the intensity potential pulling the particle towards the largest intensity: ∆y = F lat /k y . (C) The overlap of several contour lines of y 0 and ∆y. Because of the relative displacement of these contour lines in the φ direction, we generally find only two intersections per set of two measurements. The correct intersection can be retrieved using a calibration of the scattering intensity.