Molecular dynamics investigation of surface roughness scale effect on interfacial thermal conductance at solid-liquid interfaces

Non-equilibrium molecular dynamics simulations were conducted for solid-liquid-solid systems with nanometer scale grooved surfaces and an induced heat flux for a wide range of topology and solid-liquid interaction conditions to investigate the mechanism of solid-liquid heat transfer, which is the first work of such extensive detail done about the nanoscale roughness effect on heat transfer properties. Single-atom molecules were used for liquid, and the solid-liquid interaction was varied from superhydrophobic to superhydrophilic, while the groove scale was varied from single atom to several nanometers, while keeping the surface area twice that of a flat surface. Both Wenzel and Cassie wetting regimes with a clear transition point were observed due to the capillary effect inside larger grooves that were more than 5 liquid molecule diameters, while such transition was not observed at smaller scales. At the hydrophobic state, large scale grooves had lower interfacial thermal conductance (ITC) due to the Cassie regime, i.e., having unfilled grooves, while at the hydrophilic state, grooved surfaces had ITC about twice that of a flat surface, indicating an extended heat transfer surface effect regardless of the groove scale. At the superhydrophilic state, crystallization of liquid at the surface occurred, and the packing of liquid molecules had a substantial effect on ITC regardless of the groove scale. Finally, both potential energy of solid-liquid interaction and work of solid-liquid adhesion were calculated and were shown to be in similar relations to ITC for all groove scales, except for the smallest single-atom scale grooves, due to a different heat transfer mechanism. © 2019 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5081103


I. INTRODUCTION
As manufacturing technology has improved, electronic device size has greatly reduced, while their performance vastly increased. One of the challenges that have to be solved is an efficient cooling solution, as with size reduction, heat dissipation density increases, while also putting restriction on the heat sink. This is an acute problem with microchips, 1,2 as well as power modules inside electronic vehicles. [3][4][5] The heat transfer efficiency between the heat source and sink is greatly affected by the properties of their interface because of the small device dimensions.
In general, the actual contact area between two solid surfaces is a small fraction of the total apparent contact area. This essentially restricts heat flow to the contact area as the thermal conductivity of air is too low for surface-gap-surface heat transfer to have a significant contribution. Thermal interface materials (TIMs) are often placed between surfaces to reduce the gap amount, where the major types are particle-laden polymers, metal, and carbon nanotube array TIMs. 6 Much experimental work has been done on finding TIM with good thermal conductivity, low interfacial thermal resistance, while also fulfilling other requirements such as presence/lack of electric conductivity, high viscosity at liquid state, or low young modulus at solid state, as well as economical viability and regulatory compliance. 6,7 Reducing the interfacial thermal resistance is an important requirement for high performance TIM, which becomes especially critical as TIM thickness decreases, 8 and a better theoretical understanding of the effect would be greatly beneficial. At the basic level, solid-TIM interaction can be likened to solid-liquid or solid-solid interactions. Two theoretical models, "acoustic mismatch model" (AMM) and "diffusive mismatch model" (DMM), 9 assuming either specular phonon reflection and refraction, or complete diffusive scattering, respectively, are widely used to predict interface heat transfer properties. Unfortunately, actual interfaces usually do not fully conform to either of the models, and molecular dynamics (MD) simulation has been found as an effective tool to investigate interfacial properties at greater detail. [10][11][12][13][14] At the macroscopic scale of solid-liquid interfaces, interfacial thermal conductance (ITC) can be increased by attaching fins or creating grooves at the surface 15 as this increases the heat transfer surface area. Additionally, it has been demonstrated via molecular dynamics and lattice Boltzmann methods that the surface roughness at the microscopic scale has profound effects on momentum and heat transfer properties, [16][17][18] and there have also been several studies specifically investigating if this extended heat transfer surface effect is valid for surface topology of the nanometer scale. [19][20][21][22][23] A general trend was shown that a grooved or rough surface increases ITC as long as there is good contact between surface and liquid, implying interaction that results in good wettability. Unfortunately, previous studies either dealt with a single surface geometry or a single type of solid-liquid interaction strength, making it difficult to uncover a general trend between the improvement in ITC, groove scale, and the nature of solid-liquid interactions.
In this work, the surface groove scale at the interface was varied while keeping the surface area constant, while also changing the solid-liquid interaction from superhydrophobic to superhydrophilic to determine the effect on interface heat transfer properties and investigate the heat transfer mechanism. Non-equilibrium molecular dynamics (NEMD) simulations of solid-liquid-solid systems were conducted with an induced heat flux, where single-atom molecules were used to model a liquid, for simplicity. Simulation method details are given in Sec. II A, while the actual details of systems and their creation are described in Sec. II B. Sections III A 1 and III A 2 give details on the system states at various groove scales and solid-liquid interaction strengths, showing the presence of temperature jump and capillary effect at sub-nanometer groove scales. Section III B contains the obtained ITC values for all conditions, while Sec. III C obtains the work of solid-liquid adhesion for each of the systems and discusses the relationship with ITC.

II. SIMULATION METHODS
The Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 24 package was used to conduct molecular dynamics simulations with velocity Verlet integrator at a time step of 5 fs.

A. Potentials
Single atom molecules were used to represent the liquid, with 12-6 Lennard-Jones (LJ) interactions between them, where rij is the distance between atoms i and j, while ε and σ are the potential well depth and distance parameters of the LJ potential, for which the parameters of argon were used. 25 Specific values are displayed in Table I. The solid surface was represented by a FCC lattice, with the Morse potential acting between lattice atoms, where the "ss" superscript represents the solid-solid interactions and D, α, and r 0 are the potential well depth, well width, and equilibrium bond length parameters for the Morse potential, respectively. The parameters derived by Pamuk and Halicioǧlu that reproduce gold crystal properties at room temperature were used and are shown in Table I. 26 Note that because this potential does not express free conduction electrons, heat is transferred strictly via phonons, making it more similar to the non-metallic lattice with respect to thermal properties. The Morse potential was also used to represent the interaction between solid and liquid atoms, where the "sl" superscript represents the solid-liquid interactions and the η parameter is the solid-liquid interaction coupling parameter, with values ranging from 0.01 to 1. The Morse potential parameters were set according to the work done by Grenier et al., 27 where they obtained pair potentials for interaction between the gold surface and argon atoms via quantum chemistry computations. The pair potential obtained from CCSD(T)/AV5Z computations was used to determine the Morse potential well depth and bond length parameters as D sl = 13.583 meV and r sl 0 = 4.1024 Å, respectively. The potential well width parameter was obtained via nonlinear least-squares Marquardt-Levenberg fitting in the range of 3-10 Å provided by the Gnuplot plotting program, 28 which resulted in α sl = 1.3787 Å −1 . The values are summarized in Table I.
The cut-off distance was set to 12 Å for all the interactions.

Bulk lattice systems
Bulk simulations of the solid surface lattice were conducted to determine the lattice constant a of the FCC crystal. A constant The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp number of particles, constant pressure, and constant temperature (NPT) simulation of 4000 atoms in a FCC lattice at 1 atm and 110 K with a periodic boundary condition in all directions was conducted. Nosé-Hoover style equations were used, 29 with the thermostat and barostat each having three chains and damping coefficients set to 0.5 and 5 ps, respectively. The initial system dimensions were determined from a gold density of 19.3 g/cm 3 . A simulation of 10 ns was conducted, and the last 9 ns were used to determine the lattice constant to be a = 4.064 01 (0.000 01) Å, where the standard error of the mean is shown in brackets. This lattice constant was used when constructing grooved surfaces described in Sec. II B 3.

Liquid layer systems
Liquid layer simulations were conducted to determine the saturated vapor pressure and surface tension. A total of 2000 liquid atoms were positioned in a FCC lattice of approximately 40 × 40 × 50 Å 3 size, with the lattice constant set to 2 2 3 σ. These atoms were placed inside a vacuum system of 40 × 40 × 300 Å 3 size, with an all periodic boundary condition. A constant particle number, constant temperature, constant volume simulation (NVT) was conducted with a Langevin thermostat 30 with a damping factor of 100 fs and a control temperature of either 100, 110, or 120 K. A total of 150 ns of simulation were performed for each control temperature. Liquid layers formed over the xy plane creating a system with coexisting liquid and vapor phases. The last 100 ns of the simulations were used for pressure and surface tension computation. The vapor pressure resulted in approximately 0.39 (0.03), 0.79 (0.03), and 1.34 (0.04) MPa for the 100, 110, and 120 K systems, respectively. The surface tension was computed from the pressure tensor, 31 where the pairwise interaction from the virial pressure was used, 32 and resulted in approximately γ lv = 17.36(0.17), γ lv = 12.77(0.14), and γ lv = 8.62(0.17) mN/m for the 100, 110, and 120 K systems, respectively, where γ lv is the surface or liquid-vapor interface tension. The standard error of the mean is shown in brackets for both pressure and surface tension values and was calculated by taking statistical inefficiency into account. 33

Grooved surfaces
The concept of a grooved surface is displayed in Fig. 1(a), where on top of 17 initial layers with the FCC (100) crystal plane perpendicular to the z axis, grooves are created according to a set pitch value. Grooves with depth and width at half-pitch were created, thus effectively doubling the surface area. This scheme allowed creation of grooves at various scales, while maintaining the surface area the same, thus enabling the evaluation of groove scale effect on heat transfer. Surfaces with pitches of 1, 3, 5, 10, and 20 were crated in addition to a flat surface, used as a reference, where the pitch values are in lattice constant units described in Sec. II B 1. The exact cross section dimensions for each system type are displayed in Table II. Because the liquid atom van der Walls diameter is σ = 3.405, 25 groove widths can be interpreted as approximately 0.5, 1.5, 3, 6, and 12 liquid molecule diameters for pitch 1, 3, 5, 10, and 20, respectively.
Note that for pitch 1 systems, groove width is half of the lattice constant, i.e., the depth and width is a single atom. Because of the FCC (100) crystal face, this results in a surface populated by Side-view snapshots of all system types used in this work, where the pitch size of each system is given below in lattice constant units and will be used to describe the system types. The liquid particles are half-transparent to give a better view of the surface grooves. [(c) and (d)] Snapshot of a solid-liquid interface with poor and good wettability, where the groove type is pitch 20 and the solid-liquid interaction coupling parameter is η = 0.1 and η = 0.5, respectively. The visualisations were produced by the PyMOL package. 34 dangling atoms, and pitch 1 surfaces cannot be strictly considered grooved.
Only periodicity in the x and y directions was assumed.

Solid-liquid-solid systems
A flat surface and the grooved surfaces described in Sec. II B 3 were placed at the bottom of the systems together with their mirror images at the top and liquid atoms were placed in between, as displayed in Fig. 1(b). The periodic boundary condition was set in the x and y directions, while no restrictions were set in the z direction. The second outermost layers of both of the solid walls were coupled with Langevin thermostats 30 with damping coefficients of 100 fs and control temperatures set to 100 and 120 K for the bottom  and top thermostats, respectively, which is a common method to create a stable temperature gradient. 13,14 This created a heat flow in the system from the top to the bottom, and by computing the cumulative energy added and subtracted by the two thermostats, the energy flux was also easily obtained. The bottommost layer of the bottom solid wall was fixed, while for the topmost layer of the top solid wall, only the relative positions of the atoms were fixed, and a constant downward force was applied to them. This resulted in the top grooved surface acting as a piston to maintain a steady system pressure in the z direction. The pressure applied to the system was set to 1 MPa, which is in between the saturated pressure of 110 and 120 K for the argon model used in this work, as described in Sec. II B 2, i.e., this produced a slightly compressed liquid. The liquid atoms were initially positioned in an FCC lattice with the lattice constant set to 2 2 3 σ and were given an initial temperature of 60 K. The number of liquid atoms for each system type is shown in Table II. Simulations of constant particle number, constant pressure, and constant cross section were run for 100 ns with the solid-liquid interaction coupling parameter set to η = 1 and reached a stable state with a temperature gradient in the z direction. It was confirmed that the average temperature of the liquid phase was approximately 110 K and there was no gas phase for all systems. The resulting systems are displayed in Fig. 1(b). The solid-liquid interaction coupling parameter η was then set to values in the range of 0.1-0.9 with a step of 0.1, an additional range of 0.12-0.28 with a step of 0.02, and also a η = 0.01 value for very weak solid-liquid interactions. An additional run of 100 ns was conducted for each new system. Snapshots demonstrating the effect of weak and strong solid-liquid interactions are displayed in Fig. 1(c). Afterwards, a run of 150 ns for each system was conducted, from which various properties were calculated. This resulted in 19 η values for each of 6 system types, resulting in 114 unique systems. Only the bottom grooved surface was used for analysis as the top grooved surface was in constant motion. It should be noted that for several pitch 1 systems, the dangling atom structure was disrupted in the top grooved surface due to higher temperature, but this had only negligible effect on the results obtained from the bottom grooved surface.
Because of the need to compute work of solid-liquid adhesion of the grooved surface, equilibrium systems were also created in addition to the heat flow systems. The thermal equilibrium systems were created in an identical manner, by removing the thermostat from the bottom solid wall and setting the control temperature of the top one to 110 K. Both top and bottom grooves were used for analysis as only their energetic properties were needed.

III. RESULTS
A. Grooved surfaces

One-dimensional density and temperature profiles
The density and temperature distribution profiles of all system types for solid-liquid interaction coupling parameter values η = 0.1 and η = 0.5 are displayed in Fig. 2. Linear fittings of temperature profiles inside liquid and solid surface bulk areas are also displayed. The liquid bulk region was set as the area that is further than 25 Å from the mean position of the solid surface innermost layers, while the solid bulk was set as the 4th to 13th solid wall layers, when numbered from the outermost layer. Note that because the one-dimension density distribution of liquid does not take surface bump volume into account, liquid density is essentially halved inside the grooves. None the less, we can observe that at weak solid-liquid interaction of η = 0.1, almost no adsorption layer existed at the FIG. 2. One-dimensional density and temperature distributions at grooved surfaces for weak solid-liquid and strong solid-liquid interactions, where the solid-liquid interaction coupling parameter is set to η = 0.1 and η = 0.5, respectively, and the origin of z is set as the position of the bottommost solid wall layer. Black lines represent the density distribution of liquid. The temperature of each wall layer is displayed as red points and that of liquid is displayed as red lines, where temperature is only shown when liquid density is over 0.01 g/cm 3 . Linear fittings of the temperature distributions are shown as hollow blue circles over blue dashed lines.  35 although the low density indicates the presence of gas phase, as shown in the example in Fig. 1(c). A clear adsorption layer became apparent at stronger solid-liquid interactions with η = 0.5, and all grooves were filled with liquid, i.e., the Wenzel wetting regime, 36 such as displayed in Fig. 1(d).
The effect of different solid-liquid interaction strength can also be clearly observed from the liquid temperature distribution illustrated as red lines in Fig. 2. At η = 0.1, liquid temperature distribution was almost horizontal for all surfaces, indicating very low heat flux inside the liquid. Because the systems were in a stable state, heat flux between the two heat baths was uniform, which indicates that there was also only miniscule heat transfer between grooved surfaces and liquid molecules. On the other hand, at stronger solid-liquid interactions of η = 0.5, a clear gradient in liquid temperature distribution can be observed. This is also reflected in the temperature jump between surface and liquid, where stronger solid-liquid interaction resulted in smaller temperature jumps.
Contrary to the liquid, the temperature profile of the solid surface indicated by red dots in Fig. 2 is mostly horizontal regardless of the η value, indicating a much higher thermal conductivity than the liquid. Only at the solid-liquid interface, the temperature distribution is no longer linear, indicating surface effect. 13 The nonlinear effect was strongest in pitch 1 systems as can be observed from comparing pitch 1 and other systems at η = 0.5 in Fig. 2. The temperature point closest to the solid-liquid interface in the pitch 1 systems represents the temperature of the dangling atoms at the surface, which were described in Sec. II B 3, and this indicates

Two-dimensional density profiles
Because at groove scales larger than pitch 1 it is no longer possible to faithfully represent system states by one-dimensional density profiles, two-dimensional density profiles for each system type at η = 0.1 and η = 0.5 are displayed in Fig. 3. It is possible to determine system wetting regimes by investigating liquid molecule density inside the grooves: a very low density would indicate the Cassie regime, while high density comparable to that of liquid bulk would indicate the Wenzel regime, as illustrated in Figs. 1(c) and 1(d). For a flat surface, we observe that the same trends as seen in Fig. 2, i.e., at low solid-liquid interaction with η = 0.1, almost no adsorption layer existed, while at stronger interaction with η = 0.5, an adsorption layer could be observed, with multiple layers continuing into the bulk.
Density distribution of the pitch 1 system in Fig. 3 at η = 0.5 reveals a very unique interface, where high density patches of liquid were concentrated between the dangling atoms, acting as pillars, at the grooved surface and indicates packing of liquid molecules at the solid surface. These high density patches themselves acted in a similar manner and caused creation of slightly lower density patches above them. The effects can be seen continuing for several iterations, in a similar manner observed for the adsorption layers of the flat surface.
Looking at other systems with weak solid-liquid interactions at η = 0.1, we can observe that none of the grooves were filled. For large scale grooves, at pitch 10 and 20, a clear transition between the Cassie and Wenzel regimes existed at 0.2 < η < 0.24 and 0.18 < η < 0.2, respectively and can be observed from the liquid density inside grooves shown in Fig. 4. This can be explained by the capillary effect, i.e., the position of the solid-liquid interface inside the grooves is decided by the balance between interfacial tension due to interface curvature and the pressure difference between gas and liquid phases. Strictly speaking, current steady state simulations are not enough to determine the presence of capillary effect, but previous numerical simulation studies have observed capillary action in nanochannels and carbon nanotubes that are of comparable scale to larger grooves; 39,40 therefore, we can apply this macroscopic model to our systems. As the pressure difference should be negligible due to lack of gravity and small system size, only interfacial tensions, i.e., contact angles, must be considered. Because weak solid-liquid interactions translate to bad wettability, i.e., contact angle more than 90 ○ , capillary action dictates that liquid should not enter the grooves at all. On the other hand, when good wettability, i.e., contact angle less than 90 ○ , is achieved, the liquid should fill the grooves completely as is observed for all systems at η = 0.5 in Fig. 3. The exact transition happened at larger η for pitch 10 than pitch 20. As a similar relation was confirmed for equilibrium systems, this is not due to temperature difference, but likely due to the groove scale at the solid surface, which would change the transition point from Cassie to Wenzel regimes.
For grooves below pitch 5, there was no clear transition between the Cassie and Wenzel regimes, as indicated by liquid molecule density inside the grooves in Fig. 4. It appeared that a gas phase existed in the grooves, but instead of a clear wetting regime transition, the liquid density inside the grooves gradually increased with the increase in η, indicating that at this scale, the capillary action no-longer occurs. The decrease in density at high η values is because packing of liquid molecules occurred inside the grooves, and the area used to compute liquid density, shown in Fig. 3 at pitch 5, η = 0.1 mainly encompassed low density regions, while that of pitch 10 and 20, also displayed in Fig. 3 at η = 0.1, encompassed both low and high density regions due to a larger size. The packing of liquid molecules will be discussed in more detail in Sec. III B.
It is interesting to note that the current pressure control scheme was selected to prevent the formation of the vapor phase in the system, but it was present none the less inside the grooves during Cassie wetting regimes. Because the temperature of the gas molecules in the groove was no more than 110 K, the saturation vapor pressure of these molecules was below 0.79 MPa, as calculated in Sec. II B 2, which, in turn, was below the control pressure of 1 MPa. We set a control volume at 15 Å ≤ x ≤ 25 Å and 55 Å ≤ z ≤ 65 Å in the center of the groove inside the pitch 20 system at η = 0.1 so that it would be out of direct interaction with the groove walls and solid-liquid interface. In this area, the mean density and temperature were approximately ρ gas = 0.04 g/cm 3 and T gas = 108.7 K, respectively, where the "gas" superscript indicates the gas phase values inside the groove. Because in the gas phase, the greatest contribution to pressure comes from the kinetic pressure term and the virial term can be ignored, 41 we can estimate the gas pressure inside the groove by using the ideal gas law P gas ≈ ρ gas kT gas m liq. ≈ 0.92 MPa, where k is the Boltzmann constant. This is almost the same as the control pressure and clearly above the saturation vapor pressure for the temperature of the gas molecules inside the groove. This is slightly counter-intuitive, but a consequence of the need for a force balance between the liquid and gas phases in the z direction and the lack of curvature at the liquidgas interface. Therefore we conclude that during the Cassie wetting regime, the grooves were filled with over-saturated vapor, and it can The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp be speculated that given a large enough groove scale, condensation of droplets inside the grooves may begin. Finally, we note that for all systems at η = 0.5 with pitches larger than 1, liquid molecules fully entered the grooves and formed high density adsorption layers, indicating high wettability. It is notable that for pitch 3, only a single liquid atom appeared to be able to fit in the horizontal x direction in the groove, resulting in two high density patches aligned vertically. For the pitch 5 system, two liquid molecules appeared to be able to line up at the bottommost groove positions, but appear to have more leeway, which resulted in a more disordered structure distribution above them. Pitch 10 and 20 system grooves appeared to be wide and tall enough to create a uniform adsorption layer at the groove walls, except for the corners. Due to the nature of single-atom molecule liquids, we can see that adsorption layers filled up the grooves even in pitch 20 systems, preventing the formation of liquid bulk in the strict sense.

B. Interfacial thermal conductance
The interfacial thermal conductance (ITC) of an interface can be obtained by where Q is the heat flux and ∆T is the observed temperature jump at the interface. The heat flux was obtained from computing the kinetic energy that was added to the system by the two Langevin heat baths, as described in Sec. II B 4, where E top and E bottom are the total kinetic energy added by the top and bottom heat baths over time t, while A is the xy cross section area. The one-dimensional density and temperature profiles shown in Fig. 2 were used to obtain the temperature jump at the interface. The exact definition of the interface position has substantial arbitrariness for grooved surfaces. In this work, the midpoint between the bottommost point where liquid density is more than 0 and the topmost layer position of the bottom grooved surface was used, which resulted in the excluded volume area in the solid-liquid interface for flat and pitch 1 surfaces, and approximately the middle of the groove for other surfaces with larger pitch values. This interface definition allowed us to obtain consistent results. As was discussed in Sec. III A 1, the temperature profiles in the vicinity of the interface in Fig. 2 are non-linear; therefore, values extrapolated to the solidliquid interface from the linear fits are used instead of the actual values. 13 The obtained ITC is displayed in Fig. 5. The error bars display the standard error of the mean, which was obtained via propagation of error from the standard errors of mean of the temperature jump ∆T and heat flux Q, which were obtained by taking statistical inefficiency into account. 33 We can roughly divide the observed tendencies into several phases based on the solid-liquid coupling parameter. First, for superhydrophobic solid-liquid interactions at 0.01 ≤ η ≤ 0.1, very low ITC is obtained, regardless of the surface groove structure. As can be observed from Figs. 2 and 3 at η = 0.1, a high density adsorption layer is absent, which prevented effective heat transfer between solid and liquid. It has been reported that at low wettability, i.e., weak solid-liquid interaction, flat surfaces have higher ITC compared to grooved ones 19,20 because of more contact surface but is not yet observed at this stage due to the small values of ITC.
The second phase is at approximately 0.1 ≤ η ≤ 0.2, where the solid-liquid interaction is still hydrophobic, but the differences can be observed between different surface types, as is emphasised in the internal panel of Fig. 5. Pitch 10 and pitch 20 systems show a sharp jump in the ITC values at approximately η ≈ 0.2, which corresponds to the transition from hydrophobic (Cassie wetting) to hydrophilic (Wenzel wetting) stages, i.e., liquid entered inside the grooves, improving heat transfer, which can be confirmed from showing that only at this scale, the groove size starts to have an effect. Curiously, even at this stage, flat surfaces did not display higher ITC, unlike what was reported in previous studies, 19,20 but were rather on par with that of pitch 10 and pitch 20 systems. The reason for this is not completely clear but might be explained by the difference in the pressure control scheme, which allows the liquid phase in our systems to move away from an unfavorable solid, thus negating any advantage that flat surfaces might have due to the large solid-liquid contact area.
The third phase, corresponding approximately to 0.25 ≤ η ≤ 0.3, is the hydrophilic interaction phase, where the liquid fully fills the grooves, but the contact angle with the flat surface is more than 0 ○ , i.e., not fully wettable, as will be discussed in Sec. III C. Grooved surfaces all appear to have similar ITC values, which are roughly twice that of the flat surface, and demonstrate that even at very small The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp groove scales, the effect of the extended heat transfer surface can be observed. This is somewhat surprising for the pitch 1 system as looking at the density profiles in Fig. 3, the heat transfer surface does not appear to be substantially increased. It will be demonstrated later in Sec. III C that pitch 1 systems appear to use a different mechanism for heat transfer improvement compared to systems with larger groove scales. Finally, at the phase of superhydrophilic interactions at 0.5 ≤ η, we see deviations of ITC for different system types. First, for the grooved surface systems, pitch 3 appears to have the lowest conductance, followed by pitch 20. Pitch 5 and pitch 10 systems maintain the highest ITC values, approximately twice the amount of flat surfaces. Pitch 1 systems appear more erratic, having ITC values on par with pitch 5 and pitch 10 at η ≤ 0.7 and falling in between pitch 5/10 and pitch 20 values at 0.8 ≤ η. For pitch 1 systems, the behavior stems from the fact that it can not be strictly considered a grooved surface, as described in Sec. II B 3. The dangling atoms on the surface are thought to act similar to chemical groups, attached through chemical modification. 37,38 This effect can be observed as higher temperature of the top surface layer, corresponding to dangling atoms, observed in Fig. 2 for pitch 1 systems at η = 0.5. On the other hand, it is not immediately clear why pitch 20 systems exhibit a lower ITC than pitch 5 and 10 systems. A closeup display of density profiles at surface grooves with very strong solidliquid interaction is shown in Fig. 6. A clear crystallization of liquid molecule in the grooves can be observed. When comparing pitch 20 systems with pitch 5 and pitch 10 systems that have higher ITC, it is apparent that the liquid molecules at the bottom of the groove were less restricted although still being more restricted than that of flat systems. Therefore, it is apparent that even at the largest groove scales of pitch 20, liquid molecules inside the grooves were still constrained by the groove geometry, which affected packing of liquid molecules inside the grooves and changed heat transfer properties. This is further reinforced by the temperature distribution in Fig. 2(l), where the temperature gradient inside the groove is shown to decrease for pitch 20 systems at η = 0.5, whereas this is not the case for pitch 5 and 10 systems, which could indicate that less energy is being transferred. Note that these density profile tendencies still remain even with weaker solid-liquid interaction strength, such as η = 0.3, but this does not appear to affect heat transfer as much.

C. Interfacial thermal conductance and work of adhesion
Previous computational and experimental studies have demonstrated a proportional relation between work of solid-liquid adhesion or solid-liquid interaction strength and interfacial thermal conductance (ITC). 19,[42][43][44] On the other hand, it was also reported by Wang et al. that in a constrained channel, there is an optimal solidliquid affinity that achieves the highest heat conductance. 45 In this section, we investigate these relations in detail for our wide range of groove types and solid-liquid interaction conditions.
The left panel of Fig. 7 shows the relation between the solidliquid coupling coefficient η and mean value of solid-liquid potential energy per unit area ⟨u sl ⟩ for the equilibrium systems, described in Sec. II B 4. Several interesting trends can be observed from this graph, depending on the wettability regime. At hydrophobic interactions at approximately η ≤ 0.2, there is a trend that grooved surfaces appear to have larger ⟨u sl ⟩, than that of flat surfaces, which is mostly due to the difference in the solid-liquid contact area. Pitch 10 and 20 systems have clear jumps due to the change in the wetting regime, as discussed in Sec. III B. For hydrophobic interactions, at approximately 0.3 ≤ η, there is a general trend that systems with larger groove scale have smaller ⟨u sl ⟩, which appears to saturate at pitch 10. Apparently, groove corners and edges effect the solid-liquid potential unfavorably, and pitch 20 systems with least of groove corners and edges have the lowest potential energy. These energetic tendencies do not correlate well with what was observed in Fig. 5, where pitch 5 and pitch 10 systems showed greater ITC. Additionally, the pitch 1 system appears to have higher solid-liquid potential energy than even flat systems, which is surprising as pitch 1 systems clearly exhibited higher ITC, as shown in Fig. 5. In conclusion, a simple relation cannot be drawn between interfacial potential and ITC values.
We use the dry-surface method described in earlier studies 46,47 to obtain the work of solid-liquid adhesion. In brief, the work of solid-liquid adhesion for a system with solid-liquid coupling parameter η ′ can be obtained by the following thermodynamic integration: where the 1 2 coefficient on the right-hand side is to account for two interfaces created by the bottom and top grooved surfaces. In principle, free energy change due to change in system volume, P∆V, should also be considered but was omitted due to its small value.
The work of solid-liquid adhesion for each system at various η values is displayed in the right panel of Fig. 7. Overall, the same tendencies as seen for solid-liquid potential in the left panel of Fig. 7 can be also observed for work of solid-liquid adhesion: systems with greater groove scales have larger W sl values, while pitch 1 shows an anomaly, having lower work of adhesion than even flat systems. If of Chemical Physics we use the Young-Dupré equation, where θ is the contact angle, precise hydrophobic, hydrophilic, and superhydrophilic regions can be determined. We use the contact angle to define hydrophobic, hydrophilic, and superhydrophilic regions via work of adhesion values of flat surfaces, where 90 ○ ≤ θ, 0 ○ ≤ θ ≤ 90 ○ , and θ ≤ 0 ○ , respectively. When using the γ lv value for 110 K for simplicity determined in Sec. II B 2, we obtain that the transition from hydrophobic to hydrophilic state should occur in the 0.3 < η < 0.4 region, while the transition to superhydrophilic state should occur at the 0.5 < η < 0.6 region. This is close to what we stipulated in Sec. III B, but fact that the filling of grooves in pitch 10 and pitch 20 systems occurs at different points and earlier than predicted by the Young-Dupré equation suggest that the capillary action is affected by the groove scale. This is reaffirmed by the fact that smaller pitch systems do not show any capillary action effect. The jumps seen in solid-liquid potential energy are not present in work of adhesion for pitch 10 and pitch 20, but the change in wettability is reflected in the change of line gradients.
The relation between solid-liquid potential energy or work of solid-liquid adhesion and ITC is shown in the left and right panels of Fig. 8, which display similar tendencies, except for small potential energy or work of adhesion values. This similarity is expected  as it has been reported that change in work of solid-liquid adhesion and solid-liquid potential energy is the proportional ratio once the liquid adsorption layer at the solid-liquid interface is saturated. 47 A clear difference gradient in is observed for pitch 1 systems, compared to the rest of grooved and flat surfaces. This again points that pitch 1 systems are quite different from the rest and should be considered as having a different mechanism to increase ITC. On the other hand, although the rest of the systems have somewhat similar gradients, pitch 20 systems have a lower gradient compared to the other grooved surfaces, indicating that this particular groove configuration tends to result in lower ITC, while pitch 5 and pitch 10 systems have a higher gradient, indicating that their groove configurations tend to produce higher ITC. Pitch 3 systems also appear to have a gradient comparable to that of pitch 10, but the largest achievable ITC is smaller because the groove scale limits the amount of solidliquid interaction that can achieved when compared to larger groove scales.
Overall, we have confirmed that increasing liquid affinity does increase ITC, and the ratio is mostly similar regardless of the surface shape, except in the case of pitch 1 systems, which can not be strictly considered a grooved surface and is thought to be closer to a flat surface with attached chemical groups. As ITC appears to be strongly correlated with the solid-liquid potential energy, it can be concluded that as long as the groove width is larger than single liquid molecule diameter, the increase in solid-liquid interaction area results in an almost proportional increase in ITC, and the observed difference between different groove scales was due to difference in solid-liquid potential energies due to geometric constraints, which is also closely related to the difference in packing of liquid molecules at the superhydrophilic wetting regime discussed in Sec. III B.

IV. CONCLUSION
Non-equilibrium molecular dynamic simulations for solidliquid interfaces with grooved surfaces at various interaction strengths under an induced heat flux were conducted to investigate the mechanism of heat transfer at the solid-liquid interface under a wide range of conditions. Both Cassie, i.e., grooves not being filled with liquid, and Wenzel, i.e., grooves filled with liquid, wetting regimes were observed.
At the Cassie, i.e., hydrophobic, wetting regime, large scale grooves had lower interfacial thermal conductance (ITC) than small scale grooves because of the reduced solid-liquid contact area due to the capillary effect. At the Wenzel, i.e., hydrophilic, wetting regime, all grooved surfaces showed higher ITC compared to a flat surface, which was proportional to the increase in the solid-liquid contact area due to grooves.
At superhydrophilic solid-liquid interactions, liquid molecules started to crystallize inside the grooves and groove dimensions influenced packing of liquid molecules even at the largest groove scale of over 10 liquid molecule diameters. As none of the grooved systems could obtain ITC that was clearly more than twice that of a flat surface, the molecular packing only appeared to have an effect of decreasing ITC although the exact tendencies were inconsistent. This hints that for systems with highly ordered liquid molecules at the surface, surface topology should have a significant effect on ITC.
The solid-liquid potential energy and work of solid-liquid adhesion were also investigated, where, except for the smallest scale grooves, ITC was similarly proportional to both solid-liquid potential energy and work of solid-liquid adhesion for hydrophilic and superhydrophilic regimes, regardless of the groove scales. This indicates that the effect the grooves had on ITC was due to the increase in the solid-liquid interaction, i.e., increase heat transfer area, where the differences among the different groove scales was due to different patterns of molecular packing at the superhydrophilic regime, which changed the heat transfer properties of liquid inside grooves. On the other hand, in the limit of the groove scale of single atom, vastly different tendencies were observed, which indicated a different heat transfer mechanism. A proposed explanation was that the semi-dangling surface atoms acted similar to chemical groups attached to the surface.