Tracer dynamics in polymer networks: generalized Langevin description

Tracer diffusion in polymer networks and hydrogels is relevant in biology and technology, while it also constitutes an interesting model process for the dynamics of molecules in fluctuating, heterogeneous soft matter. Here, we study systematically the time-dependent dynamics and (non-Markovian) memory effects of tracers in polymer networks based on (Markovian) implicit-solvent Langevin simulations. In particular, we consider spherical tracer solutes at high dilution in regular, tetrafunctional bead-spring polymer networks, and control the tracer-network Lennard-Jones (LJ) interactions and the polymer density. Based on the analysis of the memory (friction) kernels, we recover the expected long-time transport coefficients, and demonstrate how the short-time tracer dynamics, polymer fluctuations, and the viscoelastic response are interlinked. Further, we fit the characteristic memory modes of the tracers with damped harmonic oscillations and identify LJ contributions, bond vibrations, and slow network relaxations, which enter the kernel with an almost linear scaling with the LJ attractions. This procedure proposes a reduced functional form for the tracer memory, allowing for a convenient inter- and extrapolation of the memory kernels. This leads eventually to highly efficient simulations utilizing the generalized Langevin equation (GLE), in which the polymer network acts as an additional thermal bath with tuneable intensity.

9][30][31][32][33][34] In particular, improved experimental single particle tracking methods, for example, by fluorescence microscopy, have elucidated the stochastic motion and complex physical nature of tracer motion on a wide range of microscopic timescales with high resolution.The goal of these experiments, as well as the numerous accompanying theoretical works, is to learn about the dynamical features of both, the underlying medium as well as the coupled motion of the tracer particles in these usually crowded, disordered and highly fluctuating environments.a) Electronic mail: sebastian.milster@physik.uni-freiburg.deb) Electronic mail: tanja.schilling@physik.uni-freiburg.dec) Electronic mail: joachim.dzubiella@physik.uni-freiburg.de   Occurring phenomena, like anomalous (sub-) diffusion, non-Markovianity, and memory, are often intertwined processes, and far from being well understood.Especially in polymeric systems, such as melts and networks, slow relaxation, memory and subdiffusion are commonplace, [35][36][37][38][39][40][41][42][43][44][45] and can connect to the tracer dynamics for strong coupling. 33,46,479][50][51][52][53][54][55][56] It was demonstrated that such dynamics sensitively depends on microscopic details, such as polymer density, mesh-size (or crosslink density), and in particular the tracer-polymer interactions, 33 defining effective energy landscapes for diffusion. 51,55,57hile the long-time dynamics has been well studied, we are now interested in the physics of the full timescale dependent tracer dynamics in a fluctuating polymer network, and how it is governed ('slaved') by the physical coupling to the polymer network dynamics through interactions.Intuitively, repulsive tracers in dilute networks should be weakly coupled to the polymer fluctuations, while more attracted tracers in dense networks are expected to pick up much more of the dynamic features of the matrix.How the short-time tracer dynamics is controlled by the microscopic interactions in detail is still challenging to categorize. 33 popular tool to model and interpret anomalous dynamics in (macro-) molecular liquids revolves around memory (friction) kernels in the framework of the non-Markovian generalized Langevin equation (GLE), 45,[58][59][60][61][62][63][64][65][66][67][68][69][70] leading, for example, to recent GLE descriptions for anomalous polymer 42,71,72 and tracer dynamics. 475][76][77] Nevertheless, a amenable quantification of the coupled tracer-polymer dynamics for different interaction strengths on the full time-scale is still lack-

Microscopic Model
Generalized Langevin Equation friction memory colored noise tracer solvent ing, and this work aims to fill this gap.
Our strategy and outline is as follows: In Sec.II A we start from a relatively simple model of a regular bead-spring network consisting of crosslinked monomer-resolved beads interacting with spherical tracer particles.We model the solvent implicitly by employing a Markovian Langevin thermostat to the molecular dynamics simulation.From such a 'microscopic' model, we then dynamically coarse-grain to the level of the tracer particles by integrating out the degrees of freedom of the polymeric medium as theoretically outlined in Secs.II B-II F. Particularly, we calculate and analyze the appropriate tracer and polymer memory (friction) kernels in the time and frequency domain.
In Secs.III A and III B we examine systematically the velocity and force autocorrelations as well as the kernels' behavior under the influence of the changing tracer-network interactions, and for different densities.We further present various consistency checks for the long-time dynamics based on the autocorrelations and the memory in Sec.III C. We discuss in detail how the time-dependent friction depends on the coupling to the polymer dynamics (Secs.III D and III E), where we clearly observe large correlations for larger attractions (a few times the thermal energy).
In Sec.III F we finally demonstrate that only a few modes of the kernel are necessary for a massively reduced GLE description of the tracer dynamics, where the magnitude of the friction modes can be described by simple, almost linear scaling laws as a function of the tracer-polymer interaction strength.This enables efficient predictions of the tracer dynamics in different simulation setups in which the polymer plays the role of an additional while tuneable heat bath.

A. Microscopic model
In our microscopic simulations, we consider spherical solute molecules (tracers) diffusing inside a polymeric medium, where the latter is modeled by a bead-spring network similar to previous works, 51,55,78 and as visualized in Fig. 1(a).The polymer network consist of tetrafunctional crosslinkers connecting the polymer chains (ten monomers long) in a regular fashion.
The temporal evolution of the ith particle's velocity is governed by the Markovian Langevin equation, reading where the deterministic forces F int i (t) result from bonded and non-bonded particle interactions, and where −γ 0 mv i (t) + ξ i (t) are the implicit solvent's contributions, namely a viscous drag, −γ 0 mv i (t), and a Gaussian white noise force with zero mean, ⟨ξ i (t)⟩ = 0, and intensity where ν, ν ′ ∈ {x, y, z} are the Cartesian coordinates.Here, γ 0 denotes the friction coefficient in the absence of any other polymer and tracer particles, i.e. for F int i (t) = 0 ∀t.It determines the free equilibrium diffusion via the Einstein relation, D 0 = k B T /(mγ 0 ), where k B T = β −1 is the thermal energy, and it also defines the Brownian timescale τ 0 = γ −1 0 , which we use as reference.All (non-bonded) particles interact via Lennard-Jones potentials, with distance r ij = |r i − r j |, and particle diameter σ, which we chose identical for the network and the tracer beads.The potential strength ε ij depends on the particle pair, resulting in network-network (ε nn ), tracer-tracer (ε tt ), and tracer-network (ε tn ) interactions.
The polymer bonds are modelled by harmonic springs, with strength κ b = 100 k B T /σ 2 and length r 0 = σ.The angle potentials of the network are flexible springs, described by with strength κ a = 1 k B T /rad 2 , and equilibrium angle Ψ 0 mer = π or Ψ 0 x = arccos(−1/3) (tetrahedral angle), dependent on the vertex particle type j, which can be either chain monomer ('mer') or crosslinker ('x').
We fix the bonded network parameters and employ different network-network interaction strengths, βε nn ∈ {0.35, 0.5}, resulting in typical polymer volume fractions of roughly ϕ ∈ {0.16, 0.32} after N pT equilibration. 6,51,55he tracer particles interact mutually repulsively, i.e., the tracer-tracer interaction strength is fixed, βε tt = 0.1.However, we choose a low density, c ≈ 0.005σ −3 , such that the forces F int i acting on a single tracer are known to be practically only governed by the interaction with the polymer network. 51,55More details on the microscopic simulations and a tabulated summary of all parameters are provided in the Appendix A.

B. Generalized Langevin equation
In this work, we reduce the high-dimensional Markovian model for the tracer dynamics by coarse-graining the polymer network (cf.Fig. 1).The Markovian Langevin equation [Eq.(1)] can be then replaced by a non-Markovian integrodifferential equation, namely the generalized Langevin equation (GLE), reading where the memory kernel, K(t − t ′ ), and the correlated fluctuations ζ i (t), comprise all polymer and implicit solvent bath effects of the underlying model.In fact, K carries also some tracer-tracer cross-correlations, which, in principle, can be corrected. 79,80Nonetheless, with the choice of a relatively low tracer concentration (c = 0.005σ −3 ), we may regard K as a single-tracer memory.A comparison with a test simulation with an even lower tracer concentration (c ≈ 0.001σ −3 , not presented in this work) suggests that the effect of the tracertracer interactions on the correlations and the memory is of the order of a few percent.In Eq. ( 5), the components of the fluctuating force ζ iν (t) with ν ∈ {x, y, z} and the memory kernel K(t−t ′ ) are related via the fluctuation-dissipation theorem (FDT), Note that the kernel K(τ ) in the present work does not change over time and depends only on the delay τ = t − t ′ (with arbitrary t and t ′ ).If the memory decays practically instantaneously, K(τ ) → γδ(τ ), we recover the structure of the Markovian Langevin equation with friction −mγv i (t) and delta-correlated fluctuations.More details on the derivation of the GLE and the formal definitions of the quantities appearing in the GLE are summarized in the Appendix B.

C. Relation to the autocorrelation function
By multiplying the GLE [Eq.( 5)] with v i (0) and taking the ensemble average, we obtain the equation of motion for the autocorrelation dynamics Here, the orthogonality of the fluctuating force and the initial velocity is used, namely ⟨ζ iν (t)v iν (0)⟩ = 0, and we introduce the velocity autocorrelation function (VACF) The prefactor results from averaging over the Cartesian coordinates due to the isotropy of the system.Hence, at zero delay it takes the value By taking the Laplace transform of Eq. ( 7), and using the 1D Green-Kubo relation for the diffusion constant D, we find Employing the Einstein relation, we obtain a relation between the effective friction constant and the integral over the memory kernel, This friction constant is compatible with the correct diffusion constant by definition.Further, one can use the relation between Laplace transforms and one-sided Fourier transforms, F(ω) = F(−iω), to define the complex diffusivity The complex diffusivity, which is proportional to the generalized dynamic mobility μ(ω), 60,81,82 is related to the kernel via [c.f.Eq. ( 9)] The one-sided Fourier transform of the memory kernel, K(ω) = K(−iω), provides the frequency spectrum | K(ω)| of the fluctuating forces ζ iν (t) in Eq. ( 5), 83 and can also be interpreted as the frequency-dependent friction entering the viscoelastic response, which we discuss later.
In our analysis of the microscopic system we will also evaluate its force autocorrelation function (FACF), defined as with F(t) = m v(t) from the microscopic velocities, which is related to the VACF via C F (τ ) = −m 2 Cv (τ ).

D. Calculation of the memory kernel
To calculate the memory kernel we use the equation of motion for the autocorrelation dynamics, Eq. ( 7), which does not involve the fluctuating force anymore.In principle, if one knows the autocorrelation function of the observable, which can easily be calculated from simulation data, Eq. ( 7) can be used to obtain the memory kernel.As Mori already showed, 59 formally this equation can be solved via a Laplace transform.Eq. ( 13) allows, in principle, to compute K(ω) from the VACF in the frequency domain, 45,65 however, due to numerical subtleties and possible artefacts of the Fourier treatment, 65 we calculate the kernel in the time domain.In the latter the memory kernel is usually extracted via an iterative or recursion ansatz. 64,67,84,85In this article we use a modified version of the matrix-inversion method introduced in ref. 68 adapted to the stationary case.The details of the adaptations are given in the Appendix C.However, some properties of the memory kernel can be inferred directly.Calculating the time derivative of Eq. ( 7), we see that the memory kernel carries some characteristics of the correlation functions.At short correlation times we find K(τ ) ≈ (β/m) C F (τ ), and from Eq. ( 13) we see that in the low frequency limit, when |ω| ≪ | K(ω)|, we recover the long time Einstein relation,

E. Effective long-time transport coefficients
In this section we briefly recall common methods to calculate the effective long-term equilibrium tracer diffusion, D, and friction, γ, which are connected via the Einstein relation [first equality in Eq. ( 11)].
A very common approach to calculate D from simulation trajectories and single-particle tracking experiments is the use of the mean-squared-displacement (MSD), 86 Equivalently to the MSD method, the diffusion coefficient can be computed from the integral of the velocity autocorrelation function [first equality in Eq. ( 10)], 60 coinciding with the zero-frequency limit of Eq. ( 12), i.e., D VACF = D(0).
As shown in Sec.II C, this is also consistent with calculating the friction constant via the integral over the memory kernel Eq. (11).From the FDT we also know that the same result could be obtained by integrating over the autocorrelation of the fluctuating force However, here it is important to note that this relation holds for the fluctuating forces ζ iν (t) from the projected dynamics and not the total forces F iν (t) = m viν appearing in the experiment or simulation.The two following expressions are not the same: Here, L is the Liouvillian and Q M is the orthogonal Mori projection operator (c.f.Appendix B).The fluctuating forces from the projected dynamics are usually much more difficult to obtain and, hence, the autocorrelation of the forces from the experiment or simulation are often used to estimate friction constants.Only for a tracer particle with infinite mass or under specific constraints, the Green-Kubo relation for FACF yields the correct friction coefficient, [87][88][89][90][91][92][93][94] β ∞ 0 dτ lim m→∞ C F (τ ) = mγ FACF .Finite masses, however, are not covered by the Green-Kubo relation and lead to the plateau problem and the integral vanishes.Hence, the integral needs to be cut off at τ cutoff ≈ γ −1 , which depends on the effective friction itself.In practice, cutting off the integral at its maximum provides results comparable to the other methods.We assume that some confusion regarding the plateau problem has been caused in the past by the fact that ζ iν (t) is called a fluctuating force in the literature on projection operator formalisms, while F iν (t) is a force, which fluctuates.

F. Generalized Stokes-Einstein relation
The Fourier transform of the memory kernel K(ω) plays an important role under mild nonequilibrium conditions, e.g., in the presence of weak sinusoidal driving fext (ω), because it provides information about the elastic and viscous contributions.In the frequency domain, K(ω) directly enters the mobility or, equivalently, the mechanical (viscoelastic) impedance, reading μ−1 = m(−iω + K), and one can already conclude that the viscous part, Re( K), is related to the (frequency-dependent) dissipation. 86,95The imaginary part Im(μ −1 ) contains inertial effects mω, and elastic interactions with the surrounding stored in Im( K).

A. Tracer dynamical correlations
In Fig. 2  The decay time of velocity correlations [panel (a)] decreases with higher interactions strength and system density, and deviates significantly from the simple exponential decay C 0 (τ ) ∝ exp(−γ 0 τ ) of the freely diffusing particle (gray line).The interactions with the dense surrounding thus reduce the mean free path and also the probability of the particles to perform free diffusion, hence reducing the effective diffusivity.
The interaction strength directly affects the magnitude of the deterministic forces acting on the tracer particles, as apparent in the force correlations [panel (b)].High ε tn create a sticky environment, i.e., tracers are likely to adsorb to polymer beads and pick up the network dynamics, evident by the emergent oscillations on the intermediate timescale, τ γ 0 ≈ 0.1, in C v as well as in C F .A more precise analysis of the frequencies and the source of such vibrations is provided later.
In the case of strong adsorption, the apparent decay in the velocity autocorrelation is increasingly determined by network oscillations and their damping times, and the density ϕ appears to have only little effect on the velocity and force correlations.We may conclude that the free diffusion barely occurs for large βε tn .The motion is rather governed by hopping events along the vibrating polymer chains.
Note that the oscillations in C v and C F have the same frequencies, but are shifted by roughly one quarter of the oscillation period [compare the pronounced minima for strong interactions in Fig. 2(a) and (b)].This is expected due to the relation between force and velocity correlations, C F = −m 2 Cv , and can also be rationalized in the harmonic oscillator picture, in which the velocity lags behind the forces by π/2 due to the inertia m.

B. Tracer memory kernel
The memory kernels K presented in Fig. 2(c) are related to C v and C F .With τ → 0, K converges to (β/m) C F .Note that due to the logarithmic scaling of the delay time τ , the delta-peak stemming from the solvent noise ξ(t) in the microscopic model is not visible.Although the shape of K at intermediate times, τ γ 0 ≈ 0.1-1, significantly differs from C F , we find again that the magnitude of the friction directly scales with the interaction strength ε tn , as similarly observed for C F , leading to a common intersection point of all presented kernels at τ γ 0 ≈ 0.1.As already observed for C v and C F , the increased dynamical coupling of the tracers with the network vibrations at higher ε tn is also obvious in the memory, which we discuss in more detail in Secs.III D-III F.

C. Long-time diffusion
In Fig. 3 we compare the effective diffusion coefficients obtained from the different methods, introduced in Sec.II E, for seven different interactions strengths ε tn and two different polymer volume fractions ϕ.The methods based on the MSD, VACF and the kernel represent the very same physics yielding the exact same diffusion coefficients, while the FACF estimate only reproduces the overall trend.
The essentially exponential scaling, D ∝ exp(−Aβε tn ), with scaling A, is known from similar coarse-grained simulations of diffusive tracers in flexible as well as in rather stiff polymer networks. 51,55At low interaction strength, volume exclusion effects D ∝ exp(−Cϕ/(1 − ϕ)) dominate the diffusion process. 48,49In this regime, the higher densities reduce the probability for the particles to perform free diffusion, which is in line with the observations in the VACF [Fig.2(a)] and the kernel [Fig.2(b)].Precisely, for βε tn = 0.1, C v and K are almost purely decaying functions, where C v decays faster with increasing ϕ, and where K shows increased magnitude.
Nonetheless, with increased interactions, we see that the scaling A is smaller for ϕ ≈ 0.32, and, eventually, the diffusion for tracers at very strong attractions (βε tn = 3) is higher in the denser networks.This is usually rationalized by the overlaps of the numerous LJ potentials essentially smoothening the energy landscape and thus creating pathways through the network speeding up the diffusion. 55,103,104

D. Network dynamics
The dynamical coupling between tracers and the network increases with stronger LJ interactions ε tn .Analysing the memory kernel of the polymer beads thus reveals the origin of the observed oscillations in the tracer memory.In Fig. 4 we depict the single particle memory kernel of the (tracer-free) network (with ϕ ≈ 0.32) averaged over all polymer beads (black solid line) in comparison with the tracer memory with moderate (βε tn = 1) and very strong interactions (βε tn = 3) in such a network [already presented in Fig. 2(c)].We also present the memory of the network particles for three auxiliary simulations for which we disabled the non-bonded interactions (bond-only network), turned off all bonds (LJ gas), and removed all crosslinkers (melt of polymer strands with comparable density ϕ ≈ 0.29), respectively.
The first auxiliary simulation represents a network that only interacts via bond and angle potentials (black dashed line in Fig. 4) and exhibits the same underdamped oscillations apparent in the network with the full force field (black solid line).We conclude that the vibrations due to the bonded potential dominate the polymer dynamics.The frequency ω bonds ≈ 25γ 0 is of the order of a single bond oscillation, e.g., a single harmonic potential with spring constant κ b = 100k B T /σ 2 yields ω * harm.= 2κ b /m ≈ 14γ 0 , and the two-bond stretch vibration gives ω * stretch = 4κ b /m = 20γ 0 . 63,105We also observe that the amplitude of the oscillation is more pronounced in the absence of LJ interactions, hence, the non-bonded interactions slightly inhibit the bond vibrations.
In the second auxiliary simulation (black dotted line in Fig. 4) all bonds were removed, resulting in a pure LJ system (with βε nn = 0.5), and one observes the typical signature of a LJ gas with weak attractions. 62,64,69We will demonstrate later that the LJ contributions are described by strongly damped oscillation, with ω LJ /γ 0 ≈ 10-15.In highly diluted systems, the pairwise frequency from a harmonic approximation at the LJ minimum (at r = 2 1/6 σ) provides the lower limit.It scales with ω * LJ ≈ σ −1 25ε/(2m), which yields ω * LJ (βε = 0.5) = 2.5γ 0 in our case.The actual frequency is higher, since, on the one hand, the LJ minimum is considerably steeper than the second order approximation, and, on the other hand, due to collective effects, such as the superposition of multiple potential wells.
Considering finally a polymer melt (11664 strands of ten monomers each and no crosslinks), which is our third auxiliary simulation (depicted as orange solid line in Fig. 4), the memory kernel is almost identical to the one for network (with bonded and LJ interactions) at small timescales, attributed to the similarities in the local environment and interactions.
At large delay times we find long subdiffusive tails in the memory kernels for the 'bonded networks' and the polymer strands (see inset in Fig. 4), while the LJ gas converges quickly to linear diffusion.In fact, we find strongest subdiffusion for the network with the full force field, for which the tail of the memory kernel scales with roughly τ −1/3 .Hence, the very long subdiffusion is a network effect, which is significantly magnified by LJ interactions. 35,37This conclusion is further substantiated by the observed long-time memory of the polymer strands.While the magnitude of the subdiffusive tail is increased by the LJ interactions (compare with the bond-only network), it also decays considerable faster then the tails of both 'bonded networks'.
Inspecting the tracer memory, we find a clear LJ signature for the moderate interaction strength (βε tn = 1), and only at higher interactions (e.g.βε tn = 3) the tracers pick up the bond vibrations and also exhibit a subdiffusive regime.The subdiffusion of the tracers, however, decays notably faster than the ones of the 'bonded networks', and the polymer strands.Details on the dynamical coupling become evident in the frequency domain provided below.the tracer memories from the original simulations for comparison (as in Fig. 2).The inset depicts four memory kernels with clear subdiffusive regimes at large delay times.The long-time tail of the network memory kernel (with full force field) scales roughly with τ −1/3 .See main text for more details and interpretation.
which is related to the dynamic relaxation modulus G(ω) of the whole network.

E. Frequency domain
Analyzing the memory in the Fourier space allows us to identify the spectral density of the correlated fluctuations ζ, interpret the frequency-dependent friction, and extract the dominant vibration modes contributing to the memory kernel.In Fig. 5(a) we depict K(ω) for different ε tn and for tracers in the dense network only.The low-frequency limit, Re[ K(0)] = γ kernel , is identical with the integral of the memory kernel [Eq.( 11)], and confirms the above presented longtime limit.At the low, yet non-zero frequencies, the slow convergence of Re[ K(ω)] to γ kernel signifies the existence of very long subdiffusive cross-over regimes, corresponding to the long tails in K(τ ) .
At high frequencies Re[ K(ω)] shows low-pass characteristics, but the friction is limited by uncorrelated solvent friction γ 0 (indicated by gray horizontal line).Note that γ 0 is not perfectly reached since the polymer beads are also subject to uncorrelated solvent fluctuations, which are passed on to the tracers with increasing coupling strengths βε tn .
At intermediate frequencies (ω/γ 0 ≈ 10-50) we find the superposition of several damped oscillations, whose magnitudes and resonating behaviors increase with ε tn .The maximization of the imaginary part indicates the presence of elastic interactions and corresponds to the emergent oscillations observed in the correlation functions and in the memory kernel (Fig. 2).The plateau in the imaginary part at lower frequencies originates from slower, elastic network deformations [cf.network memory in frequency domain in Fig. 5(b)], yet only plays a minor role due to the overall large friction in the low frequency regime.The real part, in contrast, is dominant and approaches slowly to the long-time limit (ω → 0), since the adsorbed tracers adopt the network subdiffusion for relatively long timescales ω ≈ 0.1γ 0 , and only detach from the polymer at very long times, which also coincides with the vanishing Im( K).
Inspecting the network kernel and the auxiliary simulations, depicted in Fig. 5(b), allows us to deconvolute the effects again.We see a strong resonating behavior stemming from the bond vibrations, a rather damped signature from the LJ interactions, and low-frequency scaling typical for networks.The latter can, to some extend, be compared to viscoelastic experiments of gels.For instance, Winter 74 evaluated the frequency scalings of the storage and loss modulus during crosslinking polymerization in order to detect the gelation point (fully connected polymer network) and proposed the general scaling Re(G) ∝ ω n and Im(G) ∝ ω n with 0 < n < 1, which relates the storage and loss, reading tan(nπ/2) = Im(G)/Re(G).Given the GSER is valid in the low-frequency domain, one may assume n ≈ (2/π) arctan[−Re( K)/Im( K)], resulting in n ≈ 2/3 for ω ≲ 0.1 in our simulations, well in line with typical hydrogels (n ≈ 0.5 -0.8). 38is further provides the correct scaling (n − 1) for the memory kernel, i.e., Re( K) ∝ ω −1/3 and Im( K) ∝ ω −1/3 , as depicted in Fig. 5(b), reflecting the τ −1/3 -scaling of the subdiffusive tails (cf.inset in Fig. 4).We hence may conclude that the low-frequency scaling is dominated by very slow network dynamics.In fact, it decouples completely from the fast bond vibrations, evident in the comparison of K for networks with varying bond stiffness (see.Fig. 7 in the Appendix D).

F. Damped harmonics approximation
We now approximate the memory kernel with a series of damped oscillations, reading summarizing the leading modes and their dependence on tracer-network coupling ε tn .In fact, the magnitude for constant ϕ scales with ε tn , and we extend Eq. ( 21) with where α k is the mode-dependent scaling.][108] With this extension we aim to fit as many memory kernels simultaneously as possible with the smallest number of modes.For demonstration, we use the dense system (ϕ ≈ 0.32) and only attractive interactions (βε tn ∈ {1.0, 1.5, 2.0, 2.5, 3.0}), i.e., the regime in which tracer adsorption to the polymer plays a significant role.The fitting focuses on the short and intermediate timescales and only considers the kernel data up to τ = 10γ −1 0 .The fitting results with N k = 4 modes are summarized in Table I, ordered by magnitude λ k , and confirm our observations from the previous sections.
The dominant contribution (k = 1) at small delay times stem from oscillations close to critical damping with a frequency of ω 1 ≈ 10γ 0 .We have already identified this mode as a clear LJ signature, and the almost linear scaling (α 1 ≈ 1) with βε tn confirms this interpretation.The second strongest mode with ω 2 ≈ 25γ 0 is most obvious at intermediate delay times, τ γ 0 ≈ 0.1-1.0, and quantifies the coupling of the tracer dynamics to the polymer bond vibrations with a scaling of α 2 ≈ 1.3.The third mode has rather small impact and cannot be directly attributed to a single interaction type.It is rather a result of the assumption that only the magnitude scales with the LJ interactions.In fact, this mode shares the same decay time as the dominating ones, and due to ω 3 ≈ 44γ 0 it essentially performs a small shift of the extrema to smaller τ with TABLE I. Memory kernel fit parameters for ϕ ≈ 0.32 and attractive tracer-network interactions εtn ∈ {1.0, 1.5, 2.0, 2.5, 3.0} with N k = 4 modes.The kernels with different εtn were fit simultaneously by Eq. ( 21) and introducing the scaling [cf.Eq. ( 22)] with parameter α k .The delta peak stemming from the Langevin thermostat [cf.Eq. (1)] was not included for fitting and needs to be added manually to the memory kernel to account for the solventrelated friction.magnitude decay frequency scaling origin  21) and the εtn-scaling Eq. ( 22), which were used to inter-and extrapolate kernels for the GLE simulations for the test data.
increasing ε tn , indicating that also ω k and τ k for the first two modes slightly scale with ε tn .Mode four has the smallest magnitude and only plays a minor role for small τ .However, it has a significantly longer decay of τ 4 γ 0 ≈ 2.5 and is decisive for the long subdiffusive tails observed for strong LJ attractions.We verify our fitting approach by comparing the timedependent diffusion coefficients, D(t) = t 0 C v (t ′ )dt ′ , of the microscopic simulations with the GLE simulations presented in Fig. 6.Details on the GLE simulations are provided in the Appendix B. In panel (a) we compare the two simulation techniques for the fit data and only observe significant deviations for larger times and very strong coupling.For the results depicted in panel (b) we have inter-and extrapolated the GLE model and compare the results to those of microscopic simulations that were not used for fitting.We find a very good shorttime agreement up to t ≈ γ −1 0 .At relatively long timescales, up to t ≈ γ −1 0 , the interpolated kernels perform well, and only at very large times we see significant deviations, particularly for extrapolated kernels at the largest and smallest tested ε tn values.

IV. CONCLUDING REMARKS
As motivated in the introduction, tracer diffusion in polymer networks is not only important for applications but provides valuable insight into the physics of fluctuations and transport in soft, viscoelastic media.Long-time relaxation modes, describable by memory kernels, are key to the rationalization of the viscoelasticity of the polymer matrix but a systematic study of the kernel behavior depending on the tracermatrix coupling had been missing, yet.Here, we presented a well-controlled tracer-network simulation model, in which memory kernels could be calculated with unprecedented statistical quality.This allowed us to identify various, frequencydependent contributions of the kernel, and to discuss the origins of the emerging physics separately, depending on the tracer-polymer coupling.
In particular, with increasing attraction strength between tracer and polymer, the tracer is more likely to adhere to the polymer matrix and pick up its dynamics; in here, we identified the signature of the LJ liquid, emergent bond vibrations, and the slow subdiffusion of the whole network on larger scales.The magnitude of these effects entering the kernel scales almost linear with the LJ interaction strength.With that we demonstrated how the single solute kernel analysis may be employed to retrieve dynamic and viscoelastic information about the solute transport as well as about the surrounding polymeric matrix.
The analysis allowed us further to construct a GLE in which the polymer network enters effectively only through the essential modes (using an auxiliary variable method) in the memory kernel and the corresponding fluctuations.In this reduced description, the polymer network serves now as an additional bath with controlled coupling to the tracer system, enabling the efficient inter-and extrapolation of the friction effects.This will be useful to predict the outcome of other couplings and modified system setups in efficient GLE-based simulations in future.
The recognized modes may also be useful to construct or optimize empirical kernels to describe and interpret specific experiments, e.g., for the mechanical tracer response in in the strong tracer-polymer coupling limit, 44,63,73 or for the nonequilibrium response in biological networks, as demonstrated recently. 1094][115] Importantly, adding external forces to the underlying microscopic model may result in an additional terms describing nonequilibrium response forces in the GLE, which can only be neglected in the zero-force limit. 116

APPENDIX Appendix A: Microscopic simulation details
The creation of the regular polymer network and the simulations of the Markovian Langevin model, Eq. ( 1), were performed with the LAMMPS 117 software package, similar to our previous work. 55The network consist of 122472 polymer beads and has a crosslink ratio of 5% yielding a chain length of ten monomers between two crosslinker beads.We fixed the bond and angle potentials (see Table II) and tuned the polymer volume fraction ϕ with the LJ network-network interactions ε nn by equilibrating the network in N pT ensemble with an isotropic Berendsen barostat (τ p = 0.1τ 0 , p target = 0) and the Langevin integration scheme.The equilibration was carried out until convergence of the cubic simulation box volume was reached, yielding ϕ ≈ 0.16 for βε nn = 0.35 and ϕ ≈ 0.32  for βε nn = 0.5.This corresponds to box lengths of L = 73.3σand L = 59.5σ for the moderate and the high polymer density, respectively.
After the network equilibration, the simulation box volume was fixed, and 1956 (for ϕ ≈ 0.16) and 1041 tracer particles (for ϕ ≈ 0.32) were randomly added to the system, yielding a tracer concentration of c ≈ 0.005σ −3 .The tracer positions were subsequently minimized and the whole system was again equilibrated in N V T ensemble (Langevin integration) for τ EQ = 100τ 0 with a timestep of dt = 0.001τ 0 .The production run was carried out with identical timestep for τ PR = 10 3 τ 0 and the tracer trajectories were output every timestep.
In order to evaluate the dynamics of the network beads we ran additional simulations without tracer solutes.Instead of storing the single particle trajectories, we made use of LAMMPS' on-the-fly calculation of the velocity autocorrelation function (VACF), providing the VACF averaged over all network beads.We performed 40 consecutive production runs for τ PR = 10 3 τ 0 with dt = 0.005τ 0 , yielding 40 individual VACFs, where the initial timestep of each run was used as reference for the delay.The final result was obtained by averaging over these individual VACFs.This procedures was employed for the network with the full force field as well as for the auxiliary simulations with partial force fields [as presented in Fig. 4 and Fig. 5(b)].From the VACFs the memory kernel was computed as described in Sec.II D.
The ρ 0 (Γ) in Eq. (B5) denotes the stationary phase-space distribution sampled from the simulations or experiments.From these definitions it is easy to see that the Liouvillian is an anti-self-adjoint operator, i.e., a = 0, and that the fluctuationdissipation theorem ⟨ζ iν (t)ζ iν (t ′ )⟩ = m 2 K(t − t ′ ) v 2 iν (B7) holds true.Note that the entire formalism can also be used for simulations that involve a Langevin thermostat as shown in ref. 118.
. Each mode has delta-correlated noise ⟨ξ k (t)ξ k (t ′ )⟩ = k B T 4τ −2 k δ(t − t ′ ), satisfying the FDT.To account for the solvent-related fluctuations, we added −mγ 0 v + ξ 0 (t) with ⟨ξ 0 (t)ξ 0 (t ′ )⟩ = 2k B T mγ 0 δ(t − t ′ ) to Eq. (E2), where γ 0 actually also scales with ε tn [see high frequency limit in Fig. 5(b)], but we approximated it by γ 0 = γ 0 for simplicity without noticeable impact.In fact, even dropping the solvent-related terms completely would not significantly change the simulation results (not presented) since the tracer dynamics are dominated by the polymer bath in the attractive regime.

FIG. 1 .
FIG. 1. (a): Simulation snapshot of the microscopic (Markovian, implicit-solvent) model of diffusive tracer particles (blue) in a (periodic) polymeric bead-spring network (red: chain segments, yellow: crosslinkers).While the tracer particles are displayed with correct size (σ), the chain monomer and crosslinker diameters are reduced for illustration.(b): Schematic representation of the microscopic model displaying classical Langevin solvent contributions and pairwise deterministic forces acting on the tracer particles.(c) Coarse-grained generalized Langevin equation (GLE) model that incorporates solvent and polymer bath effects through the memory kernel K and correlated fluctuations ζ.
we depict the velocity autocorrelation C v (τ ) [panel (a)], the force autocorrelation C F (τ ) [panel (b)], and the memory kernel K(τ ) [panel (c)] as a function of the delay time τ for different tracer-network interaction strengths ε tn [color-coded, see legend in panel (a)], and two different polymer densities [see line pattern legend in panel (b)].

FIG. 3 .
FIG. 3. Long-time diffusion coefficients for different methods (cf.Sec.II E) represented by different symbols and colors (see legend in panel (a)).The panels depict the results for different densities, (a) ϕ ≈ 0.16 and (b) ϕ = 0.32.While the diffusion calculated from the MSD, VACF, and the kernel perfectly match, the FACF method only captures the rough scaling.

FIG. 4 .
FIG. 4.Comparison of the tracer memory kernel with the network memory kernel.Black solid line corresponds to the kernel averaged of all polymer beads with full force field (networknetwork LJ interactions βε = 0.5 ⇒ ϕ ≈ 0.32 and bonded potentials) as used in the simulations (as in Fig.2).The black dashed line depicts K for the same network at same density but without the LJ interactions, hence extracting pure bond contribution to the network memory.The black dotted line depicts the kernel of the polymer beads after the removal of all bonded potentials, hence representing a LJ gas with βεnn = 0.5.The solid orange line reflects the memory kernel of monomers in a polymer melt with strand length of ten monomers, no crosslinks, and ϕ ≈ 0.29.The purple and the green solid lines ('tracer') are

FIG. 5 .
FIG. 5. (a)Real and imaginary parts of the frequency-dependent tracer friction K(ω) in the dense network (ϕ ≈ 0.32) for different interaction strengths βεtn, color-coded as in the Fig.2[from dark blue (βε = 0.1) to purple (βε = 3.0)].The real part is decisive for viscous loss, where the low and high ω limits correspond to effective long-time friction K(0) = γkernel and the uncorrelated solvent-related friction K(ω → ∞) ≈ γ0, respectively.The imaginary part reflects elastic collisions with the polymer and rises with increasing coupling to bond vibrations (ωbonds ≈ 25γ0), but also slower network modes in the low-frequency domain (ω/γ0 ≈ 0.1-1).For ω → 0, we observe a vanishing imaginary part, indicating that the adsorbed particles eventually detach from the polymer for very long times.(b) Real (black) and imaginary (red) parts of network memory in the frequency domain for the network (ϕ ≈ 0.32) with the full force field and the two auxiliary simulations, the bond-only network and the LJ gas.The low-frequency scaling Re( K) ∝ Im( K) ∝ ω n−1 with n − 1 ≈ −1/3 can be computed via n ≈ (2/π) arctan[−Re( K)/Im( K)], which is related to the dynamic relaxation modulus G(ω) of the whole network.

TABLE II .
Parameter summary for the microscopic model.