To correlate the terahertz response with temperature, calibration experiments were carried out by plasmonic heating of the GNR dispersions at various power densities between 0.95 and 60.4 W/cm. The temperature was recorded by inserting a thermocouple probe close to the heating/sensing spot in the cuvette once the temperature reached the steady-state. The change of amplitude following the increase in temperature was observed by monitoring the terahertz echo reflected from the cuvette/sample interface. This is due to the heat-induced change in the refractive index of the water solution at terahertz frequencies. The temperature was then assigned to the corresponding measured terahertz amplitude, resulting in a linear relationship with a coefficient of determinationof ∼0.99, shown in. The slight difference in the calibration curves might be due to the different dimensions of the GNRs, which lead to small variations in the refractive indices. The slopes of the calibration curves were recorded for each dispersion to directly translate terahertz amplitudes into temperature values. Following this postprocessing procedure, a temperature resolution of 1 °C was achieved and a well reproducible measurement can be guaranteed.

A robust temperature measurement relies on an accurate extraction of the amplitude associated with the second terahertz echo, the latter being reflected from the cuvette/sample interface. As shown in, the first echo has an enduring signal or “ringing” due to the ambient water vapor, which interferes with the second echo. Therefore, simply ignoring this interference will lead to an inaccurate measurement of the amplitude associated with the second echo. In order to isolate the first and second echoes, a sparsity-based deconvolution methodwas used to extract the distinct waveforms of the two terahertz reflections, as shown in. Moreover, typical THz-TDS systems suffer from amplitude fluctuations and phase drifts, which can impede reliable and precise measurements.However, the fluctuations in amplitude and phase of the first and second terahertz echoes are correlated since they arise from intrinsic system noise. These effects were mitigated by using a self-reference calibration method.Hence, following the amplitude variations extracted from the first terahertz echo, the second is adjusted accordingly. Therefore, combining the deconvolution and the self-reference calibration method, an accurate evaluation of the second terahertz echo peak value can be obtained. Detailed information about the terahertz sparse deconvolution and the self-calibration method can be found in the supplementary material

Terahertz temperature measurements were carried out based on the THz-TDS system set in reflection mode. A sealed standard 10-mm path length PE cuvette containing 3 ml of each GNR dispersion was placed in the terahertz focusing spot, as depicted in. A continuous-wave NIR laser emitting at 786 nm was employed to excite the sample inside the cuvette. A typical terahertz time trace reflected from the sample is shown in. Three terahertz echoes are clearly identified, which correspond to the terahertz reflections from the air/cuvette (1st echo), cuvette/GNR dispersion interfaces (2nd echo), and a second-round reflection (3rd echo), respectively. The NIR excitation laser spot (∼2 mm diameter) is overlapped with the terahertz focal spot, in order to record the temperature-induced terahertz changes exactly within the heating spot at the interface between the cuvette and the GNR dispersion. The terahertz echo arising from the cuvette/GNR dispersion interface (i.e., the second echo) was utilized to determine the temperature transient during the heating and cooling cycles. The first echo arising from the air/cuvette interface was instead used for self-referencing in the following postprocessing procedure.

B. Quantification of the photothermal effect via THz-TDS

η as well as MHR of the three GNR dispersions, transient temperature measurements were carried out for the heating and cooling cycles by using THz-TDS. During the heating cycle, a NIR illumination with a power intensity of 31.8 W/cm2 was maintained for 12 min. Thereafter, the laser was switched off, and the temperature decay was recorded for 15 min to acquire the temperature transient of the cooling cycle. Each temperature transient was obtained by recording the peak amplitude of the second terahertz echo (after postprocessing) every 10 s and directly translating it into a temperature value, according to the calibration curves. Figure 3 In order to calculateas well asof the three GNR dispersions, transient temperature measurements were carried out for the heating and cooling cycles by using THz-TDS. During the heating cycle, a NIR illumination with a power intensity of 31.8 W/cmwas maintained for 12 min. Thereafter, the laser was switched off, and the temperature decay was recorded for 15 min to acquire the temperature transient of the cooling cycle. Each temperature transient was obtained by recording the peak amplitude of the second terahertz echo (after postprocessing) every 10 s and directly translating it into a temperature value, according to the calibration curves.shows the complete temperature transients of the three GNR dispersions. Their temperature increased under laser illumination and then decreased toward the ambient value after the laser source was switched off, following the trend predicted by Eq. (8) . After heating for 12 min, GNR10 and GNR25 reached a similar steady-state temperature of ∼53 °C, while the temperature of GNR50 increased to ∼45 °C. The fluctuations in the measurement of the temperature transients might be associated with thermal drifts in the system and surrounding environment. The impact of the fluctuations was mitigated by fitting the experimental data with our theoretical model.

η of the GNR dispersions, the rate of energy absorption A and the rate of heat dissipation B were extracted by fitting the experimental temperature transients with Eq. B between the GNRs and the external environment can be determined from the cooling cycle of the temperature transients. The theoretical model for the temperature transient of the heat dissipation can be found by setting A = 0 and T initial = T max in Eq. T(t) = T 0 + (T max − T 0 )e−Bt, where T max is the maximum temperature achieved during photothermal heating. Based on this equation, a linear relationship is identified between ln[(T(t) − T 0 )/(T max − T 0 )] and time t. The heat dissipation rate B can be determined from the slope of this curve. It indicates a first-order decay as the heated system reaches equilibrium with the surrounding medium. By fitting the experimental data with this linear model, the parameter B for each GNR dispersion was determined to be around 3.0 × 10−3 s−1. The fitting results for the cooling curves are shown in Figs. 4(a1)–4(c1) B is retrieved, the energy absorption rate A can be obtained by fitting the heating cycles of the temperature transients with Eq. Figs. 4(a2)–4(c2) R2 is ∼0.97 and indicates an adequate fit. To evaluate the photothermal conversion efficiencyof the GNR dispersions, the rate of energy absorptionand the rate of heat dissipationwere extracted by fitting the experimental temperature transients with Eq. (8) . The rate of heat dissipationbetween the GNRs and the external environment can be determined from the cooling cycle of the temperature transients. The theoretical model for the temperature transient of the heat dissipation can be found by setting= 0 andin Eq. (8) ) =+ (, whereis the maximum temperature achieved during photothermal heating. Based on this equation, a linear relationship is identified between ln[() −)/()] and time. The heat dissipation ratecan be determined from the slope of this curve. It indicates a first-order decay as the heated system reaches equilibrium with the surrounding medium. By fitting the experimental data with this linear model, the parameterfor each GNR dispersion was determined to be around 3.0 × 10. The fitting results for the cooling curves are shown in. The heat dissipation rate does not depend on the GNR dimensions, as shown in Eq. (6) . Once the heat dissipation rateis retrieved, the energy absorption ratecan be obtained by fitting the heating cycles of the temperature transients with Eq. (8) , as shown in. The coefficient of determinationis ∼0.97 and indicates an adequate fit.

Fig. 5(a) Fig. 5 T) at the steady-state with respect to the initial ambient value. The image acquisition is essential to estimate the effective mass of the GNR dispersion. We underline that neither the temperature transients nor their spatial distributions lead to the determination of the photothermal conversion efficiency of the GNRs, but only their mutual correlation. The terahertz thermal images from the front and back window of the cuvette display a similar temperature distribution and clearly reveal that the GNR10 dispersion close to the cuvette bottom was not heated to the same temperature as that at the cuvette top due to heat convection. Therefore, the effective mass of the heated GNR dispersion in the cuvette can be calculated by the following equation: m e ff = ∬ T ( x , y ) d x d y S c m s o l , (11) where T(x, y) is the normalized temperature distribution of the front surface, dx and dy are the dimensions associated with one pixel of the terahertz image, S c is the area of the cuvette surface, and m sol is the entire mass of the GNR dispersion in the cuvette (3.0 g). By taking an average of the values achieved from the front and back windows, the effective mass of the GNR10 dispersion was determined to be 1.74 g. Moreover, the effective mass of the PE cuvette was taken into account as (m eff /m sol ) × m cuv = 0.9 g, where m cuv accounts for the mass of the cuvette surrounding the GNR dispersion. It is important to note that, since no magnetic stirrer was employed in the cuvette, only a fraction of the total mass was affected by the photothermal heating. This led to a nonuniform temperature distribution inside the sealed PE cuvette, which must be accounted for by determining the effectively heated mass of the sample during laser excitation. To this end, terahertz thermal imaging based on raster-scans with spatial resolution of 0.25 mm was performed along the entire front and back sides of the cuvette prior to heating, as well as after having reached the steady-state temperature level. In this case, the sample cuvette was mounted on a motorized X-Y translation stage in order to facilitate the terahertz raster-scans. The NIR beam was fixed to the center of the cuvette [as shown in] during the imaging process (i.e., both the NIR beam and the cuvette were moving across the terahertz beam). The resulting terahertz thermal image maps the temperature distribution at the internal cuvette/GNR dispersion interface when the heating cycle reaches the steady-state. For the GNR10 dispersion, the terahertz thermal images acquired from both the front and back window of the cuvette are shown in, exhibiting the temperature variation (Δ) at the steady-state with respect to the initial ambient value. The image acquisition is essential to estimate the effective mass of the GNR dispersion. We underline that neither the temperature transients nor their spatial distributions lead to the determination of the photothermal conversion efficiency of the GNRs, but only their mutual correlation. The terahertz thermal images from the front and back window of the cuvette display a similar temperature distribution and clearly reveal that the GNR10 dispersion close to the cuvette bottom was not heated to the same temperature as that at the cuvette top due to heat convection. Therefore, the effective mass of the heated GNR dispersion in the cuvette can be calculated by the following equation:where) is the normalized temperature distribution of the front surface,andare the dimensions associated with one pixel of the terahertz image,is the area of the cuvette surface, andis the entire mass of the GNR dispersion in the cuvette (3.0 g). By taking an average of the values achieved from the front and back windows, the effective mass of the GNR10 dispersion was determined to be 1.74 g. Moreover, the effective mass of the PE cuvette was taken into account as () ×= 0.9 g, whereaccounts for the mass of the cuvette surrounding the GNR dispersion.

Fig. 6(a) A, the MHR was also calculated based on Eq. Fig. 6(a) η achieved with these GNRs, especially with the sample containing the particles with the smallest volume [see Fig. 6(a) η of 90%, a value similar to other reports for comparable gold nanostructures. 10 Understanding the photothermal conversion efficiency of gold nanocrystals ,” Small 6, 2272– 2280 (2010). 10. H. Chen, L. Shao, T. Ming, Z. Sun, C. Zhao, B. Yang, and J. Wang, “,” Small, 2272–(2010). https://doi.org/10.1002/smll.201001109 Fig. 6(a) η and MHR follow an opposite trend with respect to the increasing average volume of the GNRs. Whereas η is highest for the GNRs with the smallest diameter, the MHR is maximum for the largest ones. These results follow the behavior that is expected for plasmonic resonators, supporting the reliability of our method and also illustrating the role of the electromagnetic response of the GNRs in the heating process, for which a schematic depiction of some relevant properties can be found in Fig. 6(b) Figure 6(c) σ abs and scattering σ scat ) cross sections corresponding to the GNRs used in the experiments. Importantly, it is the increase in their absorption cross sections that leads to a larger amount of energy absorbed in each particle, thus driving the positive trend for the MHR shown in Fig. 6(a) Therefore, by combining the temperature transients and the thermal images achieved via THz-TDS, the photothermal conversion efficiencies of the three GNR dispersions were calculated, as shown in(red bars), without the need of stirring the sample. Using the known particle concentration (see Table S1 in the supplementary material ) together with the retrieved energy absorption rate, thewas also calculated based on Eq. (10) and illustrated in(orange bars). It is worth noticing the high values ofachieved with these GNRs, especially with the sample containing the particles with the smallest volume [see]. In particular, the GNR10 sample achieves aof 90%, a value similar to other reports for comparable gold nanostructures.Furthermore,clearly shows thatandfollow an opposite trend with respect to the increasing average volume of the GNRs. Whereasis highest for the GNRs with the smallest diameter, theis maximum for the largest ones. These results follow the behavior that is expected for plasmonic resonators, supporting the reliability of our method and also illustrating the role of the electromagnetic response of the GNRs in the heating process, for which a schematic depiction of some relevant properties can be found in. First, the interaction cross sections of a plasmonic nanoparticle do increase with its polarizability, which in turn scales with the volume of the particle.shows the theoretical results for the extinction (sum of absorptionand scattering) cross sections corresponding to the GNRs used in the experiments. Importantly, it is the increase in their absorption cross sections that leads to a larger amount of energy absorbed in each particle, thus driving the positive trend for theshown in. It is relevant to recall that the GNR samples have a comparable optical density, meaning that the GNR50 sample has a significantly lower particle number concentration than the GNR10.