Nanoscale Quantum Thermal Conductance at Water Interface: Green's Function Approach Based on One-Dimensional Phonon Model.

We have derived the fundamental formula of phonon transport in water for the evaluation of quantum thermal conductance by using a one-dimensional phonon model based on the nonequilibrium Green's function method. In our model, phonons are excited as quantum waves from the left or right reservoir and propagate from left to right of H 2 O layer or vice versa. We have assumed these reservoirs as being of periodic structures, whereas we can also model the H 2 O sandwiched between these reservoirs as having aperiodic structures of liquid containing N water molecules. We have extracted the dispersion curves from the experimental absorption spectra of the OH stretching and intermolecular modes of water molecules, and calculated phonon transmission function and quantum thermal conductance. In addition, we have simplified the formulation of the transmission function by employing a case of one water molecule (N=1). From this calculation, we have obtained the characteristic that the transmission probability is almost unity at the frequency bands of acoustic and optical modes, and the transmission probability vanishes by the phonon attenuation reflecting the quantum tunnel effect outside the bands of these two modes. The classical limit of the thermal conductance calculated by our formula agreed with the literature value (order of 10 - 10 W/K) in high temperature regime (>300 K). The present approach is powerful enough to be applicable to molecular systems containing proteins as well, and to evaluate their thermal conductive characteristics.


Introduction
In recent years, advances in the fabrication and characterization of nanoscale systems allow for a better understanding of the heat flow at the microscopic level [1]. Regarding the researches for inorganic materials, studies on carbon nanotube systems have provided a prototypical example for nanoscale thermal conduction experimentally [2][3][4][5][6] and theoretically [7,8]. Both theoretical and experimental progresses have thus been reported extensively in each of these studies [1]. As for other nanoscale systems, Lervik et al. [9] analyzed via classical molecular dynamics simulations the heat transfer through nanometer-scale interfaces consisting of n-decane (2-12 nm diameter) droplets in water. Tanaka et al. [10,11] studied molecular dynamics of water by microwave heating. In contrast to these classical-mechanics-based approaches, Fisher [12] showed that the Landauer transport formalism can be applied to the formulations of the quantum thermal conductance of heat flux carried by phonons between hot and cold reservoirs on the basis of nanoscale models of inorganic materials. Furthermore, using the Landauer formulation, Rego et al. [13] studied the quantized thermal conductance of dielectric quantum wires at low temperatures.
On the other hand, attention has recently been focused on the nanoscale thermal conduction in biological systems as well. Concerning experimental reports on endogenous heat production in single cells, Baffou et al. [14] argued that the experimentally observed temperature rise of ∆T (1-2 K) of a whole cell [15][16][17] would be much larger than the theoretical ∆T estimated from the relation ∆T = P/κL, where κ is the thermal conductivity, L is the diameter of the heat source and P is the power. This equation is derived from the macroscopic heat diffusion equation for continuum [18]. Baffou et al. [19] estimated ∆T within a cell to be around 10 −5 K. However, the 10 5 gap issue (the difference between the measured temperature increase of 1 K in a single cell using cellular thermometry as compared to the theoretically calculated increase estimate of 10 −5 K) may disappear [15] if once the uncertainties of variable heat sources in stimulated cells (10 1 ), length scales (10 1 -10 2 ), and microand nanoscale thermal parameters (10 1 -10 2 ) are considered. A forthcoming challenge in the series of debates [14,15] is thus to establish a relevant theoretical framework for thermal conduction by setting up nanoscale models of biological cells containing liquid water and proteins. Regarding this issue, Lervik et al. [20] calculated the thermal conductivity at the protein-water interface in terms of classical molecular dynamics. Pandey and Leitner [21,22] quantum-mechanically evaluated the thermal energy transport through a trehalose layer between water and protein, and between gold, such as a gold nanoparticle, and its cellular environment.
In this paper, we aim at constructing a theoretical formulation to describe the quantum thermal conduction at nanoscale, paying special attention to the roles by water in biological contexts. We here derive the fundamental formula of phonon transport in water for the first step of evaluations of quantum thermal conductances in intracellular molecular environment. In the nanoscale space relevant to intracellular environment, we expect that the quantum characteristics of heat transfer may manifest themselves, which should be compared to the classical descriptions. The model of this study, which is given by a novel combination of the Landauer formulation and spectroscopic data, is powerful enough to be applicable to molecular systems containing proteins, and to quantitatively evaluate the thermal conductive characteristics in realistic systems, since experimental spectroscopic features are taken into account in the model. In the present formulation, we have applied the Green's function method (GFM) to water with a one-dimensional phonon model. This model then takes into account not only intramolecular vibration modes but also intermolecular ones. We have extracted the dispersion curves from experimental absorption spectrum of water and evaluated quantitatively quantum heat conduction characteristics such as the phonon transmission function and the thermal conductance. Based on the GFM theoretical framework and the extracted parameters, we have calculated the quantum thermal conductance and compared it with the literature values quantitatively. We show detailed mathematical formulations based on equilibrium Green's functions, non-equilibrium Green's functions and phonon transmission function in the Appendix A, and describe simplified formulas used in the calculations in the following sections of main text. Figure 1 illustrates the definitions of displacements, masses, and spring constants of the elements H 2 and O in the present one-dimensional phonon model. Figure 1 gives a background of the theoretical formulations. The j-th unit cell of water molecule with length l j consists of two elements, H 2 and O, blue and yellow colored circles, which have displacements u j,1 and u j,2 , masses M j,1 and M j,2 , and spring constants K j,1 and K j,2 , respectively. By applying the Bloch theorem u j+1,n = exp (ikl j )u j,n to Equations (A1a)-(A2b) in Appendix A, we can obtain a phonon dispersion relation of liquid water between the frequency ω and the dimensionless wavenumber:

Extraction of Phonon Dispersion Curves of Liquid Water from Experimental Results
where j and n in K j,n and M j,n are indices of unit cells and atomic components, respectively. In the model of Figure 2, j and n take values from (−N L + 1) to (N + N R ) and 1, 2, respectively, where N L , N R , and N represent the number of unit cells in the left and right reservoirs and the water, respectively, and n = 1, 2 mean two-hydrogen (H 2 ) and one-oxygen (O), respectively. The phonon dispersion relation (1) has corresponding dispersion curves with real k branches of acoustic and optical modes. Hereafter, assuming that M j,1 and M j,2 do not depend on j, we have approximated them by M 1 and M 2 , respectively, and similarly, K j,1 and K j,2 were also approximated by K 1 and K 2 , respectively. If we assume the acoustic and optical modes have angular frequency bands of 0 ≤ ω ≤ ω r and ω q ≤ ω ≤ ω p , respectively, we can obtain equations for K 1 , K 2 and ω p as follows:  One-dimensional phonon model of reservoir-water-reservoir structure. Each unit cell j has length l j , (j = −N L + 1, −N L + 2, · · · N + N R ) along the x-direction and temperature θ j , (j = 0, 1, 2, · · · N + 1), where θ −N L +1 = θ −N L +2 = · · · = θ 0 and θ N+1 = θ N+2 = · · · = θ N+N R [23].
As illustrated in Figure 2, we have assumed the model in which phonons are excited as quantum waves from the left or right reservoir and propagate from left to right of water layers or vice versa.
These reservoirs have periodic structures of solid in the red and blue areas of Figure 2, whereas water may have an aperiodic structure of liquid in the rainbow colored area of Figure 2; the reason for employing the aperiodic structure is to generally describe disordered liquid water. We have assumed a one-dimensional phonon model with intervals of water molecule as l j in the j-th unit cell, where j = 1, 2, · · · N as shown in Figure 2. Each unit cell of water is modelled to consist of two components, H 2 and O. The present model takes into account not only intramolecular vibration modes but also intermolecular ones. Phonon has two kinds of modes, acoustic and optical modes. The components of water molecules at acoustic and optical modes vibrate with the same and inverse directions, respectively. We have extracted the dispersion curves from the experimental absorption spectrum of H 2 O in the liquid phase [24] by the following procedure: 1.
Extracting wave numbers k b and k e at yellow and red broken lines of absorbance maxima in the experimental absorption spectrum of water in the liquid phase shown in Figure 3.

2.
Calculating angular frequencies ω q and ω r as 2πk b c and 2πk e c, respectively, with c being the light velocity, so that the one-dimensional phonon model is consistent with the experimental spectrum. 3.
Calculating spring constants K 1 and K 2 by using Equations (2a) and (2b) above.

5.
Calculating the wave number k of phonon by using the dispersion relation Equation (1) above.

6.
Obtaining dispersion curves in the upper part of Figure 4 below with the frequency on the abscissa and the normalized wave number kl j /π on the ordinate. . Absorption spectrum of water in the liquid phase [24]. Yellow and red broken lines represent the vibration modes (b) and (e), respectively. Figure 5 shows liquid-phase water motions [25,26] of (a) anti-symmetric stretching, (b) symmetric stretching, (c) bending, (d) libration, and (e) intermolecular vibration, respectively. In step 1, we determined wave numbers k b and k e as modes (b) and (e) from the maximum value of the experimental absorption spectrum shown in Figure 3, where we selected k b as maximum value of intramolecular vibration in the experimental spectrum because the absorbance of mode (b) was higher than that of mode (a) [27]. It is also noted that some other modes such as those associated with hydrogen bond stretch may be contained in the "libration" region (d) in Figure 3, while they are supposed to play an insignificant role in the present simplified model. This procedure is applicable not only to liquid water, but also to water in a wide range of solid, liquid, and gas temperatures [28]. Table 1 shows the wavenumbers and angular frequencies of various molecular vibrational modes of liquid water, where k a , k b , k c , k d , and k e are wavenumbers of liquid-phase water motions [25,26] of (a)-(e), respectively. Because of small difference between k a and k b , we have calculated k a as k a − k b + k b , where k a and k b are wave numbers of isolated water molecule 3756 cm −1 and 3657 cm −1 , respectively [27].

Formulation of the phonon transmission function through H 2 O
To obtain the physical outlook without losing physical essence, we have derived a simple expression for one H 2 O molecule on the basis of the formula derived from the equation of the phonon transmission function for a region consisting of N water molecules, described in Appendix A.
In this section, we derive equations of the phonon transmission function tr(T) for H 2 and O components of one water molecule, that is, the N=1 phonon model of Figure 6 for simplicity, and perform test calculations using the formula based on the N=1 water layer model without losing physical essence. We set N = 1 in Equations (A12), (A14c), (A14d), (A15a), (A15b), (A23a), (A23b), (A24), and obtained the following equations: where tr means trace of matrix; Γ L,jn , Γ R,jn , G D,jn , g L,jn and g R,jn (j, n = 1, 2) are jn elements of matrices Γ L , Γ R , G D , g L and g R , respectively; ε, ε j,n and γ j,n (n = 1-4) represent ω 2 , squared on-site energy and squared hopping energy, respectively, as defined in Equations (A2c)-(A2i) and (A6a)-(A6c) in Appendix A; Im(Z) is the imaginary part of complex number Z, and η is the imaginary part of the energy, which is a positive infinitesimal. T represents the transmission matrix: as defined in Equation (A24) in Appendix A. The complex dispersion curves shown in the upper part of Figure 4 has three modes of ( I ) acoustic wave, (II) attenuation and (III) optical wave as shown in Table 2. In this study, we have employed an interval of water molecule l 1 as 0.319 nm [29] in Figure 2.
We can calculate the transmission function by Equation (3a), and ( I ), (II) and (III) modes have the transmission probabilities of approximately 1, 0, 1 over angular frequency ranges of 0 < ω < ω r , ω r < ω < ω q , and ω q < ω < ω p , respectively, as shown in Table 2. Then we can depict the frequency characteristics of the transmission probability as shown in the lower part of Figure 4. The transmission probability takes not only discrete values of 0 or 1, but continuous values between 0 and 1 near the boundary angular frequencies of ω r and ω q .

Validation of Calculated Thermal Conductance
Since phonons are bosons, they accord to the Bose-Einstein distribution as follows: whereh, ω, k B , θ j are reduced Planck constant, angular frequency, Boltzmann constant, and temperature in the jth unit cell of water molecules in the one-dimensional phonon model of Figure 2. We can thus formulate the Bose-Einstein distribution function at the unit cell j in one-dimensional phonon model of reservoir-H 2 O-reservoir structure shown in Figure 2. Moreover, we can define the heat flow J by the following equation [1,23,31]: where the functions N 0 (ω) and N N+1 (ω) are the Bose-Einstein distributions at left and right reservoirs, respectively. Then, thermal conductance σ is defined by an equation that divides the heat flow J by the temperature difference θ 0 − θ N+1 as follows: Figure 7 shows temperature dependencies of thermal conductances by this calculation with Equation (7). The temperature dependences show fair agreement with the previous results [32] for homogeneous and heterogeneous atomic chains, where the mass of a 'Device' atom in the homogeneous case is 4.6 × 10 −26 kg (Si atom assumed), and the masses of 'Device' atoms in the two heterogeneous cases are 9.2 and 2.3 × 10 −26 kg, respectively. It is observed that the agreement is in particular better when the temperature difference |θ 0 − θ N+1 | is smaller. It is noted that the conductance curves in Figure 7 are not so sensitive to the atomic masses employed in the calculations.  (7). The parameters 1, 2, . . . 10 with unit of degrees Kelvin in the graph mean temperature differences θ 0 − θ N+1 between the left and right reservoirs.

Classical Limit of Thermal Conductance
Because the higher-order terms (hω/k B θ j ) 2 /2! + (hω/k B θ j ) 3 /3! + · · · can be neglected in the classical limithω/k B θ j 1, the Bose-Einstein distribution function N j (ω) of Equation (5) takes the following limiting form: For example, in a case of temperature θ j =300K and angular frequency range ω ≤ ω r , we havē hω/k B θ j ≤ 0.912. In this case we may use the approximation of Equation (8). {N 0 (ω) − N N+1 (ω)} of Equation (6) can then be approximated to (k B /hω)(θ 0 − θ N+1 ). The phonon transmission probability ∑ 2 k=1 ∑ N j=1 {T} j,k;j,k can be approximated to unity in the bands of 0 < ω < ω r and ω q < ω < ω p , and zero for the other bands of ω according to the lower part of Figure 4. From these two approximations, the integral with respect to angular frequency of Equation (6) can be approximated by the following equation, where we have transformed the 3rd expression in Equation (9) by using −ω q + ω p ω r . Equation (9) can be used for the quantum heat flux due to phonon propagation at the molecular size in the cell, which may be formulated according to the actual size of cell biology. In the classical limit, thermal conductance σ is thus expressed by Equations (7) and (9) as follows: By using Equation (10), we can evaluate the thermal conductance in the classical limit as σ = 7.87 × 10 −11 W/K, which agrees well with the values in the high-temperature limit in Figure 7 showing the order of 10 −10 W/K; this evaluation is also similar to the estimated value of about 100 MW/K/m 2 by Lervik et al. [20] with inclusion of the effective cross section (1 nm) 2 . The quantum thermal conductance evaluated in the present study thus agrees with the thermal conductance obtained via classical molecular dynamics simulation in the high-temperature, classical limit. One of significant results in the present study is that the magnitude of thermal conductance σ in the classical limit calculated from ω r is similar to that by classical molecular dynamics for protein-water system, whereas the former would give somewhat lower estimate due to the neglect of the contributions from other modes. Baffou et al. argued [19] that the observed temperature rise of ∆T (1-2 K) of a whole cell [16,17] would be much larger than the theoretical ∆T estimated from the relation ∆T = P/κL, where L is the diameter of the heat source, P is power, and κ is thermal conductivity estimated to be around 1 W/K/m [19]. They estimated ∆T within a cell to be around 10 −5 K, assuming that the heat source size is 10 nm and that the average heat production per hour based on glucose production is 100 pW (360 nJ/h). However, one may expect that once three uncertainties below are supposed, the 10 5 gap issue may disappear [15]: 1.
Corresponding to these issues, we can make use of analytical model dicussed above. In order to assess these issues, in this paper, we constructed the following models in the framework of nonequilibrium Green's function method: 1.
reservoir model representing heat and bath reservoirs; 2.
one-dimensional phonon model with water molecules sandwiched between the reservoirs; 3.
phonon transports at atomic scale in water.
Thus Equation (7) can provide quantitative suggestions for the discussions on intracellular temperature distributions [14,15,19,33]. At physiological temperature, our quantum-mechanical estimation for the thermal conductance above is consistent with the values evaluated through classical molecular dynamics simulations for protein-water systems [20]. Then, if we would employ the evaluations of thermal conductivity κ by Lervik et al. [20] as well, the values of κ would fall around 0.1-0.3 W/K/m, which could make the discussions [14,15,19,33] more quantitative in nanoscale intracellular conditions. Of course, further theoretical investigations would be required to resolve the gap issue.

Conclusions
We have derived the fundamental formula of phonon transport in water for the first step of evaluations of quantum thermal conductance by using a one-dimensional phonon model. We have calculated the dispersion curves reproducing OH stretching mode of water molecule and inter-molecular mode, and reproduced experimental values of the absorption spectra by fitting the parameters of the phonon modes such as ω r , ω q and so on. We also have formulated the phonon transmission function of N water molecules connected to the left and right heat baths by using the non-equilibrium Green's function method. We have employed a model in which phonons are excited as quantum waves from the left or right reservoir and propagate from left to right of H 2 O layer or vice versa. We have modelled these reservoirs as having periodic structures, whereas we can regard the water layer part as aperiodic structure of liquid. In this formulation, we have applied the Green's function method to H 2 O with a one-dimensional phonon model. The model of this formulation is powerful enough to be applicable to molecular system containing proteins, and to evaluate thermal conductive characteristics thereof.
We have extracted dispersion curves from experimental absorption spectrum of H 2 O and quantitatively evaluated quantum heat conduction characteristics such as the phonon transmission function and the thermal conductance. In addition, we have simplified the formulation of the transmission function by employing one water molecule (N=1) in the aperiodic layer. From this calculation, we have obtained the characteristic that the transmission probability is almost unity in the frequency bands of acoustic and optical modes, whereas the transmission probability vanishes reflecting the quantum tunnel effect in the frequency band between these two modes. The classical limit of the thermal conductance calculated by our formula and the literature value [32] in high temperature regime (>300 K) have agreed with each other on the order of 10 −10 W/K. The present model has also provided the calculated results for thermal conductance that are consistent with those evaluated via classical molecular dynamics simulations for protein-water systems [20] in the classical limit.
As future issues, in order to extract more information from the experimental results of water absorption spectra [34], we will expand our one-dimensional phonon model to higher-dimensional models including other vibration modes [35]. We expect this refinement of our model to provide more accurate reproduction of heat conduction characteristics. We are aiming at multi-scale heat conduction analysis by modeling not only water but also intracellular molecules such as proteins. We would like to apply the present method also to the analysis of heat conduction characteristics of biomolecules containing surrounding molecules such as flexible sugar chains [36], which would give more insights into the thermal conductive phenomena in intracellular crowding environments.

Conflicts of Interest:
The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

GFM
Green's function method Trad/s unit of tera radians per second tr trace of matrix Im imaginary part of complex number

Appendix A. Phonon Transport Equations in a One-Dimensional Chain Model of Water
We perform derivations in Appendix A.1: transport equations of phonon in water, Appendix A.2: dispersion equations in left and right reservoirs, Appendix A.3: non-equivribrium Green's functions, and Appendix A.4: quantum thermal conductance in water in the following.

Appendix A.1. Derivation of Transport Equations of Phonon in Water
We can formulate phonon transport equations in a one-dimensional model of water in terms of atomic displacements u j,n (t), j = −N L + 1, · · · − 1, 0, 1, · · · N + N R based on Figure 1 as follows: where n = 1, 2 in u j,n (t) mean two-hydrogen (H 2 ) and one-oxygen (O) atomic components, respectively. This model can mimic not only intramolecular vibration modes but also intermolecular ones. We can obtain time-independent equations of Fourier transforms of Equations (A1a) and (A1b) as follows.
The 2 × 2 matrices above consist of on-site energy and hopping energy of atoms in each water molecule of the j-th unit cell which are defined by following equations:

. Derivation of Eigenvalue Equations in Left and Right Reservoirs
We define κ and λ as Fourier transforms of atomic displacements of ±x propagation waves in the left reservoir, and µ and ν as those in the right reservoir, respectively. We assume that the Fourier transforms of atomic displacementsû j,n are linear combinations of ±x propagating waves in the left and right reservoirs as follows: where T j and R j (−N L + 1 ≤ j ≤ 0 and N + 1 ≤ j ≤ N + N R ) are amplitudes of + and -x propagating waves, respectively. Assuming that the phonon wavenumbers in the left and right reservoirs are k ± L and k ± R , respectively, we derive the following equations from the periodic boundary conditions of amplitude based on Bloch's theorem for the left and right reservoirs: In the left reservoir, applying R j = 0 and T j = 0 as the condition of no phonon propagating to ∓x directions, we obtain eigenvalue equations of plane waves of phonon propagating to ±x directions.
Similarly, in the right reservoir, applying R j = 0 and T j = 0 as the condition of no phonon propagating to ∓x direction, we obtain eigenvalue equations of plane waves of phonon propagating to ±x directions.
We can calculate dispersion curves as solutions to the eigenvalue Equations (A9a)-(A10b) in the left and right reservoirs, respectively, and these dispersion curves are the same as the solutions of the dispersion relation (1).

Appendix A.3. Non-Equilibrium Green's Functions
The purpose of this section is to formulate phonon waves propagating in the water excited from left and right reservoirs, where we use the non-equilibrium Green's functions defined by arranging Equation (A4a) [23]: where G L , G D , and G R are retarded Green's functions in the left reservoir, water molecules, and right reservoir, respectively; G LD , G DL , G RD , and G DR are retarded Green's functions at the interfaces between each reservoir and water. η is the imaginary part of the energy, which is a positive infinitesimal. Eliminating G L , G R , G LD , G DL , G RD , and G DR from Equation (A11), we can get an expression of G D as follows: where Σ L and Σ R are self-energies defined as g L and g R are equilibrium Green's functions defined as follows: g j,j = g j,1;j ,1 g j,1;j ,2 g j,2;j ,1 g j,2;j ,2 , j, j = −N L + 1, −N L + 2, · · · , 0, N + 1, N + 2, · · · , N + N R . (A14e) Substituting Equations (A4e)-(A4h), (A14c) and (A14d) into (A13a) and (A13b), we can obtain the following equations for the self-energies: We can define the equilibrium Green's function equivalent for each unit cell at each atomic position (n, n ) in the left and right heat reservoirs as follows: g L n,n = g j,n;j ,n , j, j = −N L + 1, −N L + 2, · · · , 0, (A16a) g R n,n = g j,n;j ,n , j, j = N + 1, N + 2, · · · , N + N R .

. Quantum Thermal Conductance in Water
Let us formulate equilibrium Green's functions in reservoir, g L , g R , by the the following procedure consisting of basis transformations and so on: Transforming wave functions consisting of the tight binding (TB) basis into those consisting of the plane wave (PW) basis.

2.
Finding the PW basis solutions in semi-infinite (N L , N R → ∞) heat reservoir at both ends of water with using the periodic characteristics of the heat reservoir in thermal equilibrium.

3.
Inverse transformation from PW basis wave functions to TB basis ones.
Based on the components of the j=1st and Nth unit cells of Equation (A4a), we can obtain the following equations governing the plane waves propagating to x direction in the left and right heat reservoirs, respectively: We define P L and P R in Equations (A18a) and (A18b) as matrices in which (1,1) and (2,2) components of the diagonal terms are lined with propagation constants along x direction in the left and right heat reservoirs, respectively: Next, we define the matrices U L and U R in Equations (A18a) and (A18b) as unitary matrices consisting of the eigenvectors of the eigenvalue Equations (A9a)-(A10b) as follows: We can simplify g L and g R as matrices with the same number of components as the number of the lattice points connected to the left and right heat reservoirs by using Equations (A18a) and (A18b), respectively: We can obtain spectral functions A L and A R proportional to local state densities A L /2π, and A R /2π in left and right reservoirs, respectively: We can obtain functions Γ L and Γ R by the following equations, which are matrices representing the phonon lifetime in the left and right reservoirs, respectively: We can then formulate the transmission matrix T of phonon waves traveling from the left to the right heat reservoirs via the interfaces of the water as follows [1,23,31]: We can thus express the heat flow in water J by the following equation by integrating the product of the transmission function and the Bose-Einstein distribution [32]: {T} j,n;j,n is the transmission function expressed by the trace of T; N 0 (ω) and N N+1 (ω) are Bose-Einstein distribution functions at j = 0 and N + 1 expressed as follows [37]: Thermal conductance in water σ can then be expressed as