Entanglement Characteristic Time from Complex Moduli via i-Rheo GT

Tassieri et al. have introduced a novel rheological tool called “i-Rheo GT” that allows the evaluation of the frequency-dependent materials’ linear viscoelastic properties from a direct Fourier transform of the time-dependent relaxation modulus G(t), without artifacts. They adopted i-Rheo GT to exploit the information embedded in G(t) derived from molecular dynamics simulations of atomistic and quasi-atomistic models, and they estimated the polymers’ entanglement characteristic time (τe) from the crossover point of the moduli at intermediate times, which had never been possible before because of the poor fitting performance, at short time scales, of the commonly used generalized Maxwell models. Here, we highlight that the values of τe reported by Tassieri et al. are significantly different (i.e., an order of magnitude smaller) from those reported in the literature, obtained from either experiments or molecular dynamics simulations of different observables. In this work, we demonstrate that consistent values of τe can be achieved if the initial values of G(t), i.e., those governed by the bond-oscillation dynamics, are discarded. These findings have been corroborated by adopting i-Rheo GT to Fourier transform the outcomes of three different molecular dynamics simulations based on the following three models: a dissipative particle dynamics model, a Kremer–Grest model, and an atomistic polyethylene model. Moreover, we have investigated the variations of τe as function of (i) the ‘cadence’ at which G(t) is evaluated, (ii) the spring constant of the atomic bone, and (iii) the initial value of the shear relaxation modulus G(O). The ensemble of these results confirms the effectiveness of i-Rheo GT and provide new insights into the interpretation of molecular dynamics simulations for a better understanding of polymer dynamics.


Introduction
The stress relaxation function, G(t), of the entangled polymer melts in the linear regime is one of the most fundamental and important dynamic properties in the tube theory. Since the efficient multiple-tau correlator method was proposed [1,2], the linear stress relaxation data with good quality become accessible from the molecular dynamics (MD) simulations. In order to reveal the linear viscoelastic behaviors, G(t) is generally transformed into the storage and loss moduli G (ω) and G (ω) by fitting G(t) with the multi-mode Maxwell function. However, this procedure may suffer from some uncertainties and artifacts. Recently, Tassieri et al. proposed an effective analytical tool called "i-Rheo GT" to convert the time-dependent G(t) into the frequency-dependent complex moduli, G (ω) and G (ω), via the oversampling technique, wherein the raw G(t) data are interpolated to produce oversampled data [3]. Since i-Rheo GT does not require any preconceived models, the artifacts, especially in the high-frequency regime which may exist in the early works, could be ruled out. For example, for the generic Kremer-Grest (KG) model, a crossover of G (ω) and G (ω) at the high frequency like in the vicinity of 1/τ e with τ e = N e 2 ζb 2 / 3π 2 k B T where τ e is the entanglement time, N e is the entanglement length, ζ is the monomeric friction coefficient, b is the Kuhn segment length, and T is the temperature, is observed for the first time [3]. Thus, a specific problem/challenge, that is, the qualitative disagreement between the fitted G (ω) and G (ω) from the multi-mode Maxwell function fittings of G(t) with the experimental data is overcome, indicating that a direct determination of τ e from the crossover becomes possible with i-Rheo GT for the KG model. Furthermore, when applying i-Rheo GT to the atomistic polyethylene (PE) model, a dependence of τ e on the molecular weight, M, (or the chain length, N) is found, for the first time, to be in a stretched exponential form which approaches a constant value from above. Note that the convergence performance of τ e (M) (or τ e (N)), which approaches an asymptotic constant value from bigger values with an increase in M, is qualitatively similar with the dependence of N e on the chain length N wherein N e is obtained from the modified estimators defined in the primitive path analyses (PPA) [4]. Since its emergence, i-Rheo GT has been applied in different systems (e.g., ring catenane, polymer-elastomer blends, elastomeric networks, and hierarchical simulation) [5][6][7][8] for the evaluation of the viscoelastic properties and has been implemented in the REPTATE rheology software [9].
However, the τ e -values of both coarse-grained (CG) KG and atomistic PE models estimated from i-Rheo GT seem inconsistent with previously reported results [3]. For the CG KG model, N e calculated from PPA, plateau modulus in the stress relaxation, and the monomer mean-square-displacement (MSD) is about 85 [10][11][12]. Correspondingly, τ e ≈ 10 4 τ LJ where τ LJ is the Lennard-Jones (LJ) time unit. In addition, when fitting G(t) with the slip-links model, τ e ≈ 5800 τ LJ . While τ e estimated from i-Rheo GT is about 671τ LJ , which is approximately one order of magnitude smaller than those estimated by other methods mentioned above. On the side, for the atomistic PE model, the value of τ e ≈ 0.32 ns (at 450 K) obtained from the complex moduli via i-Rheo GT is somewhat smaller than those obtained from MSD, rheology theory mapping, and experiments [3,13]. Especially with the same definition of τ e as Tassieri et al. [3], Szántó et al. [14] utilized the high-frequency rheometer to experimentally determine τ e and obtained a τ e of 3.0 +9.0 −1.8 ns (at 463 K). The lower bound of the experimental τ e is still almost four times larger than that from Tassieri et al. via i-Rheo GT [3]. Apparently, the large discrepancies of τ e estimated from i-Rheo GT and other methods, including simulations and experiments, could not be simply attributed to different physical observables. We note that N e of the KG model, which is estimated from PPA, stress relaxation, and MSD, has reached consistent results [12,13], despite the fact that PE may exhibit some complex effects of the chemical structure on entanglement characters of the melts [1,15]. Since τ e is one of the most fundamental parameters in the tube theory and all viscoelastic polymers, it is meaningful to clarify the discrepancy of τ e calculated from i-Rheo GT and other methods and find an approach to ensure consistency of τ e estimated from different methods. In this work, for a comprehensive study of the above important and fundamental issues in polymer science, we conducted the hierarchical molecular dynamics simulations at different scales and investigated the linear rheology in the entangled systems. Three typical model systems are considered, including models with high computational efficiency to obtain strongly entangled polymer melts at a reasonable computational cost and models for which extensive detailed studies and relevant reference data are available. The simulated models include the entangled dissipative particle dynamics (DPD) model, generic CG KG model [16], and atomistic PE model [17]. For the computationally efficient DPD model [18][19][20][21], the polymer melts with different chain lengths are simulated. For KG and PE models, only well-entangled melts are simulated.

DPD Simulation
In DPD, all beads interact with a pairwise conservative force (F C ), dissipative force (F D ), and random force (F R ). For example, the forces between bead i and j located at the coordinates of r i and r j , respectively, are given below: a ij is the repulsion constant, r c is the cutoff radius, γ is the friction parameter, v i is the velocity of bead i, σ is the thermal noise amplitude, ζ ij is the Gaussian random variable, and ∆t is the simulation timestep. Note γ and σ are connected according to the fluctuation-dissipation relaxation: σ 2 = 2γk B T with k B being the Boltzmann constant and T being the temperature. The sequential beads in the same chain are connected by the harmonic springs: with K B being spring constant and l being equilibrium bond length. Additionally, a bond bending potential taken as U bend (θ) = K θ (1 + cosθ) where K θ is the bending stiffness constant, and θ is the angle of consecutive bond vectors, is added to increase the stiffness of chain. The quantities are taken in the DPD reduced units, namely, the units of length, mass, and energy are set as r c , m, and k B T respectively. The corresponding time unit is τ = r c ( m k B T ) 1/2 . It is worth stressing that, different from the standard DPD, the conservative force used here is properly increased to prevent chain crossings from showing entanglement effects [22]. The parameters of the DPD model used here are the same with refs [18][19][20][21] in which a ij = 200, r c = 1, γ = 4.5, ∆t = 0.012, k B T = 1, K B = 400, l = 0.95, K θ = 2, and the density ρ = 1 in DPD reduced units.

KG Simulation
In the standard KG model [16], the bead-spring polymer chains are represented as a sequence of beads connected by the finitely extensible nonlinear elastic (FENE) spring. All beads interact by the purely repulsive Lennard-Jones (LJ) potential: where r is the distance between two beads, ε is the LJ energy scale, and σ is the LJ radius. The cutoff distance r c is set as 2 1/6 σ. The FENE potential is described as: where k is the spring constant and equals 30ε/σ 2 , R 0 = 1.5σ is the maximum bond length. All quantities are expressed in the LJ reduced units. σ, ε, and the bead mass m are equal to 1. The temperature is maintained at The corresponding time unit is τ LJ = mσ 2 /k B T 0.5 . For our MD simulations of the KG model, the simulation parameters are the same as those of Tassieri et al. [3], except for the thermostat and timestep. We note that compared to the Langevin thermostat used in the work of Tassieri et al. [3], the DPD thermostat could conserve both global and local momentums. Therefore, in the present work, the DPD thermostat is used (instead of the Langevin thermostat), and the timestep is ∆t = 0.006τ LJ (other than 0.012τ LJ ).

PE Simulation
In the case of PE, the well-known TraPPE united-atom (UA) force field is used [17]. Each terminal methyl CH 3 group and internal methylene CH 2 group is set as a single interacting point. The bonded interactions include bond stretching, bond angle bending, and torsional rotation. In the original description of the TraPPE force field, the bond lengths are fixed at 1.54. In this work, the harmonic bond potential of U = K(r − r 0 ) 2 with the spring constant of K = 120 kcal/mol and the equilibrium bond length of r 0 = 1.54 is taken. The bond angle bending potential is described by a harmonic potential: where k θ is the force constant with k θ /k B = 62,500 K rad −2 and the equilibrium angle θ 0 is 114 • . The torsional interaction is governed by the OPLS (optimized potentials for liquid simulations) potential: where the force constants c 1 , c 2 and c 3 are 355.03 k B K, −68.19 k B K, and 791.32 k B K, respectively. The non-bonded interactions are described by LJ potentials: where the parameters of CH 3 and CH 2 groups are ε CH 3 /k B = 104K, ε CH 2 /k B = 46K, σ CH 3 = 3.91 and σ CH 2 = 3.95, respectively. The interactions between unlike particles are computed by Lorentz-Berthelot mixing rules. The cutoff distance of LJ interactions is 14.

DPD
We first focus on the DPD model as an alternative to the entanglement molecular model for MD simulations, which, as noted above, has become a promising tool for studying the long-time dynamics of polymeric liquids. In Figure 1, we present the stress relaxation for the DPD polymer chains with lengths N = 129, 161, and 321. The inset shows the bond oscillations at the early time in a semi-logarithmic plot. The stress is calculated at every timestep via the multiple-tau correlator algorithm [1,2]. The stress relaxation behaviors of the DPD model are qualitatively similar to those of the KG model. Note that as indicated by the horizontal dash line in Figure 1b, the normalized stress relaxation, G(t)t 0.5 , shows a plateau as expected from the Rouse model for t after the bond length relaxation. However, the slope of G(t)t 0.5 for the standard KG model at 1 < t < 30 is about −0.13, which some deviates from the prediction of the Rouse model [1]. In this regard, the DPD model seems more suitable to describe stress relaxation. The resulting viscoelastic moduli data via i-Rheo GT by feeding the absolute values of ( ) into i-Rheo GT, as suggested in Ref 3, are shown in Figure 2. Note that the quality of the curves of = 321 are not as good as those of = 129 and 161 since the total simulation times of the latter ones are almost 200 times longer than their disentanglement time while the simulation time of the former one is less than a decade of its disentanglement time. Nevertheless, Figure 2 distinctly indicates a crossover point of the storage elas- The resulting viscoelastic moduli data via i-Rheo GT by feeding the absolute values of G(t) into i-Rheo GT, as suggested in Ref. [3], are shown in Figure 2. Note that the quality of the curves of N = 321 are not as good as those of N = 129 and 161 since the total simulation times of the latter ones are almost 200 times longer than their disentanglement time while the simulation time of the former one is less than a decade of its disentanglement time. Nevertheless, Figure 2 distinctly indicates a crossover point of the storage elastic modulus G and the loss of viscous modulus G at ω cross in the intermediate frequency scale for all systems. Furthermore, as presented in the inset of Figure 2, tan δ(10ω cross ) ≈ 2, which is in agreement with the theory and experiments [3]. We could further estimate τ e, GT from 1/ω cross where the subscript GT in τ e,GT indicates the entanglement time is obtained from i-Rheo GT. τ e,GT are about 68, 67, and 44 for N = 129, 161, and 321, respectively, much smaller (beyond an order of magnitude) than the entanglement time estimated from MSD, which is indicated as τ * e ≈ 810 or τ e = 9τ * e π ≈ 2300 given in Figure S1 in SI. This finding is similar to the results of KG and PE models reported by Tassieri et al. [3], as mentioned above. Additionally, the smaller value τ e, GT of N = 321 than those of the other two shorter counterparts indicates that τ e, GT relies on the chain length, as further evidenced by the broad chain length window in Figure S2 in SI, which qualitatively agrees with the results of PE by Tassieri et al. [3] as well. The resulting viscoelastic moduli data via i-Rheo GT by feeding the absolute values of ( ) into i-Rheo GT, as suggested in Ref 3, are shown in Figure 2. Note that the quality of the curves of = 321 are not as good as those of = 129 and 161 since the total simulation times of the latter ones are almost 200 times longer than their disentanglement time while the simulation time of the former one is less than a decade of its disentanglement time. Nevertheless, Figure 2 distinctly indicates a crossover point of the storage elastic modulus and the loss of viscous modulus at in the intermediate frequency scale for all systems. Furthermore, as presented in the inset of Figure 2, tan (10 ) ≈ 2, which is in agreement with the theory and experiments [3]. We could further estimate , from 1/ where the subscript in , indicates the entanglement time is obtained from i-Rheo GT. , are about 68, 67, and 44 for = 129, 161, and 321, respectively, much smaller (beyond an order of magnitude) than the entanglement time estimated from MSD, which is indicated as * ≈ 810 or = * ≈ 2300 given in Figure S1 in SI. This finding is similar to the results of KG and PE models reported by Tassieri et al. [3], as mentioned above. Additionally, the smaller value , of = 321 than those of the other two shorter counterparts indicates that , relies on the chain length, as further evidenced by the broad chain length window in Figure S2 in SI, which qualitatively agrees with the results of PE by Tassieri et al. [3] as well.  Experimentally, the highest sampling frequency of the rheometers is limited in the mega-hertz range and is usually quite smaller than that taken in molecular simulations, wherein the stress is often sampled every timestep, and the corresponding sampling frequency is in the tera-hertz range [3,14]. Therefore, it is interesting to check the effects of the limited dynamic window as accessed by experiment instruments on the Fourier transformation of G(t) via i-Rheo GT. We now discard the short time data of G(t), corresponding to the regime of the bond length relaxation, during the transformation from G(t) to G (ω) and G (ω) just as early works in experiments and simulations. [1,9] We note that for the DPD model used here, the bond length relaxation regime would be "naturally" discarded if the stress is calculated every hundred timesteps since the corresponding sampling time interval is 1.2(= 100∆t)τ, with ∆t being the simulation timestep and τ being DPD time unit, which is larger than the relaxation time of the bond oscillations, i.e., about 1τ as indicated by the inset of Figure 1. Taking N = 129 as an example, as demonstrated in Figure 3a, the resulting G(t) curves from stress data calculated every timestep and every hundred timesteps (denoted as "every 1" and "every 100", respectively) fall practically on top of each other, though G(t) of "every 100" is with poorer quality. As an alternative route to discarding the high-frequency domain of bond length relaxation, we can "artificially" Polymers 2022, 14, 5208 6 of 11 discard the G(t) data at an early time, i.e., G(t < 1.2), while stress and the original G(t) are calculated at every timestep, which is denoted as "discard" in what follows. As expected, in Figure 3b, the viscoelastic spectrums from G(t) of "every 100" nearly overlap with those of "discard". Since the moduli of "discard" are of better quality than the moduli of "every 100", for a better understanding of the limited dynamic window problem, we would focus on comparisons of the moduli of "every 1" and "discard" later. As indicated by the black and green lines, though the G curves of "every 1" and "discard" lie on top of each other for ω < 0.02, G of the former is smaller than that of the latter apart from ω at the terminal scale, which G is in the same magnitude as both "every 1" and "discard". Consequently, the cross point of the moduli of "discard" moves to a lower frequency regime. Interestingly, the resulting τ e is about 820, which is approximately an order of magnitude larger than that estimated from the moduli of "every 1" but very close to τ * e ≈ 810 estimated from MSD, as mentioned above. In this regard, the result of discarding the early bond oscillation regime appears to be more consistent with results from other common physical observables, independent of the "naturally" or "artificially" discarding approach used. As we see below, this also holds true for both KG and PE models. A consistent result could be anticipated when i-Rheo GT is used to analyze the experimental data since the time scale corresponding to the bond length relaxation is hardly accessible and is usually not included in the rheological experiments. Additionally, as demonstrated by the black and red lines in Figure 3, it seems very tricky that for i-Rheo GT, the sampling time interval plays a critical role in the determination of τ e , even though the time interval is still smaller than τ e by tens of times. We note that the time scale corresponding to the bond length relaxation is beyond the consideration of tube theory on the viscoelastic properties of the polymers [15], while in real polymer melts, there exists an intricate interplay between bonded and non-bonded interactions. The method of i-Rheo GT to process the data of the coupling interactions between the different dynamics mechanisms at different time regimes might be the origin of the observed deviation. This, to some extent, is qualitatively similar to the case in which single-chain models obey the stress-optical law rather well for times longer than the bond length relaxation time, whereas for many-chain models, it is valid on slightly longer time scales as a consequence of complicated coupling of the interactions. [23] Nevertheless, we note that as shown in the inset, for the moduli discarding the bond oscillation regime, tan δ(10ω cross ) ≈ 3.3, which disagrees with the theoretical and experimental expectation.  The absolute values of ( ) are input into i-Rheo GT, and the value of (0) is 43.7. The inset in panel (b) presents the loss tangent tan vs. frequency. The legends of "every 1" and "every 100" indicate the stress is calculated using every first and hundredth point, respectively. The legend of "discard" means the stress is calculated at every timestep while the data at early time ( < 1.2 ) is artificially discarded.
The results discussed above show whether the inclusion of the bond length relaxation regime or not plays an important role in the exact determination of relevant parameters in rheology via i-Rheo GT. To gain a detailed insight into the influences of the bond length relaxation regime on the determination of , we do further analysis on the "discard" data. For simplicity, we tune the spring constant, , of bonded beads and artificially replace the data of ( < 1.2) by the ones from DPD simulations with different spring constants while the data of ( ≥ 1.2) remain unchanged. Note that a smaller would result in the crossability of chains and a larger requires a smaller timestep [22]. Here, the stress is calculated at every time of 0.012 though the timestep is different for differ- presents the loss tangent tan δ vs. frequency. The legends of "every 1" and "every 100" indicate the stress is calculated using every first and hundredth point, respectively. The legend of "discard" means the stress is calculated at every timestep while the data at early time (t < 1.2τ) is artificially discarded.
The results discussed above show whether the inclusion of the bond length relaxation regime or not plays an important role in the exact determination of relevant parameters in rheology via i-Rheo GT. To gain a detailed insight into the influences of the bond length relaxation regime on the determination of τ e , we do further analysis on the "discard" data. For simplicity, we tune the spring constant, K B , of bonded beads and artificially replace the data of G(t < 1.2) by the ones from DPD simulations with different spring constants while the data of G(t ≥ 1.2) remain unchanged. Note that a smaller K B would result in the crossability of chains and a larger K B requires a smaller timestep [22]. Here, the stress is calculated at every time of 0.012τ though the timestep is different for different K B . Taking N = 129 as an example again, the early oscillations of bond relaxation for different spring constants, which range from 100 to 1600, are presented in Figure 4a. The height and number of peaks increase with the increase of the spring constant, indicating that the oscillations become more intense and complicated. Particularly, the value of G(0) increases quickly with the spring constant, as demonstrated in the inset of Figure 4a. Figure 4b presents the corresponding viscoelastic moduli. Although the G (ω < 0.3) curves of different spring constants lie on top of each other, G generally increase with the increase in the spring constant for ω > 10 −3 . The different responses of G and G to the spring constant are similar to the responses of G and G to the sampling frequency as mentioned above. Hence, the values of ω for the cross point become smaller with the increase in the spring constant. The resulting τ e depends strongly on the spring constant K B as presented in Figure 4c, even though G(t ≥ 1.2) are the same for all systems with different spring constants. Additionally, we would like to point out that τ e even shows an exponential dependence on G(0), which can be described by a power law of τ e ∝ G(0) 0.8 as demonstrated in Figure 4d. These findings convincingly demonstrate that the dynamics at the shorter time scale could intensively affect the calculation of τ e located at a longer (at least by dozens of times) scale when applying i-Rheo GT to convert G(t) into G (ω) and G (ω), while the same τ e is expected to be obtained if we "artificially" discard the data of G(t < 1.2). Altogether, our DPD results reveal the limits of applicability of i-Rheo GT to the exact determination of relevant parameters in rheology. The reason for this apparent failure lies in the method of i-Rheo GT to process the data of the coupling interactions between the different dynamic mechanisms at different time regimes.

KG
To elucidate the above issues further, we consider KG model as another test case and check the consistency of estimated from i-Rheo GT with those available reference data from other methods. The chain length is = 350, which is the same as the longest chain of Tassieri et al. [3] Note that the stress is calculated every two timesteps here to keep the

KG
To elucidate the above issues further, we consider KG model as another test case and check the consistency of τ e estimated from i-Rheo GT with those available reference data from other methods. The chain length is N = 350, which is the same as the longest chain of Tassieri et al. [3] Note that the stress is calculated every two timesteps here to keep the same sampling time interval with the work of Tassieri et al. [3] since the sampling frequency seems to play a role in i-Rheo GT analysis as mentioned above. The resulting viscoelastic moduli of the KG model are shown in Figure 5. The curves of "every 2" are generally the same as those in the work of Tassieri et al. [3] despite that the quality of the curves at the middle scale (ω ≈ 10 −5 − 10 −4 /τ LJ ) in this work is poorer. An approximate range of τ e,GT estimated from tan δ(ω cross ) = 1 (as shown in the inset of Figure 5) is 600 ∼ 900, which is comparable with the result of 671 reported by Tassieri et al. [3] but is quite smaller than those calculated by other methods as mentioned above. Additionally, the characteristic time for tan δ = 2 is about 25, which is smaller than τ e,GT by several tens of times. Consequently, the expectation of tan δ(10ω cross ) = 2 is not satisfied from the results of "every 2" for the KG model (this expectation is met for the DPD model). Like in the DPD case, we "artificially" discard the G(t) data at an early time, i.e., G(t < 1.2). As indicated by the green curves in Figure 5, the behaviors of G and G are similar to those of the DPD model described above, namely, the curve of G (ω < 0.02) of "discard" falls on top of that of "every 2" and G is on the top of that of "every 2", except for the terminal scale. τ e,GT estimated from tan δ = 1 of "discard" is τ e,GT ≈ 8400, which is comparable with the values from other estimators as mentioned above, i.e., τ * e estimated from MSD for the KG model is τ * e ≈ 2950 [24], thus τ e = 9τ * e /π ≈ 8450. Nevertheless, the expectation of tan δ(10ω cross ) = 2 is not satisfied from the results of "discard" as well.
Polymers 2022, 14, x FOR PEER REVIEW 9 of 12 Figure 5. Viscoelastic moduli of KG model. The absolute values of ( ) are input into i-Rheo GT, and the value of (0) is 67.7. The inset presents the tan vs. frequency. The legends of "every 2" and "discard" mean that the stress is calculated using every second point, and the data of ( < 1.2) of "every 2" are discarded, respectively.

PE
For a complete understanding of the limits of applicability of i-Rheo GT, we then consider the well-studied atomistic PE model. The MD simulation condition is the same as the work of Tassieri et al. [3] (at 450 K, 1 atm). The stress is calculated at every timestep (Δ = 2 fs). The molecular weight is 5.04 kg/mol, corresponding to = 360 carbon atoms per chain. Figure 6 presents the viscoelastic moduli of PE. Similar to the results of DPD and KG models, the curves of of "every 1" and "discard" lie on top of each other as < 20/ns while the curves of of these two cases overlap with each other only at the terminal regime. , are approximately estimated to be 0.55 ns and 4.2 ns from "every 1" and "discard", respectively. Note that these two , are separately compared with the value calculated from i-Rheo GT (which is 0.37 ns for the molecular weight being 5.04 kg/mol) and the value of ~5.5 ns estimated from MSD in Ref [3]. Consistent with DPD and KG models, , estimated from "discard" here is in better agreement with from other methods, however, , from "every 1" appears to be somewhat smaller. As for tan (10 ), the values are about 2.1 for "every 1"and 4.4 for "discard", which is sim- The inset presents the tan δ vs. frequency. The legends of "every 2" and "discard" mean that the stress is calculated using every second point, and the data of G(t < 1.2) of "every 2" are discarded, respectively.

PE
For a complete understanding of the limits of applicability of i-Rheo GT, we then consider the well-studied atomistic PE model. The MD simulation condition is the same as the work of Tassieri et al. [3] (at 450 K, 1 atm). The stress is calculated at every timestep (∆t = 2 fs). The molecular weight is 5.04 kg/mol, corresponding to N = 360 carbon atoms per chain. Figure 6 presents the viscoelastic moduli of PE. Similar to the results of DPD and KG models, the curves of G of "every 1" and "discard" lie on top of each other as ω < 20/ns while the curves of G of these two cases overlap with each other only at the terminal regime. τ e,GT are approximately estimated to be 0.55 ns and 4.2 ns from "every 1" and "discard", respectively. Note that these two τ e,GT are separately compared with the value calculated from i-Rheo GT (which is 0.37 ns for the molecular weight being 5.04 kg/mol) and the value of~5.5 ns estimated from MSD in Ref. [3]. Consistent with DPD and KG models, τ e,GT estimated from "discard" here is in better agreement with τ e from other methods, however, τ e,GT from "every 1" appears to be somewhat smaller. As for tan δ(10ω cross ), the values are about 2.1 for "every 1"and 4.4 for "discard", which is similar to the case of DPD.
as the work of Tassieri et al. [3] (at 450 K, 1 atm). The stress is calculated at every timestep (Δ = 2 fs). The molecular weight is 5.04 kg/mol, corresponding to = 360 carbon atoms per chain. Figure 6 presents the viscoelastic moduli of PE. Similar to the results of DPD and KG models, the curves of of "every 1" and "discard" lie on top of each other as < 20/ns while the curves of of these two cases overlap with each other only at the terminal regime. , are approximately estimated to be 0.55 ns and 4.2 ns from "every 1" and "discard", respectively. Note that these two , are separately compared with the value calculated from i-Rheo GT (which is 0.37 ns for the molecular weight being 5.04 kg/mol) and the value of ~5.5 ns estimated from MSD in Ref [3]. Consistent with DPD and KG models, , estimated from "discard" here is in better agreement with from other methods, however, , from "every 1" appears to be somewhat smaller. As for tan (10 ), the values are about 2.1 for "every 1"and 4.4 for "discard", which is similar to the case of DPD. Figure 6. Viscoelastic moduli of PE. The absolute values of ( ) are input into i-Rheo GT, and the value of (0) is 14.5 GPa. The inset shows tan vs. frequency. The legends of "every 1" and "discard" mean that the stress is calculated at every timestep, and the data of ( < 10 ns) of "every 1" are discarded, respectively. 5 GPa. The inset shows tan δ vs. frequency. The legends of "every 1" and "discard" mean that the stress is calculated at every timestep, and the data of G(t < 10 −3 ns) of "every 1" are discarded, respectively.
To further elaborate on the influence of the high-frequency dynamics, we gradually discard the data from the highest range and the obtained τ e is presented in Figure 7. For t start < 10 −4 ns, the obtained τ e remains unchanged and is quite smaller than the reported values shown by the lines and symbols on the left-hand side. The obtained τ e increases quickly with the increase of t start before the end of bond length relaxation (10 −4 < t start < 10 −3 ) and increases even faster after the bond length relaxation. Note that the experimentally accessible time range is quite larger than the relaxation time of bond length fluctuation (~10 −3 ns). Generally, the resulting τ e from i-Rheo GT would be comparable with the reported ones, which are obtained from MSD, rheology theory mapping, neutron spin-echo spectroscopy (NSE), and rheology experiments if the highfrequency dynamics is discarded, indicating that the ability to interpret the data at the high-frequency dynamics for i-Rheo GT seems insufficient in the current framework of the tube theory. It would be very interesting to look forward to improving rheology theory and experiment techniques to interpret high-frequency dynamics. To further elaborate on the influence of the high-frequency dynamics, we gradually discard the data from the highest range and the obtained is presented in Figure 7. For < 10 ns, the obtained remains unchanged and is quite smaller than the reported values shown by the lines and symbols on the left-hand side. The obtained increases quickly with the increase of before the end of bond length relaxation (10 < < 10 ) and increases even faster after the bond length relaxation. Note that the experimentally accessible time range is quite larger than the relaxation time of bond length fluctuation (~10 ns). Generally, the resulting from i-Rheo GT would be comparable with the reported ones, which are obtained from MSD, rheology theory mapping, neutron spin-echo spectroscopy (NSE), and rheology experiments if the high-frequency dynamics is discarded, indicating that the ability to interpret the data at the high-frequency dynamics for i-Rheo GT seems insufficient in the current framework of the tube theory. It would be very interesting to look forward to improving rheology theory and experiment techniques to interpret high-frequency dynamics. , which means the data of ( < ) are discarded.
from Table 2 in Ref [13] are shown on left side. The vertical dash line indicates the end of the bond length relaxation.

Conclusions
In conclusion, although i-Rheo GT could transform the stress relaxation from time to frequency domain without preconceived models and over an extremely wide range of frequencies, especially at high scales, the determination of the entanglement time from the Figure 7. Dependence of τ e of PE model on the start time, t start , which means the data of G(t < t start ) are discarded. τ e from Table 2 in Ref. [13] are shown on left side. The vertical dash line indicates the end of the bond length relaxation.

Conclusions
In conclusion, although i-Rheo GT could transform the stress relaxation from time to frequency domain without preconceived models and over an extremely wide range of frequencies, especially at high scales, the determination of the entanglement time from the transformed complex moduli could not offer a consistent result with those estimated from other common methods but depends on the sampling frequency at least for DPD, KG models and atomistic PE model in the current framework of tube theory. Interestingly, by discarding the stress relaxation data at high frequency upon applying i-Rheo GT, the consistency of the resulting entanglement time with those estimated from experiments and simulations is reached. Thus, it is conceivable that the practical approaches for determining entanglement time from G(t) might be the use of i-Rheo GT with discarding the early bond oscillations regime or the fit of G(t) with well-built models such as tube or slip-links models [1,25].
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/polym14235208/s1, Figure S1: Mean-square displacement of the middle monomers for N = 129, 161, and 321. The vertical coordinate is normalized by the power law of t 0.25 . The intersecting time of the time regimes t < τ e and t > τ e is indicated as τ * e ; Figure S2: tan δ vs. frequency ω for the chain length ranging from 64 to 321. The regime of tan δ around 1 is enlarged in the inset [26][27][28].