Thermal Transport at Solid–Liquid Interfaces: High ... - ACS Publications

Subscriber access provided by University of Newcastle, Australia

Letter

Thermal Transport at Solid-Liquid Interfaces: High Pressure Facilitates Heat Flow Through Non-Local Liquid Structuring Haoxue Han, Samy Merabia, and Florian Müller-Plathe J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/acs.jpclett.7b00227 • Publication Date (Web): 12 Apr 2017 Downloaded from http://pubs.acs.org on April 13, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry Letters is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 17

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

Thermal Transport at Solid-liquid Interfaces: High Pressure Facilitates Heat Flow Through Non-local Liquid Structuring Haoxue Han,∗,† Samy Mérabia,‡ and Florian Müller-Plathe† †Theoretische Physikalische Chemie,Eduard-Zintl-Institut f¨ ur Anorganische und Physikalische Chemie, Technische Universität Darmstadt. Alarich-Weiss-Straße 8, 64287 Darmstadt, Germany ‡Institut Lumière Matière UMR 5306 CNRS Université Claude Bernard Lyon 1, Bˆ atiment Kastler, 10 rue Ada Byron, 69622 Villeurbanne, France. E-mail: [email protected]

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract The integration of three-dimensional microelectronics is hampered by overheating issues inherent to state-of-the-art integrated circuits. Fundamental understanding of heat transfer across soft-solid interfaces is important to develop efficient heat dissipation capability. At the microscopic scale, the formation of a dense liquid layer at the solid-liquid interface decreases the interfacial heat resistance. We show through molecular dynamics simulations of n-perfluorohexane on a generic wettable surface that enhancement of the liquid structure beyond a single adsorbed layer drastically enhances interfacial heat conductance. Pressure is used to control the extent of the liquid layer structure. The interfacial thermal conductance increases with pressure values up to 16.2 MPa at room temperature. Furthermore, it is shown that liquid structuring enhances the heat transfer rate of high energy lattice waves by broadening the transmission peaks in the heat flux spectrum. Our results show that pressure is an important external parameter that may be used to control interfacial heat conductance at solid-soft interfaces.

Graphical TOC Entry

b

High pressure facilitates heat flow through non-local liquid structuring: insight beyond the first liquid layer


Page 2 of 17

Page 3 of 17

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Understanding nanoscale heat transfer across hard-soft material interfaces, e.g. solidpolymer and solid-liquid interfaces, is of great importance for improving the reliability of nanoelectronic systems where efficient heat dissipation becomes critical. 1–6 Key quantities that thermally characterize a given dissimilar interface are the interfacial thermal conductance G = q/∆Tsl , with q the heat flux and ∆Tsl the temperature difference at the interface, and the associated thermal heat (Kapitza) resistance R = 1/G. Over the past decade, considerable efforts have been devoted to directly probe G or R through both experiments 7–11 and molecular simulations. 12–24 Research groups have attempted to find a correlation between G and other quantities that characterize solid-liquid interfaces. Shenogina et al. reported a correlation between G for interfaces between water and hydrophilic and hydrophobic substrates and the work of adhesion Wadh of these interfaces. 13 Ramos-Alvarado recently discussed the limits of this result and showed that two interfaces with the same Wadh may have different G. 25 Alexeev 26 and Ramos-Alvarado 25 reported correlations between G for water on nonpolar surfaces and quantities that measure the density of the first adsorbed water layer. Hung 24 found a temperature dependence of the density of the first water layer affecting G between water and self-assembled monolayers on a gold surface. It is thus clear that the formation of adsorbed layers with enhanced density is crucial to increase interfacial heat conduction. In this Letter, we report nonequilibrium molecular dynamics (MD) simulations of nperfluorohexane (n-C6 F14 ) on a non polar-substrate at room temperature to determine the pressure (p) dependence of G. n-C6 F14 as a stable fluorocarbon-based fluid is widely used as electronics coolant liquids. 27 The reason lies in their relatively low boiling point and high evaporation latent heat. We found that increasing pressure yields a threefold to fourfold enhancement of interfacial heat conduction. This drastic enhancement is to a large extent explained by pressure induced liquid structuring which extends a few nanometers away from the interface. Our simulation results shed light on the mechanisms of interfacial heat transfer at solid-liquid interfaces, and show the necessity to consider non-local effects when the pres-



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

sure applied is larger than a few MPa. Our conclusions are different than previous simulation results which conclude that G linearly increases with the density of the first liquid layer, but do not consider the liquid layering. 25,26 We show here that a model based on the density of the first layer can only describe the pressure dependence of G at low pressure values, but underestimates the increase of G by a factor four at large pressure values. In the MD simulations reported here the all-atom model of Watkins et al. 28 was used for n-C6 F14 . The solid surface was a face-centered cubic (FCC) lattice with the lattice parameter of gold. The interaction between the gold atoms was modeled with the 12-6 Lennard-Jones (LJ) pair potential with the parameters of Heinz et al. 29 Note that this model was found to reproduce satisfactorily experimental mechanical properties such as the bulk, shear and Young’s moduli. The interaction between liquid and solid was assumed to be dominated by van der Waals interactions and modeled through the 12-6 LJ potential. The LorentzBerthelot arithmetic mixing rule was used for the distance parameters, i.e., between the gold and carbon atoms and between the gold and fluorine atoms. The energy parameters were obtained through the geometric mixing rule. The OPLS parameters were used for carbon and fluorine, while different parameters were employed for gold depending on the crystallographic plane in contact with n-C6 F14 . When n-C6 F14 interacted with the (111) crystallographic plane a value of 0.381 kJ/mol was used whereas 0.342 kJ/mol was used with the (100) plane. These two values were chosen such that both the (111) and the (100) surfaces yield a value for the wetting contact angle of 3◦ , i.e., the system has nearly fully wetting behavior. The energy parameters for gold reported above were determined through the dry-surface method to be consistent with this contact angle 30 (See section I of Supporting information). A nonequilibrium setup was used to study the heat flow across the solid-n-C6 F14 interface. Two solid-n-C6 F14 interfaces were created by coupling the n-C6 F14 liquid to the solid surfaces, as shown in Fig. 1 (See section II of Supporting information for details). Ten monolayers of solid atoms with a thickness of Lbath = 5.4 nm next to the frozen atoms at the left


Page 4 of 17

6000 5000 (a) 4000 3000 2000 1000 0 0 5 350 340 330 320 310 300 290 280 270 260

Jin -Jout 10

t (ns)

15

1.0 0.8 f(¯v) ∝ exp( − m¯v¯2 ) 2kB T 0.6 0.4 0.2 (b) 0.0 20 -2 -1 0 1 2

f(¯v) × 10 3 (s/km)

Ecum (kCal/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


T (K)

Page 5 of 17

v¯ (km/s)

Jout

Jin

(c) 50

100

150

x (Å)

200

250

Figure 1: (a) Cumulative input and output energy Ecum as a function of time t in the hot and cold heat reservoirs. (b) Histogram of the center-of-mass velocities v¯ of n-C6 F14 molecules at the center of the system. The black solid curve refers to a normalized gaussian function √ m¯ ¯ v2 f (¯ v ) = m/2πk ¯ B T exp(− 2kB T ). (c) Temperature distribution T (x) along the longitudinal direction ⃗ix of the confined system. The errorbars correspond to the standard deviations of temperatures of each bin. The red line is a linear fit to the temperature points in the n-C6 F14 liquid. The temperature of the hot and cold reservoirs is kept constant at 310 K and 270 K, respectively. The temperature drop ∆Tsl at the interface is estimated by extrapolating the linear fits to the temperature profiles in the liquid and calculating the difference at the interface. and right ends were coupled to Langevin heat baths 31 at the temperatures Th = T +∆T /2 and Tc = T −∆T /2, where T = 290 K and ∆T = 40 K. The bath time constant was tb = 1 ps, ensuring that the bath-induced mean free path of phonons satisfies Λb = vg tb ≤ 2Lb (vg = 3400 m/s is the longitudinal speed of sound in gold 32 ) so that phonons arriving at thermostats are dissipated before reflecting from the fixed ends (see section II of Supporting Information). The solid surfaces had a cross section of 3.7×3.7 nm2 and a thickness of 7.5 nm. The liquid film had a thickness of 8 nm at 1 bar. Additional simulations were carried out with a 3.2 nm-thick liquid film. The cumulative input and output energy Ecum (t) in the hot and cold heat reservoirs is shown in Fig.1(a). The linear profile demonstrates the steady-state heat flow in the system. The established steady-state heat flux is calculated as q = 1/A·∂Ecum /∂t where A is the surface area. A typical temperature profile obtained from the nonequilibrium simulation is shown in Fig. 1(c). The center-of-mass (COM) velocities v¯ (components in



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 17

the three axes of Cartesian coordinates) of n-C6 F14 molecules inside the bin at the center of the system are recorded and averaged for 500 ps. The histogram of the velocities are shown in Fig. 1(b). It can be observed that the COM velocity follows the Maxwell-Boltzmann √

distribution with the probability density of f (¯ v) =

2

m¯ ¯v m/2πk ¯ ¯ is the B T exp(− 2kB T ), where m

mass of a n-C6 F14 molecule, kB the Boltzmann constant and T the temperature. The temperature T corresponding to the velocity distribution is determined to be 289 K which agrees relatively well with the local temperature of the central bin. The temperature drop ∆Tsl at the interface is estimated by extrapolating the linear fits to the temperature profiles in the liquid and calculating the difference at the interface. We note that the measured temperature jump can be affected by the definition of the location of the solid-liquid interface. Since the temperature in the solid is nearly constant here, a solid-liquid boundary defined at the inner most solid layer in contact with liquid would result in a smaller temperature drop compared to a boundary defined at the first peak of the liquid density, and hence a lower Kapitza resistance would be obtained. Therefore, the solid-liquid boundary is chosen to be at the mid-point of the inner most solid layer (in contact with liquid) and first peak liquid density(See Supporting Information). The interfacial thermal conductance is obtained as G = q/∆Tsl , where ∆Tsl is the average of the temperature drop for the two solid-liquid interfaces. Different thermostat temperatures were applied to the systems and a linear dependence was found between q and ∆Tsl , indicating that the system is under investigation was in the linear response regime. Area and length effects were investigated and the dimensions used were those unaffected by the size of the systems. The pressure of liquid was measured on the fly during the non-equilibrium MD simulations and it is confirmed that the three components of the pressure tensor pxx (x), pyy (x), and pzz (x) are isotropic in the bulk region of the fluid in non-equilibrium state (See Supporting Information). The pressure values are calculated by using the arithmetic mean of the three main pressure components of the liquid away from the interface in non-equilibrium simulations p = (pxx + pyy + pzz )/3. 35 The pressure dependence of G up to 16.2 MPa is shown in Fig. 2(a) for the (100) and


35

(100) (1 1 1)

25 quasi-crystalline

20 15 10 5

liquid 0

5

10

p (MPa)

ρ(x) (g cm −3 )

(a)

30

G(p) (MW m −2 K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ρ(x) (g cm −3 )

Page 7 of 17

15

20

3.5 solid-like 3.0 (b) structuring 2.5 2.0 1.5 liquid 1.0 0.5 interface bulk 0.0 0 10 20 30 40 3.0 16.2 MPa 2.5 2.0 12.7 MPa 1.5 1.0 0.28 MPa 0.5 (c) 0.0 0 2 4 x (Å)

Figure 2: (a) Pressure dependence of G(p) for n-C6 F14 on the (100) and (111) surfaces. The error bars correspond to the standard deviation of five independent trajectories. (b) Liquid density distribution ρ(x) perpendicular to the interface on the (100) surface for p = 0.28 MPa, 12.7 MPa and 16.2 MPa. The quantity x refers to the distance from the surface. The region of space noted ”bulk” (blue color) is defined from the simulations at low pressure. It denotes the region where the average liquid density has its bulk value. The region of space noted ”interface” (pink color) denotes the region where density oscillations (liquid layers) are observed. The liquid films have thickness ≈ 10 nm. (c) First maximum of the spatial density distribution for the pressure values 2.6 MPa, 19.9 MPa and 81.8 MPa. The color codes in (b) and (c) are identical. (111) surfaces. It can be observed that G increases with p. An approximate fourfold increase of G is obtained in the considered pressure interval. It can be seen in Fig. 2(b) that n-C6 F14 exhibits a layer structure in the vicinity of the surfaces similarly to water 26 and other molecules like n-hexane. 33 We note that the microscopic mechanism of the formation of the liquid layers in the vicinity of the substrate is the decrease of free energy. Indeed, to minimize its free energy, the system places as many molecules as possible in the range of interaction of the solid-liquid interfacial potential. Molecules in these layers are more ordered than in the bulk of liquid and thus present structural features that tend to resemble a crystalline phase. Moreover, the magnitude of the density maxima but also the extent of the layer structure increase when p increases. We note that the slight displacement of the first liquid layer towards the solid and the emergence of the two “sub-peaks” in the first liquid layer at a high liquid pressure can also enhance the thermal conductance at the interface.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

These colorful features should originate from the complex molecular structure of n-C6 F14 , i.e. longer molecular backbone and different dihedral angles between the perfluorinated blocks, since they were not previously observed in studies on water-solid interfaces under high pressures. 25,26 These phenomena could be an important outlook of our current study and may inspire fellow researchers in related field to conduct further studies to explore the thermal transport at the interfaces between solid and complex liquids. It is also interesting to note that the (100) and the (111) surfaces lead to different values of G, despite the fact that both correspond to the same wetting angle under ambient conditions. This observation is reminiscent of water on the (111) and (100) silicon systems. 25 Qualitative understanding of that difference can be gained when comparing the group velocities of the longitudinal and transverse phonons of gold in different directions at the Brillouin zone center. In fact G depends on the phonon transmission coefficients between solid and liquid and these coefficients depend on the phonons group velocities. The group velocities of the longitudinal and transverse phonons in the Γ-X (100) direction at the zone center are obtained from lattice dynamics of gold 32 as vL100 ≈ 3400 m/s and vT100 ≈ 1500 m/s and in the Γ-L (111) direction vL111 ≈ 3900 m/s and vT111 ≈ 1400 m/s. Therefore, the (111) surface has a higher transmission coefficient than the (100) surface thus a higher thermal conductance, due to a higher average sound velocity. In addition to the difference of phonon group velocities at (100) and (111) surfaces, the larger number of interacting atom pairs between liquid molecules and the (111) surface compared to the (100) surface is also a key factor in the resultant higher interfacial thermal conductance of the former surface. As mentioned above, both G and the density of the first liquid layer ρmax increase with pressure (Fig. 2). In the case of water on non-polar surfaces, Alexeev et al. found a simple one-to-one relationship between ρmax and G. 26 To understand whether this correlation applies for the present systems, we define the reduced liquid density peak ρr = ρmax /ρb , where ρmax is the initial density peak close to solid, as shown in Fig. 2(c), and ρb is the bulk density far away from the interface. In our simulation, ρb ≈ 1.7 kg/l at 1 MPa, in good agreement


Page 8 of 17

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5.0 4.5

l ≈ 8 nm l ≈ 3.2 nm G r ∝ ρr

4.0

Gr = G(ρ)/G0

Page 9 of 17

3.5

16.2MPa 14.3MPa

vdW AMM

3.0

13.4MPa

2.5

6.4MPa

2.0

9.7MPa

2.5MPa

1.5

0.28MPa

1.0 0.5

1.2

1.3

1.4

1.5

1.6

1.7

1.8

ρr = ρmax /ρb Figure 3: Reduced thermal conductance Gr = G(ρ)/G0 of n-C6 F14 against the reduced liquid density ρr at different pressure values for the liquid films with thickness l ≈ 8 nm and 3.2 nm on the (100) surface. The grey dashed line refers to the linear correlation between Gr and ρr . The black solid curve refers to the prediction of the van-der-Waals AMM. with experimental result. 34 In Fig. 3, we plot the reduced interfacial thermal conductance Gr = G(ρ)/G0 of n-C6 F14 against the reduced liquid density ρr at different p for liquid films of thickness l ≈ 8 nm and 3.2 nm on the (100) surface. It is clearly observed that the correlation between the thermal conductance and ρmax goes beyond a simple linear relationship. This means that the higher density of the first liquid layer cannot fully account for the enhancement of the interfacial thermal conductance. It can be seen in Fig. 2 that the enhancement of p also leads to a change in the liquid density far from the interface (bulk). To possibly account for the changes in both the material bulk densities ρb and in the number of molecule collision pairs (density of the first liquid layer ρmax ), we use the vdW acoustic mismatch model 36 (vdW AMM) . In the vdW AMM, the change of ρb affects the sound velocity and therefore the mismatch of acoustic impedances of the two materials (See section III of Supporting information for calculation details). As can be seen in Fig.5, the vdW AMM provides a good qualitative trend of thermal conductance at low pressures < 10 MPa. This is in accordance with the work of Alexeev et al., 26 who correlated the thermal resistance with the density of the first liquid layer at ambient pressure. Nevertheless, the vdW AMM fails to predict the correlation of thermal bounday conductance at high pressures. Therefore,



we show that neither a simple density relation or a more complicated acoustic mismatch model can provide a quantitative prediction of the correlation of the Kapitza conductance for complex liquid under a high pressure condition. However, until now we have not discuss the role of liquid structuring in the interface vicinity on the Kapitza conductance. 10

D(ω) × 1e 3 (THz −1 )

Au C 6 F 14 , 0.28 MPa C 6 F 14 , 16.2 MPa

8 6

(a)

4 2

0

D(ω) × 1e 3 (THz −1 )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 17

5 4 3 2 1 0

0

5

10

(b) 0.28 MPa

0

1

2

15

ω/2π (THz)

3

ω/2π (THz)

4

5

5 4 3 2 1 0

20 (c) 16.2 MPa

0

1

2

3

ω/2π (THz)

4

5

Figure 4: (a) Vibration spectra D(ω) of the solid and liquid at 0.28 and 16.2 MPa for the (100) surface. (b) Vibration spectrum of n-C6 F14 on the (100) surface at 0.28 MPa. (c) Vibration spectrum of n-C6 F14 on the (100) surface at 16.2 MPa. We now discuss the enhancement of G in the light of the vibrational properties of the interface. We show in Fig. 4a the vibrational density of states (VDOS) D(ω) for the (100) surface and the adsorbed liquid film at 3.7 and 72.9 MPa. The spectra are calculated from the Fourier transform of the mass-weighted velocity auto-correlation functions D(ω) = 1 kB T

∑ i

mi

∫∞

−∞

dτ eiωτ ⟨vi (τ ) · vi (0)⟩, with mi the mass of atom i and vi its velocity. For

n-C6 F14 , a few broad peaks are present in the spectra below the Debye frequency of the solid (ωD ≈ 4.5 THz for gold). For ω > ωD , the numerous sharp peaks correspond to the intramolecular eigenmodes and are highly force-field dependent. The spectrum for the solid shows that vibrations above the Debye frequency ωD do not contribute to the interfacial heat flow, although surface modes could be present at the solid-liquid interface. 19 The effect of pressure on the VDOS of n-C6 F14 is shown in Figs. 4b and 4c where it can be seen that the bands broaden due to increased pressure. We show below that the corresponding lattice waves contribute most of the interfacial heat conduction and its pressure dependence.


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


q(ω) [MWm −2 THz)]

Page 11 of 17

500 400 300 200 100 5000 400 300 200 100 5000 400 300 200 100 0

(a) l ≈ 8 nm

25

p =0.28 MPa p =12.7 MPa p =16.2 MPa

(c) l ≈ 3.2 nm

(b) Gcum (ω) [MW m −2 K)]

20

p =0.46 MPa p =6.9 MPa

15

10 (d) l ≈ 8 nm, 16.2 MPa

0

1

2

1 st +2 nd layer 2 nd layer

3

ω/2π (THz)

4

5

5

0

0

1

2

3 4 5

ω/2π (THz)

Figure 5: (a) Frequency-dependent heat flux q(ω) from the (100) surface to the liquid at 0.28, 12.7 and 16.2 MPa for the liquid film of thickness l ≈ 8 nm. (b) Cumulated thermal conductance for the same system as in (a) at 0.28, 12.7 and 16.2 MPa. (c) q(ω) for p = 0.46 and 6.9 MPa for the film with l ≈ 3.2 nm. (d) Propagation of heat flux through the first and second liquid layers (red) compared with the propagation through only the second layer (black) bypassing the first layer for p = 16.2 MPa and l ≈ 8 nm. The frequency-dependent heat flux was calculated following: ∑ ∫ ∞ 2 q(ω) = ℜ dτ ⟨Fij (τ ) · vi (0)⟩eiωτ A j∈l,i∈s −∞

where τ is the correlation time between the velocity vi of the solid atom i, the force Fij from the liquid atom j and ℜ denotes the real part (See section IV of Supporting information). The spectral heat flux from solid to liquid q(ω) delivers information about which vibration modes most contribute G. The variation of q(ω) for different pressure values is shown in Fig. 5(a). At either low pressure (p = 3.7 MPa), two major peak are present in q(ω) around 0.65 and 1.6 THz, indicating their dominant contribution to the interfacial heat transfer. At the intermediate pressure of 47.8 MPa, two additional peaks emerge in q(ω) around 0.98 and 1.33 THz, enabled by a higher pressure. As the pressure continues to increase, the peak around 1.33 THz markedly broadens and the heat flux at higher frequencies ω/2π > 1.5 THz is enhanced. This finding is in contrast to a previous analysis which stated that low frequency phonons are the most significant modes for heat transfer in a quartz-water system 11 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

at ambient pressure by using VDOS. 15 At the highest pressure considered p = 72.9 MPa, the mid-high frequency phonons with 1 < ω/2π < 2.8 THz contributes most to the heat flux, as demonstrated by the cumulated thermal conductance in Fig. 5(b). As discussed above, p enhances interfacial liquid structuring. Such a liquid structuring strengthens the coupling of collective vibrational modes of the liquid molecules with the solid lattice waves (phonons), hence facilitating the phonon transmission from the solid into the liquid. It was shown that for a monoatomic LJ liquid, the heat flux increased with pressure over the whole spectrum without any shift of phonon spectrum, 19 whereas for the molecular liquid in this work, the mid-frequency vibrations above 1 THz and below ωD are more sensitive to the pressure. This phenomenon is confirmed by the accelerated heat transfer rate for frequencies 0.8 < ω/2π < 2.3 THz by pressure for a liquid film of 3.2 nm thick in Fig. 5(c). Our study suggests that not only the first adsorbed liquid layer is the source the enhancement of G. This finding is supported by the observation that the heat flux directly through the second liquid layer bypassing the first layer is approximately 25% of the total flux (Fig. 5(d)). Therefore, the ordered molecules in the first liquid layer in contact with the substrate are not solely responsible for the enhancement of the solid-liquid interfacial heat transfer. In summary, the pressure-mediated thermal transport between liquid n-C6 F14 and a solid surface was investigated by MD simulations. We observe a drastic enhancement of G by pressure. The enhanced liquid structuring arising from the exerted pressure is the micropscopic mechanism underlying the heat transfer enhancement. We identify, through phonon spectral analysis, the contribution of different vibrational modes to the thermal transport across the solid-liquid interface. The heat transfer rate of high-frequency phonons between 1 THz and 4.5 THz was enhanced with higher pressures by broadening the transmission peaks whereas the contribution of low energy vibrations remain unchanged. Such a pressure-induced change in the energy transfer landscape originates from the increased density of high frequency vibrations due to the enhanced liquid structuring. Our results provide fundamental insight of thermal energy transport mechanisms in solid-liquid systems and enable detailed inves-


Page 12 of 17

Page 13 of 17

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tigations of energy transfer between, e.g., water and organic molecules. Previous studies have shown that surface functionalization has emerged as an effective strategy to increase interfacial heat transfer. 22,37–39 The present work strongly suggests that pressure may also be an external parameter that can be used to enhance liquid structuring beyond a few liquid layers and produce an extended solid-like structure with its inherent conductive properties. Acknowledgements- This work was supported by the German Research Foundation (DFG) within the framework of the collaborative research center SFB-TRR75. We thank Kimmo Säa¨skilahti and Frederic Leroy for helpful discussions. We thank Alireza Gholijani, Christiane Schlawitschek, Tatiana Gambaryan-Roisman and Peter Stephan of the Chair of Technical Thermodynamics in TU Darmstadt for providing experimental contact angle. We thank the Lichtenberg-High Performance Computer at TU Darmstadt for the computing resources. Supporting Information Available: Description of the material included.

References (1) Murad, S.; Puri, I.K., Molecular simulations of thermal transport across interfaces: solidvapour and solid-solid. Molecular Simulation, 2012, 38, 642-652. (2) Marable, D. C.; Seungha, S.; Y. Nobakht. Investigation into the microscopic mechanisms influencing convective heat transfer of water flow in graphene nanochannels. Int. J. Heat and Mass Transfer 2017, 109, 28-39. (3) Vo, T. Q.; Kim, B. Transport Phenomena of Water in Molecular Fluidic Channels. Sci. Rep. 2016, 6, 33881 (4) Saleman, B.; Rafeq, A.; Chilukoti, H. K.; Kikugawa, G.; Shibahara, M.; Ohara. T. A molecular dynamics study on the thermal transport properties and the structure of the solid-liquid interfaces between face centered cubic (FCC) crystal planes of gold in contact with linear alkane liquids. Int. J. Heat and Mass Transfer 2017, 105, 168-179.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(5) Chilukoti, K. H.; Kikugawa. G.; Shibahara, M.; Ohara. T. Local thermal transport of liquid alkanes in the vicinity of quartz solid surfaces and thermal resistance over the interfaces: A molecular dynamics study. Phys. Rev. E 2015, 91, 052404. (6) Vo, T.Q.; Kim, B. Interface thermal resistance between liquid water and various metallic surfaces. International J Precision Eng. and Manufacturing 2015, 16, 1341-1346. (7) Harikrishna, H.; Ducker, W.A.; Huxtable, S.T. The influence of interface bonding on thermal transport through solid-liquid interfaces. App. Phys. Lett. 2013, 102, 251606. (8) Ramiere, A.; Volz, S.; Amrit, J. Thermal resistance at a solid/superfluid helium interface. Nat. Mater. 2016, 15, 512-516. (9) Ge, Z.B.; Cahill, D.G.; Braun, P.V. Thermal conductance of hydrophilic and hydrophobic interfaces. Phys. Rev. Lett. 2006, 96, 186101. (10) Park, P.; Cahill, D.G. Plasmonic Sensing of Heat Transport at Solid-Liquid Interfaces. J. Phys. Chem. C 2016, 120, 2814-2821. (11) Zhang, T.; Gans-Forrest, A.R.; Lee, E.; Zhang, X.; Qu, C.; Pang, Y.; Sun, F.; Luo, T., Role of Hydrogen Bonds in Thermal Transport across Hard/Soft Material Interfaces. ACS App. Mater. & Interf. 2016, 8, 33326-33334. (12) Ge, S.; Chen, M. Vibrational Coupling and Kapitza Resistance at a Solid-Liquid Interface. Int. J. Thermophys. 2013, 34, 64-77 (13) Shenogina, N.; Godawat, R.; Keblinski, P.; Garde, S. How Wetting and Adhesion Affect Thermal Conductance of a Range of Hydrophobic to Hydrophilic Aqueous Interfaces. Phys. Rev. Lett. 2009, 102, 156101. (14) Barrat J.L.; Chiaruttini, F. Kapitza resistance at the liquid-solid interface. Mol. Phys. 2003, 101, 1605-1610.


Page 14 of 17

Page 15 of 17

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(15) Goicochea, J. V.; Hu, M.; Michel, B.; Poulikakos, D. Surface Functionalization Mechanisms of Enhancing Heat Transfer at Solid-Liquid Interfaces. J. Heat Transfer 2011, 133, 082401. (16) Mérabia, S.; Shenogin, S.; Joly, L.; Keblinski, P.; Barrat, J.L. Heat transfer from nanoparticles: A corresponding state analysis. Proc. Natl. Acad. Sci. U.S.A. 2009 106, 15113-15118. (17) Murad, S.; Puri, I. K. Thermal Transport Across Nanoscale Solid-Fluid Interfaces. Appl. Phys. Lett. 2008, 92, 133105. (18) Pham, A.; Barisik, M.; Kim, B. Pressure dependence of Kapitza resistance at gold/water and silicon/water interfaces. J. Chem. Phys. 2013, 139, 244702. (19) Sääskilahti, K.; Oksanen, J.; Tulkki, J.; Volz, S. Spectral mapping of heat transfer mechanisms at liquid-solid interfaces. Phys. Rev. E 2016, 93, 052141. (20) Hu, M.; Keblinski, P.; Schelling, P.K. Kapitza conductance of siliconamorphous polyethylene interfaces by molecular dynamics simulations Phys. Rev. B 2009, 79, 104305. (21) Hu, L.; Desai, T.; Keblinski. P. Determination of interfacial thermal resistance at the nanoscale. Phys. Rev. B 2011, 83, 195423. (22) Lin, S.; Buehler, M.J. The effect of non-covalent functionalization on the thermal conductance of graphene/organic interfaces Nanotechnol. 2013, 24, 165702. (23) Kikugawa, G.; Ohara, T.; Kawaguchi, T.; Kinefuchi, I.; Matsumoto, Y. A Molecular Dynamics Study on Heat Transfer Characteristics Over the Interface of Self-Assembled Monolayer and Water Solvent J. Heat. Trans. 2014, 136, 102401. (24) Hung, S.W.; Kikugawa, G.; Shiomi, J., Mechanism of Temperature Dependent Thermal Transport across the Interface between Self-Assembled Monolayer and Water. J. Phys. Chem. C 2016, 120, 26678-26685. 15 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(25) Ramos-Alvarado, B.; Kumar, S.; Peterson, G.P. Solid-liquid Thermal Transport and its Relationship with Wettability and the Interfacial Liquid Structure. J. Phys. Chem. Lett. 2016, 7, 3497-3501 (26) Alexeev, D.; Chen, J.; Walther, J.H.; Giapis, K.P.; Angelikopoulos, P.; Koumoutsakos, P. Kapitza resistance between few-layer graphene and water: Liquid layering effects. Nano Lett. 2015, 15, 5744-5749. (27) “3M Manufacturing and Industry: 3M Fluorinert Electronic Liquids”. 090519 products3.3m.com (28) Watkins, E.K.; Jorgensen, W.L. Perfluoroalkanes: Conformational analysis and liquidstate properties from ab initio and Monte Carlo calculations. J. Phys. Chem. A 2001, 105, 4118-4125. (29) Heinz, H.; Vaia, R.A.; Farmer, B.L.; Naik, R.R. Accurate simulation of surfaces and interfaces of face-centered cubic metals using 12-6 and 9-6 Lennard-Jones potentials. J. Phys. Chem. C 2008, 112, 17281-17290. (30) Leroy, F.; M¨ uller-Plathe, F. Dry-Surface Simulation Method for the Determination of the Work of Adhesion of SolidLiquid Interfaces. Langmuir 2015, 31, 8335-8345. (31) Schneider, T.; Stoll, E. Molcuelar dynamics study of a 3-dimensional one-component model for distortive phase transitions. Phys. Rev. B 1978, 17, 1302-1322. (32) J. W. Lynn, H. G. Smith, and R. M. Nicklow, Phys. Rev. B 1973, 8, 3493 (33) Ardham, V.R.; Deichmann, G.; van der Vegt, N.F.A.; Leroy, F. Solid-liquid work of adhesion of coarse-grained models of n-hexane on graphene layers derived from the conditional reversible work method. J. Chem. Phys. 2015, 143, 243135.


Page 16 of 17

Page 17 of 17

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(34) Pineiro, M.M.; Bessieres, D.; Gacio, J.M.; Saint-Guirons, H.; Legido, J.L. Determination of high-pressure liquid density for n-perfluorohexane and n-perfluorononane. Fluid phase equilibria, 2004, 220, 125-134. (35) Fujiwara, K.; Shibahara, M. Local pressure components and interfacial tension at a liquid-solid interface obtained by the perturbative method in the Lennard-Jones system. J. Chem. Phys. 2014, 141, 034707. (36) Prasher, R. Acoustic mismatch model for thermal contact resistance of van der Waals contacts. Appl. Phys. Lett. 2009, 94, 041905. (37) Han, H.; Zhang, Y.; Wang, N.; Samani, M.K.; Ni, Y.; Mijbil, Z.Y.; Edwards, M.; Xiong, S.; Sskilahti, K.; Murugesan, M; et al. Functionalization mediates heat transport in graphene nanoflakes. Nat. Commun. 2016, 7, 11281. (38) Zhang, L.; Bai, Z.; Liu, L. Exceptional Thermal Conductance across Hydrogen-Bonded Graphene/Polymer Interfaces. Adv. Mat. Interf. 2016, 3, 1600211. (39) Wang, Y.; Zhan, H.F.; Xiang, Y.; Yang, C.;Wang, C.M.; Zhang, Y.Y. Effect of Covalent Functionalization on Thermal Transport across Graphene-Polymer Interfaces J. Phys. Chem. C 2015, 119, 12731-12738.


Thermal Transport at Solid–Liquid Interfaces: High ... - ACS Publications

Recommend Documents