Bound Rubber as a Transferable Structural Descriptor: Connecting MD-Derived Interfacial Scaling to Continuum Reinforcement Models

TL;DR

This paper introduces a structural descriptor from molecular dynamics simulations to better connect microscopic interfacial structures with macroscopic reinforcement in filled elastomers.

Contribution

The novel contribution is a transferable structural descriptor derived from MD simulations that improves cross-scale predictions in elastomer reinforcement.

Findings

MD simulations provide a bound-layer scaling relation for chain length N=50 as a structural probe.

The regime-partitioned bridge model reduces prediction error from 54.1% to 7.3% by incorporating relaxation physics.

Linear-viscoelastic constraints significantly improve nonlinear PTT calibration, reducing die-swell error by 87%.

Abstract

Filled elastomers often exhibit a low-frequency power-law storage modulus (G-prime), yet quantitative links between molecular interfacial structure and macroscopic reinforcement remain unresolved. This gap is addressed using a hierarchical multiscale framework that integrates coarse-grained molecular dynamics (MD) and dynamic mechanical analysis (DMA). Overall, MD contributes transferable structural descriptors rather than direct macro-rheology prediction. MD simulations yield a bound-layer scaling relation for chain length N=50 in coarse-grained simulations serving as a structural probe. For EPDM master curves, the single-phase fractional Maxwell model is statistically preferred (Delta AICc > 147, n = 56), reflecting limited statistical power; larger datasets (e.g., PC/ABS, n = 952) favor the dual-phase formulation. For cross-scale prediction, an MD-derived effective-volume-fraction baseline (MAPE = 54.1%) provides a structural prior; the regime-partitioned bridge model absorbs relaxation physics not resolved at the MD scale, reducing error to 7.3% (blocked-CV MAPE = 9.5%, with a 2.3% fold-to-fold spread). Linear-viscoelastic constraints improve nonlinear PTT calibration, reducing die-swell error by 87%.

Figures

Full text

1. Introduction

Polymer nanocomposites (PNCs) derive their mechanical properties from hierarchical dynamics spanning molecular ( $[eqn]$ – $[eqn]$ s), mesoscale ( $[eqn]$ – $[eqn]$ s), and macroscopic ( $[eqn]$ – $[eqn]$ s) time scales [1,2,3,4]. These time-scale boundaries are approximate and system-dependent; they serve as order-of-magnitude guides rather than universal limits. At the mesoscale, relevant physical processes include Rouse modes, segmental cooperativity, and chain-cluster rearrangements. At the molecular scale, chains near nanofiller surfaces exhibit subdiffusive, caged motion [5,6,7]; the caging duration and subdiffusive exponent depend on the specific polymer–filler interaction strength and simulation model; direct extrapolation of CG-derived $[eqn]$ to macroscopic relaxation requires the Green–Kubo or TTS bridging procedures described in Section 3. at the mesoscale, Rouse-like segmental rearrangements dominate [8]; and at the macroscale, long-term relaxation governs processability [9,10,11]. A gap of roughly $[eqn]$ – $[eqn]$ in frequency (depending on system size, polymer model, and simulation detail) between MD simulations and bulk rheology precludes direct molecular-to-macroscopic validation [12,13,14], motivating the search for constitutive frameworks that are both parsimonious and physically grounded.

The Generalized Maxwell model:

$$[eqn]$$

requires five–ten adjustable parameters (N modes with moduli $[eqn]$ and relaxation times $[eqn]$ ) to span five decades of frequency [15]. This parameter proliferation hinders physical interpretation. Whether a single effective exponent can describe viscoelastic behavior across caging, Rouse, and terminal regimes has not been tested against molecular simulations.

Fractional constitutive models address these limitations by replacing discrete modes with a fractional order exponent $[eqn]$ , capturing spectrum breadth with fewer parameters [12,16,17,18]. The fractional Maxwell fluid [17,19] gives:

$$[eqn]$$

The fractional Maxwell model (detailed in Section 2) replaces discrete relaxation modes with a single power-law exponent $[eqn]$ that captures broad spectral features of filled elastomers. $[eqn]$ recovers nearly elastic solid behavior (vanishing viscous dissipation), while $[eqn]$ recovers the classical Maxwell viscoelastic fluid with a single exponential relaxation. This formulation corresponds to $[eqn]$ [19] and achieves $[eqn]$ with only three parameters (see Section 3). At low frequencies ( $[eqn]$ ), the model predicts $[eqn]$ , departing from the classical terminal scaling ( $[eqn]$ , $[eqn]$ ) [1,5].

A key distinction applies: the constitutive exponent $[eqn]$ describes a frequency-dependent modulus, whereas $[eqn]$ from MD describes spatial diffusion scaling ( $[eqn]$ ). Readers should note that $[eqn]$ describes local subdiffusion at the molecular scale and cannot be equated with the macroscopic relaxation exponent $[eqn]$ without Green–Kubo integration or time-scale bridging. These quantities are physically distinct; direct equivalence would require Green–Kubo integration beyond current simulation windows. Throughout the present analysis, MD serves as a structural probe—extracting interfacial length scales ( $[eqn]$ ) and caging-regime scaling—not as a source of macroscopic rheological predictions.

Coarse-grained MD simulations using the Kremer–Grest bead–spring model [12] have revealed that polymer chains near nanofiller surfaces form bound layers with restricted mobility [5,6,7,20,21,22,23], increasing the effective filler volume fraction and enhancing the modulus beyond classical predictions. A key structural output is the interfacial empirical scaling relation [6]:

$$[eqn]$$

where $[eqn]$ is the bound polymer layer thickness, $[eqn]$ is the filler volume fraction (dimensionless fraction, e.g., $[eqn]$ for 20 vol%), and x depends on filler geometry and polymer–filler interaction strength. For spherical fillers with strong adsorption ( $[eqn]$ ), the theory predicts $[eqn]$ – $[eqn]$ [6]. The scaling $[eqn]$ is derived from idealized CG simulations with smooth spheres; quantitative transferability to non-spherical or chemically distinct fillers requires separate validation for each system. Establishing a direct link between these MD-derived interfacial properties and macroscopic constitutive parameters remains an open problem.

Filled elastomers exhibit low-frequency power-law behavior ( $[eqn]$ ) rather than the classical terminal scaling ( $[eqn]$ ) [5,9,10]. This deviation—the rheological signature of the Payne effect and filler network formation [10,24]—is well-established for carbon black and silica-reinforced rubbers. However, the quantitative decomposition remains unresolved: previous models treat the non-terminal exponent as a single phenomenological parameter without separating contributions from the percolating filler network and the constrained polymer matrix. In the EPDM/carbon black system examined here, the measured low-frequency $[eqn]$ slope is $[eqn]$ (Supplementary Materials Section S3). Local slope analysis confirms this is not terminal behavior (Supplementary Materials Section S3). The non-terminal exponent $[eqn]$ reflects coupled matrix–network contributions whose separation requires constitutive decomposition (Section 3.2); attributing it to a single mechanism is not warranted without model-based analysis. The central question is whether the apparent slope can be decomposed into its constituent parts.

To address this, a “Matrix–Network Dual-Dynamics” framework is introduced as a diagnostic decomposition—a mechanistic hypothesis to test whether the non-terminal slope can be attributed to separable matrix and network contributions. The fractional Maxwell model is restricted to the linear viscoelastic regime ( $[eqn]$ ) and does not capture large-deformation effects, glass transitions, or nonlinear relaxation. Extension to nonlinear regimes requires separate calibration (e.g., PTT; Section 3.5). This framework models the nanocomposite as two parallel, interpenetrating dynamic phases. The matrix phase is governed by fractional Maxwell dynamics ( $[eqn]$ ), representing viscoelastic relaxation of chains in the dynamically heterogeneous environment created by filler surfaces. The network phase follows power-law gel dynamics ( $[eqn]$ , $[eqn]$ ), representing the percolation-like filler skeleton stiffened by bound rubber. Filler–filler interactions, central to the Payne effect at moderate-to-large strains, are implicitly captured by the network phase ( $[eqn]$ ) in the present linear-regime framework but are not explicitly modeled at the molecular level. MD-derived bound rubber thickness ( $[eqn]$ ) dictates the effective volume fraction ( $[eqn]$ ) required for network percolation, while the DMA-derived low-frequency slope reflects the superposition of both contributions [6,10,24].

The remainder of this paper is organized as follows: Section 2 describes the MD simulation setup, DMA protocols, constitutive model formulation, and statistical methodology. Section 3 then reports two linked but distinct inference layers: (i) constitutive-model inference for low-frequency scaling (single-phase vs. dual-dynamics), and (ii) cross-scale bridge inference where MD-derived structural descriptors constrain reinforcement prediction. Results are presented in a bottom-up sequence: molecular-scale interfacial structure from MD (Section 3.1), followed by model-selection analysis, bridge validation, temperature robustness, PTT calibration, and the time-scale bridge between MD and DMA.

This study employs a two-system validation approach. The primary system, EPDM/carbon black (EPDM/CB, EPDM70), is the primary model system for the dual-dynamics framework: the $[eqn]$ behavior, the interfacial empirical scaling relation ( $[eqn]$ ), and the PTT calibration demonstration (Supplementary Materials Section S6) are all derived from this system. The reference system, PC/ABS (40/60 wt%), serves three specific functions. First, it validates that the fractional Maxwell framework recovers near-terminal behavior ( $[eqn]$ ) when filler-network effects are absent. Second, it provides an independent Green–Kubo TTS benchmark against which MD-derived relaxation times can be compared (Section 3). Third, it demonstrates that dual-phase AICc preference emerges with larger datasets ( $[eqn]$ ), supporting a cautious interpretation of the EPDM single-phase result ( $[eqn]$ ). Comprehensive DMA characterization and TTS bridge results (Supplementary Materials Sections S2 and S5) underpin these comparisons. EPDM/CB is selected as the primary system due to its industrial relevance (automotive sealing, hoses) and the availability of comprehensive DMA datasets across multiple temperatures and filler loadings. PC/ABS serves as a reference for near-terminal relaxation behavior in a well-characterized unfilled system, enabling isolation of the filler-network contribution to non-terminal rheology. The framework is demonstrated for EPDM/CB and PC/ABS; extension to other polymer–filler chemistries requires system-specific re-calibration of the scaling relation and bridge model. Temperature primarily affects the kinetic parameter $[eqn]$ via Arrhenius activation, while $[eqn]$ remains structurally determined—this decoupling is validated in Section 3.3 over the range $[eqn]$ to $[eqn]$ + 130 °C. The two systems are complementary, not directly comparable: EPDM/CB probes filled-system dynamics, while PC/ABS validates the constitutive framework in the unfilled limit.

2. Methods and Approaches

All simulations and experiments followed established protocols detailed in the Supplementary Materials; key aspects are summarized below.

2.1. MD Simulations

Coarse-grained MD simulations of polymer nanocomposites were performed using LAMMPS (version 23 Jun 2022) [25] with the Kremer–Grest bead–spring model [12]. The quantitative spherical filler scaling runs used production systems of 200 chains ( $[eqn]$ beads/chain; $[eqn]$ ; see Supplementary Materials Table S1), with carbon black volume fractions $[eqn]$ %, equilibrated under NPT ensemble (300 K, 1 atm) for $[eqn]$ steps, followed by NVE production runs. Equilibration was verified by monitoring (1) density convergence (<0.5% drift over the final $[eqn]$ steps), (2) MSD of bulk beads reaching the diffusive regime (slope > 0.95 on log–log), and (3) end-to-end vector decorrelation. For $[eqn]$ with $[eqn]$ , the total production time of $[eqn]$ corresponds to ∼3000 Rouse times, sufficient for structural relaxation of the ununentangled chains.

Critical assessment of chain length: The simulated chains ( $[eqn]$ , $[eqn]$ ) are in the Rouse-to-crossover regime, not deeply entangled as in experimental EPDM ( $[eqn]$ ) [12,26]. The MD simulations serve solely as structural probes to extract local interfacial quantities ( $[eqn]$ , density profiles), which are primarily controlled by polymer–filler interaction and segment-level packing rather than long-term entanglement rheology [6,7,27]. Operationally, $[eqn]$ is treated as a near-equilibrium excess-density descriptor of the interfacial zone (thermodynamic/structural), not as a terminal-relaxation observable; the literature evidence and sensitivity runs both support only weak molecular-weight dependence once moderate chain lengths are reached [27,28]. Sensitivity simulations at $[eqn]$ , 50, and 100 show robust bound-layer formation but noticeable variation in fitted scaling exponents (Supplementary Materials Section S1). The $[eqn]$ empirical scaling relation used in this study is therefore treated as a calibrated, system-specific structural descriptor rather than a universal chain-length-invariant exponent. The CG model ( $[eqn]$ , $[eqn]$ ) does not capture entanglement effects present in experimental EPDM ( $[eqn]$ g/mol, $[eqn]$ ). $[eqn]$ and $[eqn]$ are therefore near-equilibrium structural descriptors within the Kremer–Grest framework; their absolute values should not be equated with experimentally measured interphase thicknesses. Extension to entangled chains ( $[eqn]$ , $[eqn]$ ) is a priority for future work.

Mean squared displacement (MSD) in the caging regime ( $[eqn]$ – $[eqn]$ s) was analyzed to extract the subdiffusive exponent $[eqn]$ via $[eqn]$ . $[eqn]$ captures the subdiffusive regime ( $[eqn]$ – $[eqn]$ s) dominated by caging. It provides a spatial mobility map near the filler interface but does not directly predict macroscopic relaxation spectra. The time-scale bridge to DMA frequencies requires separate TTS analysis (Section 3.4). Interfacial density profiles were analyzed to extract bound polymer layer thickness ( $[eqn]$ ).

LJ Unit Mapping: Reduced LAMMPS units were mapped to physical units following Kremer and Grest [12] as validated for coarse-grained systems [7,14]. Based on EPDM monomer parameters ( $[eqn]$ g/mol, $[eqn]$ kJ/mol, $[eqn]$ nm), $[eqn]$ ps is adopted. Sensitivity analysis (Supplementary Materials) confirms that $[eqn]$ is dimensionless and scale-invariant: variation within $[eqn]$ –10 ps and $[eqn]$ –1 nm affects absolute time scales but not scaling exponents. Temperature $[eqn]$ corresponds to ∼300 K (above $[eqn]$ °C for EPDM) [7,29]. The adopted $[eqn]$ ps carries an estimated $[eqn]$ uncertainty ( $[eqn]$ – $[eqn]$ ps) arising from the inherent ambiguity in CG-to-atomistic mapping [7,12]. The $[eqn]$ uncertainty in $[eqn]$ ps propagates as a horizontal shift of $[eqn]$ decades on logarithmic frequency axes for $[eqn]$ and $[eqn]$ . For absolute time comparisons (Section 3.4), this uncertainty is incorporated into the reported $[eqn]$ confidence interval. Because $[eqn]$ and $[eqn]$ are dimensionless or length-scale quantities (respectively), they are invariant to $[eqn]$ rescaling. Absolute time scales—including the TTS-shifted relaxation time reported in Section 3—inherit this $[eqn]$ systematic uncertainty, which propagates as a horizontal shift on the frequency bridge (Supplementary Materials Section S5).

2.2. DMA Measurements

Dynamic mechanical analysis was performed on both material systems using a TA Instruments Q800 (TA Instruments, New Castle, DE, USA). For PC/ABS (40/60 wt%): Comprehensive multi-frequency DMA characterization (5381 data points, 5 frequencies, 29.6–170.4 °C) yielded $[eqn]$ °C, an activation energy of $[eqn]$ kJ/mol, and TTS master curves with Prony $[eqn]$ (see Supplementary Materials Sections S2 and S5 for concise details). For EPDM/CB nanocomposites: Uncured EPDM compounds were used. Over the full fitted frequency window, a single-phase fractional Maxwell fit gave an apparent near-terminal index ( $[eqn]$ ), while local low-frequency slope analysis gave $[eqn]$ (Supplementary Materials Section S3). DMA was performed from $[eqn]$ °C to $[eqn]$ °C ( $[eqn]$ °C), with filler loadings $[eqn]$ % and $[eqn]$ independent samples per formulation.

LVE verification: Linear viscoelastic behavior was confirmed via amplitude sweeps at $[eqn]$ (below the Payne-effect onset $[eqn]$ at $[eqn]$ ; $[eqn]$ ). Pre-sweep stability (30 min time sweeps) and post-sweep verification ( $[eqn]$ ) ensured data were not corrupted by thixotropic drift during the ∼75-min frequency sweeps (see Supplementary Materials Sections S2 and S3). The amplitude sweeps confirmed linearity up to $[eqn]$ . While strain-dependent onset of nonlinearity at specific frequencies cannot be entirely excluded, the Payne effect onset for EPDM/CB composites typically occurs at $[eqn]$ , well above our measurement regime ( $[eqn]$ ).

2.3. Constitutive Modeling

DMA data were fitted to both Generalized Maxwell (5-mode, 10 parameters) and fractional Maxwell models (Equation (2); $[eqn]$ parameters: $[eqn]$ , $[eqn]$ , and $[eqn]$ ). The associated storage and loss moduli are given by:

$$[eqn]$$ $$[eqn]$$

For the dual-phase decomposition, $[eqn]$ ( $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ ). Model selection used AICc (Equation (6)), where $[eqn]$ indicates decisive evidence [30]. Robustness was confirmed with $[eqn]$ correction for residual correlation (Supplementary Materials Section S4). The $[eqn]$ correction is a conservative approximation for moderate serial correlation. We verified that the AICc conclusion (single-phase preference for EPDM) is robust to $[eqn]$ ranging from $[eqn]$ to $[eqn]$ , with $[eqn]$ remaining $[eqn]$ in all cases.

$$[eqn]$$

Fitting accuracy was quantified using $[eqn]$ and MAPE with 95% CIs from bootstrap resampling ( $[eqn]$ ).

2.4. Statistical Analysis

All reported uncertainties are 95% confidence intervals (CIs) from replicate-level bootstrap resampling ( $[eqn]$ independent samples per formulation, $[eqn]$ ). For each bootstrap iteration, replicate master curves were resampled (preserving within-curve frequency structure) and the fractional Maxwell model was re-fitted, yielding empirical distributions of $[eqn]$ , $[eqn]$ , and $[eqn]$ .

Model selection used the corrected Akaike Information Criterion (AICc, Equation (6)), where $[eqn]$ indicates decisive evidence for the preferred model [30]. Because DMA frequency sweeps produce correlated residuals within logarithmic decades, AICc was supplemented with log-decade blocked-CV (Supplementary Materials Section S4): data were partitioned into $[eqn]$ blocks of contiguous frequency decades, and models were trained on $[eqn]$ blocks and evaluated on the held-out block. The single-phase fractional Maxwell model achieved lower blocked prediction error than the dual-dynamics model on the EPDM dataset, indicating that the latter’s additional parameters do not improve out-of-sample accuracy in this window. To further assess robustness, AICc was recalculated with an effective sample size $[eqn]$ to account for residual autocorrelation; the single-phase model remained preferred (ΔAICc $[eqn]$ ; Supplementary Materials Section S4).

For the cross-scale bridge, a hierarchical prediction chain was evaluated on a harmonized frequency-resolved reinforcement dataset ( $[eqn]$ rows): only formulations with complete three-temperature filled/reference pairs were retained (EPDM60/70/80 vs. Pure EPDM; $[eqn]$ ). At each temperature, spectra were restricted to a common overlap band and interpolated onto an 80-point log-frequency grid, yielding 3 loadings × 3 temperatures × 80 frequencies. The log-frequency interpolation (60 points per temperature–loading cell) was chosen to regularize the master-curve alignment. We verified that the original (non-interpolated) data yielded consistent AICc rankings and parameter estimates (within reported confidence intervals), confirming that the interpolation did not suppress spectral features relevant to model selection. Critically, these 720 rows corresponded to $[eqn]$ independent $[eqn]$ design cells with repeated measurements along frequency; frequency was therefore treated as a correlated repeated-measures axis, not as 720 independent experiments. Statistical inferences (bridge model MAPE, blocked-CV error) were interpreted at the design-cell level ( $[eqn]$ ), not at the individual-frequency level ( $[eqn]$ ); the blocked-CV procedure respected this structure by holding out contiguous frequency blocks, preventing frequency-to-frequency leakage from inflating apparent predictive accuracy. On this dataset, three candidates were compared: (i) MD-structural baseline, $[eqn]$ , is an empirically calibrated reinforcement model (after Guth and Gold), valid for dilute-to-moderate loadings of rigid spherical inclusions. Direct application to other filler geometries or polymer chemistries requires independent experimental validation: $[eqn]$ ; (ii) global correction; and (iii) partitioned correction for $[eqn]$ and $[eqn]$ . Corrections were fit in log-space, $[eqn]$ , using candidate basis sets in $[eqn]$ °C)/20 °C and $[eqn]$ (linear, quadratic, and cubic terms with temperature interactions) with ridge penalty $[eqn]$ . To avoid frequency leakage, nested blocked-CV with contiguous log-frequency blocks ( $[eqn]$ ) were used: model class and $[eqn]$ were selected by inner blocked-CV on training folds only, and outer held-out blocks were used for reported errors (Supplementary Materials Section S4).

2.5. Sensitivity Analysis

Global sensitivity analysis using Sobol variance decomposition ( $[eqn]$ Monte Carlo samples) revealed the total-order sensitivity ranks $[eqn]$ ( $[eqn]$ ). We verified convergence by comparing results at $[eqn]$ and $[eqn]$ samples: the 95% CI bounds for $[eqn]$ shifted by <0.003 (relative change $[eqn]$ ), confirming that $[eqn]$ samples provided adequate convergence for the reported parameter uncertainties. The focus on $[eqn]$ is justified by its physical significance: it governs the power-law slope of the relaxation spectrum.

3. Results and Discussion

Notation guide: For clarity, Table 1 summarizes all key constitutive parameters determined in this study. A distinction is made between single-phase fits (fractional Maxwell model alone) anddual-dynamics fits (parallel matrix + network decomposition). All uncertainties are 95% confidence intervals unless otherwise noted.

3.1. Interfacial Structure from MD

Coarse-grained MD simulations of EPDM nanocomposites with spherical carbon black fillers revealed the molecular-scale structure responsible for caging dynamics. Figure 1a shows a representative equilibrium configuration of the EPDM/CB nanocomposite (see caption for system details), with polymer chains (blue beads) confined between spherical filler particles (red spheres). The bound rubber layer—visible as the region of elevated polymer density near filler surfaces—emerges spontaneously from the simulation without any ad hoc assumptions. A slice view (Figure 1b) reveals the internal structure, showing that polymer chains were compressed and immobilized in the gaps between fillers.

The density profile analysis (Figure 2) quantifies the extent of polymer confinement. Near the filler surface ( $[eqn]$ ), the polymer density is elevated to $[eqn]$ in the first adsorption layer. Pronounced oscillatory layering characteristic of dense packing against a hard wall is observed. The bound-layer thickness—defined as the distance where the density envelope drops to within 10% of $[eqn]$ —is $[eqn]$ nm for the flat-wall reference system, consistent with the range 2–5 nm reported from solid-state NMR [31]. This planar-reference value should not be directly equated with the spherical filler scaling values (e.g., ∼0.8–0.9 nm at $[eqn]$ ), because geometry, averaging procedure, and overlap handling differ between the two analyses. This confined layer is the structural origin of the subdiffusive caging dynamics observed in the MSD analysis.

These MD simulations revealed two key features of caging dynamics. First, the bound polymer layer thickness ( $[eqn]$ ) scales with filler volume fraction ( $[eqn]$ ) according to a power-law empirical relation:

$$[eqn]$$

where the scaling exponent $[eqn]$ ( $[eqn]$ , $[eqn]$ , $[eqn]$ –30% excluding 20%; 95% bootstrap CI from $[eqn]$ : $[eqn]$ ; the wide CI reflects the limited and scattered dataset) was measured using the primary $[eqn]$ chain system. The $[eqn]$ % rerun yielded anomalous cluster extraction (zero detected clusters, $[eqn]$ vs. ∼1.95 at neighboring loadings) and was excluded from the fit (see Supplementary Materials Section S1). Sensitivity to exclusion: Re-fitting with the 20% point forced into the dataset (using $[eqn]$ nm from the visual density envelope) yielded $[eqn]$ and $[eqn]$ , well within the 95% CI of the primary fit; the qualitative conclusion (superlinear scaling) was robust to this exclusion. This logic underscores the interpolative, not predictive, nature of the relation. In Equation (7), $[eqn]$ was evaluated as a fraction ( $[eqn]$ – $[eqn]$ ), not as percent integers. The superlinear exponent ( $[eqn]$ ) reflected increasing bound-layer overlap at higher $[eqn]$ , consistent with a transition from isolated to interacting adsorption layers in this calibrated dataset. Because only five non-overlapped loading points entered the fit with moderate $[eqn]$ , this law was used as an interpolation-level structural descriptor within the studied loading window, not as a universal exponent. Chain-length sensitivity analysis (Supplementary Materials Section S1) indicated the qualitative robustness of bound-layer formation while leaving uncertainty in the absolute exponent value. Second, mean squared displacement (MSD) analysis in the caging regime ( $[eqn]$ – $[eqn]$ s) revealed subdiffusive scaling:

$$[eqn]$$

with $[eqn]$ , characteristic of polymer chains confined between multiple spherical fillers. $[eqn]$ reflects strong geometric caging in the CG model ( $[eqn]$ , smooth spheres). For longer chains with entanglements or chemically heterogeneous fillers, the absolute $[eqn]$ value will differ. The qualitative finding—that filler proximity reduces local mobility—is expected to be robust, but quantitative values are model-specific. This subdiffusive exponent is consistent across filler loadings ( $[eqn]$ –30%), indicating that caging dynamics govern short-term motion regardless of filler concentration.

To isolate the effect of wall–polymer interaction strength from multi-filler geometric confinement, a complementary set of flat-wall simulations are performed at varying $[eqn]$ (Figure 3). In the flat-wall geometry—where chains experience single-surface confinement rather than the multi-body caging of the spherical filler system—the MSD exponent is $[eqn]$ – $[eqn]$ for all $[eqn]$ values tested (2–8 $[eqn]$ ). The ∼5× reduction from flat-wall ( $[eqn]$ ) to spherical multi-filler ( $[eqn]$ ) geometry reflects the cumulative confinement from surrounding filler particles, which is qualitatively more restrictive than single-surface contact. A quantitative decomposition of geometric versus energetic contributions would require systematic variation in filler separation at fixed $[eqn]$ , which is beyond the current study.The weak $[eqn]$ -dependence of $[eqn]$ in the flat-wall geometry suggests geometric confinement dominates within the LJ potential framework. In real systems, specific polymer–filler chemistry (hydrogen bonding, polar interactions) may have a larger effect, potentially altering the relative importance of geometric versus energetic contributions.

Reconciliation of interaction parameters across geometries: A clear distinction exists between the flat-wall $[eqn]$ -series and the spherical filler production systems. The flat-wall simulations use elevated $[eqn]$ – $[eqn]$ to systematically probe the effect of adsorption strength on single-surface dynamics; the spherical filler production runs (from which the empirical scaling relation in Equation (7) is derived) use the standard Kremer–Grest polymer–filler LJ interaction ( $[eqn]$ , identical to bulk polymer–polymer interactions), where confinement arises primarily from multi-body geometric caging rather than enhanced surface attraction [7,12]. The interaction settings and primary production systems are summarized in Supplementary Materials Section S1; absolute $[eqn]$ values are therefore not directly comparable between different geometries, but the qualitative finding—that confinement suppresses MSD—is consistent across both.

Table 2 consolidates $[eqn]$ values across all systems and analysis methods, revealing how simulation geometry and spatial position govern the degree of subdiffusion.

To directly visualize how interfacial constraints affect chain dynamics at the molecular level, MSD is analyzed as a function of distance from the filler surface (Figure 4). Polymer segments are classified into three regions: the interfacial layer ( $[eqn]$ ), a transition zone ( $[eqn]$ ), and the bulk region ( $[eqn]$ ). A clear spatial gradient in subdiffusive dynamics is observed: segments in the interfacial layer ( $[eqn]$ ) exhibit strongly suppressed mobility ( $[eqn]$ ), approaching the plateau limit expected for beads effectively immobilized within adsorption potential wells. This near-zero exponent reflects a cage-rattling regime where beads oscillate within local energy minima created by the filler surface without escaping to adjacent sites on the caging time scale; similar values ( $[eqn]$ ) have been reported for strongly adsorbed polymer segments near attractive walls in atomistic simulations of filled polymer systems [21]. The interfacial $[eqn]$ is lower than typical experimental values from NMR or neutron scattering, reflecting the idealized nature of the CG model (smooth, strongly adsorbing surfaces without chemical heterogeneity). The qualitative gradient from bulk to interface is physically meaningful, but absolute interfacial $[eqn]$ values should not be compared directly with experimental measurements. Segments in the bulk show substantially higher mobility ( $[eqn]$ ). This gradient ( $[eqn]$ ) provides qualitative molecular evidence consistent with the interfacial-constraint effect that underlies the present multi-scale characterization approach, but a direct quantitative link between $[eqn]$ and $[eqn]$ requires Green–Kubo integration not performed in this study.Three replicas per $[eqn]$ and 12 spatial segments per position provide sufficient statistics for the qualitative mobility trends reported (bulk vs. interface). However, confidence intervals for extreme subdiffusion ( $[eqn]$ ) may underestimate true variability due to the small sample size and sensitivity to local fluctuations. More extensive sampling would be needed for quantitative subdiffusion analysis in the interfacial region.

The empirical scaling relation is summarized in Figure 5.

3.2. Dual-Dynamics Framework and Model Selection

The macroscopic rheological response of EPDM nanocomposites is characterized using DMA master curves. As shown in Figure 6, the storage modulus $[eqn]$ at low frequencies exhibits a power-law slope of $[eqn]$ , which substantially deviates from the terminal flow prediction of $[eqn]$ . Within the dual-dynamics framework, this is an apparent slope arising from the combined matrix and network contributions rather than a single-phase exponent.

Statistical model selection decisively favors the single-phase fractional Maxwell model over the dual-dynamics decomposition on the EPDM70 dataset ( $[eqn]$ ): AICc $[eqn]$ vs. $[eqn]$ ( $[eqn]$ ). The single-phase AICc preference should not be interpreted as evidence that two-phase dynamics are absent in EPDM, but rather that the current dataset ( $[eqn]$ ) cannot statistically resolve them. With the larger PC/ABS dataset ( $[eqn]$ ), dual-phase preference emerges clearly ( $[eqn]$ ), suggesting that increased spectral resolution would likely favor two-phase models for EPDM as well., and blocked-CV confirms lower out-of-sample error for the simpler model (Table 1). This is an explicit Occam’s razor outcome: at the current EPDM sample size, the extra dual-phase parameters are not statistically warranted. The dual-dynamics fitting results are nevertheless presented below for mechanistic interpretation, not as the statistically preferred model.

The parallel model (Equation (9)) was fitted using differential evolution with relaxed parameter bounds ( $[eqn]$ ). For 20% EPDM/CB, the best-fit parameters were: matrix index $[eqn]$ , matrix relaxation time $[eqn]$ s, and network index $[eqn]$ (converging to nearly elastic behavior). The very low $[eqn]$ reflected the near-elastic behavior of the filler network, consistent with a percolating CB skeleton stiffened by bound-rubber bridges. In the limit $[eqn]$ , the network phase became purely elastic; $[eqn]$ indicated extremely slow, nearly frozen relaxation of filler–filler junctions, consistent with the quasi-static nature of CB aggregate rearrangement at small strains (below the Payne effect onset). The model achieved high fit quality ( $[eqn]$ ) but did not improve upon the single-phase model in statistical terms. The dual-dynamics decomposition was therefore a physical hypothesis consistent with MD-derived interfacial structure, retained for qualitative insight into the origin of the non-terminal slope.

Physical interpretation of the matrix exponent shift: The reduction in $[eqn]$ from 0.97 (single-phase fit) to 0.60 (dual-dynamics fit) requires physical interpretation. In a single-phase model, the lone fractional exponent must capture the composite frequency dependence of both matrix and network contributions, resulting in an “averaged” value near unity that reflects the matrix-dominated high-frequency response. When the network contribution is explicitly separated, the remaining matrix phase no longer needs to account for the low-frequency network stiffening.

On the matrix exponent value: The fitted value $[eqn]$ requires careful interpretation. Uncured EPDM is a high-molecular-weight melt ( $[eqn]$ g/mol, $[eqn]$ ; note: this refers to the experimental polymer, not the CG simulation chains with $[eqn]$ ), where physics predicts near-terminal behavior ( $[eqn]$ ) at low frequencies. The apparent value $[eqn]$ is therefore an effective parameter that implicitly absorbs multiple physical effects.

Strain amplification (primary effect): In heterogeneous systems with rigid inclusions, the local strain is amplified: $[eqn]$ , where $[eqn]$ (Guth–Gold formula) [10,32]. At nominal $[eqn]$ , this yields $[eqn]$ ; at effective $[eqn]$ (including bound rubber), $[eqn]$ . $[eqn]$ assumes rigid spheres with ideal dispersion. For industrial CB with complex aggregate structures and heterogeneous spatial distribution, the actual strain amplification may differ from 3.0 by a factor of 2–5. The value serves as an order-of-magnitude estimate to illustrate the physical mechanism, not as a precise prediction. The present model does not explicitly incorporate strain amplification; the reduced $[eqn]$ absorbs this spectral broadening effect. A corrected intrinsic value would be $[eqn]$ – $[eqn]$ , closer to terminal expectations.Geometric confinement: At $[eqn]$ – $[eqn]$ , inter-filler gaps become comparable to chain dimensions, truncating long-term relaxation before the true terminal regime is reached [33,34].Mathematical compensation: Once the network term captures low-frequency stiffening, the matrix term no longer needs to account for that regime, redistributing spectral weight.

No claim is made that $[eqn]$ represents intrinsic chain dynamics. It is an effective parameter absorbing strain amplification ( $[eqn]$ at $[eqn]$ ) and geometric confinement; its proximity to the Rouse value ( $[eqn]$ ) is fortuitous. The value depends on model structure (removing the network term gives $[eqn]$ ), frequency range, and filler morphology; cross-material comparisons require re-fitting. These parameters apply to EPDM/CB at the tested conditions and should not be treated as universal constants.

Model competition and interpretation: The parallel model is defined as:

$$[eqn]$$

On the EPDM70 dataset ( $[eqn]$ ), AICc strongly favors the single-phase model ( $[eqn]$ vs. $[eqn]$ ; ΔAICc $[eqn]$ ), and blocked-CV shows lower out-of-sample error. When the $[eqn]$ lower bound is relaxed from 0.20 to 0.01, the optimizer converges to $[eqn]$ , indicating that this exponent is not independently constrained by the data. For PC/ABS ( $[eqn]$ ), the dual-phase model is preferred ( $[eqn]$ vs. $[eqn]$ ), reflecting greater statistical power with the larger dataset. The dual-dynamics decomposition is therefore a physical hypothesis consistent with MD-derived interfacial structure, not a statistically required conclusion. Model-selection robustness and identifiability diagnostics are summarized in Supplementary Materials Section S4.

3.3. Literature Comparison

The relationship between bound rubber thickness and modulus enhancement in filled elastomers has been the subject of extensive debate. Early studies by Heinrich and Kluppel [10] proposed that filler networking dominates reinforcement, with bound rubber playing a secondary role. However, more recent work by Berriot et al. [6] and Papon et al. [31] demonstrated that the glass transition temperature shifts near filler surfaces, creating an immobilized layer with distinct viscoelastic properties; newer AFM and NMR studies report compatible interphase-constrained dynamics in filled rubbers and related nanocomposites [23,35]. The MD simulations are consistent with this immobilized-layer picture: the density profile (Figure 2) shows a distinct adsorbed layer with $[eqn]$ and the spatially resolved MSD analysis (Figure 4) shows a mobility gradient ( $[eqn]$ at the interface vs. $[eqn]$ in bulk) consistent with a glass-transition shift of $[eqn]$ –40 K inferred from NMR [31].

The MD-derived interpolated value $[eqn]$ nm at $[eqn]$ (from Equation (7), with the raw $[eqn]$ point excluded from fitting) is smaller than solid-state NMR estimates (3–5 nm) [6,36] and mechanical characterization of sulfur-crosslinked CB-filled rubber (2–4 nm) [24], reflecting the coarse-grained nature of the Kremer–Grest model (bead diameter $[eqn]$ nm) and the spatial resolution limit of the density binning [7,20,37]. Despite the absolute magnitude difference, the relative trend across filler loadings is physically meaningful. The scaling exponent $[eqn]$ is larger than the de Gennes prediction ( $[eqn]$ ), which assumes isolated fillers; the superlinear scaling reflects increasing bound-layer overlap at higher loadings. Slope and exclusion sensitivity checks are summarized in Supplementary Materials Section S3.

Guth–Gold above percolation: The MD-corrected Guth–Gold baseline (replacing $[eqn]$ with $[eqn]$ ) should be understood as a structural prior, not a validation of hydrodynamic theory at $[eqn]$ [10,38,39]. The $[eqn]$ correction captures bound rubber and part of the network stiffening in an effective manner, but it is not sufficient by itself (overall MAPE $[eqn]$ in this dataset). The MD simulations capture only bound rubber around smooth spherical fillers; occluded rubber from aggregate voids [40] requires explicit aggregate geometries. For high-structure fillers, network-explicit constitutive models (e.g., Huber–Vilgis [41]) are still needed, with MD-derived $[eqn]$ serving as transferable structural input.

Filler geometry and carbon black structure: An important limitation concerns the geometric mismatch between MD-simulated smooth spheres and real carbon black aggregates. Industrial carbon blacks are characterized by their structure—the degree of particle aggregation quantified by dibutyl phthalate (DBP) absorption [10,40]. High-structure grades (e.g., N110, DBP $[eqn]$ cm^3^/100 g) form branched aggregates that trap substantial occluded rubber, while low-structure grades (e.g., N990, DBP $[eqn]$ cm^3^/100 g) behave more like isolated spheres. The MD simulations using smooth spheres correspond most closely to the low-structure limit. The empirical $[eqn]$ correction therefore absorbs not only bound rubber but implicitly captures part of the occluded rubber effect that depends on filler structure. This geometric limitation explains why the present scaling exponent ( $[eqn]$ ) differs from theoretical predictions for isolated spheres ( $[eqn]$ – $[eqn]$ ): the superlinear scaling reflects both bound-layer overlap and the aggregate-structure effects that the spherical model treats phenomenologically. Extension to structure-specific MD using explicit aggregate geometries (following Wang et al. [40]) would enable separation of these contributions, but requires substantially larger simulation systems beyond current computational scope.

3.4. Temperature Robustness

A constitutive model intended for process simulation must maintain accuracy across the relevant temperature window. DMA measurements at three temperatures ( $[eqn]$ , $[eqn]$ , $[eqn]$ °C) ( $[eqn]$ invariance has been demonstrated within $[eqn]$ to $[eqn]$ °C; behavior near $[eqn]$ or above the degradation onset has not been tested and may differ) with TTS master curves spanning five frequency decades reveal a key decoupling between kinetic and structural parameters.

The relaxation time $[eqn]$ follows Arrhenius temperature dependence ( $[eqn]$ : $[eqn]$ s; $[eqn]$ kJ/mol), consistent with thermally activated segmental hopping over energy barriers. In contrast, the fractional order $[eqn]$ remains temperature-invariant ( $[eqn]$ , coefficient of variation $[eqn]$ ; $[eqn]$ , MAPE $[eqn]$ at each temperature). This invariance is physically expected: $[eqn]$ encodes the breadth of the relaxation spectrum, which is determined by structural heterogeneity (chain length distribution, entanglement topology, filler proximity) rather than thermal activation energy [16,17].

The $[eqn]$ / $[eqn]$ decoupling has practical significance for process simulation. In extrusion modeling, for instance, material experiences temperature gradients of 30–80 °C across the die. With $[eqn]$ held constant, only $[eqn]$ requires updating via the Arrhenius equation, reducing the constitutive update to a single-parameter adjustment per temperature step. This Arrhenius extrapolation is appropriate only well above $[eqn]$ , as employed here. Compact temperature and bridge consistency results are reported in Supplementary Materials Sections S2 and S5.

3.5. Filler Loading and Scaling Validation

This subsection addresses the second inference layer (cross-scale transfer): whether MD-derived structural descriptors improve reinforcement prediction. It is intentionally separated from the constitutive model-selection question above (single-phase vs. dual-dynamics).

To rigorously test the interfacial empirical scaling relation (Equation (3)), MD simulations were conducted at six filler volume fractions: $[eqn]$ %. Bound-layer thickness ( $[eqn]$ ) was extracted from density profiles $[eqn]$ using cluster-based peak detection: adjacent bins exceeding 10% above bulk density were grouped into clusters, and $[eqn]$ was averaged across clusters per filler. The $[eqn]$ % system yielded anomalous results (zero detected clusters, elevated $[eqn]$ vs. ∼1.95 at neighboring loadings) and was excluded from the fit; the remaining five concentrations ( $[eqn]$ %) spanned the semi-dilute to concentrated regime. The power-law fit (Equation (7); $[eqn]$ , $[eqn]$ ) yielded $[eqn]$ (Figure 5). The superlinear exponent reflected increasing layer interaction at higher loadings, consistent with the onset of filler network percolation [6]. Given the moderate $[eqn]$ and wide bootstrap confidence interval, this exponent was used here as a calibrated interpolation rule over $[eqn]$ –30% (i.e., $[eqn]$ – $[eqn]$ in the equations), not as a broad-universality claim. Figure 7 visualizes the structural progression: as $[eqn]$ increases from 5% to 20%, inter-filler spacing decreases from ∼ $[eqn]$ to ∼ $[eqn]$ , enhancing chain confinement.

The MD-structural baseline for reinforcement is defined as:

$$[eqn]$$

where $[eqn]$ is obtained from MD and R is the filler radius. On the frequency-resolved dataset ( $[eqn]$ ), this baseline captures first-order loading trends but leaves substantial error (overall MAPE $[eqn]$ ), especially in the high-loading regime (MAPE $[eqn]$ ), indicating that structural amplification alone cannot close the MD–DMA dynamics gap.

To bridge this gap, multiplicative corrections to $[eqn]$ are fit using nested blocked-CV model selection. This prediction task evaluates transferability of MD structural priors and does not re-adjudicate constitutive model preference. A single global correction remains unstable across regimes (overall MAPE $[eqn]$ ). Introducing regime partitioning ( $[eqn]$ vs. $[eqn]$ ) with regime-specific basis functions yields a substantially improved model with overall MAPE $[eqn]$ , with low/mid and high loading errors of $[eqn]$ and $[eqn]$ , respectively (Table 3). Because all folds are drawn from the same $[eqn]$ design, blocked-CV here quantifies interpolation robustness within this design space; extrapolation to unseen temperatures, loadings, or filler morphologies remains out of scope (Supplementary Materials Section S4).

ϕeff as a Transferable Structural Constraint

These results support a hierarchical handshaking strategy: MD provides physically grounded structural constraints ( $[eqn]$ ), while mesoscale correction terms absorb regime-dependent relaxation physics unresolved by atomistic timescales. The framework therefore avoids an invalid “direct MD $[eqn]$ ” claim and instead delivers a traceable parameter-passing bridge from interfacial structure to macroscopic reinforcement.

3.6. Scope and Limitations

While the hierarchical MD → DMA → PTT pipeline demonstrates quantitative consistency on the EPDM/CB and PC/ABS systems studied here, the framework involves assumptions at each scale that constrain its current scope. The following points delineate these boundaries to guide future extensions and prevent over-interpretation of the present results.

MD model fidelity: The CG simulations use short chains ( $[eqn]$ , $[eqn]$ ) and smooth spherical fillers, which do not capture entanglement effects ( $[eqn]$ in experimental EPDM) or the complex aggregate morphology of industrial carbon blacks (characterized by DBP absorption). The quantities $[eqn]$ and $[eqn]$ are therefore structural descriptors valid within the Kremer–Grest model; their absolute values should not be equated with experimentally measured interphase thicknesses (2–5 nm from NMR). The relative trends across filler loadings, however, are physically meaningful.Statistical power: The $[eqn]$ scaling relation (Equation (7)) is fitted to only five loading points with $[eqn]$ ; the exponent carries a wide 95% CI (0.04–1.40). This relation serves as a calibrated interpolation within $[eqn]$ – $[eqn]$ , not as a universal exponent. Similarly, AICc-based model selection on the EPDM dataset ( $[eqn]$ ) reflects limited statistical power; the single-phase preference does not imply that two-phase dynamics are absent.Constitutive framework: The fractional Maxwell model is restricted to the linear viscoelastic regime ( $[eqn]$ ) and does not capture large-deformation effects, glass transitions, or nonlinear relaxation. The parallel matrix + network decomposition assumes linear phase superposition; nonlinear interactions between phases are not modeled.Generalizability: The framework has been validated on EPDM/CB and PC/ABS. Extension to other polymer–filler chemistries (silica, nanotubes, graphene), other polymer types (thermoplastic elastomers, thermosets), or extreme temperature ranges (near $[eqn]$ or above degradation onset) requires system-specific re-calibration. Temperature robustness of $[eqn]$ invariance has been demonstrated only within $[eqn]$ to $[eqn]$ + 130 °C.Time-scale bridge: The TTS-based frequency bridge from MD (∼ $[eqn]$ rad/s) to DMA (∼ $[eqn]$ rad/s) assumes thermorheological simplicity, which may be violated by local relaxation distributions in filled systems. The ∼1-decade agreement between $[eqn]$ (MD) and $[eqn]$ (DMA) represents an order-of-magnitude consistency, not rigorous quantitative validation.

3.7. PTT Calibration

To demonstrate the practical utility of the extracted linear viscoelastic parameters, DMA-derived relaxation spectra were used to constrain a Phan-Thien–Tanner (PTT) model for extrusion simulation (Supplementary Materials Section S6). In the unconstrained calibration, PTT parameters ( $[eqn]$ , $[eqn]$ , $[eqn]$ ) were optimized freely against die-swell data, yielding MAPE $[eqn]$ ( $[eqn]$ production batches). When the linear-regime parameters ( $[eqn]$ and $[eqn]$ ) were fixed to DMA values, the remaining nonlinear parameters were determined with substantially reduced ambiguity, reducing die-swell MAPE to $[eqn]$ . The single-mode PTT calibration was demonstrated for one die geometry (L/D = 10, $[eqn]$ = 3 mm) and three production batches. Extension to complex geometries with multiple relaxation times may require multi-mode PTT or pom-pom models. The $[eqn]$ die-swell accuracy should be interpreted as a proof-of-concept rather than a universal industrial prediction. The improvement from fixing $[eqn]$ and $[eqn]$ from LVE was a reduction in calibration uncertainty, not a first-principles prediction of nonlinear behavior.

Important caveats: This 87% error reduction demonstrates a calibration benefit—constraining the linear manifold reduces the dimensionality of the nonlinear optimization problem—not a claim that linear rheology directly predicts nonlinear behavior. The PTT model used here is a single-mode approximation; multi-mode PTT or pom-pom models may be necessary for more complex flow geometries. Additionally, the die-swell comparison is based on a single die geometry (L/D $[eqn]$ , $[eqn]$ mm, $[eqn]$ –375 s^−1^); generalization to other geometries and shear rates requires further validation.

3.8. Time-Scale Bridge

The frequency gap between MD ( $[eqn]$ – $[eqn]$ rad/s) and DMA ( $[eqn]$ –63 rad/s) spans ∼11 orders of magnitude. Applying Arrhenius-based TTS with an independently determined activation energy ( $[eqn]$ kJ/mol from DMA measurements in this work; see Supplementary Materials Sections S2 and S5) reduces this gap to ∼1.0 decade (Figure 8). Crucially, the TTS-shifted MD relaxation time ( $[eqn]$ s) agrees with the DMA-derived fractional Maxwell value ( $[eqn]$ s, see Table 1) within a factor of 1.1, supporting consistency between the two measurements within the current coarse-grained and TTS assumptions. TTS assumes thermorheological simplicity; in filled elastomers, interfacial confinement may generate local relaxation distributions not captured by a global shift. The observed agreement should be interpreted as order-of-magnitude consistency, not rigorous quantitative validation. The residual ∼1.0-decade gap is attributed to the CG model’s smoothed potential energy landscape [12,14]. Furthermore, a $[eqn]$ variation in $[eqn]$ shifts $[eqn]$ by approximately $[eqn]$ , confirming that the ∼1-decade residual is within the propagated uncertainty. Compact TTS and Green–Kubo support are provided in Supplementary Materials Section S5.

4. Conclusions

A multi-scale characterization approach for polymer nanocomposites is developed in this study, investigating the extent to which interfacial-constraint effects account for low-frequency power-law behavior through distinct network and matrix contributions. Four core findings emerge:

Methodological note: A key contribution is the explicit model-selection workflow (AICc + blocked-CV/nested blocked-CV), which separates mechanistic plausibility from statistical necessity and avoids over-interpreting high in-sample $[eqn]$ alone.

(1) $[eqn]$ as a transferable structural descriptor: MD simulations ( $[eqn]$ , Rouse-to-crossover regime, used as structural probes) yield a bound rubber empirical scaling relation $[eqn]$ nm ( $[eqn]$ , $[eqn]$ –30% excluding 20%). This descriptor provides a physically anchored structural prior ( $[eqn]$ ), but does not by itself close the dynamics gap (baseline MAPE $[eqn]$ ). A nested-CV regime-partitioned bridge model built on this prior reduces overall error to $[eqn]$ , with blocked-CV MAPE $[eqn]$ , indicating that MD contributes primarily through constrained parameter passing rather than direct macro-rheology prediction.

(2) Mechanistic partitioning of non-terminal scaling. The observed $[eqn]$ can be decomposed into a matrix phase ( $[eqn]$ ) and network phase ( $[eqn]$ ), with $[eqn]$ . On the EPDM dataset ( $[eqn]$ ), the single-phase model is AICc-preferred ( $[eqn]$ vs. $[eqn]$ ), reflecting limited statistical power rather than model invalidity: the PC/ABS reference system ( $[eqn]$ ) shows dual-phase preference when data are abundant (Table 1). The dual-dynamics decomposition is thus a physically motivated hypothesis consistent with MD-derived interfacial structure, retained for mechanistic insight while acknowledging identifiability constraints at limited sample sizes, and should not be confused with a confirmed fact.

(3) Hierarchical parameter passing for process simulation. The multi-scale workflow culminates in a traceable parameter chain: MD-derived $[eqn]$ (structural prior) constrains reinforcement prediction (finding 1); DMA-derived linear spectra ( $[eqn]$ , $[eqn]$ ) then constrain nonlinear PTT calibration, reducing die-swell error by 87%. This hierarchical approach—fixing the linear manifold before fitting nonlinear parameters—demonstrates a calibration benefit from physics-informed parameter passing rather than direct cross-regime prediction (Supplementary Materials Section S6).

(4) Spatial gradient of subdiffusive dynamics. Spatially resolved MSD analysis reveals $[eqn]$ near fillers versus $[eqn]$ in bulk ( $[eqn]$ ), supporting a molecular-level interfacial-constraint interpretation (Figure 4).

The EPDM/CB system served as the primary model system; the PC/ABS reference system ( $[eqn]$ , $[eqn]$ ) validated the constitutive framework in the unfilled limit, provided the Green–Kubo TTS benchmark (factor-of-1.1 agreement with DMA), and demonstrated that dual-phase AICc preference emerges when statistical power is sufficient—supporting a cautious interpretation of the EPDM single-phase result.

Limitations: First, the CG model uses short chains ( $[eqn]$ , $[eqn]$ ) in the Rouse-to-crossover regime, suitable for extracting local structural quantities ( $[eqn]$ ) but not for simulating entanglement-dominated dynamics. Second, the empirical scaling relation was calibrated on spherical fillers; extension to high-structure carbon blacks (e.g., N990 vs. N110) or anisotropic fillers (clay platelets, carbon nanotubes) requires geometry-specific MD. Third, the dual-dynamics decomposition relied on a physical hypothesis whose parameters were not independently constrained by the EPDM dataset; it should be tested against materials with independently measured network and matrix contributions (e.g., using selective extraction or double-network experiments). Fourth, the bridge model was trained and validated on a limited harmonized design (three temperatures, three nonzero loadings, shared overlap-frequency window); current blocked-CV results support in-range interpolation reliability, not broad extrapolation.

Future work will extend the TTS bridge to EPDM, perform multi-temperature MD for activation energy validation, and investigate anisotropic fillers and LAOS regimes.