Size-dependent position of a single aerosol droplet in a Bessel beam trap

An ideal linearly polarized Bessel beam can be described by its amplitude, order, semi-apex angle θ₀ and wavenumber k. For a zeroth-order Bessel beam, the intensity profile in the transverse plane will consist of a bright, circular core surrounded by concentric rings of approximately equal power [30, 31]. The Bessel beam itself can be thought of as a superposition of plane waves travelling at an angle θ₀ to the longitudinal axis [32]. These plane waves possess equal amplitude and are evenly distributed around all azimuthal angles and, as such, are often viewed as propagating on a cone [30, 33, 34]. The longitudinal and transverse components of the wavevector of the Bessel beam are then defined as k_z = k cos θ₀ and k_⊥ = k sin θ₀, respectively.

When a Bessel beam is generated using a Gaussian beam of waist w_c passing through an axicon with an internal angle γ and refractive index of n_ax, the semi-apex angle will be defined as [35]

$$\begin {equation} \theta _0 = \mathrm {arcsin} \, \left (\frac {n_{\rm ax}}{n_0} \mathrm {sin} \, \gamma \right ) -\gamma \label {apex} \end {equation} \tag{ 1 }$$

where n₀ is the refractive index of the medium. For small values of θ₀, equation (1) will be θ₀ ≈ (n_ax/n₀ − 1)γ. A beam generated in this way will not be diffraction free and will have a maximum propagation distance z_max defined as [36]

$$\begin {equation}z_{\max }=\frac {w_c}{\mathrm {tan} \, \theta _0}. \end {equation} \tag{ 2 }$$

With these parameters, the expression for the electric field of a zeroth-order Bessel beam that is linearly polarized along the x axis can be written as [23, 35–37]:

$$\begin {eqnarray} \textbf {E(r)}&=&E_0 \sqrt {\frac {2z}{z_{\max }}}{\rm e} ^{{\rm i} k_z z-z^2/z_{\max }^2 +(1-\pi {\rm i} )/4}\nonumber \\ &&\times \,\left \{ \begin {array}{l} \lbrack J_0(k_\perp \rho ) + \beta ^2 J_2(k_\perp \rho ) \mathrm {cos} \, 2 \phi \rbrack \hat {\textbf {x}} \\ \beta ^2 J_2(k_\perp \rho ) \mathrm {sin} \, 2 \phi \, \hat {\textbf {y}} \\ -{\rm i} 2 J_1(k_\perp \rho ) \mathrm {cos} \, \phi \, \hat {\textbf {z}} \end {array} \right \} \label {Bessel-fields} \end {eqnarray} \tag{ 3 }$$

where β = sin θ₀/(1 + cos θ₀), $\rho =\sqrt {x^2+y^2}$, and ϕ = arctan (y/x). The unit vectors $\hat {\textbf {x}}$, $\hat {\textbf {y}}$, and $\hat {\textbf {z}}$ are in the Cartesian coordinate system. A time dependence of e^−iωt has been omitted. The factor E₀ in (3) is chosen so that $|\textbf{E}(0,0,z_{\max }/2)|^2=E_0^2$. The relationship between E₀ and the power of the Gaussian beam P_T is

$$\begin {equation}E_0^2=4{\rm e}^{-1/2}\frac {P_{\rm T} k_\perp }{w_c c \epsilon _0 n_0}. \label {E0} \end {equation} \tag{ 4 }$$

To obtain beam widths that are suitable for experimental use, the Bessel beam needs to be reduced with a telescope. For a telescope consisting of two lenses with focal lengths f₁ and f₂ the magnification will be M = f₂/f₁. The transformation of the Bessel beam parameters are then described by [38]

$$\begin {eqnarray} \mathrm {sin} \, \theta _0^\prime =\frac {\mathrm {sin} \, \theta _0}{M}, \end {eqnarray} \tag{ 5 }$$

$$\begin {eqnarray}w_{\rm c}^\prime =Mw_{\rm c}, \end {eqnarray} \tag{ 6 }$$

where $\theta _0^\prime$ and $w_c^\prime$ are the semi-apex angle and beam waist, respectively, after passing through the telescope.

For the setup described in section 2, the calculated values of θ₀ and z_max are 0.780^° and 550.9 mm, respectively. The values of w_c and n_ax are known to be 7.50 mm and 1.52, respectively. After the Bessel beam is passed through a telescope of magnification M = 1/9 the values of θ₀, w_c, and z_max become 7.04^°, 0.833 mm, and 6.75 mm, respectively. For a beam with an input power of 768 mW and a wavelength of 532.0 nm, the calculated value of E₀ after the telescope is 1.10 × 10⁶ V  m⁻¹.

Optical forces were calculated using the equations derived by Barton et al [39]. All beams considered here are propagating in the + z direction so the time averaged longitudinal force 〈F_z〉 is

$$\begin {eqnarray} \langle F_z \rangle & = & -\frac {\alpha ^2 a^2 \epsilon _0 E_0^2}{2} \sum _{l=1}^{\infty } \sum _{m=-l}^{l} \mathrm {Im} \left ( {\vphantom {\sqrt {\frac {(l-m+1)(l+m+1)}{(2l+3)(2l+1)}}}} l(l+2) \right .\nonumber \\ & & \times \, \sqrt {\frac {(l-m+1)(l+m+1)}{(2l+3)(2l+1)}} \nonumber \\ & & \times \, (2 n_0^2 a_{l+1,m}a_{lm}^* + n_0^2 a_{l+1,m}A_{lm}^* +\, n_0^2 A_{l+1,m}a_{lm}^* \nonumber \\ & & + 2 b_{l+1,m}b_{lm}^* + b_{l+1,m}B_{lm}^* + B_{l+1,m}b_{lm}^*)\nonumber \\ & & \left . {\vphantom {\sqrt {\frac {(l-m+1)(l+m+1)}{(2l+3)(2l+1)}}}} +\, n_0 m (2 a_{lm} b_{lm}^* + a_{lm} B_{lm}^*+ A_{lm} b_{lm}^*) \right ) \label {force-z} \end {eqnarray} \tag{ 7 }$$

and the time averaged forces in the transverse plane 〈F_x〉 and 〈F_y〉 are

$$\begin {eqnarray} \langle F_x \rangle + {\rm i} \langle F_y \rangle & = & \frac {{\rm i} \alpha ^2 a^2 \epsilon _0 E_0^2}{4} \sum _{l=1}^{\infty } \sum _{m=-l}^{l} \left ( {\vphantom {\sqrt {\frac {(l+m+2)(l+m+1)}{(2l+3)(2l+1)}}}} l(l+2) \right . \nonumber \\ & & \times \, \sqrt {\frac {(l+m+2)(l+m+1)}{(2l+3)(2l+1)}} \nonumber \\ & & \times \, (2 n_0^2 a_{lm}a_{l+1,m+1}^* + n_0^2 a_{lm}A_{l+1,m+1}^*\nonumber \\ & & +\, n_0^2 A_{lm}a_{l+1,m+1}^*+ 2 b_{lm}b_{l+1,m+1}^* \nonumber \\ & & +\, b_{lm}B_{l+1,m+1}^*+ B_{lm}b_{l+1,m+1}^*) \nonumber \\ & & + \, l(l+2) \sqrt {\frac {(l-m+1)(l-m+2)}{(2l+3)(2l+1)}} \nonumber \\ & & \times \, (2 n_0^2 a_{l+1,m-1}a_{lm}^* + n_0^2 a_{l+1,m-1}A_{lm}^* \nonumber \\ & & +\, n_0^2 A_{l+1,m-1}a_{lm}^* + 2b_{l+1,m-1}b_{lm}^* \nonumber \\ & & +\, b_{l+1,m-1}B_{lm}^* + B_{l+1,m-1}b_{lm}^*)\nonumber \\ & & -\,n_0\sqrt {(l+m+1)(l-m)} (-2a_{lm}b_{l,m+1}^* \nonumber \\ & & + \, 2b_{lm}a_{l,m+1}^* -a_{lm}B_{l,m+1}^* + b_{lm}A_{l,m+1}^* \nonumber \\ & & \left . {\vphantom {\sqrt {\frac {(l+m+2)(l+m+1)}{(2l+3)(2l+1)}}}} +\, B_{lm}a_{l,m+1}^*- A_{lm}b_{l,m+1}^*) \right ), \label {force-xy} \end {eqnarray} \tag{ 8 }$$

where the size parameter is α = ka and a is the radius of the sphere. In practical calculations the series in l cannot be evaluated up to l = ∞ and must eventually be truncated. Here, the series was truncated at the integer closest to α + 4α^1/3 + 2 [28]. The coefficients a_lm, b_lm, A_lm, and B_lm are defined as

$$\begin {eqnarray} a_{lm}=\frac {\psi _l^\prime (\bar {n} \alpha ) \psi _l (\alpha )-\bar {n} \psi _l(\bar {n} \alpha ) \psi _l^\prime (\alpha )}{\bar {n} \psi _l (\bar {n} \alpha ) \xi _l^{(1)\prime } (\alpha )-\psi _l^\prime (\bar {n} \alpha ) \xi _l^{(1)} (\alpha )} A_{lm},\label {alm} \end {eqnarray} \tag{ 9 }$$

$$\begin {eqnarray} b_{lm}=\frac {\bar {n} \psi _l^\prime (\bar {n} \alpha ) \psi _l (\alpha )-\psi _l(\bar {n} \alpha ) \psi _l^\prime (\alpha )} { \psi _l(\bar {n} \alpha ) \xi _l^{(1)\prime } (\alpha )- \bar {n} \psi _l^\prime (\bar {n} \alpha ) \xi _l^{(1)} (\alpha )} B_{lm},\label {blm} \end {eqnarray} \tag{ 10 }$$

$$\begin {eqnarray} A_{lm}=\frac {1}{l(l+1)\psi _l(\alpha )} \int _0^{2 \pi } \int _0^{ \pi } (E_r(a,\theta ,\phi )/E_0) \nonumber \\ \quad \times \, \mathrm {sin} \, \theta \, Y_{lm}^*(\theta , \phi ) \,{\rm d} \theta \,{\rm d} \phi ,\label {A} \end {eqnarray} \tag{ 11 }$$

$$\begin {eqnarray} B_{lm}=\frac {1}{l(l+1)\psi _l(\alpha )} \int _0^{2 \pi } \int _0^{ \pi } (H_r(a,\theta ,\phi )/H_0) \nonumber \\ \quad \times \, \mathrm {sin} \, \theta \, Y_{lm}^*(\theta , \phi )\, {\rm d} \theta \,{\rm d} \phi , \label {B} \end {eqnarray} \tag{ 12 }$$

where $\xi _l^{(1)}=\psi _l-{\rm i} \chi _l$, ψ_l and χ_l are the Ricatti–Bessel functions, $\bar {n}=n/n_0$ is the relative refractive index of the particle, Y_lm are the spherical harmonics, primes indicate differentiation with respect to the argument of the given function, H₀ = E₀/(μ₀c), and E_r and H_r are the radial components of the incident electric and magnetic fields, respectively, at the surface of the sphere. Equations (11) and (12) were evaluated numerically [40]. The program used to calculate the forces was tested and validated against previously published results for Gaussian [39] and Bessel beams [36, 41]. Additionally, calculations involving a plane wave were performed using Mie theory [28]. All calculations shown in this work were performed on spherical glycerol droplets (n = 1.4746) [29] in air (n₀ = 1).

For a spherical droplet in a gas stream the drag force F_d from Stokes' law is

$$\begin {equation}F_{\rm d}=6 \pi \mu a v, \label {drag} \end {equation} \tag{ 13 }$$

where μ is the dynamic viscosity of the medium and v is the undisturbed speed of the gas that would be found at the location of the droplet.

The Poiseuillian speed distribution of a free laminar jet is [42]

$$\begin {equation}v=\frac {Q^2}{2 \pi ^2 \eta z_i r_0^2 (1+(\frac {Q r}{4 \pi r_0 \eta z_i} )^2)}, \label {speed} \end {equation} \tag{ 14 }$$

where z_i is the distance from the inlet, Q is the volumetric flow rate through the inlet, η is the kinematic viscosity of the gas, r₀ is the radius of the inlet, and r is radial distance from the z axis (the jet centerline). In all calculations, the medium was assumed to be air at a temperature of 298 K and a pressure of 1 atm, yielding viscosities of μ = 1.8 × 10⁻⁵ Pa s and η = 1.5 × 10⁻⁵ m² s⁻¹ [43]. The values for Q and r₀ were given in section 2.

In the longitudinal direction, the net force F_t that a particle experiences will be

$$\begin {equation}F_{\rm t}=\langle F_z \rangle -F_{\rm d}. \label {netforce} \end {equation} \tag{ 15 }$$

Equilibrium positions will exist when F_t = 0. Figure 3 shows the calculated net force for a droplet with a = 2.0 μm in a Bessel beam. There are two equilibrium positions: one that is unstable (z₀ = 0.14 mm) and one that is stable (z₀ = 9.38 mm) with the negative slope of the force curve characteristic of a stable restoring force that acts to immobilise the particle at the zero force crossing. By calculating the latter position across a range of radii, the retention behaviour as determined by the variation in equilibrium position with droplet size can be quantified.

Zoom In Zoom Out Reset image size

While there is no drag force in the transverse direction, gravity will be present with the chosen experimental configuration of a horizontally propagating laser beam and gas flow. However, preliminary calculations and testing found the effect of gravity not to be significant and it will subsequently be ignored here. Therefore, in the transverse direction only optical forces need to be considered.

Figure 4(a) shows the transverse force that droplets of several different radii experience as they are displaced along the x axis (the polarization axis of the zeroth-order Bessel beam) while remaining at z₀ = z_max/2 and y₀ = 0. The corresponding potential energy curves are shown in figure 4(b). Washboard-like potentials [22, 36, 44] are observed in figure 4(b) where the barrier to escape a local energy minimum is always smaller when moving towards the center of the beam (x = 0) than it is moving away from the center. Therefore, in the transverse direction of the Bessel beam, droplets will be guided towards the center of the beam [22, 36, 44].

Zoom In Zoom Out Reset image size

Minima in the potential energy curves are located at positions that approximately correspond to intensity maxima in the rings of the Bessel beam, suggesting that droplets will preferentially adopt equilibrium positions in regions of maximum light intensity. However, not every intensity maximum in the Bessel beam has an associated minimum in the potential energy curve. For instance, in figure 4(b) the only equilibrium position for the a = 1.5 μm curve is at x₀ = 0 which coincides with the center of the core of the Bessel beam. If a droplet with a radius that is initially greater than 1.5 μm is trapped in the Bessel beam and its radius were to slowly decrease (e.g. through evaporation), then it will eventually be confined along the x axis at x₀ = 0 (as its radius approaches 1.5 μm). Once located at x₀ = 0, the droplet will unlikely be able to access any of the local energy minima seen in the curves for a = 0.5 and 1.0 μm due to the high energy barrier to transverse motion. Furthermore, the low energy barriers in the a = 2.0 μm curve mean that the droplet will probably be trapped at x₀ = 0 for radii well above 1.5 μm (for this curve a droplet approaching x₀ = 0 must only overcome barriers that are less than 10 k_BT in height).

Zoom In Zoom Out Reset image size

The intensity and Poynting vector of the Bessel beam have circular symmetry around the propagation axis of the beam (the z axis). For a droplet in this beam, however, the resulting transverse forces and potential energy surface will not be circularly symmetric. For example, the potential energy curves shown in figure 4(b) will not have the same energies if the displacements were along the y axis instead of the x axis. This lower symmetry originates from the linear polarization of the Bessel beam. The main consequence of this is that the potential energy curves for displacements in other direction in the transverse plane will have smaller barriers than those along the x axis. However, this change is not significant enough to affect the main conclusion that large droplets (a > 1.5 μm) will always tend to move onto the z axis as their radius decreases. Finally, even though the beam is linearly polarized, in the Rayleigh limit (a ≪ λ) the transverse forces calculated with the method described in section 3.2 will become approximately circularly symmetric around the z axis. This result is consistent with the well-known electrostatic expression that describes the gradient force in this limit [2, 45, 46].

Figure 5 compares the longitudinal (z direction) force that a spherical glycerol droplet experiences while in (a) a plane wave and (b) a zeroth-order Bessel beam. Note that drag forces are not considered in this section. In the Bessel beam, the droplet is located at x₀ = y₀ = 0 and z₀ = z_max/2. For this comparison, the amplitude E₀ of the plane wave is the same as that of the Bessel beam. For very small radii (a < 0.2 μm), the force that a droplet experiences is similar in both beams. In this size regime, the transverse intensity profile of the Bessel beam will change very little across ρ = −a to + a and even less across z = z_max/2 − a to z = z_max/2 + a. As an approximation, the intensity can therefore be set to its value at the droplet center. This treatment yields a profile identical to that of a plane wave (which has a constant intensity at every position) and explains the similarities between the two results in this size regime.

As a increases (0.2 μm < a < 1.1 μm) the calculated forces when in a Bessel beam and in a plane wave differ considerably, as shown in figure 5. The energy flux crossing droplets of identical radii will always be greater for the plane wave than for the Bessel beam due to the constant transverse intensity profile of the plane wave and its chosen value of E₀. Consequently, the difference between the two calculated forces becomes very large. Sharp peaks corresponding to morphology-dependent resonances or WGMs [47–58] also appear in this size regime. Despite the difference in the magnitude of the calculated forces, these peaks appear at identical radii. The type of beam used in the calculation has no effect on the radii at which WGMs occur.

The most noticeable difference between the two force curves in figure 5 occurs in the next size regime (1.1 μm  < a < 1.6 μm). The sharp peaks that are observed in the plane wave force curve are either absent or greatly diminished in the Bessel beam force curve. In this range of radii, the excitation efficiency of the WGMs by the Bessel beam appears to be very low. In a spherical droplet, WGMs are excited by rays that pass through a region just outside the radius of the sphere. Known as van de Hulst's localization principle [59], excitation occurs through coupling between the incident beam and the evanescent field of the WGM. This phenomenon has been studied extensively for Gaussian beams that are incident on spherical particles [60–68]. Coupling between the incident beam and WGMs will be optimal at a radial distances b_l from the center of the sphere (note, do not confuse the symbol b_l with the Mie coefficient b_lm). The distance b_l has been referred to as the impact parameter [67–69] and for mode number l is defined as [59]

$$\begin {equation}b_l=\frac {(l+1/2)}{k}. \label {impact} \end {equation} \tag{ 16 }$$

The first order modes that exist in the 1.1 μm < a < 1.6 μm regime are from l = 15 to 23. This leads to values of b_l that range from 1.3124 μm (l = 15) to 1.9898 μm (l = 23). The first minimum in the transverse intensity profile of the Bessel beam occurs at 1.67 μm, which falls approximately in the middle of this range of b_l. This correspondence between the minimum in the intensity profile and the impact parameters suggests a straightforward explanation for the absence or diminishment of the peaks in figure 5(b): the excitation efficiency of the WGMs is greatly reduced as the availability of the number of rays from the incident beam is at a minimum.

The excitation of WGMs was further examined by plotting the internal and near-surface fields of spheres of different radii when illuminated by a plane wave (figures 6(a)–(c)) and a Bessel beam (figures 6(d)–(f)). For the plane wave, field amplitudes characteristic of WGMs are seen for radii that correspond to sharp peaks in figure 5(a). Figures 6(a) and (b) show the calculated fields for a = 1.3333 and 1.7067 μm, respectively. The WGM excitations in figures 6(a) and (b) are the ${\rm TE}_{19}^1$ mode and ${\rm TE}_{25}^1$ mode, respectively. When the radius is not near a sharp peak, no resonance effects are observed (figure 6(c)).

Zoom In Zoom Out Reset image size

Unlike in figure 5(a), a = 1.3333 μm does not correspond to a sharp peak in figure 5(b). The associated field plot in figure 6(d) shows that no WGM is seen despite satisfying the resonant conditions for the ${\rm TE}_{19}^1$ mode. In fact, there is little difference between the fields in figures 6(d) and (f) (a non-resonant case where no excitation is observed). The impact parameter for the ${\rm TE}_{19}^1$ mode is b₁₉ = 1.6511 μm which is very close to the first minimum in the transverse intensity profile of the Bessel beam (1.67 μm). This example illustrates the extent to which the excitation efficiency of a WGM can be diminished in a Bessel beam: it can fall to essentially zero if the impact parameter is near an intensity minimum in the Bessel beam cross-sectional profile.

For a = 1.7067 μm, the impact parameter is not near an intensity minimum of the Bessel beam and excitation of the ${\rm TE}_{25}^1$ mode by the Bessel beam is seen in figure 6(e). This field pattern is very similar to that seen in figure 6(b) when the same mode is excited by the plane wave. Figure 6(f) shows that, similar to the plane wave, no resonance effects are seen in the fields when the radius does not correspond to a sharp peak in figure 5.

In the final size regime (1.6 μm < a < 2.0 μm), all of the sharp peaks in figure 5(a) have corresponding peaks in figure 5(b). The impact parameter is now in the vicinity of the intensity maximum of the first ring of the Bessel beam and, as expected, the excitation efficiency of the WGMs drastically increases.

An interesting effect involving the longitudinal force is observed when a droplet is no longer centered on the propagation axis of the Bessel beam at x = y = 0. In figure 7 the longitudinal forces are shown for droplets located at the intensity maximum of the first ring of the Bessel beam. For the displacement along the x axis, the excitation efficiency of TE modes is much greater than TM modes and, for the displacement along the y axis, the excitation efficiency of TM modes is much greater than TE modes. An example of the difference in mode excitation that can occur is shown in figure 8 for first order l = 11 modes. When the droplet is displaced along the x axis, excitation of the TE mode is seen (figure 8(a)) but not the TM mode (figure 8(b)). Conversely, excitation of the TE mode is not seen when the droplet is displaced along the y axis (figure 8(c)) but the TM mode is excited (figure 8(d)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In contrast to figure 5 where the droplet is centered on the propagation axis, figure 7 does not contain a large region over which the excitation efficiency of all WGMs is poor. The overall lower symmetry of the system drastically changes the excitation efficiency of the WGMs and demonstrates that the coupling between the Bessel beam and the droplet is highly dependent on the position of the droplet within the transverse profile of the beam.

We now include both optical and drag forces to characterise the equilibrium positions of droplets captured by the Bessel beam. Figure 9 shows retention distances for glycerol droplets calculated using equation (15). Results for droplets in a Bessel beam (figures 9(a) and (b)), a plane wave (figure 9(c)), and a Gaussian beam (figure 9(d)) are presented. Sharp peaks are seen throughout all four plots. These are the result of the increase in the longitudinal force that occurs when the radius and refractive index of a droplet are such that the incident beam excites a WGM, displacing the droplet to regions of lower light intensity and giving rise to large excursions in droplet position on resonance.

Zoom In Zoom Out Reset image size

As was discussed in section 4.2, when a droplet is located on the propagation axis of the Bessel beam there will be a range of radii over which the excitation efficiency of the WGMs is poor. In figure 9(a) this occurs from approximately a = 1.1–1.6 μm. Here, peaks in the retention distance are either absent or not as sharp as they are in figures 9(b)–(d). Of course, excitation efficiencies are dependent on the position of the droplet in the Bessel beam and figure 9(b) demonstrates that when a droplet is no longer on the propagation axis of the Bessel beam it is possible to excite the same number of WGMs as a plane wave (figure 9(c)). Finally, figure 9(d) shows the retention distance curve calculated for a droplet in a Gaussian beam. This curve is provided here for reference as previous calculations involving retention distances have all been performed for spheres in Gaussian beams [18, 20, 70]. Given the loose focus of the Gaussian beam it is not surprising that the shapes of figures 9(c) and (d) are the most similar of the four plots.

A single glycerol droplet was trapped using the Bessel beam and gas flow system described in section 2. Figure 10 shows the measured retention distances as the radius of the droplet decreases due to evaporation. When initially trapped, the radius of the droplet was 1.549 ± 0.010 μm. As the radius decreases, the retention distance of the droplet varies. This change was tracked until the radius was equal to 0.823 ± 0.010 μm, at which point the droplet fell out of the trap. Notably, the droplet position varies by ∼0.2 mm over the size range.

Zoom In Zoom Out Reset image size

The sharp peaks that were predicted in figure 9 are observed throughout figure 10. Based on the resonant behaviour in the measured positional data and its correspondence to the calculated curves, there is no doubt that the origin of the observed peaks is WGMs. However, according to calculations performed in section 4.1, the droplet should be located on the propagation axis of the Bessel beam and the excitation efficiency of the WGMs should be greatly diminished across a large portion of the size regime studied. This is not observed. The likely explanation is that because the measured retention distances in figure 10 are all much larger than z_max, the droplet will be outside of the Bessel zone and the beam will be far from ideal. The intensity in between the rings of the Bessel beam will no longer drop to nearly zero and this will allow for an increased excitation efficiency of WGMs whose impact parameters fall within this region. It is also possible that the non-ideality will remove the restriction that the droplet will become trapped at x₀ = y₀ = 0, allowing for equilibrium positions that are located off of the propagation axis. The calculated retention distance curves presented in figure 10 were chosen to be representative of these two cases. For the calculated and observed curves, the correspondence between both the positions and widths of the peaks is quite good as they match extremely well and the WGMs can be easily identified. The two deficiencies of the fit are that the amplitudes in the extent of the excursions of droplet position when resonant are not well represented and the more general location about which the droplet moves is not accurately calculated. For the current experimental setup, proper quantitative modelling of the retention distance appears to be incredibly difficult due to the presence of imperfections in the beam and the inability to determine the position of the droplet in the transverse plane with sufficient precision.

Given the similarities between the retention distance curves for the three different beam types shown in figures 9(b)–(d) it is worth briefly discussing why a Bessel beam was chosen in this work. On a simple first examination of figure 9, it may be thought that simple Gaussian or plane wave illumination may be just as effective at demonstrating aerosol optical chromatography, without the complexities of establishing a Bessel beam trap. However, while the Bessel beam is well-known for the properties of being diffraction free and self-healing, of particular interest here are the strong transverse forces that allow for the confinement of sub-micron droplets within the bright central core. For a droplet in a plane wave, there are no optical forces acting in the transverse direction and such a beam has no practical use when attempting to trap particles. For a droplet in a Gaussian beam, the transverse forces are much weaker than those of the Bessel beam (for the beams and droplet radii considered here). Furthermore, if the transverse forces are enhanced by greater focussing, the degree of motion that can be achieved in the z direction is similarly reduced. Therefore, when working with droplets with radii around one micrometer, the Bessel beam is the most suitable of the three beams, allowing resolution of particles of different size at different longitudinal positions but ensuring tight transverse confinement.

Previous theoretical work on optical chromatography has presented the relationship between retention distance and radius as one where the retention distance increases monotonically with increasing radius (for a dynamic range that spans about three orders of magnitude) [18]. This relationship is not seen in either the observed or calculated retention distance curves presented in figures 9 and 10. The reason for the discrepancy is that, to our knowledge, all previous work in the field of optical chromatography has used ray optics when calculating optical forces [18, 20, 21, 70]. The approximations in such a treatment ignore WGMs, which are the origin of the retention distance peaks in the results presented here.

Further, the influence of WGMs on measured retention distances has not been observed experimentally. It is reasonable that the effect has gone unnoticed as measurements have focussed on coarse separations of systems containing particles with at most four different sizes or refractive indices and a significant variation in each parameter [17, 18, 20, 21, 71, 72]. In contrast, the slow rate of evaporation of a glycerol droplet means that the droplets studied here are a unique probe of retention distance. Changes to the radius of the droplet are recorded on such a fine interval that the measured retention distance contains all of the features of the calculated curves. Even peaks whose widths are less than a few nanometers are observed in figure 10. If such a chromatographic technique is to be used for aerosols, it will be important to ensure that the instrument operates in a regime similar to that identified in figure 9(a), where a more monotonic change in retention distance with size (and refractive index) is observed with variation in particle size. Only in this limit would particles of differing size become segregated and, thus, separated, in different spatial locations.