In order to solve equation ( ), a spherical computational domain of radius 15 times the nanoparticle length was used with a perfectly matched layer (shell) specified with a scattering type boundary condition to emulate an infinite medium. Upon finding the absorption cross section or Q from equation , equation ( ) to (3) were solved simultaneously to determine the velocity and temperature fields. To obtain the temperature fields in the case of unconfined nanospheres and nanorods, the convective term in equation ( ) can be neglected and the axial symmetry can be exploited. For the confined cases, axial symmetry is valid only for nanospheres and vertically oriented nanorods (axis parallel to the direction of gravity). For the cases of horizontally oriented (axis perpendicular to the direction of gravity) nanorods in spherical and hemispherical confinements, the computation domain was reduced to 1/4 th of the total by implementing symmetry conditions on the dividing walls (domain independence test, Figure ). For the plasmonic analysis by solving equation ( ), the perfectly matched layer was discretized by prismatic elements in five layers and the physical domain was discretized by tetrahedral elements as shown in Figure (a). Triangular and tetrahedral elements were used to construct the numerical mesh for the 2D axisymmetric and 3D cases respectively (Figure (b) and (c)). The mesh elements were constructed relatively finer close to the walls where no-slip condition is specified. To ensure the accuracy of the numerical results at the expense of optimum computational resources, grid independence tests for each configuration were conducted by comparing local electric field and velocity magnitudes for test grids with varying number of elements, Figure and S7. A tolerance limit of 1% error was set as a criteria for the grid independence test. The BDF (backward difference formulae) time stepping scheme predicts each time step according to the dynamics of the system in order to optimize the computation. Generalized alpha scheme is more robust in that it does not allow time steps larger than a specified value. Upon agreement between the two time-stepping schemes, the BDF scheme was used in order to minimize computational time.