Investigating Interfacial Thermal Resistance and Effective Thermal Conductivity in Epoxidized Natural Rubber Composites

Muhammad Uzair Riaz ^†,*, Cornelia Breitkopf ^†

1. Institute of Power Engineering, Faculty of Mechanical Engineering, Technical University Dresden, 01069 Dresden, Germany

† These authors contributed equally to this work.

* Correspondence: Muhammad Uzair Riaz

Academic Editor: Frédéric Addiego

Special Issue: Molecular Dynamics Simulations: Bridging Theory and Real-World Applications

Received: January 14, 2026 | Accepted: March 03, 2026 | Published: March 09, 2026

Recent Prog Sci Eng 2026, Volume 2, Issue 1, doi:10.21926/rpse.2601003

Recommended citation: Riaz MU, Breitkopf C. Investigating Interfacial Thermal Resistance and Effective Thermal Conductivity in Epoxidized Natural Rubber Composites. Recent Prog Sci Eng 2026; 2(1): 003; doi:10.21926/rpse.2601003.

© 2026 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.

1. Introduction

Effective thermal transport in polymeric materials is essential for the advancement of wearable electronics, soft robotics, and sustainable composite systems. However, most polymers exhibit inherently low thermal conductivity κ, often in the range of 0.1-0.5 W/m·K [1,2,3], owing to weak van der Waals forces, substantial chain disorder, and strong phonon-phonon scattering in their amorphous regions. This intrinsic limitation is further intensified in heterogeneous or multiphase systems, where the interfaces between dissimilar components create supplementary resistance to thermal conduction. The mismatch in vibrational spectra at such junctions produces a discontinuity in temperature, known as interfacial thermal resistance (ITR) or Kapitza resistance (R_K) [4,5]. In several polymer composites, this interfacial barrier predominates overall thermal resistance, surpassing even the inherent conductivity of the constituent phases.

In recent years, significant research has focused on minimizing ITR through interfacial engineering. Conventional approaches involve embedding high-κ nanofillers, e.g., graphene, carbon nanotubes, or boron nitride nanosheets into polymer matrices to construct continuous phonon pathways [6,7,8]. While proficient in improving bulk conductivity, these systems frequently experience inadequate filler-matrix adhesion and poor phonon coupling resulting from interfacial mismatch. The resulting phonon reflection and scattering at poorly bonded interfaces reduce thermal transport efficiency. As emphasized by [9,10,11,12,13], a novel alternative is to improve intrinsic polymer-polymer interfaces via chemical functionalization, crosslinking, and molecular-level coupling, thereby providing a more sustainable and structurally compatible approach to improving thermal transport.

Within this domain, soft-soft interfaces between polymers, i.e., epoxidized natural rubber (ENR) and poly (butylene succinate) (PBS), represent an effective model for assessing intrinsic interfacial heat transfer [14]. ENR represents a modified (epoxidized) natural rubber that possesses the ability to engage in both hydrogen bonding and covalent reactions due to the presence of polar epoxy functional groups. Conversely, PBS is characterized as a biodegradable aliphatic polyester exhibiting partial crystallinity and polar ester linkages, which enhance robust dipolar interactions with ENR [15]. This synergy consequently facilitates an optimal equilibrium of flexibility and polarity at the interface between these two distinct phases. Furthermore, these materials enable the modulation of phonon transport through interfacial chemistry rather than relying on filler loading. Furthermore, potential transesterification reactions between epoxy and ester groups may form covalent bridges, which can reduce vibrational mismatch and improve phonon transmission. This concept has been recently highlighted in polymer-polymer and polymer-oxide interfaces [16,17,18].

Molecular dynamics (MD) simulations provide a reliable atomistic approach for directly measuring heat transport and interfacial phenomena. MD has been successfully used to investigate crystalline polyethylene, crosslinked epoxy resins, and natural rubber, demonstrating how chain alignment and sulfur-bridge density affect anisotropic thermal conduction [19,20,21]. Two MD frameworks are widely used: Equilibrium molecular dynamics (EMD) and non-equilibrium molecular dynamics (NEMD). The former relies on Green-Kubo formalism, evaluating the thermal conductivity by integrating the autocorrelation of spontaneous heat flux fluctuations under equilibrium conditions. EMD has been successfully employed to study crystalline polyethylene, crosslinked epoxy resins, and natural rubber, depicting the impact of chain alignment and sulfur-bridge density on anisotropic thermal conduction [22,23,24]. Despite being computationally efficient, EMD’s global averaging renders it less useful for systems with explicit interfaces or strong local gradients, where spatially resolved heat transport is crucial [25,26]. Moreover, comparative analyses indicate that EMD results exhibit reduced sensitivity to domain-size effects, as polymer phonon mean free paths are typically shorter compared to those in crystalline solids. In comparison, NEMD explicitly establishes a temperature gradient by applying thermostats to hot and cold regions, driving steady-state heat flow across the sample. The measurable temperature discontinuity across the interface provides a direct route to determine ITR through the relation $R_{k}=\frac{\Delta T}{J}$, where J is the heat flux [27]. This makes NEMD the preferred method for investigating interfacial resistance in polymer bilayers, nanolaminates, and multilayer composites.

Despite extensive NEMD work on polymer-filler systems, e.g., graphene-epoxy, graphene-polyethylene, graphene oxide-polypropylene, and other systems cited in [28,29,30], comparable studies on purely polymer-polymer interfaces remain limited. The molecular complexity of such interfaces, where flexibility, polarity, and reactivity coexist, presents unique challenges for controlling interfacial bonding and interpreting phonon coupling mechanisms. Work by [31,32] indicates that even minor variations in chemical affinity or chain orientation can induce significant changes in ITR, underscoring the need for systematic, chemistry-dependent studies. Recent NEMD investigations on chemically bonded hard-soft interfaces, for example, epoxy-alumina systems functionalized with silane coupling agents, have demonstrated that covalent interfacial coupling can substantially suppress temperature discontinuities and enhance phonon transmission across dissimilar materials [33].

The current research uses NEMD simulations to investigate interfacial thermal resistance and effective thermal conductivity in PBS-ENR-PBS layered systems, with an emphasis on how the epoxidation level and transesterification crosslinking affect interfacial heat transport. Two ENR variants, i.e., ENR-25 and ENR-50 with different epoxy content, serve as model systems for tunable interfacial polarity [34]. Simulations explore how molecular connections affect key parameters, e.g., interfacial adhesion, phonon coupling, and vibrational mode overlaps at the PBS-ENR interface by gradually increasing crosslinking density (0-50%) via covalent transesterification bridges.

This research gives a molecular-level comprehension of how interfacial chemistry affects heat transfer in bio-based elastomeric systems, bridging the gap between filler-dominated composites and intrinsically designed polymer junctions. The insights acquired here apply not only to the design of thermally efficient soft composites but also to the broader development of interactive fibre-rubber composites (I-FRCs), whereby interfacial functioning controls macroscopic thermos-mechanical behavior.

2. Materials and Methods

2.1 Construction of the Simulation Model

To develop a layered PBS-ENR-PBS configuration for interfacial thermal transport analysis, the modelling structures of polybutylene succinate (PBS) and epoxidized natural rubber were created using Moltemplate software (version 2.20.3) [35] and are shown in Figures 1(a), (b). Two different levels of epoxidation were implemented in the starting structures to represent distinct functionalization degrees of the rubber phase. In the ENR-25 model, every third monomer unit along the cis-1,4-polyisoprene backbone was epoxidized, corresponding to ≈25 mol% epoxide substitution. In contrast, every alternate monomer unit was epoxidized in the ENR-50 model, resulting in ≈50 mol% epoxide content.

Figure 1 (a) Polymeric model of polybutylene succinate (b) Polymeric models of ENR-25 and ENR-50.

ENR-25 and ENR-50 were chosen as they are the most researched and commercially accessible grades of epoxidized natural rubber. The former predominantly exhibits a hydrocarbon character suitable for moderate interfacial compatibility, while the latter possesses significantly enhanced polarity and reactive epoxy functionality, thereby augmenting interfacial adhesion and the potential for chemical crosslinking with polar polymers such as polyester [14,36,37].

A PBS/ENR/PBS layered system was constructed in a 3PBC MD simulation box of initial dimensions (100 Å × 100 Å × 300 Å), containing 200 polymer chains in total, as shown in Figure 2. Each PBS chain comprised four monomer units, whereas the ENR-25 and ENR-50 chains were of 10 monomer units each. The chains were randomly packed and arranged through Packmol [38], having 80 chains of PBS on the left and right with 40 chains of ENR in the centre. This configuration generated two identical polymer-polymer interfaces for assessing thermal conduction across chemically diverse layers.

Figure 2 Atomistic model of the PBS-ENR-PBS trilayer system (before optimization) used in MD simulation, comprising 80 PBS chains on each side and 40 ENR chains at the centre.

2.2 Interatomic Potential

The OPLS-UA force field [39] was utilized to model the polymeric systems, simplifying non-polar moieties such as CH, CH₂, and CH₃ into single interaction sites to optimize computational efficiency while preserving the essential physics of long-chain polymer configurations. The system’s total potential energy was defined as

\[ E=E_{bond}+E_{angle}+E_{dihedral}+E_{non-bonded} \tag{1} \]

The first three terms describe the bonded interactions, while the non-bonded interactions were calculated using van-der-Waals interactions.

Non-bonded interactions were simulated using the Lennard-Jones (12-6) potential with a cutoff distance of 10 Å for van der Waals and 11 Å for Coulombic interactions [40]. The Lorentz-Berthelot mixing rules were applied to determine cross-interaction parameters. Bonded parameters for hydrocarbon, with polar carboxylic and alcohol segments, were taken from the OPLS-UA force field. In contrast, other polar functional groups, such as ester and epoxy sites, used OPLS-AA parameters (Table 1) to improve torsional flexibility and charge distribution. The modified dihedral parameters for −C(=O)−O− and −CH₂−CH−O− linkages were adopted from [41] to correct the over-rigidity of the carboxylate moiety. All molecular dynamics simulations were performed using the LAMMPS molecular dynamics package (KOKKOS-accelerated version, August 2022) [42].

Table 1 Force field parameters for MD simulations. The usual factor of $\frac{1}{2}$ in stretching and bending interactions is included in k_a and k_b.

2.3 Compression and Layer Stabilization Protocol

Following initial energy minimization, the preassembled PBS-ENR-PBS structure faced a multi-phase compression-equilibration protocol aimed at attaining the desired density and eliminating voids at the interfacial regions, as shown in Figure 3(a). To preserve the geometric integrity of the trilayer arrangement during densification, harmonic reflective walls were placed on both sides of each PBS-ENR interface [43]. These virtual reflective walls served as a soft boundary that prevents mixing of the two phases while still permitting localized relaxation and energy dissipation.

Figure 3 Polymeric model undergoing different steps to reach the equilibrated state. (a) Energy minimization and compression with reflective walls to prevent mixing of layers of ENR and PBS (b) Five annealing cycles using NVT and Nose and Hoover’s barostat and thermostat. (c) Final equilibration in the NPT and NVT ensembles to obtain a realistic equilibrated system.

The system was initially equilibrated under the canonical (NVT) ensemble at a temperature of 700 K, utilizing a Nose-Hoover thermostat [44,45], with a damping constant of 100 fs to facilitate segmental mobility and enhance molecular reconfiguration. Thereafter, an isotropic compression phase was implemented to elevate the system's density to approximately ≈1.18 g/cm³. A compression was maintained for 200 ps at 700 K, allowing the polymer chains to accommodate the applied volume reduction while preserving the interface geometry.

After the densification phase, the system was cooled gradually from 700 K to 298.15 K over 20 ps, maintaining constant volume and temperature control via the same Nose-Hoover thermostat. During this process, the total potential and kinetic energies, as well as instantaneous density, were constantly monitored to make sure that the energies were converging and the temperature was stable. The final configuration obtained from this stage represented the densified equilibrium state used for subsequent annealing cycles.

2.4 Cyclic Annealing and Equilibration Scheme

A five-cycle thermal annealing strategy was executed to obtain complete structural relaxation and realistic chain entanglement [43]. Each cycle consisted of four steps with a timestep of 0.2 fs: (i) heating from 343.15 K to 363.15 K under NVT conditions for 50 ps, (ii) holding at 363.15 K for 100 ps, (iii) cooling from 363.15 K to 343.15 K for 50 ps, and (iv) relaxing at 298 K and 1 atm under the isothermal-isobaric (NPT) ensemble for 300 ps, as described in Figure 3(b).

Temperature and pressure were regulated using Nose-Hoover thermostat and barostat algorithms, with damping constants of 0.1 ps and 1.0 ps, respectively. This cyclic annealing made it easier to sample different conformations by periodically breaking and reconstructing weak intermolecular connections. This lets the chains get over local energy barriers and pack more realistically.

After five annealing cycles, a long-term equilibration was carried out for 1400 ps at 298.15 K and 1 atm under NPT conditions to ensure volumetric and pressure stability, followed by a 600 ps NVT run at the same temperature to equilibrate the system under constant-volume conditions, as depicted in Figures 3(c), (d). The final equilibrated structure (density ≈ 1.20 g·cm^-3) was then used for five independent NEMD simulations with different random seeds to evaluate thermal transport across PBS-ENR interfaces.

2.5 Non-Equilibrium Molecular Dynamics (NEMD)

Thermal transport across the PBS-ENR-PBS system was evaluated through direct non-equilibrium molecular dynamics (NEMD) simulations employing the Langevin thermostat method [46], which induces a steady-state temperature gradient along the longitudinal (z) direction. Following annealing and equilibration, the system was divided into 20 equidistant slabs, as depicted in Figure 4, with the two terminal regions of 4.0 Å regarded as fixed. Adjacent to them, PBS phases are coupled to independent Langevin heat baths maintained at 320 K (hot) and 280 K (cold), respectively. The intermediate region, including the ENR interlayer, evolved under microcanonical (NVE) conditions, allowing unrestricted phonon transport between the thermostatted boundaries. After an initial equilibration period of 3 ns, the system reached steady state, characterized by a constant net heat flux and temporally stable temperature profile [47].

Figure 4 NEMD setup for the PBS-ENR-PBS system with 4.0 Å fixed boundary layers at both ends and 5.0 Å hot and cold slabs used to establish a steady heat flux across the ENR layer for evaluating interfacial thermal resistance.

The steady-state heat flux (J) was derived from the cumulative energy exchanged between the thermostats and the system. The rates of energy addition and removal, $\frac{dE_{hot}}{dt}$ and $\frac{dE_{cold}}{dt}$, were extracted from the energy logs over the final half of the production trajectory. Under steady-state conditions, these rates converge in magnitude and opposite sign, allowing the net transport rate to be determined according to eq. (2).

\[ \frac{dT}{dE}=\left(\frac{1}{2}\right)\left(\frac{dT}{dE_{hot}}-\frac{dT}{dE_{cold}}\right) \tag{2} \]

The steady-state heat flux was then computed by normalizing this rate with the cross-sectional area of the simulation box, A = L_xL_y, as seen in eq. (3).

\[ J=A_1\times\frac{dT}{dE} \tag{3} \]

expressed in W·m^-2 after conversion of all quantities to SI units.

The spatial temperature distribution, T(z), was ascertained by partitioning the cell into uniform bins along the z-axis and subsequently calculating the mean instantaneous kinetic temperature of atoms within each bin across the steady-state trajectory. To eliminate boundary artifacts, the immobile atoms near frozen walls were excluded from the analysis, retaining only the mobile domain between the hot and cold reservoirs. The resulting temperature profile exhibited two linear regions corresponding to the PBS layers and one central region corresponding to the ENR phase. Linear regression of these regions yielded local slopes m_i and intercepts c_i for each phase, fitted by eq. (4).

\[ T_{i(z)}=m_iz+c_i \tag{4} \]

where i represents PBS₁, ENR, or PBS₂. The temperature discontinuities at the ENR-PBS interfaces were determined from the difference between the extrapolated bulk temperature lines at the respective interface positions z_int,1 and z_int,2, as mentioned in eq. (5).

\[ \Delta T_1=\left|T_{PBS1(z_{\{int,1\}})}-T_{ENR(z_{\{int,1\}})}\right|,\quad\Delta T_2=\left|T_{ENR(z_{\{int,2\}})}-T_{PBS2(z_{\{int,2\}})}\right| \tag{5} \]

The effective thermal conductivity (k_eff) of the entire multilayer structure was subsequently obtained from Fourier’s law (rearranged eq. 6)

\[ R_{k,1}=\frac{\Delta T_1}{J},R_{k,2}=\frac{\Delta T_2}{J},R_{k,total}=R_{k,1}+R_{k,2} \tag{6} \]

The interfacial thermal resistance (R_K), or Kapitza resistance, quantifying the degree of phonon impedance at the junction [48], was computed from the ratio of the temperature jump to the steady-state flux according to eq. 7

\[ k_{eff}=-\frac{J}{\Delta T} \tag{7} \]

where ΔT = dT/dz represents the linear temperature gradient in the ENR layer. As a result, we obtain the k_eff that shows both the bulk and interfacial contributions to the overall heat conduction.

A slab-averaging procedure was employed to compute the local mass/density distribution, facilitating the identification of interfacial positions with atomic precision [49]. The positions of maximum gradients in the normalized density, |dρ/dz|, indicated the boundaries separating PBS and ENR domains. This made it possible to visualize both thermal and structural transitions across the system. To confirm that the interfaces were spatially distinct, void-free, and dynamically stable under steady-state conditions, density-temperature overlay plots were used.

The verification of steady-state convergence and accuracy in NEMD simulations was achieved by monitoring the temporal evolution of the cumulative energy exchange between the hot and cold regions. Similarly, the linearity and symmetric slopes of these curves validated uniform heat transfer and thermal equilibrium throughout the simulation region. This validation ensures that the computed R_K and k_eff values accurately represent true steady-state transport properties.

2.6 Algorithmic Implementation of Transesterification-Based Crosslinking

A reactive transesterification algorithm was utilized to model the interfacial chemical interactions (i.e., covalent bond formation and cleavage) between ENR and PBS in classical MD, using the heuristic framework established by [50]. This automated reaction-mapping method efficiently combines non-reactive and fully reactive potentials.

The algorithm specifies local reactive sites based on defined geometric and energetic requirements. Each reactive event includes two separate groups of atoms:

• The atoms that participate in bonding (illustrated by red squares in Figure 5) engage directly in the exchange of bonds.
• Edge atoms (shown by blue triangles) that maintain local valence and enable the propagation of the reaction.

Figure 5 Algorithm for reaction mapping of the transesterification reaction between poly (butylene succinate) and epoxidized natural rubber. Red squares denote bonding atoms, whereas blue triangles signify edge atoms. The atom indices and bond connections are updated through the reaction map file post-local reaction, thereby facilitating the formation of a cohesive polymer network.

The algorithm also constantly evaluates probable reactive pairings using the reaction map file, which contains data on the topological transformations between pre-reaction and post-reaction states. A cutoff distance of 10.0 Å is set for a bond-exchange event (i.e., C-O bond formation) to achieve optimal yield, converting the ester linkage of PBS and the epoxy moiety of ENR into a trans-esterified junction (as seen in Figure 5). This reaction proceeds via local energy minimization to maintain thermodynamic continuity between bonded and non-bonded states.

The reaction scheme was adjusted to obtain a crosslinking degree of 10-50% in the interfacial region by controlling reactive sites and restricting iterations. This ensured that the ratio of formed interfacial bonds to total epoxy sites matched the desired conversion rate, optimizing material performance.

\[ X_{cross}=-\frac{N_{formed}}{N_{available}}\times100\% \tag{8} \]

To maintain structural stability during bond formation, stabilization was initiated with a damping coefficient of 0.03 to gradually dissipate local energy. In addition, adjacent atoms were temporarily constrained for 200 integration steps to facilitate localized relaxation after the post-reaction event. The effective crosslinking density is governed by limiting the number of reaction events, rather than by the reaction probability. Reactive pairs were identified within a 10 Å cutoff distance, and all geometrically feasible interfacial encounters were assigned a reaction probability of 1.0, while the target epoxy conversion (10-50%) was achieved by capping the number of reactions. Independent random seeds were used to introduce stochastic variability in the selection and spatial distribution of reaction events, thereby enhancing the statistical representativeness of the interfacial crosslinking configurations.

After each iteration, the system follows a thermal equilibration protocol. First, an NVT ensemble was applied at 500 K for 40 ps with a 0.2 fs timestep to increase molecular mobility and facilitate reactive encounters, effectively corresponding to a heating phase. The system was then cooled from 500 K to 298.15 K over another 40 ps with the same timestep and ensemble to allow chain relaxation and network stabilization around the newly formed covalent bonds.

This reaction-heating-cooling protocol, conceptually aligned with the REACTER methodology [51], provides a physically consistent route for simulating progressive interfacial crosslink formation without explicit reactive force fields. It maintains full compatibility with classical OPLS-based potentials while enabling direct control over interfacial crosslink density, thereby allowing systematic evaluation of how bonding degree influences thermal transport and mechanical integrity across ENR-PBS interfaces.

3. Results and Discussion

3.1 PBS-ENR25-PBS System

The thermal transport properties of the PBS-ENR-PBS trilayer system, with ENR-25 as the central layer, were evaluated to set a reference configuration prior to crosslinking reactions. Figure 6 shows the temperature distribution across three distinct linear regions corresponding to the bulk PBS and ENR domains, confirming a uniform steady-state heat flux under the thermal gradient. The PBS layers are represented in red, the ENR core in blue, and the interfacial transition zones in violet. The hot and cold thermostat regions are depicted in light red and light blue, respectively.

Figure 6 Temperature and normalized density profiles of the PBS-ENR-PBS (ENR-25) trilayer system without cross-linking at steady state. Red-shaded areas represent PBS bulk layers, the blue-shaded zone corresponds to the ENR core, and violet bands indicate interfacial regions. ΔT₁ and ΔT₂ mark temperature discontinuities at the PBS/ENR boundaries, while the grey line shows normalized density variation.

The PBS-ENR interfaces were identified by normalizing the mass-density profile and determining the interface sites at ρ_norm = 0.5. These positions correspond to the maxima of the density gradient (|∂ρ/∂z|), signifying regions of considerable compositional change across the phases. Here, the interfacial temperature discontinuity ΔT is defined as the difference between the extrapolated linear temperature profiles of the adjacent bulk regions evaluated at the detected interface position, rather than from raw slab-to-slab temperature differences The detected interface position was found to be insensitive to moderate variations in density bin size and smoothing parameters, with deviations below 1-2 Å, which do not affect the extracted ΔT or interfacial resistance within statistical uncertainty.

The temperature discontinuities (ΔT₁ and ΔT₂) across these interfaces were then determined by extrapolating the linear temperature fits of adjacent bulk regions to the positions of the detected interfaces. In this way, each ΔT value corresponds exactly to the temperature difference at the midpoint of the density transition, ensuring a direct structural-thermal correlation. The resulting interfacial jumps were ΔT₁ ≈ 1.35 K at the PBS₁/ENR interface and ΔT₂ ≈ 2.65 K at the ENR/PBS₂ interface.

The layer-specific conductivities were 0.603 ± 0.230 W/m·K for PBS₁ region, (0.221 ± 0.014) W/m·K (ENR), and (0.221 ± 0.013) W/m·K (PBS₂). The slightly lower thermal conductivity of PBS₂ compared with PBS₁ can be attributed to asymmetric interfacial phonon scattering near the cold side of the system, as reported for polymer multilayers [52,53]. However, sensitivity tests with varying thermostat widths, damping parameters, and reversed heat-flux directions indicate that the ΔT₁-ΔT₂ asymmetry can also be influenced by methodological factors inherent to NEMD simulations. Therefore, the asymmetry should not be interpreted as purely intrinsic, while the combined temperature discontinuity remains within statistical variation, supporting the stability of the extracted total interfacial resistance. The interfacial resistance was calculated from the temperature discontinuity and imposed heat flux according to eq (9)

\[ R_i=-\frac{\Delta T_i}{J} \tag{9} \]

which yielded R_k,1 = (2.747 ± 0.518) × 10^-9 m²K/W and R_k,2 = (2.428 ± 0.673) × 10^-9 m²K/W, leading to a combined R_k,total = (3.740 ± 0.494) × 10^-9 m²K/W, consistent with previously reported polymer-polymer junctions [54]. Each reported value of R_{k, total} corresponds to the mean of five independent NEMD simulations with different initial velocity seeds, and the associated uncertainty is expressed as the standard error of the mean (SEM).

The total interfacial resistance is the sum of the individual Kapitza resistances R_{k, total} = R_k,1 + R_k,2 due to two PBS-ENR interfaces in the simulation cell. The linear temperature gradient observed within the ENR layer indicates bulk-like diffusive heat transport, allowing the two interfaces to be treated as thermally independent scattering centers, consistent with prior NEMD studies of polymer multilayers.

The cumulative energy-exchange profile (Figure 7) further confirms the establishment of a steady state. The red curve corresponds to the cumulative energy added at the hot boundary and the blue curve to that removed at the cold boundary; both exhibit symmetric linear trends whose slopes yield the net heat flux, |J| = 9.50 × 108 W/m. From Fourier’s law,

\[ k_{eff}=-\frac{J}{dT/dz} \tag{10} \]

the effective composite conductivity was determined as k_eff = (0.235 ± 0.014) W/m·K lies between the reported values of neat PBS (0.18-0.22 W/m·K) [55,56,57] and neat ENR (0.13-0.14 W/m·K) [58,59]. This intermediate value reflects partial phonon coupling across PBS-ENR interfaces and moderate enhancement in heat transport relative to pure ENR. The PBS-rich composition (80 wt.%) provides semi-crystalline domains that favour phonon propagation. At the same time, amorphous ENR regions introduce scattering at interfacial boundaries, resulting in the observed conductivity consistent with steady-state heat-flux behaviour.

Figure 7 Cumulative energy exchange for the hot (red) and cold (blue) thermostat regions, with linear fits (black dashed lines) confirming steady-state heat flux and energy conservation across the PBS-ENR-PBS structure.

With 10% crosslinking, the system exhibits an increased heat flux (1.01 × 10⁹ W/m²) and a reduction in ΔT₂ to 2.14 K, whereas ΔT₁ rises slightly to 5.06 K because at low conversion, the formation of isolated transesterification bonds produces a chemically heterogeneous interphase, which locally increases stiffness and phonon back-scattering rather than improving vibrational as shown in Figure 8. The total interfacial resistance decreases modestly to R_{k, total} = (7.467 ± 1.297) × 10^-9 m²K/W, while k_eff improves to (0.247 ± 0.014) W/m·K, as mentioned in Figures 9(a), (b). This enhancement results from the formation of limited transesterification bridges, which locally stiffen the interface and facilitate better phonon transmission by creating covalent coupling sites between PBS ester oxygens and ENR epoxy carbons [17]. The ENR region’s conductivity also rises from 0.196 to 0.267 Wm^-1K^-1, signifying partial energy delocalization through new covalent channels.

Figure 8 Comparison of the interfacial temperature jumps (ΔT₁ and ΔT₂) in the PBS-ENR-PBS (ENR-25) system as a function of crosslinking density. Each horizontal bar represents the mean temperature discontinuity obtained from five independent NEMD simulations. The separation between ΔT₁ and ΔT₂ at each crosslinking level reflects the asymmetry of heat transfer across the two PBS-ENR interfaces, highlighting the influence of crosslinking on interfacial thermal resistance.

Figure 9 (a) Variation of interfacial thermal resistance (ITR) of PBS-ENR25-PBS with increasing crosslinking density. (b) Effective thermal conductivity k_eff of PBS-ENR25-PBS as a function of crosslinking density. Error bars represent the standard error of the mean obtained from five independent NEMD simulations.

At 20% crosslinking, the interfacial bonding becomes more spatially continuous, promoting improved phonon transmission across the PBS/ENR boundary. As a result, both the structural and thermal profiles stabilize, with ΔT₁ ≈ 2.00 K and ΔT₂ ≈ 1.514 K, producing nearly identical results at 30% crosslinking. The interfacial resistances (ITR) converge to R_k,1 = (3.861 ± 1.270) × 10^-9 m²K/W, R_k,2 = (2.079 ± 0.318) × 10^-9 m²K/W, and R_{k, total} = (4.490 ± 0.820) × 10^-9 m²K/W, corresponding to k_eff ≈ (0.232 ± 0.014) W/m·K for 20% and k_eff ≈ (0.231 ± 0.015) W/m·K for 30% crosslinking (Figures 9(a), (b)). The negligible difference between 20% and 30% systems indicates saturation in interfacial phonon coupling. This is due to the establishment of stable ester-ether linkages at most reactive sites. At a crosslinking density of 30%, additional crosslinking mainly hinders chain mobility without generating new heat-transfer pathways, thereby constraining further conductivity enhancements. This saturation phenomenon aligns with earlier molecular-level analyses of chemically crosslinked rubbers and epoxy-based composites, which reported optimal thermal performance at intermediate crosslink densities due to a balance between vibrational continuity and molecular flexibility.

Across all cases, the PBS domains consistently show higher thermal conductivity than ENR, ranging from 0.38-0.52 W/m·K, reflecting their semicrystalline structure and more ordered chain packing. The ENR region continues to be the principal thermal barrier owing to its amorphous nature and localized vibrational modes. The cold-side PBS₂ layer demonstrates marginally reduced conductivity relative to PBS₁ in most systems, signifying asymmetrical phonon reflection at the exit interface, a recognized phenomenon in NEMD polymer multilayers due to local density and orientation fluctuations.

The ENR-25 series shows that moderate interfacial crosslinking significantly improves thermal coupling, but the benefits plateau beyond 30%, indicating a transition to bulk-limited transport. Although the uncertainty ranges partially overlap, the observed saturation of R_{k, total} beyond 20-30% crosslinking, is reproducible and exceeds the level of statistical noise. The trends in k_eff and R_{k, total} reveal that controlled interfacial chemistry, rather than excessive covalent constraints, facilitates efficient heat transfer in PBS-ENR hybrid systems [60].

3.2 PBS-ENR50-PBS System

The thermal response of the PBS-ENR50-PBS trilayer without crosslinking is illustrated in Figure 10. This illustration confirms diffusive heat transport in the presence of a linear temperature gradient across the trilayer. The temperature discontinuities, ΔT₁ ≈ 3.36 K and ΔT₂ ≈ 0.65 K, derived from the half-density (ρ_norm = 0.5) interface criterion, indicate an asymmetric distribution of interfacial resistance, as illustrated in Figure 10. This design results in a total interfacial thermal resistance of (3.93 ± 0.68) × 10^-9 m²K/W and an effective thermal conductivity of (0.245 ± 0.010) W/m·K, as seen in Figures 11(a), (b). In contrast to ENR-25, ENR-50 exhibits a higher baseline ITR due to higher epoxide content. This intensifies interfacial polarity and non-covalent interactions with PBS ester groups, thereby partially improving vibrational mode overlaps despite covalent bonds.

Figure 10 Temperature and normalized density profiles of the PBS-ENR-PBS (ENR-50) trilayer system without cross-linking at steady state. Red-shaded areas represent PBS bulk layers, the blue-shaded zone corresponds to the ENR core, and violet bands indicate interfacial regions. ΔT₁ and ΔT₂ mark temperature discontinuities at the PBS/ENR boundaries, while the grey line shows normalized density variation.

Figure 11 (a) Variation of interfacial thermal resistance (ITR) of PBS-ENR50-PBS with increasing crosslinking density. (b) Effective thermal conductivity k_eff of PBS-ENR50-PBS as a function of crosslinking density. Error bars represent the standard error of the mean obtained from five independent NEMD simulations.

Crosslinking at 10% yields minimal alterations, with the ITR maintained at (3.91 ± 0.89) × 10^-9 m² K/W and k_eff decreasing to (0.236 ± 0.017) W/m·K (Figures 11(a), (b)). The temperature increases (ΔT₁ = 2.93 K, ΔT₂ ≈ 0.92 K) indicate that the initial transesterification reactions do not significantly enhance interfacial heat transfer in ENR-50 (Figure 12). The diminished segmental mobility and elevated local stiffness of the ENR-50 phase restrict the efficacy of isolated covalent bridges in enhancing phonon transmission, in contrast to ENR-25. At low conversion, crosslink formation enhances local rigidity instead of establishing delocalized vibrational pathways, as observed in highly functionalized epoxy-polymer networks [61].

Figure 12 Comparison of the interfacial temperature jumps (ΔT₁ and ΔT₂) in the PBS-ENR-PBS (ENR-50) system as a function of crosslinking density.

At 20% crosslinking, ΔT₁ decreases to around 0.67 K, whereas ΔT₂ increases to approximately 1.74 K, signifying a change in interfacial resistance at the PBS-ENR interface, as depicted in Figure 12. Notwithstanding redistribution, the overall ITR increases to (4.64 ± 1.04) × 10^-9 m²K/W, whereas k_eff diminishes to (0.232 ± 0.012) W/m·K (Figures 11(a), (b)). As epoxidation levels rise, crosslink density excessively constricts interfacial regions, diminishing the low-frequency vibrational modes essential for heat transfer. Densely crosslinked ester-ether networks reduce phonon mean free paths by enhancing anharmonic scattering, surpassing the effects of chemical bonding. Crosslinked epoxy resins and highly grafted polymer surfaces diminish conductivity comparably at elevated functionalization levels.

The thermal reactivity of the ENR-50 system plateaus beyond 30% crosslinking, like the ENR-25 system, with values of ITR = (4.64 ± 1.04) × 10^-9 m²K/W and k_eff = (0.232 ± 0.012) W/m·K, as shown in Figures 11(a), (b). The temperature-jump graphs (Figure 12) demonstrate equivalent ΔT₁ and ΔT₂ values exceeding 20%, signifying that most reactive sites at the PBS-ENR-50 interface are depleted. Additional crosslinking constrains molecular mobility without creating additional heat-transfer pathways, shifting from interface-dominated to bulk-restricted thermal conduction.

4. Conclusions

The present work offers a comprehensive molecular-level analysis of interfacial thermal transport (ITR) in bio-based PBS-ENR-PBS trilayer systems, emphasizing the significant influence of epoxidation level and interfacial chemistry on heat conduction at soft polymer interfaces. NEMD simulations assessed interfacial thermal resistance and effective thermal conductivity through spatially resolved temperature and density profiles, facilitating a direct link between structural changes and thermal discontinuities. In ENR-25 systems, moderate transesterification crosslinking (≈10-20%) markedly enhances interfacial phonon coupling by creating covalent bridges that promote vibrational continuity while not overly restricting chain mobility. Beyond 30% crosslinking, the transport behaviour reaches saturation, as further covalent bonds predominantly enhance interfacial stiffness rather than establish new heat-transfer channels. Conversely, ENR-50 systems demonstrate elevated baseline interfacial resistance and premature saturation, indicative of the intrinsically greater polarity and stiffness linked to increased epoxide concentration. Excessive crosslinking in ENR-50 swiftly inhibits segmental dynamics, thereby constraining further improvement in thermal transport. The results indicate that flexible, bio-based elastomers can attain adjustable thermal characteristics by controlled interfacial chemistry instead of filler-based methods, maintaining mechanical compliance while enhancing heat transfer. This study addresses an underexplored area in molecular-scale simulations of polymer-polymer interfaces, which have received considerably less systematic attention than polymer-filler systems. The insights obtained herein offer practical design principles for thermally optimized, sustainable elastomeric composites and present a comprehensive simulation framework for next investigations on interactive fiber-rubber composites and soft multifunctional materials.

Acknowledgments

The computational time and support from the Zentrum für Informationsdienste und Hochleistungsrechnen (ZIH) in Dresden are also highly appreciated (project:p_thermo_dat).

Author Contributions

Conceptualization, M.U.R. and C.B.; methodology, M.U.R. and C.B; software, M.U.R.; validation, M.U.R. and C.B.; formal analysis, M.U.R. and C.B.; investigation, M.U.R.; resources, M.U.R..; data curation, M.U.R.; writing—original draft preparation, M.U.R.; writing—review and editing, M.U.R.; visualization, M.U.R.; supervision, C.B.; project administration, C.B.; funding acquisition, C.B. All authors have read and agreed to the published version of the manuscript.

Competing Interests

The authors have declared that no competing interests exist.

AI-Assisted Technologies Statement

During the preparation of this manuscript/study, the author used ChatGPT (OpenAI, version 5) for the purpose of checking English grammar mistakes, language refinement, and support in organizing and interpretating data analysis. All authors have read and agreed to the published version of the manuscript.