Recent research has highlighted the potential to achieve high-thermal-conductivity polymers by aligning their molecular chains. Combined with other merits, such as low-cost, corrosion resistance, and light weight, such polymers are attractive for heat transfer applications. Due to their quasi-one-dimensional structural nature, the understanding on the thermal transport in those ultra-drawn semicrystalline polymer fibers or films is still lacking. In this paper, we built the ideal repeating units of semicrystalline polyethylene and studied their dependence of thermal conductivity on different crystallinity and interlamellar topology using the molecular dynamics simulations. We found that the conventional models, such as the Choy-Young's model, the series model, and Takayanagi's model, cannot accurately predict the thermal conductivity of the quasi-one-dimensional semicrystalline polyethylene. A modified Takayanagi's model was proposed to explain the dependence of thermal conductivity on the bridge number at intermediate and high crystallinity. We also analyzed the heat transfer pathways and demonstrated the substantial role of interlamellar bridges in the thermal transport in the semicrystalline polyethylene. Our work could contribute to the understanding of the structure–property relationship in semicrystalline polymers and shed some light on the development of plastic heat sinks and thermal management in flexible electronics.

## I. INTRODUCTION

Polymers have been widely used in a variety of energy-related applications, such as solar cells,^{1,2} thermoelectric devices,^{3} and batteries^{4} because they are lightweight, corrosion-resistive, stable, nontoxic, and inexpensive. However, the low thermal conductivity of pure amorphous polymers, which is typically $\u223c0.1\u20130.5\u2009W/(m\u2009K)$,^{5–7} to a large extent limits their applications in the heat dissipation area. Efforts to improve the thermal conductivity of polymers-based materials can be classified into two categories: (1) blending the polymer with materials with higher thermal conductivity, such as manufacturing polymer composites with hybrid fillers,^{8,9} or metal powders;^{10} (2) mechanically aligning the polymer chains and taking advantage of the high thermal conductivity along the molecular chain backbone.^{11–13} Various theoretical calculations and simulations have proved that the thermal conductivity of a single polymer chain or a crystalline polymer unit is considerably high along the chain direction, which can even exceed 100 W/(m K).^{12,14} In this spirit, experimentally ultra-drawn polyethylene (PE) nanofibers were reported to have as extremely high thermal conductivity as ∼100 W/(m K); thermal conductivity of the commercial ultra-drawn single PE microfibers (e.g., spectra 2000) has been measured as ∼10–16 W/(m K).^{15–17} Zhu *et al.* showed that the thermal conductivity of those commercial microfibers could increase ∼80%–150% with a moderate draw ratio 2–6.^{16} Because fibers provide limited practical benefits, researchers start to propose methods to massively fabricate large-area polyethylene films with high thermal conductivity. Ghasemi *et al.* proposed to manufacture ultra-drawn highly aligned PE polymer films through a sol-gel and ultra-drawn process.^{18} The measured thermal conductivity of these preliminary high-crystallinity highly aligned PE polymer films can be comparable to a single commercial fiber, $16\u2009W/(m\u2009K)$. Recently, an even higher thermal conductivity, ∼60 W/(m K), was achieved in scalably fabricated PE films by a roll-to-roll system.^{19} In all of these nanofiber, microfiber, and film samples, the PE is in a semicrystalline state with high structural anisotropy while the thermal transport in those materials is not well understood.

Semicrystalline polymer is a two-phase mixture system with the ordered crystalline lamellar phase and the disordered amorphous phase. Isotropic semicrystalline polymer is generally treated as a spherulite with the folded chains forming an isotropic spheroidal structure; anisotropic semicrystalline polymer can be considered as an alternative stack of the two phases to simplify the physical models. An ideal repeating unit of semicrystalline polymers is a sandwich structure with two crystalline regions at the sides called lamella and the amorphous region in the middle called interlamellar region.^{20,21} Theoretically, Choy^{22} proposed a thermal conductivity model based on the Maxwell effective medium equation to predict how the thermal conductivity of semicrystalline polymers depends on crystallinity. For highly oriented semicrystalline polymers with an orientation angle $\theta $ between the lamellar axis and the drawn direction, Choy-Young's model is

where $Ke$, $Ka$, $Kc\u2225$, $Kc\u22a5$, and *X* are the thermal conductivity of the semicrystalline polymer, the amorphous counterpart, crystalline counterpart (along the chain), the crystalline counterpart (perpendicular to the chain), and crystallinity, respectively. This model fits well with the experimental data of isotropically distributed lamellas;^{22} nevertheless, it fails in highly drawn semicrystalline polymers. Besides, the series model^{23,24} is widely used to apply the thermal resistance network concept in such a two-phase mixture system. However, it underestimates the thermal conductivity of semicrystalline polymers at high crystallinity as well. The reason is that both models treat the entire interlamellar region as one amorphous region and underestimates the role of the chain topologies. Specifically, interlamellar bridges that connect two adjacent crystalline lamellas increase their fractions and change their chain alignment dramatically with stretching, which play an important role in the heat conduction. Takayanagi^{25} proposed a model that specifically treats those bridges the same as the crystalline stems, which are aligned in the extrusion direction and have the same thermal conductivity $\u2009Kc$. The thermal conductivity of the semicrystalline region depends on a fitting parameter *b*, which was roughly assumed to be the width ratio of the region formed by these highly oriented bridges (defined as the ratio of the width of the bridge region to the width of the whole semicrystalline region). However, the Takayanagi model fails to predict the temperature dependence of thermal conductivity due to overestimating the contribution from bridges. Similar to the Takayanagi's model, Halpin-Tsai equation^{26} was adopted based on the Cox model^{27} that considered the bridges as a fiber phase connecting at least two lamellas and forming a continuous crystal region, while the thermal conductivity was fitted with the aspect ratio of the fibers (defined as the ratio of the diameter to the length of the fibers) as a free parameter.^{28} Both Takayanagi's and Cox's models treat bridges as tie chains (straight ultra-aligned chains), which could be true at ultrahigh draw ratios, however, would lose their predictive power at low and intermediate draw ratios when bridges are in a much coiled state. Moreover, all the models above could not reveal the heat transfer pathways among different chain topologies in the interlamellar region.

Recently, the enhanced Monte Carlo (EMC)^{21,29} method developed by in't Veld *et al.* was widely adopted to simulate the mechanical and structural properties of semicrystalline polymers, especially under tensile deformation.^{30,31} This method is a powerful tool that performs random Monte Carlo (MC) moves to mimic the formation of interlamellar region of semicrystalline systems and allows the investigation into different interlamellar topology. In this work, we built the ideal one-dimensional repeating unit of semicrystalline PE systems with different crystallinities and interlamellar topology using the EMC method. Their thermal conductivities were calculated by the nonequilibrium molecular dynamics (NEMD) simulations. Results showed that besides the dependence on crystallinity, thermal conductivity highly depends on bridge numbers at intermediate and high crystallinity. Our study explained why the Choy-Young's model and series model lost its predictive power at high crystallinity. Moreover, we modified the Takayanagi's model by incorporating the thermal conductivity of bridges and analyzing the heat transfer pathways along bridges, loops, and tails in the interlamellar region, providing a guide on evaluating the thermal conductivity of high-crystallinity polymers.

## II. MODEL AND METHODS

### A. Semicrystalline polyethylene model construction by the Enhanced Monte Carlo (EMC) method

Figure 1 shows the schematics of a typical semicrystalline polyethylene model in our simulations. We chose the minimum representative repeating unit, including two crystalline regions and an interlamellar region, as our model. We used the crystalline PE unit cell to represent the crystalline region, which is a different but simplified version from the switchboard model of a polymer crystalline lamella. After quenching, the interlamellar region evolves into an amorphous region and a transitional region (i.e., inter-phase regions). Typically, we take into account three types of chain topologies in the interlamellar region, i.e., bridges, loops, and tails. Among these, bridge is defined as the segments connecting the crystalline stems throughout the interlamellar region; loop is defined to connect crystalline stems at the same side; and for tail, one side of it is connected to the crystalline stems while the other side terminates.

The EMC method is adopted to build those models following the procedures described in Refs. 30–33. We briefly describe the major procedures as follows: (1) initially, a crystalline PE structure comprising $4\xd76\xd776$ unit cells was generated with the lattice constant $a=7.72\u2009\u01fa,\u2009b=4.45\u2009\u01fa,\u2009c=2.53\u2009\u01fa$, resulting in 7296 atoms and 48 chains in the simulation box. We used the united-atom (UA) force field^{34} (see Table III in Appendix A for parameters in our simulation) to describe the interatomic potential. The chains were then tilted to make sure that the {201} plane is normal to the *z*-direction, following the procedure in Refs. 21, 31, and 35, because such an orientation in semicrystalline PE has been suggested as one of the most common configurations due to its minimum interfacial energy.^{36,37} We selected and fixed 10, 15, 24, 30, 33, and 38 layers of atoms on each side in the z-direction as the crystalline region in 0%-, 18%-, 50%-, 72%-, 83%-, and 100%-effective crystallinity configurations, respectively. Each layer contains $4\xd76\xd71$ unit cells. Here, the effective crystallinity $X$ is defined as $X\u2009=\u2009Lc/Ls$, with $Lc$ and $Ls$ as the effective length of crystalline and the whole semicrystalline regions used for fitting in the following NEMD thermal conductivity calculations, respectively. We note here that even for the 0%-crystallinity configurations, there are still 20 layers crystalline region on each side because we need to reserve the heat source and heat sink regions for NEMD thermal conductivity calculations where the linear fitting region is entirely amorphous. (2) Next, 16 of the 48 chains were cut randomly and 1201, 1000, 600, 343, and 214 united atoms in these chains were removed in 0%-, 18%-, 50%-, 72%-, and 83%-effective crystallinity configurations, respectively, to match the experimental amorphous density of $\u223c0.855\u2009g/cc\u2009(Ref.\u2009$38) at $T\u2009=300\u2009K$ and $P=1\u2009atm$. This procedure generated 32 tails and 32 bridges in the interlamellar region. (3) After that, all the united atoms in the interlamellar region underwent random Monte Carlo (MC) moves, including single-site displacement, end-rotation, re-bridging, end-reptation, and end-bridging, to form the interphase and amorphous regions. Among these moves, bond formation and breakage were involved. Specifically, single-site displacement displaces a random site slightly to change the bond lengths and angles; end-rotation rotates the sites at a tail end with one to three sites involved in the rotation; re-bridging cuts three sites, which are consecutively bonded, from a bridge and then reinserts them in the bridge by searching for the possible insertion positions without altering the lengths or angles; end-reptation excises and appends a segment of one tail end to another tail end with random torsion angles, forming new tails; end-bridging move introduces the conversion between bridge-tail pairs and loop-tail pairs. During this procedure, the number of tails and the total number of bridges and loops remained the same. Detailed information can be found in Refs. 21 and 31. These random MC moves were performed at 10 000 K in order to promote a rapid randomization of the interlamellar configuration, generating the amorphous phase.^{31} (4) Finally, these configurations were equilibrated at 10 000 K for 20 000 Monte Carlo cycles, followed by quenching from 10000 K to 300 K stepwise during 600 000 MC cycles. The quenching was performed in the sequence of 5000 K, 2000 K, 1000 K, 750 K, 500 K, 400 K, 350 K, and 300 K in 20 000, 40 000, 40 000, 40 000, 60 000, 80 000, 160 000, and 160 000 MC cycles, respectively. A merit in building configurations by the EMC method is that the atoms in the interphase region are covalently bonded to the crystalline stems and amorphous regions instead of simply stacked through the van der Waals interactions. We prepared 8 to 10 independent initial configurations for each crystallinity system. A configuration with *X *=* *50% built by the EMC method is shown in Fig. 2(a) for illustration.

### B. Molecular dynamics (MD) simulations

MD simulations using the LAMMPS package^{39} were performed to equilibrate the ensembles generated by the EMC methodology and to study their thermal properties. The same UA force field and periodic boundary conditions in all three directions were applied. Time step equals 1 fs in all simulations because there are no fast-vibrating hydrogens in the model. Independent configurations for each crystallinity were equilibrated in the NPT (constant total number of particles N, pressure P, and temperature T) ensemble for 8 ns, during which the pressure *P *=* *1 atm and temperature *T *=* *300 K were maintained. The pressure and temperature damping parameters were 1000 fs and 100 fs, respectively. After 8 ns of NPT equilibration, the fluctuations of the system volume are within 0.3% and energy converges within 0.7%. Then, 0.2 ns NVE ensemble (constant total number of particles N, the volume V, and the total energy E of the system) was performed for relaxation. During MD simulations, the numbers of bridges, loops, and tails do not change.

### C. Configuration analysis

We summarized the number of bridges, loops, and tails and the simulation box parameters of each configuration after MD equilibrations in Table I. The length in the *z*-direction in different systems was roughly the same ($\u223c163\u2009\u01fa$). The fitting length $Ls$ for the effective thermal conductivity is fixed as $120\u2009\u01fa$, and thus, the lengths of crystalline (amorphous) region for fitting $Lc$ ($Li$) are 0 (120), 22 (98), 60 (60), 86.5 (33.5), and 99.2 (20.8) $\u01fa$ at 0%-, 18%-, 50%-, 72%- and 83%-crystallinity, correspondingly. These crystalline and interlamellar (inter-phase and amorphous) regions were distinguished by the density profiles, with details shown in the Appendix B. The three interlamellar topologies (bridges, loops, and tails) of independent configurations are also summarized in Table I. Note that the absolute numbers of bridges and loops depend on the random MC moves in EMC and they will reach an equilibrium number after sufficient long simulation time. We listed the numbers of bridges and loops for each crystallinity. We also note that some of the independent configurations at the same crystallinity have the same bridge number; however, the alignment and the atom number of each chain are different due to the random MC moves. Among all the configurations, the tail number is fixed as 32, and the average number of bridge (loop) is 2 (30), 4.2 (27.8), 11.7 (20.3), and 20.6 (11.4) at 18%-, 50%-, 72%- and 83%-crystallinity, respectively. Obviously, the length of interlamellar region decreases with increasing crystallinity, and thus, more bridges tend to be formed in higher crystallinity systems, which is at the cost of loop number.

Crystallinity . | X = 0%
. | X = 18%
. | X = 50%
. | X = 72%
. | X = 83%
. | X = 100%
. |
---|---|---|---|---|---|---|

Bridge No. | … | 1, 2, 4 | 2, 3, 4, 6, 7 | 9, 10, 11, 12, 14, 15 | 16, 18, 19, 20, 21, 22, 23, 27 | … |

Loop No. | … | 31, 30, 28 | 30, 29, 28, 26, 25 | 23, 22, 21, 20, 18, 17 | 16, 14, 13, 12, 11, 10, 9, 5 | … |

Tail No. | … | 32 | 32 | 32 | 32 | … |

$Lx\u2009(\u01fa)$ | $33.9\xb10.9$ | $34.1\xb10.9$ | $35.2\xb10.4$ | $35.9\xb10.5$ | $36.3\xb10.5$ | 38.5 |

$Ly\u2009(\u01fa)$ | $27.8\xb10.3$ | $27.8\xb10.2$ | $27.8\xb10.1$ | $28.1\xb10.2$ | $28.3\xb10.2$ | 27.1 |

$Lz\u2009(\u01fa)$ | $165.7\xb13.5$ | $166.5\xb14.3$ | $165.8\xb11.7$ | $163.0\xb13.0$ | $161.5\xb11.6$ | 149.4 |

$Lc\u2009(\u01fa)$ | 0.0 | 22.0 | 60.0 | 86.5 | 99.2 | 120.0 |

$Li\u2009(\u01fa)$ | 120.0 | 98.0 | 60.0 | 33.5 | 20.8 | 0.0 |

Crystallinity . | X = 0%
. | X = 18%
. | X = 50%
. | X = 72%
. | X = 83%
. | X = 100%
. |
---|---|---|---|---|---|---|

Bridge No. | … | 1, 2, 4 | 2, 3, 4, 6, 7 | 9, 10, 11, 12, 14, 15 | 16, 18, 19, 20, 21, 22, 23, 27 | … |

Loop No. | … | 31, 30, 28 | 30, 29, 28, 26, 25 | 23, 22, 21, 20, 18, 17 | 16, 14, 13, 12, 11, 10, 9, 5 | … |

Tail No. | … | 32 | 32 | 32 | 32 | … |

$Lx\u2009(\u01fa)$ | $33.9\xb10.9$ | $34.1\xb10.9$ | $35.2\xb10.4$ | $35.9\xb10.5$ | $36.3\xb10.5$ | 38.5 |

$Ly\u2009(\u01fa)$ | $27.8\xb10.3$ | $27.8\xb10.2$ | $27.8\xb10.1$ | $28.1\xb10.2$ | $28.3\xb10.2$ | 27.1 |

$Lz\u2009(\u01fa)$ | $165.7\xb13.5$ | $166.5\xb14.3$ | $165.8\xb11.7$ | $163.0\xb13.0$ | $161.5\xb11.6$ | 149.4 |

$Lc\u2009(\u01fa)$ | 0.0 | 22.0 | 60.0 | 86.5 | 99.2 | 120.0 |

$Li\u2009(\u01fa)$ | 120.0 | 98.0 | 60.0 | 33.5 | 20.8 | 0.0 |

### D. Thermal conductivity calculations

The thermal conductivity was calculated by the NEMD simulations.^{40,41} In order to establish a temperature gradient along the z-direction without introducing the coupling between thermostats, vacuum regions were added on each side, and the atoms in the boundary regions along the z-direction were fixed. Then, the Langevin thermostats were applied to the heat source and sink regions as shown in Fig. 2(a).

The non-equilibrium simulations were performed for 5 ns to make sure that the system is in steady sate. In steady state, the heat flux $Jz$ along the $z$-direction is

Here, $A$ is the cross-sectional area in the $x-y$ plane, $\Delta t$ is the time step, $\Delta Esource$ and $\Delta Esink$ are the amounts of heat added to the heat source and subtracted from the heat sink, respectively. Results show that $\Delta Esource\Delta t$ and $\Delta Esink\Delta t$ are the same, confirming the symmetry of our system.

The effective thermal conductivity $Ke$ of the whole semicrystalline region is defined as

where $\Delta T$ is the temperature difference across the fitted region, $\Delta L$ is the length of the whole fitted region. All the parameters were calculated in the z-direction. Regions used to fit the thermal conductivity are far away enough from the heat bath regions, as shown in Fig. 2(b). The amorphous thermal conductivity $Ka$ was fitted in the 0%-crystallinity configurations, as shown in Fig. 2(c), and the middle region was adopted to fit $Ka$.

## III. RESULTS AND DISCUSSION

### A. Thermal conductivity analysis

Figure 3(a) shows the dependence of the effective thermal conductivity of all the model polymers on their crystallinity *X*. Thermal conductivities of independent configurations with the same *X* were averaged statistically. To compare with the Choy-Young's model, we calculated the intrinsic thermal conductivity along the chain $Kc\u2225$ and perpendicular to the chain $Kc\u22a5$ of a crystalline PE system and plotted the prediction with the consideration of the tilted angle between the crystalline chain direction and the z-direction. To calculate $Kc\u2225$ and $Kc\u22a5$, we built configurations with the {001} plane normal to the z direction so that the chain direction is along the z-direction. Since the NEMD method has a strong size effect on calculating the thermal conductivity in the chain direction, we first calculated the thermal conductivity of this configuration with an effective length 12 nm in the chain direction, $Kc\u2225$ = 35 $\xb1$ $3\u2009$ W/(m-K), and we assume the intrinsic thermal conductivity along the chain direction $Kc\u2225$ to be larger than this value. We also found that the Choy-Young's model is not sensitive to $Kc\u2225$ when $Kc\u2225$ > 10 W/(m-K), so we used $Kc\u2225$ = 35 $\xb1$ $3\u2009$ W/(m-K) in the model. $Kc\u22a5$ was obtained using the NEMD method by applying the heat baths on the two sides in the *x*-direction with an effective length 10 nm. $Kc\u22a5$ = $0.10\xb10.02\u2009$ W/(m-K), which is comparable to that of the amorphous counterpart. For the series model, we used the calculated averaged thermal conductivity of 100%- and 0%- crystallinity systems as $Kc$ and $Ka$ and plot the series model predictions $Ke$ in the shaded area. Using the UA force field and NEMD, $Kc$ is $21\xb12\u2009$ W/(m-K), and $Ka$ is $0.12\xb10.02\u2009$ W/(m K). Those calculated thermal conductivity values agree reasonably well with other simulation results^{42,43} and are consistently lower than the experimental results^{15,20,44} due to the simplified UA force field^{45,46} and size effect.^{47} The UA force field has been used to study the thermal conductivity of amorphous polymers, such as PE^{48} and polyamide-6,6.^{49} Both simulations agree well with experimental values.^{5–7,50} Besides, the UA force field has also been used to simulate the thermal conductivity of single crystalline PE chains,^{51} which gave an underestimated value compared with those using AIREBO potential.^{14} As pointed out in Ref. 45, using the UA force field, the specific heat of single PE chain at low temperature could be overestimated, and an inaccurately predicted phonon dispersion could lead to an underestimation in the acoustic sound velocity at low reduced wave vectors. We note here that the absolute thermal conductivity values from our simulations are not critical as long as those values that we used to testify the models are consistently from the same force field and simulation method. The use of UA model enables us to study the thermal properties in a much larger system with the limited computational resources. For polymers with different crystallinity, even though a length-dependent thermal conductivity might apply for the crystalline region, both the Choy-Young's model and the series model are not sensitive to the accurate value of $Kc$ that is two orders of magnitude higher than $Ka$.

Our NEMD calculation results agree well with the Choy-Young's model and the series model in their general trends, i.e., the effective thermal conductivity increases with increasing crystallinity. Quantitatively, the Choy-Young's model provides an excellent estimation of the average value at *X *=* *18% and 50%, but loses the predictive power at intermediate and high crystallinity, i.e., 72% and 83%, especially at *X *=* *83%. The discrepancy between the calculated value and the maximum Choy-Young's model prediction increases from ∼35% at 72%-crystallinity to ∼68% at 83%-crystallinity. Our results confirm the failure of Choy-Young's model at high crystallinity in the highly anisotropic semicrystalline polymer. Similarly, using the series model, the discrepancy between the calculated value and the maximum prediction from the model increases from ∼10% at 18%-crystallinity to ∼72% at 83%-crystallinity. We attribute such a large deviation of both models to the increased amount of bridges so that a single thermal conductivity of the amorphous counterpart is not able to describe the thermal conductivity in the interlamellar region. This statement is supported by the calculated “local effective thermal conductivity,” which is defined as the ratio of the heat flux to the local temperature gradient (see Appendix C for the plots). The local effective thermal conductivity of the amorphous region is neither the thermal conductivity of the amorphous counterpart nor the crystalline counterpart.

Figure 3(b) shows the thermal conductivity of semicrystalline PE as a function of their bridge number *n*. We group the calculated results into three regions: (I) at a low crystallinity, for example, *X *=* *18%, their bridge numbers are mainly in the range of 1–4. In these configurations, the effective thermal conductivity of the entire amorphous regions can be roughly treated the same as $\u2009Ka$; (II) at an intermediate crystallinity, i.e., *X *=* *50% and 72%, and their bridge numbers are in the range of 3–15. The role of bridges and their interaction with atoms in loops and tails becomes important, and the effective thermal conductivity increases linearly with n, where predictions from the series model start to deviate; (III) at a high crystallinity (*X *=* *83%), the bridge number can even achieve 27, thus heat transfer along bridges dominates so that the effective thermal conductivity increases linearly with *n* in a much larger slope. Those results indicate that not only the crystallinity, but also the number of bridges should be considered to estimate the effective thermal conductivity of semicrystalline polymer systems with high crystallinity.

### B. Modified Takayanagi's model

Even though Takayanagi's model takes into consideration the role of bridges, the failure of this model lies in the treatment of interlamellar bridges the same as crystalline stems. In order to investigate the role of bridges at high crystallinity, we seek to modify the Takayanagi's model. Instead of treating bridges as a continuous region of crystalline stems, we treated them separated from crystalline region, with a thermal conductivity $Kb$, and we built the model based on a thermal resistance network concept. Figure 4(a) shows our modified model. In the interlamellar region, bridges are in parallel with the amorphous part formed by the rest of the chain topologies, i.e., loops and tails. And the interlamellar region is in series with the crystalline region. The modified effective thermal conductivity is

Here, *b* is the width ratio of the bridge region, which assumes to be *b *=* *n/N, where n is the bridge number, and N is the total chain number (In our model, N is 48).

The interfacial thermal resistance between the crystalline and interlamellar region was negligible at 300 K in our configurations because (1) the covalent bonding between the two regions is much stronger than the van der Waals interactions. The interfacial thermal conductance has been proved to increase dramatically with stronger bonding strength.^{52–54} (2) The thermal resistance of the interface region is much smaller compared with that in the crystalline region and amorphous region. Since thermal conductivity of the amorphous region is two orders of magnitude smaller than that of the crystalline region, the effective thermal conductivity is dominated by the interlamellar thermal conductivity. Therefore, the interfacial thermal resistance between the crystalline and interlamellar region can be reasonably neglected. This negligible interface thermal resistance can be seen from the slope of temperature profile at the interface region, which does not show a noticeable deviation from that of the crystalline region in the NEMD simulation [Fig. 2(b)], proving the estimation in Ref. 20.

### C. Contributions to heat flux

To further understand the transition from region I to III, we divided the atoms in the interlamellar region into three groups, which include: (1) the upper group, where tails and loops are connected to the upper crystalline chain stem (group $u$); (2) the lower group, where tails and loops are connected to the lower crystalline chain stems (group $l$); (3) bridges, which connect the upper and lower crystalline chain stems (group $b$). The contribution from each group to the total heat flux was then studied. The possible pathways of heat transfer among these groups due to non-bonded interactions are simplified in Fig. 4(b).

In LAMMPS, the heat flux vector $J$ in a specific group can be calculated by summing the contributions from all the atoms in this group^{39}

where $V$ is the volume of the specific group, *i* is the atom index, $ei$ is the energy, i.e., potential and kinetic energy, of atom *i*, vector $vi$ is the velocity of atom $i$, and $Si$ is the stress tensor for atom *i*.

Since heat flows from heat source to heat sink in our simulations and atoms in group $l$ are not bonded to other groups, we reasonably assume the heat obtained by this group l is mainly due to the non-bonded interactions with atoms in the other two groups. Therefore, to study the role of non-bonded interactions in our model, we excluded specific pathways between different groups, relaxed the system in 0.2 ns NVT ensemble at 300 K (temperature damping parameter is 100 fs) followed by 0.2 ns NVE ensemble, and calculated the contributions from each group to heat flux and the thermal conductance of the interlamellar region after performing 5 ns non-equilibrium simulations.

Seven representative configurations were considered as summarized in Table II: (i) 18%-crystallinity with 2 bridges (Conf. i), (ii) 50%-crystallinity with 3 bridges (Conf. ii), (iii) 50%-crystallinity with 7 bridges (Conf. iii), (iv) 72%-crystallinity with 12 bridges (Conf. iv), (v) 83%-crystallinity with 20 bridges (Conf. v), (vi) 83%-crystallinity with 23 bridges (Conf. vi), and (vii) 83%-crystallinity with 27 bridges (Conf. vii).

Configurations . | Conf. i . | Conf. ii . | Conf. iii . | Conf. iv . | Conf. v . | Conf. vi . | Conf. vii . |
---|---|---|---|---|---|---|---|

Crystallinity (%) | 18 | 50 | 50 | 72 | 83 | 83 | 83 |

Bridge No. | 2 | 3 | 7 | 12 | 20 | 23 | 27 |

Configurations . | Conf. i . | Conf. ii . | Conf. iii . | Conf. iv . | Conf. v . | Conf. vi . | Conf. vii . |
---|---|---|---|---|---|---|---|

Crystallinity (%) | 18 | 50 | 50 | 72 | 83 | 83 | 83 |

Bridge No. | 2 | 3 | 7 | 12 | 20 | 23 | 27 |

Figure 5 shows the comparison among the contributions from different groups to heat flux. When we excluded pathways P2 and P3, i.e., non-bonded interactions between groups u and l, and that between groups l and b, heat in group u can only be transferred by non-bonded interactions within group b. Figure 5(a) shows the corresponding contributions to the total heat flux from groups b and u in different configurations. When *n* = 2, *X *=* *18% (region I), the contribution from loops and tails is around 30%, demonstrating the importance of pathway P1. Therefore, the interlamellar region can be treated roughly as a pure amorphous region. Although group u is not bonded to group b, the large number of atoms within the non-bonded interaction cut-off distance ensures this pathway competitive with those along bridges, which is consistent with the results in Ref. 55. To prove this idea visually, snapshots with each chain individualized from Visual Molecule Dynamics (VMD)^{56} are shown in Fig. 6(a). In contrast, with the increase of n, contributions from bridges become dominant, especially when $n\u226520$, where contributions from bridges are even larger than 97%. Thus, the modified Takayanagi's model has to be introduced in these configurations.

Similarly, if we excluded the pathways P1 and P2—the non-bonded interactions between groups u and l, and that between groups u and b—heat from group l is totally transferred by non-bonded interactions with group b. As shown in Fig. 5(b), similar results were observed, proving the effect of the non-bonded interactions between groups u and b and that between groups u and l are almost the same. We also simulated the case when the pathways P1 and P3 (non-bonded interactions between groups u and b, and that between groups l and b) were excluded. In this case, heat from group u can be transferred to group l through non-bonded interactions, but there is no non-bonded interaction between group b and the other two groups. As shown in Fig. 5(c), when *n* is small, the contribution to heat flux from groups u and l is even larger than that from group b. However, the contribution from group b increases with n and becomes dominant at high *n* number. Especially, when $n\u226520$, the contribution from group b is larger than 96%. Such a dominant role of bridges at high crystallinity is consistent with the results in Ref. 48, where they found that the strong covalent bonds dominate in the heat transfer compared with the nonbonded interactions.

In the following, we focus on the role of bridges of configurations in regions II and III. When all the non-bonded pathways among groups l, u, and b were excluded, heat is mainly transferred along bridges, and therefore, the thermal conductivity of bridges can be roughly estimated. One thing to note here is, at *X *=* *50% and 72%, their bridges are much coiled, excluding non-bonded interactions could introduce a slight difference in the chain alignment, and thus our results only show a rough estimation of $Kb$ [see Fig. 5(d)]. Here, $Kb$ is the intrinsic thermal conductivity of group b by substituting A in Eq. (2) with the effective cross-sectional area ($n48A$). When *n* is in region II, $Kb$ fluctuates around 0.9 W/(m K), while for *n* in region III, $Kb$ increases linearly with a slope of 0.076. Thus, we calculated the effective thermal conductivity using our modified model with $Kb=0.9\xb10.1\u2009W/(m\u2009K)$ for region II and $Kb=\u22120.53+0.076n$ for region III.

We attribute such a transition of *K*_{b} from region II to region III to two reasons: (i) bridges at *X *=* *83% are much aligned compared with those at *X *=* *50% and 72%. Figure 6(b) shows the VMD snap shot. Namely, coiled bridges at *X *=* *50% and 72% offer more phonon scattering centers; the increased number of scattering centers and the increased bridge number affect $Kb$ in an opposite and competitive manner, resulting in $Kb$ fluctuating around a constant value; (ii) at *X *=* *82%, with the large increase of bridge number, the role of non-bonded interaction becomes negligible, and most of the heat is transferred along chains through bond interactions. VMD snapshot also shows the more aligned bridge chains at high *X*, which is expected since the bridge chains tend to align when the interlamellar region is short (high *X*) and more bridge chains (large *n*) are present there. As simulations^{57} have pointed out that thermal conductivity along polymer chains exponentially scales with the chain orientation, such a less-coiled configuration could greatly enhance the heat transfer capability.

Finally, we test our modified model by predicting the thermal conductivity of PE with a crystallinity *X *> 20%, which is plotted in Fig. 3(b) to directly compare with the NEMD calculated data. The thermal conductivity values used in the model were from separate simulations: $Ka=0.12\xb10.02\u2009W/(m-K)$, $Kc=21\xb12\u2009W/(m-K),\u2009Kb=0.9\xb10.1\u2009W/(m-K)$ for *X *=* *50% and 72%, $Kb$ varies with *n* as $Kb=\u22120.53+0.076n$ with $\xb110%\u2009$ uncertainty for *X *=* *83%. We found that the predictions from our modified model fit excellently well with the effective thermal conductivity calculated from NEMD, which proves that this modified model captures the heat transfer pathways in the semicrystalline PE and the important role of bridges in the interlamellar region. In all the previous experimental studies of semicrystalline polymers, due to the challenges of defining and measuring the parameter *b*, only crystallinity *X* is considered for the data fitting and explanation. Our study shows that at high crystallinity, such fittings could give large errors (e.g., ∼70%) and how to treat the interlamellar bridges is the key for the modeling.

## IV. CONCLUSION

In summary, we investigated the effective thermal conductivity of semicrystalline PE with crystallinity varying in a large range, from 18% to 83%. Our results clearly show that the predictions from Choy-Young's model and the series model for the thermal conductivity value of semicrystalline polyethylene fail, especially at high crystallinity. We attributed this deviation to the substantial role of bridges in carrying the heat by studying their contributions to the heat flux in the interlamellar region. We modify Takayanagi's model to take into account the effective thermal conductivity of bridges, and this modified model successfully predicts the effective thermal conductivity of semicrystalline PE at intermediate and high crystallinity. Our results prove that there is a transition from the intermediate crystallinity (*X *=* *50% and 72%) to high crystallinity (*X *=* *83%) due to the non-bonded interactions and the alignment of bridges, and suggest that both crystallinity and the bridge number in the interlamellar region need to be considered to explain the thermal conductivity of polymers with high crystallinity. As a preliminary study, the simple one-dimensional model in this study could contribute to the understanding of thermal transport in semicrystalline polymers. We expect a structure-property guide based on a more sophisticated model that can describe the true structure of high-crystallinity polymer fibers or films in the future work. Our work could shed some light on the development of plastic heat sinks and the thermal management in flexible electronics.

## ACKNOWLEDGMENTS

The authors would like to thank Gregory C. Rutledge and Raghavan Ranganathan for their helpful discussions on building the semicrystalline PE model. K.K., T.L., and J.L. acknowledge the financial support from the NC State Faculty Research and Professional Development Funds. T.L. also acknowledges the funding support from China Scholarship Council. X.L. acknowledges the National Natural Science Foundation of China (NSFC) Grant No. 51506062. J.Z. acknowledges the program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning Grant No. TP2014012. G.C. acknowledges DOE-BES Award No. DE-FG02-02ER45977. Computational resources were provided by the High Performance Computing Center at North Carolina State University.

### APPENDIX A: UA FORCE FIELD

The united atom (UA) Force field was used in the EMC modeling and MD simulations consistently. The parameters are listed in Table III.

Non-bond . | $En=4\u03f5\sigma r12\u2212\sigma r6,\u2009r<rc$ . | ||
---|---|---|---|

Type . | $\u03f5\u2009Kcal\u2009mole\u22121$ . | $\sigma \u2009\xc5$ . | $rc\u2009\xc5$ . |

CH_{2}-CH_{2} | 0.09339913958 | 4.009 | 10.0225 |

CH_{3}-CH_{3} | 0.226542543 | 4.009 | 10.0225 |

Non-bond . | $En=4\u03f5\sigma r12\u2212\sigma r6,\u2009r<rc$ . | ||
---|---|---|---|

Type . | $\u03f5\u2009Kcal\u2009mole\u22121$ . | $\sigma \u2009\xc5$ . | $rc\u2009\xc5$ . |

CH_{2}-CH_{2} | 0.09339913958 | 4.009 | 10.0225 |

CH_{3}-CH_{3} | 0.226542543 | 4.009 | 10.0225 |

Bond . | $Eb=Kr\u2212r02$ . | ||
---|---|---|---|

Type . | $K\u2009Kcal\xb7mole\u22121\xb7\xc5\u22122$ . | $r0\u2009\xc5$ . | . |

CH_{2}-CH_{2} | 449.5 | 1.53 | |

CH_{2}-CH_{3} | 449.5 | 1.53 | |

CH_{3}-CH_{3} | 449.5 | 1.53 |

Bond . | $Eb=Kr\u2212r02$ . | ||
---|---|---|---|

Type . | $K\u2009Kcal\xb7mole\u22121\xb7\xc5\u22122$ . | $r0\u2009\xc5$ . | . |

CH_{2}-CH_{2} | 449.5 | 1.53 | |

CH_{2}-CH_{3} | 449.5 | 1.53 | |

CH_{3}-CH_{3} | 449.5 | 1.53 |

Angle . | $Ea=K\theta \u2212\theta 02$ . | ||
---|---|---|---|

Type . | $K\u2009Kcal\xb7mole\u22121\xb7rad\u22122$ . | $\theta 0\u2009deg$ . | . |

CH_{2}-CH_{2}-CH_{2} | 60 | 112 | |

CH_{2}-CH_{2}-CH_{3} | 60 | 112 | |

CH_{3}-CH_{2}-CH_{3} | 60 | 112 |

Angle . | $Ea=K\theta \u2212\theta 02$ . | ||
---|---|---|---|

Type . | $K\u2009Kcal\xb7mole\u22121\xb7rad\u22122$ . | $\theta 0\u2009deg$ . | . |

CH_{2}-CH_{2}-CH_{2} | 60 | 112 | |

CH_{2}-CH_{2}-CH_{3} | 60 | 112 | |

CH_{3}-CH_{2}-CH_{3} | 60 | 112 |

Dihedral . | $Ed=12K11+cos\u2009\varphi +12K21\u2212cos\u20092\varphi $ $+12K31+cos\u20093\varphi $ . | ||
---|---|---|---|

Type . | $K1\u2009Kcal\xb7mole\u22121$ . | $K2\u2009Kcal\xb7mole\u22121$ . | $K3\u2009Kcal\xb7mole\u22121$ . |

CH_{2}-CH_{2}-CH_{2}-CH_{2} | 1.6 | −0.867 | 3.24 |

CH_{2}-CH_{2}-CH_{2}-CH_{3} | 1.6 | −0.867 | 3.24 |

CH_{3}-CH_{2}-CH_{2}-CH_{3} | 1.6 | −0.867 | 3.24 |

Dihedral . | $Ed=12K11+cos\u2009\varphi +12K21\u2212cos\u20092\varphi $ $+12K31+cos\u20093\varphi $ . | ||
---|---|---|---|

Type . | $K1\u2009Kcal\xb7mole\u22121$ . | $K2\u2009Kcal\xb7mole\u22121$ . | $K3\u2009Kcal\xb7mole\u22121$ . |

CH_{2}-CH_{2}-CH_{2}-CH_{2} | 1.6 | −0.867 | 3.24 |

CH_{2}-CH_{2}-CH_{2}-CH_{3} | 1.6 | −0.867 | 3.24 |

CH_{3}-CH_{2}-CH_{2}-CH_{3} | 1.6 | −0.867 | 3.24 |

### APPENDIX B: DENSITY PROFILE AND THE DEFINITION OF CRYSTALLINE, INTERPHASE, AND AMORPHOUS REGIONS

The regions of crystalline, interphase, and amorphous were distinguished according to the density profile illustrated in Fig. 7. The average density of the crystalline region is $\u223c1\u2009g/cc$ and that of the amorphous region is $\u223c0.87\u2009g/cc$, which are reasonable compared with the experimental results. The drop in the density from crystalline to amorphous occurs in a transition region which is defined as the interphase region. The thickness of the interphase is estimated as $\u223c1\u2009nm$.^{58}

### APPENDIX C: LOCAL EFFECTIVE THERMAL CONDUCTIVITY

The local effective thermal conductivity is defined as the ratio of heat flux to the local temperature gradient. The local effective thermal conductivity $KL\u2009$ of representative samples was fitted. Four adjacent data points on each side of the fitted point were adopted for fitting. While for the fitted point close to the interface, one adjacent data point on each side was used for fitting. Figure 8 shows the $KL$ in configurations with *X *=* *72% and 83%. The $KL$ in the interphase region behaves as a transitional region between the crystalline and the amorphous region. Moreover, comparing Figs. 8(a) and 8(b), the amorphous region has a higher $KL$ with higher crystallinity because of the increasing bridge number and better alignment.