Thermal conductivity of polymers: A simple matter where complexity matters
Debashish Mukherji





[]debashish.mukherji@ubc.ca
Quantum Matter Institute, University of British Columbia, Vancouver V6T 1Z4, Canada



§ ABSTRACT 
Thermal conductivity coefficient κ measures the ability of a material to conduct a heat current.
In particular, κ is an important property that often dictates the usefulness of a material over a wide range of environmental conditions. 
For example, while a low κ is desirable for the thermoelectric applications, a large κ is 
needed when a material is used under the high temperature conditions. These materials range from common crystals 
to commodity amorphous polymers. The latter is of particular importance because of their use
in designing light weight high performance functional materials. In this context, however, one of the major limitations of the amorphous polymers is 
their low κ, reaching a maximum value of about 0.4 W/Km that is 2–3 orders of magnitude smaller 
than the standard crystals. Moreover, when energy is predominantly transferred through the bonded connections, κ≥ 100 W/Km. 
Recently, extensive efforts have been devoted to attain a tunability in κ via macromolecular engineering. 
In this work, an overview of the recent results on the κ behavior in polymers and polymeric solids is presented. 
In particular, computational and theoretical results are discussed within the context of complimentary experiments. 
Future directions are also highlighted. 


Thermal conductivity of polymers: A simple matter where complexity matters
    Debashish Mukherji
    September 9, 2024
==========================================================================





§ INTRODUCTORY REMARKS



When a material is subjected to a non–equilibrium condition of temperature T, such that one end of the 
material is kept at an elevated temperature T_ high and a lower temperature T_ low is maintained at 
the other end, a heat current flows from the hot to the cold region. This is a direct consequence of the second law of thermodynamics and 
is quantified in terms of the heat flux vector j⃗. Here, the Fourier's law of heat diffusion
states that j∝( T_ high - T_ low)/ℓ across a sample of length ℓ.
The proportionality constant is thermal transport coefficient κ, which is a key material property that commonly 
dictates the usefulness of a material for a particular application <cit.>. 

Traditionally, heat flow in crystalline materials and in nano–structures have been of primary 
interest  <cit.>.
Due to the long range order in crystals, the phonon mean free path Λ are rather large and thus κ≥ 100 W/Km <cit.>. 
In the carbon–based materials, κ can even exceed 1000 K/Wm <cit.>.

A complete opposite class to the crystals is the amorphous solids, where Λ is  small, 
i.e, within the direct atom–to–atom contact. Here, κ≤ 2 W/Km and heat propagates via localized vibrations <cit.>. 

Within the class of amorphous solids, polymers are of particular importance because they usually provide a flexible platform for 
the design of advanced functional materials <cit.>.
Some examples include, but are not limited to, organic solar cells <cit.>, electronic packaging and/or heat sinking materials <cit.>, 
thermal switches <cit.> and thermoelectric applications <cit.>.
However, the typical κ values of the amorphous polymeric solids are further 5–10 times smaller <cit.> 
than the standard amorphous materials (such as amorphous silicon). This often hinders the usefulness of polymers under the high T conditions.

Most commonly known (non–conducting) commodity polymers can be categorized into the systems where non–bonded monomer–monomer 
interactions are either dictated by the van der Waals (vdW) forces or by the hydrogen bonds (H–bonds) <cit.>. 
Here, the interaction strength of vdW is about k_ BT at a temperature T = 300 K and the Boltzmann constant k_ B, 
while the strength of a H–bond is between 4–8k_ BT depending on the dielectric constant of the medium <cit.>.  

A few examples of the commodity polymers is shown in Fig. <ref>.
Note that these particular polymers are chosen because their experimental and simulation data is readily available.



Experiments have reported that κ≃ 0.1-0.2 W/Km for the vdW systems <cit.> 
and for the H–bonded polymers κ→ 0.4 W/Km <cit.>.

Table <ref> lists κ and the corresponding glass transition temperatures T_ g for the polymers in Fig. <ref>. 

It can be seen in Table <ref> that PMMA (a vdW–based polymer) has κ≃ 0.20 W/Km 
and its T_ g≃ 378 K, while PVA (a H–bonded system) has κ≃ 0.310 W/Km 
and a lower T_ g≃ 348 K <cit.>. These values indicate that T_ g and κ are not correlated, which is also visible across many polymeric systems <cit.>. Furthermore, a closer look suggests that– within a simple approximation, T_ g is directly related to the trans–to–gauche free energy barrier Δ E_ t-g (i.e., local fluctuations),
which is dictated by a delicate combination of the bonded, the angular and the dihedral interactions along a chain backbone. The higher the Δ E_ t-g, the larger the T_ g. Here, PMMA with a larger side group has a higher Δ E_ t-g than a PVA, hence the observed trend in T_ g. Within this simple discussion it becomes reasonably apparent that the exact T_ g is a completely irrelevant quantity within the context of κ in polymers.

Amorphous polymeric solids are a special case because even when their macroscopic κ values are very small, 
at the monomer level they have different rates of energy transfer. For example, energy can be transferred 
between the bonded monomers and that between the neighboring non–bonded monomers. 

In this context, a closer investigation of the polymer structures reveal that the carbon–carbon (C–C) covalent bond 
constitute the most common backbone of commodity polymers, see Fig. <ref>. 
Here, it is known that the stiffness of a C–C contact is E ≥ 250 GPa <cit.>, while E of a vdW or a H–bond 
system vary between 2–5 GPa <cit.>. Given that κ∝ E <cit.>, the energy transfer between 
the two non–bonded monomers (soft contacts) is significantly smaller than along a (stiff) bond <cit.>. A simple schematic of this scheme is shown in Fig. <ref>.

The thermal behavior in polymers are predominantly dictated by the non–bonded contacts, while the energy transfer along a bonded contact plays a lesser important role. This is particularly because a chain in a frozen configuration follows the random walk statistics <cit.>, i.e., when it is quenched to T ≪ T_ g 
from a melt configuration.

Within this picture, when heat flows along a chain contour, it experiences scattering due to the bends and the kinks 
along the path <cit.>. Energy also occasionally hops off to a non–bonded neighboring monomer. 
A combination of these two effects facilitates a knocking down of κ in the amorphous polymers <cit.>. Note also that the exact monomer structure, i.e., the side groups connected to a backbone play an additional (delicate) role <cit.>, which will be discussed at a later stage within this short overview.



Over the last 2–3 decades, extensive efforts have been devoted to study the heat flow in polymers using the experimental, theoretical, and computational approaches <cit.>. 
In particular, even when the polymers are a class of simple matter with their great potential in designing flexible materials, establishing a microscopic understanding in polymers is rather complex. Here, one of the grand challenges in this field is to attain a predictive tunability in κ (almost at will) using macromolecular engineering.
This requires a protocol that can properly account for a delicate balance between the bonded to the non–bonded interactions, the chain conformations, 
and their morphology. Motivated by the above, this manuscript aims to highlight the latest developments in the field of polymer thermal conductivity. 
For this purpose, comparative experimental and simulation results will be discussed to put forward the key concepts.



§ EFFECT OF BLENDING ON THE THERMAL CONDUCTIVITY OF POLYMERS


In a system consisting of one (linear) polymer component, κ is rather limited because of the 
restrictive monomer–level interactions <cit.>. To circumvent this problem, 
studies have suggested that κ of a polymeric solid may be enhanced by blending a second component with a relatively higher κ. 
Here, an obvious choice is the carbon–based materials, such as the carbon nanotube (CNT) and polymer composites <cit.>. 
In such a composite, a significant increase in κ requires concentration of CNT ϕ_ CNT
exceeding their typical percolation threshold. While a CNT–polymer composite certainly show a significantly higher 
κ than the bare polymers <cit.>, it also has two major drawbacks: (1) it looses the underlying flexibility (typical of polymers) because of a large ϕ_ CNT and their physical properties 
then get dominated by the CNTs present in the background polymers. 
(2) Polymers are rather cost effective, having their typical 
prices of about 2–3 orders of magnitude lower than the CNTs
and thus a CNT–polymer composite inherently becomes significantly costlier. 

A more plausible alternative is the polymer blends,
where the non–bonded interactions can be altered by changing ϕ_ second of the second polymer component.
Here, however, a prerequisite is that the two components remain fairly miscible over the full range of ϕ_ second <cit.>. 
On the contrary, when the two components in a blend phase separate, they create zones within a sample 
consisting of the individual components. These separate zones usually have very weak interfacial interaction 
and thus induce resistance for the heat flow, akin of the Kapitza resistance <cit.>.

In a symmetric polymer blends (i.e., when both polymer components have a comparable degree of polymerization N_ℓ), 
κ varies monotonically between the two pure components <cit.>. 
A recent experimental study, however, reported that an asymmetric blend (consisting of the longer PAA and shorter PAP chains) shows a larger enhancement in κ≃ 1.5 W/Km around a PAP concentration of ϕ_ PAP≃ 30% <cit.>. 
It was argued that the PAP chains act as the H–bonded cross–linkers between the neighboring PAA chains, 
forming a 3–dimensional H–bonded stiff network.

On the contrary, another set of experiments did not attain the same enhancement, instead found that PAA–PAP phase separate around ϕ_ PAP≃ 30% <cit.>.

Motivated by the above contradicting experimental results, a simulation study suggested 
that the miscibility can be enhanced when PAP is replaced with PAM, i.e., a PAA–PAM system consisting of 
a long PAA and a short PAM <cit.>. 
PAA–PAM showed a weak non–monotonic variation in κ with ϕ_ PAM, attaining a maximum κ≃ 0.4 W/Km around ϕ_ PAM≃ 30%. 
Note also that the size of PAM molecules in this simulation study was chosen to be of the order of persistence length ℓ_ p≃ 0.75 nm (or 3 monomers), 
i.e., a PAM as a stiff linker that fits perfectly between two PAA chains. 
When the cross–linker length N_ℓ≫ℓ_ p, they form flexible cross–linking. This on one hand makes a network soft <cit.>, on the other they also induce effective free volume (weak spots) within a network <cit.>. Collectively, these two effects reduce κ. Something that speak in this favor is that the experimental results of cross–linked PAA reported κ≃ 0.27 W/Km, which is about 25% smaller than κ≃ 0.37 W/Km measured in a linear PAA <cit.>. 

The above discussions suggest that there is a need to look beyond the simple amorphous polymers. Therefore, in the 
following section some analytical approaches are first presented that may provide a guiding tool for
the remaining discussions presented herein.



§ ANALYTICAL MODELS
 

In an isotropic material, κ is directly related to the volumetric heat capacity c, 
the group velocity v_ g,i(ν), and the phonon mean–free path Λ(ν) = τ(ν) v_ g,i(ν).
Here, τ(ν) is the phonon life time and ν is the vibrational frequency. 
Starting from the above description, κ can be written as,

    κ(ν) = 1/3∑_i c(ν) v_ g,i^2(ν) τ(ν).
 
For an non–conducting amorphous material a well known theoretical 
description is the minimum thermal conductivity model (MTCM) <cit.> that is discussed in the following.



 §.§ The minimum thermal conductivity model



Following Eq. <ref>, the general expression of κ for a 3–dimensional isotropic system reads <cit.>,

    κ =  (ρ_ N h^2/3k_ B T^2)
        ∑_i ∫τ(ν) v_ g,i^2(ν) 
        ν^2 e^hν/k_ BT/(e^hν/k_ BT -1 )^2
        g(ν)  dν,

where ρ_ N, g(ν), h are the total particle number density, the vibrational density of states, and the Planck constant, respectively. 
Within this description, MTCM uses the Debye model of lattice vibrations in Eq. <ref> and proposes 
that a sample can be divided into regions of size Λ(ν)/2, whose frequencies are given by the low ν sound wave velocities v_i, 
and thus approximates τ = 1/2ν <cit.> and v_ g,i = v_i. This gives,

    κ =  (ρ_ N h^2/6k_ B T^2) 
        (v_ℓ^2 + 2v_t^2 ) ∫ℐ(ν)
        g(ν)  dν,


with

    ℐ(ν) = ν e^hν/k_ BT/(e^hν/k_ BT -1 )^2.


v_ℓ = √(C_ 11/ρ_ m) and v_t = √(C_ 44/ρ_ m) are the longitudinal and the transverse sound wave velocities, respectively. 
Here, C_ 11 = K + 4 C_ 44/3, K is the bulk modulus, C_ 44 is the shear modulus, and ρ_ m is the mass density.

One key quantity in Eq. <ref> is g(ν), which can be calculated by the Fourier transform
of the mass–weighted velocity auto–correlation function c_ vv(t) = ∑_i m_i ⟨v_i(t) ·v_i(0)⟩ <cit.>
obtained from the classical simulations,

    g(ν) = 1/C∫_0^∞cos(2πν t) c_ vv(t)/c_ vv(0) dt.

The prefactor C ensures that ∫ g(ν) dν = 1. The representative g(ν) for four different commodity polymers are shown
in Fig. <ref>.

It can be appreciated that the polymers have many high ν quantum degrees–of–freedom that do not contribute to κ at T = 300 K. 
For example, a C–H bond vibration frequency in a polymer is ν≃ 90 THz, while the representative ν_ room≃ 6.2 THz at T = 300 K. 
Such a mode, together with many other stiff modes (for ν > ν_ room), remain quantum–mechanically frozen at T = 300 K <cit.>. 
If the contributions of these individual modes are not properly incorporated in Eq. <ref> via the Bose–Einstein weighted function in Eq. <ref>, 
one can easily overestimate κ within the classical simulations in comparison to the experimental data <cit.>.   



Standard analytical descriptions typically use the Debye estimate of parabolic vibrational density of states g_ D(ν) = 3 ν^2/ν_ D^3
in Eq. <ref>. Here, ν_ D is the Debye frequency and is written as <cit.>,

    ν_ D = ( 9 ρ_ N/4π)^1/3( 1/v_ℓ^3 + 2/v_t^3)^-1/3.

While g_ D(ν) is certainly a good approximation for the standard amorphous solids when
T ≪Θ_Θ, with Θ_ D = hν_ D/k_ B being the Debye temperature. 
Typical examples are amorphous silica and/or silicon, where Θ_ D≥ 480 K <cit.>.
Moreover, in the case g(ν) is complex, such as in the polymers (see Fig. <ref>), 
simplistic g_ D(ν) may lead to a wrong estimate of the low ν vibrational modes 
and thus leads to the default artifacts in computed κ. This is particularly because Θ_ D≃ 180-220 K 
for the commodity polymers listed in Table <ref>, i.e., 20–40% smaller than T = 300 K <cit.>
where typical experiments and simulations are performed.

When exact g(ν) from Fig. <ref> are used in Eq. <ref>, κ values can be reasonably 
reproduced within 5–20% error, see Fig. <ref>.

Note that this data is obtained by taking v_ℓ and v_t from the experiments <cit.>, 
while ρ_ N and g(ν) are taken from simulations <cit.>.



  §.§.§ High temperature approximation with correction of the stiff modes


Within the high T classical limit, i.e., when all modes are considered in Eq. <ref> <cit.>,
the original MTCM for amorphous polymers can be written as,

    κ_ MTCM   =  (π/48)^1/3 k_ BN^2/3
         (v_ℓ + 2v_t).

The corrections for the stiff modes in a polymer can then be incorporated in Eq. <ref> by considering an effective number of 
atoms N = 2(N - N_ H)/3 and c = 3 N k_ B. Here, N eliminates the stiff modes associated 
with the number of hydrogen atoms N_ H and other stiff backbone modes. Following this, Eq <ref> can be simply written as,

    κ_ MTCM =  (π/432)^1/3 k_ B^1/3c^2/3 (v_ℓ + 2v_t).

Eq. <ref> gives estimates that are 30% larger than the corresponding experimental data <cit.>.



  §.§.§ Accurate computation of κ via specific heat correction








An accurate computation of κ within the standard classical molecular simulation setups is a daunting
task, because polymers have quantum degrees–of–freedom whose exact contribution to the heat balance is rather non–trivial.
Furthermore, T > Θ_ D for the commodity polymers (see Section <ref>)
and thus the polymer thermal properties are dominated by the low ν classical modes that are 
dominated by the non–bonded interactions (or the localized vibrations). On the contrary, the stiff modes (i.e., for ν > ν_ room) remain
quantum–mechanically frozen and do not contribute to κ. In this context, one of the key quantities 
that dictates κ behavior in Eq. <ref> for polymer is its c.

In classical simulations, every mode in a polymer contribute equally to c, i.e., given by the Dulong–Petit classical estimate.
Here, it is well documented that the classically computed c are always overestimated in comparison to the corresponding 
experimental data <cit.> and thus
also leads to an overestimation of κ <cit.>.

Given the above discussion, if c is estimated accurately, it will automatically lead to an accurate computation of κ.
Recently a method has been proposed to compute the quantum corrected c. This method uses the Binder approach to estimate the 
contributions of the stiff harmonic modes <cit.>, which is then used to get the difference between the classical
and the quantum descriptions <cit.>,

    Δ c_rel(T)/k_ B 
        = ∫_0^∞{ 1 - 
        (hν/k_BT)^2 e^hν/k_ BT/(e^hν/k_ BT -1 )^2} g(ν)  dν.

Finally the quantum corrected estimate of c(T) is can be given by,

    c(T) = c_ cl(T) - Δ c_rel(T).

Here, the classical heat capacity is calculated using c_ cl = [H(T- Δ T) - H(T - Δ T)]/2 Δ T 
and H(T) is enthalpy. The main advantage of using Eq. <ref> is that the stiff harmonic modes 
are corrected, while the contributions from the anharmonic (low ν) modes remain unaffected.

Using Eq. <ref>, one can then calculated the quantum corrected κ(T),

    κ(T) = c(T)κ_ cl(T)/c_ cl(T).

Here, κ_ cl(T) is classically computed thermal transport coefficient using the standard equilibrium <cit.> and/or non–equilibrium <cit.> methods.
Fig. <ref> show the computed κ(T) for PMMA using c(T) <cit.>.
It can be seen that the quantum corrected κ(T) compares reasonably with the corresponding experimental data, 
while the classical estimate is about a factor of three too high.
The method proposed in Ref. <cit.> also highlighted different strategies to estimate c(T) 
accounting for the missing degrees–of–freedom (DOF) within the united–atom and/or coarse–grained models. 
A direct implication is that a certain percentage error in κ computed in the united–atom models comes from
the missing DOFs <cit.>.

It is also important to highlight that the simple scaling in Eq. <ref> works reasonably for polymers because 
only the low ν modes dominate their thermal properties. When dealing with the crystalline materials special attention
is need. For example, in a crystal, not only c that has quantum effects, rather v_ g (i.e., stiffness) <cit.> and τ 
also has quantum contributions at low T ≪Θ_ D.



 §.§ Single chain energy transfer



In the introduction, it is discussed that a C–C bond is significantly stiffer <cit.> than the typical non–bonded contacts.
A direct consequence of this microscopic interaction contrast is that the energy transfer between two bonded monomers is over 100 times faster than 
the energy transfer between two non–bonded monomers <cit.>. Taking motivation from such distinct microscopic interactions, 
experimental and computational/theoretical studies have reported a large enhancement in κ in the systems where bonded 
interactions dominate, such as in the single extended chains <cit.>, polymer fibers <cit.>, 
and/or molecular forests <cit.>. However, until recently there existed no direct theoretical 
framework that could quantitatively decouple the effects of these two separate microscopic interactions in dictating the macroscopic heat flow in polymers <cit.>.
Therefore, in this section, the key ingredients of this simple chain energy transfer model (CETM) will be discussed.



Starting from a homogeneous sample consisting of linear polymers, CETM considers the diffusion of energy along a 
chain contour, i.e., between the bonded monomers. This involves multiple hops along a chain
before infrequent energy transfers to the neighboring non–bonded monomer belonging to another chain.
Note that there may also be non–bonded contacts between the two monomers belonging to the same chain, but topologically 
far from one another. However, this will require loop–like conformations of a chain in a dense polymeric system. 
Furthermore, the free energy difference to form such a loop of segment length 𝒩 is given 
by ℱ(𝒩) = mk_ BT ln(𝒩) with a critical exponent m = 1.95 <cit.>.
Within this picture, a loop can only form when it overcomes a free energy barrier of several k_ BT,
which has a very low probability in a dense system. Note also that the CETM method does not distinguish between the
intra– and intre–molecular non–bonded hopping.

Considering the first and the second neighboring bonded monomer transfers along a chain contour,
The rate of change in the internal energy ℰ for any monomer i can then be simply written as <cit.>,

    dℰ_i/ d t = c_ m d T_i/ d t
           =G_b(T_i+1-2T_i+T_i-1)
         +G̃_b(T_i+2-4T_i+1+6T_i-4T_i-1+T_i-2)
         +nG_nb(T_bulk-T_i) ,

with G_ b/c_ m, G̃_ b/c_ m, and G_ nb/c_ m are the bonded, next nearest bonded, and non–bonded energy transfer 
rates, respectively. Here, individually G values are thermal conductances, c_ m is the specific heat of one monomer, n is the number of non–bonded neighbors, T_i is the temperature of the
i^ th monomer, and T_ bulk = 300 K. 
Following the treatment presented in Ref. <cit.>, diagonalizing Eq. <ref> along the chain 
contour will lead to an exponential relaxation of the eigen–modes,

    T̂_p(t) ∝e^-α_p t,
	
with,

    T̂_p (t) = ∑_i = 0^N-1{T_i (t) -T_ bulk}cos[pπ/N(i+ 1/2)],

and

    α_p = 4 G_ b/c_ msin^2(pπ/2N) - 16 G̃_ b/c_ msin^4(pπ/2N)
    	+ n G_ nb/c_ m.

In a nutshell, p gives the effective length scale in a system, i.e., a particular p mode corresponds to a length scale of N_ℓ/p.
Fig. <ref> shows the variation in α_p for PMMA.

Fitting the simulation data (symbols) with Eq. <ref> gives G_ b/G_ nb≃ 63. 
Note also that G_ b/G_ nb≃ 155 for a polyethylene (PE) chain <cit.>. This difference is because of the rather bulky 
side group in PMMA that act as an additional scattering center for the energy transfer. This aspect will be discussed at a later stage. 

It can also be appreciated in Fig. <ref> that α_p almost plateaus for 4sin^2(π p/2N_ℓ) ≥ 0.9 (or p ≥ 11).
This length scale is comparable to ℓ_p≃ N_ℓ/p = 2.7 monomers (or 0.7 nm for PMMA), i.e., a length scale where 
ballistic energy transfer may dominate that is not considered within the formalism of CETM.

The energy transfer rates obtained using Eq. <ref> can also be used to get a theoretical estimate of thermal transport coefficient 
within the Heuristic Random–Walk model <cit.>,

    κ_ HRW =ρ_ N/6[n G_nbr_nb^2+(G_b-4G̃_b) 
        r_b^2 +G̃_br̃_b^2].

Here, r_ nb, r_ b, and r̃_ b are the average distances between a monomer 
and its first bonded, second bonded and non-bonded first shell neighboring monomers, respectively. 
It should, however, be noted that κ_ HRW is underestimated for all investigated commodity
polymers <cit.>.

One of the central assumptions in Eq. <ref> is that the monomers surrounding the reference chain is kept at a constant 
T_ bulk = 300 K <cit.>. This is certainly a good approximation for the common amorphous polymers,
where the heat leakage between the non–bonded monomers is rather weak and mostly restricted upto the
first non–bonded neighbor. Moreover, in the polymers and lubricants 
under high pressure <cit.>, in the confined hydrocarbons <cit.>, 
and/or in the systems where π-π stacking is dominant <cit.>, 
heat leakage between the non–bonded monomers can be significantly enhanced. In these cases, the formalism within CETM may not be directly applicable 
without properly accounting for T_ bulk that will have a gradient as a function of the radial distance from the central chain.










§ THERMAL CONDUCTIVITY OF CROSS–LINKED POLYMER NETWORKS




A predictive tuning of κ purely based on the non–bonded interactions is certainly a non–trivial task, 
if not impossible. Therefore, a more plausible protocol might be to make use of the distinct microscopic interactions (i.e., bonded vs non–bonded),
chain conformation, and possibly also their morphology to understand their effects on κ. 
In this context, one of the most common classes of polymeric materials where the bonded interactions dominate their properties is the epoxies, 
commonly also referred to as the highly cross–linked polymer (HCP) networks.

In a typical HCP, an individual monomer can form more than two bonds (unlike in a linear chain) 
and thus forms a 3–dimensional bonded network. The HCPs are usually light weight high performance materials with extraordinary mechanical 
response <cit.>, attaining E values that can be 2–3 orders of magnitude 
larger than the common amorphous polymers, consisting of linear chains, and
may provide a suitable materials platform toward the enhancement of κ <cit.>. 
Therefore, recent interest has been devoted in investigating the κ behavior in HCPs. 



 §.§ Amine cured epoxies


One common example of amorphous HCP is the amine cured epoxy networks <cit.>, 
where monomers are cross–linked with different amine hardeners with varying stiffness and cross–linker bond length ℓ_ cb, 
see the insets in Fig. <ref>. It can be seen from the experimental data in Fig. <ref> that just by changing the hardener, 
κ can be tuned by about a factor of two <cit.>. Moreover, even in the best case (shown by the red circles in Fig. <ref>), 
κ only increases by a factor of 1.35 in comparison to a linear PMMA, i.e., κ≃ 0.20 W/Km <cit.>.

The specific enhancement in κ is rather small considering that κ in epoxies is expected to be dominated by the bonded interactions. 
What causes such a small variation in κ for HCPs? To answer this question, direct information about the network micro–structure is 
needed where several competing effects control their physical properties. In this context, obtaining any reasonable information regarding such 
microscopic details is a rather difficult task within the commonly employed experimental techniques. 
Therefore, simulations may be of particular interest in studying the κ behavior in epoxy networks, where 
a direct access to the  microscopic network details are reasonable available <cit.>. 

Give the chemical specificity of epoxies, one may expect that the all–atom simulations might be the best possible choice.
Moreover, creating a network structure at the all–atom level is difficult and also is computationally expensive, especially when dealing 
with a broad range of system parameters. Complexities get even more elevated because of the large system sizes coupled with 
spacial and temporal heterogeneity <cit.>. Therefore, an alternative (and possibly a better choice) is a 
bead–spring type generic simulation technique <cit.>. Broadly, generic simulations address the common 
polymer properties that are independent of any specific chemical details and thus a large number of systems can be 
explained within one physical framework <cit.>. Additionally, tuning the system parameters is rather straightforward within a generic setup,
such as the relative bond lengths, their stiffness, and/or bond orientations that are usually inspired by the underlying 
chemical specific systems <cit.>. 





Fig. <ref> shows κ as a function of ℓ_ cb for a set of model HCPs with different bond stiffness <cit.>. 
It can be appreciated that the (relatively) soft bonds (representing the amine hardeners) give reasonably consistent values as 
in the experiments <cit.>, i.e., κ/κ_ amorphous≃ 1.25-1.50 and it is about 1.35 in Ref. <cit.>
(represented by an arrow in Fig. <ref>). 

For the stiff bonds, a significant increase in κ is observed (represented by the diamond symbols in Fig. <ref>). 
This large enhancement can be understood by looking into the networks micro–structures. 
From the simulation snapshots in Fig. <ref>, it can be appreciated that there are large voids (or free volume) 
in all the cured samples. Such voids exist when the neighboring monomers form all their bonds pointing out of each other <cit.>
and the monomers along the periphery of a void only interact via vdW forces. These are usually the weak spots within a network,
hence resist the heat flow. 

The observed void sizes are larger for the tri–functional moderately cross–linked polymers (MCP),
while the tetra–functional HCP have relatively smaller free volume. 
A direct consequence is that the MCPs (open symbols) usually have lower κ than the HCPs (solid symbols) 
in Fig. <ref>. In summary, it is not only that the increasing bonded contacts can by default increase κ. 
Instead the cross–linked bond stiffness, ℓ_ cb, and their effects on the network micro–structures control κ.





 §.§ Ethylene cured epoxies


Another class is the ethylene cured epoxy networks <cit.>, where the free volume can be tuned by changing the length
of the ethylene linkers, see the schematic in Fig. <ref>(a). Here, it is important to mention 
that ℓ_ p≃ 0.65 nm for of a PE chain (or equivalent of one ethylene monomer) <cit.>.
When an ethylene linker is longer than N_ℓ≥ 4, it is soft because of its small flexural stiffness. 
The longer the N_ℓ, more flexible is the linker and thus there are also larger free volume in a sample. 
A direct consequence of such a linker is that– together the free volume and soft linker– they significantly reduce 
stiffness of a materials and as a result κ decreases with increasing N_ℓ, see Fig. <ref>(b).



§ THERMAL CONDUCTIVITY OF CHAIN ORIENTED SYSTEMS




 §.§ Extended chain configurations


A typical representing system where bonded interactions dominate κ is the polymer fibers <cit.>, 
where individual chain are extended along the direction of heat flow and thus κ is dominated 
by the energy transfer between the bonded monomers <cit.>. 
In this context, it has been experimentally reported that a PE fiber can attain κ > 100 W/Km <cit.>,
which is significantly higher than κ≃ 0.3 W/Km for an amorphous PE <cit.>. 



A closer look at an extended chain configuration reveals that it can be viewed as a quasi one–dimensional (Q1–D) 
crystalline material <cit.>. This is a direct consequence of the periodic arrangement 
of monomers along a chain backbone, see the top schematic in the inset of Fig. <ref>.
In such a system, phonons carry a heat current and the coupling strength between the lattice sites is dictated by the bonded interactions. 
Usually a pure (pristine) sample has a large Λ and thus also a high κ≃ 160 W/Km <cit.>, see the main panel in Fig. <ref>. 
However, whenever there appears a kink or a bend along a chain contour (see the bottom schematic in the inset of Fig. <ref>), 
it scatters the phonons. The larger the number of kinks along a chain, the larger the resistance for heat flow and thus a lower κ.
This picture is well supported by the simulation results of an extended PE chain with varying number of kinks, 
see Fig. <ref>.



To investigate the effects of kinks on κ, different PE configurations were specifically engineered <cit.>. 
However, a natural system where the number of kinks and the backbone stiffness can be controlled almost at will is the bottle–brush polymers (BBP) <cit.>. 
A polymer is referred to as a BBP when a linear polymer of length N_ℓ is grafted with the side chains with varying
length N_ s and grafting density ρ_ g. Here, ρ_ g is defined as the 
number of side chains grafted per backbone monomer. For example, if every backbone monomer is grafted with one side chain, then ρ_ g = 1.
BBPs are of interest because of their potential in designing one–dimensional organic nano–crystals <cit.>.

In a BBP, the backbone flexural stiffness (controlling the number of kinks along a chain) is dictated by N_ s and ρ_ g <cit.>.
The heat management in the BBPs is of particular interest because their κ is dictated by two competing effects <cit.>: 
(1) The presence of side chains act as the pathways for heat leakage that effectively reduce κ. 
(2) The side chains increase the backbone flexural stiffness and thus there exists less kinks (or defects) along the backbone, 
which effectively increases κ. To investigate the extent by which these two effects control κ, recent simulations
have been performed using a generic model. The representative data is shown in Fig. <ref>(a).

It can be appreciated in Fig. <ref>(a) that κ shows a non–monotonic variation with ρ_ g, where
two regimes are clearly visible: For ρ_ g≤ 1, scattering because of the side chains reduces κ, while the 
backbone stiffening via side chains increases κ for ρ_ g > 1 <cit.>. This backbone stiffening scenario is also 
supported by g(ν). It can be seen in Fig. <ref>(b) that there are almost indistinguishable changes in g(ν) 
for ρ_ g≤ 1. Moreover, when ρ_ g > 1 a peak becomes more prominent around ν≃ 9-10 t_∘^-1. 
This is associated with the flexural stiffness that increases with increasing ρ_ g, 
as revealed by the shift in this peak towards the higher ν values. The full–width–of–half–maxima ν_ FWHM also 
decreases with increasing ρ_ g and thus increases the phonon life time τ∝ 1/ν_ FWHM (or κ).



 §.§ Molecular forests


The knock down in κ via kinks is a concept that can also be helpful in dictating the heat flow in more complex 
molecular assemblies. One example is molecular forest, where Q1–D are grafted perpendicularly on a
surface forming a two dimensional assembly, such as the forests of CNT <cit.>, silicon nanowires <cit.>, 
and/or polymers <cit.>.
These forests often exhibit intriguing and counter–intuitive physical behavior. 
In this context, it had been experimentally reported that– while a single CNT has κ≥ 10^3 W/Km <cit.>, 
the same CNT in a forest shows a drastic reduction in κ <cit.>. 
This phenomenon is commonly referred to as the heat trap effect (HTE) in the carbon nanotube (CNT) forests <cit.>.

Even when the counter–intuitive HTE phenomenon was known for over a decade, there existed no clear understanding of this behavior.
Here, the simple concepts known from soft matter physics turned out to be reasonably useful in understanding
certain aspects of a hard matter problem of complex molecular assemblies. For this purpose, generic simulations were performed <cit.>.
The key assumption in this model is that a Q1–D is considered as a single extended polymer chain and thus
represents a molecular forest as a polymer brush. The only input parameter in such a modelling approach is ℓ_ p.
In this context, it was readily observed that a CNT can be characterized by their bending stiffness, as measured
in terms of ℓ_ p, that increases with CNT diameter d <cit.>. For example, a CNT of d = 1 nm has ℓ_ p = 50-60 μm. 
Within this picture, a CNT forest of 2 mm height can have about 30-40ℓ_ p
or as many number of kinks. Note that in this simple argument we do not discuss the effect of grafting density Γ. 

The representative data and the corresponding simulation snapshot is shown in Fig. <ref>.
As expected, a chain in a forest shows a significant reduction in κ <cit.>, see the left panel in Fig. <ref>. 
As discussed above, this knock down is direct consequence of the kinks that act as the scattering centers for the heat flow.
The kinks are also evident from the simulation snapshot in the right panel of Fig. <ref>.





§ THERMAL CONDUCTIVITY OF CRYSTALLINE POLYMERS


A somewhat different class to the amorphous polymers is the polymers with certain degree of crystalline order, 
where the long range order facilitates phonon propagation that carry a heat current and thus results in 
an enhanced rate of energy transfer. The typical examples include liquid–crystalline 
materials <cit.>, poly–peptide sequences <cit.>, 
and/or semi–crystalline polymers <cit.>. 



 §.§ Liquid crystalline polymers


Liquid crystals usually have κ≃ 0.3 W/Km <cit.>, i.e., similar to the commodity
polymers. However, one of the advantages of a liquid crystalline material, such as the azobenzene–based liquid crystals, 
is that an azobenzene undergoes a re–entrant trans–to–cis transition when they are exposed to near–ultraviolet light <cit.>. 
Such a transition also alters the molecular order in a liquid crystal and hence κ switches 
between 0.1 W/Km (in cis state) and 0.3 W/Km (in trans state) <cit.>.
This can simply be viewed as a light responsive thermal switch.

When a liquid crystalline polymer is cross–linked with the ethylene linkers, they give very interesting 
and counter–intuitive trends in κ <cit.>. For example, earlier experimental studies have reported 
that κ shows a zig–zag variation with increasing N_ℓ, varying between 1.0–0.2 W/Km <cit.>. 
For an even number of carbon atoms in a linker, κ always has a higher 
value than the next system with a linker with an odd number of carbon atoms. This behavior is commonly known as the odd–even effect in κ,
which was initially reported in a simulation study <cit.>. Moreover, a more detailed investigation is recently reported <cit.>. 
Such an odd–even effect is also well–known in various other properties of the liquid crystalline polymers <cit.>.
While these studies gave very nice insight into the κ behavior of these complex systems,
an exact molecular level understanding of such a non–trivial odd–even effect is still somewhat lacking.



 §.§ Conjugated polymers


Another polymeric system, where crystalline ordering is probably most important, is the conjugated polymer
because they are often used under the high temperature conditions. The crystalline order in such a system
is because of the π-π stacking of their backbone consisting of the aromatic structures <cit.>. 
In this context, a significantly large value of κ→ 2.0 W/Km was reported in poly(3-hexylthiophene) (P3HT) <cit.>.
This κ can also be further increased by blending with the multi–wall CNTs <cit.>. The latter
study reported a non–monotonic variation with ϕ_ CNT, reaching a maximum value of about 5.0 W/Km around ϕ_ CNT≃ 30%.



 §.§ Poly–peptide sequences


A natural soft matter that shows structural order is a poly–peptide sequence. 
Previous experimental results have shown that– by controlling the specific amino acid residues
along a poly–peptide sequence, one can significantly alter its degree of secondary structure d_ s. 
In such a system, an enhancement of up to κ≃ 1.5 W/Km is observed in a hydrated poly–peptide <cit.>. 
Note also that κ≃ 0.6 W/Km for water and hence the observed increase is a direct consequence of d_ s.

Here, however, it is important to mention that the controlled synthesis (aka precision polymerization) of a specific 
poly–peptide sequence is a grand challenge and they are commercially expensive <cit.>. 



 §.§ Semi–crystalline polymers





A more plausible alternative to poly–peptides might be the semi–crystalline synthetic commodity polymers, 
such as the PVA, PLA, and PE systems <cit.>. 

While κ of a semi–crystalline sample can be rather large because of the long range order, 
it may be inferred that if they are cross–linked, the combination of these two effects might lead to an even greater 
increase in κ than the bare amorphous polymers. Motivated by this, a recent experimental study has 
reported that a cross–linked semi–crystalline network can only achieve an increase of up to a factor of 
2.5 times than the pure PMMA sample, i.e., κ≃ 0.5 W/Km <cit.>. This rather surprising behavior was also 
investigated in a simulation study, where a similar increase in κ was observed
for a critical ℓ_ cb, see Fig. <ref>(a). 

A maxima in κ is only observed when ℓ_ cb is comparable to the lattice constant ℓ_ a, such that the 
degree of crystallinity d_ c increases with ℓ_ cb, as revealed by the peak heights in the scattering function in Fig. <ref>(b) <cit.>. 
When ℓ_ cb increases beyond a certain threshold (i.e., ℓ_ cb/ℓ_ a≥ 0.9)–
in one hand d_ c increases only slightly, on the other hand there is a large decrease in the 
bond density ρ_ cb, see Fig. <ref>(c). To summarize the data in Fig. <ref>(a), 
the observed initial increase in κ for ℓ_ cb/ℓ_ a≤ 0.9 is due to the increased d_ c and 
the decrease in κ for ℓ_ cb/ℓ_ a≥ 0.9 is due to the 
reduced ρ_ cb. This readily suggests that a delicate combination of d_ crystal and ρ_ cb controls κ <cit.>.



§ THERMAL CONDUCTIVITY OF POLYELECTROLYTES


In the preceding sections, a short overview of the κ behavior in neutral polymers are presented. 
Possible ideas are also discussed that can be used to tune κ in amorphous systems by macromolecular
engineering. However, there are systems where electrostatic interaction also plays an important role, 
examples include but are not limited to, organic (soft) electronics <cit.>, bio–inspired materials <cit.>, 
and flexible chips <cit.>. 
In these systems it is always desirable to attain a large κ that can act as a heat sink and thus improves device 
performance/durability. Because of this need, extensive efforts have been devoted in studying the κ behavior 
in electrostatically modified polymers <cit.>.

One of the classical examples of polyelectrolytes is the modified PAA with varying degree of ionization. 
Here, a recent experimental study has investigated the effect of pH on ionized PAA, which reported κ→ 1.2 W/Km <cit.>.
This is an enhancement of about 3–4 times than the neutral amorphous PAA, where κ≃ 0.20-0.37 W/Km <cit.>. 
This enhancement was also coupled with a significant increase in the materials stiffness E, i.e., consistent with the 
predictions of the MTCM that κ∝ E <cit.>. The increased E was predominantly because
electrostatic interaction stretches a PAA and thus makes bonded interaction more dominant than in the case 
of an uncharged PAA system. This characteristic extension of an ionized chain is also consistent with the earlier 
studies investigating the effective stretching and ℓ_ p of polyelectrolytes <cit.>.

The κ behavior in the polyelectrolytes suggest that the influence of the electrostatic interaction is rather indirect, i.e., 
they help stretch a chain and thus the bonded interactions become more dominant, which increases κ.
Something may speak in this favor that the electrostatics alone do not influence κ,
as in the case of the ionic liquid consisting of small molecules <cit.>, where κ≃ 0.2 W/Km <cit.>.



§ THERMAL TRANSPORT IN SMART RESPONSIVE POLYMERS


The backbone structure of the commodity polymers are commonly dominated by the C–C covalent bonds, see Fig. <ref>.
Such a bond is extremely strong with its strength of about 80k_ BT and thus these bonds live forever under the unperturbed 
environmental conditions. This creates severe ecological problems, which get even worse when dealing with water 
insoluble polymers, as shown by a few example in the top panel of Fig. <ref>. 
This is one of the main reasons why the recent interest has been diverted to water soluble (H–bonded) polymers, 
referred to as the “smart" polymers <cit.>, as shown by a few examples in the bottom panel of Fig. <ref>. 
Additionally, it is also preferred if a polymer can be bio–degradable and/or pH responsive, such as the acetal–linked 
copolymers <cit.>. 

A polymer is referred to as a “smart" responsive when a small change in the external stimuli can significantly alter 
their structure, function, and stability. These stimuli can be temperature <cit.>, 
pressure <cit.>, pH <cit.>, light <cit.>, and/or cosolvent <cit.>. 
One common example of smart polymer is PNIPAM that shows a coil–to–globule transition in water around T_ℓ≃ 305 K (or 32^∘ C) <cit.>. 
This is a typical lower critical solution (LCST) behavior <cit.>
driven by the solvent entropy <cit.>.

The fast conformational switching of PNIPAM in water may be extremely useful in the thermal applications.
Thermal switching is one such application that controls heat flow in various systems, including, but are not limited to, 
thermoelectric conversion, energy storage, space technology, and sensing <cit.>. 
In this context, the conventional thermal switches often suffer from their slow transition rates and thus also have poor performance.

Recently there has been considerable interest in studying κ in the smart polymer with a goal to attain a fast switching 
in κ <cit.>. 
In particular, experimental studies in the aqueous solutions of PNIPAM <cit.> and PNIPAM–based hydrogels <cit.> 
have shown that their κ behavior follow the same trend as the LCST transition around T_ℓ≃ 305 K 
(or 32^∘ C) <cit.>, see Fig. <ref>.

It is important to note that κ increases with T in the liquids and in the amorphous materials <cit.>
because of an increased vibrations. 
This behavior is also visible in pure water, where a weak increase in κ is observed with T, see the black data set in Fig. <ref>.
Moreover, the sudden drop in κ around T ≥ 303 K (or 30^∘ C) for the PNIPAM concentration above 5 × 10^-3 g/mL
is predominantly due to the coil–to–globule transition of PNIPAM. This drop is likely due to the loss in the number of hydrogen bonds needed 
to stabilize a PNIPAM configuration and the resultant breakage of the water caging around PNIPAM <cit.>. 
These broken water–PNIPAM H-Bonds effectively creates weak interfaces that act as resistance for the heat flow.
Contrary to these results, another experimental study has reported an opposite trend for the concentrated PNIPAM solutions, 
i.e., κ increases above T_ cloud <cit.>. This is simply because a chain under
a high concentration does not collapse into a globule, instead it remains rather expanded surrounded by the other neighboring chains 
and hence the κ behavior become dominated by the energy transfer between the bonded monomers.

Furthermore, these distinct results highlight that the polymer concentration, chain size at a given concentration, 
relative interaction/coordination (monomer–monomer, monomer–solvent, and solvent–solvent), water tetrahedrality around 
a PNIPAM and a delicate balance between these effects play key role in dictating κ behavior in the polymer solutions. 
A detailed understanding of such effects on the κ behavior is a rather open discussion.

Lastly, it might also be important to highlight another (possible) system for the thermal switching application. 
In this context, elastin–like poly–peptides (ELP) are a modern class of biomimetic polymers that also shows LCST transition <cit.>.
One important aspect of the ELPs is their proline isomerization (ProI) that can have either a cis or a trans conformation, which can 
dictate their relative conformations <cit.>. The free energy barrier of such a cis–to–trans transition is 
about 30k_ BT, which the free energy difference between these two states in only about 2k_ B T and thus can be switched via light.
Given the above discussion, ELPs with ProI may also be alternatively used as a thermal switch, similar to that in the liquid crystalline materials <cit.>.



§ CONCLUDING REMARKS


Ever since the seminal publication of Hermann Staudinger <cit.>, the field of polymer science has traversed a long journey 
with many new interesting developments for the future design of advanced functional materials. In the constant quest to find 
new polymeric materials with improved performance, significant attention has been devoted within the field thermal transport 
of polymers over the last 2–3 decades. Especially because the polymeric plastics usually have very low thermal conductivity coefficient κ, which is 
typically a few orders of magnitude smaller than the common crystals. Here, one of the grand challenges is to attain a predictive tuning
of κ via macromolecular engineering. In the context, experiments have investigated a plethora of systems that
include– linear polymers, symmetric and asymmetric polymer blends, polymer composites, cross–linked networks, polymer fibers, 
crystalline polymers and electrostatically modified polymers, to name a few. Motivated by these studies, computational 
studies have also been conducted to establish a structure–property relationship in polymers within the context of their κ behavior.

While it is certainly rather difficult to address all aspects of a huge field of research within one short overview,
in this work an attempt has been made to highlight some of the latest developments in the field of heat conductivity 
in polymers and polymeric materials. In particular, computational results are discussed within the context of the complementary 
experiments with a goal to establish a detailed microscopic understanding that dictates macroscopic polymer properties. 
Available theoretical models are also discussed that may pave the way to guide the future experimental and/or simulation studies.

Some discussions are also presented that showed that the simple concepts know from the basic polymer (soft matter) science <cit.> can be
used to understand a complex problem from an opposite class of hard matter physics <cit.>. This further highlights why polymer science
is such a vibrant and active field of research which is not only restricted within the soft matter community. Rather, it reaches 
across a wide range of interdisciplinary fields. 


Acknowledgement: The contents presented in this review have greatly benefited from the discussions with many colleagues. 
In particular, the development of two key concepts presented here would not have been possible without 
very fruitful collaborations with Martin Müser and Marcus Müller, whom I take this opportunity to gratefully acknowledge.

This draft is a contribution towards a special issue to celebrate the 40th anniversary of Max Planck Institute for Polymer Research. 
I take this opportunity to gratefully acknowledge very fruitful continual collaborations with many MPIP colleagues,
especially Kurt Kremer for numerous stimulating discussions that led to the foundation of my works in MPIP and also after.

I further thank Kyle Monkman for useful comments on this draft.


Conflict of interest: The author declares no conflicting financial interest.


Copyright permission statement: Copyright permissions are obtained for all figures used in this review.

ieeetr


