Nonlinear plasmonics at high temperatures

1 Introduction

Nanoplasmonic systems have been intensively studied in recent decades due to their unique potential for local field enhancement and subwavelength confinement and are considered as promising candidates for a wide variety of applications [1], [2]. However, the inherent absorption in the metal proves to be a substantial obstacle toward the realization of real-life applications.

Accordingly, in recent years, the applied plasmonic research focused on applications that exploit the absorption in the metal as a means to generate heat on the nanoscale [3], [4] – a research topic usually referred to as thermoplasmonics. This resulted in a wide range of emerging applications, at different ranges of temperatures, starting from photothermal (PT) imaging [5], [6], through cancer treatment [7], temperature measurement [8], plasmonic photovoltaics [9] and water boiling, sanitation, and superheating [10], [11], [12], [13], up to thermophotovoltaics [14], [15], diffusive switching [16], radiative heat transfer [17], plasmon-mediated photocatalysis [18], [19], [20], plasmon-assisted chemical vapor deposition [21], and heat-assisted magnetic recording [22], which may involve temperatures even higher than 2000 K.

In the majority of works in thermoplasmonics, the optical and thermal properties are assumed to be fixed. However, as the heat is induced by laser illumination (unlike external heating assumed in thermal emission engineering and nanoscale radiative heat transfer [23], [24]), it is necessary to account for the coupling between the electromagnetic fields, the temperature, the optical properties, and the thermal properties (i.e. heat capacity, thermal conductivity, and Kapitza resistance) of the constituents to achieve a quantitative understanding of the field and temperature distribution. To the best of our knowledge, such a systematic, self-consistent study was not performed so far in the context of thermoplasmonics. Specifically, the temperature dependence of metal permittivity was accounted for in some studies either through a (cubic) thermal nonlinearity [25], [26] or, more generally, based on a combination of the two-temperature model [27] with a complex model of the permittivity. Such models have to correctly account for a rather large number of competing effects within the metal, including electron scattering, thermal expansion, band shifting, the effect of Fermi smearing on intraband transitions [25] and interband transitions [28], [29], [30], and more. To the best of our knowledge, such a comprehensive study was done only by Stoll et al. [31]. Moreover, during intense illumination, all the effects mentioned above are modified due to the deviation of the electron distribution function, scattering rates, etc., from their equilibrium values – effects that were studied only partially [25], [32], [33], [34]. Most importantly, the studies such as [25], [26], [27], [28], [29], [30], [31], [32], [33], [34], and the many less complete ones, were all dedicated to the ultrafast regime, such that efforts were made to avoid the longer-term (i.e. few picoseconds and longer) thermal effects, for example, by looking at nanoparticles (NPs) not larger than a few nanometers in diameter, lowering the repetition rate, etc. (see discussion in [35]), to the extent that the study of these longer-term thermal effects themselves was neglected. Thus, the relative weight of the above-mentioned effects was not studied so far for continuous wave (CW) illumination. However, the study of the CW limit becomes interesting again, with the growing interest in thermoplasmonic applications.

In the few studies dedicated to longer-term thermal effects (i.e. for CW illumination) [36], [37], the model used for the temperature dependence of metal permittivity neglected some of the dominant physical effects, most prominently, the temperature dependence of the interband transitions. Even in the mid-infrared regime, where no intense sources are available, it did not include all dominant effects known from the literature. Thus, the quantitative predictions in these studies are questionable.

Finally, the temperature dependence of the thermal properties was not taken into account before, either for ultrafast or CW illumination to the best of our knowledge.

To close this knowledge gap, in this article, we perform a thorough theoretical study of the high-temperature regime of nanoplasmonic systems under intense optical illumination in the visible range and for CW illumination. Our study focuses on the (classical) interaction between the temperature, permittivity, and electromagnetic fields. We use experimentally measured data for the various optical and thermal properties to avoid the need to dwell on the details of the underlying physics, which, as explained, is only partially understood; however, where possible, we try to identify the relevant physical mechanisms by comparing the theory to experimental results. Specifically, we will show that the thermo-optical nonlinearity can be very strong, thus, allowing us to explain the experimental observations of the strong nonlinear scattering from metal NPs [38], [39], [40], [41]. More generally, this study should serve as the starting point for further experimental and theoretical studies of the underlying physics, of other regime of parameters (specifically, of pulsed illumination, different materials, geometries, etc.), and enable a quantitative study of the various applications mentioned above as well as several others such as nonlinear composites/metamaterials [35], [42], [43], [44], optical limiting [45], [46], [47], [48], plasmon lasing [49], [50], [51], and superresolution techniques based on metal NPs [6], [38], [39], [40], [41], [52], [53], [54], [55], [56].

The paper is organized as follows. We start by explaining why the temperature dependence of the optical and thermal properties is usually neglected and identify cases where the temperature dependence and mutual coupling of the Maxwell and heat equations is nonnegligible. We then solve the Maxwell and heat equations self-consistently for small metal spheres illuminated by intense visible light and elucidate the large errors in the calculations of the temperature and field distributions associated with neglecting the temperature dependence of the gold permittivity both in-resonance and off-resonance. We then show that the temperature dependence of some additional parameters, such as thermal conductivity and Kapitza resistance, is also required for a correct quantitative prediction of the temperature and field distributions. Finally, we discuss the implications of our results to a previous experimental work and specify several future measurements necessary for further studies of the strong temperature nonlinearity of metals.

2 Self-consistent calculation of temperature in metal nanostructures

We want to calculate the scattering of an incident CW off some metal-dielectric nanostructure as a function of pumping intensity and/or temperature. This requires us to understand how much does the metal temperature increase under this illumination and how much, in turn, does this temperature increase affect metal permittivity (hence the electromagnetic field distribution around the nanostructure).

In the simplest model, the metal-dielectric system is assigned a single, spatially nonuniform temperature T (i.e. we neglect the difference between the electron and lattice temperatures). Then, under CW illumination (and with no temporal pump modulation), the heat equation governing the temperature dynamics reduces to the Poisson equation:

where κ is the thermal conductivity (specifically, κ_m and κ_host for the metal nanostructure and the dielectric host, respectively). Note that, in principle, thermal conductivity can be temperature dependent. The typical boundary conditions accompanying Eq. (1) are the continuity of the temperature T and heat flux κ∇T across the interface between different materials.^[1]

The heat source p_abs represents the density of absorbed power (in units of W/cm³). Classically, the relation between the absorbed power and the local incident electromagnetic field intensity is given by

where E→ is total (local) electric field (namely, the solution of the vectorial Helmholtz equation), ω is the pump frequency, and ϵm=ϵ′m+iϵ″m is the complex (relative) permittivity of the metal, which serves as the heat source in this problem. This expression is sometimes replaced by ~αI (or ~σ_absI), where α is the absorption coefficient (absorption cross-section) and I is the (local) beam intensity. Quantum mechanically, p_abs has to be calculated as the spectral integral over all the possible transitions from electron levels and over the frequency content of the pump pulse times the respective photon energy, involving both absorption and emission to and from each level [33].

The common model in heat calculations of plasmonic systems (frequently referred to as thermoplasmonics [4]) is to solve the Maxwell equations first for ambient conditions, that is, assuming ϵ_m(ω; T=T_env), where T_env is the temperature far away from the heat-generating (metal) objects. Then, one substitutes the resulting electric field distribution into the heat source term p_abs (2) in the heat equation (1). Below, we refer to this approach as the temperature-independent permittivity (TIP) model.

The TIP model is appropriate as long as the relative change of the permittivity Δϵ_m~ΔTdϵ_m/dT is small. Typically, the thermoderivative dϵ_m/dT varies between ~10⁻⁵/K for standard dielectric materials [58] and up to 10⁻⁴–10⁻³/K for water [59] or metals [60], [61]. Thus, as long as the temperature increase (with respect to T_env) is modest (i.e. limited to a few degrees), the change of the permittivity is indeed negligible. Potentially, the opposite signs of the thermoderivatives of the dielectric material and metal may cause the overall temperature dependence of the system under consideration to be weaker than in each of the constituents [61], thus, providing further justification for treating the permittivity as temperature independent. The permittivity changes may be negligible also in the wavelength regime for which dϵ_m/dT vanishes. Peculiarly, it turns out that, for gold, this regime is about 520–550 nm [60], [61] (i.e. it coincides with the plasmon resonance wavelength of small metal nanospheres, which have been subject to extensive study) [4], [31], [38], [39], [40], [61].

However, in plasmonic nanostructures under intenseillumination (as for thermophotovoltaics [14], [15], plasmon-mediated photocatalysis [18], [19], [20], plasmon-assisted chemical vapor deposition [21], and heat-assisted magnetic recording [22]), the conditions prescribed above are typically not fulfilled. Indeed, whereas the (relative) modification of the real part of metal permittivity due to the changes of the temperature may be small, a substantial increase of the temperature (a few tens of degrees or more) may cause the imaginary part of metal permittivity ϵ″m, to change substantially.

In the case of external heating (and weak illumination), one has to use the appropriate permittivity data for the ambient temperature and solve only the Maxwell equations, as done routinely for room-temperature studies. In contrast, intense (laser) illumination will result in mutual coupling of the heat and Maxwell equations via p_abs (2), requiring them to be solved simultaneously. In these cases, the standard model described above (TIP model), which does not take into account the thermo-optical (nonlinear) response to the electromagnetic field, would have to be replaced with a temperature-dependent permittivity (TDP) model. This is essential to make the results from such high temperature applications quantitatively relevant.

Remarkably, it is a common practice to take into account the thermal nonlinearities of the host medium (e.g. for PT imaging [5], [6], [62], cancer treatment [7], and thermal lensing [63]). However, the majority of studies within the plasmonics community ignore the temperature dependence of metal permittivity. Some of the earlier studies did account for the thermal response of the metal by approximating it with a cubic nonlinearity (see [25], or [26] for a more recent review). This approach was used, however, only for cases where the pump pulse was not longer than a few nanoseconds and in the perturbative regime, where the relative permittivity changes were small such that the cubic approximation is sufficient. Moreover, these studies focused primarily on the electric field distribution and ignored the temperature itself. Similarly, studies of effective medium theories applied to media with (thermal) cubic nonlinearities also focused on the field rather than the temperature distribution (see, e.g. [25], [35], [42], [43], [44], [45], [46], [47], [48], [64] and references therein). One of the reasons for that is obviously that measuring the temperature in the near field of the NPs remains a very difficult task despite the progress made recently [65], [66].

Within a TDP model, we expect to be able to distinguish between two scenarios. In the general scenario, as the temperature (hence the imaginary part of metal permittivity) changes, the heat generation rate p_abs changes as well. In particular, if ϵ″m grows with temperature, the TIP model will provide an underestimation of the actual temperature, as could be calculated from the fully coupled (TDP) model. The field distribution, on the contrary, may differ only slightly from the field distribution predicted by the simplified TIP model since typically, ϵ″m ≪ |ϵ′m|.

In contrast, at plasmon resonance, the metal nanostructure acts as a cavity whose quality factor scales inversely with the imaginary part of metal permittivity ϵ″m. Accordingly, if ϵ″m increases with temperature, the local electric field (hence the heat power dissipation) drops, the resonance broadens and the overall power dissipation decreases with respect to the prediction of the TIP model. Accordingly, the temperature will rise more slowly, such that the simplified TIP model will provide an overestimate of the actual temperature. In this case, the TDP model will also predict a substantial change of the scattered field with respect to the prediction of the TIP model. The opposite will happen if ϵ″m decreases with temperature.

These two effects will be demonstrated analytically and numerically in Section 3. However, in the meantime, from this discussion, it is obvious that there is a strong spectral sensitivity such that the solution would almost never follow the predictions of the TIP model. Below, we demonstrate the differences between the TIP and TDP models in some specific examples, showing that they can be substantial for realistic cases and for many applications that are studied extensively these days.

To quantify these differences, one has to have available comprehensive data of the temperature dependence of metal permittivity. However, quite surprisingly, such data hardly exist, even for gold, which is the plasmonic material studied most extensively (see detailed discussion in [67]). In the absence of elaborate experimental data, theoretical models for the temperature dynamics [27], [32], [68] and metal permittivity dynamics [31] were developed. However, as mentioned above, effectively all these studies focused on the ultrafast (up to a few picoseconds) regime, and only a few of these studies accounted for all the relevant physical mechanisms [31]. Similarly, the multitude of models where the thermal response is approximated as a cubic nonlinearity [25], [26] did not consider the thermal response on time scales longer than a few nanoseconds. The quantitative predictions in the few studies of the CW nonlinear response should, as mentioned, be taken with a grain of salt due to the missing ingredients in the permittivity models employed.

Thus, to the best of our knowledge, there is no complete model for the slow thermal response as appropriate for CW illumination. In this regime, the electronic response that dominates the ultrafast response becomes negligible, and other effects, such as lattice heating and thermal expansion [31], stress and strain, band shifting [69], and indirect (i.e. phonon-assisted) interband transitions [70], take dominance.

To close this knowledge gap, we have recently performed ellipsometry measurements to retrieve the permittivity data of bulk gold at increasing temperatures [67]. Our study showed that ϵ″m increases substantially with temperature across the visible spectral range. Indeed, in the temperature regime of 300–570 K, Figure 1 shows changes of up to ~25%–30% in the visible range for selected wavelengths. In the near-infrared regime, an increasingly stronger dependence on temperature is observed [67].

Figure 1:

Au ellipsometry data.

(A) Real and (B) imaginary parts of the (relative) permittivity extracted from the ellipsometry measurements of an annealed (blue) and unannealed (red) Au film at λ=533 nm for 300–570 K. (C and D) Same data for λ=671 nm.

Our study also shows that the changes to the real part of the gold permittivity are substantially smaller with respect to the room temperature values (about 1%–3% in the temperature regime studied here). Similar results appear in two recent independent studies [71], [72] as well as for Ag [71], [73].

To simplify the modeling and discussion, we assume in what follows that the host (dielectric) material permittivity is purely real and nondispersive. This assumption has a negligible effect on our results. Indeed, the numerical examples below show that the changes of the host permittivity and of the real part of metal permittivity have a secondary effect on the temperature and field distribution. This residual temperature dependence will have to be taken into account in applications of PT imaging [5], [6], [62] and treatment [7], water boiling [11], [12], plasmonic (thermo)photovoltaics [14], [15], thermal lensing [63], and plasmon-assisted catalysis [18], [19], [20], where the level of temperature rise of the surrounding medium is critical. This is, however, left for future studies.

3 Metal spheres

Although, in general, the problem at hand requires a self-consistent numerical iteration scheme involving both the Maxwell and heat equations [36], [37], for some simple geometries, one can avoid solving the Maxwell equations and rely on known solutions. As a generic example, we now consider the temperature of a single small (with radius a ≪λ) spherical metal NP illuminated by a plane wave. Hence, the quasi-electrostatic solution of Maxwell equations [74] holds such that the electric field inside the NP (E→NP≡E→(r<a)) is uniform and given by

where E→inc is uniform (i.e. a fixed parameter). In what follows, we suppress the vector symbol. Using Eq. (3) in Eq. (2) gives a uniform power dissipation:

The solution of the Poisson equation (1) for this case is [4]

where κ_m and κ_host, the thermal conductivities of the NP and host, respectively, are assumed for the moment to be temperature independent (hence uniform) and Pabs(ω, T)≡∫pabsdV is the total power dissipated in the NP.

Because typically κ_m≫κ_host, it follows from Eq. (5) that heat diffusion is sufficiently strong to homogenize the temperature within the NP; this assumption is supported by exact simulations, showing temperature uniformity even for much larger NPs [37]. Thus, neglecting the small temperature variation, we can define T_NP≡T(r<a) so that P_abs=4πa³p_abs/3; by Eq. (5), we then get

Substituting Eq. (4) in Eq. (6) gives

where ϵ′tot≡ϵ′m+2ϵd.

Equation (7) is a simple root equation for the NP temperature that is easy to solve. However, before presenting detailed numerical examples, we discuss several general properties of the solution.

In the general (off-resonance) case, the real parts of the permittivities do not perfectly cancel, such that, typically, ϵ′tot≫ϵ″m. Then, we get

Indeed, as predicted, we see from Eq. (8) that the absorbed power (hence the overall temperature) will be higher in the TDP model compared to the TIP model.

On the contrary, at resonance, the real part of the denominator vanishes such that upon a temperature increase, the power dissipation (4) drops for increasing ϵ″m. This is due to the resonant nature of the interaction – ϵ″m is inversely proportional to the quality factor of this effective resonator. At the same time, the resonance broadens and the temperature itself is given by

We thus see from Eq. (9) that at resonance the absorbed power (hence the overall temperature) will be lower in the TDP model compared to the TIP model.

From the above discussion, the reasons for neglecting the changes of the real part of the metal (and dielectric) permittivities become apparent. Indeed, the changes of ϵ′tot are relatively small and cause only a slight shift of the plasmon resonance position. It is clear, however, that the important parameter is the shift of the real part of ϵ′tot rather than the shift of any of its constituents, as various combinations of thermoderivatives of the metal and dielectric will give rise to shifts in different directions.

Once the NP temperature is determined, one can calculate the scattered field using the quasi-static solution [74]. In the case of a single intense (pump) beam, the scattered field E→sc≡E→(r>a)−E→inc, is given by

When the intense beam is accompanied by a second, weaker (probe) beam, the scattering of the probe will be given by the same expression, where the permittivities and fields are evaluated at the probe frequency. Because this case is effectively similar to the standard linear case or to the case of external heating [67], [75], it will not be considered further.

3.1 Numerical examples

Based on the experimental data for annealed gold [67], as appropriate for metal NPs made by the pulsed laser ablation of gold films [76], we initially solve Eq. (7) for λ=533 nm (permittivity data given in Figure 1A and B). Figure 2A and B shows that, when the system is tuned away from resonance (the host permittivity is ϵ_host=5.5), for sufficiently large pumping intensity, the naive (TIP) model indeed underestimates the temperature rise in the particle. For example, we see that the TDP model predicts T=594 K, whereas the TIP gives T=552 K (i.e. an error of ~17% of the temperature increase). This error is commensurate with the corresponding change of ϵ″m, which increased by ~30% (the real part changes by ~1%). Nevertheless, the TIP model still predicts the scattered field with reasonable accuracy (Figure 2C). Indeed, in this case, the denominator for the expression for the scattering [see Eq. (10)] is only slightly affected by the change of the imaginary part of the permittivity.

Figure 2:

Temperature and intensity for intense off-resonant illumination.

(A) Calculated temperature for λ=533 nm and the off-resonance case (ϵ_d=5.5) for the TDP model (blue dots) and TIP model (solid blue line) based on annealed permittivity data. (B) Temperature difference between both models. (C) Peak intensity of the scattered field as a function of the incoming pump intensity for the two models.

At resonance, on the contrary (λ=533 nm but with ϵ_d=2.25), the naive (TIP) model overestimates the temperature rise in the NP (see Figure 3A and B). For example, when the TDP model predicts T=594 K, the TIP gives T=694 K (i.e. ~34% error in temperature rise measurement). More importantly, in this case, the TDP model predicts a 40% decrease of the scattering (Figure 3C). Figure 3 also shows the results based on nonannealed permittivity data as appropriate for metal NPs synthesized in solution [67]. One can see that, although the results are qualitatively similar, the nonannealed gold shows a much stronger sensitivity to the rising temperature. This emphasizes the need to account for the relevant permittivity data depending on the metal particle preparation method [67].

Figure 3:

Temperature and intensity for intense on-resonant illumination.

Same as in Figure 2 for the on-resonance case (obtained by setting ϵ_d=2.25). Also shown in (a) are the results for unannealed data (TDP model, red dots; TIP model, solid red line), and (b) shows the corresponding difference between the two models. Finally, (c) shows also the (normalized) experimentally measured data of scattering from a single 40 nm Au NP embedded in index matching oil under CW illumination.

Most importantly, Figure 3C also shows a comparison to measured scattering data^[2] from a single 40 nm Au NP embedded in index matching oil under CW illumination. One clearly observes very good qualitative agreement between the theory and the measurement, achieved without any fitting parameters. This agreement between the scattered fields also reveals the NP temperature – the lowest scattering levels are attained for a temperature rise of only a few hundred degrees (i.e. well below the melting temperature of the Au NPs), which is somewhat less than 1000 K [13]. Note that a quantitative agreement of the scattered field and temperature requires a further refinement of our modeling (see Section 5).

As a comparison, we show in Figure 4 the temperature and scattered field for resonant illumination at λ=671 nm (ϵ_d=2.5); the permittivity data at this wavelength are shown in Figure 1C and D. Although the trends are qualitatively similar to the case of λ=533 nm, the nonlinear response is stronger – the scattered field drops by 60% and the temperature error is up to ~250 K.

Figure 4:

Temperature and intensity for intense on-resonant illumination (λ=671 nm).

Same as in Figure 2 for the on-resonance case (λ=671 nm and ϵ_d=2.5).

Finally, we note that, when the variation of the real part with temperature is taken into account [i.e. when we solve Eqs. (7) and (10) for ϵ′tot(T)≠0], there is no substantial change of any of the results described so far. This happens because the (absolute, as well as) relative changes of the real part of metal permittivity with temperature are typically smaller than those of the imaginary part; accordingly, they have far smaller influence on the temperature and scattering.

4 Additional considerations

Below, we discuss two additional aspects of the thermo-optical problem at hand that, to the best of our knowledge, are discussed for the first time in the current context.

4.1 Role of the temperature dependence of thermal conductivity

So far, we assumed that the thermal conductivities are temperature independent. However, the temperature dependence of the thermal conductivity is well known for a wide range of materials. Remarkably, its variation with temperature is comparable to that of metal permittivity. For example, the thermal conductivity of water increases by about 10% between 300 and 400 K; beyond this temperature, the water boils. Oil exhibits comparable changes over a wider temperature range, with some oils exhibiting increased conductivity with growing temperature and some exhibiting reduced conductivity. The thermal conductivity of other materials, such as collagen [77], quartz, silicon wires, or aluminum oxide, exhibit even stronger temperature dependence. The thermal conductivity of the metal itself also varies substantially with the temperature; however, because it is typically much larger than the host conductivity, this variation plays a negligible role for our purposes [see Eq. (5)]. Thus, it is clear that this dependence has to be taken into account to accurately determine the temperature and field distributions. In general, if the host thermal conductivity increases with temperature, then the temperature rise is lower than that predicted by a model that ignores this effect and vice versa.

The exact solution (5) used so far will not hold anymore for a temperature (hence space)-dependent thermal conductivity. However, exploiting again the uniformity of the temperature inside the NP allows us to keep using the implicit relation (7). Numerical simulations (see Figure 5) show that the error associated with the change of thermal conductivity with the temperature (taken as T_env/κ_envdκ/dT~±10%) is of the same order of the temperature change itself. As expected, a similar trend is found also for the off-resonant case (not shown).

Figure 5:

Effect of the temperature dependence of the host thermal conductivity on NP temperature.

Calculated temperature for on-resonance case (λ=533 nm) for the TIP model (solid red line) and TDP model with temperature-independent thermal conductivity (blue dots), both as in Figure 3, compared to a thermal conductivity that increases with temperature (green dots) and thermal conductivity that decreases with temperature (black dots).

In the latter two cases, thermal conductivity changes by as much as 13%.

5 Discussion

The results shown above raise a clear need to take into account the temperature dependence of the optical and thermal properties of the metal (and its surroundings) in calculations of field and temperature under intense illumination conditions. In particular, the errors associated with the neglect of the temperature dependence of these quantities grow monotonically with the temperature rise and can reach several tens or even hundreds of degrees for the refractory applications (i.e. even up to 100% relative errors); for resonant illumination, there are comparable relative errors in the scattered fields. In fact, for some applications, such as PT imaging [5], [6], [62], correcting errors of even a few percent could be substantial.

More generally, our calculations provide a complete treatment of nonlinear plasmonic systems at the high-temperature regime that goes beyond the perturbative description of the thermo-optical response [25], [26], [60]. Indeed, we intentionally avoid any assumption on the functional dependence of the metal permittivity on the temperature or intensity (e.g. an assumption of a cubic nonlinear response [25], of a constant thermoderivative dϵ_m/dT [44], [60], or of an averaged response in the effective medium spirit [44]). This approach allowed us to identify the thermo-optical mechanism as being responsible for the nonlinear scattering of monochromatic waves from Au NPs that was observed experimentally [38], [39], [40], [41], showing deviations from the linear prediction (TIP) of several tens of percent (see Figure 3). Indeed, such changes of scattering are shown to be commensurate with the change of the imaginary part of metal permittivity with the temperature (see Figure 1). Remarkably, we confirm that the effect occurs on a subwavelength scale, from a single NP and potentially its immediate surrounding (via the thermal conductivity), rather than being an effect accumulated on macroscopic distances or due to interparticle interactions or aggregation, as one may conclude from previous studies of NP suspensions (see, e.g. [43], [45], [46], [47], [48], [80]).

In contrast to previous works, which relied on a theoretical model that missed some dominant physical effects [36], [37], [75], our study relies on experimentally-measuredpermittivity data [67] and focuses on the visible range. Furthermore, we show stronger effects from NPs smaller than those studied before. Yet, it is important to note that our model provides only a qualitative match to the experimental data. A quantitative agreement requires accounting for the actual size of the particles (i.e. to go beyond the quasi-static approximation employed here, as done in [37]) and for the temperature dependence of the thermal properties of the metal and host, which is currently not known. For completeness, it is also desired to develop a theoretical model for the (slow) thermal nonlinearity of gold to support the experimental results. In contrast to the common models (used in some previous studies [36], [37], [75]), a complete theoretical model will have to account for the temperature dependence of metal permittivity on both intraband and interband transitions and specifically for the effects of the temperature on the NP volume, electron scattering rates, electron distribution (Fermi smearing), lattice spacing (band shifting), and stress/strain build-up, as well as for nonequilibrium effects and multiphoton absorption, which we neglected. A model that describes the interplay and relative importance of these effects is yet to be developed. Such a model will be also particularly important to explain the nonlinear scattering under pulsed illumination, which typically involves higher intensities than for CW (up to GW/cm²), and exhibited opposite trends to those observed for CW illumination [45], [46], [47], [48], [80].

In the same vein, we should mention that the current analysis of the thermal effects may not be sufficient to address the complete intensity dependence of the scattering from metal NPs. Indeed, it was shown [38], [40] that, for sufficiently high excitation intensities, the decrease of the scattering changes to a sharp increase, occasionally and somewhat confusingly, referred to as “reverse saturation”.^[3] This effect may be related to electron population redistribution due to Fermi smearing (i.e. based on the distribution of thermalized electrons), which shows a rather complicated and nonintuitive spectral dependence with several spectral regimes where the permittivity decreases upon heating [28], [29], [30], [31]. Alternatively, the increased scattering may be related to absorption saturation (i.e. based on the distribution of nonthermal electrons) [25], [33] or to an effect associated with the host [e.g. (nonlinear) absorption and phase/structural change]. The determination of its origin also awaits the comprehensive permittivity model and thus left for a future study.

In that regard, we emphasize that the use of the permittivity data under external heating in laser illumination calculations (as adopted in the current study or in [31]) is justified only if the effects associated with nonthermal electrons, which accompany intense illumination, are negligible compared to the effects associated with thermalized electrons. This seems to be the case for gold [26] – the absorption saturation due to interband transitions, which is related to partial population inversion, is predicted [25] and experimentally verified [43], [64] to be smaller than the nonlinearity associated with heating (thermalized electrons, Fermi smearing) at least for moderately high intensities and pulsed illumination. A simple estimate based on the measured cubic nonlinearity, χAu(3)′​′∼5⋅10−9cm2/W at λ=532 nm [25], [43], [64], together with the associated field enhancement within the gold, |3ϵd/ϵ″m|2∼20 (see Eq. (3) and [81]), shows that the permittivity change for I_inc=1 MW/cm² induces a ~10% change of the imaginary part of the permittivity, as indeed we have observed experimentally (see Figures 1 and 3).^[4]

Having said that, we emphasize that the estimates above, which are based on the measurements of the ultrafast thermal response, are only partially appropriate for the current context of a CW illumination. Indeed, the ultrafast thermal nonlinearity was derived in [25] by neglecting the diffusion of heat from the NP to its surroundings (see also [58]). In our configuration, however, heat diffusion is clearly important [see, e.g. Eq. (7)], so that the overall thermal response depends also on the host properties as well as on the particle size. This may give rise to different values of the nonlinearity. In general, though, as already noted above, a complete quantitative match of the model to the experimental data will have to be deferred to a future study.

Finally, we hope that our study would motivate further studies of thermo-optical nonlinearities at the high-temperature regime of other gold NPs as well as similar studies of other metals, especially those proposed for use in refractory plasmonics applications [22].