Effect of boundary chain folding on thermal conductivity of lamellar amorphous polyethylene

Thermal transport properties of amorphous polymers depend significantly on the chain morphology, and boundary chain folding is a common phenomenon in bulk or lamellar polymer materials. In this work, by using molecular dynamics simulations, we study thermal conductivity of lamellar amorphous polyethylene (LAPE) with varying chain length (L0). For a short L0 without boundary chain folding, thermal conductivity of LAPE is homogeneous along the chain length direction. In contrast, boundary chain folding takes place for large L0, and the local thermal conductivity at the boundary is notably lower than that of the central region, indicating inhomogeneous thermal transport in LAPE. By analysing the chain morphology, we reveal that the boundary chain folding causes the reduction of both the orientation order parameter along the heat flow direction and the radius of gyration, leading to the reduced local thermal conductivity at the boundary. Further vibrational spectrum analysis reveals that the boundary chain folding shifts the vibrational spectrum to the lower frequency, and suppresses the transmission coefficient for both C–C vibration and C–H vibration. Our study suggests that the boundary chain folding is an important factor for polymers to achieve desirable thermal conductivity for plastic heat exchangers and electronic packaging applications.


Introduction
Recent studies suggest that polymers have great potential applications in many elds, such as thermoelectrics, 1,2 solar cells, 3,4 light-emitting diodes, 5 biomedical devices, 6 nonlinear optical devices, 7,8 and electronic packaging. 9 As the most commonly used material for electronic packaging, the heat transfer capability of a polymer lm directly determines the heat dissipation and operating performance of the electronic device. In this regard, a comprehensive understanding of the thermal transport properties of amorphous polymers is crucial. The morphology of the surface polymer chain in a bulk polymer is mainly the lamellar chain-folded type, which is formed during the crystal growth process from quiescent melt or solutions of exible long-chain polymers. 10,11 Actually, many polymeric materials contain surface chain-folded lamellar and intervening amorphous regions. 12,13 Due to a lower free-energy barrier compared to the inter-molecular secondary nucleation, the chain folding is attributed to intra-molecular secondary crystal nucleation at the crystal growth front. 10 During the growth of polymer, the conguration of adjacent chain-folding can limit the thickness of polymer. 14 The bulk amorphous polymers typically have a low thermal conductivity around 0.1-1.0 W m À1 K À1 at room temperature, 15,16 which is two to three orders of magnitude lower than that of crystals. 17,18 The random orientation of the polymer chains and the weak inter-chain couplings are the two main mechanisms responsible for the low thermal conductivity. 19 Proper alignment of polymer chains is found to be an effective method to enhance the polymer's thermal conductivity, which can achieve more than two orders of magnitude enhancement for amorphous polyethylene (PE). 17,20 Using molecular dynamics simulations, Henry et al. 21,22 found that thermal conductivity of a single PE chain can reach up to 300 W m À1 K À1 , which is almost as high as that of some inorganic semiconductors and metals.
The morphology of polymer chains also plays an important role on the thermal conductivity of polymer. Zhang et al. 23 used PE as a model system to systematically investigate the fundamental relationship between the molecular morphology and thermal conductivity in bulk amorphous polymers. They found that thermal transport through covalent bonds dominates the effective thermal conductivity in polymers, and the thermal conductivity of rigid rod-like polymer can be increased by more than one order of magnitude, compared to that of polymers with very so backbone. For lamellar amorphous polyethylene (LAPE), Ma et al. 24 found that stronger connement and less entanglement lead to a smaller thermal conductivity. However, the impact of boundary chain-folding on thermal transport in lamellar polymer has not been understood yet. This is particularly relevant for lamellar polymer because the surface-tovolume-ratio is rather large 7 due to the limited lm thickness, and the inuence of chain folding cannot be neglected.
In this work, we study the thermal conductivity of LAPE by using molecular dynamic simulation to elucidate the effect of chain folding on thermal transport. Our results show that the chain folding results in the non-uniform local temperature gradient in LAPE. To describe the thermal transport property of the whole system in the presence of non-uniform temperature gradient, the validity of using effective thermal conductivity is discussed. The physical origin for the non-uniform local thermal conductivity is revealed by comparing various factors in different regions, including chain morphology, orientation, and vibrational spectrum. Our study provides a solid understanding for the inuence of chain conguration on the thermal transport in LAPE, which will be benecial for the design of polymer based materials for thermal management applications.

Model and simulation
To take into account the variation of polymer chain length in experiments, we rst construct PE chains with varying lengths, which are sufficiently extended and quasi-straight by randomly changing the torsion angles (see Fig. S1(a) in ESI †). Here the chain length refers to the total length of the polymer chain when it is fully extended (straight). The chain length distribution satises the Gaussian distribution where L 0 is the average length, and s is the standard deviation. Table 1 summarizes the detailed congurations for various PE chains with different lengths with a given standard deviation of s ¼ 1 nm. As the total number of atoms in our study is xed (about 27 000), the number of chains will decrease with the increase of L 0 . We then use the PACKMOL soware 25,26 to randomly place the PE chains into a sufficiently large box (50 nm Â 10 nm Â 10 nm), which allows for the construction of the initial lamellar structure (see Fig. S1(b) in ESI †). Aer preparing the initial structure, molecular dynamics simulations are performed to study the dynamics of the system by using the large-scale atomic/molecular massively parallel simulator (LAMMPS) soware. 27 The adaptive inter-molecular reactive empirical bond order (AIREBO) potential 28 is used to describe inter-and intra-molecular interactions between carbon and hydrogen atoms, which has been widely used to study the polymer systems in literature. 21,22 The time step is 0.2 fs, and the force cutoff distance is set as 10Å. Periodic boundary conditions are adopted in all directions.
The initial structure is relaxed by performing equilibrium molecular dynamics (EMD) simulation with the isothermalisobaric (NPT) ensemble at 1 atm and 300 K for 2 ns, aer which the thermal equilibrium and structure relaxation have been achieved (see Fig. S2 in ESI †). More details about the structure generation are included in the ESI. † EMD simulations with canonical (NVT) ensemble at 300 K are then performed for another 1 ns. A typical relaxed atomic structure of LAPE is shown in Fig. 1.
Aer structure relaxation, we then carry out the nonequilibrium molecular dynamics (NEMD) simulations to compute the thermal conductivity of LAPE. As shown in Fig. 1, two Nosé-Hoover thermostats at 310 K and 290 K are applied to the center and two ends of the simulation domain, respectively, in order to establish the temperature gradient along the LAPE. The lattice thermal conductivity k of the amorphous PE is then calculated as: where VT is the temperature gradient along the length (X) direction, and J is the heat ux dened as the energy transported across unit area per unit time 29 where dQ/dt is the average rate of the energy injection and extraction in the thermostated regions, and S is the crosssectional area. The factor of two in the denominator accounts for the bidirectional heat ow ( Fig. 1) due to the use of periodic boundary condition. We record the temperature prole along the X-direction by dividing the simulation domain into slices with a uniform width of 2Å. The temperature gradient is then calculated by tting the linear portion of the temperature prole. The NEMD simulations are performed long enough (2 ns) to ensure the non-equilibrium steady state in which the temperature prole and heat ux are time-independent. The nal thermal conductivity results are averaged over three independent NEMD simulations with different velocity distribution as initial conditions. It should be noted that the thermal transport behaviors in nanomaterials are signicantly different from that of bulk materials. For instance, in one-dimensional nanotubes and nanowires, [30][31][32][33] and two-dimensional graphene, 34 the peculiar size dependent thermal conductivity has been reported. Therefore, we focus on the effect of boundary chain folding on the thermal conductivity of PE for the nite-size system, which allows us to discover some new phenomena different from bulk PE.

Effective thermal conductivity
To test the AIREBO potential for thermal conductivity of PE, we calculate the thermal conductivity of bulk PE, which is found to be 0.34 W m À1 K À1 at room temperature. This result falls in the range of 0.22-0.37 W m À1 K À1 reported by the literature study [35][36][37] for the bulk PE. Furthermore, the nal structures for different LAPE samples with varying lengths have a similar density about 0.70 g cm À3 aer NPT relaxation (Table 1), which is in good agreement with the reported density 0.68-0.73 g cm À3 at ambient condition. 38 Fig. 2 shows the typical temperature prole along the length direction for LAPE with different average length (L 0 ) and constant standard deviation (s ¼ 1 nm). For short PE chains (L 0 ¼ 15 nm), a linear temperature prole (square in Fig. 2) is established in the whole region between the heat source and heat sink. In contrast, the temperature prole for long PE chains (L 0 ¼ 30 nm) is no longer linear across the whole region, but becomes piecewise linear, exhibiting a two-segment kinked characteristic (circle in Fig. 2). Besides, the temperature gradient in the boundary region (X < 20Å) is greater than that in the internal region (X > 20Å). Furthermore, we repeat the simulations for L 0 ¼ 30 nm with Langevin heat bath, and the same non-uniform characteristics in temperature prole are also observed (see Fig. S4 in ESI †).
According to the Gaussian distribution, the length of various chains mainly falls into the range from L 0 À 3s to L 0 + 3s. When the upper bound of the chain length (L 0 + 3s), for instance L 0 ¼ 15 nm, is smaller than the length of the simulation domain (L ¼   1 nm). The results are averaged over two mirror symmetry regions. The square and circle represent the temperature profile for L 0 ¼ 15 nm and L 0 ¼ 30 nm, respectively. In the presence of non-uniform temperature gradient for L 0 ¼ 30 nm, k 1 and k 2 denote the local thermal conductivity in the boundary and internal region, respectively. The inset depicts the calculation of effective thermal conductivity k 0 for L 0 ¼ 30 nm based on the artificial uniform temperature gradient between two ends. 20 nm), each PE chain can be placed into the simulation box so that no boundary chain folding occurs. On the other hand, when the length of PE chains is longer than the simulation domain, they have to be folded at the boundary in order to t into the simulation box. The two-segment kinked temperature prole obviously suggests that the boundary chain folding has a remarkable impact on the thermal transport in LAPE.
In order to quantify the effect of boundary chain folding, we calculate the thermal conductivity of LAPE with varying average length L 0 . When a uniform linear temperature gradient (square in Fig. 2) is established, thermal conductivity can be directly calculated according to eqn (1) for the whole region. In the presence of boundary chain folding, however, the two-segment piecewise linear temperature gradient (circle in Fig. 2) occurs, and thus only local thermal conductivity in each linear segment can be calculated by eqn (1) instead. As shown in Fig. 2, k 1 and k 2 represents the local thermal conductivity for the boundary and internal region, respectively.
Here, a natural question to be asked is how to dene the effective thermal conductivity for the whole system in the presence of non-uniform temperature gradient? This problem has been seldomly emphasized in literatures as most materials are homogeneous along the temperature gradient direction. For instance, one widely used experimental technique to measure thermal conductivity is the thermal bridge method, 34,39,40 which directly relates thermal conductivity of sample to the temperature change in the heater and sensor assuming a linear temperature gradient across the sample. Since the non-uniform temperature gradient is not considered in the experimental model, the second question is the validity of such calculation (k 0 in the inset of Fig. 2) and underlying assumption.
To address these questions, let us start with the relationship between the thermal resistance (R) and thermal conductivity (k) k ¼ L/(SR), where L and S denotes the length and cross sectional area for the system of interest, respectively. Then, one can dene the local thermal resistance as: where L 1 and L 2 denotes the length of boundary and internal segment with linear temperature gradient, respectively. Similar to the electric circuit case, one has to rely on the serial or parallel connection model in order to get the total resistance. As the boundary segment and internal segment are placed adjacently along the temperature difference direction, it is reasonable to assume serial heat conduction channels rather than parallel connection. Therefore, the effective (total) thermal resistance R eff is given by the serial connection rule as R eff ¼ R 1 + R 2 . Correspondingly, the effective thermal conductivity k eff for the whole system can be expressed as: Based on eqn (4), one immediate derivation is that k eff z k 2 in the limit when L 1 ( L 2 . In other words, when the length of the boundary segment is much smaller than that of the internal segment, the assumption of a uniform linear temperature gradient is well justied as the temperature gradient for the internal segment is equivalent to the temperature gradient across the whole system in this limit case. This condition can be satised in most experimental measurements when the contact (boundary) region is much smaller than the internal suspended region.
Interestingly, we further prove rigorously in ESI † that k eff given by eqn (4) is actually identical to k 0 in which one ctitious temperature gradient is used in the calculation (inset in Fig. 2), regardless of the emergence of the non-uniform temperature gradient. That is to say, even in the presence of non-uniform temperature gradient among different segments, if one assumes serial connection between each heat dissipation segment, it is essentially equivalent to dene an effective uniform temperature gradient established across the whole system. Similarly, for thermal conductance, we can also dene an effective thermal conductance analogue to the effective thermal conductivity (see ESI † for more detail).
To verify such equivalence numerically, we compute two types of effective thermal conductivity k eff (eqn (4) with two linear ttings) and k 0 (eqn (1) with one linear tting) for different LAPE samples with varying length (Fig. 3). For short length samples (L 0 < 20 nm), only k 0 can be dened as there is only one uniform temperature gradient. Besides, thermal conductivity exhibits an increasing trend with the average chain length, which is consistent with the previous theoretical study on amorphous polymer at short length scale. 41 For longer samples (L 0 $ 20 nm) where boundary chain folding takes place, k eff and k 0 overlap with each other within the statistical error, and both of them show no obvious length dependence. Moreover, the local thermal conductivity for the internal region (k 2 ) is consistently higher than that for the boundary region (k 1 ) for the same L 0 . Obviously, the emergence of non-uniform local thermal conductivity hinders the continuous growth of effective thermal conductivity with length in LAPE. These results suggest that the boundary chain folding has a remarkable impact on thermal conductivity of LAPE.

Effect of chain length distribution
In order to explore the effect of chain length distribution on the thermal transport in LAPE, we further compute thermal conductivity k 0 of PE chains with varying s for two typical length L 0 . As shown in Fig. 4(a), when L 0 is equal to 10 nm, k 0 is almost constant ($0.4 W m À1 K À1 ) when varying s from 0 to 3 nm. For L 0 ¼ 15 nm, k 0 is higher ($0.7 W m À1 K À1 ) and remains constant when increasing s up to 2 nm, but drops to a lower value when s ¼ 3 nm, in which the boundary chain folding process takes place. By computed the local thermal conductivity for L 0 ¼ 15 nm and s ¼ 3 nm, we nd the thermal conductivity for the internal region k 2 (empty triangle in Fig. 4(a)) is still close to that before the emergence of chain folding, but thermal conductivity for the boundary region k 1 (solid triangle in Fig. 4(a)) is much lower (only $0.2 W m À1 K À1 ), leading to the reduced effective thermal conductivity according to eqn (4). Therefore, we can conclude that the variation of s value has no impact on thermal conductivity of LAPE unless the boundary chain folding occurs.
The variation of thermal conductivity with s can be understood by analyzing the distribution function of PE chain length. As shown in Fig. 4(b), the upper bound for chain distribution for the sample with L 0 ¼ 10 nm (solid symbols) does not exceed the length of simulation domain (L ¼ 20 nm) for all s values, which means that each chain has enough space to extend. As we x the average length L 0 , with the variation of s value, the chain length of some PE chains decreases while that of others increases, so that the thermal conductivity for the whole system is size independent. This is also the case for the sample with L 0 ¼ 10 nm and s < 3 nm. However, for the sample with L 0 ¼ 15 nm and s ¼ 3 nm, there is a notable portion of PE chains whose chain length is longer than the length of simulation domain. Consequently, these long PE chains have to be folded at the boundary in order to t into the simulation box, which gives rise to the non-uniform local thermal conductivity.

The inuence of morphology
As the mass density is uniform along the length direction of LAPE, the emergence of non-uniform local thermal conductivity suggests that the boundary chain folding induced morphology change could be an important factor. To directly visualize the chain morphology in these samples, we show in Fig. 5(a) a typical snapshot of the molecular conformation for LAPE with a long length (L 0 ¼ 30 nm and s ¼ 1 nm). Here three individual chains are highlighted with different colors. It is clearly shown that each PE chain is more extended in the internal region (region B), while it becomes more twisted in the boundary region (region A).
In order to quantify the morphology change, we compute the radius of gyration (R g ) of PE chains, which is one of the most important parameters to characterize the conformation of amorphous PE. The radius of gyration for an individual PE chain is dened as: 24,42 where M is the total mass, r cm is the center of mass position, and m i and r i are the mass and position for atom i, respectively. For a single polymer chain, the summations in eqn (5)- (7) include all atoms on the same chain. For a polymer chain with the xed length, a larger R g value means more extended chain morphology. As the chain folding phenomenon mostly takes place in the boundary region, in order to capture the morphology change in a specic region, we further compute the local radius of gyration R g l , which is obtained by restricting the summations in eqn (5)- (7) only to the atoms on the same chain inside a specic region and then averaging over all polymer chains in the same region.
Since the morphology of a PE chain is mainly determined by carbon atoms on the backbone, 43,44 we only include the carbon backbone in order to simplify the calculation of R g l . In order to record the structure variation, we run EMD simulation for another 1 ns aer the structure relaxation process. Fig. 5(b) shows the variation of R g l versus time for the sample with L 0 ¼ 30 nm and s ¼ 1 nm in different regions. For each region, R g l is almost time-independent, conrming that the polymer  structure has been fully relaxed. Besides, R g l in region A is consistently lower than that in region B. This difference means PE chains near the boundary are more entangled as a result of the chain folding, which is consistent with the result that the local thermal conductivity for the boundary region is lower than that for the internal region (Fig. 3). Such relation suggests that the entangled PE chains suppress the thermal transport in LAPE. Previous studies also have shown that more extended chain morphology with larger radius of gyration is favorable for more efficient thermal transport in amorphous polymers. 45 Although the radius of gyration can measure the degree of entanglement for polymer chains, it contains no information about the chain alignment direction, which is an important factor that can affect the polymer's thermal conductivity. 46 To this end, we introduce here the orientational order parameter P 2 that can be used to characterize the segment alignment along different directions dened as: 47-49 where e x is the unit vector in the direction of the X-axis, Y-axis, or Z-axis, the angular bracket means the ensemble average, and e i is the direction vector for a line segment connecting two nearest neighbors on the same chain given by A value of P 2 ¼ À0.5 means that the orientation is perpendicular to the selected direction, and P 2 ¼ 1 means that the orientation is parallel to the selected direction. Similar to the calculation of R g l , we compute the local orientational order parameter P 2 l , including only those atoms in a specic region. As shown in Fig. 5(c), P 2 l in region B exhibits the largest value along the X-axis (temperature bias direction), but P 2 l in region A has a peak value along the Y-axis. In particular, along the X-axis, P 2 l in region B is larger than that in region A. These results indicate that PE chains in the internal region are mostly aligned along the length (temperature bias) direction, while the alignment of PE chains in the boundary region deviates substantially from the thermal transport direction due to the chain folding process. As the heat conduction in polymer is mainly through the covalently bonded backbone, 37 the conformable alignment of polymer chains along the temperature bias direction will greatly promote the thermal transport in polymer systems. This explains the origin for the non-uniform local thermal conductivity in the presence of chain folding, and the lower thermal conductivity in the boundary region. Our results are consistent with the previous study which reported the improvement on the chain alignment, achieved by the mechanical stretching, can increase the polymer's thermal conductivity. 46 Therefore, in addition to the morphology change, the chain alignment is another factor responsible for the higher thermal conductivity in the internal region.

Vibrational spectrum analysis
In dielectric solids, thermal energy is mainly carried by the atomic vibrations. In crystalline materials, the atomic vibration can be described by the phonons. The reduced phonon lifetime is one of the mechanisms responsible for the reduction of thermal conductivity in crystalline solids, 50,51 due to the enhanced phonon scatterings. However, for amorphous materials, such as the amorphous polyethylene in our work, the phonon properties (such as the group velocity and lifetime) are no longer well-dened due to the lack of lattice periodicity. 52 Vibrational spectral analysis is an important tool to understand the thermal transport mechanism in various nanostructures. [53][54][55] However, few studies have analyzed the inuence of the chain morphology on vibrational properties of polymer. To probe the underlying mechanisms for the effect of chain folding on thermal conductivity, we compute the vibrational density of states (VDOS) as 56-58 where v j (t) is the time-dependent velocity for atom j, and N is the total number of atoms in a specic region. Fig. 6(a) shows the VDOS in different regions (as shown in Fig. 5(a)) for the sample with L 0 ¼ 30 nm and s ¼ 1 nm. For the extended PE chains in region B, VDOS exhibits the high frequency peak around 90 THz, which is a characteristic for CH 2 stretching vibration. 59 Moreover, the VDOS peak at $50 THz is the characteristic G-peak of carbon material, 60 and the lowfrequency peaks (<30 THz) originates from the exibility (angular bending and dihedral rotations) of the carbon-carbon bonds. 56 Compared to the carbon-carbon bond, carbonhydrogen bonds are associated with higher vibrational frequencies due to the lighter weight of hydrogen atoms. 56 As a result of the chain folding in region A, the high frequency spectrum ($50 THz) shis to the lower frequency, and the ultrahigh frequency vibrations ($90 THz) are signicantly suppressed in the PE chains. Such frequency shi is qualitatively consistent with the previous molecular dynamics study by Brayton et al., 61 in which they found that the vibration spectrum shied to the lower frequency when the structure of PE changed from crystal to semi-crystal to amorphous structure. In our study, the PE chains in region A are less extended and more conned compared to those chains in region B. The chain folding restricts the vibration of the covalent C-C bond and reduces the order of chain alignment.
In order to further understand the effect of chain folding on thermal transport properties of LAPE, we calculate the frequency spectrum of transmission coefficient G(u) dened as 62,63 GðuÞ ¼ qðuÞ where k B is the Boltzmann constant, DT is the temperature difference, and q(u) is the spectral heat ux across the imaginary interface (Fig. 5(a)). The heat ux can be computed from NEMD simulations according to 62,64 qðuÞ where A is interface cross sectional area, v i is the velocity of atom i, F ij is the force between the atom i and j, and the symbol L and R denotes the le and right side of the imaginary interface, respectively. Fig. 6(b) shows the calculation results of G(u) from NEMD simulation. Compared to the transmission spectrum for region B (solid line), the original characteristic peaks $30 THz and 50 THz have been greatly suppressed in region A (dashed line). In particular, the transmission of C-H vibration at very high frequency ($90 THz) is also forbidden in region A. Therefore, the boundary chain folding induces substantial suppression of the transmission for both C-C vibration and C-H vibration, which leads to the reduced local thermal conductivity in the boundary region.

Conclusion
In summary, we have studied the thermal transport behavior of the lamellar amorphous polyethylene with varying chain lengths by using molecular dynamics simulations. For short chain length, a uniform temperature gradient can be established and thermal conductivity increases with the average chain length. With the increase of average chain length, the phenomenon of chain folding starts to occur at the boundary, leading to the non-uniform temperature gradient (thus local thermal conductivity) along the heat transport direction. The local thermal conductivity for the boundary region is found to be notably lower than that for the internal region, and the effective thermal conductivity of the whole system shows no obvious length dependence. Compared to the boundary chain folding region, we nd that, based on the analysis of the chain morphology, the polyethylene chains in the internal regions are more extended and better aligned along the heat transport direction. Furthermore, the vibrational spectrum analysis reveals that the boundary chain folding shis the vibrational spectrum to lower frequency, and suppresses the transmission coefficient for both C-C vibration and C-H vibration, leading to the reduced local thermal conductivity in the boundary region. Our study provides a microscopic insight to the understanding of fundamental thermal transport behavior in lamellar amorphous polymers.

Conflicts of interest
There are no conicts of interest to declare.