The Interfacial Thermal Conductance (ITC) is a fundamental property of materials and has particular relevance at the nanoscale. The ITC quantifies the thermal resistance between materials of different compositions or between fluids in contact with materials. Furthermore, the ITC determines the rate of cooling/heating of the materials and the temperature drop across the interface. Here, we propose a method to compute local ITCs and temperature drops of nanoparticle–fluid interfaces. Our approach resolves the ITC at the atomic level using the atomic coordinates of the nanomaterial as nodes to compute local thermal transport properties. We obtain high-resolution descriptions of the interfacial thermal transport by combining the atomistic nodal approach, computational geometry techniques, and “computational farming” using non-equilibrium molecular dynamics simulations. We use our method to investigate the ITC of nanoparticle–fluid interfaces as a function of the nanoparticle size and geometry, targeting experimentally relevant structures of gold nanoparticles: capped octagonal rods, cuboctahedrons, decahedrons, rhombic dodecahedrons, cubes, icosahedrons, truncated octahedrons, octahedrons, and spheres. We show that the ITC of these very different geometries varies significantly in different regions of the nanoparticle, increasing generally in the order face < edge < vertex. We show that the ITC of these complex geometries can be accurately described in terms of the local coordination number of the atoms in the nanoparticle surface. Nanoparticle geometries with lower surface coordination numbers feature higher ITCs, and the ITC generally increases with the decreasing particle size.
I. INTRODUCTION
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.6
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,12 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.1 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.17
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 = Jq/ΔT, by computing explicitly the interfacial temperature jump, ΔT, and interfacial heat fluxes, Jq. 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,17 uncovering a correlation between high nanoscale curvature and high thermal conductance. This observation supports earlier studies of fluid–fluid19 and nanoparticle–fluid interfaces28 as well as theoretical models.29 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.27 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.
II. METHODS
A. Nanoparticle models and simulation details
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 , 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°.17 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,17 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.30 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 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 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 with the NP density set to . The temperature T of each particle of mass mi and velocity vi in the simulation was determined using the equipartition equation
where kB is the Boltzmann constant and Ndof = 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* = kBT/ɛf. To determine the heat flux at each position, we used a local version of the formula derived by Irving and Kirkwood,31
where Jq,i is the heat flux corresponding to particle i with associated volume Vi, 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 rij, 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, Vi, which we identify with the Voronoi volume that is calculated for each atom as a function of time. Practically, the per-atom volume, Vi, was obtained using the Voronoi construction of both fluid and nanoparticle atoms in SNEMD simulations spanning 105 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 106 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 × 104 timesteps were discarded, and the remaining 5 × 104 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
B. Bulk binning approach
Previous computational studies of the ITC17,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 × 104 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).
Illustration of BBA for (a) spherical or cylindrical (top view) shells and (b) slabs for flat surfaces. These methods use a smooth and continuous shell to approximate the surface properties.
Illustration of BBA for (a) spherical or cylindrical (top view) shells and (b) slabs for flat surfaces. These methods use a smooth and continuous shell to approximate the surface properties.
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.
C. Temperature and heat flux profiles
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 qr 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 across various radii using the Irving–Kirkwood formula and by examining its relationship with the numerical thermal gradient calculated from simulated data.
1. Continuity equation for nanorods
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 rc and the temperature T difference from a reference point Tc. The heat rate experienced by the NP and fluid can be expressed with qr and the thermal conductivity of the particle (p) or fluid (f), κ,
The heat rate experienced by the NP and fluid can be expressed with qr and is the same for both NP and fluid due to the conservation of the energy in a steady-state system.
where rc is the radial distance defining the temperature jump, (see Subsection II C 3).
2. Continuity equation for nanospheres
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
3. Using the error function at the NP–fluid interface
To model the temperature transition between the NP Tp(r) and fluid Tf(r) and account for their different thermal conductivities κp and κl and the interfacial temperature “jump” we used an error function equation,
where rc refers to the transition point and ω quantifies the steepness of the error function .34 The temperature “jump” ΔT can be found using ΔT = Tp(rc) − Tf(rc), which is defined in terms of the temperature of the particle, , and the fluid, , at rc.
D. The atomistic nodal approach
1. The philosophy behind ANA
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.
(a) Illustration of the fluid (blue) and NP atoms (yellow) and the heat flux vectors (red arrows) for a rhombic dodecahedron gold nanoparticle. The NP atoms are treated as nodes. A cubic grid of length 0.4σf is superimposed on the fluid region to compute the properties of the fluid atoms. (b) The surface (yellow) of the NP identified by the α-shape algorithm. The atoms in the nanoparticle solvation shell are highlighted in red.
(a) Illustration of the fluid (blue) and NP atoms (yellow) and the heat flux vectors (red arrows) for a rhombic dodecahedron gold nanoparticle. The NP atoms are treated as nodes. A cubic grid of length 0.4σf is superimposed on the fluid region to compute the properties of the fluid atoms. (b) The surface (yellow) of the NP identified by the α-shape algorithm. The atoms in the nanoparticle solvation shell are highlighted in red.
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 . The probe radius is an important variable in our method. Our choice () 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.37
Density profile of the fluid surrounding the NP. The high density fluid regions are represented in yellow, and the low density regions are in purple. The NP surface atoms are colored in red, and the rest of the NP atoms are in gray. Orange lines joining NP surface atoms with fluid atoms indicate the nearest fluid-neighbors identified.
Density profile of the fluid surrounding the NP. The high density fluid regions are represented in yellow, and the low density regions are in purple. The NP surface atoms are colored in red, and the rest of the NP atoms are in gray. Orange lines joining NP surface atoms with fluid atoms indicate the nearest fluid-neighbors identified.
Subsequently, the temperature gradient 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).38 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 whose radii rNP − rF are within a fixed distance M1 of the atom in the NP, where M1 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 TNP,i and TF are the per-atom NP (atomic node, i) and fluid grid temperatures, respectively.
The ITC was obtained by using the heat flux, Jq, 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 .
2. General ANA simulation approach
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 103 simultaneous independent simulations, 5 × 104 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.
Algorithmic approach used to analyze the NEMD data generated with the ANA method. The ANA method makes use of the computational geometry methods, α-shapes to identify surface atoms, and kd-trees to find nearest neighbors for the NP surface-atom nodes.
Algorithmic approach used to analyze the NEMD data generated with the ANA method. The ANA method makes use of the computational geometry methods, α-shapes to identify surface atoms, and kd-trees to find nearest neighbors for the NP surface-atom nodes.
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.
III. RESULTS AND DISCUSSION
A. Comparison between BBA and ANA approaches
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.
Dependence of the temperature [panels (a) and (c)] and heat flux [panels (b) and (d)] profiles with radial distance for cylindrical and spherical nanoparticles obtained with the ANA and BBA methods. The red and blue areas represent the location of the hot (T* = 0.92) and cold thermostats (T* = 0.72), respectively. The temperature profiles and the heat fluxes obtained with the BBA and ANA approaches agree with the profiles modeled using Eq. (19) for the fluid and the nanoparticle. We show in the insets of panels (a) and (c) a cylindrical NP of radius r* = 5.6 and a spherical NP of radius R* = 5.0.
Dependence of the temperature [panels (a) and (c)] and heat flux [panels (b) and (d)] profiles with radial distance for cylindrical and spherical nanoparticles obtained with the ANA and BBA methods. The red and blue areas represent the location of the hot (T* = 0.92) and cold thermostats (T* = 0.72), respectively. The temperature profiles and the heat fluxes obtained with the BBA and ANA approaches agree with the profiles modeled using Eq. (19) for the fluid and the nanoparticle. We show in the insets of panels (a) and (c) a cylindrical NP of radius r* = 5.6 and a spherical NP of radius R* = 5.0.
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 Jq = qr/(2πLr) for cylinders and Jq = qr/(4πR2) 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.
B. Relationship between coordination number and conductance: The spherical nanoparticle
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 (CNbulk−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).
(a) Dependence of the ITC with the local CN and radial distance for a spherical nanoparticle of R* = 20. The color maps indicate different thermal conductances, G*, or radial distances, R*. (b) The ITC results have been projected on the NP surface atoms to identify regions of high and low conductance. (c) Temperature distribution of the NP atoms with hotter regions corresponding to atoms with low thermal conductance and vice versa. The size of the points in panel (a) indicates the ratio of surface atoms with a particular CN and radial distance R* from the NP center.
(a) Dependence of the ITC with the local CN and radial distance for a spherical nanoparticle of R* = 20. The color maps indicate different thermal conductances, G*, or radial distances, R*. (b) The ITC results have been projected on the NP surface atoms to identify regions of high and low conductance. (c) Temperature distribution of the NP atoms with hotter regions corresponding to atoms with low thermal conductance and vice versa. The size of the points in panel (a) indicates the ratio of surface atoms with a particular CN and radial distance R* from the NP center.
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.
C. ITC of complex geometries: Induced edge effects on conductance
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.
(a) Dependence of the ITC with the radial distance of the atoms to the main axis of nanowire periodic in the axial direction (infinite length). The dashed horizontal lines represent the ITC obtained with independent simulations of periodic (infinite) surfaces spanning the whole simulation cell. The inset shows a cross section of the nanoparticle. The heat map indicates the ITC of the atoms at the surface of the NP. The bulk atoms in the NP are shown as gray spheres. (b) Voronoi volumes as a function of the radial distance. The inset features the fluid atoms (red) and NP atoms (blue) and their corresponding Voronoi tessellation.
(a) Dependence of the ITC with the radial distance of the atoms to the main axis of nanowire periodic in the axial direction (infinite length). The dashed horizontal lines represent the ITC obtained with independent simulations of periodic (infinite) surfaces spanning the whole simulation cell. The inset shows a cross section of the nanoparticle. The heat map indicates the ITC of the atoms at the surface of the NP. The bulk atoms in the NP are shown as gray spheres. (b) Voronoi volumes as a function of the radial distance. The inset features the fluid atoms (red) and NP atoms (blue) and their corresponding Voronoi tessellation.
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.
D. Thermal conductance of complex nanoparticle shapes
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.39
We correlate the ITC of the different NPs with the average coordination number, CNavg. Our calculation relies on the α-shape method to identify surface atoms. For the latter, we consider only atoms whose coordination number is . We provide a justification for this criterion in Fig. 15 in the supplementary material. We define the CNavg of a particle as
where the index i runs over coordination numbers 1, …, 9 and si 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 CNavg(n), we use the analytical equations derived by Kaatz and Bultheel.40 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 CNavg(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 . Using CNavg(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 limn→∞CNavg(n) = 9.0 (tetrahedron, octahedron, dodecahedron, and decahedron) consist of (111) facets only.
Analytical solutions for CNavg(n) as a function of the number of shells, n, for 9 common NP structures; n ≥ 1 in all cases. Note that for the rhombic dodecahedron, only even values of n are considered. To illustrate the NP shapes, we represent cartoons generated with the code polyhedral-viewer (https://github.com/tesseralis/polyhedra-viewer). Images of corresponding NPs synthesized in experiments are also shown (left panel next to each cartoon). The experimental images are taken from the following references and are reproduced with permission. (a) Reproduced with permission from Zheng et al., Chem. - Asian J. 9, 2635–2640 (2014). Copyright 2014 John Wiley & Sons.41 (b)–(e) and (h) Reproduced with permission from Hsu et al., Inorg. Chem. 49, 4149–4155 (2010). Copyright 2010 American Chemical Society.42 (f) Reproduced with permission from Ma et al., ACS Nano 14, 9594–9604 (2020). Copyright 2020 American Chemical Society.43 (g) Reproduced with permission from Kim et al., Angew. Chem., Int. Ed. 43, 3673–3677 (2004). Copyright 2004 John Wiley & Sons.44 (i) Reproduced with permission from Personick et al., J. Am. Chem. Soc. 135, 3800–3803 (2013). Copyright 2013 American Chemical Society.45
Analytical solutions for CNavg(n) as a function of the number of shells, n, for 9 common NP structures; n ≥ 1 in all cases. Note that for the rhombic dodecahedron, only even values of n are considered. To illustrate the NP shapes, we represent cartoons generated with the code polyhedral-viewer (https://github.com/tesseralis/polyhedra-viewer). Images of corresponding NPs synthesized in experiments are also shown (left panel next to each cartoon). The experimental images are taken from the following references and are reproduced with permission. (a) Reproduced with permission from Zheng et al., Chem. - Asian J. 9, 2635–2640 (2014). Copyright 2014 John Wiley & Sons.41 (b)–(e) and (h) Reproduced with permission from Hsu et al., Inorg. Chem. 49, 4149–4155 (2010). Copyright 2010 American Chemical Society.42 (f) Reproduced with permission from Ma et al., ACS Nano 14, 9594–9604 (2020). Copyright 2020 American Chemical Society.43 (g) Reproduced with permission from Kim et al., Angew. Chem., Int. Ed. 43, 3673–3677 (2004). Copyright 2004 John Wiley & Sons.44 (i) Reproduced with permission from Personick et al., J. Am. Chem. Soc. 135, 3800–3803 (2013). Copyright 2013 American Chemical Society.45
Non-archimedean and non-platonic shapes such as spheres and cylinders do not possess such magic numbers. Hence, the numerator of Eq. (23), , is a second order polynomial in r for spheres (4πr2) and first order for cylinders of length L, (2πrL), where r is the NP spherical and cylindrical radii. The resulting equation, , 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).
Dependence of CNavg (CN 10) with NP size (in terms of number of atoms in the NP) for 9 FCC faceted solids, cylinders, and spheres. The total number of atoms in the NP varies with the number of shells as n3 (see Sec. 1 in the supplementary material for particle sizes). A spherical NP of 3 × 104 gold atoms has a diameter of ∼10 nm. A plot with a log scale for the number of atoms is shown in Fig. 17 in the supplementary material.
Dependence of CNavg (CN 10) with NP size (in terms of number of atoms in the NP) for 9 FCC faceted solids, cylinders, and spheres. The total number of atoms in the NP varies with the number of shells as n3 (see Sec. 1 in the supplementary material for particle sizes). A spherical NP of 3 × 104 gold atoms has a diameter of ∼10 nm. A plot with a log scale for the number of atoms is shown in Fig. 17 in the supplementary material.
CNavg(n) increases in a similar fashion with the number of atoms for nanoparticles of different shapes (see Fig. 8). Nanoparticles with 8.0 < limn→∞CNavg(n) < 9.0 (cuboctahedron and truncated octahedron) have both (100) and (111) faces. The ratio of the surface area of these faceted surfaces determines limn→∞CNavg(n). The hexagonal bipyramid NPs consist of (111) and (210) facets, resulting in a limn→∞CNavg(n) = 8.5, an average of coordination numbers 9 and 8. Cubes are constructed using (100); hence, limn→∞CNavg(n) = 8.0. Finally, the shape of the rhombic dodecahedron (RD) is dominated by (110), resulting in limn→∞CNavg(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 CNavg(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 CNavg(n), while icosahedrons, decahedrons, tetrahedrons, and octahedrons possess the highest CNavg(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 CNavg(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 CNavg(n). NP with lower CNavg(n) conducts heat better than particles with high CNavg(n) with ratios between maximum ITC and minimum ITC of 19/11 . 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.
(a) Variation of ITC with average coordination number for each individual unique atomic CN across various NP shapes. (b) Dependence of the Kapitza resistance with the average coordination number of the NP surface atoms for various NP shapes. The ITC has been projected on the surface atoms for each shape with red and orange colors representing atoms with high conductance and blue representing atoms with low ITC (see the color bar on the right-hand side of the figure). The NPs used in both plots have sizes in the range of 3 × 103 to 3 × 104 atoms. The size of each point represents the number of atoms in each NP type (see the supplementary material Sec. 1.2 for more details). The gradient and intercept of this plot are 0.34 ± 0.03 and −1.2 ± 0.2, respectively.
(a) Variation of ITC with average coordination number for each individual unique atomic CN across various NP shapes. (b) Dependence of the Kapitza resistance with the average coordination number of the NP surface atoms for various NP shapes. The ITC has been projected on the surface atoms for each shape with red and orange colors representing atoms with high conductance and blue representing atoms with low ITC (see the color bar on the right-hand side of the figure). The NPs used in both plots have sizes in the range of 3 × 103 to 3 × 104 atoms. The size of each point represents the number of atoms in each NP type (see the supplementary material Sec. 1.2 for more details). The gradient and intercept of this plot are 0.34 ± 0.03 and −1.2 ± 0.2, respectively.
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).
E. Size dependence of the ITC
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.27 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
Dependence of the (a) Kapitza resistance and (b) ITC with the average coordination number of the NP surface atoms for various NP shapes and sizes. The number of atoms in the NP is visualized by the marker size of each particle type. (c) Relationship between the ITC of spheres of various sizes and the adsorption of the first fluid layer (see Sec. 1 in the supplementary material for information on particle sizes). (d) Dependence of the adsorption with the NP curvature 1/R* of the nanospheres.
Dependence of the (a) Kapitza resistance and (b) ITC with the average coordination number of the NP surface atoms for various NP shapes and sizes. The number of atoms in the NP is visualized by the marker size of each particle type. (c) Relationship between the ITC of spheres of various sizes and the adsorption of the first fluid layer (see Sec. 1 in the supplementary material for information on particle sizes). (d) Dependence of the adsorption with the NP curvature 1/R* of the nanospheres.
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 Nl 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 M1 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).
IV. CONCLUSIONS
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 nm2 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.
SUPPLEMENTARY MATERIAL
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.
ACKNOWLEDGMENTS
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).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
DATA AVAILABILITY
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.