Molecular dynamics study on thermal energy transfer in bulk polyacrylic acid

We performed non-equilibrium molecular dynamics (NEMD) simulations on bulk amorphous polyacrylic acid (PAA) with three polymer chain lengths to investigate molecular mechanism of thermal energy transfer in heat conduction. Thermal conductivity obtained by NEMD simulations increased as the polymer chain length of PAA increased, and its dependence on polymer chain length exhibited a saturation behavior. By decomposing heat flux into each contribution of molecular interactions, it was found that dominant mechanism of the thermal energy transfer in PAA was intramolecular interaction, and contribution of the intramolecular interaction to thermal conductivity increased as the polymer chain length increased, and resulted in increase in total thermal conductivity. On the other hand, coiled conformation of PAA advanced in response to elongation of the polymer chain length; and this coiled conformation inhibited further increase of thermal conductivity due to the polymer chain elongation. Consequently, the elongation of the polymer chain length had two conflicting effects: increasing and suppression of thermal conductivity, due to increase in intramolecular interaction and change in conformation, respectively. This is the reason of the saturation tendency of thermal conductivity as a function of the polymer chain length. Detailed understanding of the molecular mechanism of thermal energy transfer obtained in the present study provided the in-depth knowledge to clarify the thermal energy transfer mechanism and will lead to the characterization of thermal energy transfer in more complicated materials such as a layer-by-layer membrane. © 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.5080432


I. INTRODUCTION
Thermal management in mechanical and electric devices is important to obtain high performance of such devices. Recently, due to improvement and downsizing of power semiconductor and other electronic devices, heat generation density has been increasing, and hence, further reduction of interfacial thermal resistance between two joining components in such devices has become a requirement. In order to reduce the interfacial thermal resistance at the joining surfaces, thermal interface material 1-3 (TIM) is inserted between the two components. For TIM, various soft materials such as thermal grease with suspended microparticles, elastomer matrix with fillers particles, and gels are utilized, and development of the TIM to attain high thermal conductivity has been continuing. For example, dispersing filler particles with high thermal conductivity into the TIM raises the overall thermal conductivity of the TIM. However, in the case of the TIM with fillers, the limit on the amount of fillers that can be added and that material cost increases by the filler's cost becomes a bottleneck. Furthermore, forming continuous heat transfer path via fillers in the TIM is not sufficient to meet the industrial high requirements. Under such circumstances, search for novel materials for TIM has been continuing.
Several types of soft matter have peculiar heat conduction characteristics. For example, lipid bilayer membranes, in which phospholipid molecules are aligned in a direction normal to the membrane and form a bilayer structure, have a cross-plane thermal conductivity five times higher than the in-plane one. 4 This asymmetric thermal conductivity results from the fast thermal transport along the covalent bonds in the molecules aligned in the cross-plane direction. Among such soft matters, the layer-by-layer (LBL) membrane has become the focus of the authors' research for promising materials for the next-gen TIMs. The LBL membrane is fabricated by the sequential deposition of positively charged and negatively charged material. For the deposition materials, such electrifiable materials as polymer, 5,6 protein, 7 carbon nanotube 8 and graphene 9 are used to produce various features. 10 The LBL membrane is made by simple and low cost fabrication and it is possible to form the membrane on substrates with various surface structures. Moreover, it is easy to selectively produce a membrane with the thickness at nanoscale by the sequential deposition protocol. The thickness of the TIM is an important factor determining the thermal resistance, and hence, ease of control of the thickness is a great advantage for the TIM.
In the viewpoint of thermal energy transfer in LBL structure, it is expected that thermal energy transfer of the materials itself works for the in-plane heat conduction within the layer, while thermal energy transfer between the two materials works for the cross-plane heat conduction. Especially, LBL membranes made by polymer molecules with hydrophilic groups are expected to exhibit high thermal conductivity helped by intermolecular energy transfer due to electrostatic forces. For example, Kim et al. reported that a blend of polyacrylic acid (PAA) and poly-N-acryloyl piperidine (PAP), where PAA is hydrogen-bond donating polymer and PAP is hydrogen-bond accepting polymer, has thermal conductivity about seven times higher than pure PAA and PAP. 11 We are thinking about application of polyethylenimine/polyacrylic acid (PEI/PAA) LBL membrane. Prior to the examination of the LBL membrane, thermophysical properties and molecular mechanism of thermal energy transfer in the bulk of the two materials were analyzed. In the present paper molecular dynamics simulation of bulk PAA is performed to examine molecular mechanism of thermal energy transfer. Although thermal conductivity of dry amorphous PAA 11,12 and PAA solution in water 13 were measured experimentally, the detailed molecular mechanism of thermal energy transfer that determines the thermal conductivity is not clear. In view of this, non-equilibrium molecular dynamics (NEMD) simulations of bulk PAA with three polymer chain lengths were performed. Furthermore, the heat flux decomposition analysis, 14 which has been usefully applied to various liquids and soft matters, e.g., bulk liquid, 15 binary mixture, 16 solid-liquid interface, 17,18 polymer 19 and lipid bilayer membrane, 4 helped us to understand how the various molecular mechanisms contribute to heat conduction in PAA.

A. Molecular model
The CHARMM General Force Field (CGenFF), 20,21 which is an all-atom force field and includes the protein, nucleic acid, lipid, carbohydrate and general organic molecule parameter sets, was used for PAA in the present simulations. In this study Ligand Reader & Modeler in CHARMM-GUI 22,23 was used to create the force field parameters and structures of PAA because there are no specific parameters for PAA in CGenFF. The details of the parameters used in our simulations are given in the supplementary materials (see Table S1-S5). The chemical structure of PAA is shown in Fig. 1, where n is the number of monomeric units. PAA consists of the monomeric units which have a main chain (-C-C-) with three hydrogen and a carboxyl group. Chain ends are terminated by hydrogen. In the present study, we used three types of PAA with different polymer chain lengths, n = 5, 10 and 20, which are referred to as PAA5, PAA10 and PAA20, respectively. In addition, atactic polymers were used and dissociation of PAA was not considered.

B. Simulation system
The NEMD simulation system used in the present study is shown in Fig. 2. The simulation system is a rectangular parallelepiped unit cell with periodic boundary conditions in the x, y and z directions. The lengths of the unit cell in the x, y and z directions are 3L x , 3L y and 14L z . Here, L x , L y and L z are the characteristic lengths which are used to determine slab widths, and defined after NPT simulation. The simulation protocol is described below and the detail simulation conditions are also given in the supplementary materials.
First, one of the three type of PAA were randomly arranged and oriented in the unit cell of which the lengths in

ARTICLE
scitation.org/journal/adv the x, y and z directions were 12L x , 12L y and 56L z , i.e., of which the volume was 64 times larger than the one to give the target density, 1.2 g/cm 3 . The initial temperature of the system was set at 1000 K. The initial values of the lengths L x , L y and L z were chosen to be larger than the characteristic scale of the coiled PAA used here and all set at 18 Å (see supplementary materials). The unit cell was shrunk isotropically during 0.05 ns in the NVT simulation at 1000 K so that the density would reach the target value of 1.2 g/cm 3 . Before and after this shrinking process the NVT simulations were carried out for 0.025 ns as a relaxation run. After that the system was cooled gradually to 300 K in the NVT simulation for 2 ns, and then NPT simulation at 300 K and 1 atm was carried out for 10 ns until the system density reached a steady state. The length L x , L y and L z after the NPT simulation are 1/3, 1/3 and 1/14 of the each length of the unit cell in x, y and z directions, respectively, and these are shorter than the initial value of 18 Å because the unit cell shrank during the NPT simulation.
The unit cell obtained by the above processes was used for NEMD simulations. In the NEMD simulation, a constant heat flux J = 1000 MW/m 2 was applied in the z direction by using the enhanced heat exchange (eHEX) algorithm proposed by Wirnsberger et al.. 24 In this system, thermal energy was added and subtracted from the hot slabs and cold slab, respectively, at a constant rate by rescaling the velocities of all particles contained in each slab. The hot slabs with a length of L z in the z direction were set at both ends of the unit cell and the cold slab with a length of 2L z was set at the center of the unit cell in order to generate the heat flux with a temperature gradient in the z direction. To measure the system heat flux, two control volumes with lengths of 3L z in the z direction were placed on both sides between the hot and cold slabs and are separated by a length L z from both the hot and cold slabs to eliminate the effects of the artificial energy control. This NEMD simulation was carried out for 20 ns as relaxation run to ensure a non-equilibrium steady state. All calculations of above protocol were performed by using LAMMPS. 25 The smooth cutoff was applied between 10 and 12 Å for the Lennard-Jones (LJ) potential that was applied for the van der Waals (vdW) interaction. The electrostatic interactions were handled by particle-particle particle-mesh (PPPM) method 26 with a cutoff of 9 Å. The reversible reference system propagator algorithms (r-RESPA) method 27 was used as the time integration algorithm. The time step was 0.25 fs for the intramolecular covalent bond interactions, 0.5 fs for the intramolecular angle interactions and 1.0 fs for the others. The NVT and NPT simulations were handled by the Nosé -Hoover style equations. 28 After the relaxation run by LAMMPS, the NEMD simulation was carried out using our in-house program package for 3 ns as the production run. The simulation system by our inhouse program has several differences from that by LAMMPS, which are described as follows. The constant heat flux was applied by using the scheme of Jund and Jullien. 29 The cutoff for LJ potential was 12 Å without smoothing. The electrostatic interactions were handled by the smooth particle mesh Ewald (SPME) method 30 and the unit cell was divided into charge grids so that the charge grid spacing was almost 1 Å. The time steps of 1.0 and 0.2 fs for the inter-and intramolecular interactions, respectively, were used via r-RESPA method.
In the preset study, ten simulation systems with different initial structures were used for each case of polymer chain length to reduce dependence on the initial structure and to increase the number of samples for calculation of accurate thermal properties. For a single case, the production run was carried out for 3 ns, thus for each polymer chain length we obtained the amount of samples equivalent to that of 30 ns, which were used to calculate the thermal conductivity and the heat flux decomposition described later in Section II. C.
In the NEMD simulation, the constant temperature gradient in the z direction was formed in the control volume. The thermal conductivity λ can be obtained from the Fourier law where ∂T/∂z is the average of the temperature gradients which were calculated from a linear fitting to the local temperature profile in z direction in each control volume.

C. Heat flux decomposition
Macroscopic conduction heat flux in media is composed of thermal energy transfers due to various molecular mechanisms. The total heat flux J total in the z direction in the control volume V CV in a classical system is given by 14 (2) The first term in the right-hand side of Eq. (2) represents the contribution of the transport of internal energy of sites by the migration of them, and is composed of the kinetic and potential components. Here, E s and s,z are the total energy, i.e., sum of the kinetic and potential energies, and the z component of velocity of the site s, respectively. The second term of Eq. (2) represents the contributions of the mechanical work of the inter-and intramolecular forces on the sites defined by the n-body potentials. Here, F s,U is the intermolecular or intramolecular forces exerted on the site s according to the interaction of an n-body potential defined among a set of sites U = {s 1 , s 2 , · · · , s n }, v s and z s are the velocity vector and the z coordinate of the site s, respectively. In this study, the second term is composed of the contributions of bond stretching, valence angle bending, torsions of both dihedral and improper dihedral angles, vdW, and electrostatic interactions. in the brackets [ ] represents taking the summation over all pair of sites under the case where the line segment connecting the two sites is fully or partially located in the control volume, and (z sα − z sβ ) * is the portion of (z sα − z sβ ) located in the control volume, and hence, (z sα − z sβ ) * = (z sα − z sβ ) holds in the AIP Advances 9, 025302 (2019); doi: 10.1063/1.5080432 9, 025302-3 ARTICLE scitation.org/journal/adv case where the line segment between the two sites is fully located in the control volume. The total heat flux is composed of the various contributions of molecular interactions that have been described. Each of the components was obtained in the present simulations and its contribution to the heat flux was evaluated to elucidate the mechanisms that govern the transfer of thermal energy.

III. RESULTS AND DISCUSSION
A. Thermal conductivity The temperature and density distributions in the z direction of a single PAA5 system for the case where the average temperature and density in the control volumes are T ave = 303 K and ρ ave = 1.25 g/cm 3 , respectively, are shown in Fig. 3 as an example of the NEMD simulations. The computational domain was divided into 50 slabs in the z direction, and the local temperature and density for each slab was calculated considering atoms located in each slab. Linear temperature distributions from the hot slab to the cold slab were formed in the control volumes, while density distribution fluctuated around the average value. This fluctuation results from the fact that the simulated bulk PAA has amorphous solid structure and molecular migration is limited. Thus, the local density is mostly fixed in space over the simulation period and depends on the initial structure. Although the density gradient is expected to appear in response to the temperature gradient, such tendency is not recognized clearly due to spatial density fluctuation. The thermal conductivities λ for PAA with each of the three polymer chain lengths were calculated by Eq. (1) and are shown in Fig. 4 together with the average density for the three polymer chain lengths. Error bars in Fig. 4 represent the standard error of the values obtained in ten cases for each chain length. The thermal conductivity and average density both increase as the polymer chain length increases and have similar tendency. The density of the system increases, i.e., the system shrinks, as degree of polymerization increases, because the equilibrium distance of the covalent bond between carbon atoms is shorter than that of the vdW interaction. Therefore, the density of PAA increases  due to elongation of the polymer chain length. Higher density PAA has more heat transfer paths per unit volume compared to lower density PAA, additionally many paths in PAA with high degree of polymerization have fast thermal transport along the covalent bonds, thus longer PAA has the higher thermal conductivity in the range of the polymer chain length simulated here.
The thermal conductivities calculated in this study were λ = 0.302-0.351 W/(m·K). In case of experimental values of the thermal conductivity of bulk amorphous PAA (Mw = 450,000), Kim et al. reported λ = 0.22 ± 0.02 W/(m·K) which was measured by using a standard differential 3ω method in ambient conditions, 11 and Xie et al. reported λ = 0.37 ± 0.02 W/(m·K) which was measured by using time-domain thermoreflectance (TDTR) method. 12 Our simulated thermal conductivities fall within the rage of the experimental data, and hence, we concluded that our simulations are reasonable to clarify the molecular-level thermal transport mechanism in PAA.

ARTICLE
scitation.org/journal/adv of PAA. For comparison, the percentage of each contribution to the total heat flux for the linear alcohol liquids from ethanol to decanol, i.e., for the linear alcohols with various chain lengths, at the standard state, 298 K and 1 atm, 31 are also shown in Fig. 5. Here, the total heat flux is decomposed into eight components as follows: J total = J kin + J pot + J bond + J angle + J torsion + J nb,intra + J LJ,inter + J elec . (3) Here, J kin and J pot represent the contributions of the transport of the kinetic and potential energies of the molecules according to migration of the molecules, respectively, and is equivalent to the first term of Eq. (2). The other components belong to the second term of Eq. (2). J bond , J angle , J torsion and J nb,intra are the contributions of the intramolecular interaction, and represent the contribution of the mechanical work of the bond stretching, valence angle bending, torsions of both dihedral and improper dihedral angles, and the others intramolecular nonbonded interaction excluding the 1-5 or more separated intramolecular electrostatic interaction, respectively. J LJ,inter represents the contribution of the intermolecular interaction via vdW interaction. J elec represents the contribution of the 1-5 or more separated intramolecular and the intermolecular electrostatic interaction. We used the modified Wolf method 32 to calculate the per-atom electrostatic potential for the first term in Eq. (2) because it is technically difficult to calculate that by using the Ewald sum. On the other hand, because the method to calculate the stress due to the electrostatic potential in the finite control volume in the case of the Ewald sum has not been developed, the partial heat flux due to the electrostatic interaction in the second term in Eq. (2) was obtained indirectly by subtracting the partial heat flux of the other contributions from the applied total heat flux (J total = 1000 MW/m 2 ) as with the previous work for linear alcohols. 31,33 As indicated in Fig. 5, the contribution of the transport term for PAA is minuscule compared to that for the linear alcohol liquids, and is less than 4% to the total heat flux. This is because simulated PAA is amorphous solid material, and hence, atomic migration is limited compared to liquid.
Similar to the tendency in the case of the linear alcohols, for PAA the contributions of the intramolecular interaction increases as the polymer chain length increases, because the elongation of chain increases the path of thermal energy transfer utilizing the stiff intramolecular bonds. In the previous work for saturated liquid of linear alkanes 15,34 and normal alcohols at the standard state, 31 a shift of the dominant mechanism of heat transfer from intermolecular to intramolecular interaction as the chain length of molecule increases was reported, and the switching chain length, i.e., the number of carbon atoms in a molecule, from inter-to intramolecular interaction being dominant is around 10 to 16 for liquid alkanes and around 4 to 6 for linear alcohols. Additionally, in the previous work, 31,34 when the intermolecular interaction is dominant, increase of the contribution of the intermolecular interaction, i.e., decrease of the chain length, raises the partial thermal conductivity due to the intermolecular interaction, resulting in increase of the total thermal conductivity. On the other hand, when the intramolecular interaction is dominant, the partial thermal conductivity due to the intramolecular interaction and the total thermal conductivity and increase as the contribution of the intramolecular interaction increases, i.e., as the chain length increases. PAA has 2 carbon atoms on every monomeric unit along its main chain, and thus the number of carbon atoms along the main chain of PAA5, PAA10 and PAA20 are 10, 20 and 40, respectively. Therefore, PAA molecules used here are long enough for the dominant mechanism to switch to the intramolecular interaction, and indeed the sum of the contributions of J bond , J angle , J torsion and J nb,intra is over 60% of the total for all cases of PAA.
The partial thermal conductivity 35 due to a contribution of i is given by (4) Figure 6 shows the partial thermal conductivities due to each contribution according to Eqs. (3) and (4) as functions of the polymer chain length. Here, the subscripts of λ correspond to those in Fig. 5 with their colors in the figure. The partial thermal conductivity due to the intramolecular interactions, especially the bond stretching and angle bending, increase as the polymer chain length increases, and results in the increment of the total thermal conductivity. In focused molecules, carbon, oxygen, and hydrogen atoms are drawn as green, red, and white, respectively.
The contribution of the electrostatic interaction is almost constant around 11-15% and the partial thermal conductivity due to the electrostatic interaction is also almost constant for all PAA used here. This is contrary to the fact that for the linear alcohol liquids it decreases as the chain length increases. In the case of the linear alcohol, a hydroxyl OH group connects to a primary carbon atom. Therefore, the density of thermal transfer path due to the electrostatic interaction decreases as the chain length increases because the longer linear alcohol molecule has lower OH group density. In contrast, PAA has an OH group on every monomeric unit along the main chain, thus the OH group density is independent of the polymer chain length. Thus, the contribution of the electrostatic interaction is almost constant for all PAA. On the other hand, the OH group density of PAA is higher than that of the linear alcohol liquids. Therefore, the electrostatic interaction contributes to the heat transfer as much as the intermolecular vdW interaction in PAA, while the contribution of the electric interaction is smaller than that of the intermolecular vdW interaction in the linear alcohol liquids.

C. Morphology of polymer
In this part, we discuss the effect of polymer chain length on thermal conductivity from the aspect of the polymer conformation. The radius of gyration, R g , is often used to describe a size of coiled polymer chain, and roughly, the coiled polymer occupies a spherical space with radius R g . Its square R 2 g represents the second moment around the center of mass of the polymer as follows: Here, M and r CM are the mass and position vector of center of mass of the polymer, m i and r i are the mass and position vector of atom i inside the polymer, respectively. Figure 7 shows the average of the radius of gyration, R g , as a function of the number of monomeric unit n with a fitting curve. Error bars in Fig. 7 represent the standard error of the values obtained in the ten cases for each case of chain length. The fitting curve was calculated by the leastsquare method for log R g = α log n + β, where α and β are parameters. In general, the radius of gyration of a polymer consisting of n monomeric units has the relation, R g ∝ n α . The exponent α for an ideal linear chain polymer without an excluded volume is 0.5, and for a real chain polymer with the excluded volume is around 0.59. 36 As shown in Fig. 7, it was observed that α = 0.549 for PAA, and hence the tendency of the conformation of PAA depending on the polymer chain length can be considered as that of a typical linear polymer. The value of α ∼ 0.5 represents that coiling of the polymer advances in response to the increase of polymer chain length. Figure 8 shows snapshots of typical conformations of PAA in the central part of the control volume for the cases with the three chain lengths. The displayed region sizes are 36 × 36 × 36 Å 3 . Carbon, oxygen, and hydrogen atoms in focused molecules are drawn as green, red, and white, respectively, and the remaining non-focused molecules are colored in yellow. We can observe the highly coiled conformation of the long PAA as compared to the short PAA. In the case of amorphous polymer, highly coiled and entangled polymer structure inhibit further increase of thermal conductivity in response to increase in chain length. 37 As an extreme example, if TIM consisted of very long polymers and a single polymer molecule could connect from heat source to heatsink like molecular junction, 38 the TIM could have significantly high thermal conductivity due to the intramolecular heat path which can transfer the thermal energy directly. However, even if the polymer chain length increases, the dimension in which the thermal energy can be transferred directory via intramolecular interaction is limited due to the coiling conformation. This is the reason why the total thermal conductivity in Fig. 4 and the partial thermal conductivities due to the contributions of the intramolecular interactions in Fig. 6 as the functions of the polymer chain length eventually saturate.

IV. CONCLUSIONS
NEMD simulations of heat conduction were carried out to investigate thermal energy transfer mechanism in bulk PAA at the molecular scale. Thermal conductivities obtained in our NEMD simulations for PAA with three polymer chain lengths fall within the reasonable range compared to the experimental values. Thermal conductivity and density of PAA increases as the polymer chain length increases. To clarify the relationship between thermal conductivity and the polymer chain length, the heat flux decomposition and the morphological analysis were examined.
As a result of the heat flux decomposition, more than 96% of thermal energy is transferred by the mechanism of the ARTICLE scitation.org/journal/adv intra-and intermolecular interactions, while molecule migration only accounts for less than 4%, because bulk PAA has amorphous solid structure and the migration of molecules is limited. Especially, the contribution of the intramolecular interaction is dominant, with more than 60% of the total. The elongation of the polymer chain length raises the contributions of the intramolecular interaction, resulting in increase of the thermal conductivity. On the other hand, the contribution of the electrostatic interaction which is expected to be an effective heat transfer path in the LBL membrane is about 10% regardless of the polymer chain length, because in the present study PAA without dissociation was assumed and OH groups are distributed homogeneously in the polymer. In the actual LBL membrane, PAA is a polyanion, i.e., exists in ionized state. Therefore, how the electrostatic interaction enhances thermal energy transfer in the LBL membrane should be evaluated by using the ionized polymers, and hence, further investigation is needed.
Finally, the conformation of PAA with respect to the polymer chain length was investigated. The longer PAA has highly coiled conformation. Although the elongation of the polymer chain tends to raise thermal conductivity, increase in thermal conductivity is limited by the coiled conformation which advances in response to the increase of polymer chain length. Consequently, the relation between the thermal conductivity and the polymer chain length results in a saturation curve. Our work provides the proper understanding of the molecular mechanism of the thermal energy transfer in bulk PAA and provides the important knowledge in order to clarify the thermal energy transfer mechanism and the heat conduction character in the LBL membrane which is our on-going work.

SUPPLEMENTARY MATERIAL
See supplementary material for detail force field parameters of PAA and detail values of conditions and results of NEMD simulations.