Heterogeneous thermal conductance of nanoparticle–fluid interfaces: An atomistic nodal approach

Nanoparticles (NPs) have garnered increasing attention due to their wide range of applications in medicine, materials science, and catalysis,^1–5 stemming from the significant dependence of their physicochemical properties with size and shape. In addition, metallic nanoparticles have attracted considerable interest due to their efficient conversion of light into heat, which provides unique opportunities to generate and control temperature fields at the nanoscale.⁶ 

Nanoparticle heating is important in material characterization, for instance, in diffraction studies, where the intense illumination using x-ray probes can trigger nanoparticle melting in ps to ns timescales.^7–9 Atomistic modeling has proven particularly helpful to address the limitations of continuum theory to describe thermal transport processes involving a phase change.^10,11 The unique capabilities of nanomaterials for light/heat conversion have motivated the development of a new research field, thermoplasmonics,¹² which is impacting several areas: hot-electron chemistry, light harvesting, microfluidics, and photothermal therapies. In particular, the heat sources used in photothermal therapy treatments are nanoparticles that emit heat upon excitation by electromagnetic radiation due to local plasmon resonance.¹ Gold nanoparticles have attracted significant interest in biomedical applications due to their low toxicity and wide range of absorption wavelengths, targeting the so-called biological windows.^4,13,14 Therefore, predicting the heat fluxes and temperature fields emerging from these hot particles is an important challenge. Addressing this challenge will allow significant advances in many areas from medicine to material processing.

The interfacial thermal conductance (ITC) or its inverse, the Kapitza resistance, is crucial to quantify and explain heat transport across NP-fluid interfaces. The Kapitza thermal resistance arises from the different thermal transport properties (acoustic contrast) of two phases in contact.^15,16 The most vivid illustration of the Kapitza resistance is a significant temperature drop at subnanometer length scales, namely, at the interface between two materials or materials and fluids. Unfortunately, the ITC is often neglected in studies based on continuum theory since its experimental measurement is still scarce for nanoparticles and prone to significant uncertainties.¹⁷ 

Computer simulations provide a robust approach to compute the ITC by linking thermal transport to interatomic forces and relevant properties of solids and fluid phases, such as the thermal conductivity or the vibrational density of states. Molecular Dynamics (MD) simulation methods have evolved in different directions. Equilibrium methods (EMD) rely on the computation of cross-correlation functions, while Non-Equilibrium Molecular Dynamics (NEMD) mimics the experimental setup by explicitly simulating temperature gradients and heat fluxes across interfaces. Within NEMD, transient non-equilibrium MD (TNEMD) relies on fitting the time-dependent temperature relaxation to the solution of the heat diffusion equation. The fitting to a first-order decay law (the “lumped” method) gives direct access to the ITC. However, this approach neglects the thermal conductivity of the nanoparticle, a variable that is important at the nanoscale.

Models that include both the thermal conductivity and thermal conductance were proposed and applied to analyze the transient data of nanoparticles and biomolecules.^18–21 One caveat of the TNEMD methods is that the heat diffusion equations are geometry dependent, making it challenging to analyze the heat transport of nanoparticles with contorted geometries. This issue prevents the identification of regions in the nanoparticle that might feature very different ITCs. Specifically, current TNEMD methods do not account for local thermal conductance in nanoparticles with complex geometries, such as facets, edges and vertices, and more generally, rough surfaces.

Stationary non-equilibrium MD (SNEMD) simulations provide another route to obtain the ITC, G = J_q/ΔT, by computing explicitly the interfacial temperature jump, ΔT, and interfacial heat fluxes, J_q. SNEMD has been applied to solid–liquid,^22,23 liquid–vapour,^24,25 and nanoparticle–fluid interfaces.^17,26,27 Recently, the interplay of nanoparticle curvature and nanoparticle wettability was investigated,¹⁷ uncovering a correlation between high nanoscale curvature and high thermal conductance. This observation supports earlier studies of fluid–fluid¹⁹ and nanoparticle–fluid interfaces²⁸ as well as theoretical models.²⁹ Recent SNEMD computations of gold nanoparticles in hexane reported a dependence of the ITC on particle size also. However, the ITC increased with decreasing curvature.²⁷ That work also reported the thermal conductance of faceted nanoparticles (cuboctahedral and icosahedral shapes) and proposed a correlation between the ITC and the local atomic environments in each nanoparticle via an atomic coordination number.

Previous SNEMD and TNEMD studies on NP thermal conductance relied on “bulk binning” methods. These methods partition the simulation box into bins with simple, well-defined geometries, e.g., slabs or spheres, which are used to obtain density, temperature, and heat flux profiles. This approach works well when considering structures with well-defined geometries, such as spheres or planes. However, the bulk binning method averages out interfacial structural features that deviate from these geometries. By averaging the physical properties in each bin, the bulk binning approach (BBA) assumes a continuous atomic distribution within a sub-volume (bin) and isotropy in the interfacial properties of the regions sampled in a particular sub-volume. Real nanostructures feature roughness, facets, and edges, where the heat flux and the temperature change rapidly depending on the specific location, hence introducing additional degrees of freedom that are difficult to model, e.g., at intersections between facets. Therefore, the heterogeneous structure of nanoparticle surfaces can potentially introduce markedly different local ITCs at subnanometer length scales. This situation is difficult to map into a simple geometry; hence, resolving the ITC spatially requires a higher level of computation to account for the heterogeneity of the interfaces.

We present in this work an algorithmic workflow, the Atomistic Nodal Approach (ANA), to quantify local thermal properties, namely the heat flux, temperature “jumps,” and the ITC at angstrom level resolution. Our approach becomes particularly advantageous to investigate small NPs (few nm diameter) and anisotropic nanoparticles featuring under-coordinated atoms. Moreover, the ANA can be generalized to study a wide range of complex nanoparticle geometries without introducing specially designed binning methods to capture the nanoparticle shapes. We demonstrate the ANA method by investigating a wide range of nanoparticle structures. The ANA approach uncovers the heterogeneous character of the ITC in nanoparticle–fluid interfaces, providing key information to explain the dependence of the ITC of the nanoparticles in terms of the local coordination number of the surface atoms in contact with the surrounding fluid. This makes it possible to identify high and low ITC regions in a nanoparticle. Finally, we discuss the general dependence of the ITC on nanoparticle size and geometry.

We constructed nanoparticles of different geometries. The atoms in the nanoparticle and in the fluid interact through a Lennard-Jones (LJ) potential,

where U is the pairwise interatomic interaction potential, σ is the atomic diameter, and ɛ is the interaction strength. Hereafter, we employ reduced units (denoted with an asterisk “*”) to present our results. We define the interaction parameters as $σNP=σf=1.0,εf=1,εNP=10$⁠, where NP and f refer to the nanoparticle and fluid, respectively. The cross interaction strength, ɛ_NP−f, was obtained using

where F is the heteroatomic interaction parameter set to 0.25. This value of F = 0.25 corresponds to a full wetting state, i.e., the NP-fluid contact angle is 0°.¹⁷ The atomic interactions were truncated at r* = 2.5, and the linear and angular momentum of the nanoparticle were kept to zero at every step to maintain the position and orientation of the NP throughout the simulation. Previous studies of spherical nanoparticles employed the same LJ model,¹⁷ and the wettability of the nanoparticle was quantified as a function of ɛ_NP−f. Moreover, the LJ model has been parameterized to investigate fcc cubic metals, showing excellent performance in terms of densities, surface tensions, metal–water interfacial properties, and mechanical properties.³⁰ Therefore, this generic model represents a good starting point for our investigation and to test our numerical approach.

Our simulations involved several steps. First, we assigned random velocities to the particles, drawn from the Maxwell–Boltzmann distribution at T* = 0.90. Afterward, we equilibrated the system using the Nosé–Hoover thermostat at T* = 0.90 with a relaxation time of 200 δt* where δt* = 0.0025 is the simulation timestep. Nanoparticle heating can be achieved in an experimental setting using optical or magnetic excitation. Following the equilibration period, we simulated the heating of the nanoparticle using a simple velocity rescaling thermostat applied to the core of the nanoparticle. The thermostatting volume was defined at a radial distance $dRhot*=2$ within the NP surface identified using the α-shape algorithm, while the remaining atoms in the nanoparticle were not thermostatted and therefore moved following newtonian dynamics. Hence, the Density of States of the interfacial atoms (nanoparticle and fluid) is not directly affected by the thermostatting process. We find that the interfacial thermal conductance, heat flux, and temperature gradients are independent within the statistical uncertainty of our computations on the thermostatting frequency (see pp. 18–19 in the supplementary material and the result for thermostatting frequencies between 1 and 100 timesteps). We set a cold thermostat at a radial distance $dRcold*=7$ away from the NP surface. The hot and cold regions generate a constant heat flux in the stationary state.

For all the simulations reported in this work, the lattice parameter of the NP prior to equilibration was set to $lc=4/ρNP*1/3$ with the NP density set to $ρNP*=ρNPσ3=1$⁠. The temperature T of each particle of mass m_i and velocity v_i in the simulation was determined using the equipartition equation

where k_B is the Boltzmann constant and N_dof = 3 is the number of translational degrees of freedom for each particle, i.e., the dimensionality of the system. The reduced temperature is defined by T* = k_BT/ɛ_f. To determine the heat flux at each position, we used a local version of the formula derived by Irving and Kirkwood,³¹ 

where J_q,i is the heat flux corresponding to particle i with associated volume V_i, which we calculate using the Voronoi polyhedra approach. The first term in Eq. (4) represents the kinetic energy flux of particle i, and the second term represents its potential energy flux, where ϕ_i is the potential energy of particle i in the field of all the other particles. The third term, referred to as the collisional contribution, takes into account the flux contribution to the forces between particle i and j separated by a distance r_ij, and it is the dominant contribution in dense fluids and solids.

A similar equation can be written for a control volume (V) used in the “bulk binning” method,

The specific version of the Irving-Kirkwood derivation in Eq. (4) is used to obtain the local heat flux on each atom in the nanoparticle. The location of the NP atoms is well defined throughout the simulation (with only minor vibrations about their equilibrium positions), i.e., a given atom does not diffuse from its initial position, and it occupies the same relative position with respect to the other NP atoms throughout the simulation. This allows us to calculate the per-atom local volume, V_i, which we identify with the Voronoi volume that is calculated for each atom as a function of time. Practically, the per-atom volume, V_i, was obtained using the Voronoi construction of both fluid and nanoparticle atoms in SNEMD simulations spanning 10⁵ steps and ten independent runs [see Fig. 7(b) for a visualization of a Voronoi tessellation].

The simulation workflow was divided into three stages: initialization for 10⁶ steps and NVT equilibration under the Nosé–Hoover thermostat, followed by thermostatting using the velocity rescaling thermostat every timestep at the hot and cold regions. The initial 5 × 10⁴ timesteps were discarded, and the remaining 5 × 10⁴ production steps were used for heat transfer measurements. The resulting trajectories were then used to perform the analyses either with the bulk binning or the ANA approaches. All the trajectories were generated using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) (19 March 2020 version).^32,33

Previous computational studies of the ITC^17,20,27 have predominantly employed the bulk binning approach (BBA), where the heat flux and temperature profiles are obtained from averaging the thermal properties of all atoms inside sampling volumes located between the thermostats and in the direction normal to the interface of interest [see Figs. 1(a) and 1(b)]. We performed BBA calculations for specific particle geometries to provide a benchmark for our implementation of the ANA method. Specifically, we computed density, temperature, and heat flux profiles, using trajectories spanning 5 × 10⁴ steps. We performed five independent runs to obtain the final averages and calculate statistical uncertainties. The ITC was obtained from the analysis of the temperature and heat flux profiles by fitting them to the continuity equations. This is further elaborated in the following (see Subsection II C).

We used the BBA route to quantify the ITC of nanocylinders and nanospheres [see Figs. 5(a) and 5(c)]. The cylindrical nanoparticle spanned the whole simulation box in the axial direction; hence, we considered only the radial heat flux since there is no temperature gradient nor heat flux in the axial direction.

The validity of the continuity equations such as Fourier’s law of thermal conduction was evaluated both numerically and analytically,

Fourier’s law states that the rate of heat transfer q_r over an area A(r) is proportional to the thermal conductivity κ of the material normal to the temperature gradient. Fourier’s law can be verified by sampling the radial heat flux $qrA(r)$ across various radii using the Irving–Kirkwood formula and by examining its relationship with the numerical thermal gradient calculated from simulated data.

For a cylinder of length L, Fourier’s law can be derived analytically to obtain a logarithmic relationship between the radial distance r from a reference radius r_c and the temperature T difference from a reference point T_c. The heat rate experienced by the NP and fluid can be expressed with q_r and the thermal conductivity of the particle (p) or fluid (f), κ,

The heat rate experienced by the NP and fluid can be expressed with q_r and is the same for both NP and fluid due to the conservation of the energy in a steady-state system.

where r_c is the radial distance defining the temperature jump, $Tcp−Tcf$ (see Subsection II C 3).

Following Subsection II C 1, we use Fourier’s law to obtain an equation for the radial heat flux through a spherical surface,

The temperatures for the particle and the fluid are given by

where again r_c is a reference radius corresponding to the interfacial temperature jump. The heat rate in Eqs. (7) and (13) was calculated from the average amount of the energy exchanged at the hot or cold thermostats in time δt, q_r = |ΔU/δt|.

To model the temperature transition between the NP T_p(r) and fluid T_f(r) and account for their different thermal conductivities κ_p and κ_l and the interfacial temperature “jump” we used an error function equation,

where r_c refers to the transition point and ω quantifies the steepness of the error function $erf(π(r−rc)ω)$⁠.³⁴ The temperature “jump” ΔT can be found using ΔT = T_p(r_c) − T_f(r_c), which is defined in terms of the temperature of the particle, $Tcp$⁠, and the fluid, $Tcf$⁠, at r_c.

Unlike the bulk binning approach, the ANA makes no assumptions about the NP geometry. Instead, it projects local temperatures and heat fluxes on each atom (“node”) in the NP. In addition, it uses various computational geometry techniques to identify the atoms defining the interface and the ITC.

A cuboid grid is used to calculate the local properties of the fluid atoms around the NP. The calculations are performed on a per-atom and a per-grid basis for NP and fluid atoms, respectively. We show in Fig. 2 an example of a nanoparticle immersed in a fluid, which illustrates the definitions used for the NP nodes and fluid grids used in this work.

Following the calculation of the averaged per-atom data from the NP coordinates, we employed the α-shape algorithm,^35,36 which was used to construct a triangular surface around the NPs, using a spherical probe with α-radius $Rα*=1.1$⁠. The probe radius is an important variable in our method. Our choice (⁠$Rα*=1.1$⁠) corresponds to the nearest neighbor distance between atoms belonging to the nanoparticle. This probe radius identifies under-coordinated atoms at the nanoparticle surface and prevents the inclusion of atoms that do not belong to the surface in the ITC analysis. One advantage of the α-shape approach is that it quickly identifies the atoms protruding from the surface, allowing easy identification of the solid–liquid interfaces. This method has also been employed successfully to investigate the intrinsic surface of liquid–vapour interfaces.³⁷ 

Subsequently, the temperature gradient $∇T(x,y,z)=∂T∂x,∂T∂y,∂T∂z$ was calculated to determine the direction of heat flow at the NP–fluid interface. ∇T(x, y, z) was computed using trilinear interpolation methods, further elaborated in the supplementary material (see Sec. 3). The direction of the vector ∇T(x, y, z) was observed to be largely parallel to the normal of flat faceted surfaces. We employ in all our computations of the ITC the heat flux component in the direction of the thermal gradient.

To calculate the interfacial temperature “jump,” ΔT, i.e., the NP-fluid temperature discontinuity per atom, we used the temperature of the fluid grids within the first solvation shell of the NP’s surface atoms. To find the nearest grids, we constructed a binary tree for both the fluid and particle atoms. The resulting k-dimensional tree (KD-tree) was built by dividing the space occupied by the fluid grid into cells using the “sliding midpoint” method. Physical properties such as the fluid density and temperature of the nearest fluid neighbors surrounding a NP atom can be efficiently identified using KD-trees in combination with the sliding midpoint rule (see Sec. 9 in the supplementary material for additional information on these methods).³⁸ Each node on the tree represents a partition of the space surrounding a reference point (i.e., a surface NP atom) that allows us to quickly query the fluid (F) nearest neighbors $NNP$ whose radii r_NP − r_F are within a fixed distance M₁ of the atom in the NP, where M₁ is the first minimum in the NP-fluid radial distribution function (see Fig. 3 for an illustration of the nearest fluid-neighbors) ,

After identifying the nearest neighbors for each atom at the surface of the NP, we computed the local temperature, density, and temperature gradient ∇T. We used a density weighted average across all local neighbors to compute the fluid temperature in each grid box in Fig. 2. The density weighted average was used to ensure that each neighboring grid property is determined by the fluid density ρ_F of each grid to avoid over-sampling under-populated regions. The temperature discontinuity ΔT corresponding to each NP atom, i, at the surface is calculated from

where T_NP,i and T_F are the per-atom NP (atomic node, i) and fluid grid temperatures, respectively.

The ITC was obtained by using the heat flux, J_q, projected along the temperature gradient, and the corresponding ΔT [from Eq. (21)] is obtained for each NP surface atom,

We note that the α-shape sorting method discussed above gives access to other important local properties, such as the coordination number (CN) between atoms in the nanoparticle and the number of fluid atoms solvating a given atom i located at the surface of the nanoparticle, which can be extracted from $NNP,i$⁠.

Figure 4 shows the flow diagram employed in the ANA introduced in this work. First, we identify the core atoms and apply the hot thermostat on those atoms up to r* = 2 away from the nanoparticle surface defined by the α-shape algorithm. The calculations are performed using “computational farming,” namely, we run 10³ simultaneous independent simulations, 5 × 10⁴ steps each, and we use these simulations to obtain averages. We then calculate the local heat flux using the IK Eq. (4) and the local temperature gradients using linear interpolation (see Sec. 3 in the supplementary material). The projected heat flux along the direction of thermal gradient and the local ΔT [Eq. (21)] gives the ITC.

Except for long rectangular cuboids used for the octagonal nanowire specified later, all other simulations used cubic grids of l* = 0.4. This size was optimized to obtain the local fluid properties at a good resolution. For each NP atom, a unique identification number was assigned, which allowed the averaging of their thermodynamic properties across different simulations. Information on the particle dimensions and simulation box dimensions are provided in the supplementary material, Sec. 1.

We show in Fig. 5 the temperature and heat flux profiles of cylindrical and spherical nanoparticles obtained with the ANA and BBA approaches. For the BBA calculations, we used cylindrical binning for nanowires and spherical binning for nanospheres. We fitted the radial temperature profile data between the hot and cold thermostatting regions to Eq. (19) and the corresponding temperature equations in Sec. II C 1 (cylinders) and Sec. II C 2 (spheres). The thermal conductivities (κ_p and κ_f) obtained from the fitting of the temperature profiles to Eqs. (17) and (18) are consistent with those reported in Ref. 17 for spherical nanoparticles. Equation (19) was used to fit all the temperature profiles investigated in this work and to calculate the temperature “jump,” ΔT.

The results obtained with the BBA and ANA methods agree within the statistical uncertainty of our computations, showing that our implementation of the local heat flux [see Eq. (4)] produces consistent results with the heat fluxes obtained using a cylindrical or spherical sampling volume [via Eq. (5)]. The radial heat fluxes agree with the results of the continuity equations [see Figs. 5(b) and 5(d)]: continuity equation J_q = q_r/(2πLr) for cylinders and J_q = q_r/(4πR²) for spheres. The continuity result provides an additional consistency check for the heat flux computations using the IK method [Eqs. (4) and (5)]. We find excellent agreement between the continuity equation and the IK methods, supporting the validity of our ANA implementation. Equation (19) models accurately the temperature profiles [see Figs. 5(a) and 5(c)], reproducing the local temperature inside the nanoparticle, fluid, and interface.

We have shown above that the ANA method produces temperature profiles and heat fluxes, consistent with the calculations using common control volumes and the continuity equations. In the following, we use the ANA approach to obtain the local ITC of nanoparticle–fluid interfaces. We will establish a relationship between the ITC and the local structure of the NP. With this purpose, we use the local coordination number (CN) between nanoparticle atoms as a fingerprint characterizing their local environment. In this context, “low” values of CN indicate a more significant number of contacts between that NP atom and the fluid atoms. We note that even for shapes of high symmetry, such as spheres, the nanoscale roughness of the surface results in unique atomic environments and therefore somewhat different local coordination numbers. We illustrate this notion in Fig. 6 for a sphere of diameter 2R* = 40 (∼12 nm using σ_f = 0.3 nm, which is of the order of the water molecule diameter). The surface of this spherical nanoparticle features a heterogeneous environment with the coordination numbers of the surface atoms spanning values between 6 and 9. Under-coordinated atoms with CN values lower than bulk (CN_bulk−NP = 12) feature higher ITCs. This result is consistent with a larger number of contacts with the fluid, which reduces the thermal resistance. Indeed, NP surface atoms with low CN feature higher local fluid densities (see Fig. 13 in the supplementary material).

Figure 6(b) shows a color map of the ITC projected on the NP surface atoms. Our ANA method uncovers regions with contrasting thermal conductance. Facets with high CN = 9 feature significantly lower ITCs than regions with low CN = 6. The differences in conductance between these regions are very significant, ∼50%. The surface atoms with lower conductance experience a larger temperature “jump” and are hotter than the high conducting ones. Atoms with high ITC dissipate heat better and cool down faster. Conversely, the fluid surrounding low ITC regions remains cooler due to poorer heat transfer [see Figs. 6(b) and 6(c)]. Overall, the ANA approach uncovers important details missed by previous methods, such as bulk binning, which only provides an ITC corresponding to the average over regions with very different ITCs. We have compared the thermal conductance obtained from both methods, ANA and BB, by computing the arithmetic mean of the local ITCs obtained with the ANA approach. We obtain G* = 0.74 ± 0.02 in good agreement with the result from the radial bulk binning method, G* = 0.75 ± 0.05, showing the consistency of the ANA method.

We have shown that the ANA allows for calculating local ITCs in relatively simple geometries such as spheres. In the following, we use the ANA method to compute the ITC of complex nanoparticles and resolve the different contributions due to atoms located at facets, edges, and vertices.

First, we investigated an infinite nanowire with an octagonal cross section and an average cylindrical diameter 2R* = 11.1 (3.3 nm using σ = 0.3 nm). This NP contains under-coordinated atoms at edges and high CN atoms at the faces (100) and (110). We show in Fig. 7 the NP structure along with local ITC maps projected onto the NP surface atoms. We have represented our results as a function of the radial distance to the main axis of the NP. Note that atoms belonging to the same facet in this representation might have a different radial distance to the axis. The ITC varies significantly depending on the atom location, and the ITC decrease in the order: edges $>$ (110) facets $>$ (100) facets. Following the discussion in Sec. III B, we can explain the decrease of the ITC in terms of the local CN of the surface atoms. The CN increases in the order 6, 7 for edges, (110), and 8 for (100) facets, and we get consistently higher ITC for (110) faces. This result supports our view that lower CN is connected to higher ITC.

The ITC of the NP atoms located at (110) and (100) faces is very similar to those obtained simulating periodically extended surfaces with that symmetry [see Fig. 7(a)]. This result indicates that the ITC of the nanoscale facets can be approximated with the ITCs of macroscopic surfaces. We also calculated the average ITC of the nanoparticle by using the arithmetic mean of the ITCs and considering the fraction (fr) of surface atoms with different CN: 8 (fr = 20/36), 7 (fr = 8/36), and 6(fr = 8/36). Despite the smaller fraction of atoms with lower CN, they contributed significantly to the average due to their high ITC. We obtained G* = 0.724 ± 0.005 using the arithmetic mean, which is in agreement with the coarse-grained radial bulk-binning method, G* = 0.74 ± 0.01.

It is worth noting the slight increase in the ITC from the center of the (100) facet (lower radial distance) toward the edges (higher radial distance) [see Fig. 7(a)]. All the atoms in this facet have the same CN, 8. These small changes in ITC are therefore not associated with different coordination numbers. The Voronoi volume analysis helps rationalize the observed increase in conductance near the edges as it considers the packing effects of both fluid and NP atoms surrounding each surface atom. Therefore, it can identify subtle variations in the local atomic environment.

The Voronoi volume per atom increases as the atoms in the (100) facet approach the edges (larger radial distance) [see Fig. 7(b)], hence providing a correlation with the slightly higher conductance reported in Fig. 7(a). Consistently, the atoms in the (110) face (CN = 7) and the edges (CN = 8), which have higher ITC, do also feature higher Voronoi volumes. The Voronoi analysis is useful since it shows that atoms in the same facet with the same CN can have slightly different volumes, which ultimately indicates whether the atoms are in the core of the facet (lower volumes) or toward the edges (higher volumes). A higher atomic Voronoi volume is also observed to correlate with greater NP–fluid interactions. Since the fluid has lower density than the NP and surface NP atoms that are uncoordinated with other NP atoms bind more with the less dense fluid, the Voronoi volume of the under-coordinated NP atoms is increased.

Nanoparticles can be synthesized in various complex shapes incorporating in their structure different facets and surface atoms with different coordination numbers. We analyze in the following the ITC of different shapes corresponding to experimental structures. Some of these shapes, e.g., cylinders or particles featuring triangular faces, i.e., with “spikes,” can potentially be of interest in thermal therapy applications.³⁹ 

We correlate the ITC of the different NPs with the average coordination number, CN_avg. Our calculation relies on the α-shape method to identify surface atoms. For the latter, we consider only atoms whose coordination number is $≤9$⁠. We provide a justification for this criterion in Fig. 15 in the supplementary material. We define the CN_avg of a particle as

where the index i runs over coordination numbers 1, …, 9 and s_i is the number of atoms with coordination number i. The number of surface atoms is proportional to the surface area of the NP.

To calculate CN_avg(n), we use the analytical equations derived by Kaatz and Bultheel.⁴⁰ for archimedean and platonic solids whose size and number of atoms are linked to the number of complete shells, n, that it is composed of (see Ref. 40 for further details). The analytical solutions from Kaatz and Bultheel’s equation are collected in Fig. 8 along with a representation of the NP shapes and examples of these shapes taken from experimental studies. The CN_avg(n) of various platonic and archimedean solids can be defined in terms of a rational expression consisting of a second order polynomial in the number of surface shells such that $CNavg(n)=∑i=19isi(n)∑i=19si(n)$⁠. Using CN_avg(n) and the equations provided in Ref. 40, we get the rational expression:

where a, …, e are coefficients of second order polynomials listed in Fig. 8. The denominator in Eq. (24) quantifies the total number of surface atoms with CN ≤ 9. The ratio of coefficients a/d quantifies the highest surface coordination environment attainable for a particle in the limit n → ∞. Shapes with lim_n→∞CN_avg(n) = 9.0 (tetrahedron, octahedron, dodecahedron, and decahedron) consist of (111) facets only.

Non-archimedean and non-platonic shapes such as spheres and cylinders do not possess such magic numbers. Hence, the numerator of Eq. (23), $∑i=19isi$⁠, is a second order polynomial in r for spheres (4πr²) and first order for cylinders of length L, (2πrL), where r is the NP spherical and cylindrical radii. The resulting equation, $CNavg(r)=ar2+br+cdr2+er+1$⁠, resembles Eq. (24) where a, …, e are coefficients obtained by fitting the dependence of the average surface coordination number with the number of atoms in the nanoparticle (see Fig. 9).

CN_avg(n) increases in a similar fashion with the number of atoms for nanoparticles of different shapes (see Fig. 8). Nanoparticles with 8.0 < lim_n→∞CN_avg(n) < 9.0 (cuboctahedron and truncated octahedron) have both (100) and (111) faces. The ratio of the surface area of these faceted surfaces determines lim_n→∞CN_avg(n). The hexagonal bipyramid NPs consist of (111) and (210) facets, resulting in a lim_n→∞CN_avg(n) = 8.5, an average of coordination numbers 9 and 8. Cubes are constructed using (100); hence, lim_n→∞CN_avg(n) = 8.0. Finally, the shape of the rhombic dodecahedron (RD) is dominated by (110), resulting in lim_n→∞CN_avg(n) = 7.

We show in Fig. 9 the dependence of the average coordination number of the surface atoms for different shapes taking into account the NP size. The CN_avg(n) increases with the number of atoms in the NP. This increase is monotonic, except for some shapes (see Fig. 17 in the supplementary material). Following the results from Sec. III C, we therefore expect that the ITC will decrease as the particle size increases. We also find that the rhombic dodecahedron features the lowest CN_avg(n), while icosahedrons, decahedrons, tetrahedrons, and octahedrons possess the highest CN_avg(n). Advancing the discussion in the following and considering our analysis in Sec. III C, we expect that these structures will show the lowest (icosahedron) or highest (rhombic dodecahedron) ITCs. We conclude from this analysis that among the various NP types investigated, spherical, cylindrical, and rhombic dodecahedron shapes are likely to have higher conductances.

We used the ANA method to compute the local ITC of different NP structures and sizes and the BBA for three flat faceted surfaces [see Fig. 10(a)]. The ITC varies inversely with the average CN. This idea can be made concrete by representing the Kapitza resistance (inverse ITC) as a function of CN_avg(n) [see Fig. 10(b)]. The NP structures represented in the insets highlight regions of high and low conductance. Under-coordinated atoms at vertices and edges have much higher ITCs than highly coordinated atoms located at flat surfaces. Our results show a linear correlation between the Kapitza resistance and CN_avg(n). NP with lower CN_avg(n) conducts heat better than particles with high CN_avg(n) with ratios between maximum ITC and minimum ITC of 19/11 $∼1.7$⁠. These differences are substantial considering the particles have approximately the same average diameter (10 nm) (see Sec. 1.2 in the supplementary material for information on the number of atoms in each nanoparticle). The results presented in Fig. 9 highlight the importance of the CN as a variable to rationalize the ITC of particles with very complex shapes.

Shapes with multiple facets such as the truncated octahedron, cuboctahedron, nanowire, and nanorod show multi-modal distributions due to the different thermal properties of different facets as well as a small distribution of high conductance atoms from edge atoms located near the vertices. Shapes with a single facet such as a cube, rhombic dodecahedron, and octahedron show higher homogeneity in their conductance distributions with a dominant distribution reflecting the thermal property of the facet that constitutes the shape. For a nanosphere that features a variety of facets present from cutting an FCC lattice into a sphere, a large degree of heterogeneity can be seen from the overlapping distributions of interfacial thermal conductance (see Fig. 20 in the supplementary material).

The dependence of the ITC with particle size has attracted attention recently.^{17,18,27–29} Several of these studies showed an increase of the ITC with interfacial curvature or decreasing particle size. Interestingly, an increase in the ITC with particle size has been reported also.²⁷ We showed in Fig. 9 that the average coordination number (CN) generally increases with particle size, and in Subsection III D, we established a relationship between higher thermal conductance and lower CN. These results suggest that the ITC will generally decrease with increasing particle size.

To test the size dependence of the ITC, we performed additional simulations of the octahedron, spheres, and rhombic dodecahedron (RD) nanoparticles. Our results (see Fig. 11) indicate that the RD features the lowest interfacial thermal resistance (highest ITC), while the octahedron features a high thermal resistance. This general dependence is in line with the results discussed in Subsection III D with the NP featuring a higher average CN showing lower ITCs. The spherical particles do also show a higher resistance for higher average CN or with larger particle size. This result is consistent with previous observations.^17,18,28,29

We have further explored the correlation between fluid adsorption and thermal conductance. This correlation was studied before in Refs. 17 and 24, showing that the ITC increases linearly with interfacial adsorption. We calculated the adsorption using

where N_l is the number of particles in the first solvation shell, A* is the NP surface area, ρ*(r) is the fluid density at distance r*, and M₁ represents the first minimum of the density profile (see Fig. 16 in the supplementary material for plots of the radial density profiles). We obtain an approximately linear dependence between the ITC and the adsorption for the spherical particles investigated here (see Fig. 11).

We have introduced a computational approach to calculate the Interfacial Thermal Conductance (ITC) of nanomaterial–fluid interfaces. The Atomistic Nodal Approach (ANA) relies on combining “computational farming,” namely, using thousands of independent non-equilibrium molecular dynamics stationary simulations. We also use tools from computational geometry (KD-trees, Voronoi polyhedra, and α-shapes) to identify surface atoms in nanoparticles and their solvation environment, which is quantified through their local coordination number. Local thermal gradients were computed via multivariate trilinear interpolation. Furthermore, the atoms in the nanomaterials are used as a nodal network to project heat fluxes, local densities and temperature gradients, and compute the interfacial thermal conductance at atomistic resolution.

We have tested the performance of the ANA method against widely used binning methods that partition the system of interest in well-defined geometries, such as spheres and cylinders. These geometries are amenable for theoretical treatments using the classical heat transfer theory and the heat diffusion equation. The ANA produces results entirely consistent with the binning methods and goes beyond these methods by enabling the calculation of local properties at the atomic level. Our approach makes it possible to identify regions in the surface of the nanomaterials with contrasting ITC. Moreover, it enables us to quantify the thermal conductance of nanomaterials with complex shapes (nanocrystals with facets, vertices, and edges) and, in general, surfaces with significant roughness, which cannot be easily mapped to standard geometries. We note that our method is general, and it can be extended to more complex systems. Calculating the heat flux in molecular systems and more complex forcefields, e.g., for metallic solids, requires proper consideration of intramolecular degrees of freedom and non-pairwise potentials.

We applied the ANA method to nanoparticles of different shapes and sizes immersed in an atomistic fluid. The nanoparticles were heated at stationary conditions, and the heat transport across the nanoparticle–fluid interface was analyzed using the ANA method. We conclude the following from our investigation:

• We have established a correlation between the ITC and the coordination number (CN) of atoms in the surface of the nanoparticle, where the CN quantifies the number of nearest neighbors of surface atoms with other atoms in the nanoparticle. The ITC increases with decreasing CN. Surface atoms with smaller CN share more contacts with the solvent. These contacts enhance the thermal conduction from the hot nanoparticles to the fluid. Generally, undercoordination correlates positively with higher ITC. Hence, the coordination of the surface atoms reveals itself as an approach to design nanoheaters for which the heating performance can be tuned by changing the particle geometry and modifying at the same time their ITC.

Thermal therapies are of interest in cancer treatments. The most efficient nanostructures provide efficient light/heat conversion, absorption in the so-called biological windows, and high interfacial conductances for better heat transfer from the nanoparticle to the surrounding fluids. Our work tackles the latter point. We infer from our results that low surface coordination offers better prospects for these applications [see Figs. 10(a) and 10(b)].

• For nanoparticles featuring facets with well-defined symmetries (111), (100), and (110), the local ITC obtained for small surfaces spanning a few nm² is very similar to the CN of extended infinite surfaces. This is an interesting result that provides a design principle to estimate the ITC of nanoparticles using the ITC of macroscopic surfaces. Furthermore, the ITC of atoms at edges and vertices can be obtained readily using the ANA approach, making it possible to predict the ITC of particles for a wide range of geometries.

• Nanoparticle size and morphology act synergistically and determine the average ITC of nanoparticles. Smaller NPs have higher thermal conductances due to a more significant proportion of under-coordinated atoms at the surface. Hence, small nanoparticles are better heat conductors. Instead, larger NPs form surfaces with a higher average surface coordination number, reducing their ITC. The general dependence of the average surface coordination number with the number of NP atoms can be described using a rational expression.

• The Kapitza resistance (the inverse ITC) of nine different nanoparticle structures varies linearly with the average CN of the surface atoms. Nanoparticles mainly consisting of (100) and (110) faces (rhombic dodecahedron and nanocubes) feature the lowest thermal resistance. Those with (111) faces (icosahedron and dodecahedron) have the highest ITC with spherical and cylindrical geometries having intermediate conductances.

We expect the atomistic nodal approach proposed here will stimulate further theoretical and experimental work to advance our understanding of thermal transport at the nanoscale and its dependence on the particle size. In addition, knowledge of local ITC should be helpful to rationalize experimental observations and assist the design of nanoparticles for biomedical and, more generally, thermoplasmonic applications.

Additional data on simulation details, calculation of heat profiles, ITCs of various nanoparticles shapes, and local solvation structures of surface atoms are available in the supplementary material.

We acknowledge the Leverhulme Trust (UK, Grant No. RPG-2018-384) for financial support and the Imperial College RCS High Performance Computing facility and the UK Materials and Molecular Modelling Hub for computational resources partially funded by the EPSRC (Grant Nos. EP/P020194/1 and EP/T022213/1). F.B. acknowledges financial support from the EPSRC (Grant No. EP/J003859/1).

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request. The code and script used for generating and averaging data are publicly available in github at https://github.com/bresmegroup/Atomistic_Nodal_Approach_Nanoparticles.git.