# Heat transfer through covalent versus non-covalent bonding: a case study in the crystalline π-conjugated P3HT polymer by approach-to-equilibrium molecular dynamics

Cheick Oumar Diarra, Carlo Massobrio, Évelyne Martin

► To cite this version:

Cheick Oumar Diarra, Carlo Massobrio, Évelyne Martin. Heat transfer through covalent versus non-covalent bonding: a case study in the crystalline π-conjugated P3HT polymer by approach-to-equilibrium molecular dynamics. Physical Chemistry Chemical Physics, 2025, 10.1039/D5CP02084K. hal-05198298

HAL Id: hal-05198298

https://hal.science/hal-05198298v1

Submitted on 4 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Cite this: DOI: 00.0000/xxxxxxxxxx

# Heat transfer through covalent versus non-covalent bonding: a case study in the crystalline π-conjugated P3HT polymer by approach-to-equilibrium molecular dynamics

Cheick Oumar Diarraª, Carlo Massobrioª and Evelyne Martinª

Received Date
Accepted Date

DOI: 00.0000/xxxxxxxxxx

Heat transport in a crystalline polymer, poly(3-hexylthiophene) (P3HT), is studied using molecular dynamics. The potential energy is described by density functional theory (DFT), to which a correction is added to account for van der Waals dispersion forces. Approach-to-equilibrium molecular dynamics (AEMD) method is used to create and analyse thermal transients and target transport along or between polymer chains. As expected, the approach leads to a length-dependent thermal conductivity of the system along the covalent polymer backbone. Less expectedly, we observe heat transfer between two polymers, i.e. along a non-covalent heat path, and quantify this transfer using a thermal conductance.

## 1 Introduction

Thermal transport by conduction is intimately associated with the concept of phonons, which are derived from the force constants between atoms. In cases where bond strengths are not as substantial as in covalent bonds, the mechanisms underlying heat transfer are less well-established. Recent theoretical work has concentrated on the case of two-dimensional materials and the transfer associated with van der Waals forces between graphene and hexagonal boron nitride, for example, as cited in¹. The impact of van der Waals forces has also been considered in the total thermal conductivity calculated for a disordered system of small molecules². In the same vein, heat path was experimentally adjusted to pass through more or fewer hydrogen bonds between molecules, thereby increasing the thermal conductivity of the system³.

In order to provide a valuable addition to the existing literature on this subject, and enhance the understanding of the relationship between the nature of bonds and the intensity of heat transfer, we study in the present work a π-conjugated crystalline system, namely poly(3-hexylthiophene) (P3HT). In this system (see Fig. 1), polymers are stacked by π-π interactions, thereby combining covalent character along the polymer chains and non-covalent character between them. In contrast to the amorphous blends of polymers or molecules, where entities may be more or less distant from each other, crystallinity goes hand in hand with coupling along the entire polymer.

Thermal studies have already been conducted on crystalline polymer systems, in particular to identify a hydrodynamic transport regime⁴. In the present study, we expand upon existing models limited to an ideal crystal at T = 0 K by considering strands of varying lengths, to probe length dependence, and at finite temperature to account for thermal fluctuations, but still maintaining a π-stacked arrangement. This configuration aligns with the aligned polymer systems (P3HT in particular) employed for the development of organic thermoelectric generators⁵, for which a comprehensive understanding of thermal transport is crucial in the optimisation of the overall device performances.

The P3HT crystal is studied at room temperature using first-principles molecular dynamics (FPMD). The description of interatomic forces in this scheme, as determined by the calculation of the electronic structure, is independent on empirical assumptions as done in classical molecular dynamics. The density functional theory (DFT), augmented here with the consideration of van der Waals forces, is selected to facilitate the incorporation of interactions that are weaker than covalent or ionic forces, as is the case in the direction of the π-stacking. We resort to FPMD to create thermal transients and exploit them in the specific framework of approach-to-equilibrium molecular dynamics (AEMD).

The paper is organised as follows: in Section 2, crystalline P3HT systems and simulation method are described. In Section 3, thermal conductivities characterising heat transport along polymer chains are calculated and analysed. Section 4 focuses on heat transport between π-conjugated polymers and the calculation of the corresponding thermal conductance. Conclusions are drawn in the final section.

<footnote>ª Université de Strasbourg, CNRS, Laboratoire ICube, UMR 7357, F-67037 Strasbourg, France. E-mail: evelyne.martin@cnrs.fr</footnote>
## 2 Methods and system models

Molecular dynamics (MD) simulations are performed using the Car-Parrinello (CP) method⁶ as implemented in developers version 4.3 of the CPMD code.⁷ For the exchange-correlation component of the total Kohn-Sham energy, the formulation proposed by Perdew, Burke and Ernzerhof (PBE)⁸ is employed. The valence-core interactions are described by norm-conserving Troullier-Martins⁹ pseudopotentials. The valence states are developed on a plane wave basis with a 70 Ry cutoff, and the Brillouin zone sampling is limited to the Γ point. The fictitious electronic mass is taken equal to 400 a.u. and the time step is set to $\Delta t = 0.1$ fs (4 a.u.), to ensure that constants of motion are conserved. The temperature of the ions is controlled by a Nosé-Hoover thermostat¹⁰,¹¹, while the fictitious electronic kinetic energy is regulated by a Blöchl-Parrinello thermostat¹² to ensure that adiabaticity is preserved. The masses of thermostats (in units of corresponding frequencies) are equal to $1100 \text{ cm}^{-1}$ for the ions and $5000 \text{ cm}^{-1}$ for the fictitious electronic degrees of freedom.

The study of thermal properties is conducted through the implementation of temperature transients within the system, employing the AEMD method¹³. This approach involves the application of local thermostats to establish hot and cold blocks, constituting the “Phase 1” of the process. Subsequently, the thermostats are released, and the temperature of the blocks is monitored throughout the approach-to-equilibrium, during a “Phase 2”. The characteristic time of the transient is then utilised to extract the thermal conductivity or the interface conductance depending on the configuration, as detailed in sections 3 and 4.

In order to estimate size effects, supercells of P3HT crystals of varying dimensions were constructed and the configurations are captured at room temperature in figure 1. Periodic boundary conditions are applied in the three directions. The size and number of atoms are Specified in table 1. The size along the *x* direction is systematically adjusted so that the leftmost atom of the chain forms a covalent bond with the rightmost one. This means that an infinitely long polymer chain is obtained. Reference system (A) contains two strands of eight thiophenes each, with the hexyl chains included to allow the structure to naturally adopt the experimentally observed 'Form I' configuration. In this configuration, the backbone formed by the thiophenes is inclined with respect to the stacking direction (z). This system is obtained in accordance with the protocol that was employed in our previous work¹⁴. To summarise, this entails commencing from a planar configuration (not tilted), and subsequently introducing temperature through gradual heating to ambient temperature. Systems B is duplicated along the stacking directions, and contains four strands by unit cell. Systems C and D are analogous to system A, but with larger sizes along the polymer axis *x*. Finally, in system E, hexyl chains have been substituted by methyl groups.

For each of the aforementioned systems, the electronic structure is calculated using a convergence criterion of $10^{-6}$ a. u. The geometry is then relaxed in the corresponding potential energy, and the system is heated to $T = 10$ K, followed by $T = 100$ K, prior to being stabilised for approximately 10 ps at $T = 300$ K.

<table><caption>Table 1 Characteristics of the five models (number of atoms *N*, and size of the orthorombic supercell along *x* ($L_x$), *y* ($L_y$) and *z* ($L_z$))</caption><thead><tr><th>Label</th><th><i>N</i></th><th><i>L</i><sub>x</sub> (Å)</th><th><i>L</i><sub>y</sub> (Å)</th><th><i>L</i><sub>z</sub> (Å)</th></tr></thead><tbody><tr><td>A</td><td>400</td><td>31.0</td><td>15.7</td><td>7.6</td></tr><tr><td>B</td><td>800</td><td>31.0</td><td>15.7</td><td>15.2</td></tr><tr><td>C</td><td>400</td><td>46.6</td><td>15.7</td><td>7.6</td></tr><tr><td>D</td><td>800</td><td>62.0</td><td>15.7</td><td>7.6</td></tr><tr><td>E</td><td>160</td><td>31.0</td><td>7.3</td><td>7.1</td></tr></tbody></table>

## 3 Heat transport along polymers and thermal conductivity

### 3.1 Assessment of the thermal conductivity

A regime of thermal transport along the polymer chains is initiated by subjecting the systems to an initial non-equilibrium state, as depicted in the upper section of figure 2. The *x* position of atoms within the skeleton serves as the basis for delineating two distinct blocks, which are then subjected to thermostats at hot and cold temperatures, respectively. The alkyl chain atoms are kept, via the use of a thermostat, at the temperature of the backbone atom to which they are attached. The number of atoms in each block thus formed is identical. Local thermostats are maintained for approximately 10 ps to stabilise the hot and cold temperatures, as illustrated in figure 3. The corresponding temperature profile is depicted in ochre in figure 2. The resulting shape, which is not uniform in each block, is due to the small size of the system. In the second phase, thermostats are removed, leading to a decline in the temperature of the hot block and an increase in the temperature of the cold block (Fig. 3). This process continues until the temperature of the system is uniform. The duration of the approach-to-equilibrium is system-dependent, being specifically influenced by its thermal conductivity and size. This relationship can be substantiated by solving the heat equation under the same conditions as employed in the molecular dynamics simulation, i.e a 1D heat transport and initial periodic hot/cold conditions¹³. The result is a Fourier series, with a predominant term as follows:

$$T(t,x) = T_{\text{av}} + \frac{2\Delta T_0}{\pi} \sin \left( 2\pi \frac{x}{L_x} \right) e^{-t/\tau} \quad (1)$$

where $T_{\text{av}}$ is the average temperature reached at the end of phase 2, and $\Delta T_0$ is the initial temperature difference between the hot and cold block, equal to 200 K in the present case.

As illustrated by figure 2, the temperature profile during phase 2 assumes a sinusoidal form, as indicated by the perfect superposition of the pink curve with the averaged profile over the duration of the phase 2 (shown in blue). By calculating the temperature difference $\Delta T(t)$ between the hot and cold blocks analytically using the temperature profile from equation 1, we obtain:

$$\Delta T(t) = \frac{1}{L_x/2} \int_{x=L_x/2}^{L_x} T(t,x)dx - \frac{1}{L_x/2} \int_{x=0}^{L_x/2} T(t,x)dx = \Delta T(0)e^{-t/\tau} \quad (2)$$

As demonstrated in figure 4, the present simulations yield the same exponential decay of $\Delta T(t)$. It follows from the aforementioned that the relationship issued from the heat equation between the thermal conductivity $\kappa$ and the transient decay time
<figure><img src="image_2.png" /><figcaption>Fig. 1 Unit cells in the xy-, xz- and yz-planes of models: A, reference model containing two P3HT strands, with each strand consisting of eight thiophenes; B, model A duplicated along the stacking direction (z); C, two P3HT strands of twelve thiophenes; D, two P3HT strands of sixteen thiophenes, i. e. model A duplicated along the strand axis (x); E, two strands of eight thiophenes in which the hexyl side chains have been replaced by methyl groups (P3MT). Carbon atoms are represented in brown, sulfur in yellow, and hydrogen in white. The configurations are captured during a trajectory at T = 300 K.</figcaption></figure>

$\tau$ will be employed in order to calculate the thermal conductivity of the system studied by molecular dynamics. The relationship in question is as follows:

$$ \kappa = \frac{L_x^2 C_v \frac{1}{\tau}}{4\pi^2 V} \quad (3) $$

with the heat capacity equal to $C_v = 3\gamma Nk_B$. N denotes the number of atoms in the supercell of volume $V$, and the factor $\gamma=1.01\pm0.01$ is calculated from the dependence of the total energy with the temperature during the trajectories simulated for P3HT (0.86±0.11 for P3MT).

In order to evaluate the statistical error in the value of $\tau$ and therefore in the corresponding $\kappa$, one has to take the average over different implementations of phase 2. This amounts to extending phase 1, selecting a new starting point for phase 2 fully decorrelated from the first and extracting a new value of $\tau$. This procedure can be repeated at will for each system under consideration to obtain distinct values for the thermal conductivity, each one

of them equally contributing to the average. In our case, three uncorrelated phases 2 have been created.

## 3.2 Length dependence of the thermal conductivity

Figure 5 summarises the calculated thermal conductivities in systems A, C, D and E. The results are presented as a function of the length of the polymer chains, which also corresponds to the size of the simulation cell according to $x$ and to the period of the temperature profile. The figure shows the results of successive AEMDs as small dots, with the large dots and error bars representing the mean and standard deviation of the calculations, respectively.

A comparison of the thermal conductivities of systems A and E, which are identical in length but differ in the presence of hexyl (A) or methyl (E) chains, reveals that the mean values are in perfect agreement, with the error bar being higher in the case of P3MT (E system). This last point can be explained by the much smaller number of atoms in the E system (160 compared with 400 for P3HT). The observed agreement can be interpreted as in-
<figure><img src="image_3.png" /><figcaption>Fig. 2 Temperature profile along the polymer axis during phases 1 and 2 and fit by a sinusoidal curve. Each point in the profile is obtained by averaging the kinetic energies of the atoms belonging to the interval centred at x. P3HT, 4 polymers (system B).</figcaption></figure>

dicating a negligible role for hexyl chains in heat transport along polymers, suggesting that covalent interactions between atoms in the polymer backbone are predominantly responsible for this process.

In contrast, for the diverse P3HT systems with lengths ranging from 3 to 6 nm, a substantial variation in thermal conductivity is observed, which increases with the length of the simulation cell. This size dependence of thermal conductivity has been observed in all the other systems studied by AEMD. To comprehend the physical origin of this effect, it is necessary to revert to the origin of thermal conduction, namely heat carriers, or phonons, in crystalline materials. According to the principles of kinetic theory, thermal equilibrium is established by collisions. In this case, these collisions involve phonon-phonon scattering, and possibly scattering of phonons by defects in the material or surfaces and interfaces. Consequently, heat transport in a material is attributable, on the one hand, to the contribution of the ballistic transport mode of phonons (between two collisions) and, on the other hand, to the diffusive transport of phonons (where collisions dominate). In ballistic transport, phonons evolve without dispersion, and this evolution is limited only by the size of the crystal. Consequently, ballistic thermal conductivity increases with crystal size in the direction of heat propagation¹⁵.

In the diffusive mode, heat transport is predominantly con-

strained by collisions. Consequently, diffusive thermal conductivity is a function of the range of phonon mean free paths existing in the system. In a given material, these free paths can span a broad range of values (from 10 nm to 10 µm in silicon, for example¹⁶). Within the AEMD framework, for a simulation cell of dimension $L_x$ in the direction of heat flow, phonons responsible for transporting heat from the hot block to the cold block can be categorised into two distinct groups. The first category comprises phonons with mean free paths shorter than $L_x$, which undergo collisions as if they were in an infinite crystal. Conversely, phonons with mean free paths greater than $L_x$ will not undergo collisions within the cell and their trajectory will be ballistic. As the size of the simulation cell increases, the proportion of diffusive phonons increases until the size exceeds the largest phonon mean free path of the material. At this juncture, the thermal conductivity ceases to be contingent on the dimensions of the simulation cell, a phenomenon that has been observed, for instance, in crystalline $Ge_2Sb_2Te_5$¹⁷. In such a scenario, the thermal conductivity of the material attains a value observed at the macroscopic scale. These phenomena were modelled analytically by Alvarez and Jou¹⁸, and the resulting law was used in our previous works to extrapolate the AEMD calculations to infinite size. However, extrapolation would not make sense in the present case because the sizes studied are too small compared to the saturation value.

The thermal conductivity of P3HT in the direction of the polymer chains has been determined through experimental methods and documented in a variety of works in the existing literature. For instance, as illustrated in figure 5, the horizontal grey band corresponds to the value of $\kappa$ as determined by Degousée et al.¹⁹ in a mechanically aligned P3HT film. The experimental value is attained by our calculations for system D of length equal to $L_x = 6$ nm. In addition, the size of the crystalline domain is equal to 8-9 nm for P3HT aligned by this method²⁰. However, owing to the computational workload of the present calculations, the study of systems of this length was not possible and the observation of the saturation of the thermal conductivity at a constant value was not feasible. In any case, it should be borne in mind that comparison with experimental results should take into account the contribution of the amorphous domains of smaller thermal conductivity, which adds to that of the crystalline domains to form the average macroscopic thermal conductivity of P3HT.

### 3.3 Other dependencies

In addition to the variation in thermal conductivity as a function of system length in the direction of thermal transport previously discussed, other dependencies were investigated using system B. This system is twice the size in the z direction compared with system A. Initially, the hot and cold blocks are delineated in the x direction, as was previously employed to obtain the results presented in Fig. 5. The thermal conductivities of these systems are presented in Fig. 6, along with those of system D, which is twice as large in the x direction of heat transport. It is evident that the thermal conductivity exhibits a substantially lesser variation when the size is increased along the cross section than when the size is increased in the direction of thermal transport (B vs.
<figure><img src="image_4.png" /><figcaption>Fig. 3 Temperatures in the hot (red curve) and cold (turquoise curve) blocks during phase 1 and phase 2. System B, shown in inset with the delimitations of the two blocks.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 Semi-log plot of the temperature difference $\Delta T$ between the hot and cold block during phase 2 and exponential fit to extract the decay time $\tau$, here equal to $5.17 \pm 0.02$ ps, which corresponds to a thermal conductivity equal to $\kappa = 0.21 \pm 0.01$ W K$^{-1}$ m$^{-1}$. System B.</figcaption></figure>

D). A further comparison of the thermal conductivities of systems A and B reveals that the error bar decreases significantly as the cross section increases. This is due to the fact that there are more atoms in each of the hot/cold blocks, and therefore fewer fluctuations in their average temperature. In summary, the effect of the size of the cross-section can be summarised as an impact on the error bar, but not significantly on the thermal conductivity itself.

System B was also employed to ascertain the thermal conductiv-
tivity in the direction perpendicular to the polymer axis. In order
to do this, the top two polymers were thermalised to $T=400$ K,
and the bottom two to $T=200$ K. In this case, heat propagation
during the approach to equilibrium (phase 2) can no longer take
place via paths including covalent bonds. Instead, it must pass
through $\pi-\pi$ interactions between chains. We therefore expect
the thermal conductivity to be lower. This hypothesis is confirmed
by the results presented in Fig. 6, which show a significant de-
crease in thermal conductivity (from 0.28 to 0.028 W K$^{-1}$ m$^{-1}$)
in the system under investigation when compared with the sys-
tem A of reference. It is important to note that size effects may
also be present in this direction. Ideally, the thermal conductiv-
ity values should be compared when the $\kappa(L)$ curves saturate.
However, the scope of the present study does not permit a full
exploration of this issue. Nevertheless, the conclusion remains
unambiguous: the presence of a gap between polymers does not
impede heat transport, although it is less efficient, which provides
the first indication of the role of non-covalent interactions, in this

<figure><img src="image_6.png" /><figcaption>Fig. 5 Thermal conductivity of P3HT versus length (systems A, C and D), without hexyl chains (P3MT, system E), and range of the experimental data.</figcaption></figure>

case $\pi-\pi$ interactions. This topic is more specifically addressed
in the following section.

4 Heat transport between polymers and thermal conductance

4.1 Assessment of the thermal conductance

A regime of thermal transport between the polymer chains is initi-
ated by subjecting the systems to an initial non-equilibrium state,
as depicted in the right section of figure 7. Local thermostats
are maintained for approximately 6 ps to form a hot and a cold
block, as illustrated in figure 8. The corresponding temperature
profile is depicted in ochre in figure 7, and the evolution of blocks
temperature with time in figure 8. In phase 2, thermostats are re-
moved. The temperature profile maintains a similar shape, with
one polymer at a higher temperature than the other (Fig. 7), but
a decline in the temperature of the hottest one and an increase
in the temperature of the coldest one is observed (Fig. 8). Once
again, this is a sign that there is heat transfer across the gap sepa-
rating the two polymers, and therefore via $\pi-\pi$ interactions. The
transfer continues until the temperature of the system is uniform.
<figure><img src="image_7.png" /><figcaption>Fig. 6 Thermal conductivity of P3HT: reference (system A), impact of the section (system B) and length (system D), and of the heat transport direction (system B).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7 Temperature profiles along the z axis during phases 1 and 2 (left) and illustration of the hot and cold blocks, i.e. upper polymer is hot and lowest one is cold (right). System D.</figcaption></figure>

In what follows, we will examine how this transient can be used to characterise this transfer by means of a thermal conductance G.

Firstly, it should be noted that even in the event of transfer between the two polymers, this process is likely to be less efficient than heat diffusion along the polymer. Consequently, the assumption that the temperature remains consistent throughout the polymers is substantiated, and moreover corroborated by our simulations. This approximation is commonly referred to as the lumped capacitance approximation²¹. It has been employed in first-principles molecular dynamics to analyze heat transfer between two silicon blocks connected by a molecular layer²². The difference in the presence work is that there are two channels of exchange between the hot and cold blocks, i.e. with the polymer above and below in the direction of the π stacking. In such conditions, the transient phase 2 is governed by the equilibrium between the internal energy variation of each polymer (the left-hand term in the equation below) and the exchanges with the two neighbouring polymers (the right-hand term), i.e. for the hot polymer:

$$C_v \frac{\partial T_H}{\partial t} = 2G(T_C - T_H) \quad (4)$$

and for the cold polymer:

$$C_v \frac{\partial T_C}{\partial t} = 2G(T_H - T_C) \quad (5)$$

where G is the thermal conductance, in units of W K⁻¹, and T_H and T_C are respectively the temperature of the hot and cold polymers. By subtracting Eqs. 4 and 5, we obtain an equation governing the evolution of the temperature difference ΔT(t) = T_H(t) - T_C(t) between the hot and cold blocks:

$$C_v \frac{\partial \Delta T}{\partial t} = -4G(\Delta T). \quad (6)$$

The temperature difference therefore decreases exponentially with time according to:

$$\Delta T(t) = \Delta T_0 \exp(-t/\tau), \quad (7)$$

the decay time constant τ being related to the conductance G by:

$$G = \frac{3\gamma N k_B}{4\tau}. \quad (8)$$

As demonstrated in figure 9, the present molecular dynamics simulations yield the same exponential decay of ΔT(t). It follows from the aforementioned that the above relationship (Eq. 8) will be employed to calculate the thermal conductance G between the polymers from the transient decay time τ obtained by molecular dynamics.

## 4.2 Impact of lateral chains and length dependence

The thermal conductances are presented on Fig. 10. The results are presented as a function of the length of the polymer chains, which also corresponds to the size of the simulation cell according to x. The figure shows the results of successive AEMDs as small ochre dots, with the large black dots and error bars representing the mean and standard deviation of the calculations, respectively. The magnitude of interchain thermal conductance is found to be lower in P3MT (system E) in comparison to P3HT (system A). However, the distance d defined as the shortest distance between the polymer chains, is equivalent in both systems¹⁴. This observation indicates that the hexyl chains enhance heat transfer between the polymers, likely by introducing steric interactions that strengthen the coupling between the chains.

Furthermore, in the case where the length of the polymers is doubled (system D), the thermal conductance doubles compared with the reference system (A). This increase is simply an effect of scale: doubling the length Lₓ results in a corresponding doubling of the effective cross-section of the interchain interactions, thereby increasing the thermal conductance. Moreover, if we normalise the conductance by the cross-section, which is generally done experimentally, we obtain a value of 138 ± 20 MW K⁻¹ m⁻² for a length Lₓ = 3 nm and 135 ± 3 MW K⁻¹ m⁻² for Lₓ = 6 nm. The other effect of increasing the length is that the uncertainty in the conductance value is reduced, which can be attributed to the reduction in temperature fluctuations in systems containing more atoms.
<figure><img src="image_9.png" /><figcaption>Fig. 8 Temperatures in the hot (red curve) and cold (turquoise curve) polymers during phase 1 and phase 2. System D, shown in inset with the delimitations of the two blocks.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9 Semi-log plot of the temperature difference $\Delta T$ between the hot and cold polymers during phase 2 and exponential fit to extract the decay time $\tau$, here equal to $6.04 \pm 0.03$ ps, which corresponds to a thermal conductance equal to $G = 1.37 \pm 0.01$ nW K$^{-1}$. System D.</figcaption></figure>

According to our calculations, heat transfer is not limited to the backbone of the polymer where atoms are covalently bonded. In this scenario, it can be hypothesised that the transfer is enabled by $\pi-\pi$ interactions, with a potential contribution from alkyl chains. Here we can assume that it is the $\pi-\pi$ interactions that allow the transfer, with possibly a complementary effect provided by the alkyl chains. It is anticipated that the thermal conductances thus obtained will be lower than those reported in the existing literature, which pertain to interfaces between covalently bonded materials. For instance, a thermal conductance of 700 MW K$^{-1}$ m$^{-2}$ was measured at the TiN/Al$_2$O$_3$ and TiN/MgO interfaces by time-domain thermoreflectance (TDTR) measurement$^{23}$. At the Al/Si interface, the conductance varies between 125 and 200 MW K$^{-1}$ m$^{-2}$ depending on the roughness of the interface and the preparation conditions$^{24}$. A value of 300 MW K$^{-1}$ m$^{-2}$ has been measured at the PMMA/Si interface$^{25}$. At the CuPc/C$_{60}$ interface, the thermal conductance is $450 \pm 130$ MW K$^{-1}$ m$^{-2}$.$^{26}$ Furthermore, Losego et al.$^{27}$ have shown that the thermal conductance increases as the chemical bond becomes stronger, as in the case of the Au/quartz interface. These values remain lower than the thermal conductance at the metal-metal interface, which is, for example, 4 GW K$^{-1}$ m$^{-2}$ at the Al/Cu interface$^{28}$.

## Conclusions

In this article, we have presented the study of heat transfer along and between P3HT polymer chains. The calculations were per-

<figure><img src="image_11.png" /><figcaption>Fig. 10 Thermal conductances between chains of P3MT (system E) and P3HT (systems A and D).</figcaption></figure>

formed using thermal transients and the AEMD method. Regarding thermal conductivity along the chain axis, we observed a size effect, whereby conductivity increased with chain length. The value obtained for the longest system (6 nm) aligns well with experimental data for a mechanically aligned P3HT film. The effects of dimensions perpendicular to thermal transport, studied by stacking four polymers, are negligible. We also demonstrated that thermal conductivity remains unaffected when hexyl chains are replaced by methyls. For thermoelectric applications, the calculated thermal conductivity for the smallest size considered (3 nm) suggests that reducing the size of the crystalline domains could halve the thermal conductivity and increase the figure of merit.

We have also demonstrated that heat transfer does take place even in the absence of covalent bonding. This was first demonstrated by calculating the thermal conductivity in the $\pi-\pi$ stacking direction and confirmed in the second part of the paper by tackling heat transfer between two polymer chains thermalised at different temperatures. We evidenced that the temperature of the hot chain decreased and that of the cold chain increased, indicating heat exchange. We characterised these exchanges by
calculating the thermal conductance, which, as expected from an
exchange not involving strong bonds, has a low value compared
with measurements in the literature. Furthermore, we demon-
strated that hexyl chains significantly increase conductance in
P3HT compared to P3MT. There is no size effect when the poly-
mer chain length is increased apart an increase in exchange sur-
face area.

Author contributions

C. O. Diarra: investigation, software, validation, and visualiza-
tion. C. Massobrio: writing – review and editing. E. Martin: su-
pervision, methodology, and writing – original draft.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the
Supplementary Information.

Acknowledgements

The authors acknowledge the High Performance Computing Cen-
ter of the University of Strasbourg for supporting this work by
providing scientific support and access to computing resources.
Part of the computing resources were funded by the Equipex
Equip@Meso project (Programme Investissements d'Avenir) and
the CPER Alsac calcul/Big Data. Calculations on the larger sys-
tems were performed by using resources from GENCI (Grand
Equipement National de Calcul Intensif) (grant no. 0910296 and
0905071). This work of the Interdisciplinary Institute HiFunMat,
as part of the ITI 2021-2028 program of the University of Stras-
bourg, CNRS and Inserm, was supported by IdExUnistra (ANR-
10-IDEX-0002) and SFRI (STRAT'US project, ANR-20-SFRI-0012)
under the framework of the French Investments for the Future
Program.

References

1 H. M. Dong, H. P. Liang, Z. H. Tao, Y. F. Duan, M. V. Milosevic and K. Chang, *Phys. Chem. Chem. Phys.*, 2024, **26**, 4047–4051.

2 T. Zhou, Z. Li, Y. Cheng, Y. Ni, S. Volz, D. Donadio, S. Xiong, W. Zhang and X. Zhang, *Phys. Chem. Chem. Phys.*, 2020, **22**, 3058–3065.

3 G.-H. Kim, D. Lee, A. Shanker, L. Shao, M. S. Kwon, D. Gidley, J. Kim and K. P. Pipe, *Nature Mater.*, 2015, **14**, 295.

4 Z. Zhang, Y. Ouyang, Y. Guo, T. Nakayama, M. Nomura, S. Volz and J. Chen, *Phys. Rev. B*, 2020, **102**, 195302.

5 N. J. Pataki, S. Guchait, B. Jismy, N. Leclerc, A. Kyndiah, M. Brinkmann and M. Caironi, *Advanced Energy Materials*, 2025, 2404656.

6 R. Car and M. Parrinello, *Physical Review Letters*, 1985, **55**, 2471-2474.

7 CPMD, Copyright IBM Corp. 1990-2022, copyright MPI für Festkörperforschung Stuttgart 1997-2001. Available on GitHub under MIT License. https://github.com/OpenCPMD.

8 J. P. Perdew, K. Burke and M. Ernzerhof, *Phys. Rev. Lett.*, 1996, **77**, 3865-3868.

9 N. Troullier and J. L. Martins, *Physical Review B*, 1991, **43**, 1993-2006.

10 S. Nosé, *The Journal of Chemical Physics*, 1984, **81**, 511–519.

11 W. G. Hoover, *Physical Review A*, 1985, **31**, 1695-1697.

12 P. E. Blöchl and M. Parrinello, *Physical Review B*, 1992, **45**, 9413-9416.

13 E. Lampin, P. L. Palla, P.-A. Francioso and F. Cleri, *Journal of Applied Physics*, 2013, **114**, 033525.

14 C. O. Diarra, M. Boero, E. Steveler, T. Heiser and E. Martin, *Phys. Chem. Chem. Phys.*, 2023, **25**, 15539–15546.

15 R. Saito, M. Mizuno and M. S. Dresselhaus, *Phys. Rev. Appl.*, 2018, **9**, 024017.

16 A. S. Henry and G. Chen, *Journal of Computational and Theoretical Nanoscience*, 2008, **5**, 1–12.

17 I. Bel-Hadj, M. Guerboub, A. Lambrecht, G. Ori, C. Massobrio, E. Martin and ADynMat consortium, *Journal of Physics D: Applied Physics*, 2024, **57**, 235303.

18 F. X. Alvarez and D. Jou, *Applied Physics Letters*, 2007, **90**, 083109.

19 T. Degousée, V. Untilova, V. Vijayakumar, X. Xu, Y. Sun, M. Palma, M. Brinkmann, L. Biniek and O. Fenwick, *J. Mater. Chem. A*, 2021, **9**, 16065–16075.

20 M. Brinkmann, *Journal of Polymer Science Part B: Polymer Physics*, 2011, **49**, 1218–1233.

21 T. L. Bergman, *Fundamentals of heat and mass transfer*, John Wiley & Sons, 2011.

22 T.-Q. Duong, C. Massobrio, G. Ori, M. Boero and E. Martin, *The Journal of Chemical Physics*, 2020, **153**, 074704.

23 R. M. Costescu, M. A. Wall and D. G. Cahill, *Phys. Rev. B*, 2003, **67**, 054302.

24 P. E. Hopkins, L. M. Phinney, J. R. Serrano and T. E. Beechem, *Phys. Rev. B*, 2010, **82**, 085307.

25 M. D. Losego, L. Moh, K. A. Arpin, D. G. Cahill and P. V. Braun, *Applied Physics Letters*, 2010, **97**, 011908.

26 Y. Jin, C. Shao, J. Kieffer, K. P. Pipe and M. Shtein, *Journal of Applied Physics*, 2012, **112**, 093503.

27 M. Losego, M. Grady, N. Sottos, D. G. Cahill and P. V. Braun, *Nature Mater*, 2012, **11**, 502–506.

28 B. C. Gundrum, D. G. Cahill and R. S. Averback, *Phys. Rev. B*, 2005, **72**, 245426.