A review of heat transport in solvated gold nanoparticles: Molecular dynamics modeling and experimental perspectives

Nanoscale View Article Online

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. REVIEW Cite this: Nanoscale, 2025, 17, 20803

View Journal | View Issue

A review of heat transport in solvated gold nanoparticles: Molecular dynamics modeling and experimental perspectives Md Adnan Mahathir Munshi,a Emdadul Haque Chowdhury,a Luis E. Paniagua-Guerra,a Jaymes Dionne,b Ashutosh Girib and Bladimir Ramos-Alvarado *a Turning gold nanoparticles (AuNPs) into nanoscale heat sources via light irradiation has prompted significant research interest, particularly for biomedical applications, over the past few decades. The AuNP’s tunable photothermal effect, notable biocompatibility, and ability to serve as vehicles for temperaturesensitive chemical linkers enable thermo-therapeutics, such as localized drug/gene delivery and thermal ablation of cancerous tissue. Thermal transport in aqueous AuNP solutions stands as the fundamental challenge to developing targeted thermal therapies; thus, this review article surveys recent advancements in our understanding of heat transfer and surface chemistry in AuNPs, with a particular focus on thermal boundary conductance across gold- and functionalized-gold–water interfaces. This review article highlights computational advances based on molecular dynamics simulations that offer valuable insights into nanoscopic interfacial heat transfer in solvated interfaces, particularly for chemically functionalized AuNPs. Additionally, it outlines current experimental techniques for measuring interfacial thermal trans- Received 12th June 2025, Accepted 21st August 2025

port, their limitations, and potential pathways to improve sensitivity. This review further examines computational methodologies to guide the accurate modeling of solvated gold interfaces. Finally, it concludes

DOI: 10.1039/d5nr02521d with a discussion of future research directions aimed at deepening our understanding of interfacial heat rsc.li/nanoscale transfer in solvated AuNPs, crucial to optimize thermoplasmonic applications.

1. Introduction

The utilization of gold nanoparticles (AuNPs) in the biomedical field has gained significant momentum in recent years due to their unique characteristics, such as tunable optical properties and remarkable biocompatibility.1,2 Upon light irradiation, free electrons in noble metal nanoparticles undergo oscillations,3,4 creating localized surface plasmon resonance (LSPRs), which enables efficient light absorption, scattering, and focusing.5,6 Metallic NPs with LSPR, known as plasmonic NPs, act as potent nanoscale heat sources. By tuning the absorption rate relative to the cooling rate, their thermal response can be optimized to meet specific requirements, giving rise to thermoplasmonics.7 As a result, their applications include photothermal reaction acceleration,8 chemical catalysis,9 solar energy harvesting,10 additive manufacturing,11 and thermal sensing at solid–liquid interfaces.12 a

Department of Mechanical Engineering, The Pennsylvania State University, University Park, PA 16802, USA. E-mail: bzr52@psu.edu b Department of Mechanical Engineering, University of Rhode Island, Kingston, RI 02881, USA

This journal is © The Royal Society of Chemistry 2025

By enabling precise interactions with biological systems at the molecular level, thermoplasmonics enable therapies such as photothermal ablation of cancerous tissue,13 bacterial eradication,14,15 and advanced drug delivery systems capable of co-delivering anticancer agents,16,17 such as RNA, DNA,18–20 and proteins.21,22 Consequently, a successful thermoplasmonic implementation requires a deep understanding of the thermal response of heated plasmonic NPs. Effective plasmonic NP-based therapy and drug delivery require optimal size and morphology for cellular uptake, high biocompatibility, and efficient heat generation while minimizing damage to healthy tissues. Effective bioincorporation requires precise NP delivery, ensuring selective accumulation. In passive targeting, the enhanced permeability and retention facilitates NP accumulation in tumors with leaky vasculature.23–25 When passive targeting is insufficient, active targeting is employed by functionalizing NPs with specific ligands for selective receptor binding.26–29 Post-delivery, adverse effects are mitigated by utilizing laser wavelengths within biological windows (700–980 nm and 1000–1400 nm),30,31 where tissue absorption and scattering are minimized due to enhanced optical transparency. Consequently, plasmonic NPs in

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Review biomedical applications primarily operate within the visible and near-infrared spectrum.32 Among the materials explored for photothermal effectbased biomedical applications, gold, silver, and copper stand out due to their LSPR-spanning wavelengths, which allow for adequate tissue penetration.33 While silver has a higher photothermal conversion efficiency than gold, both silver and copper exhibit significant toxicity and lack the chemical stability required for in vivo applications.34,35 As a result, AuNPs are preferred due to their chemical stability, low cytotoxicity, and versatile functionalization capabilities,36 i.e., their ability to undergo surface functionalization via strong sulfur–gold bonds, enhancing biocompatibility.37,38 Thiol functionalization enables AuNPs to conjugate therapeutic molecules, targeting ligands, and passivating agents, thereby improving in vivo stability and biological interactions.13,39 Thiol-functionalized AuNPs have been proposed as drug delivery vehicles, enabling targeted drug release via bond cleavage under reducing conditions.26,40–42 They have been extensively studied for nucleic acid delivery via click chemistry,42–46 a method that provides precise temporal drug release while maintaining biocompatibility and minimizing biological disruption.47,48 Diels–Alder (DA) reactions (click chemistry) are widely used to form stable cyclohexene derivatives through a reaction between a conjugated diene and a dienophile.49,50 At elevated temperatures, the reaction can reverse via the retro-Diels–Alder (rDA) pathway, regenerating the original diene and dienophile products (see Fig. 1).50,51 The thermal response of the diene/dienophile linkers can be fine-tuned to achieve controlled drug release, enabling temporal delivery of multiple drugs. Thus, understanding the thermal behavior of plasmonic NPs becomes crucial for optimizing drug delivery systems.

Md Adnan Mahathir Munshi is a PhD student in the Department of Mechanical Engineering at Penn State University. He completed his B.Sc. in Mechanical Engineering from the Bangladesh University of Engineering and Technology (BUET), Dhaka, Bangladesh. He previously served as a lecturer at the Department of Textile Engineering in Primeasia Md Adnan Mahathir Munshi University, Dhaka, Bangladesh. Now, he is on leave from his position as an Assistant Engineer (Mech) at Bangladesh Water Development Board to pursue his higher studies. His research interests lie at the intersection of Mechanics, Materials Science, Additive Manufacturing, and Machine Learning, with a particular focus on leveraging data-driven techniques to predict material behavior.

Fig. 1 Near-infrared laser light irradiation-induced photothermal heating triggers retro Diels–Alder cleavage of the surface of the PEG-DA-modified gold nanorods, releasing PEG and causing gold nanorod aggregation. (Reprinted (adapted) with permission from ref. 51. Copyright 2011 American Chemical Society).

Advancing thermoplasmonics requires a deep understanding of heat transport at the NP-solvent interface, which is quantified by the thermal boundary conductance (TBC). Thermal transport between NPs and their surrounding environment is influenced by several factors such as composition,52 size,53–56 surface properties,12,57,58 the dynamics of the solvent molecules,59 solid–liquid affinity,12 their mobility,60 and the density of covalent bonds at the NP surface61 (see Section 2 for an extended discussion). Ligand functionalization of plasmonic NPs adds complexity due to the formation of a three-component interface comprising metal, ligands, and solvent. Independent of the interface’s composition, the TBC is essential for temperature control in thermoplasmonics. Accurate local temperature measurements in nanomaterials remain difficult, underscoring the need for advanced tools to characterize heat dissipation in plasmonic NP systems. Earlier

Emdadul Haque Chowdhury is currently pursuing his PhD in Mechanical Engineering at The Pennsylvania State University. He holds a Bachelor of Science (BSc) in Mechanical Engineering from the Bangladesh University of Engineering and Technology. His research focuses on the fundamental study of nanoscale thermal, physical, and mechanical properties through atomistic level modeling, as well as the study of solid–liquid interfacial phenomena.

This journal is © The Royal Society of Chemistry 2025 View Article Online

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Nanoscale investigations used continuum heat transfer models to correlate time-dependent NP temperature changes with heat flow; these models have shown limited success.62 In contrast, atomistic simulations have become a powerful alternative, offering high-resolution insights and greater flexibility for investigating nanoscale thermal transport.63,64 As previously indicated, the spatiotemporal temperature control of plasmonic NPs depends on fine-tuning the photothermal effect and the particle-solvent heat dissipation. This review article surveys the recent literature on the latter. Section 2 reviews the literature on the fundamental physics and mechanisms governing heat transfer across solid–liquid interfaces, with a focus on AuNPs systems and interfaces functionalized with a variety of ligands and self-assembled monolayers (SAMs). Section 3 discusses recent advances in experimental techniques for characterizing thermal transport across solid– liquid interfaces, their limitations, and emerging alternatives to improve measurement sensitivity. Section 4 outlines the computational methodologies and models available for simulating functionalized gold–water interfaces, highlighting their role in predicting interfacial thermal transport. Finally, the review article concludes in Section 5 by first highlighting major challenges and research gaps, followed by a summary of key findings and an outlook on future research directions. Fig. 2 illustrates a knowledge map outlining the structure and scope of this review.

2. Thermal transport in solvated metal nanoparticles The study of thermal transport across solid–liquid interfaces and the concept of TBC trace back to Pyotr Kapitza’s investigation of the thermal conductivity of helium capillaries.65 Kapitza observed a temperature discontinuity between heated metallic surfaces and liquid helium, leading to the idea that

Dr Luis E. Paniagua-Guerra is a Postdoctoral Scholar at the Mechanical Engineering Department of The Pennsylvania State University. He received his PhD and master’s degree in mechanical engineering from The Pennsylvania State University and a bachelor’s in mechanical engineering from the University of Guanajuato in Mexico. His major interests include the fundamental study of Luis E. Paniagua-Guerra nanoscale thermal, hydraulic, and mechanical properties via atomistic level modeling, the study of solid–liquid interfacial phenomena, and the thermal characterization of liquid cooling solutions for electronics.

This journal is © The Royal Society of Chemistry 2025 Review Fig. 2 Knowledge map of the literature review on solvated AuNPs.

this discontinuity is proportional to the heat flux across any heated surface. Heat transfer across the interface between different materials inherently encounters resistance due to the abrupt change in thermal properties; this is known as thermal boundary resistance (TBR). The inverse of the TBR is the TBC, defined as J = GΔTint, where J represents the heat flux across the interface, ΔTint is the temperature discontinuity at the interface, and G denotes the TBC. During the latter half of the 20th century, the research focus shifted from solid–liquid to solid–solid interfaces, driven by advancements in microelectronics.66 In Cahill et al.’s seminal review on nanoscale heat transfer, solid–liquid heat transfer was not explicitly addressed;67 however, a follow-up review indicated a growing

Dr Jaymes Dionne is a researcher at the Naval Undersea Warfare Center (NUWC) Newport. He previously received his PhD in Mechanical Engineering from the University of Rhode Island. His research interests include the study of nanoscale heat transport phenomena through computational modeling to tune thermal transport in microelectronics and battery devices, and the development of next-generation battery materials.

Nanoscale, 2025, 17, 20803–20830 | 20805 View Article Online

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Review interest in solid–liquid interfaces between 2002 and 2012.68 Similarly, Luo and Chen69 recognized their increasing significance in biomedical applications, catalysis, energy generation, and colloidal suspensions. The early 21st century saw a renewed interest in solid–liquid TBC, particularly following the work of Ge et al.,70 who measured thermal transport at hydrophilic and hydrophobic interfaces, underscoring the need for a deeper understanding of interfacial heat transfer mechanisms. Extensive research indicates that the TBC across solid– liquid interfaces is governed by a complex interplay of factors: (i) the nature of the bonds at the interface, (ii) the interfacial liquid structure, (iii) the strength of the atomic interactions across the interface, and (iv) the vibrational mismatch between the solid and liquid phases, see Fig. 3 for a graphical depiction of these factors. Molecular dynamics simulations confirm that stronger solid–liquid coupling and surface nanostructuring enhance interfacial heat transfer by promoting the adsorption of fluid molecules.71,72 This forms an ordered, “solid-like liquid layer” that reduces the vibrational mismatch between the solid and bulk liquid.71,72 For instance, the silica–water interface has an interfacial thermal conductance (ITC) over 10 times higher than gold–water because its surface hydroxyl groups create hydrogen bonds and act as a vibrational bridge; of these factors, vibrational coupling is the most influential.73 The system temperature also plays a complex role, as the interfacial thermal resistance (ITR) for some systems is non-monotonic.74 For example, the ITR for graphene–liquid interfaces can reach a minimum value at a liquid-specific temperature (e.g., 285 K for water vs. 335 K for ethylene glycol). This behav-

Dr Ashutosh Gir is an Associate Professor of Mechanical Engineering at the University of Rhode Island (URI). He received his PhD in Mechanical Engineering from the University of Virginia in 2016. His research group utilizes laser-based optical techniques and state-of-the-art computational tools to study energy transport, conversion, and storage at the nano/macro scales. He has received several Ashutosh Giri awards, including the Office of Naval Research Young Investigator award, ASME’s BerglesRohsenow Young Investigator Award in Heat Transfer, and the Presidential Early Career Award for Scientists and Engineers (2025).

Nanoscale ior is attributed to temperature-dependent changes in the near-wall liquid density and interfacial binding energy, which alter phonon coupling at the interface.74 Despite these advancements, the mechanisms governing solid–liquid thermal transport remain only partially understood, necessitating interdisciplinary efforts to uncover the underlying principles. This complexity is further heightened in solvated NPs, where the morphology of the NP must be considered, and clustering can occur. Moreover, when organic ligands or polymers are used to functionalize solid surfaces, the dynamics of the atoms at the interface and thus the heat transfer mechanisms are modified. Accordingly, this Section is organized into three subsections: Section 2.1 reviews the fundamentals of interfacial thermal transport at solid–liquid interfaces. Section 2.2 focuses on specific effects observed in solvated NP systems, and Section 2.3 examines the additional complexities associated with characterizing thermal transport across functionalized interfaces. 2.1

Fundamentals of solid–liquid interfacial thermal transport

A consensus from early research is that the solid–liquid affinity, often quantified by the equilibrium contact angle (θc), is the key parameter controlling interfacial heat transfer. A pioneering experimental study by Ge et al.70 provided the first direct measurements of this notion, reporting TBC values in the range of 50–60 MW m−2 K−1 for hydrophobic surfaces, compared to values of 100–180 MW m−2 K−1 for hydrophilic surfaces. These observations led to the argument that hydrophilic surfaces attract a higher density of liquid molecules near the interface, thereby increasing the energy carriers avail-

Dr Bladimir Ramos-Alvarado is the Kenneth Kuan-Yun Kuo Early Career Professor, an Associate Professor of Mechanical Engineering, and the principal investigator of the Interfacial Phenomena Lab (IPHEL) at Penn State University. He is a member of the Materials Computation Center at Penn State, and an Associate of the Institute for Computational and Data Sciences (ICDS) at Penn State. Bladimir Ramos-Alvarado He also serves as a Journal Content and Subject Matter Editor for Applied Thermal Engineering (Elsevier). Dr Ramos obtained a Bachelor’s and a Master’s degree, both in Mechanical Engineering, from the University of Guanajuato, Mexico, in 2009 and 2011, respectively. He continued his education at the Georgia Institute of Technology (Georgia Tech), where he obtained a PhD and a Master’s degree in Mechanical Engineering in 2016. After brief Postdoctoral and Instructor appointments at Georgia Tech, Dr Ramos joined the Department of Mechanical Engineering at Penn State University as an Assistant Professor in May 2017.

This journal is © The Royal Society of Chemistry 2025 View Article Online

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Nanoscale Review

Fig. 3 Interfacial parameters affecting thermal transport across solid–liquid interfaces and their adaptation to AuNP–water interfaces. The central panel (with green background) demonstrates a spherical model of a functionalized AuNP surrounded by water to study interfacial heat transfer using the heat flux control method. The central panel depicts a thermal resistance network with nodes: Au, self-assembled monolayers (SAMs), and water, for the calculation of the TBC (G). The surrounding panels summarize the TBC governing parameters: surface chemistry (top-right), nanoparticle size/shape (top-left), phonon modal mismatch (left-center), interfacial bonding (bottom-left), and short-range interfacial liquid structuring (bottomright).

able for heat transfer near the interface. Consequently, TBCwettability relationships were sought after and developed. Early molecular dynamics (MD) investigations75–77 showed that the TBC increases sharply as the solid–liquid bonding strength increases, approaching a finite value in the limit of complete wetting, corresponding to a contact angle of θc = 0°. Shenogina et al.78 further identified a scaling law of the form G ∼1 + cos(θc), relating it to the work of adhesion (Wad), where Wad = γlv[1 + cos(θc)], with γlv being the liquid–vapor surface tension (see Fig. 4(a)). Subsequent investigations79–81 confirmed the G ∼1 + cos (θc) scaling law, and Alexeev et al.79 further suggested that it could be general across different interfaces. The applicability of quasi-universal TBC-wettability relationships has been challenged due to the oversight of the complex interplay of interfacial mechanisms dictating solid–liquid heat transfer. For example, Acharya et al.82 highlighted the limitations of using the contact angle on superhydrophilic surfaces and at interfaces with large curvatures, where an accurate determination of the contact angle is difficult, TBC-wettability relationships become limited (see Section 2.2 for an extended discussion on curvature effects). Furthermore, the generality of the scaling relationship G ∼1 + cos(θc) has been challenged by Ramos-Alvarado et al.,83–85 revealing two paradigm shifts: first, more wettable crystallographic planes can be less conductive; This journal is © The Royal Society of Chemistry 2025

second, in such cases, G ∼1 + cos(θc) holds independently for each plane but lacks universality (see Fig. 4(b)). Subsequent work demonstrated that factors such as chemical composition,86 crystallographic structure of the solid surface,87,88 and interface curvature56 can undermine the quasi-universal nature of the TBC-wettability relationship. These contradictions to the early TBC-wettability relationship highlight the intricate nature of heat transfer across solid–liquid interfaces, emphasize the need to look beyond the solid–liquid interaction strength, and underscore the importance of exploring additional factors influencing TBC behavior. The nature of interfacial interactions, particularly noncovalent forces, plays a crucial role in determining the TBC. These interactions are primarily governed by electrostatic interactions, including Coulombic attraction and polarization effects, as well as the formation of hydrogen bonds between the solid and liquid. For instance, several investigations have reported an increased TBC across various solid–liquid interfaces due to polarization effects,89–93 which are crucial for metallic NPs immersed in biological environments.94 In these cases, polarization causes negligible changes in wettability and interfacial free energy, but it influences the molecular ordering of the liquid phase, which has been attributed to enhancing the TBC of polar interfaces due to favoring the formation of hydrogen bonds (H-bonds).86,91,92,95 The main- Nanoscale, 2025, 17, 20803–20830 | 20807

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Review Nanoscale

Fig. 4 Descriptions of the TBC using wettability. (a) Early contributions suggested a quasi-universal law of the form G ∼1 + cos(θc), from a relationship to the work of adhesion Wad = γlv[1 + cos(θc)] (Reproduced from ref. 81. With the permission of AIP publishing). Challenges to the TBC-wettability notion. (b) TBC across silicon- and graphene-coated–silicon–water surfaces. (Reprinted (adapted) with permission from ref. 83. Copyright 2016 American Chemical Society.).

stream perspective suggests ordered H-bonds pull liquid molecules closer to the solid surface, enhancing the TBC. Alternatively, the TBC enhancement in polarizable interfaces has been explained by the excitation of additional degrees of freedom, such as vibrational modes in polar solvents.89 Nevertheless, it has been reported that polarizability negligibly modifies the vibrational density of states (vDOS) of the metal or liquid phases.92 In such cases, the enhancement in TBC has been attributed to a rise in the phonon transmission probability at the interface, which is sensitive to molecular ordering and interatomic spacing. Vibrational mode mismatch between solid and liquid particles is a key factor governing phonon-mediated thermal transport across interfaces, where the overlap of the phonon density of states (DOS) between the two phases quantifies this effect.96–100 Giri and Hopkins101 used simple Lennard-Jones (LJ) solid–liquid MD models to show that stronger solid–liquid bonding enhances low-frequency phonon coupling, broadens the interfacial DOS, and introduces new phonon modes; thus, increasing the TBC. In contrast, weak or hydrophobic interfaces act like free surfaces, limiting phonon transmission. Han et al.102 reported a similar shift of vibrational modes to higher frequencies in perfluorohexane, though driven by increased liquid pressure rather than bonding. Surface functionalization, like self-assembled monolayers (SAMs) and chemical passivation, can reduce the vibrational mismatch by introducing buffer interfacial modes, improving phonon overlap even when the bulk DOS differs.103–106 However, the relationship between modal overlap and TBC is not universal; interfacial liquid structuring and the directionality of heat flux (in-plane vs. out-of-plane) also influence transport. Out-of-plane modes dominate at low-affinity interfaces, while strong bonding and ordered structuring (e.g., via hydrogen bonding or electrostatics) enhance in-plane contributions.107–111 These findings highlight the need to consider vibrational mismatch, interfacial chemistry, and liquid ordering for a comprehensive understanding of nanoscale heat transfer across solid–liquid interfaces.

The role of the liquid’s molecular interfacial organization in the TBC has been extensively investigated. Several authors have shown that liquid layering at the interface is essential in determining the TBC. The prevailing hypothesis is that since molecules are the primary energy carriers in liquids, their availability and proximity to the interface significantly impact the energy transfer probability. Furthermore, microcalorimetry and heat capacity measurements have revealed that absorbed water on metal oxide surfaces exhibits distinct thermodynamic properties compared to bulk water,112 suggesting enhanced thermal properties for interfacial liquids. Early MD investigations linked the TBC to the height and location of the first hydration layer near the interface,79,90,91,113 showing that higher TBC values correlate with enhanced liquid layering, although this relationship is non-universal.91 Subsequent works investigated the complex liquid layering that extends beyond the first hydration layer,114,115 while also accounting for interfacial pressure effects115 and the formation of solidlike structures in the liquid phase.116 Building on the observed dependence of TBC on interfacial liquid layering, RamosAlvarado et al.83 used the density depletion length δ as a parameter to reconcile the anisotropic TBC calculations of silicon– water interfaces as illustrated in Fig. 5(a). More recently, Motokawa et al.118 introduced the radial density depletion length (RDDL) to account for single-atomic structures on solid surfaces, demonstrating that liquid ordering significantly modulates interfacial thermal transport as illustrated in Fig. 5(b). The concept of δ, which quantifies the deficit or surplus of liquid molecules near the interface, has also been employed to describe hydrodynamic slip.119–122 Subsequent contributions validated δ as a reliable parameter for describing the TBC across different interfaces.86,88,111,117 Recent work by Paniagua et al.117 further underscores the importance of interfacial liquid organization, showing that the TBC is significantly enhanced when liquid molecules near the interface organize into cluster-like structures, as opposed to uniform, layered arrangements. These cluster-like structures refer to localized, high-density regions of water that lack long-range

This journal is © The Royal Society of Chemistry 2025 View Article Online

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Nanoscale Review

Fig. 5 (a) Reconciliation of the anisotropic TBC calculated for silicon surfaces in contact with water using the density depletion length δ. (Reprinted (adapted) with permission from ref. 83. Copyright 2016 American Chemical Society). (b) TBC as a function of the radial density depletion length ξ. (Reprinted (adapted) with permission from ref. 118. Copyright 2024 American Chemical Society.) Density contours of flat Au–water interfaces, comparing (c) layered and (d) cluster-like organizations; white regions indicate interfacial zones where water is fully excluded. (Reproduced from ref. 117, with the permission of AIP Publishing.).

lateral order and are instead confined to irregular, spatially heterogeneous domains, often influenced by variations in ligand chemistry or surface affinity. A visual comparison between layered and cluster-like organization is provided in Fig. 5(c) and (d), respectively. Recent studies have expanded the understanding of heat conduction at the nanoscale beyond purely diffusive mechanisms. The observation of second sound in graphite at temperatures exceeding 100 K demonstrates that collective phonon transport becomes relevant even in systems where diffusive models have traditionally been assumed to apply.123 This challenges the universal applicability of Fourier’s law and suggests that non-Fourier effects may also emerge in systems with confined geometries and strong vibrational coupling. In parallel, investigations into ultrathin coatings have revealed that even a single atomic layer can substantially impact interfacial thermal resistance. For example, the presence of a monolayer of graphene at a Cu–water interface was shown to increase the Kapitza length by a factor of 2.5, despite preserving the macroscopic wettability of the surface.124 This finding underscores

This journal is © The Royal Society of Chemistry 2025 the importance of considering vibrational mismatch and interfacial structure, even when coatings are only a single molecule thick. In functionalized AuNP systems, several features may influence these non-classical transport regimes, e.g., the presence of structured water layers, ligand–water hydrogen bonding, modal mismatch attenuation, and temperature-sensitive chemistry at interfaces. Additionally, transport in ligands and across adsorbed water layers may exhibit dominant ballistic or hydrodynamic phonon transport. Surface functionalization further modulates phonon scattering processes primarily through bonded interactions at the gold–ligand interfaces and non-bonded interactions, depending on the chemistry and structure of the interface. These factors highlight the need for future studies that integrate non-equilibrium and nondiffusive frameworks to better characterize heat transport in solvated nanoparticle systems. 2.2

Thermal transport measurements across individual NPs are challenging; thus, several experimental contributions have Nanoscale, 2025, 17, 20803–20830 | 20809 View Article Online

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Review focused on macroscale insights into the thermal relaxation of colloidal NP solutions, also known as nanofluids. These works have revealed dependencies on the intensity of optical excitation125 and NP concentration.126,127 Additionally, enhancements in the thermal conductivity of nanofluids have been shown to depend on NP size.128 To understand the mechanisms driving these enhancements, MD models have focused on computing the effective thermal conductivity (ETC) of nanofluids. Sarkar and Selvam129 demonstrated that while the diffusion of NPs is slower than that of liquid atoms, the interfacial liquid atoms surrounding the NPs exhibit enhanced movement compared to bulk liquid atoms, which contributes to increased ETC. Further research has examined the impact of NP aggregation, reporting that aggregation enhances the ETC of nanofluids, with chain-like NP aggregates providing a greater increase in ETC compared to spherical aggregates.130,131 Liquid layering at NP interfaces has also been explored, suggesting that the local ETC of liquid adsorption layers increases with NP wettability,132 shedding light on the role of interfacial phenomena in thermal transport within nanofluids. A deeper understanding of thermal transport across curved interfaces can be achieved by focusing on the TBC at the NP interface rather than the ETC. Fundamentally, the cooling dynamics of rapidly heated NPs can be estimated by assessing the TBC and the NP’s size. In their experimental work, Ge et al.59 showed that for sufficiently large spherical NPs or high TBC interfaces, the NP’s temperature decay is limited by thermal diffusion in the surrounding fluid. The characteristic diffusion time (τd) can be estimated by equating the particle’s heat capacity with that of the fluid within a thermal diffusion length. Conversely, for sufficiently small NPs or low TBC interfaces, the cooling rate is limited by the TBC, with the characteristic decay time (τi) determined by the ratio of the particle’s

Nanoscale heat capacity to the total interfacial thermal conductance. Based on this, Ge et al.59 proposed a critical TBC value, TBCc ¼ 3Cf kf rp Cp ð1Þ

where kf is the thermal conductivity of the fluid, Cf and Cp are the volumetric heat capacities of the fluid and NP, respectively, and rp is the NP radius. This formulation effectively demarcates two distinct cooling regimes: when TBC ≫ TBCc, the cooling is diffusion-limited; when TBC ≪ TBCc, it is interfacelimited. More recently, Wilson et al.133 proposed an alternative definition of the critical TBC, grounded in the concept of the Kapitza length (effective length creating the same thermal resistance as an interface). They defined the critical conductance for water-solvated particles as, Gc ¼

where kw is water’s thermal conductivity, and d is the NP’s diameter. Similar to Ge et al.’s59 conclusions: If G ≫ Gc, the NP’s cooling process is dominated by water diffusion; alternatively, if G ≪ Gc, the NP’s cooling is controlled by the interface. Wilson et al.133 defined the diffusion-dominated regime as G > 10Gc, the interface-dominated regimes as G < 0.1Gc, and a mixed regime as 0.1Gc < G < 10Gc. These regions, along with a survey of experimental and computational data, are plotted in Fig. 6(a), where it can be observed that most AuNP systems exist in the mixed regime and skew towards the interface-dominated area; thus, supporting the need for further research on interfacial heat transfer in solvated AuNPs. A 10 nm AuNP in water can be considered as an example to contrast the two critical TBC models. In this situation, critical TBCs of 300 MW m−2 K−1 and 100 MW m−2 K−1 are obtained

Fig. 6 (a) Impact of the TBC on the nanoparticle’s spatiotemporal temperature regulation (Reproduced from ref. 131, with the permission of AIP Publishing). (b) Normalized interfacial thermal conductance of nanoparticles as a function of the nanoparticle radius. The normalized conductance ðmaxÞ (G̃k) corresponds to the conductance at a given radius and contact angle, Gk(Rnp,θ), divided by the maximum conductance value, Gk (θ), for that specific wetting condition. (Reproduced from ref. 135, with the permission of AIP Publishing.).

This journal is © The Royal Society of Chemistry 2025 View Article Online

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Nanoscale using eqn (1) and (2), respectively, and while these numbers differ by a factor of three, they exist in the same cooling regime per Fig. 6(a). The difference lies in the fact that eqn (1) was derived by obtaining the ratio of the thermal time constants in the particle and surrounding liquid, i.e., transient heat transfer parameters, and eqn (2) was derived using a Kapitza conduction length analogy. Notably, Ge et al.’s59 model is more conservative if one were to follow the same mapping strategy as Wilson et al.133 depicted in Fig. 6(a), underscoring the need for a deeper fundamental understanding of interfacial heat transfer in solvated NP systems. NP size and curvature are synonyms of the same parameter that can be computationally investigated. Merabia et al.134,135 used MD simulations of solvated AuNPs to demonstrate that curvature significantly alters the thermodynamic properties of the interfacial liquid. Due to their curved geometry, spherical AuNPs could be heated above their melting temperature without causing a phase change in the adjacent liquid. Additionally, a vapor layer, which typically forms on flat interfaces under similar heating conditions, was notably absent at curved spherical interfaces. The delay in liquid phase change and AuNP melting was attributed to the extremely high pressure near the curved interface, i.e., the Laplace pressure generated by the NP’s curvature. Later, in their computational work on nanoscale boiling around AuNPs, Gutiérrez-Varela et al.136 demonstrated that the formation of a low-density liquid layering during heating transiently reduces the TBC and delays vapor nanobubble onset at the AuNP–water interface, which also explains the absence of interfacial water phase change reported by Merabia et al.134 Notably, when evaporation was reached, it was reported that nanobubbles nucleate more rapidly on hydrophilic nanoparticles, contrary to predictions from isothermal classical nucleation theory. Merabia et al. and Gutiérrez-Varela et al.134–136 contributions demonstrate the potentially devastating effects of poor spatiotemporal temperature control of AuNP therapies. While particle melting and water nucleation could be delayed due to the large Laplace pressure around spherical NPs, evaporation is still plausible, and its subsequent effects are lower TBCs and eventually potential AuNP melting. Thus, particle size effects and interfacial liquid structure properties must be better understood to engineer AuNPs’ temperature controls. The computational work by Tascini et al.56 was one of the first to systematically establish a direct relationship between NP curvature and the TBC. Using a generic nanoparticle-fluid model, they demonstrated that the TBC increases with interfacial curvature across a wide range of fluid-solid interaction strengths. Their findings revealed an empirical relationship described by G = G∞ + c/r, where 1/r is the NP’s curvature, G∞ is the TBC in the limit of r going to infinity or a flat surface, and c is a fitting parameter. They observed that stronger interfacial interactions lead to larger values of c. Building on this, Gutiérrez-Varela et al.137 investigated the impact of curvature and size on the TBC of AuNPs across three distinct wetting regimes: strong, intermediate, and weak, as illustrated in Fig. 6(b). Their calculations matched the curvature effect of

This journal is © The Royal Society of Chemistry 2025

Review Tascini et al.’s56 empirical correlation, but using a realistic metal–liquid system. Gutiérrez-Varela et al.137 explained the enhanced conductance of smaller AuNPs using two arguments. First, smaller AuNPs have higher solid–liquid coordination numbers (a greater number of water molecules per surface atom), which creates a higher water–Au potential energy. Second, the NP’s curvature alters the interfacial vibrational spectrum: as the AuNP size decreases, the high-frequency van Hove peak fades while the low-frequency peak strengthens, aligning Au and water vibrations more closely and enhancing the TBC. Additionally, they observed that for smaller NPs, the amplitude of the first peak in the water density profile increases, consequently enhancing the structuring of interfacial water. Yet, they caution that the correspondence between interfacial conductance and fluid density is not universal. Expanding on these insights, Paniagua-Guerra and RamosAlvarado117 investigated interfacial heat transfer at AuNP– water interfaces, emphasizing the role of the density depletion length (δ). Their MD simulations demonstrated that curved interfaces consistently exhibited higher TBC than flat surfaces. This enhancement was attributed to the larger availability of water molecules at the interface, which facilitated energy transfer. Additionally, they identified an exponential relationship between the TBC and δ, indicative of the transferability of the TBC-δ relationship to curved interfaces, where traditional wettability metrics are difficult to compute. The influence of NP morphology on TBC has been further explored by studying NPs with various shapes. Neidhart and Gezelter138 dispersed bare AuNPs, icosahedral, cuboctahedral, and spherical, in solvent and systematically examined how NP morphology influences the TBC by quantifying the density of undercoordinated sites on the solid surface. They observed higher TBC values for particles with a greater fraction of exposed undercoordinated atoms. Building on this concept, Jiang et al.139 showed that TBC can vary locally across an NP’s surface: solid atoms with lower coordination numbers i.e., fewer neighboring atoms make more contact with the solvent, enhancing local heat transfer. Similarly, Gutiérrez-Varela137 quantified the number of water molecules interacting with a surface gold atom, via the water–gold potential energy, and demonstrated that decreasing the NP’s size increases this number, thereby enhancing the TBC. These findings emphasize the critical role of NP shape and atomic coordination in determining interfacial thermal transport properties. In summary, the strong dependence of the TBC on morphological factors highlights the complexity of describing interfacial thermal transport across solid–liquid interfaces. This complexity extends beyond a simple characterization of interfacial bonding strength, as the energy landscape at the interface is influenced by both the solid surface morphology and the strength of interfacial interactions. This complexity is compounded by the choice of computational models used to simulate interfaces. For example, the treatment of surface polarization in molecular simulations can dramatically impact the calculated TBC.140 It has been demonstrated that polarizable

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Review force fields, such as those using Drude oscillators, can introduce artificial vibrational couplings between the model’s internal modes and the librational modes of water, leading to a significant overestimation of the TBC.140 Crucially, two different gold–water models predicting the same interfacial tension (i.e., the same wettability) can yield vastly different TBCs, demonstrating that ITC is not directly correlated with interfacial free energy alone and depends heavily on the specific vibrational landscape of the model.140 The surface affinity-TBC analysis becomes increasingly intricate as the interface structure grows more complex. However, the underlying mechanisms and physics governing interfacial thermal transport remain consistent. Consequently, much of the knowledge gained from thermal transport across bare solid–liquid interfaces can be applied to functionalized solid–liquid interfaces, as will be explored in detail in the following Subsection. 2.3 Thermal transport at functionalized solid–liquid interfaces Interfacial heat transport in functionalized AuNPs is strongly influenced by the ligand type and surface morphology and can be easily represented by the resistance network panel in Fig. 3. The dependence of interfacial thermal transport on ligand– solvent interactions was initially highlighted by Ge et al.59,70 who demonstrated that the affinity between the solvent and the terminal groups of ligands significantly influences the TBC. More recent investigations have further demonstrated that the pronounced enhancement in TBC arises from tight coupling at two interfaces: from the metal core to the ligand shell, and from the ligands to the surrounding fluid. For instance, Au surfaces functionalized with polyethylene glycol (PEG) exhibit significantly higher TBCs in water compared to those functionalized by citrate or cetyltrimethylammonium bromide (CTAB) ligands.141 This enhancement arises from the strong Au–S bonds that couple the Au core to the PEG ligands, adding to the increased physical contact between those ligands and the surrounding solvent. Similarly, the presence of a ligand layer that structurally/chemically matches the solvent can create a buffer layer that reduces the vibrational modal mismatch at the interface.142 Lastly, the ligand surface coverage could increase the availability of channels for interfacial conduction.143–145 These three individual effects on the TBC will be discussed further in this section. For heat to flow from a functionalized solid to a solvent, efficient transport occurs from the solid to the ligands due to strong covalent chemical bonds. This is followed by heat transfer from the ligands to the liquid solvent, which could be enabled via vibrational coupling.142,145–148 The ligand layer may act as an intermediary, bridging the solid and liquid phases, which typically exhibit significant vibrational mismatches.143,145,147 For instance, Kikugawa et al.145 demonstrated via vibrational analysis that self-assembled monolayers on Au significantly reduce the interfacial thermal resistance compared to bare gold–solvent interfaces. Similarly, Hannah and Gezelter143 showed that hexylamine ligands enhanced vibrational overlap between CdSe and hexane (the surrounding

Nanoscale solvent), thereby facilitating improved interfacial heat dissipation. Contrariwise, Hung et al.104 reported a negligible phonon spectral overlap effect in SAM-coated gold and water, where better vibrational coupling did not correlate with higher TBC. Instead, they found that thermal transport is primarily facilitated by the aggregation of water molecules around the terminal atoms of SAM. Thus, depending on the vibrational properties of both the liquid solvent and the ligand layer, the ligand–liquid interface can exhibit either the largest146 or the smallest145 thermal resistance within the three-component solid–ligand–liquid interface as depicted in Fig. 3. This “phonon bridge” effect was demonstrated in simulations of a gold–pentacene (organic semiconductor) interface functionalized with SAMs.149 It was found that SAMs effectively connect the low-frequency phonon density of states of gold with the disparate vibrational modes of the organic material, creating new energy transport channels that are absent in a bare interface.149 Similar to bare solid–liquid interfaces, non-bonded interactions between the ligands and solvents affect the TBC.146,150 For polar interfaces, electrostatic forces promote the formation of stable hydrogen bonds at the ligand–liquid interface.151 Stronger hydrogen bonding draws polar organic solvent molecules closer to the interface, resulting in tighter molecular packing. This intermolecular attraction requires that the organic molecules, either on the solid side or in the solvent, contain the necessary functional groups with highly electronegative atoms like oxygen or nitrogen.152 The increased proximity, along with a higher atomic number density of the organic liquid near the interface, facilitates thermal energy transport. The significant impact of terminal group chemistry was systematically shown for Au–pentacene interfaces, where the TBC was enhanced by 6–7 times using SAMs with nonpolar –CH3 and –NH2 groups, but by 11 times when using a highly polarized -COOH terminal group.149 This superior enhancement was attributed to stronger interfacial affinity, evidenced by higher adhesion energies and the formation of hydrogen bonds at the interface, which pull the adjacent molecules closer and provide additional pathways for energy transfer.149 The strength of ligand–liquid interactions can be tailored by modifying the chemical composition of the functional groups on the ligand layer.82,153 Shavalier and Gezelter154 investigated the influence of ligand-to-solvent hydrogen bonding on heat transfer and demonstrated that PEG-capped AuNPs in water exhibit enhanced thermal conductance. Their analysis of vibrational power spectra revealed an increased population of low-frequency heat-carrying modes (0–70 cm−1) for the thiolated PEG. Because of the Bose–Einstein weighting of lower frequency modes, improved thermal transport was observed. Their findings also suggest that solvent penetration and ligand configuration-specifically, the orientational ordering of ligand chains-play crucial roles in interfacial heat dissipation. Alternatively, experimental measurements by Tian et al.153 demonstrated that the TBC is insensitive to the ligand’s chain length, suggesting that interfacial transport at the ligand–water interface is primarily dictated by the chem- This journal is © The Royal Society of Chemistry 2025

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Nanoscale istry of the terminal groups on Au surfaces exposed to water molecules. These findings, highlighting the influence of chain length and solvent penetration, have been further corroborated by MD simulations, as demonstrated in the work by Stocker and Gezelter, examining thiolate-capped gold surfaces.155 Computational investigations on surface ligand coverage have shown that partly covered surfaces enhance the TBC, as illustrated in Fig. 7.143,148 This enhancement is attributed to the increased number of thermal exchange pathways and improved vibrational coupling between the hexylamine ligands-passivated CdSe surfaces and the surrounding hexamine solvent. However, at near 100% surface coverage, a critical turning point is reached at which TBC begins to decrease.156 This reduction occurs because excessive ligand coverage prevents effective penetration of liquid molecules into the ligand layer. As surface coverage continues to increase, the reduced mobility of liquid molecules within the densely packed ligand layer hinders interfacial heat transfer, leading to diminished TBC. Alternatively, Zhang et al.103 demonstrate that decorating interfaces with high-coverage polymeric SAMs can significantly enhance the TBC even between materials with considerable vibrational mismatch. Specifically, they reported a 430% increase in the TBC after coating both sides of graphene with 7.14% polyethylene (PE), using polymethyl methacrylate (PMMA) as the surrounding medium. This enhancement was attributed to three key factors: (i) the formation of extended and well-aligned polymer chains within the PE/PMMA blending region, (ii) strong vibrational coupling between PE and PMMA, and (iii) covalent bonding between graphene and PE chains. The role of liquid mobility at functionalized solid–liquid interfaces in determining TBC remains a topic of ongoing debate. Some studies have shown that reduced mobility of interfacial liquid molecules, particularly water, can enhance the TBC in systems such as AuNPs functionalized with organic ligands.157 This enhancement is often attributed to improved vibrational coupling between the ligand layer and liquid mole-

Review cules, which facilitates phonon transmission across the interface.142,155 For example, when hexene molecules align with thiolate chains, thermal transport improves due to stronger vibrational overlap. However, other studies have reported the opposite effect. Low liquid mobility can hinder heat transfer when molecules become trapped or immobilized at the interface, limiting energy exchange through molecular diffusion.155 These conflicting findings indicate that the influence of liquid mobility is not yet fully understood, and it is likely solvent–ligand pair specific. Therefore, further research is needed on how the properties of the thiolate layer, such as surface coverage and chemical affinity with the solvent, affect liquid mobility and interfacial heat transfer. In addition to mobility, the organization of liquid molecules at the interface and their proximity to both the solid surface and the ligand layer play a critical role in interfacial thermal transport. Liquid molecules in closer proximity to the solid surface and the ligand layer facilitate more effective energy exchange at the interface.150 However, the role of liquid layering and structuring in thermal transport across functionalized solid–liquid interfaces remains debated. Neidhart and Gezelter158 found that a higher solvent density peak within the penetration region correlated with an increased TBC. Conversely, Sun et al.147 observed a weaker dependence of the TBC on liquid layering for gold slabs coated with SAM. They concluded that liquid layering effects are more pronounced in bare solid–liquid interfaces than in functionalized ones. This highlights the nuanced and context-dependent role of liquid structuring in thermal transport across functionalized interfaces. Analyzing heat transport across functionalized solid–liquid interfaces is inherently complex due to the coexistence of a three-component interface. Unlike bare solid–liquid interfaces, additional factors must be considered when evaluating the solid–ligand–solvent layer. (i) Ligands exhibit localized thermal motion due to vibrational and conformational fluctuations, unlike rigid solid atoms; however, they lack the transla-

Fig. 7 (a) Interfacial thermal conductance, G, as a function of ligand density on various CdSe surfaces. (b) Different CdSe surfaces considered in ref. 141. (Reprinted (adapted) with permission from ref. 141. Copyright 2015 American Chemical Society.).

This journal is © The Royal Society of Chemistry 2025 Nanoscale, 2025, 17, 20803–20830 | 20813 View Article Online

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Review Nanoscale

tional mobility of liquid molecules and do not undergo diffusion.144 (ii) The ligand–water interface is not well-defined because water molecules can penetrate the ligand layer.143,157 This penetration significantly influences the ligands’ temperature profile.157 For bare or hydrophilic ligand-coated interfaces, the temperature profile typically shows a single steep descent at the solid–ligand interface. However, for hydrophobic ligand-coated interfaces, the temperature profile becomes more complex, exhibiting an initial drop at the solid– ligand interface, followed by a plateau along the ligand, and finally, a second drop at the liquid–ligand interface.157 This intricate behavior underscores the need for detailed analysis to understand thermal transport in such systems. Computing the TBC at a complex solid–ligand–liquid interface requires a simplified model to account for its intricacies. A common approach involves calculating a global TBC by considering the temperature change from the solid surface to an idealized sharp ligand–liquid interface.143,145,155 This method reduces the three-component interface into two independent interfaces: the well-defined solid–ligand interface and the diffuse ligand–water interface, which is approximated as a sharp boundary. Alternatively, some authors employ an effective thermal resistance model (the inverse of TBC), which represents the interface as a network of smaller thermal resistances (as illustrated in Fig. 3). These resistances are defined by the discrete temperature jumps observed at different points across the interface.148,155 This approach allows for a more detailed characterization of thermal transport mechanisms at the interface. Major findings related to TBC modeling efforts and their implications in biomedical fields are summarized in Table 1.

3. Experimental measurements of thermal transport across solid–liquid interfaces In the past two decades, there has been considerable advancement in the theoretical understanding of interfacial thermal

Table 1 transport across solid–liquid interfaces, which has mainly been fueled by the tremendous progress in the atomistic modeling based on MD simulations of several solid–liquid interfaces,75,78,113,116,159–164 and the advances achieved in the analytical description utilizing the phonon theory of liquid thermodynamics.165–168 Although comparatively there have been fewer experimental works focusing on understanding thermal transport across solid–liquid interfaces,70,81,153,169–174 a handful of these contributions have provided crucial validation to the theoretical advancements. For instance, Harikrishna et al.81 have shown that by varying the terminal of the alkane–thiol monolayers on the gold thin film surface, the thermal conductance values monotonically increased in the range of 60–190 MW m−2 K−1 as the work of adhesion increased (Fig. 4(a)), and the solid–liquid contact angles spanned from 25° to 118°. The measurements were carried out with the time-domain thermoreflectance technique (TDTR), which utilizes a femtosecond pulsed laser system to monitor (in real-time) the temperature changes on a metallic surface induced by the laser pulse absorption.81,175,176 In fact, such pump–probe laser-based techniques are particularly well suited for this application and have been the popular choice for investigating interfacial heat transfer across solid–liquid interfaces. In the pump–probe-based thermoreflectance techniques, the laser pulses are absorbed by the solid (usually Au thin films deposited on a transparent substrate), and a bidirectional heat flow model is used to back out the TBC across metal–liquid interfaces (as schematically represented in Fig. 8(a)). The popular choice for the liquid has been water, and to vary the interfacial adhesion (hydrophobicity), the wellknown gold–thiol chemistry (as utilized by Harikrishna et al.81) is utilized. However, the use of the traditional thermoreflectance techniques, which have been the current standard for measuring thermal boundary conductance, lacks sufficient sensitivity to accurately quantify the interfacial heat conduction across solid-liquid interfaces.169 The experimental insensitivity to solid–liquid TBC, in general, originates from the large thermal resistance posed by

Overview of research on TBC across Au-functionalized interfaces pertinent to biomedical applications Source Shavalier et al. 141 Methodology Key findings Implications Computational – MD simulations

• Surface morphology significantly impacts interfacial heat transport in functionalized AuNPs. • Stronger ligand–water interactions promote better coupling, enhancing TBC at the gold–water interface. • Ordered structuring of interfacial water molecules enhances thermal conductivity.

• Ligand selection can improve heat dissipation and drug delivery efficiency. • Tailoring ligand chemistry can fine-tune interfacial heat transfer properties for biomedical applications. Li et al.95 Computational – MD simulations

Tascini et al.56 Computational – MD simulations Experimental • Curved AuNPs exhibit altered vibrational modes that influence interfacial thermal conductance. • Ligand-to-water hydrogen bonding enhances thermal coupling, enhancing TBC.

Computational – MD simulations • Higher surface ligand density initially improves TBC, but excessive coverage inhibits effective energy transfer. Ge et al.59 Hannah et al.143 20814 | Nanoscale, 2025, 17, 20803–20830

• Advanced interfacial models should incorporate liquid structuring effects for accurate TBC predictions. • Curvature effects should be considered in AuNP design for targeted heat transport applications. • Functionalization strategies should focus on maximizing ligand-water interactions to enhance heat transfer. • Balancing ligand surface coverage is crucial for optimizing thermal response without hindering drug release.

This journal is © The Royal Society of Chemistry 2025 View Article Online

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Nanoscale Review

Fig. 8 (a) Schematic of the three-layer system (usually incorporating glass/Au/water layers) utilized in typical TDTR measurements, where the pump and probe beams are incident on the gold surface after passing through the transparent substrate. (b) Measurement of the TBC (hk) lower bounds for a thin film Au/H2O and Au/ethanol interfaces. (c) Sensitivity of the solid–liquid TBC in comparison to the thermal conductivity of the liquid in typical TDTR analyses.

the liquids relative to that of the interfacial heat flow.177 This has been shown quantitatively through calculations of the sensitivity (based on TDTR sensitivity analysis carried out by Costescu et al.178) of the typical TDTR signal (representing the ratio between the in-phase and the out-of-phase signals) to the various parameters in the thermal model used to back out the thermal boundary conductance (Fig. 8(c)). The relative sensitivity of the solid–liquid TBC is significantly lower in comparison to the thermal conductivity of the liquid due to their low thermal diffusivities. Recently, Tomko et al.169 highlighted the lack of sensitivity of the typical pump–probe experiments to solid–liquid interfacial heat flow by using gold films in contact with several different liquids. In this work, the pump and probe beams were passed through transparent glass substrates and focused on the surface of gold films in contact with various liquids at the other end. Similar to the conventional approach, a bidirectional heat flow model was used to back out the thermal boundary conductance. The results from the measurements were compared with TBC measured across other gold–substrate interfaces. In this regard, it is instructive to compare the values as a function of the ratio of the longitudinal sound velocities for the two media comprising the interface, as illustrated in Fig. 8(b). Assuming a simple acoustic mismatch model (AMM) for thermal boundary conductance,179 increasing overlap of the sound velocities, ν, of the two media is expected to enhance the heat conduction and reduce the temperature drop that occurs at the interface. Although the comparison of the measured values suggested that the TBC across gold–liquid interfaces can be as high as those of other solid–solid interfaces associated with gold films, these measurements represented lower bounds (with the error bars representing a 5%

This journal is © The Royal Society of Chemistry 2025 error in the film thicknesses). For these gold–liquid interfaces, it was not possible to obtain a nominal value for the upper bound of the measured TBC with the conventional TDTR technique alone, as the experimental insensitivity to solid–liquid TBC originated from the large thermal resistance posed by the liquids. Although it has been difficult to accurately determine the TBC with the typical thermoreflectance techniques, the work by Tomko et al.169 showed that alternative pump–probe experiments to quantify nanoscale energy transport at solid–liquid interfaces can support the traditional measurements and provide the much-needed validity. Namely, they probed the damping of acoustic phonon modes (commonly referred to as picosecond acoustics) in the solid layer upon interaction with the solid–liquid interface. This technique is based on the optical detection of the propagation of acoustic modes through the piezo-optic effect and can provide information on the transmissivities of acoustic phonons across the interface between the solid film and contacting layers.169 In other words, the ultrafast pump pulse excitation of metal films produces an oscillatory strain wave that travels in the thin film and interacts at the metallic film interface, which attenuates the oscillatory strain, and thus allows for the measurement of the phonon mode transmissivity across that interface.180–182 Tomko et al.169 showed that the transmissivities increased with the increase in the work of adhesion at the solid–liquid interfaces. They further supported their measurements with experiments that monitored the ablation threshold for the various samples, which served as a metric for changes in thermal transport at the gold–liquid interfaces. More specifically, the ablation threshold for gold thin films in contact with different liquids was shown to correlate well with TDTR measurements of TBC across the Au–liquid interfaces.183 Note

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Review that the ablation threshold is a quantitative measure of the minimum laser fluence required to remove mass from the thin gold films in the experiments. Notably, in Tomko et al.’s work, this threshold was highly dependent on the TBC between the Au films and contacting layers, where the ablation threshold increased linearly with increasing TBC. By correlating this linear increase with the ablation thresholds measured at the various gold–liquid interfaces, it was shown that the TBC across the Au–liquid interfaces could be determined with lower uncertainties as compared to the regular TDTR measurement analysis procedure with large uncertainties.169 Jiao et al.184 introduced an innovative experimental approach using the 3ω technique to evaluate the TBR between water and superhydrophobic surfaces in the Cassie state, where water sits atop structured surfaces with air gaps underneath. Their findings reveal that the presence of air at the interface significantly increases the TBR due to its poor thermal conductivity. To accurately measure this effect, they employed a combined differential and bi-directional 3ω method, enabling precise characterization of heat transfer at the solid–liquid interface. This approach offers enhanced sensitivity compared to conventional thermoreflectance techniques, which are generally inadequate for such interfaces. Therefore, although current thermoreflectance techniques lack the required sensitivity to accurately determine the TBC across solid–liquid interfaces, a combination of such pump–probe experiments with other sensitive techniques like the 3ω method can provide a fundamental platform for advancing our fundamental understanding of interfacial heat flow in these systems.

4. Computational modeling of solvated gold nanoparticles Molecular dynamics (MD) is a powerful computational tool used to model atomic-level interactions using a classical physics framework. In the field of heat transfer, MD facilitates the calculation of the thermal conductivity and TBC. Furthermore, the atom-level resolution offered by MD enables a deeper understanding of the underlying mechanisms of thermal transport. Achieving meaningful results from MD simulations requires atomic structures to ideally be derived from experimental data or principles of crystallography to accurately represent a solid. Additionally, it is of paramount importance to select or develop a suitable set of force fields (FFs) that faithfully capture the interactions within the system and align with the specific research objectives. This section provides guidance for proper MD modeling of functionalized gold–water interfaces, which is critical for advancing the application of plasmonic nanoparticles in the biomedical field. Accordingly, this section is organized into three parts: Section 4.1 reviews various MD methods for calculating interfacial thermal transport at solid–liquid interfaces. Next, Section 4.2 delves into relevant literature covering MD foundations on modeling different components of bare and functionalized gold water interfaces. Finally, Section 4.3 offers a concise over- 20816 | Nanoscale, 2025, 17, 20803–20830

Nanoscale view of the models in two subsections: Subsection 4.3.1 summarizes the extensive work on Au–water interactions, and Subsection 4.3.2 covers models of thiolate adsorption on Au surfaces, including mathematical representations of the Au– sulfur bond. 4.1 Molecular dynamics methods for thermal boundary conductance calculations The TBC is crucial for characterizing nanoscale heat transfer, and MD simulations offer several methods to calculate it. Rajabpour et al.185 compared four primary MD techniques for evaluating the TBC at nanoparticle–water interfaces: the transient non-equilibrium molecular dynamics (TNEMD) method using both lumped capacitance and finite internal resistance models, the steady non-equilibrium molecular dynamics (SNEMD) method, and the equilibrium molecular dynamics (EMD) approach. TNEMD exists in two versions: (i) the lumped capacitance model, which assumes negligible internal conduction resistance; and (ii) the finite internal resistance model, which accounts for internal temperature gradients. The lumped model closely mirrors transient experimental setups but struggles with accurately defining temperatures due to rapid cooling. Conversely, the finite internal resistance model captures internal temperature gradients but risks inaccuracies when nanoparticles are smaller than phonon mean free paths. In the lumped capacitance model, the NP is selectively heated at a temperature Ti and then allowed to cool down in a large constant temperature fluid reservoir at T∞. The NP’s temperature decay is fitted using eqn (3): t

where TNP(t ) is the NP’s temperature at time t, and τ is the thermal relaxation time. The TBC is then calculated as: G¼ CNP Aτ ð4Þ

where CNP is the NP’s heat capacity and A is its surface area. In contrast, the finite internal resistance model accounts for spatial-temporal variations of temperature inside an NP, while assuming that the macroscopic heat conduction equation (eqn (5)) governs this problem: @TNP ðr;tÞ ¼ α∇ 2 TNP ðr;tÞ @t

ð5Þ where TNP(r,t ) is the local temperature at radius r and time t, and α is the thermal diffusivity of the NP. By tracking the temperature of discrete spherical shells over time and fitting to analytical or numerical solutions, the TBC is calculated from the NP’s surface boundary condition, which is a convection or Robin boundary condition. Both methods enable quantitative evaluation of the TBC, but the choice between them depends on the particle’s Biot number, where the TBC is not known a priori. The differences between the two TNEMD approaches, manifested as less accurate temperature relaxation

This journal is © The Royal Society of Chemistry 2025 View Article Online

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Nanoscale Review

profiles and subsequently less precise estimates of interfacial thermal conductance,186 for the interface between an alkane nanodroplet and water. The SNEMD method is widely used to determine the TBC by imposing a continuous temperature gradient across an NPfluid system. In this approach, the system is divided into distinct thermal regions: a central NP or solid region is maintained at a higher temperature (heat source), while distant fluid regions are kept at a lower temperature (heat sink), see the inset in Fig. 9(a). The temperature gradient can be created either using thermostats185 or heat input/output regions.117 Once a steady-state temperature profile is established (Fig. 9(b)), the TBC is calculated using eqn (6): G¼

where q is the heat transfer rate across the solid–liquid interface, A is the interfacial area, and ΔT is the temperature drop at the interface, see Fig. 9(a). If the thermostat method is used, q can be computed from the rate of energy added to (or removed from) the thermostatted region as: 1 ΔE q¼ Δt A

ð7Þ where ΔE is the total energy added or extracted over time Δt, see Fig. 9(b). To obtain ΔT, the solvent’s and NP’s temperature profiles are spatially averaged, and a linear fit is performed on the solid regions away from the interface; the extrapolated temperature discontinuity at the interface yields ΔT. A detailed understanding of interfacial thermal coupling can be obtained from EMD methods.75 Unlike non-equilibrium simulations, EMD does not generate temperature gradients; instead, it relies on the fluctuation-dissipation theorem to relate heat transport properties to equilibrium energy exchanges between the solid and liquid regions. Once the system reaches thermal equilibrium at a given temperature T0,

the EMD method uses fluctuations in the interfacial heat power to calculate G via the Green–Kubo formalism:75 G¼ 1 AkB T0 2 ð1 hPðtÞPð0Þidt ð8Þ 0

where A is the interfacial cross-sectional area, kB is the Boltzmann constant, T0 is the equilibrium temperature of the system, P(t) P(0) is the instantaneous heat power/flux across the interface, 〈〉 denotes an ensemble average for the time autocorrelation function of the fluctuating heat power. Barrat and Chiaruttini75 initially proposed to calculate the fluctuating heat power across the solid–liquid interface according to eqn (9): X X PðtÞ ¼ Fij  νi ð9Þ i[liquid j[solid

where Fij is the force vector produced from the interaction between atoms i and j across the interface, and νi is the velocity vector of atom i. Similarly, the heat power across the interface can be calculated as:187 PðtÞ ¼

dEi ðtÞ dt ð10Þ where, Ei (t ) is the internal energy on one side of the interface at time t. Barrat and Chiaruttini75 noted that the conventional Green–Kubo relation shown in eqn (8) is strictly valid only in the thermodynamic limit, where the volumetric heat capacity cv → ∞. Since MD simulations inherently deal with finite systems, they proposed a modified Green–Kubo expression (eqn (11)) for solid–liquid interfaces in finite domains, calculating G from the long-time integral of the interfacial heat power autocorrelation function.75 ats

Fig. 9 SNEMD simulation setup for an AuNP–water system: (a) Temperature distribution from the heated AuNP to the cooled water at the boundaries of the simulation domain. (b) Time-dependent accumulation of energy input into the AuNP and extraction of water; the linear slope indicates the established heat flux.

This journal is © The Royal Society of Chemistry 2025 Nanoscale, 2025, 17, 20803–20830 | 20817 View Article Online Review Nanoscale

Recently, Rajabpour and Volz188 developed a new Green– Kubo expression for the TBC of finite systems, in which G is derived by integrating the time autocorrelation function of the instantaneous temperature difference between the solid and the surrounding liquid:

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. 1 A ¼ G kB T0 2 ð1 hΔTðtÞΔTð0Þidt

ð12Þ 0 where ΔT (t ) is the instantaneous temperature difference between the solid and the first solvation shell of the surrounding liquid. To enhance statistical reliability, the temperature autocorrelation function is averaged over time and multiple configurations, as shown in Fig. 10. While EMD avoids artifacts introduced by artificial temperature gradients, it requires long simulation times and careful statistical averaging to obtain converged results due to the inherently noisy nature of equilibrium fluctuations. Nevertheless, it is particularly suited for systems where applying external gradients would be physically unrealistic or introduce unwanted nonlinearities. The SNEMD method requires large and unphysical temperature gradients to compute a TBC with minimal statistical noise. Similarly, it may introduce artificial artifacts in the simulation due to nonlinear temperature effects from thermostating.185 In contrast, the EMD approach avoids such artifacts by relying on equilibrium temperature fluctuations, making it more suitable for evaluating intrinsic interfacial properties. Rajabpour et al.185 systematically compared four MD-based techniques and reported that while all predicted TBCs within the same order of magnitude, SNEMD and EMD differed by less than 10%. In contrast, TNEMD with a lumped thermal resistance approximation underestimated the TBC by approximately 25% compared to EMD, whereas TNEMD with a finite internal resistance overestimated it by a similar margin. This discrepancy is likely due to the transient nature of TNEMD,

Fig. 10 Equilibrium MD results for an AuNP in water. Left Y-axis: Temperature difference autocorrelation function over time. Right Y-axis: Interfacial thermal resistance (Inverse of TBC). (Reprinted from ref. 117. With the permission of AIP Publishing.).

20818 | Nanoscale, 2025, 17, 20803–20830 where the NP is initially heated and allowed to cool, making it challenging to define a precise interfacial temperature and, therefore, accurately quantify the TBC. More recently, the TBC calculations produced by the three different MD techniques (lumped capacitance TNEMD, SNEMD, and EMD) were reported by Paniagua et al.117 for the TBC across Au–water interfaces. It was reported that although the EMD method provides fewer artificially imposed conditions on the system, the implementation is computationally expensive and prone to instabilities during the computation of the time autocorrelation function. Therefore, due to its smaller temperature gradients and affordability, the SNEMD method was appropriate to compare the TBC of the different interfaces.117 Thus, careful consideration of these factors is crucial for accurate TBC evaluation in molecular-scale thermal transport investigations. 4.2 Computational modeling of solvated functionalized AuNPs Gold is a highly valuable material with diverse applications in fields such as materials science, nanotechnology, and catalysis.189 Consequently, the literature on MD modeling of Aubased systems is extensive. Common FFs for modeling gold include Lennard-Jones (LJ),190 Morse, Embedded Atom Method (EAM), and Quantum Sutton-Chen (QSC). Pairwise FFs, such as LJ and Morse, offer computational efficiency and are parameterized to match key physical properties. LJ potentials effectively capture the lattice constant, cohesive energy, bulk modulus, and free surface energy of metals.190–194 In contrast, Morse potentials, featuring tunable parameters, are commonly used to represent covalent bonds195,196 but can also represent non-covalent interactions, though LJ potentials perform better in this role.197 Morse FFs accurately reproduce the bulk properties of Au,191 and Au surfaces, providing useful modeling of thiol adsorption on Au substrates198 (see Subsection 4.3.2 for more details). While pairwise FFs are widely used in MD simulations, they have notable limitations, especially in modeling metallic solids. Specifically, pairwise FFs cannot adequately represent electron delocalization and charge transfer, essential for accurately modeling metallic bonds.193,199 To address these shortcomings, more sophisticated FFs like EAM, first introduced in the seminal work by Daw and Baskes, have been developed.200,201 For Au, EAM has successfully described the properties of both bulk systems202 and NPs,203–205 but at higher computational costs compared to pairwise FFs. Another sophisticated approach is the QSC FF, which also accurately describes the metallic interactions.206 However, both the EAM and QSC FFs do not account for polarization effects, limiting accuracy in simulating the interfaces between a metal substrate and a polar adsorbate. To address this gap, Bhattarai et al.92,207 developed the density readjusting embedded atom method (DR-EAM), explicitly incorporating polarization effects by assigning partial charges to atoms, thereby adjusting valence electron densities. Interestingly, this contribution demonstrated that polarization effects had minimal influence on the interfacial vibrational characteristics and thermal trans- This journal is © The Royal Society of Chemistry 2025

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Nanoscale port properties at metal–water interfaces. Thus, the choice of FF for modeling Au should be guided by specific interaction properties, system complexity, and computational efficiency needs. Water is widely used in solvated AuNP models due to its similarities to biological fluids208 and well-documented properties. Its models typically use Coulomb’s law for electrostatics and LJ potentials for dispersion and repulsion. Charges may be placed at atomic cores or on dummy sites, with the LJ term frequently reserved for oxygen–oxygen interactions. Water models are distinguished by factors such as interaction point count, rigidity or flexibility, and polarization inclusion. Threesite models, such as SPC (Simple Point Charge) and TIP3P (Transferable Intermolecular Potential 3-Point), are among the earliest water models.209–212 These models apply partial charges and LJ parameters to the oxygen atom and the two hydrogen atoms. While their simplicity allows for efficient large-scale MD computations, they regard the water molecule as rigid and nonpolarizable, limiting their capacity to effectively represent temperature-dependent characteristics. Foursite models, notably TIP4P and its improved versions such as TIP4P/2005, add a dummy atom to refine the electrostatic distribution.209,213 Although they outperformed three-site models in terms of predicting structural properties and phase diagrams, they lacked explicit polarizability, limiting their accuracy in strong electric fields or inhomogeneous situations. To further enhance the accuracy, five-site models like TIP5P and six-site models have been developed,214,215 but they are less frequently employed due to high computational demand. Sirk et al.216 conducted a comprehensive comparison of the thermal conductivities of rigid and flexible water models using MD. Flexible models, such as TIP3P/Fs, SPC/Fw, and SPC/Fd, demonstrated 15–25% greater conductivities than their rigid counterparts due to having more degrees of freedom, allowing for more efficient energy transmission. On the other hand, thermal conductivities for rigid models like SPC, SPC/E, TIP3P-Ew, and TIP4P-Ew ranged between 0.776 to 0.816 W m−1 K−1, with an average of 0.799 W m−1 K−1, while the experimental value is 0.609 W m−1 K−1 at 300 K. Recent advancements include polarization-corrected and flexible water models. Polarization enhancements better represent water properties in polarizable environments.212 Rigid water models use position constraints to treat bonded interactions implicitly, while flexible models capture anharmonic O–H bond stretching. Researchers demonstrated that the choice of water model has a considerable impact on the strength and heat flux dependency of TBC at the nanoscale Au–water interface.217 They found that due to increased phonon coupling at low frequencies, the rigid TIP3P model provides higher and more consistent TBC values. The flexible TIP3P model, on the other hand, yields somewhat lower and temperature-dependent TBC, with increased conductance at greater heat fluxes. Given the impact of the water model choice on TBC calculations, careful model selection is essential in solvated AuNP simulations. As shown by Sirk et al.216 and Munjiza et al.,217 flexible water models (e.g., TIP3P/Fs, SPC/Fw) consistently

This journal is © The Royal Society of Chemistry 2025

Review yield higher bulk thermal conductivities (up to 25% more than their rigid counterparts) due to additional vibrational degrees of freedom. These models also exhibit stronger TBC-heat flux dependence, capturing non-linear effects that may be relevant in high-temperature applications. Alternatively, rigid models such as SPC/E or TIP3P-Ew are computationally efficient and often produce more stable TBC values, particularly in nearequilibrium conditions. Munjiza et al.217 found that the TBC between AuNPs and water remained nearly constant for the rigid TIP3P model, while the flexible model produced increasing TBCs with higher heat fluxes. Consequently, researchers could prioritize rigid water models for low-flux or screening studies, where computational cost is a constraint, and consider flexible models when studying flux-dependent interfacial behavior. Similarly, authors should consider adding to their study water models widely used in the field for benchmarking purposes. Ultimately, model selection should reflect the study’s objectives, the expected thermal regime, and the level of accuracy required for interfacial water structuring and dynamics. Thiolates (R–S–) are sulfur-based radicals bonded to organic groups, commonly found in thiols (R-SH) and as ligands in Au–SAM interfaces or thiolate-protected AuNPs. Their interactions in functionalized metallic nanostructures are often modeled using all-atom (AA) and united-atom (UA) FFs. Examples include the OPLS-AA and OPLS-UA FFs,138,157,218,219 the TraPPE-UA FF,158,220 and FFs developed for organic molecules and proteins, such as CHARMM,221 AMBER,222 and GROMOS.223 All-atom FFs explicitly account for interactions between each atom intramolecularly, while united-atom FFs group certain atoms, such as hydrogen bonded to carbon atoms, into a single interaction center. Both FF types incorporate bonded (covalent) and non-bonded (van der Waals and electrostatic) interactions. Bonded interactions include terms for bond stretching, angle bending, dihedral angle torsion, and improper torsion, providing a detailed framework for modeling thiolates. In practice, these different force fields are often combined to model multi-component systems. For instance, in a comparative study of silica–water and gold–water interfaces, the CHARMM potential was used for the hydroxylated silica surface, a Morse potential for gold, and the TIP3P model for water, all within a single simulation framework to elucidate the mechanisms of heat transfer.73 Unlike traditional FFs, which need individually parameterized models for gold, water, and thiol ligands, and frequently rely on empirical combination rules to address cross-interactions, ReaxFF offers a unified reactive framework capable of expressing any pertinent interactions under a single parameter set.224 ReaxFF can dynamically describe bond formation and breaking, allowing for reliable modeling of Au–S chemisorption, interfacial water structuring, and ligand reconfiguration under temperature gradients.224,225 This avoids the need to mix diverse non-reactive force fields (e.g., EAM for Au, SPC for water, and OPLS-AA for ligands), minimizing compatibility issues and enabling more adaptable and chemically consistent simulations under different circumstances.

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Review Nanoscale

While, in general, ReaxFF is more computationally demanding than classical, non-reactive force fields, this cost is a necessary trade-off for its ability to model complex chemical reactions. Simulations of hydrocarbon systems show that ReaxFF can be approximately 50 times slower per time-step than a classical MD code like GROMACS that uses a simplified united-atom model.226 This increased cost stems directly from the method’s complexity, which includes the dynamic calculation of bond orders and geometry-dependent atomic charges at each MD step to accurately simulate bond formation and breaking.224,227 ReaxFF’s cost can be justified by significant gains in accuracy and predictive power in systems where chemistry is paramount. Studies demonstrate that ReaxFF excels at predicting experimental data and quantum mechanics (QM)-level calculations for a wide range of properties. ReaxFF accurately reproduces experimental heats of formation for dozens of hydrocarbons with an average deviation often less than 4 kcal mol−1, a level of accuracy comparable to or better than the semiempirical PM3 method.228 Furthermore, it reliably predicts molecular geometries, such as bond lengths and angles, that are in excellent agreement with experimental data and ab initio QM results.226,228 Crucially, ReaxFF has proven highly effective at modeling the energetics of chemical reactions. It accurately describes the potential energy surfaces for various bond dissociation events, showing strong agreement with high-level QM calculations.227,228 This capability extends to complex processes like phenolic pyrolysis, where calculated bond energies and reaction barriers show a reasonable match with DFT and high-level CCSD(T) results.72,229 This predictive accuracy is achieved at a fraction of the computational cost of QM methods. Performance comparisons show that for a system of approximately 450 atoms, ReaxFF is a million times faster per iteration than DFT.228 This enormous performance advantage is amplified in larger systems, as ReaxFF scales nearly linearly with the number of atoms, Table 2

whereas the computational cost of QM methods scales much more poorly, typically ranging from O(N3) to O(N7).226 This favorable scaling allows for reactive simulations of thousands of atoms over nanosecond timescales—a regime entirely inaccessible to more fundamental methods.226,228 Therefore, ReaxFF occupies a vital position, providing the accuracy needed to reliably model reaction chemistry while retaining the computational efficiency required to study the large-scale, long-timescale dynamics relevant to AuNP–ligand interactions and drug release. Table 2 summarizes the FFs for solvated functionalized AuNPs modeling. 4.3

Modeling interfaces in MD is inherently more complex than modeling bulk materials or isolated molecules. Interfaces involve the coexistence of two dissimilar materials, requiring detailed analysis of molecular adsorption onto the interface and diffusion across the interface. This is particularly true when modeling interfaces with engineered geometries, such as the nanostructured or finned surfaces used to enhance heat transfer. These systems require careful analysis of phenomena like fluid adsorption within nano-grooves and the formation of a “solid-like liquid layer,” which are critical for accurately predicting thermal transport.72 Interfacial FFs are typically optimized to capture surface properties such as wetting behavior, surface tension, and adsorption energies.194,230 However, FF parameters tailored for interface interactions are not always readily available. In such cases, simple approaches, such as empirical combination rules for LJ FF parameters, are often employed. Nonetheless, for interfaces exhibiting both physisorption and chemisorption, non-bonded interactions alone lack the precision required to accurately describe interfacial structures. For strongly interacting interfaces, such as Au–S systems, the FF must be developed with a focus on the underlying adsorption physics to ensure accuracy. The literature on ligand–solvent interfaces is relatively limited. Consequently, it is common practice to model ligand–

Summary of FFs for solvated functionalized AuNPs models Force field Principal characteristics Pros/Cons Au modeling Pairwise FFs (LJ, Morse) Simple two-body potentials

• Computationally efficient • Poor at capturing metallic bonding • Better accuracy for metals • Higher computational cost Many-body FFs (EAM, QSC, DR-EAM) Include many-body and/or polarization effects Water modeling Rigid models (SPC, TIP3P, TIP4P/2005)

Fixed bond lengths and angles Flexible models (TIP3P-Fs, SPC/Fw) Allow bond vibrations

Ligands modeling All-atom FFs (OPLS-AA, CHARMM, AMBER) United-atom FFs (OPLS-UA, TraPPE-UA) Unified Option: ReaxFF

Explicitly model each atom; detailed bonded and non-bonded interactions Groups non-polar H atoms with heavy atoms for efficiency Reactive FF handling bond formation/ dissociation across Au, water, and ligands

20820 | Nanoscale, 2025, 17, 20803–20830 • Efficient for large simulations • Limited accuracy for temperature-dependent properties • Captures thermal effects better • Costlier than rigid models • High accuracy for organics • Tedious parameterization and high cost • Lower computational demand • Less detail for H bonding • Chemically consistent unified modeling • Relatively computationally expensive

This journal is © The Royal Society of Chemistry 2025 View Article Online

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Nanoscale solvent non-bonded interactions using Lorentz–Berthelot combining rules, which merge parameters from all-atom FFs (see Section 4.2) with the non-bonded interaction parameters of water models.158,189,220,231 For instance, in systems of functionalized metallic nanoparticles immersed in water, several works have employed the SPC/E water model combined with the OPLS FF.157,220,232 The body of work on MD models for Au– water interfaces is extensive, and concise details are provided in the following Subsection. 4.3.1 Metal–solvent interfaces. Most MD models of metallic–fluid interfaces define the interfacial interactions between the metal and the adjacent fluid using a truncated 12–6 LJ FF, with Lorentz–Berthelot combining rules commonly applied to determine the FF parameters. Berg et al.233 examined the limitations of using simple LJ FFs and combination rules for describing interfacial interactions in Au–water systems. Their findings revealed that LJ FFs poorly replicate adsorption energy curves obtained through density functional theory (DFT) simulations, particularly when compared to more sophisticated pairwise FFs, such as Buckingham or Morse. Furthermore, Berg et al.233 were unable to identify a suitable set of FF parameters that could successfully use combination rules to match DFT-calculated adsorption energies, highlighting the need for more accurate approaches to model such interfaces. Alternatively, FF parameters for metal–fluid interfaces have been determined by optimizing interfacial properties of interest, such as experimental contact angles,135 adsorption energies,230 surface tension,194 or DFT-derived adsorption energy curves.233,234 However, FF parameters optimized using different methods often exhibit discrepancies when calculating interfacial properties.111,122 Additionally, these parameters are rarely transferable across different systems. It is then imperative to shift the modeling strategy from generic empirical FFs to new FFs based on the physics and chemistry of specific interfaces. ReaxFF,224 a bond order-based force field capable of

Review accounting for reaction energy barriers, chemical bonding, and non-bonded interactions of solid-water interfaces, is needed to model Au–water interactions for a better understanding of thermal transport. There is substantial research focusing on MD models for thiolate adsorption on Au surfaces and mathematical representations of the Au–sulfur bond. These models are summarized in the following Subsection. 4.3.2 Gold–SAM interfaces. Despite extensive research on SAM–Au interfaces and their numerous applications, the precise arrangement and configuration of adsorbed thiolates have been a topic of debate for decades.235–241 Notably, most investigations on SAM–Au interfaces focus on systems with an exposed Au (111) surface.242 Early electron diffraction characterizations suggested that thiols organize into a hexagonal (√3 × √3) R30° lattice,243,244 commensurate with the underlying Au (111) surface, as depicted in Fig. 11(a). Theoretical and computational investigations of this (√3 × √3) R30° model propose that thiolates preferentially occupy threefold-coordinated hollow sites, twofold-coordinated bridging sites, or positions directly above Au surface atoms,236–238,245–250 as illustrated in Fig. 11(b). The traditional (√3 × √3) R30° lattice arrangement has been increasingly challenged over the years. X-ray measurements revealed the existence of a centered (4 × 2) superlattice derived from the (√3 × √3) R30° structure. This centered (4 × 2) superlattice consists of four atoms, where two adsorbed thiols are equivalent, and the other two occupy distinct lateral and vertical positions relative to the underlying Au (111) surface251 (see Fig. 11(a)). Furthermore, recent scanning tunneling microscopy (STM) visualizations have identified the presence of Au adatoms at the SAM–Au interface.235,252–254 These observations demonstrated that thiolates form RS–Au– SR complexes on a reconstructed Au (111) surface, particularly at low thiol coverage.235,253,255,256 STM has revealed that the Au–SAM interface is better described as a complex assembly of thiolates bonded to Au

Fig. 11 Adsorption of sulfur on the Au (111) surface. (a) Early experiments identified that adsorbed sulfur atoms organized in a hexagonal (√3 × √3) R30° lattice above the Au (111) surface (solid black line), and the existence of a c (4 × 2) superlattice (dashed blue line). Gold, teal, and white spheres represent Au, carbon, and hydrogen atoms, respectively. (b) Theoretical and computational calculations concluded three possible adsorption sites for the sulfur atoms: hollow sites, bridge sites, and atop sites. Gold, orange, and ochre spheres represent Au atoms in the outermost, 2nd and 3rd layer of an Au (111) surface, respectively, while yellow spheres represent adsorbed sulfur atoms.

This journal is © The Royal Society of Chemistry 2025 Nanoscale, 2025, 17, 20803–20830 | 20821 View Article Online

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Review adatoms, rather than thiolates directly adsorbed onto an atomically flat Au (111) surface.252 Additionally, the Au–SAM interface has been shown to exhibit dynamic behavior, including the diffusion of thiolates on the Au surface and the exchange of adsorption sites.257,258 For instance, DFT combined with ab initio MD simulations has demonstrated that adatom-based structures can emerge from the reconstruction of the interface when a system initially organized under the (√3 × √3) R30° lattice model is allowed to relax.255 However, STM measurements confirming the presence of RS–Au–SR complexes have not yet been achieved for intermediate to full thiolate coverage on Au–SAM interfaces.259,260 The presence of adatoms is particularly prominent in AuNPs, which exhibit a core–shell-like structure, where the gold core atoms are surrounded by shell-like Au adatoms bonded to the sulfur head groups of thiolates.261 The formation of RS–Au–SR complexes was first observed in AuNPs by Jadzinsky et al.,262 who identified dimeric and monomeric RS– Au–SR staples (see Fig. 12), with their distribution and ratios varying depending on the nanoparticle size. More recent research has suggested the existence of trimeric SR(–Au–SR) staples,263 bridging thiolates,264,265 and cyclic –Au–SR structures.263 However, distinguishing the latter from the more abundant dimeric and monomeric RS–Au–SR staples remains challenging.266 DFT has been widely employed to investigate the structural properties of thiolated Au surfaces. However, the omission of dispersion forces in some DFT simulations has introduced uncertainties in the adsorption behavior of thiols on Au surfaces,268 with results showing a strong dependence on the chosen functional.237 Several investigations have focused on elucidating the fundamental aspects of AuNP-ligand bonding. For instance, Reimers et al.269 used DFT calculations to show that the stability of sulfur-stabilized AuNPs arises from local gold–sulfur covalent interactions rather than from complete electron shell closure. Complementing this, Tang and Jiang270 systematically evaluated the binding strengths of various ligands on gold surfaces, revealing clear trends where bulky N-heterocyclic carbenes and alkynyl groups form particularly robust bonds. In a similar vein, Pensa et al.260 reviewed the complexity of the sulfur–gold interface, highlighting the multiple coordination modes and dynamic surface reconstructions that challenge simplified models of AuNP functionalization. Methodological advances have also played a crucial role in deepening our understanding of these systems. Fusaro et al.271

Nanoscale combined DFT with solvent models to accurately predict adsorption energetics and ligand exchange processes, while Berg et al.233 optimized force fields for water–gold interactions, improving the reliability of MD simulations in reproducing experimental observations. Several reviews272–274 have synthesized these computational approaches, addressing challenges such as weak intermolecular forces, multiscale phenomena, and the dynamic nature of the bio interface. Nevertheless, the high computational cost of DFT makes it impractical to model the length and time scales found in heat transfer processes in complex Au–SAM interfaces. Consequently, MD simulations have become the preferred approach for investigating thiol-protected AuNPs275–279 and flat Au–SAM interfaces.280,281 Unfortunately, FF parameters available for describing the Au–S bond are often developed focusing on specific adsorption models. As a result, these FFs are frequently non-transferable between Au–SAM and thiolated AuNP systems, or even between AuNPs of different sizes.282 However, reactive FFs, such as ReaxFF, may be trained using the high-fidelity insights from DFT calculations, which include surface reconstructions, bonding configurations, and adsorption energies to make them more accurate and transferable across different Au–S interfaces. An early attempt to model the Au–sulfur interaction in SAMs utilized an LJ FF.283,284 However, this approach had notable limitations, including the oversimplification of the Au surface as flat, as the Au–S potential energy was treated as dependent solely on the perpendicular distance from the surface. Furthermore, LJ FFs have been shown to inadequately capture the chemisorption behavior of sulfur on Au surfaces.197 To address these limitations, Perstin and Grunze285 developed a modified LJ FF that incorporated a surface corrugation function and explicitly accounted for Au–S–C angle bending in thiolates, providing a more accurate representation of the Au–SAM interface. Morse FF parameters have been developed as an alternative to better describe the chemisorption of sulfur on Au surfaces.286–288 However, these parameters are often tailored to specific sulfur binding sites, such as the threefold hollow site on Au (111) surfaces,286,287 limiting their transferability to other Au surface configurations or sulfur binding models. To address this limitation, more sophisticated functional forms and models have been proposed to accommodate different binding site scenarios. For instance, Longo et al.289 developed

Fig. 12 Monomeric and dimeric RS–Au–SR staples in the shell of AuNPs.267 The green and yellow spheres represent gold and sulfur atoms, respectively. (a) Rectangular staple formation present in Au38, Au102, and Au144 clusters. The V-shape staple formation in (b) appears in Au25 and Au38 clusters and in (c) in Au102 clusters (together with (a)). (Reprinted (adapted) with permission from ref. 267. Copyright 2016 American Chemical Society.).

20822 | Nanoscale, 2025, 17, 20803–20830 This journal is © The Royal Society of Chemistry 2025 View Article Online Nanoscale Table 3 Summary of key insights from DFT and experimental studies on thiolated Au interfaces

Aspect Key findings Method/reference SAM lattice structure

Hexagonal (√3 × √3)R30° and c(4 × 2) superlattices observed depending on coverage and reconstruction243,244 Occupation of hollow, bridge, and atop sites;236–238,245–250 sensitive to the local environment and surface relaxation RS–Au–SR staples (monomeric, dimeric) were observed especially at low coverage and in AuNPs267 Covalent Au–S bonding dominates; dispersion effects are important; binding strength varies by ligand type269 Thiol adsorption induces surface rearrangements and dynamic site switching260 Fixed-site or LJ-based FFs fail to reproduce adsorption flexibility and reconstruction Use of Morse,286–288 Gupta,289 GolP,290 ReaxFF228,291–293 FFs, and DFT-calibrated approaches for specific binding modes

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. Review

Ligand binding stability Surface reconstruction Limitations of traditional FFs Advanced modeling strategies

a modified Gupta FF to account for Au vacancies and adatoms at the Au–sulfur interface, while the GolP FF290 was designed to accurately represent the atop binding site for sulfur. These advancements offer improved flexibility in modeling complex Au–sulfur interfaces. The non-transferability of Au–SAM FFs poses a significant challenge in modeling thiolate-protected AuNPs. Due to the highly anisotropic surface of AuNPs, FFs developed for specific binding sites on flat Au surfaces cannot be directly applied. To address this, robust elastic network models with optimized force constants have been developed to describe the AuNP core–shell structure (including adatoms) and Au–sulfur interactions.157,218 In particular, ReaxFF228,291–293 is well-suited to capture complicated chemisorption and surface reconstruction events in thiolated AuNPs due to its capability to dynamically represent bond dissociation and formation. Similarly, Pohjolainen et al.282 developed a transferable all-atom FF for thiolated AuNPs, with parameters optimized to account for various staple units, enabling more versatile modeling of complex Au–sulfur interfaces. Table 3 summarizes the key insights from DFT and experimental studies on Au–S interfaces.

Although substantial progress has been made in computationally and experimentally elucidating interfacial heat transport in solvated AuNPs, several critical challenges remain unsolved. Highlighting these limitations is essential to contextualize existing findings and to inform more reliable design strategies for spatiotemporal control of solvated AuNPs. The following are itemized research challenges: • Lack of comprehensive insight into solid–liquid thermal transport across pristine interfaces, further challenged by morphology-dependent behavior, and aggregation tendencies in solvated nanoparticles. • Complex interfacial structures hinder accurate analysis of heat transfer across bare solid–liquid interfaces. • Complexity in analyzing heat transport across functionalized solid–liquid interfaces due to the three-component solid- This journal is © The Royal Society of Chemistry 2025

DFT, STM STM, XPS, DFT DFT studies with functionaldependent results Ab initio MD, STM MD modeling comparisons Computational FF development studies

ligand–solvent structure, ligand localized thermal motion due to vibrational and conformational fluctuations, water penetration into the ligand layer, distinct temperature profiles depending on ligand hydrophobicity, and contrasting interfacial solvent mobility effects. • Experimental quantification of solid–liquid thermal conductance is limited by the low sensitivity of traditional thermoreflectance techniques, primarily due to the dominant thermal resistance of liquids. • Equilibrium MD (EMD) avoids artifacts from artificial temperature gradients but requires long simulations and extensive averaging due to noisy equilibrium fluctuations. • Transient heating in TNEMD makes interfacial temperature definition difficult, limiting accurate TBC quantification. • EAM and QSC force fields neglect polarization effects, reducing accuracy for metal–polar adsorbate interface simulations. • Rigid and nonpolarizable three-site water models (e.g., SPC, TIP3P) fail to capture temperature-dependent properties despite being computationally efficient. While four-site water models (e.g., TIP4P) improve structural predictions, but lack explicit polarizability, reducing accuracy in strong electric fields or inhomogeneous systems. • Non-transferability of Au–SAM force fields limits accurate modeling of thiolate-protected AuNPs due to anisotropic nanoparticle surfaces. • Reactive force fields able to couple chemistry and heat transfer are mostly not optimized for both phenomena. Additionally, these force field simulations are computationally expensive, but faster than QM-level calculations. 5.2

AuNP-enabled therapy has the potential to advance the medical field, particularly in applications, such as drug/gene delivery and photothermal therapy. This review highlights their unique properties, including plasmonic behavior, biocompatibility, and versatile surface functionalization, which collectively and combined with precise thermal modulation, could deliver highly-localized drug release and tissue ablation therapies. The consideration of rDA reactions in drug delivery systems further exemplifies how the thermal properties of

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Review AuNPs can be leveraged for spatiotemporal therapeutic interventions. A key focus has been the interfacial thermal transport mechanisms of AuNPs, particularly at gold–water interfaces, where functionalized ligands significantly influence heat dissipation. Computational works have provided critical insights into the interplay between ligand chemistry, interfacial water structuring, and nanoscale thermal transport, shedding light on the modulation of the TBC. Despite these advancements, challenges persist in optimizing solvated functionalized AuNP systems with adequate spatiotemporal temperature control in complex biological environments. Such optimization requires a proper description of interfacial heat transfer across functionalized Au–water systems, which is governed by a complex interplay of multiple mechanisms. In the experimental field, TDTR measurements have dominated in obtaining measurements of TBC across solid–liquid interfaces, though they have critical limitations in terms of sensitivity. Improving or developing new experimental techniques with higher sensitivity to solid–liquid TBC is essential to enhance the accuracy of measurements. For example, combining TDTR with alternative methods like picosecond acoustics and ablation threshold measurements can provide a more comprehensive understanding of TBC and reduce uncertainties. Given the current experimental limitations, it is not surprising that computational modeling remains the favored approach for understanding heat transfer across functionalized Au–water interfaces. In the computational realm, MD models have strongly dominated over ab initio efforts due to the high computational cost of the latter, which limits their applications to the large size and time scales required to model complex functionalized Au interfaces. Numerous MD works have extensively investigated the individual mechanisms governing interfacial thermal transport in functionalized Au–solvent interfaces, including the atomic interaction strength, the role of hydrogen bonding, vibrational mismatch, and the mobility and molecular organization of the liquid phase near the interface. However, the isolated descriptions fail to fully capture the coupled nature of interfacial thermal transport. Additionally, notable discrepancies exist in the literature on the role of mobility and interfacial structuring in defining the TBC of functionalized interfaces. These discrepancies arise from the use of different frameworks to characterize the liquid layering effect, overlooking the broader molecular organization of the interfacial liquid. Notably, for bare Au–solvent interfaces, the role of interfacial structuring is better understood, as TBC calculations have been explained using the liquid depletion layer parameter δ. To further explore interfacial thermal transport, a continuum model incorporating AuNPs and their surrounding medium can be developed. This model must integrate molecular level granularity as temperature-dependent chemistry, and interfacial liquid property variations into a temperaturedependent TBC calculation. By capturing the transient thermal response of AuNPs and the surrounding water, continuum models would ensure interfacial temperature continuity, 20824 | Nanoscale, 2025, 17, 20803–20830

Nanoscale providing a comprehensive framework for nanoscale heat transfer analysis. Moreover, a quantitative analysis of energy exchange contributions at each sub-interface could provide a more precise characterization of interfacial heat transfer in functionalized Au–solvent interfaces. Developing comprehensive transient thermo–chemical models will lay the foundation for investigating complex biomedical applications, such as drug delivery systems using rDA reactions. Nevertheless, as demonstrated in MD simulations, accurately describing interfacial atomic interactions is crucial for predicting interfacial thermal transport. Precise parameterization of FF parameters is essential for reliably reproducing interfacial properties. Reactive force fields, such as ReaxFF, offer a realistic representation of metal–liquid interfaces by capturing bond formation and dissociation; however, their implementation involves high computational costs that must be carefully managed. In this context, ab initio models provide a cornerstone for deciphering the complex surface chemistry of functionalized Au interfaces. Additionally, the advent of machine learning FFs,294 trained to provide the accuracy of ab initio methods in an efficient classical framework, has opened new avenues for the modeling of interfaces.295 Addressing these challenges will unlock the full potential of AuNPs, establishing them as indispensable tools for advancing biomedical technologies. These efforts will enhance the precision and effectiveness of therapeutic interventions while contributing to the broader vision of personalized medicine, where treatments can be tailored for maximum efficacy and minimal side effects. As a final remark, it is important to note that several points discussed throughout this review are not exclusive to functionalized Au–solvent systems and can be extended to other solid–liquid interfaces of interest with appropriate considerations. Nevertheless, for the sake of brevity and coherence, the discussion primarily addressed functionalized Au–water interfaces, as they are the main system of interest in biomedical applications of plasmonic NPs and the central focus of this review.

Author contributions Md Adnan Mahathir Munshi: data curation (experimental and numerical part, equal), formal analysis (experimental and numerical part, equal), investigation (experimental and numerical part, equal), writing – original draft (experimental and numerical part, equal). Emdadul Haque Chowdhury: data curation (experimental and numerical part, equal), formal analysis (numerical part, equal), investigation (numerical part, equal), writing – original draft (numerical part, equal), writing – review & editing (equal). Luis E. Paniagua-Guerra: data curation (experimental and numerical part, equal), formal analysis (experimental and numerical part, equal), investigation (experimental and numerical part, equal), writing – original draft (experimental and numerical part, equal), writing – review & editing (equal). Jaymes Dionne: data curation (experimental part, equal), formal analysis (experimental part, This journal is © The Royal Society of Chemistry 2025

Open Access Article. Published on 21 August 2025. Downloaded on 5/31/2026 1:29:35 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence.

Nanoscale equal), investigation (experimental part, equal), writing – original draft (experimental part, equal). Ashutosh Giri: data curation (experimental part, equal), formal analysis (experimental part, equal), investigation (experimental part, equal), writing – original draft (experimental part, equal), writing – review & editing (equal), funding acquisition (equal), supervision (equal). Bladimir Ramos-Alvarado: data curation (experimental and numerical part, equal), formal analysis (experimental and numerical part, equal), investigation (experimental and numerical part, equal), writing – original draft (experimental and numerical part, equal), writing – review & editing (equal), funding acquisition (equal), supervision (equal).

Conflicts of interest There are no conflicts to declare.

Data availability No primary research results, software, or code have been included, and no new data were generated or analyzed as part of this review.