Tailored optical propulsion forces for controlled transport of resonant gold nanoparticles and associated thermal convective fluid flows

Noble metal nanoparticles illuminated at their plasmonic resonance wavelength turn into heat nanosources. This phenomenon has prompted the development of numerous applications in science and technology. Simultaneous optical manipulation of such resonant nanoparticles could certainly extend the functionality and potential applications of optothermal tools. In this article, we experimentally demonstrate optical transport of single and multiple resonant nanoparticles (colloidal gold spheres of radius 200 nm) directed by tailored transverse phase-gradient forces propelling them around a 2D optical trap. We show how the phase-gradient force can be designed to efficiently change the speed of the nanoparticles. We have found that multiple hot nanoparticles assemble in the form of a quasi-stable group whose motion around the laser trap is also controlled by such optical propulsion forces. This assembly experiences a significant increase in the local temperature, which creates an optothermal convective fluid flow dragging tracer particles into the assembly. Thus, the created assembly is a moving heat source controlled by the propulsion force, enabling indirect control of fluid flows as a micro-optofluidic tool. The existence of these flows, probably caused by the temperature-induced Marangoni effect at the liquid water/superheated water interface, is confirmed by tracking free tracer particles migrating towards the assembly. We propose a straightforward method to control the assembly size, and therefore its temperature, by using a nonuniform optical propelling force that induces the splitting or merging of the group of nanoparticles. We envision further development of microscale optofluidic tools based on these achievements.

The optical forces exerted on a colloidal particle is determined by Maxwell's stress tensor and their analysis requires the application of numerical calculation methods. The easily interpreted analytical expression for time-averaged force F, exerted by monochromatic field, is only obtained in the dipole approximation [1][2][3] where E(r) = (E x , E y , E z ) is an electric field vector with E q ≡ E q (r) = |E q (r)| exp[iΨ q (r)], ∇ is a gradient operator expressed in the chosen coordinate system and q is a placeholder for the spatial coordinates x, y, z. Here ε 0 is the permittivity of vacuum, ε m = n 2 m is the relative permittivity of surrounding medium (usually considered as a real parameter) and α = α + iα is the particle polarizability. The first term in Eq. (1) can be rewritten in the vector form as F ∇ = n m α ∇I/2c, where I = n m ε 0 c |E| 2 /2 is the intensity (irradiance) of the incident wave in the medium and c is the speed of light in vacuum. This conservative force, often called intensity-gradient or dipole force, pulls the particle toward the position of the intensity maximum (minimum) if the dispersive (real) part of polarizability α > 0 (α < 0) and allows for the particle confinement in optical tweezers. The second term describes nonconservative forces arising from light absorption and scattering by the particle and depends on * jarmar@fis.ucm.es the dissipative (imaginary) part of the polarizability α . It is determined by phase gradients of field components and takes into account the field polarization [4,5]. Since the α is always positive then this so called scattering force, F scat , propels the particle in the direction of the phase gradients. In particular, in the case of the freestyle optical trap [6] the axial component of the scattering force has to be compensated by the intensity-gradient force to achieve stable three-dimensional (3D) trap. While the transverse component of the scattering force could help in radial confinement and propels the particle along the designed trajectory (e.g., a ring) in the plane transverse to the beam propagation.
Our goal is to consider the application of such propulsion scattering force for versatile motion control of metal NPs. Strictly speaking the expression Eq. (1) is valid for a Rayleigh particle with radius a < 0.1λ, with λ = λ 0 /n m being the light's wavelength in the medium (where the wavevector is k = k 0 n m = 2π/λ). In this case α = σ ext /k = (σ abs + σ scat )/k, where σ ext , σ abs and σ scat are cross sections for extinction, absorption and scattering, correspondingly. As a rule of thumb it is often assumed that the scattering force for a larger particle can also be found as where σ ext depends on the relative permittivity ε p of the particle, its form, size, etc. and can be calculated using Mie theory [7]. The extinction cross section is σ ext = 0.455 µm 2 for the considered gold NP (radius a = 200 nm), at the illumination wavelength λ 0 = 532nm.
In order to study the light-induced motion of resonant metal NPs we have used the socalled freestyle laser trap which allows for both the optical confinement and transport of the particle along arbitrary 3D trajectory with independent control of the optical propulsion force [8,9]. Notice that the scattering optical force is significantly increased when the incident laser wavelength is near the plasmon resonance of the particle, in which case Re [ε p ] ≈ −2ε m holds. This resonance increases the particle speed along the optically defined trajectory, however, stable 3D trapping is not possible due to strong axial component of the scattering force pushing the particle in the beam's propagation direction, which cannot be compensated by the axial component of the intensity gradient force. Nevertheless, stable 2D confinement and transport of a resonant metal NP is possible against a substrate such as the glass cover-slip. We have applied a freestyle laser trap created by a polymorphic beam [10] which is focused in form of 2D laser curve in the plane transverse to the beam propagation allowing the optical transport of the particles along the target curve. Note that in the experiment the laser curve can be slightly defocused onto the substrate. In this work, we are interested in the optical control of the particle dynamics along a fixed trajectory which, as an example, is a circumference (ring) of radius R. Here we have used circular polarization that prevents from an anisotropy in the optical forces responsible for persistent NP velocity oscillations when a linear polarized trapping beam is used [11,12]. Specifically, the circular polarized polymorphic beam associated to the considered ring trap is expressed in Cartesian coordinates (in the input aperture of the objective lens with focal distance f) as the vector where ± = (1, ±i) is the circular polarization vector. This expression can be rewritten in polar coordinates as a superposition of helical Bessel functions The complex weight function allows designing the beam's amplitude (through |g (t)|, in the unit of electric field, V m −1 ) and phase (through Ψ (t)) distributions along the target curve [10].
In paraxial approximation with unlimited objective aperture, the shape of the beam E in the focal plane is described by a ring of radius R as it is easy to prove by calculating the Fourier transform of E 0 (i.e., E = FT(E 0 )) [10]. We will further consider the uniform intensity distribution along the ring that corresponds to |g (t)| = A 0 2πR 2 λf , with A 0 being a constant which depends on the light power of the incident beam at the input aperture of the objective lens. Thus the optical force propelling the particle around the ring trap is described only by the azimuthal component of scattering force, which in considered approximation is directly proportional to the phase derivative Ψ (ϕ) as it follows: where we have introduced the polar coordinates (r, ϕ) in the focal plane. For further analysis it is convenient to express the phase as with S(ϕ) being an arbitrary real function describing the phase distribution along the curve, while the parameter m = Ψ (2π)/2π defines the global phase accumulation along the entire curve, which can be associated with beam generalized topological charge. The azimuthal component of scattering force tangential to the curve (circumference), further refereed to as optical propulsion force, is given by where In the particular case of the uniform ξ-trap −its uniform phase distribution (S(ϕ) = ϕ) yields ξ(ϕ) = m/Rk 0 − the associated polymorphic beam Eq. (4) is reduced to the following helical Bessel beam of order m: and the corresponding optical propulsion force is uniform along the whole curve.
When the helical Bessel beam is focused by a high numerical aperture objective lens, the optical propulsion force can be derived by applying the Eq. (2) and the Richards-Wolf expressions [13] further generalized for the circular polarized vortex beams [14,15]. Thus we where: with β = arcsin(NA/n imm ) and NA being the numerical aperture of the objective, where n imm is the refractive index of immersion medium, whereas symbols ± correspond to the circular polarization vector ± = (1, ±i). The irradiance is given by As we observe from the

SIMULATION METHOD OF 2D LIGHT-DRIVEN MOTION OF COLLOIDAL PARTICLE
The transverse optical propulsion force F = F ϕ (R, ϕ)u ϕ , Eq. (8) with u ϕ being the unit vector tangent to the curve, acting upon the NP induces a dynamics described by the 2D

Langevin equation of motion
where r = r(t) = (x N P (t), y N P (t)) is the position of the NP of mass M , v =ṙ = dr/dt is the speed of the particle, ν is the Stokes drag friction coefficient, and ζ(t) is the stochastic thermal noise term responsible for Brownian motion of the particle [18][19][20][21]. The noise term ζ(t) follows a Gaussian probability distribution such that ζ p (t)ζ q (t ) = 2νk B T δ p,q δ(t − t ), where: the angle brackets denote an average over time, δ p,q is a Kronecker delta function over the coordinate indices and δ(t − t ) is a delta function of time. The drag friction coefficient is given by the expression ν = 6πaη where η is the dynamic viscosity of the medium. Since the considered spherical particle is transported near the glass cover-slip there exists an increase in the hydrodynamic drag approximated by the expression [22] that depends on the particle radius a and the distance h between the particle and the substrate (glass cover-slip) where the particle is trapped.
The equation of motion Eq. (14) can be solved by using a splitting-method time-integration scheme (BAOAB) [18][19][20]: Eq. (14) is split into three parts labelled A, B, and O corresponding to the kinetic term, potential term, and friction-thermal fluctuation term, respectively. In the BAOAB method the integration sequence is [18]: where the index j denotes the temporal dependence as r (j) = r(t), thus r (j+1) = r(t + ∆t), while F (j) = F (r(t)). The friction-thermal fluctuation term has the following parameters: c 1 = exp(−γ∆t), c 2 = (1 − c 1 )/γ and c 3 = k B T (1 − c 2 1 ) with γ = ν/M and β (j) being a white noise function (a set of normally distributed random numbers, ranging from 0 to 1, thus with mean 0 and variance 1). Note that such a noise contribution can be very small if |F| > |ζ(t)|. The integration problem for the non-conservative force field F is studied in detail in [20].
In our case the medium is water and the NP is strongly confined by the laser curve against the glass cover-slip such that h ∼ a can be assumed, see [12]. Therefore, the particle trajectory can be well approximated by the 2D curve's shape R(ϕ) = (R cos ϕ, R sin ϕ). All these facts allow the simplification of Eq. (16) to a 1-dimensional equation motion along the curve R(ϕ). Indeed, the set of expressions Eq. (16) provide the velocity v(ϕ(t)) and position r(t) = R(ϕ(t)) of the particle confined in the curve. In the experimental results the NP experiences a strong radial confinement in the ring trap such |R(ϕ) − R| ≤ 100 nm, which supports the assumed 1-dimensional NP motion along the curve. In the simulation of the light-driven particle motion we have considered a very small time-step value ∆t = 0.1 ms, which is a requirement when the force field F is stiff [18][19][20].
In practice, it is often assumed that the motion of the NP is over-damped,r = 0, thus the equation of motion could be simplified aṡ that in turn allows for the estimation of the force field F (propulsion force) from the measurement of the speed obtained from the particle tracking data.

TEMPERATURE INCREASE OF METAL NPS AND SURROUNDING LIQUID
The temperature increase ∆T N P experienced by a single spherical NP embedded in a homogenous medium (e.g., a fluid) of thermal conductivity κ is given by the expression [23] ∆T N P = Iσ abs 4πκa .
When the NP is located close to the substrate (glass) the thermal conductivity is given as an average κ = (κ w + κ s )/2, where κ w = 0.6 W/(m·K) and κ s = 1.38 W/(m·K) is the thermal conductivity of water and glass, correspondingly, considered in our case. The absorption cross section σ abs = 0.137 µm 2 of the considered gold NP (radius a = 200 nm) has been calculated applying the Mie theory. In the experiment (e.g., in the case of the ξ(ϕ)-trap) we have used an irradiance I = 0.73 mW/µm 2 and therefore ∆T N P ∼ 40 K. The uniform-temperature approximation inside the metal NP can be assumed [24].
Outside the hot NP the steady-state (equilibrium) distribution of temperature increase of the surrounding water can be calculated by using the expression ∆T (r) = ∆T N P a/r [24], see Fig. S3(a), where r > a is the radial distance from the centre of the particle. Since thermal processes are much faster in the metal, the surrounding fluid (water) governs the value of the time scale to reach the steady-state regime, which according to the study carried out in [24] is t = a 2 /D water with D water = 1.43 · 10 −8 m 2 /s being the thermal diffusivity of water. In our case t ∼ 2.8 µs and therefore the equilibrium distribution of temperature increase ∆T (r) is reached almost instantaneously as the NP moves (driven by any phase gradient) around the ring trap. Moreover, due to the low Rayleigh number of water, the thermal-induced fluid convection is not supposed to distort the temperature distribution within the liquid, regardless of the temperature increase or the size of the heat source [25]. Then the spatial distribution of temperature around the moving NP is the same as around a static one.
The dynamic viscosity of the water surrounding the NP varies as a function of the temperature [26] and therefore of the position as η(298 K + ∆T (r)). In Fig. S3(b) the spatial distribution of η(298 K + ∆T (r)) around the considered NP is displayed. The most rapid NP moves with a speed of 20 µm/s therefore in a time ∼ 2.8 µs (required to reach the temperature steady-state regime) the particle only travels a distance of 0.056 nm. Thus, the temperature rise of the surrounding water at this distance can be approximated as ∆T N P = 40 K, see As it has been mentioned in the main text, the local increase of the temperature in the interval ∆T ∼ 70 − 250 K produces a convective fluid flow arising from the temperatureinduced Marangoni effect at the liquid water/superheated water interface due to nanobubbles formation [27]. This fluid flow could explain the motion of tracer NPs toward the G-NP (heat source) observed in our experiments. Since the G-NP is a dynamic assembly of NPs (identical nano-spheres) with varying number of NPs and inter-particle distances an exact calculation of its temperature rise is challenging. Thus, in order to estimate the temperature rise, we simplify the problem considering only the NPs that comprise the G-NP core. They are in close proximity but not attached to each other, in part, because the NPs are coated with charged surfactants that prevent from particle contact and aggregation. It is supported by the fact that NP clustering has not been observed in our experiments. The temperature of the G-NP could be estimated following the expressions given in [23], however, in such a case all the inter-particle distances have to be known. Alternatively, if the inter-particle distances are unknown, the temperature rise ∆T S at the centre of a structure comprising such NPs can be estimated by using the expression derived in [28] for dispersed NPs in a fluid where a S is the radius of the heated region, ρ N is the concentration of hot NPs. In our case, this heated structure corresponds to the core of the G-NP comprising N NPs of radius a = 200 nm grouped in form of disk with radius a S and heigh 2a. Thus, the concentration in the core region of the G-NP is ρ N = N/V with V = 2πa 2 s a being the disk volume. Therefore, the temperature increase ∆T GN P at the centre of the G-NP (at its core) can be estimated from Eq. (19) as it follows which is consistent with the fact of superposition of temperature fields from all the hot NPs.
This collective thermal behaviour can create a homogenous heat source distributed throughout the entire core region of the G-NP [28]. The measured average radius of the disk-like core of the G-NP is a S ∼ 0.6 µm, which gives a number of N ∼ 5 hot NPs comprising the core. Thus, the G-NP transported in the ξ(ϕ)-trap (with I = 0.73 mW/µm 2 ) experiences a temperature rise ∆T GN P ∼ 200 K. Note that the plasmonic coupling between the NPs leads to a red shift of the plasmonic resonance and therefore the temperature rise (∆T N P ) of a NP can be lower in the core than in the case of a single one. However, even if the ∆T N P value is decreased by a factor of two the resulting estimated temperature rise of the G-NP, 100 K, is compatible with appearance of nanobubbles formation [27,29]. On the other hand, the estimated temperature of the superheated water 400 − 500 K surrounding the G-NP is less than the reported threshold for microbubbles formation (580 ± 20 K) [30,31].
In the experiments we have also observed that some NPs are eventually expelled from the G-NP's core when there are more than N ∼ 6 particles. This could be explained by an increased value of the core's temperature rise (i.e., ∆T GN P ∼ 250 K) that leads to a temperature of the superheated water near the threshold for microbubbles formation, thereby creating a stronger water flow able to expel NPs from the core. We think this could be an evidence of a temperature rise limit for the formation of a stable G-NP core. Taking into account the temperature rise limit ∼ 250 K for the nanobubbles regime we infer that a threshold number of NPs comprising a stable core could be N ∼ 250 K/∆T N P , which coincides with N ∼ 6 for ∆T N P ∼ 40 K. Nevertheless, a further study is required to completely confirm this hypothesis.

OPTICAL TRANSPORT ALONG ARBITRARY TRAJECTORIES
In this work we have applied a freestyle trap [8,9] whose shape and propulsion force can be easily adapted to the standing application. As an example, in the main text, we have studied a ring trap with tailored phase-gradient propulsion force to optically transport single and multiple resonant NPs. In this section, we demonstrate that indeed the shape of the optical trap can be arbitrary while preserving the independent control of the ξ(ϕ)-tailored optical propulsion force. In particular, we have considered a knot (eight-shaped) circuit [32] for the optical transport of the gold NPs as shown in Fig. S4. Specifically, the intensity and phase distribution of the trapping beam are displayed in Fig. S4(a), where the phase gradient prescribed along the circuit corresponds to the linear increasing ξ(ϕ)-trap as shown in Fig. S4(b). The corresponding phase-gradient propulsion force drive the transport of the optically transport NPs and G-NP optothermal convertors along complex trajectories. Note that this experiment has been also recorded at 100 frames per second, the time lapse image series has been created from the raw data by using a time lapse interval of 100 ms. Further information about the generation of the applied freestyle laser traps for programmable optical transport can be found in [9,32].