1. Introduction

The pioneering work of Arthur Ashkin [1–4], who discovered optical tweezers, has paved the way for further research and potential applications of the optical manipulation of nanoparticles (NPs). Optical manipulation has crucial applications, particularly in the field of biological systems [5–8], nanotechnologies [9–14], biological sensing [5,6], drug delivery [7,8], and nano-and-microswimmers [13,14]. Qualitatively, the optical forces can be categorized into three categories: (1) the radiative pressure [15–21], (2) the gradient force [15–21], and (3) the spin-curl force [22]. The first 2 components are the most well-known and the former arises due to the scattering of the incident photons by NP, which is along the direction of the propagation of the incident beam, while the latter arises due to the gradient of the optical intensity or phase of the incident light [1,2]. Finally, the spin curl force arises from the curl of the spin angular momentum of the photons [22] due to the non-uniform distribution of the spin density associated with beams that have helical phase fronts, e.g. circular polarization. In the case when the incident beam is linearly polarized, the spin curl force is essentially zero [22]. Thus, the characteristics of optical force on an NP depend on the profile of the incident light (intensity or phase profile along the beam axis or on a plane perpendicular to the beam axis) and the scattering properties of the NP. Several strategies for sculpting the profile of the incident beam have been proposed, such as the solenoid beam [20], the tractor beam [22], and the Bessel beam [22]. With these special beams, manipulations of NPs have been demonstrated, but at relatively short working distances of ∼ 10-30 µm.

A plane wave has a very-long working distance and has a wavefront that is easily accessible in the laboratory scale. The plane wave does not have a gradient in its intensity profile along the beam-propagating axis, i.e. it can only exert scattering forces on an NP. In addition, due to the momentum conservation rule, the plane wave usually results in an optical force on NP in the wave-propagating direction. However, it has recently been found that a single plane wave can exert an optical force on an NP against the wave-propagating direction when the NP is encapsulated with a lower refractive index spherical cavity than the surrounding medium [23–25]. In an experiment, the surface plasmonic resonance (SPR) of the metallic NP has been leveraged to form the nanobubble (NB) to encapsulate the NP (i.e., lower index spherical cavity), and a loosely focused Gaussian beam has been used to mimic the wavefront of a plane wave [20]. The experiment showed that the single beam can pull the metallic NP for a distance of ∼ 0.1 mm against the beam-propagating direction besides the anticipated pushing motion. This finding has opened an opportunity for a plane wave to be used to manipulate NPs in water.

Qualitatively, one can explain that the pulling of a supercavitating NP can be from the reflected stream of the incident photons by the water-NB interface like a macroscale mirror [25,26]. In a more quantitative term, an analysis with electromagnetic multipoles tried to explain that the negative optical force appears when unique multipole interaction terms dominate due to the lower refractive index of the spherical cavity [26]. Recently, a study has illustrated the vector field lines of the optical force on a TiN NP in a nanobubble [27]. This study uses the analytical form of optical force calculated from the coefficients of vector spherical harmonics (VSH). It can visualize the optical force field and has determined the size ranges of NB and NP required to initiate the negative motion of the supercavitating NP. However, these studies assume the first-order scatting by the NP in the NB to estimate the optical force, thus has not fully explored the higher order scattering effects of the optical force on a supercavitating NP. The first-order scattering method neglects any reflected lights toward the NP by the water-NB interfaces in calculating the optical force on NP. Thus, in a situation where the portion of reflected lights increases, such estimated optical force may become inaccurate. For example, those cases include the wavelength regime close to the SPR peak or when the spatial location of the NP is near the water-NB interface. Rings et al. developed a Markov model to explore the hydrodynamics focusing on the kinetic effects arising from the non-uniform temperature and viscosity fluctuations of the Brownian motion of hot nanoparticles under laser illumination [28]. However, the study lacks consideration of the optical forces and the effect on the NP movement under the influence of a cavity.

The supercavitating NP is an optical system consisting of two spheres with different origins (i.e., an eccentric spherical system). To accurately describe the optics of such a system, the multiple events of the reflected light by the NB/water interfaces toward the NP must be considered. It can be achieved by obtaining a completed set of VSH of the NP, which satisfies the boundary conditions of Maxwell's equations at the NB/water interfaces. In the concentric spherical system case, Aden and Kerker [29] have shown that the VSH solution includes all the degrees (i.e., angular modes) from the first-order harmonics (i.e., the first azimuthal mode), as the symmetry leads to the cancellation of the Transverse Magnetic and Transverse Electric modes [30], i.e. the higher order azimuthal modes. However, an eccentric spherical system lacks such a symmetry, and thus all the possible interactions between each mode in the VSH should be considered. Mackowski [31] has provided a robust formulation of scattering between multiple spheres at arbitrary locations, but in this case, the scatterers are outside of each other. Fikioris et al. [32] have presented a general formulation for an eccentrically stratified dielectric sphere but focused only on on-axis translation. Gouesbet et al. [33] have attempted to formulate a generalized Lorentz-Mie theory of an eccentrically located spheroidal NP inside a sphere but under Gaussian beam illumination. Thus, to understand the supercavitating NP, it is necessary to formulate a complete set of VSH of an eccentric spherical system under the incident of a plane wave. It is noted that a finite element method (FEM) can accurately obtain the electromagnetic profiles in such an eccentric spherical system, but it is computationally challenging to fully explore vector field lines of the optical force on the NP inside the NB.

In this study, we have developed a generalized analytical model based on VSH to estimate the electromagnetic field distribution of an eccentric spherical system, under a plane wave illumination, which can calculate the optical force on any arbitrary NP inside a spherical nanocavity surrounded by a medium. As a representative example, we use a solid Au NP encapsulated by a vapor NB surrounded by water and calculate the optical force using the developed analytical model. We have confirmed the accuracy of the developed model using the FEM at conditions where the first-order scatting method fails to accurately calculate the optical forces. We have investigated vector field lines of the optical force on the Au NP and studied its possible moving paths in the NB. We have found that the Au NP can have three kinds of trajectory outcomes: trapped, reaching the light-incoming surface of the NB, or reaching the light-outgoing surface of the NB. We have also investigated the probability of each trajectory outcome at various kinds of system parameters (e.g., the sizes of NP and NB).

2. Theory

In this study, we focus on the optical force on the NP in a spherical nanocavity at the frequency domain. Specifically, we are interested in the possible moving paths of the NP, driven by the optical force, enclosed by the nanocavity surrounded by an arbitrary medium. Due to the strong photo-thermal process, a plasmonic NP suspended in water (medium) can be instantaneously encapsulated by an NB (nanocavity) with the size of O (∼100 nm), which happens at the time scale of ∼0.1 ns [34]. In an NB, the refractive photon stream by the incident plane wave can exert optical force on the NP, and it can direct the NP to move (e.g., toward the surface of the NB) in the bubble. The location of the NP in the NB in turn influences the optical force direction and thus determines if it will result in an optical pushing or pulling motion. For example, if the optical force drives the NP close to the NB surface facing the incident of a plane wave, it can induce a negative force on the NP and thus its backward motion. Therefore, accurately estimating the optical force field is imperative.

For a general analytical formulation, we assume an NP is arbitrarily located in a spherical-shaped cavity (e.g., NB) (see Fig. 1). The particles are spherical objects, optically homogeneous, and the scattered light interacts only with the elements (e.g. NP and NB) defined in the system. The plane wave can exert an optical force on the NP, and the order of magnitude of the optical force in experiments is 10⁻¹² ∼ 10⁻¹⁴ N. In this condition, it is possible to assume that the motion of the NP is limited by the viscosity of the steam (i.e., low Reynold number of ∼ 10⁻³), and the vector field line of the optical force can represent the moving path of the NP in the NB [23,26].

 

Fig. 1. (a) Schematic of the optical configuration of an Au NP enclosed by an NB formed within a medium of water. The separate coordinate systems of the NB and the NP are [X, Y, Z] and [x, y, z], respectively. R_NP and R_NB are the radii of NP and NB, respectively. r₀, θ₀, and φ₀ are the distance, polar angle, and azimuthal angle between the origins of NP and NB coordinate systems, respectively. The plane wave propagation direction is denoted as k_z.

Download Full Size | PDF

The vector field line of the optical force can be calculated using the following recursion relation [27]:

The supercavitating NP is an optical system where a spherical NP is in a spherical NB with different origins (i.e., eccentric spherical system). In this system, the electromagnetic fields can be expressed by VSHs as:

The analytical formula of the optical force ($\textrm{F}$) in Eq. (4) is a result of integrating the time-averaged Maxwell’s stress tensor over the surface of the NP [36,27].

On the other hand, in Eq. (2), to obtain the VSH coefficients, boundary conditions of the tangential electric and magnetic fields [29] need to be used at both the surfaces of the NP (NP/NB interface) and the NB (NB/water interface). The complete set of the tangential boundary conditions of an eccentric spherical system is as follows:

Here, ${Y_{mn}}({\alpha ,\gamma } )$ is the Spherical Harmonics, where α and γ are the end-fire incidence angle and the incident polarization, respectively. The coordinate translation coefficients can be formulated with the translation theorems for VSHs provided by Cruzan [37] or the addition theorem for scalar spherical harmonics [38], and they can be expressed as:

In obtaining Eq. (8), we use the spherical Bessel function of the 1^st kind to interpret the VSHs for the NB with those for the NP, or vice versa. It is also noted that the VSH of the incident plane wave is with the origin at the center of the NB.

The analytical formulation has followed the way suggested by Lorenz - Mie theory. Therefore, in terms of the size of the NP or nanocavity, the validity domain of the multiple-harmonics model is bound to those of the Mie and Lorenz model. Specifically, the multiple-harmonics model uses the vector spherical harmonics (VSH) originally used by Lorenz and Mie. The translation coefficients [37] developed by Cruzan are also derived directly from the VSHs of the generalized Lorenz-Mie Theory. Furthermore, for larger particles above the Mie scattering limit (λ<radius), a large number of harmonics need to be considered which becomes inefficient, and in this range, other methods like ray-tracing, etc., are more desirable.

3. Results and discussions

To check the accuracy of the theoretical model, which we refer to as the ‘multiple-harmonics model, we assume an NP encapsulated by a NB with a refractive index of ${n_{NB}}$=1 (steam) which is surrounded with a refractive index of ${n_{med}}$=1.33 (water), and their optical indices are assumed to be a constant for the wavelength regime (400 nm ∼ 980 nm). The NP-NB system is illuminated by an x-polarized plane wave propagating in the z-direction. We consider an Au NP, and set the sizes of the NB and the NP to be 500 nm and 80 nm, respectively. Our interest is the regime where the NP is close to the NB/water interface. In this regime, the reflections between the NP and the interface become intensive, and thus the higher order harmonics of the VSH are necessary to include such intensive reflections. We selected two wavelength cases, where in Fig. 2(b) the NP-NB system is illuminated by a monochromatic plane wave with λ = 800 nm, which is far away from the SPR peak of the Au NP, and in Fig. 2(c) the NP-NB system is illuminated by a monochromatic plane wave with λ = 560 nm, which corresponds to the SPR of the NP. The wavelength-dependent complex refractive indexes of the Au NP for λ = 800 nm and λ = 560 nm are 0.15352 + 4.9077i and 0.28496 + 2.7390i, respectively [35]. The intensity of the illuminating plane wave considered is taken from the experiment [23], and is 12 mW/µm².

 

Fig. 2. (a) The normalized z-component of the optical force as a function of the number of angular modes (R_NP= 80 nm, R_NB= 500 nm, r₀= 400 nm, θ₀ = 150° and φ₀ = 45°). (b, c) The calculated optical force as a function of φ₀ using the multiple harmonic model, FEM, and the first-order scattering method: (b) the z-component of the optical force for λ = 800 nm and (c) λ = 560 nm (R_NP= 80 nm, R_NB= 500 nm, r₀= 400 nm, θ₀ = 150°, and φ₀ varies from 0°-360°)

Download Full Size | PDF

Practically, one can investigate the threshold number of VSH harmonics by investigating the convergence in the magnitude of optical properties (e.g., optical force or absorption cross-sections). We have investigated F_z as a function of the number of angular modes (n) of VSH at two different NP locations (e.g., r₀ = 100 nm and r₀ = 400 nm). When r₀ is 100 nm, the calculated-F_z converges at n = 4. However, at r₀ = 400 nm, where the Au NP is very close to the surface of NB, it demands higher-order harmonics (n = 16) for the convergence (see Fig. 2(a)).

The optical force components (F_x, F_y, and F_z) are calculated and compared with the results estimated using FEM. We have also included the results estimated by the first-order scattering method [27]. In Figs. 2(b) and 2(c), it is clear that the results from the multiple-harmonics model agrees well with those from FEM for both cases. They have similar optical force values at the examined spatial points, with a ∼1-5% difference, which may originate due to numerical (precision) errors from either of the methods. In the meantime, the first-order scattering model shows the trend of optical forces similar to both the multiple-harmonics model and the FEM, but the amplitude deviate from both models significantly (at a maximum of ∼45%).

It is to be noted the NP can reach temperatures [23] up to ∼ 1000 K due to strong plasmonic photo-thermal heating effect. At high temperatures, the electron-phonon scattering rate increases which in turn makes the Drude damping larger. As a result, the refractive index increases as a function of temperature [38,39]. The change in refractive index can yield a change in the optical force. In the meantime, for the vector field line, the main factor is the normalized directional vector of optical force, as the NP is at a low Reynold number (the direction of optical force at each spatial point is more dominant than the magnitude of optical force itself). As a result, we found that the normalized optical directional vector of optical force is not sensitive to the change of refractive index due to the temperature increment (see Supplementary 1).

To further investigate the validity of the multiple-harmonics model, it is important to review the forces in cases where the resonances are more pronounced. For example, Silicon (Si) NPs have quadrupole resonances in the near-infrared region [40]. In Fig. 3, we have considered a spherical Si NP (${R_{NP}}$ = 220 nm), with a refractive index of 3.48 [40] corresponding to the resonance wavelength of 1105 nm. Here, ${R_{NB}}$ = 500 nm, r₀ = 200 nm. We calculated optical force with the multiple-harmonics model by including the harmonic number up to n = 13, to compare to the forces estimated using FEM. To validate the results, we also considered all the 3 components of the optical forces (F_x, F_y, and F_z) in Fig. 3, where θ₀ has been varied (φ₀ = 45°) and the insets show the φ₀ variation (θ₀ = 150°). The multiple-harmonics model shows good agreement with the forces calculated using FEM, at the examined spatial points with a deviation in the range of 0.25% ∼ 1.4%. However, when the multiple-harmonics model includes the harmonic number up to n = 4 (or n = 1), it has been found that a deviation increases to 3% ∼ 31% (or 43% ∼ 63%). This fortifies the necessity to consider higher orders of scattering in an eccentric Si NP-nanocavity system. It is also noted that, in terms of computational cost, the multiple-harmonics model can calculate an optical force at a spatial point in ∼0.7 sec, which is 1,618 times faster than FEM, showing its efficiency.

 

Fig. 3. (a, b, c) The calculated optical force acting on Si NP (R_NP= 220 nm) inside the NB, with r₀= 200 nm, ${n_{NP}}$ = 3.48 corresponding to λ = 1105 nm as a function of θ₀ and φ₀ using different angular modes (n) of the multiple harmonic model (n = 1,4,13) and FEM. (a) the x-component of the optical force, (b) the y-component of the optical force, and (c) the z-component of the optical force. Here, θ₀ is varied from 0°-180° with φ₀ = 45°, and φ₀ has been varied from 0°-360°, in the inset, with θ₀ = 150°.

Download Full Size | PDF

Before moving on, it is important to mention that Cruzan’s translation coefficient introduces some singularity at higher azimuthal modes (m indexes). This is in the range of 0.3∼1.4% for Au NP (see Fig. 2(c)). We also found that this instability descales with the omission of the complex part of the refractive index as evident from the ∼0.03% discrepancy for the Si NP-NB system (see Fig. 3(c)). Even though the discrepancy is minimal, if need be, our formulism can be decoupled from Cruzan’s translation coefficient, and substituted with any relevant methods of VSHs translation.

In the experiment, it has been observed that the long-range motion of a supercavitating NP is heterogeneous, i.e., some move in the positive direction, some in the negative direction, while some drift slowly. These motions depend on a variety of factors, including the optical and topological heterogeneity of the NPs, the inherently stochastic nature of bubble dynamics (nucleation and growth), and the relative position of NP inside the NB [23]. Among these, the relative position of the NP w.r.t. the NB is the key factor determining the magnitude and direction of the exerting optical force [26]. The previous study [27] with the first-order scattering method has shown that there are three kinds of NP moving paths, which lead to the either NP touching the northern (positive z-direction) or southern hemispheric (negative z-direction) surface of NB or being trapped in the NB. Photo-thermal heating makes the NP hot (>850 k), and as soon as the NP comes in contact with either hemispheres (NB-water interface) it can instantaneously evaporate the water molecules and develop a new NB-water interface. For example, a NP coming in contact with the northern hemisphere under the influence of a positive F­_z, will vaporize and extend the NB-water interface along the beam propagation direction. Again at any moment, the relative position of the NP in the newly extended NB will determine the magnitude and direction of the optical force on the NP. As a result, the long-range motion of the NPs is an accumulation of such stochastic events, calculated at each moment.

In these analyses, if the NP reaches the northern (or southern) hemisphere of NB, the NP may have the positive (or negative) motion along (or against) the incident-beam propagating direction. However, as seen in the comparison studies in Fig. 2, the first-order scattering method can be inaccurate when the NP is closer to the NB/water interfaces, and thus it is meaningful to revisit the moving paths of the NP in the NB with the multiple-harmonics methods.

When the NP is close to the NB-water interface, due to heat transfer, the interfacial vaporized region can have a lower temperature than that near the NP surface on the other side. This temperature difference of vapor region across the NP can lead to thermophoretic forces [41]. The thermophoretic force can potentially impact the trajectory of the NP inside the NB. Thus, it is important to estimate the magnitude of the thermophoretic forces. We found that for our case (Fig. 2(c)), the magnitude of thermophoretic force has the order of 10⁻¹⁵ (see Supplementary 2), while that of optical forces has the order of 10⁻¹². Therefore, the optical force mainly determines the trajectory of the NP. This analysis is well aligned with the previously reported result of Janus NPs (half-dielectric/half-metallic) [41] which are purposely designed to have a strong temperature difference (∼100 K) across the NP, where it is exerted by the thermophoretic force with the order of 10⁻¹⁴ N.

We use the multiple-harmonics methods to study the moving paths of an NP in the NB with the vector field line analysis in Eq. (1) when it is illuminated by a plane wave. For this study, the NP-NB system is illuminated by an x-polarized plane wave propagating in the z-direction with the wavelength of λ = 560 nm, which is the average SPR peak of the NPs with a radius in the range of 20-100 nm. It is noted that the precise SPR wavelength may not alter the conclusion of our study in any significant manner. The northern (or southern) hemisphere of NB is defined by the polar and the azimuthal angle as θ₀ = 0°-90° (or θ₀ = 90°-180°) and φ₀ of 0°-360°, and if the NP, that reaches the NB surface of the northern (or southern) hemisphere, has a positive (or negative) F_z, the case can be considered to induce the positive (or negative) motion of the supercavitating NP. To systematically investigate the vector field lines, we set the starting points as the equal-gridded points (36 points) of θ₀ and φ₀ spaces at a distance r₀ = 300 nm (∼ 0.6R_NB).

Figure 4(a) shows the representative vector field lines of the NP with R_NP = 55 nm and R_NB = 500 nm, and it shows that there are two distinctive characteristics of the lines. One is the vector fields line ending at the northern surface of the NB, which can induce the positive motion of the NP, and we refer to these lines as the ‘positive motion lines’. The other lines finish at a point on the central axis of the NB (i.e., the z-axis), meaning the NP on this vector field line can be trapped in the NB (it is referred to as ‘trapped line’). The trapped point is where all the components of the optical forces are ‘zero’ and this point can be termed the ‘nodal point’. Also, we count the number of the positive motion lines and the trapped lines to approximate their probabilities. In this case, the probability that the NP shows a positive motion line is ∼ 41%, while the NP being trapped at has a probability of ∼ 59%.

 

Fig. 4. (a and b) The vector field lines of the NP trajectories starting from different θ₀: (a) R_NP= 55 nm, R_NB= 500 nm, r_o= 300 nm, φ₀= 45° (b) R_NP= 30 nm, R_NB= 500 nm, r_o= 300 nm with, φ₀= 45° (c and d) The probabilities of the three different trajectory outcomes (trapped, positive motion, negative motion) as a function of R_NP: (c) R_NB= 500 nm and (d) R_NB= 250 nm.

Download Full Size | PDF

In the meantime, for the 500 nm NB, we also investigate the vector field lines for a smaller NP with R_NP = 30 nm (Fig. 4(b)). In this case, we can see another kind of vector field line that can bring the NP to the southern surface of the NB, which can lead to the negative motion of the supercavitating NP against the incident beam propagation (referred to as ‘negative motion line’). The probability of this negative motion line is ∼58%. Interestingly, the probability of the trapped lines is ∼ 30%, and that for the positive line is ∼15%. In comparison to the case with the larger R_NP = 55 nm, the probability of positive lines is reduced significantly. This result agrees with the previous theoretical study [27] with the first-order scattering method, which shows that the smaller NP in a given size of NB can have negative motion.

We further investigate the NP size effect on the probabilities of those lines in the 500 nm NB (Fig. 4(c)). For the NP with radius smaller than 35 nm, it has a higher probability of achieving negative motion than positive motion. In the meantime, for the NP with radius larger than 35 nm, it is seen that the NP has a higher probability of taking the positive motion lines. Also, the probability of trapped lines stays around ∼37% for 20 nm < R_NP < 45 nm and rapidly drops to zero at R_NP > 55 nm. On the other hand, if the radius of the NB is reduced to 250 nm, the trapped lines are not found at all examined sizes of NP (20 nm < R_NP < 100 nm, see Fig. 4(d)). Also, the probability curves of the negative or positive motion lines, reveal some obvious characteristics. For example, the reduction of the NP size leads to an increment of the portion of negative motion lines, and at R_NP ∼20 nm, our prediction indicates that there are only negative motion lines. Conversely, the positive motion lines increase with a larger size of NP, reaching 100% at R_NP > 50 nm.

These results indicate that the probability of positive or negative motion or being trapped strongly depends on the size of NB and NP. To further study the NP-in-NB system for the trajectory outcome, we investigate the regime of R_NP of 20-100 nm and R_NB of 150-500 nm. We have considered r₀= 0.6R_NB for the starting position of the NP and 18-equal-gridded starting points. For some cases where the NP starts at a close physical proximity to the NB surface, it is excluded since we imposed the condition R_NP+ r₀+ 20 nm ≥ R_NB. These excluded regions are painted ‘white’ in Fig. 5. Figure 5(a) shows the probability distribution of the NP getting trapped for the aforementioned parameter space. The trapping event is apparent when R_NB is relatively large (∼300-500 nm) and the R_NP is relatively small (20-60 nm), with the highest probability of ∼60%. The probability distribution of the positive motion lines has been shown in Fig. 5(b). This case has 100% probability when the size of the NP is relatively bigger (> 60 nm) all across the parametric space and when it is 45-60 nm with R_NB< 300 nm. Figure 5(c) shows the probability distribution of the negative motion lines. This case shows that R_NP of 20-45 nm can have probabilities from ∼40-100%. The probability increases to 100% with smaller NB (R_NB ∼160-275 nm).

 

Fig. 5. The probability contour map of the NP for three different trajectory outcomes (trapped, positive motion, negative motion) as a function of R_NB and R_NP (a) trapped (b) positive motion (c) negative motion. Here, the R_NB is varied from 150 nm to 500 nm, and the R_NP is varied from 20 nm to 100 nm. In all the plots, the white triangle indicates exclusion of the points from calculation, due to the starting point of the NP being in close physical contact with the NB, i.e., R_NP+ r₀+ 20 nm ≥ R_NB, where r₀ = 0.6 R_NB.

Download Full Size | PDF

4. Conclusion

In summary, we have developed the generalized multiple-harmonics model based on VSH for the NP-in-nanocavity system under plane wave illumination. The multiple-harmonics model is computationally efficient and accurately estimates the optical force on an NP inside an NB surrounded by water. The optical forces calculated by our model match well with that of the FEM, even when the NP is very close to the surface of NB, where the first-order scattering method is not accurate. The study of vector field lines shows that the solid Au NP in an NB can be trapped in certain locations inside the NB or has light-induced positive or negative motions. The probability of each trajectory outcome as a function of R_NP and R_NB is studied. The results show that an Au NP larger than ∼60 nm prefer to have positive motions when enclosed by a supercavitation. In the meantime, one should use an Au NP smaller than ∼40 nm to enable the negative motion if the NB size is smaller than ∼270 nm, which may be achievable by adjusting the fluence of the laser that excite the NB through the plasmonic heating effect. Our fundamental study may provide useful guidance on designing various NP manipulation applications.