Omnidirectional gradient force optical trapping in dielectric nanocavities by inverse design
Beñat Martinez de Aguirre Jokisch Benjamin Falkenberg Gøtzsche Philip Trøst Kristensen Martijn Wubs Ole Sigmund Rasmus Ellebæk Christiansen


Omnidirectional gradient force optical trapping in dielectric nanocavities by inverse design]Omnidirectional gradient force optical trapping in dielectric nanocavities by inverse design

[1,2]Beñat Martinez de Aguirre Jokischbmdaj@dtu.dk

2,3]Benjamin Falkenberg Gøtzschebfago@dtu.dk


2,3]Philip Trøst Kristensenptkr@dtu.dk

2,3]Martijn Wubsmwubs@dtu.dk

1,2]Ole Sigmundolsi@dtu.dk

1,2]Rasmus Ellebæk Christiansenraelch@dtu.dk

*[1]Department of Civil and Mechanical Engineering, Technical University of Denmark, Nils Koppels Allé Building 404, Kongens Lyngby, 2800, Denmark

[2]NanoPhoton - Center for Nanophotonics, Technical University of Denmark, Ørsteds Plads 345A, Kongens Lyngby, 2800, Denmark

[3]Department of Electrical and Photonics Engineering, Technical University of Denmark, Ørsteds Plads Building 343, Kongens Lyngby, 2800, Denmark





Optical trapping enables precise control of individual particles of different sizes, such as atoms, molecules, or nanospheres. Optical tweezers provide free-space omnidirectional optical trapping of objects in laboratories around the world. As an alternative to standard macroscopic setups based on lenses, which are inherently bound by the diffraction limit, plasmonic and photonic nanostructures promise trapping by near-field optical effects on the extreme nanoscale. However, the practical design of lossless waveguide-coupled nanostructures capable of trapping deeply sub-wavelength particles in all spatial directions using the gradient force has until now proven insurmountable. In this work, we demonstrate an omnidirectional optical trap realized by inverse-designing fabrication-ready integrated dielectric nanocavities. The sub-wavelength optical trap is designed to rely solely on the gradient force and is thus particle-size agnostic. In particular, we show how a nanometer-sized trapped particle experiences a force strong enough to overcome room-temperature thermal fluctuations. Furthermore, through the robust inverse design framework, we tailor manufacturable devices operating at near-infrared and optical frequencies. Our results open a new regime of levitated optical trapping by achieving a deep trapping potential capable of trapping single sub-wavelength particles in all directions using optical gradient forces. We anticipate potentially groundbreaking applications of the optimized optical trapping system for biomolecular analysis in aqueous environments, levitated cavity-optomechanics, and cold atom physics, constituting an important step towards realizing integrated bio-nanophotonics and mesoscopic quantum mechanical experiments.




 
[
    [
    23 April 2024
=================










Optical tweezers are versatile tools for interdisciplinary research due to their precise control and manipulation of micron-sized objects <cit.>. With applications ranging from studies in microbiology  
<cit.>
to fundamental physics research <cit.>, the optical tweezer is a decorated research tool. 
Undoubtedly, the success of the optical tweezer is owed to its ability to omnidirectionally trap objects by use of a single laser; however, since it is based on free-space optics, it is inherently bound by the diffraction limit. This limit can result in prohibitive power requirements to trap nanometer-scaled particles. As shown in <cit.>, a micrometer-sized polystyrene sphere can be stably trapped with a fraction of a milliwatt, while a 10 nm sphere requires 1.5 W. To overcome this power limitation, nanostructured metallic devices have been employed <cit.>, which can achieve stable trapping of nanometer-scaled particles at much lower powers, utilizing the strongly localized optical forces of plasmonic resonances. Such plasmonic devices have successfully been applied to the optical trapping of single proteins <cit.>, and single particle Raman spectroscopy <cit.>, among others. However, inherent losses in metals lead to heating and potential stability concerns of the electromagnetic resonator or adverse effects on the trapped objects <cit.>. To solve these issues, near-lossless nanostructured dielectric optical traps have been proposed <cit.>. Experimental realizations of dielectric traps include the trapping of single quantum dots <cit.> and the trapping of particles utilizing Fano resonances <cit.>, to name a few. 

Conventionally, due to lower field strength, most dielectric traps show lower trapping stiffness than their plasmonic counterparts. However, recent developments in the design and fabrication of dielectric bowtie-based nanocavities suggest access to hitherto unexplored field strength in dielectrics with deeply sub-wavelength light confinement <cit.>. Bowtie nanocavities have previously been used as efficient optical traps <cit.>, by utilizing the strong field enhancements at the material interfaces <cit.>. However, this field enhancement inevitably makes the particle stick to the resonator walls <cit.>, impeding omnidirectional trapping with gradient forces. For mesoscopic quantum mechanical experiments, any contact with the resonator material is detrimental to the coherence of the prepared state <cit.>, and in biomolecular analysis, it may lead to undesired charge reconfiguration <cit.>.

In this work, we directly address the issue of omnidirectional trapping in nanostructures by inverse-designing dielectric nanocavities that trap particles with sizes significantly below the diffraction limit. A recent study <cit.> addresses the issue of omnidirectionality for a limited range of particle sizes, by utilizing gradient forces and self-induced back-action (SIBA) effects <cit.>. Here, by only relying on gradient forces, we deterministically tailor the device geometry of a nanostructure to feature a particle-size independent omnidirectional trapping potential. 
The inverse design process relies on topology optimization (TO), a design optimization method widely used in the design of optical applications like waveguides <cit.>, cavities <cit.>, demultiplexers <cit.>, microresonators <cit.>, and more. Experimental evidence <cit.> shows that topology-optimized structures can directly incorporate manufacturing constraints <cit.>, ensuring precise fabrication of optimized device blueprints. 
In a recent work, a plasmonic optical trap was designed using TO by maximizing the electric-field strength at the center of a cavity <cit.>, yielding the well-known bowtie-like structure, that can trap particles in two dimensions within the nanoscale gap. In this work, we employ a novel TO scheme based on fitting the electromagnetic field profile to a desired Gaussian shape, enabling omnidirectional trapping within a specified volume. 
Applying the TO procedure at different wavelengths and background materials, we design waveguide-coupled and fabrication-ready devices, that may pave the way for future on-chip levitated optomechanics, on-chip biomolecular analysis, and cold-atom-based quantum many-body systems <cit.>. 



§ INVERSE-DESIGNING OMNIDIRECTIONAL TRAPPING



To trap particles with a radius R ≪λ, where λ is the wavelength of light, we may apply the dipole approximation <cit.>. To show that our results apply to a wide range of lossless particles, we consider a non-resonant spherical particle with refractive index n=2 and radius R=15 nm, representative for both proteins (n ≃ 1.6, R ∈ [1 nm- 100 nm]) <cit.> and semiconductor quantum dots (n ≃ 2.4,  R ∈ [2 nm- 50 nm]) <cit.>. 
For our choice of lossless and non-resonant particles, scattering forces are negligible and the force may thus be completely described by the gradient force. This conservative force can be described as the gradient of the trapping potential <cit.>, 

    U(𝐫)=-α_R/4[𝐄^*(𝐫)·𝐄(𝐫)] ,

where 𝐄 is the electric field and α_R is the real part of the spherical particle's polarizability, as given by the Clausius-Mossotti relation <cit.>. Therefore, for a given particle in the nanostructure, the trapping potential is directly described by calculating the electric field distribution of the empty cavity. The dipole approximation may, however, break down for particles that are larger, exhibit loss, are resonant, or possess a higher refractive index, where additional effects like SIBA or scattering forces (e.g., radiation pressure or spin-curl forces) <cit.> may need to be accounted for. In post-analysis, we verify that our devices are unaffected by these effects.


To model optical trapping in the dipole approximation, we calculate the electric field of the empty cavity, by solving Maxwell's equations. The trapping device consists of two silicon waveguides connected to a central design domain, in which TO is applied, see [fig:1]Fig. 1a. The central region contains a cylindrical air exclusion region with a radius of  R_exc=300 nm and thickness of 800 nm, in which the stable trap is located. One of the waveguides is excited with the fundamental mode at an input power of P_in=15 mW and at a free-space wavelength of λ. The wavelength is a freely selectable design parameter, determined by the desired application of the trap. Here, as a demonstration, we inverse design a device in the near-infrared regime (λ = 1.55µm), and in one of the following sections, we design an optical trap in the optical regime (λ=775 nm). The material distribution in the design domain is optimized to obtain an electric field distribution that results in an omnidirectional trapping potential. Accordingly, we formulate our design problem as a continuous optimization problem where the material distribution is controlled by our design parameters and where we seek to minimize a Figure of Merit (FOM) that defines the difference of the electric-field norm with respect to a reference field 𝐄_ref. 
This expression may be written as:

    FOM≡Φ= √(∫_Ω[ Θ(𝐄(𝐫)/𝐄(𝐫_0)-𝐄_ref(𝐫)/𝐄_ref(𝐫_0))]^2  dΩ) ,

where 𝐫_0 is the center point in the design domain, Θ(x) is a smoothed Heaviside threshold function <cit.> and Ω is the optimization domain defined by the exclusion region, see [fig:1]Fig. 1a. To ensure sub-wavelength omnidirectional trapping, we select the reference field to be a three-dimensional Gaussian potential with standard deviations σ_x=σ_y=300 nm and σ_z=400 nm, to ensure that it features a stable trapping minimum in all spatial directions. In this expression, the Heaviside projection Θ(x) ensures that the FOM promotes only electric field distributions as steep as, or steeper than the target field. The standard deviations are chosen to match the dimensions of the cylindrical exclusion nanocavity with a volume V_cav below the diffraction limit: V_cav=0.22 µm^3< (λ/2)^3=0.465 µm^3. For more information on the inverse design framework, see Methods.



 §.§ Inverse-designed nanocavity in the near-infrared


The inverse-designed omnidirectional trap and its key characteristics are presented in <ref>. In [fig:1]Fig. 1.a we show the electric-field intensity for the optimized structure in the near-infrared, which has a well-defined bell-shaped curve in all spatial directions. This is confirmed by observing the different plane- and line-cuts for the trapping potential in [fig:1]Fig. 1.b. One can identify a minimum of the potential at the optimization region's center, resulting in an omnidirectional trapping potential. Interestingly, the potential has a similar shape in all directions, meaning that the Heaviside projection Θ(x) in the FOM allowed the potential to become steeper than the reference Gaussian along the z-axis. For the 15 mW of optical input power and the reference particle with R=15 nm and n=2, the trapping minimum is stable against thermal fluctuations at room temperature (T=300 K), which conventionally requires a characteristic trapping depth of U ≃ 10  k_B T <cit.>, where k_B is the Boltzmann constant. Remarkably, to achieve this trapping depth, the optimized design efficiently couples and confines light into the cavity while using an input power 2 orders of magnitude less than conventional tweezers (∼1 W <cit.>). Notably, we find that devices optimized for smaller exclusion radii R_exc and associated standard deviations of the reference field, yield deeper trapping potentials. However, this gain in trapping depth crucially comes at the expense of achieving omnidirectional trapping. This sets an upper bound on the minimum length scale of R_exc before the omnidirectionality of the trap is lost, in our case at around R_exc≃ 250 nm. For more details on this study refer to section S2.3 in the Supplementary Information (SI).



From the potential in the dipole approximation, we calculate the axial components of the force experienced by the particle in [fig:1_2]Fig. 2.a. Specifically, we show that close to the origin, the force becomes linear to a good approximation, and thus, the particle-trap system behaves as a linear spring-mass system. The force in the linear region can then be written as, F_i(r_i)=κ_i· r_i, where i ∈{ x,y,z } is an axis-index, κ is a vector with the axial components of the trapping stiffness and r_i is the position on the axis. Fitting the curves in the central region with a linear function yields trapping stiffnesses of κ_x=0.53 fN/nm, κ_y=0.46 fN/nm, and κ_z=0.51 fN/nm, which are comparable to state-of-the-art (SOTA) dielectric devices <cit.>. To validate the dipole approximation, we solve the full Maxwell's equations with the sphere (R=15 nm and n=2) present in the optimized geometry, and calculate the force using the Maxwell Stress Tensor (MST) formalism <cit.>. In [fig:1_2]Fig. 2.b we show the axial components of the force experienced by the particle when displaced in increments of 35 nm in the three directions. By comparing the force derived from the MST calculation and the dipole approximation, we confirm that there is excellent agreement between the results and that there are no other dominant effects, such as SIBA <cit.>, or scattering forces <cit.>.



We have shown that once the particle is at the center of the trap, it will be omnidirectionally trapped. However, loading the particle into the trap may potentially be jeopardized by other physical effects, such as the field enhancement observable at the bottom interface of the structure in [fig:1]Fig. 1.a or the Casimir-Polder (CP) forces near the material interfaces.  Analyzing the trapping potential landscape near material interfaces reveals that CP forces dominate over optical forces. These forces, however, do not compromise the loading of the trap, as long as the particle remains further away than 20 nm from the device walls. To quantify this, we calculate that the worst-case probability of a particle getting stuck at the device walls during the loading of the trap is less than 12%. Moreover, we find that for traps with deeply subwavelength exclusion region radii R_exc the CP-force dominates optical forces, making omnidirectional trapping of levitated particles unattainable. See section S3 in the SI for details. 



 §.§ Benchmarking trapping performance




In the optical trapping community, different metrics are used to benchmark optical trapping platforms of lossless particles; among others, power normalized trapping depth  <cit.>, trapping force <cit.>, and trapping stiffness <cit.> are commonly used. These measures may depend on the physical properties of the trapped particle and its environment, the device footprint, or the input power, making comparisons difficult. Therefore, we propose a normalized trapping stiffness metric (η_i) for optical trapping that normalizes trapping stiffness to particle volume, to the contrast between particle and background dielectric permittivity, and to input power,

    η_i = κ_i  ε_0/α_R  P_in ,

where ε_0 is the vacuum permittivity and i ∈{ x,y,z } is an axis-index. This metric, which characterizes the power efficiency of optical forces per distance, enables one-to-one comparisons of optical trapping among different platforms. For our trap it yields efficiencies of η_x = 1.67 pN/(µm^4·µW), η_y = 1.44 pN/(µm^4·µW) and η_z = 1.6 pN/(µm^4· µW) along the three axial directions. We note that SOTA plasmonic designs, in general, show larger trapping stiffnesses but do not provide omnidirectional trapping. In addition, plasmonic devices come at the expense of optical losses to heating, which can compromise trapping characteristics <cit.>. Compared to other dielectric platforms  <cit.>, the present design shows similar normalized trapping stiffnesses while simultaneously providing omnidirectional trapping through gradient forces. Compared to optical tweezers, which provide omnidirectional trapping, we achieve similar trapping stiffnesses at significantly lower input power, thus yielding higher normalized trapping stiffnesses. This demonstrates the usefulness of miniaturizing the free-space optics by means of a waveguide-coupled nanostructured device. 
For more details on the comparison with SOTA devices, refer to section S4 in the SI.




 §.§ Inverse-designed nanocavities for omnidirectional trapping applications




The miniaturized integrated circuit shown in [fig:1]Fig. 1.a bridges the omnidirectional trapping of optical tweezers with the near-field optics of ultra-compact nanostructures, thereby enabling a range of chip-scale applications. As an example, we demonstrate its potential in levitated cavity optomechanics. Owing to the harmonic potential close to the cavity center, we calculate the natural frequencies of the harmonic oscillator for the different axes as Ω_0,i = √(κ_i/m), where m is the mass of the particle and i ∈{ x,y,z } is an axis-index. Using a density of ρ=2 g/cm^3, which represents particles like proteins (ρ∈ [1.4 g/cm^3 - 1.5  g/cm^3]) <cit.> and quantum dots (e.g. ρ_Si=2.33 g/cm^3, ρ_GaAs=5.32 g/cm^3), we find Ω_0,x=4.33 rad·MHz, Ω_0,y=4.03 rad·MHz and Ω_0,z=4.25 rad·MHz. In the quantum-mechanical limit, the mean thermal occupancy of the mechanical energy states is given by ⟨ n⟩=k_B T / ħΩ_0 <cit.>, where ħ is the reduced Planck constant. Resolving the quantum ground state requires ⟨ n ⟩ < 1, which with the natural frequencies of this system, requires center-of-mass equilibrium temperatures of T_x = 0.21 mK, T_y = 0.19 mK and T_z= 0.20 mK. These temperatures are two orders of magnitude higher than for conventional optical tweezers <cit.>, which combined with the particle-size agnostic omnidirectional trapping, facilitates mesoscopic quantum optomechanical experiments with optically trapped and cooled nanoparticles. One could reach these temperatures by lowering the ambient temperature or employing methods like parametric feedback cooling <cit.>. A representation of the optomechanical system is shown in [fig:2]Fig. 2.a.





To demonstrate the versatility of the design framework, we move towards visible optical frequencies and design a device suited for integrated biological sensing. For biological characterization in aqueous environments, it is crucial to operate within the biological window <cit.>, a spectral range minimizing the optical excitation of vibrational modes in water molecules. This minimizes heating the environment, which could otherwise compromise the trapping setup  <cit.>. To this end, we scale all geometric design parameters and the wavelength by half, shifting towards the visible optical frequencies at λ=775 nm. Accordingly, we set the background refractive index to water and adjust the refractive index of silicon. 
The optimized trap in [fig:2]Fig. 2.b yields an omnidirectional stable trap for a single particle at the center of the design, with trapping stiffnesses of κ_x=7.49 fN/nm, κ_y=4.49 fN/nm, κ_z=4.77 fN/nm, for a particle with radius R=15 nm and for an input power of P_in=15 mW.  The normalized trapping stiffnesses for this device are η_x = 39.78 pN/(µm^4·µW), η_y = 23.84 pN/(µm^4·µW) and η_z = 25.33 pN/(µm^4·µW). These values are higher than those of the device operating in the near-infrared due to the compactness of the device, given the same input power and particle size. The inverse-designed optical trapping device is thus an integrated optical component capable of stably trapping biomolecules in all spatial dimensions in an aqueous solution. One could envision using the platform, represented in [fig:2]Fig. 2.b, integrated with an optical detection scheme to trap and detect particles in situ in biological environments, opening a path for the development of new experiments and technologies in integrated biophotonics systems.




§ DISCUSSION

The omnidirectional trapping of sub-wavelength nanoparticles in integrated optical devices is central to many applications in the field of optical trapping, such as microbiology <cit.>, biophysics <cit.>, or fundamental physics <cit.>. For deeply sub-wavelength particles, where the dipole approximation is valid, we demonstrate that it is possible to tailor the electric field distribution and the trapping potential, by careful nanostructuring of the dielectric environment. The target distribution of the electric fields results in a particle size agnostic omnidirectional optical trap based on gradient forces. To test the geometric limits of omnidirectional trapping, in the SI, we demonstrate how CP forces set a lower size limit for the exclusion radius (R_exc), while a parametric study indicates omnidirectional trapping only for exclusion radii above R_exc≥ 250 nm. To compare the optimized device to other nanophotonic platforms we propose the metric of normalized trapping stiffness (η), which shows unprecedented values compared to other SOTA omnidirectional traps <cit.>, highlighting the deeply sub-wavelength nature of the trapped particles and the low input power required for trapping. Moreover, our framework enables the design of manufacturable <cit.> optical traps for a given application by rescaling the wavelength and the spatial dimensions of the device and adjusting the material parameters accordingly, before applying the inverse design process. We have demonstrated this in the near-infrared (λ=1.55 nm) and optical frequencies (λ=775 nm). Since the designs can be fabricated by standard electron-beam lithography, we anticipate experimental realizations of the optical traps, with applications in levitated cavity optomechanics and integrated biophotonics technologies.

We also foresee developments in our inverse design framework. For instance, modifications to the FOM and the optimization framework can enable the design of optical traps based on different materials, traps for lossy or resonant particles, traps with multiple trapping spots, or novel optical traps for multiple quantum emitters, such as quantum dots or cold atoms. These networks of trapped emitters may, in turn, be used to study and develop quantum many-body systems <cit.>.



§ METHODS




 §.§ Inverse design framework


We simulate the electromagnetic fields using Maxwell’s equations in the frequency domain for the inverse design process, assuming time-harmonic behavior <cit.>. The model is discretized and solved using the finite-element method with first-order Nedelec elements <cit.>. From the electric field, we calculate the FOM in <ref>, which is sought to be minimized through TO. Additionally, we control the minimum depth of the trapping potential by adding a constraint to the optimization problem, which prescribes a minimum electric-field norm value at 𝐫_0. By enforcing the constraint, the optimizer avoids local minima where the trapping potential has the correct shape but is not deep enough for stable trapping. To ensure the manufacturability of the device, we add a constraint to connect the design to the two waveguide ends <cit.> as well as two minimum length scale constraints that act on the solid and void regions respectively <cit.>. To minimize the constrained FOM, we apply TO on the design domain, introducing one design variable per finite element in our discretized design domain, which is used to interpolate between air (0) and silicon (1). The design variables are fixed to only vary in the (x,y) plane and are linked in the z direction <cit.>. We apply a filtering and thresholding procedure to regularize the design <cit.>. The filtered and thresholded design variables are related to the material's refractive index through a material interpolation scheme <cit.>. The optimized design is obtained by solving the topology optimization problem from a single uniform initial guess, using the globally convergent method of moving asymptotes as the optimizer <cit.>. For more information on the forward problem, the inverse design framework, and the full optimization problem see the SI.





Acknowledgements

We would like to thank Jesper Mørk for the useful discussions and support. BFG, PTK, and MW would also like to thank Kishan Dholakia for the useful discussions. We gratefully acknowledge financial support from the Danish National Research Foundation through NanoPhoton - Center for Nanophotonics, grant number DNRF147. 










 
 
 
 
 
 

 
Fig.





    Supplementary information for
    

    Omnidirectional gradient force optical trapping in dielectric nanocavities by inverse design
    






§ FORCES IN OPTICAL TRAPPING


To design a nanostructure that can omnidirectionally trap single particles we need to perform optical force calculations. For particles with a radius R ≪λ, where λ is the wavelength of light, one may apply the dipole approximation. Calculating the cycle-averaged force acting on a point-dipole for monochromatic electromagnetic fields yields <cit.>:

    ⟨𝐅⟩ = 𝐅_grad + 𝐅_rad + 𝐅_SC ,

where 𝐅_grad is the gradient force, 𝐅_rad is the radiation pressure force and 𝐅_SC is the spin-curl force <cit.>. The gradient force is given by:

    𝐅_grad = α_R/4∇ [𝐄^*·𝐄] ,

where 𝐄 is the electric field, α_R and is the real part of the polarizability. The radiation pressure force is given by:

    𝐅_rad = σ_p/c⟨𝐒⟩ ,

where  σ_p = α_Ik/ε_0 is the particle's total cross-section, α_I is the imaginary part of the polarizability, k is the wave-number, ε_0 is the vacuum permittivity and c is the speed of light. The spin-curl force is given by:

    𝐅_SC = σ_p c [ ∇×⟨𝐋⟩] ,

where ⟨𝐒⟩ = 1/2Re{𝐄×𝐇^*} is the cycle-averaged Poynting vector, 𝐇 is the electric field, ⟨𝐋⟩ = ε_0/4iω 𝐄×𝐄^* is the cycle-averaged spin density of the electromagnetic field and ω is the angular frequency.

For spheres in a background medium with permittivity ε_back, the total polarizability is given by


    α = α_0/1-ik^3/6πε_0α_0 = α_R + iα_I ,   where  α_0 = 4 πε_0 R^3 ε-ε_back/ε+2ε_back ,

where α_0 is the Clausius-Mossotti polarizability <cit.> and  ε(ω) is a frequency-dependent dielectric permittivity. For lossless and non-resonant materials α_R≫α_I and the force can be completely described by the conservative gradient force, which is generated by a potential:

    U(𝐫)=-α_R/4[𝐄^*(𝐫)·𝐄(𝐫)] ,

and the real part of the polarizability is described entirely by the Clausius-Mossotti equation α_R≃α_0. Note that to obtain an omnidirectional trapping potential, given that the force is 𝐅(𝐫)=-∇ U(𝐫), it requires that the expression in <ref> has a single minimum in all spatial directions. The particle will experience a net zero force when it is at the minimum of the potential landscape, while if the particle is displaced from there in any direction, it will feel an attractive force pushing it back to the minimum, yielding omnidirectional trapping.

Calculating the cycle-averaged force acting on a particle directly from Maxwell's equations is also possible, without employing the dipole approximation. For a particle enclosed by the surface ∂ V  the force is given by <cit.>:

    ⟨𝐅⟩=∫_∂ V⟨↔𝐓(𝐫, t)⟩·𝐧(𝐫) d a

where 𝐧 defines the vector normal to the particle surface and ↔𝐓 is known as the Maxwell stress tensor and can be written as

    ↔𝐓=[ε_0 ε𝐄 𝐄+μ_0 μ𝐇𝐇-1/2(ε_0 ε E^2+μ_0 μ H^2) ↔𝐈] ,

where 𝐇 is the magnetic field, ε and μ are the relative dielectric permittivity and permeability of the medium surrounding the particle, and ε_0 and μ_0 are the free-space permittivity and permeability. The expression in <ref> has been used to calculate the force acting on a particle introduced into the optimized optical trap to validate the results in the dipole approximation.




§ INVERSE DESIGN FRAMEWORK




 §.§ The forward problem


To inverse design the optical trap, we first need to define and solve an appropriate forward problem in order to model the nanophotonic system. This means solving Maxwell’s equations in the frequency domain, assuming time-harmonic field behavior <cit.>. The model in [fig:S1]Fig. S1.a is discretized and solved using the finite element method with first-order Nedelec elements <cit.>. The model is based on a simulation domain consisting of two silicon optical waveguides (in dark gray) connected to a design region Ω_D (in blue) and an air cladding (in light gray). For the design in the near-infrared (λ=1.55µm), the dimensions (x, y, z) for the simulation domain are (L_sim, w_sim, h_sim)=(9.1µm, 4µm, 1.5µm), which are large enough to allow the fields to decay away from the cavity, avoiding artificial numerical boundary reflections[This is also the case at the waveguide ends. The waveguide is excited with the fundamental mode at a port located a wavelength (λ) away from the edge of the simulation domain. This allows the fields reflected in the waveguide to decay away from the input port.]. The waveguide has dimensions of (L_wg, w_wg, h_wg)=(2.8µm, 275  nm, 400 nm). The height of the waveguide was chosen to be at least half of the wavelength since this was the smallest value yielding omnidirectional trapping. The dimensions for the design domain are (2L_Ω_D,L_Ω_D, h_wg)=(2.5µm, 1.25µm, 400 nm) and are chosen to have a footprint of around 5(λ/n)^2, where n is the refractive index of silicon. This device size is large enough for the topology optimization to design couplers and mirrors. Lastly, in the center of the design domain, there is a cylindrical exclusion region to trap the particle, with radius R_exc=300 nm. This radius has been carefully chosen to obtain omnidirectional optical traps since there is a trade-off between the compactness of the exclusion region radius and the possibility of obtaining omnidirectional trapping, as further discussed in <ref>. Note that for the integrated biophotonic design the wavelength and all dimensions are halved.


To solve the problem in a computationally efficient way, we assume that the device is symmetric around the (x,y) and (x,z) planes, leaving only a quarter of the total simulation domain to be solved. To impose the symmetry we apply perfect electric conductor (PEC) boundary conditions on the (x,z) plane and perfect magnetic conductor boundary (PMC) conditions on the (x,y) plane. On the rest of the boundaries, we apply first-order absorbing boundary conditions. Lastly, PML regions with a length of λ are defined at the ends of both waveguides to avoid reflections at the ends of the simulation domain. 





 §.§ The optimization problem


Having computed the solution to the forward problem for a given material distribution in the modeling domain, we compute the Figure of Merit (FOM) Φ, which is to be optimized. As outlined in the main document, the FOM defines the difference of the electric-field norm with respect to a reference field, which is directly related to the shape of the potential and is given by:

    FOM≡Φ= √(∫_Ω[ Θ(𝐄(𝐫)/𝐄(𝐫_0)-𝐄_ref(𝐫)/𝐄_ref(𝐫_0))]^2  dΩ) ,

where 𝐫_0=(x_0, y_0, z_0) = (0,0,0) is the center of the design domain, 𝐄_ref is a reference electric field, Ω is the optimization domain defined by the exclusion region in [fig:S1]Fig. S1.a, and Θ(x) is a smoothed Heaviside threshold function <cit.>. The target reference field is chosen to be a three-dimensional Gaussian of the form:

    𝐄_ref(𝐫)/𝐄_ref(𝐫_0) = exp(-((x-x_0)^2/2 σ_x^2+(y-y_0)^2/2 σ_y^2 + (z-z_0)^2/2 σ_z^2)) .

where σ_x = σ_y = 300 nm and σ_z= 400 nm are the standard deviations in all directions, chosen to match the exclusion size. In [fig:S1]Fig. S1.b we show a projection of the field described by <ref> for the (x,y) plane, which showcases that the Gaussian has a single maximum at 𝐫=𝐫_0. As pointed out in <ref>, if the electric-field norm has a single maximum, this will result in a global minimum of the trapping potential, which is necessary to trap single particles omnidirectionally. For a discussion on the choice of the standard deviations σ_i, where i ∈{ x,y,z }, refer to <ref>.

The FOM is minimized by optimizing the material distribution of silicon (Si) and air in the design region (Ω_D). We formulate the design problem as a continuous optimization problem, where the material distribution is controlled by a design field ξ, which is discretized into a piecewise constant field coinciding with the finite elements used to discretize the physics model, with one design variable controlling the value of the design field in each element. To regularize the design and enable length-scale control, we adopt a filtering and
thresholding scheme. For the filter, we use a Helmholtz-based filter <cit.>:

    -(r_f/2√(3))^2∇ξ̃+ξ̃=ξ ,

where ξ̃ is the filtered design field, and r_f is the filter radius. The filter operation is followed by a smoothed Heaviside threshold (Θ) <cit.>:

    ξ̅̃̅ = Θ(ξ̃
        ) = tanh (β·η)+tanh (β·(ξ̃-η))/tanh (β·η)+tanh (β·(1-η)),   β∈[1, ∞), η∈[0,1] ,

where ξ̅̃̅ is the filtered and thresholded design field, and β and η control the threshold sharpness and value respectively. To translate the design field into the material distribution in the physics model we employ a non-linear material
interpolation <cit.>:

    ε_r(ξ̃)=(n(ξ̃)^2-k(ξ̃)^2)-i(2 n(ξ̃) k(ξ̃))-iαξ̃(1-ξ̃), 
        n(ξ̃)=n_air +ξ̃(n_Si -n_air ), 
        k(ξ̃)=k_air +ξ̃(k_Si -k_air ),

where ε_r is the relative dielectric permittivity, n is the refractive index, k is the extinction coefficient, α is a problem-dependent parameter known as the artificial attenuation, and  “Si" stands for silicon. With this setup we optimize the design using the method of moving asymptotes <cit.>  with a single homogeneous design guess (ξ = ξ_0), leaving the optimizer free to tailor the device geometry. The optimizer solves the following optimization problem:


    min _ξ:      FOM≡Φ= √(∫_Ω[ Θ(𝐄(𝐫)/𝐄(𝐫_0)-𝐄_ref(𝐫)/𝐄_ref(𝐫_0))]^2  dΩ) , 
     s.t. :    𝐒(ε_r(ξ̅̃̅, 𝐫)) 𝐄(𝐫) = 𝐅(𝐫)  , 
    
      :     log_10(||𝐄(𝐫_0)||)≥γ ,
    
      :     ∇· (-c   Θ(ξ̅̃̅(𝐫)) ∇ C(𝐫) ) =f Θ(ξ̅̃̅(𝐫)),    C=0    ∀𝐫∈Γ_i,    i∈(1,2) , 
    
      :       g^s = 1/n∑_i ∈ℕ I^s_i [min{(ξ̃_i - η_e), 0}]^2 ≤ϵ , 
    
      :       g^v = 1/n∑_i ∈ℕ I^v_i [min{(η_d - ξ̃_i ), 0}]^2 ≤ϵ , 
    
    
    
    
     :      0<ξ(𝐫)<1  ,
    
     :     ξ=0   ∀𝐫∈Ω_D   ,


where the FOM in <ref> is subject to the constraints given by the individual subequations:

    
  * <ref> is the discretized form of Maxwell's equations, which is equivalent to solving a linear algebraic system that yields the electric field for the system matrix 𝐒 and the excitation term 𝐅.
    
  * <ref> is a constraint for the electric-field norm in the center of the domain, where γ is a problem-dependent parameter. By selecting γ=1.15 log_10(||𝐄(𝐫_0)||_i=0), where i=0 refers to the initial design, we ensure that in the initial steps of the optimization, the optimizer avoids local minima where the trapping potential has the correct shape but is not deep enough for stable trapping. Once the constraint in <ref> is fulfilled, it ensures a stable trapping potential that overcomes thermal diffusion, and the optimizer directly targets the FOM in <ref> to achieve the correct shape of the trapping potential.
    
  * <ref> refers to two connectivity constraints formulated using an artificial heat-transfer problem, that ensures a design connected to the two waveguides. In this expression, c is a material interpolation for the artificial conductivity, C denotes the artificial temperature field, f denotes the artificial heat generated by materials and  Γ_i describes the heat sink boundaries given by our two (i=1,2) waveguide ends <cit.>. The rest of the parameters are selected as detailed in <cit.>.
    
  * The constraints g_s and g_v, in <ref> and <ref>, are the solid and void connectivity constraints, where η_e and η_d are the eroded and dilated thresholds, n denotes the number of elements, I^s_i=ξ̅̃̅_i ·e^-c·|∇ξ̃_i|^2 and I^v_i=(1-ξ̅̃̅_i) ·e^-c·|∇ξ̃_i|^2 are the solid and void indicator functions, c is a problem-dependent parameter and ϵ is the length scale error <cit.>. To ensure the minimum length scale of the design features we select c = r_f^4 <cit.>, where r_f is the filter radius, and η_e=0.75, η_d=0.25 which in our optimization problem ensure that no features in the design have a radius of curvature smaller than 60 nm <cit.>, making the design variable using e.g. standard electron beam lithography.


We run the optimization problem in [eq:optimization:FOM]Equations S11a-S11h with a single homogenous initial guess of the design field ξ=0.6 for all the design variables. Additionally, we use a continuation scheme to exploit the continuous nature of the design field in solving the optimization problem while achieving a final, physically realizable binarized design. This consists of increasing the parameter that controls the threshold sharpness (β) and the artificial attenuation (α) every 50 iterations of the optimization, pushing the design field toward binary values and thus physical realizability. We also gradually reduce the length scale error (ϵ) introduced in [eq:optimization:FOM]Equations S11e-S11f to ensure that the final design fulfills the minimum length scale requirement. The continuation scheme parameters are summarized in <ref>. 





 §.§ Tuning the exclusion radius: from omnidirectional trapping to strong light confinement


One of the parameters that requires the most careful tuning in the inverse design procedure is the radius of the exclusion region (R_exc) defined in <ref>. In principle, to achieve the highest possible trapping stiffness, which is defined as the force exerted over a distance, one would try to make the exclusion region as small as possible. This would allow for stronger light confinement on a sub-wavelength scale, enhancing the trapping stiffness. However, through systematic optimizations[For each optimization the intensity constraint in  <ref> was normalized to the first iteration by choosing γ=1.15 log_10(||𝐄(𝐫_0)||_i=0), while the rest of the parameters remain the same as in the original optimization problem presented in <ref>.] selecting different exclusion radii, we find a trade-off between the increased trapping stiffness and the possibility of achieving omnidirectional trapping with gradient forces. This is summarized in <ref>, where for decreasing exclusion region radii, we optimize six nanocavity designs and calculate their trapping potential. The cavities optimized for smaller exclusion region radii can achieve strong field enhancements below the diffraction limit, providing deep trapping potentials. In fact, and based on the bowtie-like cavity design for the smallest exclusion region (R_exc=50 nm), we see that if dielectric material was allowed in the exclusion region, the optimizer would create extreme dielectric confinement bowtie-like cavities, similar to ones in <cit.>. However, similar to bowtie-like cavity designs <cit.>, when checking the trapping potential distribution for exclusion radii below R_exc=200 nm in <ref>, it is clear that the strong field confinement located at the material interfaces does not allow for an omnidirectional trap in the center of the cavity, due to the spatial overlap of the field from the tip-like structures close to the center. This is similar to the double nanohole structure with a deeply sub-wavelength exclusion region found by another inverse design work <cit.> where the electric field was maximized at the center of the cavity, yielding a bowtie-like non-omnidirectional trapping potential. Interestingly, even for an exclusion radius of R_exc=250 nm the trapping potential in the y direction starts to lose its concavity, which is further accentuated for smaller radii. In other words, there is a minimum length scale of the exclusion radius close to R_exc≃ 250 nm, which sets a geometrical limit on omnidirectional trapping. By observing the projection of the trapping potential onto the cartesian axes it seems however, that one recovers stable trapping when going down to the smallest exclusion radius (R_exc=50 nm), but by looking at the two-dimensional profile of the trapping potential it is evident that a particle would be trapped at the nanocavity tips, where the field enhancement is strongest. In our work, we choose an exclusion radius R_exc=300 nm, since it simultaneously allows for an omnidirectional trapping without compromising the depth of the trapping potential or the trapping stiffness. Moreover, and as shown in the main document, for an exclusion region radius  R_exc=300 nm, we are still able to trap deeply sub-wavelength particles in a sub-wavelength exclusion region while overcoming thermal fluctuations.





§ SURFACE FORCES

A particle located at the center of our optimized device will be omnidirectionally trapped. There are, however, two physical effects that could potentially compromise loading the particle into the optical trap: optical lightning-rod effects <cit.> and vacuum fluctuation-induced Casimir-Polder (CP) forces <cit.>. We compare these two surface forces and show how the CP forces have a stronger probability of resulting in particles sticking to the interface than any lightning-rod effects present in the design. Finally, we show that, in the worst-case scenario, there is only a 12 % percent chance of particles sticking to the walls due to the CP force.



 §.§ Lightning-rod effects

By taking a close look at the field distribution at the center of the nanocavity in Fig. 1.a in the main document, one can see a significant local field enhancement at the material interface for z = -400 nm. This is an inherent consequence of the boundary conditions of Maxwell's equations <cit.>, where the field is locally enhanced
at kinks and corners <cit.>. This so-called lightning-rod effect is strongly confined to the interface but still creates an attractive force on the particle if the particle is sufficiently far from the center, and sufficiently close to the interface. 

To quantify this effect, we investigate the trapping potentials for the y line at different out-of-plane (z axis) heights in <ref>[These results have been validated for a 15 nm radius particle using the MST formalism.]. At the bottom of the device (z=-400 nm), we see that there is a region y ∈[-220 nm, 220 nm] where the particle will be stably trapped, feeling a restorative force towards the center. In the other region (in gray) the particle will feel an attractive force towards the interface, which will lead to the particle sticking to the interface. If we go up inside the structure (z=-330 nm), or down out of the structure (z=-450 nm) we see the sticking effect disappearing. This means that if we assume an initial random position of the particle above the hole in the simulation domain in <ref>, we can compute the probability that it will get stuck at the interface as it travels into the omnidirectional trap. We do this by assuming a simple positional averaging for the particle as it travels into the trap and then calculating the volume ratio between the volume where the particle feels an attractive force to the interface and the total volume. As we go up or down the structure, the region where the particle will get stuck to the interface (in gray in <ref>), becomes smaller until it disappears entirely as we reach the limits z=-450 nm and z=-330 nm. To simplify the calculations and provide a worst-case estimate, we assume that the width of this region in the y axis remains constant as we move up and down in the structure. With this assumption the probability of sticking is

    P_stick^LR  [%]= V_stick/V_total =  h_stick(R^2_exc - R^2_stick )/ h_sim R^2_exc≃ 4.5 % ,

where “LR" denotes lightning-rod, h_stick=120 nm is the height of the region where the particle can stick to the surface, in the interval  z∈[-450 nm, -330 nm]. Therefore, there is a small probability that the particle might stick to the interface due to the lightning-rod effect at the top and bottom interfaces of the structure, as it is being loaded into the trap. It is noted that this is a conservative approximation since it assumes that the unstable effect does not fade out as we go into or out of the structure and does not account for the inherent velocity of the particle, or other forces that could pull the particle out of the unstable region. Thus, our estimate constitutes a worst-case sticking probability P_stick^LR, and a significantly smaller value is expected in reality.  




 §.§ Casimir-Polder forces


Aside from the optical forces, the other key surface forces to consider are the CP forces <cit.>, which are vacuum forces that also exist in the absence of a light source. In the same way that the lightning-rod effects might be problematic and decrease the probability of correctly loading the particle into the omnidirectional trap, the CP forces result in a net attractive force to the surface that scales non-linearly with the distance between the particle and the surface. To calculate the probability of the particle getting stuck at the interface due to CP forces, we approximate the nanosphere as a polarizable point-dipole, close to an infinite dielectric plane, which is a good approximation since the CP forces are known to dominate only at short length scales. Given the high symmetry of this system, the CP energy landscape of this system reduces to a function of the distance d between the
center of the particle and the plane. Following <cit.>, we calculate this function by integration of the Green tensor as formulated in <cit.>. We show the CP energy landscape as a normalized trapping potential in <ref> at room temperature (where k_B T = 0.025 eV), for a particle with a refractive index n=2 and a radius of R=15 nm, close to an infinite silicon interface with a refractive index n=3.48.



In the simplified model, the trapping potential diverges toward infinity as we approach the interface. To calculate the probability of the particle sticking to the interface due to the CP forces, we again consider the stable and unstable regimes in a cylinder extruded from the exclusion hole. In contrast to the lightning-rod effect, the CP forces will affect all interfaces. Given that the potential at the edges of our exclusion domain is around  U≃ -7.5  k_BT, the CP potential will start to dominate the gradient force potential around d≃ 20 nm from the interface. For these values the probability of sticking to the surface is

    P_stick^CP  [%]= V_CP/V_total =  (R^2_exc. - R^2_CP )/  R^2_exc.≃ 11.75 % .

Compared to the lightning-rod effect, CP forces will more significantly affect the omnidirectionality of the trap, as it will result in the particle sticking to the dielectric interface with a higher probability, given that P_stick^CP≥ P_stick^LR. It is important to note that the CP force is inherent to all air-mode cavities, and thus, it is not possible to avoid this source of particles sticking to the surface, only to minimize it by reducing the surface-to-volume ratio near the trapping region. For instance, if one reduced the exclusion region radius as in <ref> the surface-to-volume radius would increase, resulting in even more sticking events due to the CP forces. For example, from the results in <ref>, an exclusion region radius of 20 nm yields a sticking probability of P_stick^CP≃ 90 %, which will approach unity as the exclusion becomes even smaller. Therefore, for deeply sub-wavelength design features, such as the ones in other inverse design dielectric nanocavities <cit.>, the CP forces would entirely dominate the optical trapping, making levitated omnidirectional optical trapping unattainable. In other words, vacuum forces set a minimum length scale on the exclusion region radius (R_exc) for levitated omnidirectional trapping in nanostructures, which for this device is around R_exc≃ 20 nm. For omnidirectionally trapping nanostructures, the minimum length scale set by CP forces can be reduced by increasing the input power, so that optical forces dominate vacuum forces. This would be unfeasible for plasmonic structures, where losses would be enhanced with increased input power, but is only viable for lossless dielectric trapping devices, such as the ones presented in this work.



§ COMPARING THE DEVICE PERFORMANCE WITH STATE-OF-THE-ART DEVICES


Our work constitutes the first dielectric nanostructure design featuring an omnidirectional gradient force trapping potential. To see how its trapping characteristics compare to other state-of-the-art trapping devices, we have gathered and compared the main trap parameters across several platforms in <ref>. Although another work achieved omnidirectional trapping in nanostructures using SIBA forces <cit.>, here we mainly focus on comparing devices that rely on gradient forces to create the optical trap, since this is the operational principle of the optimized device. None of the reference nanostructures can provide sub-wavelength stable trapping in all spatial directions, which is highlighted in the stable trapping directions column in <ref> and with the row coloring. Additionally, the proposed device is an integrated photonic device, meaning that it is waveguide-coupled, in contrast to many other devices that are excited out-of-plane by the use of lasers <cit.>. For these devices, where the laser excitation intensity is reported, we use the laser spot area to compute the total power in <ref>. As noted in the main document, the main comparison proxy we have used to compare optical trapping devices has been the normalized trapping stiffness (η_i), which is calculated as:

    η_i = κ_i  ε_0/α_R  P_in ,

where α_R is the real part of the polarizability of the particle, ε_0 is the vacuum permittivity, κ_i is the trap stiffness for different axes denoted by the index i ∈{ x,y,z }, and P_in is the operating input power of the device. The metric in <ref> measures the trapping stiffness (force per unit distance) normalized by particle size, the contrast between particle and background dielectric permittivity, and input power, making comparisons across platforms possible. Note that this metric metric is not applicable for trapping platforms with no light source, such as vacuum force trapping structures, since there is no input power (P_in) and the normalized trapping stiffness η diverges. The metric could be modified for these devices by redefining <ref> without including the input power.

Using the normalized trapping stiffness as in indicator, we observe that devices relying upon plasmonic resonances can achieve stronger normalized trapping stiffnesses <cit.> than our inverse-designed device. However, plasmonic devices lack omnidirectional trapping, which leads to particles sticking to the surface, and suffer from heating problems due to the optical losses in metals, which can be detrimental in many trapping settings <cit.>. We have highlighted this problem in <ref> by coloring the platforms that lose energy to heating in red and the ones that do not, in blue. Regarding dielectric platforms <cit.>, we see that our platform offers normalized trapping stiffnesses in the same order of magnitude, while being the only proposed solution offering gradient force optical trapping in all spatial directions. As discussed in <ref>, for optical trapping nanostructures, there seems to be a trade-off between achieving stability in the spatial directions on the one hand, and the normalized trapping stiffness achievable by the optical platforms on the other hand. If we compare our platform to optical tweezers <cit.>, which do indeed provide diffraction-limited omnidirectional trapping, we find that they can achieve similar trapping stiffnesses as our platform, but have much lower normalized trapping stiffnesses. This is mainly due to the lack of integration and waveguide coupling of the optical tweezers, which results in larger input power requirements. More importantly, our device operates below the diffraction limit by relying on optical near-field effects, meaning that we can, in principle, achieve higher trapping stiffnesses while relying on a much smaller footprint than conventional optical tweezers. This is also exemplified by comparing the normalized trapping stiffness η in <ref> for our vacuum device (λ=1.55 µm) and water device (λ=775 nm), which for the same particle size and input power, reveal more than an order of magnitude increase of η by operating at smaller scales.






