Guillaume Baffou
email: guillaume.baffou@fresnel.fr

Pierre Bon

Julien Savatier

Julien Polleux

Min Zhu

Marine Merlin

Hervé Rigneault

Serge Monneret

Thermal Imaging of Nanostructures by Quantitative Optical Phase Analysis

Keywords: Microscopy, thermodynamics, phase imaging, plasmonics, Gold nanoparticles

We introduce an optical microscopy technique aiming at characterizing the heat generation arising from nanostructures, in a comprehensive and quantitative manner. Namely, it permits to i) map the temperature distribution around the source of heat, ii) map the heat power density delivered by the source and iii) retrieve the absolute absorption cross section of light-absorbing structures. The technique is based on the measure of the thermal-induced refractive index variation of the medium surrounding the source of heat. The measurement is achieved using an association of a regular CCD camera along with a modified Hartmann diffraction grating. Such a simple association makes this technique straightforward to implement on any conventional microscope with its native broadband illumination conditions.

We illustrate this technique on gold nanoparticles illuminated at their plasmonic resonance. The spatial resolution of this technique is diffraction limited and temperature variations weaker than 1 K can be

Probing temperature at the nanoscale is fundamentally a complicated task mainly because heat has a non-propagative nature, unlike light. Yet, the ability to investigate thermal-induced phenomena on the nanoscale would help developing unprecedented applications in many active current areas of research, [START_REF] Govorov | Generating Heat with Metal Nanoparticles[END_REF] namely photothermal cancer therapy, [START_REF] Lal | Nanoshell-Enabled Photothermal Cancer Therapy: Impending Clinical Impact[END_REF][START_REF] Jain | Au Nanoparticles Target Cancer[END_REF] drug release, [START_REF] Han | Drug and Gene Delivery using Gold Nanoparticles[END_REF] nano-surgery, [START_REF] Urban | Single Step Injection of Gold Nanoparticles through Phospholipid Membranes[END_REF] opto-acoustic imaging, [START_REF] Eghtedari | High Intensity of in vivo Detection of Gold Nanorods Using a Laser Optoacoustic Imaging System[END_REF][START_REF] Beard | Biomedical Photoacoustic Imaging[END_REF] nanochemistry, [START_REF] Jin | Localized Temperature and Chemical Reaction Control in Nanoscale Space by Nanowire Array[END_REF][START_REF] Nitzan | Theoretical Model for Enhanced Photochemistry on Rough Surfaces[END_REF] micro-fluidics, [START_REF] Donner | Plasmon-Assisted Optofluidics[END_REF] magnetic recording [START_REF] Challener | Heat-Assisted Magnetic Recording by a Near-Field Transducer with Efficient Optical Energy Transfer[END_REF] and thermonics. [START_REF] Wang | Thermal Memory: A Storage of Phononic Information[END_REF] Twenty years ago, as a temperature microscopy technique, the use a sharp hybrid tip as a nano-thermocouple was proposed. [START_REF] Pollock | Micro-Thermal Analysis: Techniques and Applications[END_REF][START_REF] Sadat | Nanoscale Thermometry using Point Contact Thermocouples[END_REF] Even though such an approach achieved a spatial resolution of 50 nm, it appeared to be very invasive due to the proximity of the tip to sample. More recently, several thermal microscopy technique have been developed that are much less invasive since they are based on far-field optical measurements. These optical techniques rely either on fluorescence intensity, [START_REF] Löw | High-Spatial-Resolution Surface-Temperature Mapping Using Fluorescent Thermometry[END_REF][START_REF] Ross | Temperature Measurement in Microfluidic Systems Using a Temperature-Dependent Fluorescent Dye[END_REF] fluorescence polarization anisotropy, [START_REF] Baffou | Temperature Mapping Near Plasmonic Nanostructures Using Fluorescence Polarization Anisotropy[END_REF][START_REF] Baffou | Mapping Heat Origin in Plasmonics Structures[END_REF] fluorescence spectroscopy, [START_REF] Vetrone | Temperature Sensing Using Fluorescent Nanothermometers[END_REF][START_REF] Li | Single Quantum Dots as Local Temperature Markers[END_REF][START_REF] Aiguouy | Scanning Thermal Imaging of Microelectronic Circuits with a Fluorescent Nanoprobe[END_REF] Raman spectroscopy, [START_REF] Dang | Femtosecond stimulated Raman Scattering Picosecond Molecular Thermometry in Condensed Phases[END_REF] infrared spectroscopy [START_REF] Sade | Scanning Near Field Infrared Radiometry for Thermal Imaging of Infrared Emitters with Subwavelength Resolution[END_REF] or X-ray spectroscopy [START_REF] Van De Broek | Temperature Determination of Resonantly Excited Plasmonic Branched Gold Nanoparticles by X-ray Absorption Spectroscopy[END_REF] measurements. However, all these techniques suffer from drawbacks that make them suited only for particular cases. Namely, they suffer from a slow read-out rate, [START_REF] Vetrone | Temperature Sensing Using Fluorescent Nanothermometers[END_REF][START_REF] Li | Single Quantum Dots as Local Temperature Markers[END_REF][START_REF] Aiguouy | Scanning Thermal Imaging of Microelectronic Circuits with a Fluorescent Nanoprobe[END_REF][START_REF] Dang | Femtosecond stimulated Raman Scattering Picosecond Molecular Thermometry in Condensed Phases[END_REF] a low resolution, [START_REF] Sade | Scanning Near Field Infrared Radiometry for Thermal Imaging of Infrared Emitters with Subwavelength Resolution[END_REF] a lack of reliability, [START_REF] Löw | High-Spatial-Resolution Surface-Temperature Mapping Using Fluorescent Thermometry[END_REF][START_REF] Ross | Temperature Measurement in Microfluidic Systems Using a Temperature-Dependent Fluorescent Dye[END_REF] a poor temperature sensitivity [START_REF] Dang | Femtosecond stimulated Raman Scattering Picosecond Molecular Thermometry in Condensed Phases[END_REF] or the necessity to modify (tag) the medium with fluorescent molecules or nanocrystals. [START_REF] Löw | High-Spatial-Resolution Surface-Temperature Mapping Using Fluorescent Thermometry[END_REF][START_REF] Ross | Temperature Measurement in Microfluidic Systems Using a Temperature-Dependent Fluorescent Dye[END_REF][START_REF] Baffou | Temperature Mapping Near Plasmonic Nanostructures Using Fluorescence Polarization Anisotropy[END_REF][START_REF] Baffou | Mapping Heat Origin in Plasmonics Structures[END_REF][START_REF] Vetrone | Temperature Sensing Using Fluorescent Nanothermometers[END_REF][START_REF] Li | Single Quantum Dots as Local Temperature Markers[END_REF][START_REF] Aiguouy | Scanning Thermal Imaging of Microelectronic Circuits with a Fluorescent Nanoprobe[END_REF] In this paper, we propose a thermal microscopy technique sensitive to the thermal-induced variation of the refractive index of the medium of interest. This optical technique circumvents the limitations of all the previous techniques mentioned above and is straightforward to implement on any conventional optical microscope. We explain how it leads to a quantitative real-time retrieval of the distributions of both the temperature and the source of heat with a sub-micrometric resolution. The theory and the retrieval algorithms are detailed. We chose to illustrate the capabilities of this technique on gold nanoparticles (NPs), which act as ideal nanosources of heat when illuminated at their plasmonic resonance frequency. Moreover, we show how a quantitative measurement of the actual absorption cross section can be achieved whatever the nature of the absorbing structure. We chose to name this technique TIQSI for Thermal Imaging using Quadriwave Shearing Interferometry.

Results and discussion

We consider in this work a two-dimensional distribution of absorbers (gold nanoparticles) distributed at the interface between a solid (glass substrate) and a liquid environment (like water). When the absorbers are illuminated, heat is generated and a temperature steady state profile appears over a time scale that depends on the spatial extension of the source of heat (typically below 10 µs for structures smaller than 1 µm, in water 25 ). This temperature increase is responsible for a variation of the refractive index of the surrounding liquid over a typical length much larger than the size of the particle itself. This effect has been recently exploited independently by the groups of Lounis and Orrit as a localization technique of NPs or single molecules. [START_REF] Boyer | Photothermal Imaging of Nanometer-Sized Metal Particles among Scatterers[END_REF][START_REF] Berciaud | Photothermal Heterodyne Imaging of Individual Nonfluorescent Nanoclusters and Nanocrystals[END_REF][START_REF] Cognet | Photothermal Methods for Single Nonluminescent Nano-Objects[END_REF][START_REF] Gaiduk | Room-Temperature Detection of a Single Molecule's Absorption by Photothermal Contrast[END_REF] NPs and molecules of interest in these works were non-luminescent and far too small to scatter any incident light, making them invisible using any conventional optical technique. To circumvent this limitation, the idea developed in these groups was to exploit the light absorption capabilities of the nano-objects to thermally induce an extended variation of the refractive index of the surroundings, much easier to detect than the nano-objects themselves. However, this localization technique is based on the modulation of the heating laser to achieve an heterodyne detection. For this reason, it cannot be used to measure any steady-state temperature distribution. Moreover, this technique remains slow and requires a complex experimental configuration comprising an acousto-optical modulator and a heterodyne detection.

Recently, our group has developed an optical phase microscopy technique capable of mapping quantitatively the optical path difference (OPD) through a refractive object in real time with a sensitivity of ∼ 1 nm. [START_REF] Bon | Quadriwave Lateral Shearing Interferometry for Quantitative Phase Microscopy of Living Cells[END_REF] This technique was illustrated by imaging microscopic dielectric beads and living organisms. The experimental approach is based on the use of a quadriwave lateral shearing interferometer as a wave front analyzer (WFA) (see Supplementary Information for more details). Basically, it consists of a modified Hartmann grating (MHG) associated with a regular CCD camera. The device, named Sid4Bio ® , was purchased from Phasics SA. 31 The interferogram produced by the MHG and recorded by the CCD camera can be processed in real time to retrieve the OPD. Interestingly, this interferometric technique is achromatic and can be mounted on any conventional white-light transmission microscope. In the present work, we demonstrate how this experimental approach can be the basis of an efficient and easy-to-implement thermal imaging technique. The experimental setup is detailed in Figure 1. All the samples we used consist of a sandwich structure where a liquid layer of thickness h is limited by two glass coverslips in the z direction. The gold nanoparticles acting as nanosources of heat are lying on the bottom coverslip.

When experiencing a temperature variation δ T (r) = T (r) -T ∞ , a given medium undergoes a variation The sample is illuminated using a Köhler configuration with a reduced numerical aperture (D2) to increase the spatial coherence of the light. A laser is used to excite the gold structures at their plasmonic resonance and heat them (either at 532 or 808 nm, depending of the sample). It is focused at the objective entrance pupil by L1 to perform a wide field illumination. A diaphragm D3 is used to control the beam diameter and to achieve a uniform illumination at the sample location. This diaphragm is optically conjugated with the sample by the objective/L1 association The sample is imaged onto a CCD camera by L2 through a modified Hartmann grating. This creates an interferogram on the CCD that is processed to retrieve the optical path difference generated by the thermalinduced refractive index variation.

δ n(r) of its refractive index. This variation is usually described by a Taylor development of the N th order: 32

δ n = N ∑ j=1 β j δ T j . ( 1
)
where β j are empirical parameters (see Supplementary information or Ref. [START_REF] Thormählen | Refractive Index of Water and its Dependence on Wavelength, Temperature and Density[END_REF] for values). For the sake of simplicity, we shall consider in this theoretical part only the first (linear) term of this Taylor series:

δ n = β 1 δ T .
In water, this assumption is valid if the temperature increase over the ambient temperature does not exceed ∼ 20 • C. [START_REF] Thormählen | Refractive Index of Water and its Dependence on Wavelength, Temperature and Density[END_REF] It makes the equations linear and simplifies the algorithm used to retrieve the temperature profile. When higher temperature increases are considered, a non-linear iterative numerical procedure is required. The approach is more complex but still simple to implement as detailed in supplementary information.

To begin with, let us consider a point-like source of heat located at the solid-liquid planar interface at (x, y, z) = (0, 0, 0). In the surrounding liquid, the temperature distribution is governed by the Poisson equation:

κ ∇ 2 T (r) = P 0 δ (r) (2)
where κ is the thermal conductivity of the surrounding liquid (0.6 W/m/K for water), T (r) the temperature and P 0 the power absorbed (i.e. delivered) by the point-like source of heat. In spherical coordinates, the solution reads:

T (r) = T ∞ + P 0 4π κ r = T ∞ + P0 r , (3)
T (r) = T ∞ + P0 G T (r) (4)
where T ∞ is the ambient temperature, P0 = P 0 /4πκ is the normalized power and G T (r) = 1/r the standard Green's function associated to the Poisson equation.

Experimentally, we use the WFA to measure two phase images of a given area: One reference image acquired (once far all) without heating and one image under heating. The subtraction of these two images provides a measure of the spatial distribution of the optical path difference δ ℓ(x, y) through the liquid layer specifically due to thermal-induced variation of the refractive index distribution δ n(x, y, z). It can be written as:

δ ℓ(x, y) = h 0 δ n(r)dz = h 0 β 1 δ T (r)dz. (5)
According to Eq. ( 3), the profile of the temperature variation δ T reads in radial coordinates (ρ, z):

δ T (ρ, z) = P0 ρ 2 + z 2 . ( 6
)
Using Eqs. ( 5) and ( 6) and after integration, the OPD reads:

δ ℓ(ρ) = P0 β 1 ln h/ρ + 1 + (h/ρ) 2 , (7)
δ ℓ(ρ) = P0 β 1 sinh -1 (h/ρ) = P0 G ℓ (ρ), (8)
where G ℓ (ρ) = β 1 sinh -1 (h/ρ) stands for the Green's function for the phase distribution, since the source of heat is a Dirac distribution. By eliminating the radial coordinate ρ in Eqs. ( 6) and ( 7), one ends up with a bijective relation between the optical phase difference and the temperature increase at z = 0:

δ T (ρ, 0) = P0 h sinh δ ℓ(ρ) P0 β 1 . (9)
Note that this formula is only valid for a point-like source of heat, since it was derived using formula (6).

In the most general case (extended and non-uniform source of heat), such a bijective relation does not exist between the OPD and the temperature and an inversion procedure is required, which we shall detail now.

Consider now an extended non-uniform source of heat characterized by a delivered heat power density (HPD) p(x, y) = 4πκ p(x, y). The subsequent OPD and temperature distributions can be expressed as convolutions with the associated Green's functions:

δ ℓ(x, y) = [ p ⊗ G ℓ ](x, y), ( 10
) δ T (x, y) = [ p ⊗ G T ](x, y), (11)
which turns into simple multiplications in the Fourier domain:

δ ℓ(k) = p(k) Gℓ (k), (12)
δ

T (k) = p(k) GT (k), (13)
where f stands for the Fourier transform of f . From the experimental phase image δ ℓ(x, y), the heat source distribution p(x, y) can be simply retrieved using Eq. ( 12): p(k) = δ ℓ(k)/ Gℓ (k) followed by an inverse Fourier transform. However, such a basic approach is likely to cause divergence problems in case of low signal to noise ratio in the OPD images. This standard deconvolution issue is well-known and different approaches have been proposed in the literature. Among them, one of the simplest is the Tikhonov regularization method, [START_REF] Tikhonov | Solution of Ill-posed problems[END_REF] which we used here. Its principle is detailed in supplementary information. Then, the temperature distribution T (x, y) is retrieved using Eqs. (13) followed by an inverse Fourier transform.

The spatial resolution of the TIQSI technique is given by λ /NA (since the illumination is spatially coherent) where λ is the characteristic wavelength of the Köhler illumination and NA the numerical aperture of the objective. Our set of measurements corresponds to λ = 650 nm and NA = 1.4, which leads to a theoretical spatial resolution of around 450 nm. The frame rate of our CCD camera is 10 Hz (100 ms per 11). e) Measured radial profile of the OPD distribution in image b (solid line), fitted using Eq. (9) (dashed line). f) Delivered power P 0 , as function of the laser irradiance, measured by fitting the OPD radial profile using Eq. ( 9) with P 0 and h as the unknown parameters. image). We usually average a series of 10 to 30 frames to improve the signal to noise ratio. Hence, one usually needs a few seconds per image. This is much faster than most of the other thermal microscopy techniques that need a few minutes or even a few hours to record an image. This readout rate is not a fundamental limitation of the TIQSI technique and it could be improved with a faster and better cooled CCD camera.

We shall first present and discuss experimental results obtained on a point-like source of heat, i.e. a gold disc, 300 nm in diameter and 40 nm thick illuminated by a laser beam at λ 0 = 808 nm (Figure 2a).

Figure 2b displays the raw image of the OPD measured by the WFA. From this image, we have applied the Tikhonov algorithm to retrieve the actual HPD and temperature distributions, as displayed in Figure 2cd. Interestingly, the temperature extension (Figure 2d) is much narrower that the extension of the OPD distribution. This observation is consistent with theln(ρ) decrease of the OPD Green's function (Eq. ( 7)), compared to the much steeper 1/ρ decrease of the temperature Green's function (Eq. ( 4)). In this measurement, 30 frames of 100 ms were averaged. The standard deviation of the noise in the raw OPD image (Figure 2b) is 0.15 nm and in the temperature image (Figure 2d), 1.1 K.

In experimental optics, determining the actual absorption cross section of a unique absorber is not straightforward mainly because the absorption process is precisely not associated with any light reemission that could contain information and be detected in the far-field (unlike extinction or scattering cross-sections measurements). Optical techniques to measure absorption cross sections have been proposed recently. [START_REF] Arbouet | Direct Measurement of the Single-Metal-Cluster Optical Absorption[END_REF][START_REF] Celebrano | Single-Molecule Imaging by Optical Absorption[END_REF] However, all of them suffer from at least one of these drawbacks: i) They consist in fact in measuring extinction instead of absorption. This assimilation is valid only for small particles, i.e. when scattering is negligible (extinction=absorption+scattering), ii) They do perform absorption spectroscopy but they do not yield a quantitative measurement of the cross section. In this second case, a quantitative ACS can be eventually estimated indirectly by comparing the photothermal signal with the one obtained on a reference absorber with a supposedly known ACS. [START_REF] Gaiduk | Room-Temperature Detection of a Single Molecule's Absorption by Photothermal Contrast[END_REF] Our approach circumvent these two limitations at the same time: It allows to get a direct quantitative measurement of the actual absorption cross section σ abs without any assumption and independently of the nature of the absorber. One just needs to estimate the absorbed power P 0 and then calculate σ abs = P 0 /I where I is the irradiance of the laser beam (power per unit area).

There are two possibilities to estimate P 0 . The first one consists in fitting the radial profile of the OPD image (as shown in Figure 2e) using formula (7) with P 0 (and eventually h) as the unknown parameter.

This first approach gives in the case of Figure 2: σ abs = 14930 ± 530 nm 2 . The second approach consists in estimating P 0 by a spatial integration of the HPD image (Figure 2c). It naturally yields a very similar result: σ abs = 15300 ± 600 nm 2 . Calculations using the Green's dyadic tensor method 18 on a gold disc, 288 nm in diameter, 39 nm thick, discretized using 2066 cells and illuminated at lambda=808 nm, led to σ abs = 17200 nm 2 , which is the good order of magnitude. The difference of 14% with the experimental measurement can be mainly explained by the fact that the real structure is not a perfect disc. The accuracy of the measurement can be optimized by acquiring a set of images at different irradiances, as illustrated in Figure 2f. The first approach is handy since it does not require the computation of the HPD distribution, while the second approach is more general and applies even for extended (not point-like) sources of heat.

Note that one does not need to optically resolve the NP to measure its ACS. Indeed, the ACS retrieval procedure is based on the fit of the logarithmic OPD profile around the particle that extends to infinity, even for nanometric/point-like particles (see Figure 2a-c). We shall now present experimental results that address the situation of an extended distribution of heat sources. In Figure 3, experiments have been carried out on an array of quasi-hexagonally organized NPs.

The sample has been done by block-copolymer lithography. [START_REF] Aydin | Micro-Nanostructured Protein Arrays: A Tool for Geometrically Controlled Ligand Presentation[END_REF] It consists of spherical gold NPs, d = 30 ±3 nm in diameter and separated by p = 100 ± 10 nm. The NP interdistance is small enough compared to the resolution of the microscope to consider the HPD as a continuous physical quantity. The heat is gen-erated by a laser beam at λ 0 = 532 nm, close to the resonance wavelength of the gold particles. In the first case, the NP array covers the whole sample (Figure 3a) and the laser illuminates uniformly a circular area of diameter 8.0 µm. In the second case, NP distributions form cross-like structures (Figure 3f). One of them is uniformly illuminated by the heating laser. From the optical path difference measured by the WFA (Figure 3b-g) the distribution of the HPD is computed using the Tikhonov inversion algorithm. As expected, the heat source distributions (Figure 3c-h) are well delimited while the temperature distributions tend to spread out (Figure 3d-i). In order to verify the consistency of the temperature distribution, we performed numerical simulation of the temperature profile in Figure 3d using the Green's function formalism 25 (Figure 3e). The input parameters regarding the geometry and the illumination conditions were I = 105 µW/µm 2 , D = 8.0 µm, p = 105 nm, d = 28 nm. The only variable parameters were d and p and remain very close to the measured values on the SEM image. The matching is very good as observed in Figure 3e, which reinforces the validity of our post-processing numerical model. A last result that we wish to share is presented in Figure 3j. It shows that, as expected, the temperature increase is proportional to the laser irradiance, even over a wide range of temperature. The non-linearities of the problem are thus well taken into account in our numerical procedure (see Supplementary information).

In this last paragraph, we wish to discuss the temperature sensitivity that the TIQSI method can achieve. Figure 4 presents a measurement performed on an extended quasi-hexagonal array of nanodots. The structure is illuminated by a uniform circular beam of diameter 11 µm and and power P = 0.35 mW. Such a low irradiance (I = 3.7 µW/µm 2 ) led to a fairly noisy OPD image (Figure 4a) due to a weak temperature increase. Processing of the OPD image led to a temperature image presented in Figure 4b and revealed that the temperature increase was about 1 K. One can see that such a temperature increase is perfectly visible on the OPD image, even though it leads to a very noisy temperature image. The standard deviation of the noise in the OPD image (Figure 4a) is 0.044 nm and in the temperature image (Figure 4b), 0.36 K. The higher noise level in the temperature image compared with the OPD image is inherent to the deconvolution procedure, and mostly features low frequencies. Below a temperature increase of 1 K, it is hard to obtain a well-defined temperature map. However, since most of the noise observed on the interferogram recorded by the CCD image is due to dark current, it could be possible to increase the signal to noise ratio and reach even better temperature sensitivity by increasing the number of accumulated frames, or by using a better cooled CCD camera. Note that no physical limitation exists on the temperature accuracy, unlike the spatial resolution that is bound to remain diffraction-limited.

Conclusion

In summary, we introduced a thermal microscopy technique sensitive to thermal-induced refractive index variation. The measurement is achieved using a simple experimental configuration consisting of a Hartmann diffraction grating attached to a CCD camera. This technique cumulates many valuable advantages compared with previous thermal imaging techniques: i) It is fast. According to the CCD frame rate, it takes around one second to acquire an acceptable image in most cases. ii) No scanning is required.

It is a one-shot measurement of a 2D image. iii) It has a diffraction limited resolution (less than 500 nm). iv) No modification of the sample is required, like using fluorescent markers. It works whatever the sample as far as there is a liquid environment. v) No need to use sophisticated devices such as heterodyne detection, acousto-optic modulator, spectrometer, etc, like previous thermal imaging techniques. All the thermal measurements can be performed by plugging in the CCD/grating device on a conventional optical microscope, just like a regular CCD camera. vi) The TIQSI techniques features an un precedented versatility since it can quantitatively measure temperature, heat source density and absorption cross sections of nanoparticles and nanostructures. We chose to illustrate the abilities of the technique on gold NPs, which stand for ideal nano-sources of heat. However, we wish to emphasize that the method is not restricted to metal NPs and can be extended to any 2D-distribution of heat sources, independently on their nature. Moreover, this technique is straightforward to implement on any conventional microscope. For these reasons, pending important applications are envisioned not only in plasmonics but also in areas such as nano-electronics, micro-fluidics or cell biology.

Methods

Fabrication of the quasi-hexagonal arrays of gold nanodots: In a typical synthesis, polystyrene(1056)-block-poly(2-vinylpyridine)(495) (PS1056-b-P2VP495) from Polymer Source Inc. was dissolved at room temperature in anhydrous toluene (Sigma-Aldrich) with a concentration of 5 mg/mL and stirred for 2 days. The quantity of gold precursor was calculated relative to the number of P2VP units with a loading parameter equal to 0.5, i.e., 1 molecule of HAuCl 4 for 2 P2VP monomers. Hydrogen tetrachloroaurate (III) trihydrate (HAuCl 4 .3H 2 O, Sigma-Aldrich) was added to the dBCP solution and stirred for 2 days in a sealed glass vessel. Glass cover slips (Carl Roth) were cleaned in a piranha solution for at least 5 hours and were extensively rinsed with MilliQ water and dried under a stream of nitrogen. Micellar monolayers were prepared by dip-coating a glass cover slip into the previously prepared solutions with a constant velocity equal to 24 mm/min. To remove the organic template and to form inorganic nanoparticles, the dip-coated glass slides were exposed to oxygen plasma (150 W, 0.15 mbar, 45 min, PVA TEPLA 100 Plasma System). In the case of the micropatterned arrays of gold nanodots, a 1-µm thick layer of AR-P5350 photoresist (Allresist GmbH) was spin-coated onto the arrays. Illumination was performed under a maskaligner (Suss Microtec GmbH) with a dose of 175 mJ/cm 2 from a HBO 350 mercury lamp and the microstructures were developed according to the manufacturer's instructions. [START_REF] Aydin | Micro-Nanostructured Protein Arrays: A Tool for Geometrically Controlled Ligand Presentation[END_REF] The enlargement of gold nanodots was performed in an aqueous solution containing hydroxylamine (0.02 mM) and HAuCl 4 (0.1 wt%) for 10 min. Scanning electron measurements were performed with a Dual BeamTM (FIB/SEM) instrument (Quanta 3D FEG, FEI, Hillsboro).

Acknowledgment:

We acknowledge fruitful discussions with Philippe Réfrégier. We thank the Institute of Photonic Sciences (ICFO) for providing the lithographic sample (Figure 2a) . We acknowledge financial support from the region Provence-Alpes-Cote d'Azur, the Fonds Unique Interministériel (QuITO project), the Centre National de la Recherche Scientifique, the

Figure 1 :

 1 Figure1: Schematic of the experimental setup. The sample is illuminated using a Köhler configuration with a reduced numerical aperture (D2) to increase the spatial coherence of the light. A laser is used to excite the gold structures at their plasmonic resonance and heat them (either at 532 or 808 nm, depending of the sample). It is focused at the objective entrance pupil by L1 to perform a wide field illumination. A diaphragm D3 is used to control the beam diameter and to achieve a uniform illumination at the sample location. This diaphragm is optically conjugated with the sample by the objective/L1 association The sample is imaged onto a CCD camera by L2 through a modified Hartmann grating. This creates an interferogram on the CCD that is processed to retrieve the optical path difference generated by the thermalinduced refractive index variation.

Figure 2 :

 2 Figure 2: a) Scanning Electronic Microscope (SEM) image of the gold structure under study. This structure was illuminated at a wavelength λ 0 = 808 nm. b) Raw image of δ ℓ(x, y), I = 12.8 mW/µm 2 . c) Heat power density p(x, y) processed from image b using the NL algorithm. d) Temperature distribution δ T (x, y) processed from image c using Eq. (11). e) Measured radial profile of the OPD distribution in

Figure 3

 3 Figure 3: a) SEM image displaying a quasi-hexagonally ordered array of gold nanoparticles. b) Raw image of the measured optical path difference. Heating is performed by a circular beam of diameter D = 8.0 µm and uniform irradiance 105 µW/µm 2 at λ = 532 nm. c) Processed heat power density. d) Processed temperature increase. e) Radial profile of the temperature distribution. f) Optical image of the micropatterned array of gold nanodots. A SEM image in the inset displays their organization. This whole structure is illuminated by a uniform beam of irradiance 250 µW/µm 2 . g) Raw image of the measured optical path difference. h) Heat source density. i) Temperature increase. j) Maximum temperature of the microstructure as function of the laser irradiance.

Figure 4 :

 4 Figure4: a) Raw thermal-induced OPD image obtained on a quasi-hexagonally ordered array of gold nanoparticles (see Figure3a). b) Processed temperature distribution (smoothed using a 10-pixel Gaussian function). c) Cross section of image (a). Cross section of the raw image (b) without smoothing (dash line) and with smoothing (solid line).

7

  th Framework Program of the European Community (FP7) through the Marie Curie Action for career development (IEF -"Multi-PGNAs") and the Max Planck Society.


Guillaume Baffou
email: guillaume.baffou@fresnel.fr

Hervé Rigneault

Femtosecond-pulsed optical heating of gold nanoparticles

We investigate theoretically and numerically the thermodynamics of gold nanoparticles immersed in water and illuminated by a femtosecond-pulsed laser at their plasmonic resonance. The spatiotemporal evolution of the temperature profile inside and outside is computed using a numerical framework based on a Runge-Kutta algorithm of the fourth order. The aim is to provide a comprehensive description of the physics of heat release of plasmonic nanoparticles under pulsed illumination, along with a simple and powerful numerical algorithm. In particular, we investigate the amplitude of the initial instantaneous temperature increase, the physical differences between pulsed and cw illuminations, the time scales governing the heat release into the surroundings, the spatial extension of the temperature distribution in the surrounding medium, the influence of a finite thermal conductivity of the gold/water interface, the influence of the pulse repetition rate of the laser, the validity of the uniform temperature approximation in the metal nanoparticle, and the optimum nanoparticle size (around 40nm) to achieve maximum temperature increase.

I. INTRODUCTION

Gold nanoparticles (NPs) can act as efficient nanosources of heat under visible or infrared illumination at the plasmonic resonance due to enhanced light absorption. 1 The ability to locally heat at the nanoscale opens the path for new promising achievements in Nanotechnology and especially for nanoscale control of temperature distribution, 2 chemical reactions, 3 phase transition, 4 material growth, 5 photothermal cancer therapy [6][7][8] and drug release. 9,10 The use of femtosecond-pulsed illumination on gold NPs expands the range of applications compared to continuous (cw) illumination. First, it can lead to nonlinear optical processes like two-photon luminescence or second harmonic generation with applications mainly in bio-imaging. 11,12 Then, it can trigger a sudden temperature increase at the sub-nanosecond scale and subsequent effects such as acoustic waves used for opto-acoustic imaging 13,14 or bubble formation for nanosurgery. 15 A sharp and brief temperature increase of a NP generated by a femtosecond laser can also contribute to confine the temperature increase at the close vicinity of the NP to avoid extended heating of the whole medium when not desired. 16 Several experimental and numerical approaches aimed at studying the internal processes of heat generation under pulsed illumination and the subsequent effects observed in the surrounding medium, e.g. temperature and pressure variations, [17][18][19] acoustic wave generation, 18 vibration modes, [20][21][22] , cell apoptosis, 11 drug release 9,23 nanosurgery, 15,24 bubble formation, [25][26][27][28] NP shape modification 29 and melting, [29][30][31] nanosecondpulses for biomedical applications, 32,33 extreme thermodynamics conditions. [33][34][35][36][37] In this paper, we present and use a versatile numerical framework to investigate theoretically and numerically the evolution of the temperature distribution of a gold nanoparticle immersed in water when shined by a femtosecond-pulsed laser. Various degrees of complexity exist to describe theoretically such a problem. We chose a progressive approach consisting in going from simple to more sophisticated: We shall start with the more ba-sic description of the problem, a point-like source of heat to model the NP, and then refine the description of the system by taking into account successive refinements, namely a finite-size structure, an gold/water interface conductivity and a non-uniform NP inner temperature. In each case, the physics and the associated constitutive equations are detailed and the approximations discussed. The aim is to answer all the questions related to characteristic time, space and temperature increase in fs-pulsed optical heating of gold NP.

Details regarding the numerical algorithm we developed are given in appendix.

II. RESULTS AND DISCUSSION

A. Physical system

We consider a system with spherical symmetry consisting of a gold nano-sphere of radius R immersed in water (Fig. 1). This nanoparticle is uniformly illuminated by a laser light at its plasmonic resonance angular frequency ω = 2π c/λ 0 = k c/n w , c is the speed of light and n w the optical index of water. No mass transfer like fluid convection or bubble formation is considered. We shall focus on moderate temperature increase (typically a few tens of degrees) and not consider extreme thermodynamic conditions. For this reason, all the parameters describing the materials (water and gold) are assumed to remain constant within the temperature range investigated. Whenever a femtosecond pulse is mentioned, it has to be understood as a pulse, the duration of which is smaller than the characteristic time of electron-phonon scattering τ e-ph ∼ 1.7 ps. For instance, this can correspond to the use of a Ti:Sapphire laser, which usually provides a pulse duration around 100 fs.

In the following, any mentioned temperature T will have to be understood as a temperature increase above this initial ambient temperature.

Anywhere in the system, we can define the thermal energy density and thermal current density that read respectively u th (r, t) = ρ c T (r, t) and j th (r, t) = -κ∇T (r, t) where ρ is the mass density, c the specific heat capacity at constant pressure and κ the thermal conductivity of the system at the position r. From the energy conservation equation ∂ t u th (r, t) + ∇ • j th (r, t) = p(r, t), on obtain the heat diffusion equation:

ρ c ∂ t T (r, t) = κ∇ 2 T (r, t) + p(r, t) (1)
where p(r, t) is the heat power density (nonzero only inside the NP (r < R), where the light is absorbed), For the system under consideration in this work, this yields a set of two differential equations, one for each medium (gold and water), along with two boundary conditions at the gold/water interface:

                 Diffusion equations: ρ Au c Au ∂ t T (r, t) = κ Au ∇ 2 T (r, t) + p(r, t) for r < R ρ w c w ∂ t T (r, t) = κ w ∇ 2 T (r, t) for r > R Boundary conditions at r = R: κ w ∂ r T (R + , t) = κ Au ∂ r T (R -, t) T (R + , t) = T (R -, t) (2)
Au and w subscripts refer to gold and water. The first boundary condition ensures heat flux conservation at the NP interface.

It has been demonstrated experimentally that an interface resistivity at the gold/water interface exists and can play a significant role in the heat release. [38][39][40][START_REF] Alper | 42 Handbook of Physics[END_REF] The interface resistivity can reach appreciable values when the liquid does not wet the solid. The wetting depends on the nature of the interface, and in particular a possible molecular coating. Namely, hydrophobic coatings are associated to poor thermal conductivities. The direct consequence of a finite interface conductivity g (or resistivity 1/g) is a temperature drop/jump/discontinuity ∆T at the NP interface such as:

P(t) = 4πR 2 g ∆T (t) ( 3
)
where ∆T is defined such as:

∆T (t) ≡ T (R -, t) -T (R + , t) (4)
(in this work, the symbol ≡ symbolizes a definition).

The released heat power P(t) is also related to the temperature gradient on the NP surface through the energy conservation equation:

P(t) = -4πR 2 κ w ∂ r T (R, t). (5)
Equations ( 3) and ( 5) yield a modification of the second boundary condition of system (2) at the nanoparticle interface r = R:

-∂ r T (R + , t) = 1 l K ∆T (t) (6)
where l K = κ w /g is named the Kapitza length and 1/g is the associated Kapitza resistivity.

Let us define from now two dimensionless constants that we shall often use in the following,

β ≡ ρ w c w ρ Au c Au ≈ 1.680 (7) γ ≡ κ Au κ w ≈ 512 (8)
and one dimensionless parameter:

λ K ≡ κ w g R = l K R (9)
that is the Kaptiza length normalized by the NP radius R. We also introduce dimensionless space ρ and time τ variables defined such as:

ρ ≡ r/R (10) τ ≡ a w t/R 2 (11)
where a = κ/ρc is called the thermal diffusivity. R and R 2 /a w are indeed the natural space and time units associated to the system. Within this work, we will make an extensive use of dimensionless variables and constants, first because it yields simpler, more natural formulae and more general results, and then because it shows how the algorithm can be properly written, i.e. working with numbers close to one and not unnecessary powers of ten. However, for the sake of simplicity, even when using normalized variables, we shall keep the same function names -e.g. T (r, t) and T (ρ, τ ). But this mathematical digression should not cause any clarity issue. Using these dimensionless variables, parameters and constants, the set of equations ( 2) along with the additional boundary condition (6) read then:

                     Diffusion equations: ∂ τ T (ρ, τ ) = γ β ρ 2 ∂ ρ ρ 2 ∂ ρ T (ρ, τ ) + p(ρ, τ ) for ρ < 1 ∂ τ T (ρ, τ ) = 1 ρ 2 ∂ ρ ρ 2 ∂ ρ T (ρ, τ ) for ρ > 1 Boundary conditions at ρ = 1: ∂ ρ T (1 + , τ ) = γ ∂ ρ T (1 -, τ ) = - 1 λ K ∆T (τ ) (12)
where the Laplacian operator ∇ 2 has been reformulated using spherical coordinates.

From now on, the use of dimensionless formalism will not be systematic, but preferred when it simplifies the notations and make the results more general.

B. Numerical method

The set of equations (12) has no simple analytical solution. We chose to solve it numerically by developing a finite difference method (FDM) and a Runge-Kutta (RK) algorithm. [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF] Basically, it is based on a spatio-temporal discretization of the system of equations (12) according to:

ρ i ≡ i × δρ i ∈ [0, N ] τ j ≡ j × δτ j ∈ [0, M ] T i,j ≡ T (ρ i , τ j ) ∂ ρ T (ρ, τ ) → T i+1,j -T i,j δρ ∂ τ T (ρ, τ ) → T i,j+1 -T i,j δτ . (13)
This discretization procedure is associated with a RK algorithm of the fourth order (RK4) that ensures a higher accuracy -compared to regular Euler algorithms of the first order -in the estimation of T i,j at each spatiotemporal step. Further details regarding the numerical algorithm are given in the appendix.

C. cw illumination

Before studying what occurs under femtosecond (fs) pulsed illumination, it is worth describing first what happens under continuous (cw) illumination. We consider in this paragraph a uniform cw illumination of irradiance I and wavelength λ 0 .

Under cw illumination, the establishment of the steady state temperature profile will be preceded by a transient evolution. By dimensional analysis of the two diffusion equations of system (2), one find that two time scales come into play:

τ w d = ρ w c w κ w R 2 = R 2 a w , (14)
τ Au d = ρ Au c Au κ Au R 2 = R 2 a Au . ( 15
)
τ w d is the characteristic time associated to the evolution of the temperature profile in the surrounding water while τ Au d characterizes the temperature evolution inside the gold NP. Since a Au ≫ a w , the thermalization inside the NP occurs much faster. Consequently, one can consider that the global establishment of the temperature profile of the overall system is governed by the time scale τ w d .

We consider now the final steady state regime. The set of equations reads now:

         κ Au ∇ 2 T (r) = -p(r) for r < R κ w ∇ 2 T (r) = 0 for r > R κ w ∂ r T (R + ) = κ Au ∂ r T (R -) = -g ∆T (16)
where ∆T ≡ T (R -)-T (R + ). If one consider an average power density p 0 = P 0 /V , the solution has a simple form and reads:

T cw (r) = P 0 4π κ w r for r > R (17)
T cw (r) = P 0 4π κ w R 1 + 1 2γ 1 - r 2 R 2 + λ K for r < R (18)
where P 0 is the heat power dissipated in the NP. Note that the temperature profile outside the NP does not depend on the NP surface conductivity. [START_REF] Baffou | [END_REF] Since γ ≫ 1, one can usually consider -whatever the NP size -that the inner temperature of the NP is uniform and equals:

T cw NP = P 0 4π κ w R (1 + λ K ) = σ abs I 4π κ w R (1 + λ K ) (19)
where σ abs is the optical absorption cross section of the NP. 2,45 For spherical nanoparticle smaller than 2R ∼ 30 nm, a good approximation can be used: 46

σ abs = k Im(α) - k 4 6π |α| 2 (20)
where

α = α 0 1 -(2/3)ik 3 α 0 ( 21
)
and

α 0 = 4π R 3 ε Au -ε w ε Au + 2 ε w ( 22
)
with ε Au the gold permittivity and ε w the water permittivity.

From equation (19), we obtain thus a simple and very useful formula that gives the steady-state temperature of small spherical NP (R < 15 nm) under cw illumination:

T cw NP = k I R 2 κ w Im ε Au -ε w ε Au + 2 ε w , (23)
T cw NP ≈ 2.00 × k I R 2 κ w . ( 24
)
The latter formula applies for λ 0 = 520 nm. No interface resistivity is assumed in this formula (λ K = 0). If the NP is not necessarily smaller than 2R = 30 nm, the previous formalism becomes inappropriate and the more sophisticated and general Mie theory has to be used. [START_REF] Bohren | Absorption and scattering of light by small particles[END_REF] Within this model, the absorption cross section reads:

σ abs = 2π k 2 ∞ j=1 (2j + 1) (|a j | 2 + |b j | 2 )
where

a j = m ψ j (v) ψ ′ j (u) -ψ j (u) ψ ′ j (v) m ψ j (v) ξ ′ j (u) -ξ j (u) ψ ′ j (v) , b j = ψ j (v) ψ ′ j (u) -m ψ j (u) ψ ′ j (v) ψ j (v) ξ ′ j (u) -m ξ j (u) ψ ′ j (v)
and m 2 = ε Au /ε w , u = k R and v = m u. The primes indicates differentiation with respect to the argument in parenthesis. ψ j and ξ j are Ricatti-Bessel functions defined such as:

ψ j (x) = π x 2 J j+ 1 2 (x), ξ j (x) = π x 2 J j+ 1 2 (x) + i Y j+ 1 2 (x) .
J ν and Y ν are the Bessel functions of first and second order respectively. The derivatives can be expressed as follows:

ψ ′ j (x) = ψ j-1 (x) - j x ψ j (x), ξ ′ j (x) = ξ j-1 (x) - j x ξ j (x).
In the following, the Mie theory will be used whenever a calculation of the absorption cross section is required.

D. Pulsed illumination & initial temperature increase

We consider now a fs-pulsed illumination of average irradiance I , pulsation rate f and wavelength λ 0 .

The absorption of a fs-pulse by a metal nanoparticle can be described as a three-step process [START_REF] Grua | [END_REF]48 that involves different time scales:

Electronic absorption. During the interaction with the fs-pulse, part of the incident pulse energy is absorbed by the gas of free electrons of the NP, much lighter and reactive than the ions of the lattice. The electronic gas thermalizes very fast to a Fermi-Dirac distribution over a time scale τ e ∼100 fs. 48 This leads to a state of nonequilibrium within the NP: the electronic temperature T e of the electronic gas has increased while the temperature of the lattice (phonons) T p remains unchanged. The absorbed energy E 0 reads:

E 0 = σ abs I /f = P 0 /f. (25)
Electrons-phonons thermalization. Subsequently this hot electronic gas relaxes (cools down), through internal electron-phonon interaction characterized by a time scale τ e-ph to thermalize with the ions of the gold lattice. This time scale is not dependent on the size of the NP except for NP smaller than 5 nm due to confinement effects. 49 Above this size and for moderate pulse energy, the time scale is constant and equals τ e-ph ∼ 1.7 ps. [50][51][52] At this point, the NP is in internal equilibrium at a uniform temperature (T e = T p ), but is not in equilibrium with the surrounding medium that is still at the initial ambient temperature.

External heat diffusion. The energy diffusion to the surroundings usually occurs at higher characteristic time scale τ d , which leads to a cooling of the NP and a heating of the surrounding liquid. The time scale of this process depends on the size of the NP and ranges from 100 ps to a few ns. For small NP (< 20 nm), this third step can overlap in time with the electron-phonon thermalization 19 (as discussed hereafter).

If one consider that the electron-phonon thermalization (step 2) occurs much faster than the external heat diffusion (step 3), the NP temperature reaches an initial maximum temperature T 0 NP that is straightforward to estimate by energy consideration. It is related to the absorbed energy through the relation:

E 0 = V ρ Au c Au T 0 NP (26)
where V is the volume of the NP and V ρ Au c Au its heat capacity. Using equation (25), we find that the maximum initial NP temperature is:

T 0 NP = σ abs I V ρ Au c Au f . ( 27
)
This formula is not restricted to spherical nanoparticles.

As an example, for a gold nanorod, 50 nm long and 12 nm in diameter, at the plasmonic resonance (λ 0 =800 nm), considering a random orientation, f = 86 MHz and I =1.0 mW/µm 2 , we obtain T 0 NP ≈ 30 • C. Note that for a given laser power, the temperature increase does not depend on the pulse duration, but only on the pulse energy I /f . It is worth comparing the instantaneous temperature increase T 0 NP after a single fs-pulse and the steady-state temperature T cw NP achieved under cw illumination. From Eqs. ( 19) and ( 27), we obtain a dimensionless number η 0 that quantifies the gain obtained when using pulsed illumination:

η 0 ≡ T 0 NP T cw NP = 3 β a w f R 2 (1 + λ K ) . (28)
Figure 2 illustrates what η 0 , T 0 NP and T cw NP represent on an particular example. This useful formula is only valid for large values of η 0 . When R or f tends to be high, successive pulses may overlap (as explained later in section II H), which makes the assumption of a single pulse illumination wrong. Furthermore, for very small particles, a temperature damping effect occurs (by a typical factor 2) and T 0 NP is not the initial temperature, as explained hereafter. Consequently, this formula is a good approximation of the temperature gain achieved under fs-pulsed illumination compared to cw illumination as far as η 0 remains large and the NP radius not too small (R > 10 nm). The true values of η 0 -i.e. whatever the NP radius R and surface conductivity g and without approximations -is numerically computed and discussed in section II G.

E. Subsequent evolution of the temperature profile

We discuss now the subsequent evolution of the temperature field T (r, t) after a single-pulse illumination, i.e. after an initial temperature increase T 0 NP . In the ideal case consisting of a point-like NP (R → 0), the heat power density can be described by a Dirac distribution:

ρ w c w ∂ t T (r, t) = κ w ∇ 2 T (r, t) + E 0 δ(r)δ(t). (29)
This ideal problem has a simple analytical solution that reads: Then, the envelope T max (r) ≡ max t (T (r, t)) of the temperature profile over time can be easily obtained by using formula (30) and calculating the time t for which ∂ t T (r, t) = 0 for any position r. It yields a temperature envelope:

T (r, t) = E 0 c w ρ w 1 (4πa w t) 3/2 exp - r 2 4πa w t . (30)
T max (r) = 1 3 √ 3 E 0 c w ρ w 1 r 3 .
Interestingly, the temperature profile under pulsed illumination follows thus a 1/r 3 spatial decrease, which makes a real difference compared to continuous illumination characterized by a steady state profile in 1/r (see Eq. ( 17)). This is a first clue that pulsed illumination achieves a much higher temperature confinement around the NP.

We consider now the more general and realistic case of a spherical NP defined by a finite radius R = 0. In this case, no simple analytical solution exists but some approximations can be done to simplify the problem. First, the initial temperature profile T (r, 0) can be considered as uniform (equals T 0 NP ) inside the NP since the electronphonon thermalization usually occurs much faster than the external heat diffusion. Then, one can also suppose that the NP temperature remains uniform during the evolution of the system since κ Au ≫ κ w . Hence,

∀t, T (r, t) = T NP (t) for r < R. (31)
The validity of these two approximations will be investigated and discussed in more details in section II G. Under these hypotheses, the system of equations ( 2) can be simplified:

                
Diffusion equation:

ρ w c w ∂ t T (r, t) = κ w 1 r 2 ∂ r r 2 ∂ r T (r, t) for r > R Boundary conditions: V ρ Au c Au dT NP (t) dt = κ w 4πR 2 ∂ r T (R, t) = -g 4πR 2 ∆T (t) (32)
The first equation is the heat diffusion equation outside the NP. The two other equations come from considerations of energy conservation and will control the boundary condition at the NP interface (r = R).

Interestingly, two new characteristic times arise from the boundary equations. They read:

τ NP d = R 2 ρ Au c Au 3 κ w , (33)
τ s d = R ρ Au c Au 3 g . ( 34
)
These characteristic times are associated to the evolution of the nanoparticle average temperature -i.e. the heat energy stored in the NP. When the surface resistivity 1/g is high (resp. small),the evolution is governed by τ s d (resp. τ NP d ). These new time scales differ from τ Au d and τ w d that were respectively the characteristic times associated to the establishment of the internal temperature equilibrium inside the NP and temperature diffusion outside the NP. If the surface conductivity g is small enough, the evolution of the NP temperature is governed by τ s d . If the surface resistivity is negligible, the evolution is governed by τ NP d . The temperature evolution inside and outside the NP is thus governed by a subtile interplay between four time scales. Three of them are linked by the relation:

τ w d = γ β τ Au d = 3 β τ NP d , (35)
which yields:

τ w d ≈ τ NP d ≫ β τ Au d . (36)
Regarding the fourth time scale τ s d , it can be dominant or negligible depending on the values of R and g.

Using dimensionless time τ and space ρ variables, the system of equations ( 32) can be recast into this simpler form:

               Diffusion equation: ∂ τ T (ρ, τ ) = 1 ρ 2 ∂ ρ ρ 2 ∂ ρ T (ρ, τ ) for ρ > 1 Boundary conditions: dT NP (τ ) dτ = 3 β ∂ ρ T (1, τ ) = - 3 β λ K ∆T (τ ) (37)
where ∆T (τ ) = T NP (τ ) -T (1, τ ). Using the RK4 algorithm, system (37) was resolved numerically. The discretization parameters are (δρ, δτ, N, M ) = (15 × 10 -3 , 4 × 10 -6 , 400, 10 6 ). The interface resistivity 1/g was set to zero (non-zero values are discussed in the next section). The result is presented in Fig. 3 that displays the universal normalized evolution of the temperature of a NP. This profile applies for any particle size since it uses the normalized variables (ρ, τ ): For a given particle size R, the normalized coordinate ρ has to be multiplied by R and the normalized time τ by R 2 /a w to recover the actual coordinate r and time t. Hu and Hartland 19 have shown experimentally that the NP temperature can be conveniently fitted using a stretched exponential function:

F (τ ) = e -(τ /τ0) n . ( 38
)
We used this function to fit the evolution of the NP temperature as represented in the inset of Fig. 3 (dashed line). The optimized fit parameters are n = 0.39 and τ 0 = 0.041. This yields a useful formula giving the normalized NP inner temperature evolution for any particle radius R:

F R (t) = exp - a w t 0.041 R 2 0.39 . ( 39
)
Note that this useful formula assumes a zero interface resistivity 1/g. For a finite value of g, the values of the fit parameters τ 0 and n are different and have to be recalculated using the FDM.

Figure 4 aims at comparing the temperature profiles under pulsed and cw illuminations. It replots the series of temperature profiles of Fig. 3 along with the temperature envelope and the temperature profile of the steady state (cw illumination) given by Eq. (17). As derived in the previous section, when considering a point-like (R → 0) source of heat, the envelope of the temperature profile follows a 1/r 3 profile outside the NP. When considering now a finite-size NP, it appears that such a simple law does not exist. Instead, a stretched exponential function can also be conveniently used to fit the envelope of the spatial temperature profile in the surrounding water:

F (ρ) = exp - ρ -1 ρ 0 n . ( 40
)
The fit parameters are n = 0.45 and ρ 0 = 0.060. This results illustrate to what extent pulsed illumination achieves a much higher degree of temperature confinement compared to cw illumination.

Note that the evolution of the system is characterized by an energy conservation law. At any time t, the energy of the system is constant and reads:

E 0 = 4 3 πR 3 ρ Au c Au T NP (t) + ∞ R ρ w c w 4πr 2 T (r, t)dr, (41)
or using the normalized variables and constants defined above, the normalized energy reads:

ǫ 0 = T NP (τ ) T 0 NP + ∞ 1 3 β ρ 2 T (ρ, τ ) T 0 NP dρ = 1. (42)
This conservation law can be conveniently used in numerical simulations as a verification of the consistency of the result. For example, in the simulations shown in figure 3, it varied by less than 0.2%.

F. Finite conductivity of the gold-water interface

In this section, we shall go one step further into the refinement of the analytical description of the problem. We still consider the NP temperature as uniform, but we take into account a finite interface conductivity g. The set of equations describing the system is given by (37). Fig. 5 displays the evolution of the temperature profiles for a set of different normalized Kapitza length (i.e. gold/water interface conductivity). Usual values of the interface conductivity g range from 50 to ∞ MW/(m 2 K). 38,[START_REF] Alper | 42 Handbook of Physics[END_REF] As an example, gold nanorods coated with CTAB molecules are endowed with a typical surface conductivity of 130 MW/(m 2 K). 38 For NP radii ranging from 5 to 50 nm, typical normalized Kaptiza length λ K are thus ranging from approximatively 0 to 2.

For large Kapitza length (high interface resistivity), the heating of the surrounding fluid can be highly inefficient as observed in Fig. 5a. This is due to the fact that the high resistivity of the NP interface tends to slow down the heat release, which yield a weaker maximum temperature in water. The same amount of energy is released in the surroundings but more softly.

G. Beyond the approximation of instantaneous temperature increase

In the previous section, we used two approximations regarding the NP temperature to simplify the problem.

First, we have considered that the NP temperature increase T 0 NP subsequent to the pulse of light was instantaneous and uniform inside the NP. This was because the characteristic time of electron-phonon scattering (τ e-ph ∼ 1.7 ps) is usually smaller than the diffusion time into the surrounding (see Eq. ( 33)), but this may become non-valid for small NP.

Second, we assumed that the NP temperature remained uniform inside the NP during the overall evolution. That was because the thermal conductivity of gold is much larger than the one of water.

Yet, the validity of these two approximations may depend a priori on the size R and the interface conductivity g of the NP. The numerical algorithm used in this work allows one to consider the inner temperature of the NP as not necessarily uniform and to investigate the step where the temperature increases while the NP is being illuminated by the pulse. We can thus investigate the validity of the two approximations mentioned above by simulating the evolution of the inner temperature profile, both during the electron-phonon thermalization and during the subsequent heat diffusion in the surrounding medium. This is the purpose of this section.

During the electron-phonon thermalization, one can assume that the heat power density p(r, t) is uniform all over the NP.

p(r, t) = p 0 (t) = E 0 V τ e-ph exp(-t/τ e-ph ). ( 43
)
Indeed, we have seen in section II D that the electronic thermalization (step 1) occurs much faster than the energy transfer from the electrons to the lattice (step2). The heating of the lattice is thus performed by a uniformly hot electronic gas. The dimensionless form reads:

p 0 (τ ) = T 0 NP τ w p τ e-ph exp(-τ τ w p /τ e-ph ).
Note that under cw illumination, the heat generation density is on the contrary highly non-uniform within the NP. 53 Within this new approach characterized by an initial zero temperature profile and a heat power density in the NP, the energy conservation law reads now, at any . The ideal temperature profile that the NP would achieve if there were no heat release in the surrounding is represented by a dashed line. c) Evolution of the temperature T (r = 0, t) for the two above-mentioned cases along with the normalized heat power p0(t) (dark line) delivered by the excited electronic gas to the phonons. The normalized temperature does not reach 1 due to heat release during the NP heating.

normalized time τ , under cw illumination:

3β T cw NP τ = 1 0 3ρ 2 T (ρ, τ )dρ + ∞ 1 3ρ 2 β T (ρ, τ )dρ ( 44
)
and under pulsed illumination:

T 0 NP 1 -exp -τ τ w d τ e-ph = 1 0 3ρ 2 T (ρ, τ )dρ + ∞ 1 3ρ 2 β T (ρ, τ )dρ. ( 45
)
Figure 6 presents the results of the numerical simulations for two NP sizes R = 5 and R = 50 nm, and for 1/g = 0. It shows the temperature evolution during the heating phase from the initial zero uniform temperature up to the time when the NP temperature reaches its maximum value. We can see that when considering a possible heat leak into the surrounding during the pulse illumination, the NP temperature does not reach necessarily T 0 NP . For large NP, the inner temperature profile suffers from some distortion when the temperature reaches its highest value (Fig. 6a). However, the maximum temperature reaches practically the ideal maximum temperature T 0 NP . Indeed, one can see that the heat diffusion into the surroundings is almost absent. For small NP (Fig. 6b), the temperature diffusion into the surrounding medium during this initial step is much more visible. The direct consequence is that the NP temperature does not reach its maximum ideal value T 0 NP represented by dashed line in Fig. 6b. The evolution of the NP inner temperature (T (r = 0, t)) for the two above mentioned cases (R = 5 and R = 50 nm) is plotted in Fig. 6c. The temperature profiles are displayed along with the evolution of the heat power p 0 (t) provided by the free electron gas. For small particles, it is clear that the time scale τ e-ph becomes of the same magnitude as the characteristic time of the heat diffusion in the surrounding. This explains why the NP temperature cannot reach the maximum ideal temperature T 0 NP . Note that the temperature T 0 NP used to normalized the curves in Fig. 6 depends on R (cf. Eq. ( 27)) and is thus not the same for the two cases.

In any case, the temperature remains quasi-uniform inside the NP, which validate the usual assumption (31).

The maximum temperature achieved in the system is discussed in more details in Fig. 7. Temperatures are plotted as function of NP size R and interface conductivity g. Figs. 7a-b represent the maximum temperature inside and outside the NP:

T max NP = max t [ T (r = 0, t) ] , (46)
T max (R + ) = max t T (r = R + , t) (47)
normalized by the ideal temperature increase T 0 NP (Eq. ( 27) and Fig. 7g). When T max NP /T 0 NP ≪ 1, it means that the heat release outside the NP is too fast. The characteristic time τ w d (Eq. ( 14)) becomes close to the characteristic time of the inner thermalization and step 2 and 3 occur almost simultaneously. This happens when the particle is too small and when the surface conductivity is not weak. When T max w /T 0 NP ≪ 1, it means that the water temperature increase is not optimum. This happens when the NP is too small, or when g is too small, giving rise to a large temperature drop at the NP interface.

Figures 7c-d represent the maximum temperature inside and outside the NP for a given irradiance I = 0.1 mW/µm 2 . Interestingly, we evidence an optimal NP size. NPs around 2R = 40 nm turn out to be the most efficient nano-source of heat for a given laser irradiance. Below this size, the cooling of the NP is too fast and the temperature has no time to reach T 0 NP , as explained above. Above this size, the absorption cross section is no longer proportional to the volume of the nanoparticle, as seen in Fig. 7f and in Eq. (20), which tends to damp the temperature increase. This trend is also observed in Fig. 7d that plots the associated temperature increase at the vicinity of the NP under the same illumination conditions. Moreover, we can see that to achieve a high 17)). f) Absorption cross section of a spherical gold NP as function of radius obtained from Mie theory (Eq. ( 25)) (solid line). The deviation from the R 3 law is represented by a dashed line. g) Ideal temperature increase T 0 NP as function of NP radius. h) Temperature on the NP surface under cw illumination (calculated from Eq. ( 17)). The deviation from the R 2 law (see Eq. ( 24)) is represented by a dashed line.

temperature increase inside (resp. outside) the NP, a low (resp. high) surface conductivity g is preferred.

Figure 7e aims at comparing the efficiency of pulsed versus cw illumination. In Eq. ( 2), we defined the gain η 0 with the approximation of an instantaneous temperature increase and no interface resistivity. Now, we can define the exact gain η that takes into account the initial temperature damping and a finite interface conductivity:

η = T max (R + )/T cw (R + ). (48)
We observe that for small NP and high surface conductivity (upper left corner of Fig. 7e), a fs-pulsed illumination achieves a temperature rise two orders of magnitude higher than cw illumination of the same average irradiance. However, depending on the NP size and surface conductivity, a pulsed illumination does not necessarily achieve a more pronounced temperature increase in the surrounding medium. Namely, for big NP around 2R = 100 nm in diameter, it appears that the use of pulsed illumination becomes inefficient since the gain in temperature reach a maximum value of around 3. For low surface resistivity, the surrounding maximum temperature can even be higher when using cw illumination (lower right corner of Fig. 7e). Consequently, the use of a fs-pulsed laser is not necessarily efficient if one want to achieve the highest temperature possible for a given laser power. This result questions the utility of fs-pulse illumination on lithographic structure, since usually bigger than 100 nm.

H. Influence of the pulse repetition rate

When the relaxation time of the system is sufficiently long, the temperature of the NP has no time to return to zero between two successive pulses. This effect has not been considered so far and is the purpose of this section. The parameter that controls the appearance of such a regime is:

ξ ≡ f τ d = f (τ NP d + τ s d ) = f R 2 ρ Au c Au κ w (1 + λ K ). ( 49
)
Let G R,g (t) be the NP temperature evolution after a single-pulse illumination for a NP of radius R and interface conductivity g. Note that a fit formula of G R,∞ (t) for g → ∞, assuming a square initial temperature profile, is given by Eq. (39). Due to the linearity of all the equations governing the heat release and diffusion in the system, the NP temperature at any time t after a series of N pulses at the repetition rate f is:

T NP (t) = N -1 j=0 G R,g (t -j/f ). ( 50
)
Figure 8 plots the results of numerical simulations for two cases corresponding to two different NP radii. The interface resistivity 1/g is assumed to equal zero. For a radius of R = 50 nm (Fig. 8a), ξ is close to unity, which yields a temperature offset: the NP temperature has no time to return to zero between two successive pulses. For even larger NP, the temperature evolution would be even more smoothered and would tend to T cw NP . For weaker values of ξ, as shown in Fig. 8b, the successive temperature pulses do not overlap and a regime of time localization can be achieved.

I. Extension to non-spherical nanoparticles

The numerical technique we have developed could be extended to 2D systems modeling a NP with axial symmetry, like nanorods. It would require a 2D mesh and longer computation times. However, most of the results presented in this work can be easily extended to nonspherical NP without carrying out more sophisticated calculations: While the optical properties of NP particles depend very sensitively on the NP geometry -in particular the resonance frequency-, thermal processes of NP are only slightly dependent on the geometry. 2 For non-spherical NP, we can define an effective radius R eff such that the NP volume equals:

V = 4 3 πR 3 eff . (51)
Most of the reasonings we have done are also valid when replacing R by this effective radius R eff . 2 In particular, the equation ( 27) giving the ideal temperature increase T 0 NP remains exact whatever the shape of the NP using R eff . For gold nanorods, large aspect ratios would tend to make the cooling of the NP slightly faster by increasing the surface to volume ratio. However, all the orders of magnitude of time and space will remain identical as long as the aspect ratio remains moderate (less than ∼ 4).

III. SUMMARY AND CONCLUSION

To summarize, we have developed a versatile and yet simple numerical framework to investigate femtosecondpulsed optical heating of spherical gold nanoparticles (NPs). This approach enabled us to address different models, from simple to more sophisticated, and to discuss the validity of their approximations. Most of the reasonings and calculations are made using dimensionless variables, parameters and constants, which enabled us to compute universal behaviors. The constitutive equations governing the system are derived and explained. We show how such a physical system turns out to be governed by a subtile interplay between four characteristic time scales.

We investigated the influence of the size of the NP: For small NP (diameter < 30 nm), the temperature rise is not as high as expected due to fast heat release. For bigger NP, the temperature rise is damped because large NP are not efficient absorbers. This leads to a particle size compromise around 40 nm that optimizes the temperature increase (for a given laser irradiance).

The role of a possible molecular coating is also investigated in details based on considerations on interface thermal resistivity. For high interface resistivity, a temperature damping is observed outside the NP while an enhanced temperature increase is observed inside. A high interface resistivity tends in parallel to slow down the heat release in the surrounding medium making the heating of the surroundings less efficient.

A detailed comparison between cw and pulsed illumination was drawn and two main results came out: i) While a temperature profile in 1/r in the surrounding medium is observed under cw illumination, a much more confined temperature envelope in 1/r 3 characterizes a pulsed illumination. A refined model even further demonstrates a higher degree of confinement with a spatial profile following a stretched exponential. ii) Unexpectedly, pulsed illumination does not necessarily achieve much higher temperature increase in the surroundings compared to cw, especially for nanoparticles bigger than 100 nm (typically lithographic plasmonic structures). It can even be worse when the gold particle is endowed with a poor thermal surface conductivity (due to an hydrophobic molecular coating for example).

Finally, the influence of the repetition rate is discussed and two regimes are identified depending on the NP radius R and the pulsation rate f : One time-localization regime, where the temperature increase is confined spatially and temporally and one regime that tends to resemble to the regime observed under cw illumination.

Within this work, we restricted ourselves to gold nanoparticles with spherical geometry (radius R), but most of the results obtained herein are also valid for non-spherical particles when considered as particles of characteristic size R. The numerical techniques we developed could also be refined to address problems with 2D symmetries requiring a longer computation time. This versatile numerical technique could also take into account other materials than gold and water and various pulse durations, from femto-to nanosecond-scales.

We consider in this part a spherical NP endowed with a surface conductivity g and characterized by a uniform inner temperature T NP (τ ). This problem is described by the set of equations (37). Since the NP temperature is assumed to be uniform, only the surrounding medium has to be meshed (see Fig. 9). The spatio-temporal meshing of the system is such as:

ρ i ≡ 1 + i × δρ i ∈ [0, N ] τ j ≡ j × δτ j ∈ [0, M ] T i,j ≡ T (ρ i , τ j )
To simplify the notations and the explanations, we consider here a regular spatio-temperal mesh. However it would be wise to use for instance a refined mesh close the NP interface and a rougher mesh further.

The boundary condition (37) at the NP interface reads:

- T 2,j -T 1,j δρ = T 0,j -T 1,j λ K . (A1)
The temperature drop ∆T occurs between the coordinates ρ 0 and ρ 1 . The temperature gradient on the NP surface is calculated between the coordinates ρ 1 and ρ 2 (see Fig. 9). The initial temperature profile is set to:

T 0,0 = T 0 NP (A2) T 1,0 = T 0 NP 1 + δρ/λ K (A3) T i,0 = 0 for i > 1 (A4)
This initial temperature profile corresponds to a zero temperature outside the NP and fulfills the boundary condition (A1). The numerical algorithm consists in a spatial loop over the position i inside a temporal loop over the time j. At each time j + 1, here is the procedure.

First, the NP temperature T 0,j+1 is calculated from the second equation of system (37): T 0,j+1 = T 0,j + 3 β δτ T 2,j -T 1,j δρ . (A5)

Then, using Eq. (A1), one can compute the temperature on the NP surface (i = 1):

T 1,j+1 = l K T 2,j + δρ T 0,j+1 λ K + δρ . (A6)
Finally, the computation of the temperature profile in the surrounding medium (i > 2) is performed according to the RK4 procedure: 43

T i,j+1 = T i,j + k 1 + 2 k 2 + 2 k 3 + k 4 6 (A7)
where k 1 = δt K(T i,j ) k 2 = δt K(T i,j + k 1 /2) k 3 = δt K(T i,j + k 2 /2) k 4 = δt K(T i,j + k 3 ) and K(T i,j ) = 2 ρ T i+1,j -T i,j δρ + T i+1,j -2 T i,j + T i-1,j δρ 2 .

(A8) Note that the increments δρ and δτ cannot be chosen arbitrarily and independently. To ensure a proper convergence, δρ and δτ have to be chosen much smaller than unity and respecting as well the convergence criteria:

δρ 2 /δτ ≫ 1. (A9)
This ensures that the right hand member of equation (A7) -the variation of the NP temperature -is small compared to the NP temperature T i,j . It is worth estimating the total energy of the system and its conservation at each time step to check the consistency of the The normalized energy at time τ j expressed from equation (41) reads:

ǫ j = T 0,j T 0,0 + 3 β δρ N i=1 ρ 2 i T i,j T 0,0 (A10)
and should remain close to unity at any step j of the numerical procedure.

Model including a non-uniform NP temperature

We assume in this section that the inner temperature of the NP is not necessarily uniform. We shall also explain how both the heating and cooling of the NP can be investigated numerically.

In this context, the normalized coordinate reads now:

ρ i ≡ i × δρ i ∈ [0, N 1 +N 2 -1]
and in particular ρ N1 = N 1 × δρ = 1. N 1 is the number of mesh points in the NP and N 2 in the surrounding medium (see Fig. 10).

To simplify the notations and the explanations, we consider here a regular mesh. However it would make sense to use for instance a rough mesh inside the NP and a refined mesh in the surrounding medium, close to the NP interface.

The discretized boundary conditions of system (12) yields this time: T N1+2,j -T N1+1,j δρ = -T N1,j -T N1+1,j λ K = γ T N1,j -T N1-1,j δρ (A11)

At time τ 0 , the initial temperature profile is set to zero: ∀i, T i,0 = 0. At time τ j+1 , the calculation of the temperature profile inside the NP (i < N 1 ) is performed using the RK procedure (Eq. (A7)) and replacing the RK function K(T i,j ) with K Au (T i,j ) = β γ K(T i,j ) + p 0 (τ j ). This source term p 0 (τ j ) has to be added when one wants to take into account a heat generation in the NP. For a pulsed illumination, it reads:

p 0 (τ j ) = T 0 NP τ p τ e-ph
exp(-τ j τ p /τ e-ph ).

And for a cw illumination:

p 0 = 3β T cw NP .
The convergence criteria is now:

δρ 2 /δτ ≫ β γ ≈ 860. (A12)
This means that while investigating the dynamics of the NP temperature, the convergence criteria is more drastic. However, this is expected and not an issue since the associated time scale of the temperature diffusion inside the NP is much faster. Consequently, δτ can be chosen much smaller.

Then the temperatures at the nanoparticle interface T N1,j and T N1+1,j have to be calculated using the boundary conditions (A11). It yields: T N1,j+1 = γ (δρ + λ K ) T N1-1,j + δρT N1+2,j γ λ K + δρ + γ δρ , T N1+1,j+1 = T N1+2,j + γ (T N1-1,j -T N1,j+1 ).

Finally the temperature profile within the surrounding medium (i > N 1 + 1) is computed using the RK4 procedure (Eq. (A7)) and the regular RK function K(T i,j ) (Eq. (A8)).

When considering an initial zero temperature profile, the new energy conservation laws ( 44) and ( 45) reads for cw illumination:

FIG. 1 .

 1 FIG.1. Nature of the system investigated: A spherical gold nanoparticle immersed in water shined by a femtosecondpulsed laser.

FIG. 2 .

 2 FIG.2. Evolution of the temperature of a NP a radius R = 50 nm under cw illumination (dashed line) and after a singlepulse illumination (solid line) for a given irradiance I = 0.1 mW/µm 2 . This figure aims at defining visually the initial temperature increase T 0 NP under pulsed illumination and the steady state temperature T cw NP under cw illumination. For this particular case, T 0 NP = 31 • C while T cw NP = 9.2 • C.

FIG. 3 .

 3 FIG. 3. Result of a numerical simulation showing the universal evolution of the temperature profile for a NP of radius R at different normalized time τ = awt/R 2 . No interface resistivity is considered (1/g = 0). The inset shows the evolution of the NP inner temperature as function of time.

FIG. 4

 4 FIG. 4. a) Comparison between the steady state temperature profile under cw illumination and the envelope of the temperature evolution under pulsed illumination. The inset shows the fit of the NP temperature profile in the case of a pulsed illumination. b) Temperature profile around the NP under cw illumination. c) Envelope of the temperature profile over time around the NP during subsequent to a fs-pulse illumination.

FIG. 5 .

 5 FIG. 5. a, b, c) Temperature envelope of the temperature evolution for three different values of normalized Kapitza length λ K . Temperature profiles are also represented at different times (in gray). d) Evolution of the NP inner temperature corresponding to the three previous situations.

FIG. 6 .

 6 FIG.6. Temperature rise subsequent to a femtosecond-pulsed illumination for a NP 100 nm in diameter (a) and 10 nm in diameter (b) (interface resistivity 1/g = 0). The ideal temperature profile that the NP would achieve if there were no heat release in the surrounding is represented by a dashed line. c) Evolution of the temperature T (r = 0, t) for the two above-mentioned cases along with the normalized heat power p0(t) (dark line) delivered by the excited electronic gas to the phonons. The normalized temperature does not reach 1 due to heat release during the NP heating.

NP surface conductivity g (x10 6 W

 6

FIG. 7

 7 FIG. 7. a) Numerical calculations of the temperature maximum of the NP normalized by the ideal temperature increase T 0 NP for a set of different NP radii and surface conductivities. b)Numerical calculations of the temperature maximum achieved in the surrounding medium normalized by the ideal temperature increase T 0 NP . c) NP maximum temperature when illuminated by a pulsed laser, at a pulsation rate f = 86 MHz and an irradiance I = 0.1 mW/µm 2 . d) Maximum temperature of the surrounding medium under the same illumination conditions. e) Ratio between the maximum temperature on the NP surface achieved under pulsed illumination and the temperature on the NP surface under cw illumination (Eq. (17)). f) Absorption cross section of a spherical gold NP as function of radius obtained from Mie theory (Eq. (25)) (solid line). The deviation from the R 3 law is represented by a dashed line. g) Ideal temperature increase T 0 NP as function of NP radius. h) Temperature on the NP surface under cw illumination (calculated from Eq. (17)). The deviation from the R 2 law (see Eq. (24)) is represented by a dashed line.

FIG. 8 .

 8 FIG. 8. NP temperature during a series of 10 pulses, f = 86 MHz, I = 1 mW/µm 2 , for a NP radius R = 50 nm (a) and R = 10 nm (b). In figure a, the NP temperature corresponding to a cw illumination is T cw NP = 1.78 • C and is represented by a dashed line. In the second case, the NP temperature under cw illumination is T cw NP = 0.30 • C (not represented).

T NP (τ j ) ρ 2 ∆TFIG. 9 .

 29 FIG.9. Schematic of the temperature profile at time τj around a spherical gold nanoparticle as modeled in the FDM-RK4 algorithm. In this first model, the NP temperature T NP (τj) is assumed to be uniform and a temperature discontinuity ∆T occurs at the nanoparticle interface to take into account a finite interface conductivity.

FIG. 10 .

 10 FIG.10. Schematic of the temperature profile at time τj around a spherical gold nanoparticle as modeled in the FDM-RK4 algorithm. In this second model, the temperature inside the NP is also meshed.

ρ 2 ii=N1+1 3β ρ 2 i

 22 T i,j δρ And with pulsed illumination:T 0 NP 1 -exp -τ j τ w d τ e-ph = N1 i=0 3ρ 2 i T i,j δρ + N1+N2 T i,j δρ.

TABLE I .

 I Physical constants used in this work associated to gold and water. a .

	Name	Gold	Water	Unit
	Thermal conductivity κ	317	0.60	W/m/K
	Specific heat capacity b c	129	4187	J/kg/K
	Mass density ρ	19.32	1.00	×10 3 kg/m 3
	Thermal diffusivity a	127	0.143	×10 -6 m 2 /s
	a Values at approx. 25 • C taken from ref. [42]
	b at constant pressure

ACKNOWLEDGMENTS

We thank Christian Girard, Philippe Réfrégier, Damien Riedel and Jérôme Wenger for helpful discussions.

Appendix A: Numerical algorithm

In this section, we explain and detail how the physical system was modeled using a finite difference method (FDM) and in particular what the RK4 algorithm consists in.

We shall specifically use in this section the equations and the formalism based on dimensionless space and time variables (ρ and τ ) and dimensionless constants and parameters (β, γ and λ K ).

1. Model assuming a uniform NP temperature