Unraveling the Regimes of Interfacial Thermal Conductance at a Solid/Liquid Interface

The interfacial thermal conductance at a solid/liquid interface (G) exhibits an exponential-to-linear crossover with increasing solid/liquid interaction strength, previously attributed to the relative strength of solid/liquid to liquid/liquid interactions. Instead, using a simple Lennard-Jones setup, our molecular simulations reveal that this crossover occurs due to the onset of solidification in the interfacial liquid at high solid/liquid interaction strengths. This solidification subsequently influences interfacial energy transport, leading to the crossover in G. We use the overlap between the spectrally decomposed heat fluxes of the interfacial solid and liquid to pinpoint when “solid-like energy transport” within the interfacial liquid emerges. We also propose a novel decomposition of G into (i) the conductance right at the solid/liquid interface and (ii) the conductance of the nanoscale interfacial liquid region. We demonstrate that the rise of solid-like energy transport within the interfacial liquid influences the relative magnitude of these conductances, which in turn dictates when the crossover occurs. Our results can aid engineers in optimizing G at realistic interfaces, critical to designing effective cooling solutions for electronics among other applications.


Finite Size Effects Study
We model our system (Figure 1 of the main article) based on that of Sääskilahti et al., 1 with adjustments to ensure that there are no finite-size effects.Sääskilahti et al. 1 demonstrated that interfacial thermal conductance (G) is not impacted by the fluid channel length (L liquid = 55 Å).However, we choose to increase L liquid to 155 Å so as to ensure confinement effects are avoided, as we test for higher wettabilities.The total thickness of the solid walls L wall is 40 unit cells, with L heat = L cool = 25 unit cells thick.The unit cell constant is 1.56σ .The walls comprise a face-centred cubic (FCC) crystal with its [100] crystallographic plane oriented in the x-direction.A finite-size effects study was conducted to confirm that G is independent of wall thickness beyond 40 unit cells, as seen in Figure S1, for ε SL = 10.3 meV.Sääskilahti et al. 1 found that G is independent of crosssectional dimensions beyond 8×8 unit cells.Nonetheless, we model a larger cross-sectional area of 15×15 unit cells to improve the statistical accuracy of our results.To quantify the relationship between ε SL and wettability when the surface is partially wetted, the contact angles (θ ) of quasi-2D cylindrical droplets are calculated at various magnitudes of ε SL .A portion of the simulation domain implemented to achieve this is shown in Figure S2, where the inset depicts a three-dimensional view of the droplet in the x,y, and z directions.The domain boundaries are periodic except along the y direction, for which we used the 'fixed' boundary condition in LAMMPS.This prevents condensation along the periodic image of the wall that would result at top of the domain.A fictitious reflective wall is placed at the top of the domain using the 'fix wall/reflect' command to prevent the loss of atoms through the fixed boundary.The dimensions of the domain are 300 × 45 × 5 unit cells in the x,y, and z directions, respectively.This ensures the droplets can spread sufficiently without interacting across the x periodic boundary.The wall size in the y direction is set to 3 unit cells, which ensures a thickness larger than the interaction cut-off distance (r cut ).The droplets contain approximately 10,000 atoms.

Contact Angle Simulations
The system is equilibrated under the micro-canonical ensemble at 100 K for 250 ns.The spreading of the droplets is assessed by tracking their center of mass.Following equilibration, the twodimensional density profiles of the droplets are sampled for 50 ns.Subsequently, a circular fit is made through the droplet interface.This is identified by locating the "equimolar" points, where the local density is equal to the average of the bulk liquid and bulk vapour densities. 5The angle that the circular fit through equimolar points makes with the wall is the resulting contact angle θ , as depicted in Figure S2.The resulting values of θ for several magnitudes of ε SL can be seen in Table 1 of the main text.investigated, as shown in Figure S3.However, as Sääskilahti et al. 1 provided data solely for these wettabilities, they did not observe the cross-over in G.
Xue et al. 4 were the first to observe a transition in G from an exponential regime to a linear one for a Lennard-Jones (LJ) system.Modelling alternating slabs of solid and liquid atoms each 10 unit cells wide, for a total of four solid slabs and four liquid slabs, they observed the cross-over threshold to occur at ε SL = ε LL = 10.3 meV.Giri and Hopkins 3 also studied the relation between ε SL and G at an LJ solid-liquid interface.However, they chose to simplify the set-up used by Xue et al., 4 instead modelling a solid wall in the center of the domain encased by two liquid slabs to its left and right, with values of G reported in Figure S3.Giri and Hopkins 3 also mentioned the transition from an exponential regime to a linear one, but did not specify a cross-over threshold.In a subsequent publication, Giri et al. 2 utilised the same all-LJ system, reporting significantly different values of G (see Figure S3).It is unclear where these differences arise from as the parameters are identical in both setups.However, Giri et al. 2 observed the cross-over in G at the threshold of ε SL = ε LL = 10.3 meV, as originally discovered by Xue et al. 4 While an exact cross-over threshold was not reported in their earlier study, 3 it is clear from Figure S3 that this would be very different from their subsequent work.These differences are not explicitly acknowledged or addressed in Giri et al. 2 Comparing our results to Giri et al. 2 in Figure S3, we find that our G values are lower in magnitude when ε SL < 10.3 meV, and higher when ε SL > 10.3 meV.This discrepancy makes our cross-over threshold further to the right at ε SL = 50 meV instead of the threshold of ε SL = 10.3 meV.It is important to note here that Giri et al. 2 equilibrated their system at 170K and 0 MPa, while we equilibrated our systems at 100K and 50 MPa.As a reduction in pressure would produce a decrease and not an increase in G, 1,6 the higher G values obtained by them for ε SL < 10.3 meV cannot be attributed to a different system pressure.
Figure S3: Dependence of interfacial thermal conductance (G) on solid/liquid interaction strength ε SL .Our MD results (50 MPa) are compared to the values of G computed by Giri and Hopkins, 3 Sääskilahti et al., 1 and Giri et al. 2 for similar Lennard-Jones (LJ) systems.Excellent agreement is found with results by Sääskilahti et al. 1 In an attempt to replicate the results by Giri et al., 2 we equilibrate our system at 1 MPa and compute G for the same range of ε SL , as presented in Figure S3.We did not equilibrate the system at 0 MPa because the low magnitudes of G and our fixed heat-flux condition would cause the temperature of the interfacial fluid to exceed its saturation temperature, resulting in phase change.The additional 1 MPa of pressure ensures that phase change does not occur, while still having minimal influence on the magnitude of G. 6 While G values vary due to the reduction in pressure from 50 MPa to 1 MPa, we find that we still cannot reproduce the values of G and the cross-over threshold obtained by Giri et al. 2 Our values of G are still significantly lower than those reported by Giri et al. 2 across the entire range of wettabilities studied.Given our excellent agreement with Sääskilahti et al. 1 and the unusual choice of system pressure by Giri et al., 2 we speculate that the pressure control scheme implemented by Giri et al. 2 could have yielded inconsistent pressures across the different values of ε SL .Additionally, as noted earlier, the LJ material being modelled cannot exist in its liquid phase at the 170K at 0 MPa state 7 specified by Giri et al. 2 As such, without explicitly computing the actual pressure within their system at each wettability, their results cannot be replicated.
Considering the differences in the thermodynamic state with Giri et al. 2 and the discrepancies between their two studies, as well as the difference in geometry compared to the system modelled by Xue et al., 4 we rely on our validation with Sääskilahti et al. 1 for our simulation setup.

Interfacial Liquid Density Layering
The density layering was computed in the interfacial liquid by dividing it into rectangular bins of width 0.125 Å along the x axis (i.e.direction of heat flow, see Figure 1 of the main article) and computing the mass density ρ within each bin.The density layering is shown for select values of ε SL in Figure S4.We observe increasing peak ρ with increasing ε SL , consistent with literature. 8,9e observe that full adhesion of the first liquid layer occurs around ε SL = 30 meV, meaning that liquid density drops to zero in the first valley between adsorbed liquid layers, as annotated in Figure S4.This finding demonstrates that full adhesion of the first liquid layer cannot fully explain the cross-over in G, as it occurs prior to the cross-over threshold of ε SL = 50 meV.Full adhesion of the second liquid layer does not occur even when ε SL = 100 meV (ρ does not fall to zero in the second valley, see the inset of Figure S4), and again cannot be directly linked to the cross-over threshold.Note that the pronounced density layering and large magnitudes of peak ρ are a consequence of our chosen values for ε SL .
Figure S4: The mass density (ρ) of the interfacial liquid computed at select values of the solid/liquid interaction strength ε SL .The inset shows that ρ ̸ = 0 in the valley of the second liquid layer for all ε SL , demonstrating the lack of adhesion of the second layer for all wettabilities.

Interfacial Liquid Radial Distribution Functions
To further investigate if interfacial liquid structuring is linked to the cross-over, the radial distribution function (RDF) was computed within the first, second, and third interfacial liquid layers, as presented in Figure S5.We observe the onset of long-range ordering in layers 1 and 2 around ε SL = 30 meV, as demonstrated by the presence of multiple peaks within 15 Å, as seen in Figure S5.As this occurs prior to the cross-over threshold, we conclude that in-plane structuring cannot fully explain the crossover.No long-range ordering can be seen in layer 3, as demonstrated by the lack of further peaks beyond the first solvation shell across all wettabilities.Note that the significant degree of long-range ordering in layers 1 and 2 is a consequence of our chosen range of values for ε SL .6 Spectral Decomposition of Heat-flux

Methodology
The spectral decomposition of heat-flux (SDHF) is defined as the Fourier Transform of the forcevelocity cross-correlation function (FVCF). 1 To bypass the time-consuming computation of the FVCF, the computational efficiency of the Fourier Transform can be exploited.The cross-correlation theorem states that the Fourier Transform of a cross-correlation function such as the FVCF is equivalent to the dot product between the complex conjugate of the Fourier transform of one time series and the Fourier Transform of the second time series. 10sing the cross-correlation theorem, the SDHF of the interfacial solid can be computed via: where A is the system's cross-sectional area, M is the length of the time series (i.e.number of samples), and ∆t s is the sampling interval.ℜ represents the real part of the resulting dot product.F * i j (ω) is the complex conjugate of the Fourier Transform of F i j (t), where F i j (t) is the time series of the cumulative force on a solid atom i due to all liquid atoms j within the cut-off distance.V i (ω) represents the Fourier Transform of the velocity time series v i (t) of the solid atoms i.
To compute the SDHF of the interfacial liquid, we similarly use: Here, F * i j (ω) remains the complex conjugate of the Fourier Transform of F i j (t).However, F i j (t) is now the time series of the cumulative force on a liquid atom i due to all solid atoms j within the cut-off distance.Similarly, V i (ω) is now the Fourier transform of the velocity time series v i (t) of the liquid atoms i.
We compute the SDHF of both the interfacial solid and liquid by dividing the simulation into sequential time segments.Each time segment comprises 20,000 samples (M = 20000), sampled every 25 timesteps, with a timestep of 0.002 ps (∆t s = 25 × 0.002 = 0.05 ps).This yields time segments each 1 ns long.The SDHF of both media are computed from the independent time segments, and are then averaged to yield the spectra presented.

Validation
To validate our SDHF calculation methodology, we replicate the system simulated by Sääskilahti et al. 1 and compute the SDHF of the interfacial solid.Dividing the SDHF by the interfacial temperature discontinuity ∆T yields the spectral decomposition of the interfacial conductance.This is then plotted against the results presented by Sääskilahti et al. 1 Reasonable agreement can be seen in Figure S6, demonstrating the validity of our approach.Statistical noise in SDHF predictions can be diminished by computing more time windows and further averaging.The vibrational density of states (VDOS) is computed from the Fourier transform of the velocity autocorrelation function (VACF) for a group of atoms, where the time-varying VACF is defined as 11 Here, v i (t) represents the velocity of an atom i at time t, while v i (0) represents the initial velocity.Similar to the computation of the SDHF, the computational efficiency of the Fourier Transform can be made use of to avoid the direct computation of the VACF.The autocorrelation theorem is a specific case of the cross-correlation theorem where the two time series are the same, and states that the Fourier Transform of an autocorrelation function is equivalent to the squared modulus of the Fourier Transform of the original time series. 12Thus, the VDOS can be obtained using: where i is the atom belonging to a group of atoms of size N.
In the case of the LJ solid, we compute the interfacial VDOS from the outermost layer of the solid in contact with the liquid.For the LJ fluid, we compute the interfacial VDOS from the first adhered liquid layer, defined as first peak in the interfacial liquid density layering, as presented in Section 3. To account for the diffusion of the liquid, we only compute the VDOS for atoms that remain within the interfacial region for the entirety of the calculation, discarding atoms that diffuse in/out of this region. 13Similar to the computation of the SDHF, we divide the simulation into sequential time segments.Each segment is 500 samples long (M = 500), sampled at intervals of 25 timesteps.For a timestep of 0.002 ps, this results in time segments each 25 ps long.The resulting spectra are then to yield the VDOS

Validation
To validate our methodology, we compute the VDOS of bulk aluminium at 300K and validate our normalised results against those reported by Korotaev et al. 14 Excellent agreement can be seen in Figure S7, thus validating our VDOS calculations.
Figure S7: The vibrational density of states (VDOS) of bulk aluminium at 300K.Excellent agreement can be seen when comparing to results by Korotaev et al. 14

Figure S1 :
Figure S1: Interfacial thermal conductance (G) computed for the solid/liquid interaction strength ε SL = 10.3 meV at various wall thicknesses.No dependence on wall thickness is found beyond 40 unit cells, indicating the absence of finite size effects.

Figure S2 :
Figure S2: An illustration of the method employed to compute the contact angle (θ ) of quasi-2D cylindrical droplets.The inset depicts a three-dimensional view of the droplet.

Figure S5 :
Figure S5: Radial Distribution Function (RDF) computed for the first, second, and third interfacial liquid layers.

Figure S6 :
FigureS6: The spectral decomposition of G at the interfacial solid, defined as SDHF(ω)/∆T , computed for the solid/liquid interaction strength ε SL = 10.3 meV.Reasonable agreement can be seen when comparing to results by Sääskilahti et al.1 for the same system.