Temperature and curvature-dependent thermal interface conductance between nanoscale-gold and water from molecular simulation

1 Department of Mechanical Engineering, The University of Texas at Dallas, Richardson, TX 75080, USA 2 Department of Chemistry and Biochemistry, The University of Texas at Dallas, Richardson, TX 75080, USA 3 Institut de Thermique, Mécanique, Matériaux EA 7548, University of Reims Champagne-Ardenne, Reims, France 4 Department of Bioengineering, The University of Texas at Dallas, Richardson, TX 75080, USA 5 The Center for Advanced Pain Studies, The University of Texas at Dallas, Richardson, TX 75080, USA 6 Department of Surgery, The University of Texas Southwestern Medical Center, Dallas, TX 75390, USA


Introduction
Interfacial heat transfer can be characterized by the thermal interface conductance (G), which quantifies the rate 11 at which heat flows through an interface between two adjacent media at different temperatures. Experimental studies 12 indicate that the thermal interface conductance of nanoparticles is affected by the chemical composition of the 13 nanoparticle, its ligand coating, and the surrounding medium. [20,[22][23][24][25][26] Molecular dynamics simulation studies of 14 solid nanospheres immersed in mono-atomic fluids and nanodroplets of n-decane in water suggest that interface 15 curvature (inversely proportional to particle size) and temperature can also affect nanoparticle thermal interface 16 conductance. [27][28][29] Recent molecular dynamics studies also report that AuNPs immersed in organic solvent exhibit 17 curvature dependent changes in thermal interface conductance. [30,31] Although previous experiments indicate that 18 characteristic times for AuNP cooling in water scale with particle surface area, [32] it remains unclear if interface 19 curvature or temperature affect the thermal conductance at nanoscale gold-water interfaces. These factors are of 20 practical importance in biomedical applications that utilize plasmonic heating of AuNPs because the AuNPs are 21 embedded in primarily aqueous media and may experience significant temperature variations. 22 In this article, we use atomistic molecular dynamics simulations to investigate the effects of temperature and 23 curvature on heat transfer through bare nanoscale-gold interfaces with water during continuous plasmonic heating. 24 Molecular dynamics simulations allow us to assess trends in the thermal interface conductance that can be difficult to 25 explore experimentally due to limits in measurement sensitivity and particle monodispersity. [33] We simulated four 26 nanoscale-gold surfaces with different curvature under two wetting conditions (moderate vs poor), and subjected each 27 gold structure to various heating intensities to account for the effects of temperature increases during plasmonic 28 heating. Altering the gold wettability allowed us to account for how differences in the gold-water interaction potential 29 could affect the heat transfer, while also approximating the effects of coating the gold surface with ligands of different 30 hydrophilicity. First, we confirmed that our results for the interface thermal conductance fall within the range of 31 values reported from experiments. By comparing our results along with those from experiments to the critical value 32 of the interface thermal conductance, we confirmed that the thermal interface conductance is an important element of 33 AuNP heat dissipation, particularly, for small AuNPs (≤ 10 nm in diameter). Next, we analyzed the heat flux 34 dependence and curvature dependence of the interface thermal conductance at the two wetting conditions. Our 35 analysis indicates that the thermal interface conductance increases with curvature under both wetting conditions, but 36 it increases non-linearly with heat flux only under moderate wetting of the gold surface. We found that the curvature 37 dependence of the interface conductance coincided with a curvature dependent increase in water adsorption at the 38 gold surfaces, consistent with trends observed previously for a simpler model of solid nanospheres immersed in 39 monoatomic fluid. [27] We also observed a temperature dependent shift in the water vibrational density of states that 40 improved the VDOS overlap with that of gold. This improved vibrational overlap may account for temperature 41 dependent changes in the interface thermal conductance under the moderate wetting condition. This work advances 42 our understanding of how temperature and curvature affect the thermal conductance across nanoscale gold-water 43 interfaces. 44 2 Results and Discussion 45 2.1 Analysis of the thermal interface conductance 46 In order to evaluate the impact of interface curvature on the heat transfer between nanoscale-gold and water, we 47 modeled in atomistic detail four nanoscale gold-water interfaces of different curvatures (Figure 1a): three spherical 48 AuNPs with diameters (d) of 3, 5, and 10 nm, as well as a planar gold-water interface to represent the zero curvature 49 limit. To estimate the thermal interface conductance we performed steady-state non-equilibrium molecular dynamics 50 simulations in which the gold structure was subjected to a constant heat power, Q, while an outer region of water 51 was maintained at 300 K as a heat sink (Figure 1b), approximating plasmonic heating of gold under continuous-wave 52 laser stimulation. Under these conditions, the steady-state temperature profile in the direction normal to the 53 interface between gold and water exhibits a temperature jump (∆T ) (Figure 1c), which is inversely proportional to 54 the thermal interface conductance, G, Figure 1. Simulation design. a) Images depicting the simulated gold surfaces. b) Simulation design for the spherical AuNP geometry. Heat is applied uniformly to the gold atoms while the dark blue region of water acts as a heat sink, maintained at T=300 K. c) Representative temperature profile for a AuNP under continuous heating, which demonstrates the temperature jump at the interface.

3/20
where q = Q/A is the heat flux across the interfacial area A, which is taken as the gold surface area: A = πd 2 for 56 spherical AuNPs with diameter d and A = 2L x L y for the planar gold surface. L x and L y are the lateral dimensions 57 of the gold surface within the simulation cell, and the factor of two accounts for the fact that two sides of the planar 58 surface are exposed to solvent. To consider the effects of temperature, we conducted simulations under various 59 heating powers (Q, Table S1, Supporting Information). Note that in some cases, the heating caused the AuNPs to 60 thermally expand (Table S2, Supporting Information), which was accounted for when computing the thermal 61 interface conductance.

62
G is strongly affected by solvent adsorption (i.e., wetting) at the interface, scaling with the strength of the 63 interaction between the solid and liquid atoms. [27,34] To account for how differences in wetting at the gold-water 64 interface may affect the curvature and temperature dependence of G, we modeled two wetting conditions: a moderate 65 wetting case using the gold-water interaction parameters developed by Merabia et al. to reproduce a gold-water 66 contact angle < 30°, [35] as well as a poor wetting case in which we reduced the gold wettability by decreasing the 67 strength of the gold-water interaction by roughly a factor of 4. There are three reasons to consider multiple wetting 68 conditions for the gold-water interface. First, there are various parameterizations of the gold-water interaction 69 available in the literature with different parameter values and energy functions, [36][37][38][39] the particular choice of which 70 can affect both the wetting and thermal interface conductance. [40] Second, the standard 12-6 Lennard-Jones 71 potential we used to model the gold-water interaction (Equation 5) does not explicitly include image charge effects 72 arising from polarization of the metal atoms at an interface, [41,42] which can also affect adsorption at the metal 73 surface. Polarization can also affect the thermal interface conductance. [43] Lastly, the hydrophilicity of solid-water 74 interfaces can be tuned by ligand functionalization, [44] with hydrophobic ligands yielding lower thermal interface 75 conductance than hydrophilic ones. [25,45] Altering the gold wettability approximates the effect of coating the gold 76 surface with ligands of varying hydrophilicity. 77 Our simulations fall within the previous measured range for G. The steady-state temperature profiles of gold and 78 water for each case are shown in Figure S1 and S2 (Supporting Information). The thermal interface conductance (G) 79 was determined according to Equation  (Figure 2), [20][21][22]26] while the poor wetting results fall within the range of values (G ∼ 5 − 50 83 MW m −2 K −1 ) reported from experiments of spherical AuNPs in organic solvents. [20,22,26] G values for CTAB 84 coated gold nanorods in water have also been reported, with estimates ranging from ∼ 50 − 450 MW m −2 K −1 . [24,48] 85 G of the planar surface in the moderate wetting case was ∼110 MW m −2 K −1 at all heating powers, which is close to 86 the previously reported values of 105 ± 15 MW m −2 K −1 for 100 nm and 110 ± 10 MW m −2 K −1 for 18.5 ± 3.0 nm 87 citrate-stabilized AuNPs in water. [21,26] The planar surface had a G ∼ 25 − 26 MW m −2 K −1 in the poor wetting 88 case. For spherical AuNPs the relative importance of G in regulating their heat dissipation can be gauged by the 91 Kapitza number, [1] where l K is the Kapitza length, [49,50] d is the particle diameter, and κ W is the thermal conductivity of the 93 surrounding water. If λ K ≪ 1, the heat transport is limited by heat diffusion through water and interfacial heat 94 transfer has little effect on the process. In contrast, if λ K ≫ 1 the heat transport is limited by the interface. By 95 setting λ K = 1, we can define a critical G value, Analogously, when G ≫ G c AuNP heat dissipation is limited by heat diffusion through water and when G ≪ G c it is 97 limited by the interface. Heat transfer through the interface only has a negligible impact on the kinetics of AuNP 98 Figure 2. Thermal interface conductance. The simulated nanoscale gold-water thermal interface conductance values (this study) are plotted versus particle diameter, along with experimental data for spherical nanoparticles in water. Note that the different values of G shown for the simulated AuNPs at a given size correspond to simulations at different heating power. Experimental data were collected from Wilson et al., [20] Plech et al., [21] Ge et al., [22] and Stoll et al. [26] Error bars were included for both G and d of the experimental data points if specified in the given references. Note that we used the lower limit value of 20 MW m −2 K −1 reported by Wilson et al. in the plot for the Au-thiol-toluene interface. [20] 5/20 heat dissipation in the former case (i.e., when G ≫ G c ). To gauge the relative importance of heat transfer through 99 the gold-water interface to AuNP heat dissipation, we compared the G values from our simulated cases as well as 100 those reported from experiments of AuNPs in water with G c evaluated across a range of AuNP sizes ( Figure 3).

101
From this analysis, we confirm that for AuNPs within the size range we have simulated (d ≤ 10 nm) interfacial heat Figure 3. Critical thermal interface conductance. The simulated nanoscale gold-water thermal interface conductance values (this study) are plotted versus particle diameter and compared to the critical value of the thermal interface conductance G c (solid line), along with experimental data for spherical AuNPs in water. Experimental data were collected from Plech et al., [21] Ge et al., [22] and Stoll et al. [26] The boundaries between the thermal diffusion limited regime (dark gray region) and the interface limited regime (light gray region) were taken as G = 10G c ≫ G c (dashed line) and transport is a particularly important factor in their heat dissipation, with G ≲ G c under both wetting conditions. For 103 the simulated poor wetting case, the G values of the 3 and 5 nm AuNPs even approach the boundary for interface 104 limited heat dissipation. Considering the experimental estimates of G for AuNPs in water, [21,22,26] it appears that G 105 only approaches the boundary for thermal diffusion limited heat dissipation around d ∼ 100 nm, suggesting that the 106 interface slows AuNP heat dissipation for diameters d < 100 nm. However, coating AuNPs with a sufficiently 107 hydrophilic ligand such as PEG may lead to much higher G, [25] which may make heat dissipation out of particles 108 with d < 100 nm thermal diffusion limited. Conversely, G for CTAB stabilized gold nanorods in water may be as low 109 as ∼ 50 MW m −2 K −1 in some cases, [48] so more hydrophobic ligands will promote interface limited heat dissipation 110 out of AuNPs with d ≤ 5 nm.  : Supporting Information), an increase in the thermal interface conductance with heat flux suggests an underlying 116 temperature dependence. However, in the poor wetting case, G did not exhibit a clear dependence on the heat flux 117 (Figure 4b), suggesting that wetting also plays a role in how G is affected by temperature.

118
Temperature dependent thermal conductance has been observed in other systems. A linear dependence on 119 temperature has been observed for thermal conductance at solid-solid interfaces at high temperature, which was 120 attributed to an increase in inelastic scattering of phonons at the interface. [51,52] Various simulation studies have 121 reported that the thermal interface conductance of nanoscale solid-fluid interfaces tends to increase with 122 temperature, [53][54][55] typically following a power law expression. However, others indicate a decrease in the thermal 123 interface conductance with temperature. [28,35,56,57] Chen et al. reported a dependence of the thermal interface 124 conductance on heat intensity for simulations of a 3 nm spherical AuNP immersed in a flexible water model that was 125 similar to the one we observed for the AuNPs, [40] but they also found that the thermal interface conductance was 126 higher overall and did not change with heating power with a rigid water model. [40]  132 where d is particle diameter and 2/d the interface curvature, δ is the change in interface conductance with curvature, 133 and G 0 is the interface conductance at the zero-curvature limit (i.e., lim 2/d→0 ) which should match the conductance 134 of the corresponding planar surface. Similarly, we found that at low and high heat flux the thermal interface

149
Curvature dependence of the interface thermal conductance has previously been reported from simulations of other 150 systems including simpler solid nanospheres in mono-atomic Lennard-Jones fluids, [27,28] and nanodroplets of 151 n-decane in water, [29] as well as solid and liquid ZnO nanoparticles in tetradecane. [58] In those cases, as with our 152 simulations of the gold-water interface, δ > 0 and the interface conductance increased with curvature (decreases with 153 particle size) (Table S3, Supporting Information). In contrast, G from simulations of bare spherical AuNPs in hexane 154 exhibited the opposite trend with δ < 0 (G decreases with curvature). [30,31] In this case, the curvature dependence 155 was attributed to size dependent changes in the structure, coordination, and vibrational modes of the AuNP surface 156 atoms. In the case of n-decane nanodroplets, the increase in G with curvature (i.e., δ > 0) was attributed to a 2.2 Analysis of water adsorption at the nanoscale-gold surfaces 161 The structuring and adsorption of solvent at a solid surface is an important factor in the interfacial heat transfer 162 and can have a large effect on thermal interface conductance. [34,45,59] For solid spherical nanoparticles, Tascini et al. 163 found that smaller nanoparticles with larger curvature had higher surface adsorption density [27] , which coincided 164 with the increases in G that they observed with the particle curvature. To determine if curvature dependent changes 165 in solvent adsorption could explain the curvature dependence of the AuNP-water G, we examined the structuring and 166 adsorption of water at the surfaces of nanoscale gold (Figure 5a). For each wetting condition, we first examined the 167 equilibrium density profiles of gold and water at 300 K (i.e., no heating) (Figure 5b,c). In the moderate wetting case, 168 the peak in the mass density corresponding to the adsorption layer was 1.7 g mL −1 for the AuNPs and 2.4 g mL −1 at 169 the planar interface (Figure 5b), whereas in the poor wetting case it was 1. context of proteins, [60,61] as a metric to help quantify the role of hydration in protein structure, folding and binding. 177 We computed the solvent accessible surface area, A SASA (see Section 4.3.5), of the initial configuration of each gold 178 structure and normalized it relative to its surface area (A = πd 2 for spherical AuNPs and A = 2L x L y for the planar 179 slab) to quantify the relative accessibility: A SASA /A. The relative accessibility A SASA /A increased linearly with the 180 interface curvature (Figure 6b). Since A SASA effectively quantifies the maximum amount of water that can adhere at 181 the gold surfaces independent of wetting strength and temperature/heating, the trend in A SASA /A suggests a 182 geometric effect in which there is more space for water to adsorb at the curved gold surfaces. This effect explains why 183 the more highly curved gold surfaces sustain higher ϱ in both the moderate and poor wetting conditions. 184 We also analyzed the steady-state non-equilibrium molecular dynamics simulations and evaluated the steady-state 185 density profiles of gold and water under each simulated heating power ( Figure S4 and S5, Supporting Information). 186 In most cases, there were only minor changes in the height and width of the peak corresponding to the water 187 adsorption layer with heating. However, in the poor wetting case the distinct peak for the adsorption layer near the 5 188 and 10 nm AuNPs disappeared at heating powers of Q = 4000 nW and Q = 8000 nW, respectively ( Figure S5  In the absence of electron-phonon coupling, interfacial heat transport is mediated by the transmission and 199 scattering of phonons (vibrational energy carriers) across the interface. [18,62] In particular, models of thermal 200 interface conductance based on diffuse mismatch emphasize coupling between the two materials' vibrational modes, 201 or phonon density of states, as a critical factor in heat transport through the interface. [62,63] In molecular dynamics 202 simulations, the vibrational coupling between nanomaterials and solvent can be qualitatively characterized by 203 computing their respective vibrational density of states (VDOS) and comparing the overlap between the two 204 distributions. [27,31,35,40] We therefore computed the VDOS (Equation 10) of each gold structure and the surrounding 205 water at 300 K (Figure 7a). The water distribution was broad, exhibiting a low intensity peak in the low frequency 206 domain and a broad peak at higher frequency, consistent with previous reports. [30,40] Only the low frequency water 207 peak has significant overlap with the VDOS of gold. The VDOS of the gold was concentrated in the low frequency 208 regime with two distinct peaks similar to the tranverse (TA) and longitudinal (LA) acoustic modes in bulk gold, [64] 209 Figure 5. Water adsorption at the surface of nanoscale-gold a) Representative snapshot of a cross-section view of the AuNP-water interface (without heating, moderate wetting) showing the first layer of water adsorbed onto the gold surface. b,c) Mass density profiles of gold and water in the absence of heating (T =300 K) for the moderate (b) and poor wetting (c) cases; the layer of water adsorbed onto the gold surface is shaded in purple. consistent with previous reports (Figure 7a). [40,65] Increasing the interface curvature leads to lower peak intensity 210 and broader outer tails in the gold VDOS, especially for the LA mode (Figure 7a, inset). Since the temperature of 211 gold and nearby water are elevated (above 300 K) during heating ( Figure S5, Supporting Information), we also 212 computed the vibrational density of states at 370 K and compared the results to those at 300 K. The 3 nm AuNP 213 exhibited the most obvious shift in VDOS (Figure 7b), but overall the temperature dependent changes in the gold 214 VDOS were minor (Figure 7b and Figure S6, Supporting Information). In contrast, the water density underwent a 215 noticeable redistribution to the lower frequency regime at 370 K (Figure 7c). The percentage of the water VDOS in 216 the range 0 to 6 THz increased from 11% at 300 K to 14% at 370 K, which suggests there may be better vibrational 217 coupling between gold and water at elevated water temperatures. Therefore, temperature dependent shifts in the 218 water VDOS may contribute to the increases in G with heat flux observed in the moderate wetting case.

219
To further characterize the level of overlap between the gold and water VDOS we computed the Bhattacharyya 220 coefficient (BC, Equation 11), [66] a measure of overlap between probability distributions which ranges from 0 (no 221 overlap) to 1 (perfect overlap), between the VDOS of each gold structure and water (Figure 7d). At 300 K, the BC 222 values ranged from 0.28-0.29. At 370 K, the BC values increased to 0.30-0.32, indicating improved overlap between 223 the gold and water VDOS. To further confirm this effect, we computed the cross-temperature BC values (Figure 7e). 224 Similarly, the BC only increased when the water temperature was 370 K, confirming that shifts in the water VDOS 225 are the primary factor improving the overlap between the gold and water VDOS. 226

227
In this work, we used atomistic molecular dynamics simulations to investigate the heat transfer through nanoscale 228 gold-water interfaces. We simulated four nanoscale-gold surfaces with differing curvature, subjected to various 229 heating intensities, and under two wetting conditions. An analysis of the critical value of the thermal interface 230 conductance confirmed that the thermal interface conductance is an important factor in the aqueous heat transfer of 231 particles in the size range we investigated. We found that for a given AuNP size, the thermal interface conductance 232 increased non-linearly with the applied heat flux under moderate wetting of the gold surface, and did not change 233 under poor wetting conditions. Analysis of the gold and water vibrational density of states showed a temperature 234 dependent shift in the water distribution that improved the overlap with the gold density of states. This shift in the 235 water VDOS may contribute to the increases in thermal interface conductance with applied heat flux under the 236 moderate wetting condition. Our results also showed that the thermal interface conductance between nanoscale-gold 237 and water increased linearly with interface curvature at both low and high heat flux under both wetting conditions. 238 Our analysis indicated that the curvature dependence of the interface conductance coincided with a curvature 239 dependent increase in water adsorption at the gold surfaces, consistent with that observed previously with a simpler 240 model of solid nanospheres immersed in monoatomic fluid. [27] 241 This work helps elucidate the temperature and curvature effects on the thermal conductance across nanoscale 242 gold-water interfaces. Our work here focused on bare spherical particles immersed in water. Although we considered 243 two conditions for the gold wettability which approximates the effect of coating the gold with ligands of different 244 hydrophilicity, there are additional effects that ligand functionalization can have on interfacial heat transfer that were 245 not directly accounted for in our study. [30,46,47] We also did not directly investigate how the observed temperature 246 and curvature dependent changes in the thermal interface conductance impact the spatiotemporal temperature 247 dynamics during pulsed laser stimulation of AuNPs, or how they might affect any particular applications of 248 photothermally heated AuNPs. Future work will be required to address these issues, as well as to experimentally 249 validate our findings. 250 Figure 7. Analysis of the vibrational density of states (VDOS). a) VDOS of gold and water at 300 K. The grey region denotes the frequency range 0 to 6 THz which is shown in more detail in the inset. b) VDOS of the 3 nm AuNP at 300 K and 370 K c) Comparison of the water VDOS at 300 K and 370 K. The shaded areas highlight the density in the range 0 to 6 THz with the labeled percentage by area. Note that in each case when plotting the VDOS we re-scaled the density by the maximum value of the 10 nm AuNP so that y-axis ranges from 0 to ∼ 1. d) Comparison of the BC values for each gold surface computed from the VDOS at 300 K and 370 K (gold and water are at the same temperature in each case). e) Cross-temperature BC values for the two combinations of gold and water at different temperature (300 K or 370 K).

252
We constructed four atomistic models of nanoscale-gold immersed in water: three spherical AuNPs with diameters 253 of 3, 5, and 10 nm, as well as a slab of gold with a planar gold-water interface. The nanoscale-gold structures were 254 constructed and solvated using VMD software. [67] The planar gold slab was constructed as a rectangular block of  All molecular dynamics simulations were carried out with the LAMMPS software (http://lammps.sandia.gov) 265 using periodic boundary conditions and a timestep of 1 fs. [68] Note that under periodic boundary conditions, the 266 planar gold slab extends semi-infinitely in the xy plane. 267 We modeled nanoscale-gold using the embedded atom model, [69] which is a many-body potential that more 268 accurately represents metallic bonding than standard pairwise potentials. Water was modeled using the flexible 269 SPCE-F model, [70] chosen primarily for its ability to reasonably estimate the experimental surface tension of water 270 at temperatures between 300-500 K. With the SPCE-F model, the van der Waals interactions are described by a 271 standard 12-6 Lennard-Jones potential, defined as for atoms i and j separated by distance r, with a cutoff distance of r c . ϵ ij is the energy coefficient (kJ mol −1 ) for bond distance b with force constant K b = 2322.20 kJ mol −1Å −2 and equilibrium bond distance b 0 = 1.0Å. The 277 intramolecular H-O-H angle is also maintained with a harmonic potential for angle θ with force constant K a = 191.7549 kJ mol −1 rad −2 and equilibrium angle θ 0 = 109.4°. The interaction 279 between gold and water was also modeled using a standard 12-6 Lennard-Jones potential. Only the oxygen-gold 280 interactions were included in the description of the gold-water interaction. In the moderate wetting case, the 281 parameters for the water oxygen and gold interaction were ϵ O−Au = 2.47 kJ mol −1 and σ O−Au = 3.6Å as developed 282 by Merabia et al. [35] to reproduce a moderate wetting condition with a gold-water contact angle < 30°. In the poor 283 wetting case, the parameters were adjusted to ϵ O−Au = 0.59 kJ mol −1 and σ O−Au = 3.383Å, yielding a weaker 284 gold-water interaction. The model parameters are summarized in Table S4 (Supporting Information). A global cutoff 285 of r c = 10Å was used for the Lennard-Jones interactions. The particle-particle particle-mesh solver was used for 286 electrostatics with a desired relative error in the forces of 0.0001 with a switching cutoff of 8Å for smoothly 287 transitioning to long range electrostatics. [71] The gold structures were prevented from drifting during the simulations 288 by tethering their center of mass to the origin using a harmonic spring with force constant 960 kJ mol −2Å −1 .

289
Initially, each nanoscale gold-water system was energy minimized and then equilibrated under NVT conditions for 1 290 ns at 300 K, followed by an additional equilibration phase under NPT conditions for 1 ns at 300 K and 1 atm. The 291 non-equilibrium steady-state simulations were carried out under NPH conditions at 1 atm. The gold structure was 292 subjected to continuous heating at a given heat power Q, distributed uniformly across the gold atoms. In order to 293 create a non-equilibrium heat flux an outer region of water was maintained at 300 K with a Langevin thermostat, [72] 294 serving as a heat sink. In the planar gold-water system, the heat sink was defined as the region encompassing z > 200 295 A and z < −200Å within the unit cell; note that these two regions are joined into a single region by the periodic 296 boundary conditions. In the AuNP systems, the heat sink region was defined using a radial distance criterion, namely 297 the region encompassing r > 65Å, r > 85Å, and r > 110Å from the origin, for the 3, 5, and 10 nm AuNPs, 298 respectively. The non-equilibrium steady-state simulations were run for 2, 3, 4 and 6 ns, for the 3, 5, and 10 nm 299 AuNPs and planar slab, respectively. For the poor wetting case, the simulation times were all 4 ns. Simulation analysis was carried out using a combination of LAMMPS functionality, Python code 302 (https://www.python.org), and VMD scripts (https://www.ks.uiuc.edu/Research/vmd). [67] Our analysis in 303 Python incorporated various libraries, including MDanalysis, [73,74] SciPy, [75] NumPy, [76] Matplotlib, [77] 304 seaborn, [78] Numba, [79] pandas, [80] and Jupyter. [81] Simulation snapshots were rendered from VMD using the 305 Tachyon ray tracer. [82] 306 In the analysis of the non-equilibrium steady state simulations, we discarded the first 1 ns of the AuNP system 307 trajectories and the first 3 ns of the planar slab system. In the poor wetting case, we adjusted the time to 2 ns for 308 the 5 and 10 nm AuNP systems. We visually inspected the time course of the gold temperature ( Figure S8    The average values and standard errors of the thermal interface conductance G values were estimated by block 320 averaging of the temperature profile. [83] For each block, ∆T was estimated as the temperature difference between the 321 outermost gold spatial bin and the innermost water spatial bin. To account for the possibility of thermal expansion of 322 the gold, the gold size (diameter for AuNPs and thickness for planar slab) was estimated from the outermost gold 323 spatial bin position (Table S2, Supporting Information). The corresponding G for the block was estimated according 324 to Equation 1. Block sizes of 200 ps, 250 ps, 500 ps, and 750 ps were used for the block averages for the 3, 5, and 10 325 nm AuNPs and planar slab, respectively. In the poor wetting case, the block sizes were adjusted to 200 ps, 250 ps, 326 400 ps, and 300 ps for the 3, 5, and 10 nm AuNPs and planar slab, respectively. The gold temperature was examined 327 using blocked standard error analysis with varying block sizes ( Figure S10 and Figure S11, Supporting Information) 328 to confirm that these block sizes resulted in uncorrelated data for block averaging. The reported error in G values 329 was taken as 1.96 times the standard error estimate (from block averaging) and thus represents an estimate of the 330 95% confidence interval. Instantaneous (spatial) mass density profiles for gold and water were computed separately with a spatial binning of 333 width 0.25Å (radial shells for AuNPs and rectangular blocks along the z-direction for the planar surface). At each 334 15/20 time step, the coordinates of atoms were shifted so that the center of mass of the gold structure was at the origin. In 335 the case of the planar surface, the profile was averaged over the positive and negative z-directions. The profiles for 336 systems in the absence of heating were computed by averaging over the final 500 ps of the equilibration simulations, 337 which were under NPT conditions at 300 K and 1 atm. For systems under heating, the steady-state profile was 338 computed by time averaging the instantaneous profiles over the final 1 ns of each non-equilibrium steady-state 339 simulation. The time averaged profiles were then further smoothed by applying a Gaussian filter with σ = 2 to 340 generate the final equilibrium or steady-state profile. The surface adsorption density (ϱ) was computed according to the following equations adapted from Tascini et 343 al.: [27] for spherical AuNPs, and the planar gold surface, where N w is the number of water molecules adsorbed on the surface of gold with surface area A, N A is Avogadro's 346 number, M w the molar mass of water, ρ w (r) is the radial mass density of water at radius r, and ρ w (z) is the mass 347 density of water along the direction perpendicular to the planar surface, z. The upper limits of the integrals, r max 348 and z max , were taken as the profile position corresponding to the minimum following the peak of the adsorption band 349 (for example, see Figure 5b,c). The integrals were estimated using numerical quadrature. Solvent-accessible surface area estimates were computed using the measure sasa function in VMD with a probe 352 radius of 4.04Å, [67] which corresponds to the minimum energy in the Lennard-Jones potential for the gold-water 353 interaction in the moderate wetting case.
where ω is the vibrational frequency and ⟨⃗ v(t) · ⃗ v(0)⟩ is the velocity autocorrelation function computed and averaged 358 over all atoms of given type (either gold or water) and over multiple time origins (denoted by angular brackets). The 359 velocity autocorrelation function was computed by a separate set of simulations executed with a time step of 0.25 fs 360 and NVE conditions, initialized from the final configurations of each system after equilibration under NPT conditions 361 at the specified temperature (300 or 370 K) and 1 atm. The simulations were executed for a total of 40000 time steps, 362 and the velocity autocorrelation function was computed for two time origins taken over two consecutive 20000 time 363 step intervals. The Bhattacharyya coefficient (BC) quantifying the level of overlap between the gold and water vibrational density 366 of states was computed using numerical estimation of the following integral: [66] 367 BC = VDOS Au (ω) × VDOS W ater (ω)dω , where ω is the vibrational frequency and VDOS Au (ω) and VDOS W ater (ω) are the vibrational density of states 368 (VDOS) of gold and water, respectively. The integral was taken on the entire range of frequencies over which the 369 VDOS was computed (0-32 THz). 370

16/20
Supporting Information 371 Supporting Information is available online or from the author.