Concurrent two-scale model for the viscoelastic behavior of elastomers filled with hard nanoparticles

A dynamic two-scale model is developed for describing the mechanical behavior of elastomers filled with hard nanoparticles. Using nonequilibrium thermodynamics, a closed system of evolution equations is derived, coupling continuum mechanics with a fine-scale description on the level of filler particles. So doing, a constitutive stress–strain relation emerges that is applicable to transient situations. In addition to the number density of filler particles, the particle arrangement is captured by the distribution of the difference vector between two representative interacting particles, which makes this model efficient in comparison with many-particle models. The two-particle model presented here is analyzed numerically in oscillatory deformation, for two purposes. First, the nonlinearity of the model is studied in detail, in terms of the Payne effect, that compares favorably with the literature. And second, the two-particle model is compared with a corresponding many-particle model in the literature.

than on the effect of clusters and aggregates. Therefore, silica-filled elastomer nanocomposites (in the sequel simply called "nanocomposite" to emphasize the role of the microstructure) serve as a prototype example for the approach presented here.
For the description of the mechanical behavior of nanocomposites, adequate modeling is cumbersome. This is due to the intricate phenomena on the level of the nanosized filler particles that lead to the reinforcement effect on the macroscopic engineering scale (see, e.g., [1]). A feature of particular interest in this paper is the following. While the elastomeric matrix can to a good approximation be described by the theory of elasticity [2][3][4], filling the elastomer with hard silica particles adds significant dissipative contributions to the material behavior. Although filled elastomers have been studied for many years, the physical origin of these significant dissipative effects was not understood until recently. Key are the dynamical properties in the vicinity of interfaces, as they have been probed by studies such as NMR in [5]. More recently, model systems have been studied with well-dispersed spherical silica particles with a diameter of about 50 nm and typical distance between the particles of about 20 nm [6,7]. By studying the mechanical properties of the samples as functions of frequency and temperature, it was shown [6,7] that the dynamical properties are modified as compared to those of the pure elastomer. The interpretation proposed is that the presence of interfaces, if the interaction between the polymer and the substrate is sufficiently strong, results in an immobilization and thus in an increase in the glass transition temperature T g of the polymer matrix close to the interface. This increase of T g can be described equivalently by a glassy layer around the filler, given by the region in which T g is larger than the actual temperature T . If the glassy layers of different particles do not overlap, the effect on reinforcement (defined as the ratio between the modulus of the filled sample and that of the unfilled sample) is essentially geometrical: The increase in the elastic modulus can be quantitatively explained by an increase in the effective filler volume fraction due to the presence of the glassy layer. Upon lowering the temperature, the thickness of the glassy layer increases to the point that the glassy layers overlap, leading to so-called glassy bridges. This happens at temperatures close to T g,b + (10-15)K in the considered samples (where T g,b is the glass transition temperature of the elastomer in the bulk), resulting in very high reinforcement up to about 100, depending on the filler volume fraction. In the case that so-called glassy bridges are formed between adjacent particles, the phenomenology of the resulting particle interaction can be described similar to that of bulk glassy polymers [8][9][10]. Correspondingly, the strength of a glassy bridge depends on both temperature and applied load. In other words, the behavior of bulk polymer glasses is helpful to rationalize typical effects in nanocomposites. Particularly, the yielding of glassy bridges under applied load macroscopically leads to decreasing the elastic modulus G , i.e., the Payne effect [11]. Furthermore, inter-particle glassy material may also age physically after cessation of large-amplitude oscillatory deformation, leading to the Mullins effect [12,13].
Many efforts have been made to investigate the properties and mechanical behavior of silica-filled nanocomposites. Among the numerical approaches, both molecular dynamics (MD) and Monte Carlo (MC) simulation techniques have been performed on the smallest scales. For example in [14,15], MD simulations have been performed of thin films of noncross-linked polymer melts confined between substrates. In both studies, the simulations revealed a decrease in mobility close to the substrate and an increase in the thickness-averaged glass transition temperature upon decreasing the substrate spacing; these results corroborate the existence of glassy layers close to the substrate. In [16], the properties of a polystyrene matrix filled with nanoparticles were investigated by coarse-grained MC sampling, and a substantial segmental ordering close to the particle surface was found. On coarser scales, the very broad study of filled (either carbon black or silica particles) elastomers from the microscopic perspective was made in [17][18][19][20], where the authors developed a numerical scheme akin to dissipative particle dynamics (DPD) technique. The dynamics of the filler particles has been studied, taking into account soft spring forces, representative of the rubbery matrix, and temporary hard spring forces, representative of the transient network of glassy bridges. The essential idea in their approach is that each glassy bridge has got a finite lifetime that depends on the local glass transition temperature, which in turn is affected by the particle distance and the local load (i.e., the lifetime is changing in the course of time). This concept proved to be key for appropriately simulating [18][19][20] various nonlinear effects of filled elastomers, e.g., the Payne effect [11] and the Mullins effect [12,13]. On the other hand, fully macroscopic numerical approaches, based on providing a multi-level finite element method (MLFEM) simulation for heterogeneous polymeric systems, have been employed to study specifically polymeric matrices with either voids or filled with rubber inclusions [21][22][23][24]. There, a homogenization method has been proposed to account for the large deformation of viscoelastic media, considering both microscopic and macroscopic levels, and implemented to MLFEM afterward.
An analytical approach to describe filled systems is given in [25][26][27], where concepts from fractal geometry were applied. In [25], colloidal carbon black or silica aggregates were considered as fractal objects, since their structures have scale invariance. These authors have established power-law relationships between stress and elongation of filled systems, with an exponent that depends on the fractal dimension. A broad and detailed overview on the filled system is given in [26] and [27], where the viscoelastic properties of filled elastomers were considered. A power law between elastic modulus and filler volume fraction has been established [26], and a constitutive micro-mechanical model of stress softening of filled elastomers was developed, up to large imposed strain [27].
In this article, our main goal is to develop a two-scale model that describes the mechanics of silica-filled elastomers. Particularly, we strive to explicitly account for the mutual coupling of the macroscopic engineering scale with the underlying filler-particle scale. Since transient effects are hallmarks of filled elastomer mechanics, our approach should establish the coupling of both levels also in transient situations. Finally, the approach should be applicable to studying macroscopically inhomogeneous deformation, in order to enable studies of wear and tear experiments. This last criterion in combination with the requirement for dynamics asks for alternatives to the approaches mentioned in the previous paragraphs. The work of Long et al. [18][19][20] discussed above shows the significance of the local glass transition temperature and dynamics that depend on the local structure and load, and this complexity helps to understand why until recently no mechanical (viscoplastic) model was able to satisfactorily explain and reproduce mechanical properties of these systems, and in particular their dissipative properties. However, since their model considers a many-particle system to obtain the effective constitutive behavior (at each macroscopic material point), it is computationally rather demanding. Therefore, the aim of this paper is to propose a drastically simplified description of the model developed in [18][19][20] and to combine that with macroscopic elastoviscoplasticity (e.g., see [4,28]) in a thermodynamically consistent manner.
Before we start, let us give a word about notation used in this manuscript. Throughout the entire manuscript, Latin indices i, j, k, . . . denote components of a set of dynamical variables, while Greek indices α, β, γ , . . . are used for the Cartesian components of vectors and tensors. Einstein's summation convention is used for indices that occur twice. Furthermore, subscripts α and (α, β) imply contraction with any vector A α and any matrix A αβ multiplied from the left, respectively, while subscripts γ and (γ , ε) imply contraction with the vector A γ and matrix A γ ε multiplied from the right, respectively. Finally, it is pointed out that boundary terms are neglected throughout the entire paper. This simplification is justified because this paper is concerned with modeling bulk material properties, in terms of local evolution equations. It is therefore implicitly assumed that the boundary does not affect the dynamics in the inside over a distance (apart from the fact that any bulk dynamics must obviously be supplemented with proper boundary conditions when a numerical solution is sought after). For completeness, it is mentioned that systems that interact with the environment have been studied earlier, e.g., [29,30].
The manuscript is organized as follows: In Sect. 2, the GENERIC formalism of nonequilibrium thermodynamics is introduced, and in Sect. 3 this technique is applied to filled systems in order to formulate evolution equations and constitutive relations. The general two-scale model is summarized in Sect. 3.5. In Sect. 4, specific exemplary choices are made to concretize the model and translate it into stochastic differential equations that are amendable to numerical simulations. In Sect. 5, the model is solved numerically for a concrete examplenanocomposite under oscillatory deformation. Specifically, the two-particle model presented here is compared with the corresponding many-particle model developed in [18][19][20], and the Payne effect is examined. The paper closes with a conclusion, Sect. 6.

Methods: nonequilibrium thermodynamics
As a guideline for developing the two-scale model, nonequilibrium thermodynamics is used. In the end of nineties, the framework called General Equation for the Non-Equilibrium Reversible-Irreversible Coupling formalism (GENERIC) has been developed, which is used to describe closed nonequilibrium systems [31][32][33]. This method has two important advantages over other nonequilibrium thermodynamics techniques. First, it is useful for setting up multi-scale models. And second, GENERIC provides a concrete procedure for connecting different levels of description, i.e., for coarse graining.
The first step in setting up the GENERIC for a concrete system is to specify a set of dynamic variables, denoted by x. The time evolution equation is then written as The derivatives of the energy E and entropy S, that both are functionals of the dynamical variables x, drive the reversibleẋ| rev and irreversibleẋ| irr dynamics, respectively, through the Poisson operator L and the friction matrix M, which in general also depend on x. The symbol means not only summations over discrete indices but may also include integration over continuous variables; this issue will be discussed in the next Section when applying the GENERIC. The operator δ/δx represents functional rather than partial derivatives. The Poisson operator is equivalent to a Poisson bracket { , }, and the friction matrix is associated with a dissipative bracket [ , ], for two arbitrary functionals A and B. The Poisson operator and the friction matrix must satisfy a number of conditions. First, there are the degeneracy conditions. Namely, the functional derivative of the entropy lies in the null space of the Poisson operator L, and the functional derivative of the energy lies in the null space of the dissipative matrix M, As a second set of conditions, L must be antisymmetric, whereas M is positive semi-definite and Onsager-Casimir symmetric 1 where the sign * denotes the adjoint operator. The final condition is that the Poisson operator must satisfy the Jacobi identity (see [33, pp. 14-16]), for all functionals A, B and C of the variable x. It is pointed out that, rigorously speaking, the conditions (3)(4)(5) on the operators L and M must be checked in terms of the corresponding brackets, {, } and [, ], respectively. In the model formulation in Sect. 3, these brackets are integrals, and the verification of the conditions (3)(4)(5) involves integrations by parts. As explained in Sect. 1, due to our interest in bulk material behavior, one may neglect the boundary terms that arise as a result of the integrations by parts. Using the above conditions, one can show that the reversible dynamicsẋ| rev possesses the Hamiltonian structure, while the irreversible dynamicsẋ| irr is a generalization of the Ginzburg-Landau equation. Furthermore, the above stringent conditions (3)(4)(5) imply the desired energy conservation and nonnegative entropy production,Ė = 0,Ṡ ≥ 0.
For more details about the GENERIC framework, the reader is referred to [33]. The time evolution Eq. (1) for x contains a reversible contribution,ẋ| rev , and an irreversible part,ẋ| rev . Usually, reversibility means invariance upon time reversal, t → −t. However, in this study, it is rather a split between contributions related to affine deformation, that will be captured by L (δ E/δx), and others, that will be captured by M (δS/δx). We will come back to the difference between the time-reversibility-split and the affine-nonaffine-split in Sect. 3.4.2.

Choice of variables
The envisioned two-scale model should carry the following features. First, the model should account for macroscopic continuum mechanics in order to allow for the description of macroscopically inhomogeneous deformations. Second, the effects of the nanofillers are to be included in such a way that the understanding of the nanofiller-dynamics developed in [18][19][20] (e.g., the relevance of glassy bridges) can be incorporated, i.e., new dynamic "microscale" variables must be included. And third, the model should properly couple the two levels, not only in stationary but also in transient situations.
While the mechanical behavior of the nanocomposite as a whole is viscoelastic, as evidenced by experimental data for the storage and loss moduli, we choose in this paper to distinguish between the matrix and the filler-particle contributions. In terms of variables, separate descriptors are used to quantify the state of the elastomer matrix and that of the network of particles connected by glassy bridges, respectively. The elastomer matrix is considered as purely elastic with a low modulus, while the stiff network has viscoelastic, i.e., ratedependent, behavior with a high elastic modulus and relaxation processes due to the transient nature of the glassy bridges.
It has been shown [4,28] that macroscopic finite-strain nonisothermal elasticity and viscoplasticity can in principle be formulated within the GENERIC framework. Here, we make use of these earlier efforts to design a properly adjusted model for describing silica-filled nanocomposites. To describe the macroscopic finite-deformation elasticity of the elastomer matrix in an Eulerian setting, we choose the momentum density u(r), the absolute temperature T (r) and the total deformation gradient F(r). The vector r is a macroscopic position vector. The velocity field can be obtained from the momentum density as v = u/ρ, where mass density, ρ, can be derived from the deformation gradient, namely where ρ 0 denotes the mass density in the undeformed state. In addition to the variables described above, another variable must be incorporated on the level of the filler particles to account for the effect of glassy bridges adequately. For that, we take inspiration from the many-particle model developed in [18][19][20], in which several force contributions are included between filler particles, e.g., a force contribution to the glassy bridges (see left box in Fig. 1). However, for reasons discussed in the Introduction, we are here looking for a drastic simplification of that model. Particularly, rather than Fig. 1 Left System of many particles (green) with glassy layers around them (gray), and a representative particle pair highlighted in yellow. Right In the representative particle-pair model, only two particles are considered explicitly (green with yellow layer), while the other surrounding particles are taken into account only implicitly (see text for details) (color figure online)

Fig. 2
Representative particle pair under applied shear deformation. Imposed deformation moves the particle configuration from an initial state (light green and gray) to another state (green and dark gray) (color figure online) studying a collection of many particles, we consider a representative particle pair only (see right box in Fig. 1), since the mechanics of the glassy-bridge network roots in these nearest-neighbor interactions.
The representative particle pair is characterized as follows. The current separation vector between the two particles shall be denoted by R. For the further arguments in this paragraph, we concentrate on the glassy-bridge network, neglecting other effects of the matrix material. If deformation is imposed on the entire nanocomposite, the glassy-bridge network will deform, and with it also the representative pair. If the imposed deformation is so rapid that the glassy bridges cannot relax (i.e., yield), then we expect that the network will return to its original state upon rapid unloading, due to the interactions between all particles. However, considering only a single representative pair, this restoring mechanism will not occur automatically. To nevertheless allow for the restoring toward the unloaded state to occur, we include another vector Q, which describes the particle separation vector in the unloaded state (Fig. 2). Therefore, upon rapid loading and unloading, the 'reference' Q will stay constant, while the actual state R departs from Q upon loading and is pulled toward Q upon unloading. If, however, the loading happens at a rate on which the glassy bridges can yield, unloading will not lead to the original state. We will capture this by setting up Q-dynamics that is of pure relaxational form toward R. So, for example, in a step-strain experiment, R immediately attains the new actual conformation, while Q slowly relaxes toward R by re-shaping the glassy network. This implies that (i) the force experienced in the deformed state will decrease with time (since Q relaxes toward R), and (ii) upon sudden unloading during this relaxation process the actual R does not relax to its original state, but rather to the new, relaxed state Q.
The use of the particle-based vectors R and Q is akin to, in the continuum setting, the multiplicative decomposition of the total deformation gradient F into elastic F e and plastic F p parts, respectively, F = F e · F p [34]. To illustrate this analogy, let us denote the particle separation vector in the undeformed state by Q 0 . We then can say R = F · Q 0 and Q = F p · Q 0 , which leads to the following intuitive relations. 2 First, the dynamics of the unloaded state Q is related to the dynamics of the plastic deformation gradient F p , which is of purely irreversible origin. And second, the difference R − Q is related to the elastic deformation gradient, namely R − Q = (F e − 1) · F p · Q 0 . Since, in the absence of strain hardening and related effects, F e − 1 is known to be the origin of elastic stresses, so is R − Q the origin of the structural restoring force. This implies that the force for the particle separation vector depends on R and Q only through the difference R − Q. It must be noted that the strength of that particle-level restoring force depends on the connectivity in the entire glassy-bridge network. However, since we consider only a representative particle pair, such connectivity information is not available explicitly, but can enter in the representative pair idea only through an appropriate choice of the force strength. This aspect of the model is discussed in detail in Sect. 4.1 and "Appendix D".
For practical purposes, it is easier to consider the statistical distribution function p of the vectors R and Q as dynamic variable rather than the vectors themselves, p(r, R, Q), which depends on the macroscopic position r. To be more precise, p is the product of the (normalized to unity) probability distribution function of R and Q, denoted byp, and the number density of particle pairs, n, p(r, R, Q) = n(r)p(r, R, Q) , where n, using the normalization ofp, is given by The integral is performed over the set Q 2 , that is a Cartesian product of two subsets of Euclidean space, i.e., It should be noted thatp(. . . , R, . . .) is closely related to the pair correlation function and hence in principle can be measured experimentally, since it is the Fourier transform of the structure factor. However, in contrast to the pair correlation function, we are here rather looking at the first peak only, since we are interested in representative pairs with glassy bridges, i.e., nearest-neighbor arrangements. While the use of the distribution functionp seems in close analogy to usual practice in polymer kinetic theory [35,36], it should be emphasized that the concept of a variable to quantify the unloaded state, Q, is quite different. In polymer kinetic theory, an entire polymer chain can be represented by a single end-to-end vector. This end-to-end vector of the polymer chain, R, is not only affected by the imposed deformation, but it also enters into an entropic spring force, f S (R). For purely statistical reasons, i.e., maximization of the configurational entropy, that force drives the polymer toward the spherically symmetric coil conformation, which is achieved for f S (R) = 0 at R = 0. Therefore, in these systems there is no need for an additional variable to denote the unloaded state. The situation is quite different for the case of nanocomposites discussed in this paper. The current position separation vector R of the representative pair reacts immediately to imposed deformation; however, this vector does not carry enough information for the elastic restoring to take place upon unloading. For a representative pair in the nanocomposite, the force balance is achieved not by R going to zero, but rather by R approaching the previous unloaded state, Q.
In summary, according to GENERIC, in a Eulerian setting, the set of independent state variables is where r ∈ , ⊂ R 3 , and R, Q ∈ Q, Q ⊂ R 3 . For later convenience, we introduce the variable ξ , called internal configuration.

Generating functionals: energy and entropy
For the specification of the total energy E and entropy S, we refrain from making explicit choices since this is (i) not required at this point, and (ii) would disguise the structure of the equations. Let us consider the system occupy an area B r ∈ , ξ ∈ Q 2 . It is reasonable to only make the following ansatz for the energy and entropy, with v = u/ρ for the velocity, ρ ,F = −ρ F − , and F − the transpose of F −1 . For the partial derivatives, we use the notation f ,F ≡ ∂ f /∂F. The functionals U and S will be concretized in Sect. 4.

Reversible dynamics: Poisson operator L
The reversible part of the x-evolution equation is presented by the first term on the right-hand side of general evolution Eq. (1). For specification of the Poisson operator, we note that in continuum mechanics the term reversible deformation is synonymous with affine deformation, which is the basics of finite-deformation elasticity theory [3,4,37]. In an Eulerian setting, each evolution equation must contain convective terms, i.e., a contribution proportional to the velocity field v. Since the u-derivative of energy, is the macroscopic velocity field and u is the first variable in the set (10), one may use the first column of the Poisson operator to represent exactly those convective terms. In addition to the macroscopic convective terms, the momentum balance equation must contain the divergence of the stress tensor, and therefore, the first row (i.e., the transpose of the first column) of the Poisson operator should reproduce that contribution. On the particle level, reversible dynamics for the structural variable, p(r, ξ ), occurs due to the dynamics of R, which we assume 3 to be deformed affinely with a macroscopic flow field v. To study this in detail, let us consider two particles, given by their position vectors R 1 , R 2 , respectively. If both particle positions are following the flow field v, and if we assume that the particles are rather close to each other compared to the characteristic length scale for v-variations, one finds for the evolution of the connector vector This result is completely analogous to the dumbbell model in complex fluids, and it is already formulated within the GENERIC framework (see [33, pp. 134-136]). We note that the variable Q is irrelevant at this point, since it has relaxational dynamics only. To concretize the above arguments in terms of specific forms for the elements of L, one can make use of earlier results for finite-strain elasticity [4], together with the dumbbell model in complex fluids [33], formulated both within the GENERIC framework. So doing, the Poisson operator assumes the form The elements L (uu) (r), L (Fu) (r), and L ( pu) (r, ξ ) are a direct consequence of the (known) kinematics of affine deformation [4,33]. The elements L (u F) (r) and L (up) (r, ξ ) follow due to the anti-symmetry property (4a) of the Poisson operator. L (uT ) (r) can then be obtained from the degeneracy condition (3a), and again using the antisymmetry property (4a) leads to L (T u) (r). All elements being determined (see 6 for the explicit expressions), the Jacobi identity (5) can be verified by a lengthy, but straightforward calculation. For further details on the calculation of the elements of L, the reader is referred to "Appendix A". The construction of the Poisson operator is thus complete. Calculating L (δ E/δx) with the explicit forms (A.2, A.3, A.7, A.8) for the Poisson operator, and (13) for the energy derivatives, one obtains the following reversible contributions to the evolution equations, whereσ αγ andσ S αγ are given by the following relationŝ It should be noted thatσ αγ andσ S αγ have a form closely related to the stress tensor and its entropic part, respectively. However, the complete form for the stress tensor expression might be calculated only after including the irreversible parts of evolution equations as well (particularly, see Sect. 3.4.2). In the relations above, and in the sequel, the abbreviations Note that has the dimension of temperature, and if E and S derive from a Helmholtz free energy, then = T .

Irreversible dynamics: friction matrix M
In this section, the irreversible contribution to the evolution equations, namely the second term on the right-hand side of (1), is constructed. This is the place where one implements the viscoelastic glassy-bridge dynamics, which will be discussed in Sect. 3.4.1. However, there is an additional irreversible contribution. It has to do with the fact that instead of a many-particle system we only consider a representative particle pair, which brings about nonaffine deformation effects, as will be analyzed in Sect. 3.4.2. It is noted that all the conditions for the friction matrix discussed in Sect. 2 are linear in M. Therefore, if one writes the friction matrix as a sum of two contributions, and if both contributions satisfy all GENERIC conditions independently, then also the sum M satisfies all conditions. In the following,M stands for the glassy-bridge dynamics, i.e., the relaxation of Q toward R, and M represents nonaffine deformation of R. It will turn out thatM is symmetric and positive semi-definite, whilẽ M is antisymmetric. Heat conduction is not described in this paper, since its GENERIC formulation has already been given in Öttinger [33] and several other papers.

Dissipative matrixM
By analyzing many examples of dissipative dynamics in the literature, a split of the dissipative matrixM in the following way has been proposed [38],M which is essential to establish a link [33] to linear irreversible thermodynamics [39]. In that factorization, the matrix D contains all dynamic material information (e.g relaxation times, viscosity, thermal conductivity), while the operator C represents the so-called mechanical part. It has been noted that is a useful condition to construct C * , and hence C. Not only is the degeneracy conditionM δ E/δx = 0 automatically satisfied, but also it has been found that the thermodynamic driving forces d td for the dissipative processes can be written in the form The corresponding thermodynamic fluxes j td are then given by In other words, C * translates the entropy gradient into the thermodynamic driving forces (24), while C determines how the thermodynamic fluxes j td enters the evolution equations, By requiring the symmetry and positive semi-definitness of the matrixM are guaranteed.
It has been explained in Sect. 3 that the glassy bridges yield under applied load, and thus, rate-dependent effects come into play. As discussed, this is to be described by a relaxation process of Q toward R. Since this Q-dynamics occurs in three dimensions, D will be a 3 × 3-matrix. In view of (26), C therefore maps a 3-vector j td onto the irreversible contributionsẋ| irr . While the detailed derivation of C is described in "Appendix B," the main steps are summarized here. The Q-relaxation leads to no contributions to the evolution equations of the momentum density u and the total deformation gradient F, which implies that some rows in C vanish. For the other elements, the p-and T -related rows, one proceeds as follows. Due to the conservation of probability in ξ -space and the number density of representative particle pairs, the evolution equation of the variable p is that of a density of a conserved quantity, and therefore for the relaxation dynamics of the glassy bridge captured in j γ . The form (28) can be used to determine the p-component of C, C ( p) , and due to the symmetry condition (4b) follows C * ( p) . That latter result combined with the degeneracy condition (3b) results in C * (T ) , and again with (4b) one obtains C (T ) . Collecting all contributions and using (24)(25)(26), one finds in addition to the p-dynamics (28) the evolution equation for the temperature, with the current related to the driving force according to (25) with a yet unspecified matrix D. For the driving force (24) for Q-relaxation, one obtains which in turn enters the constitutive relation for the flux, (25), with a yet unspecified 3 × 3-diffusion matrix that depends on both ξ and ξ .

Antisymmetric contributionM
It has been discussed in Sect. 3 that the force to restore the structure to the unloaded state must depend on the difference R − Q. If this force was derived from a potential, , it means that the potential will depend on |R − Q|. The ramifications for the dynamics discussed up to this point are the following. In the momentum balance (17a), the quantityσ plays the role of the stress tensor. In the last contribution toσ in (18a), δ A/δp is to be identified with the restoring potential (see also Sect. 4). If the latter depends on R and Q only through |R − Q|, this implies that ∇ R is parallel to R − Q, and hence,σ is in general not symmetric. This is problematic since the related torques result in an unwanted rotation of the system. While in a many-particle system (e.g., [18,19]) torques may arise locally in small clusters of particles, they may cancel out when summing up the contributions of all interacting particles. However, in the approach presented in this paper, we are actually considering a single representative pair; therefore, no such cancelation of nonsymmetric contributions is possible. Since for silica-filled elastomers there is no reason for expecting a nonsymmetric stress tensor, the cancelation of torques needs to be included explicitly in the model. The Poisson operator is representative of affine deformation; see Appendix B in [33]. However, as shown in Sect. 3.3, implementing affine deformation automatically results in a nonsymmetric stress tensor, due to the rigid connection between the dynamics (15, 17d) and the expression forσ . As a result, any modification of the model to arrive at a symmetric stress tensor must be implemented through the irreversible dynamics in (1). However, since this contribution to the dynamics is related to the kinematics and thus must not be dissipative, we look for an antisymmetric contributionM in (21), i.e.,Ṡ| ∼ = (δS/δx) M (δS/δx) = 0. As a matter of fact, a rather analogous case of nonaffine kinetics is known in the field of complex fluids, for which the Gordon-Schowalter derivative is used (see [40][41][42], or in [33, pp. 119-120] for details), representing the relative slippage of mesoscopic objects relative to the rotational components of the macroscopic background deformation.
The guiding principle for constructingM is the elimination of the antisymmetric contribution in the stress tensor. Following the above discussion and in analogy to [33, pp. 119-120], the antisymmetricM can be calculated, using the conditions of degeneracy (4b) and anti-symmetry forM, see "Appendix B" for details. The resulting modifications of the model arising from the nonaffine contribution can be written in the form withσ non-a αγ = Both of these stress-like expressions are antisymmetric. This underlines that the rotational part in the flow field, see (32c), goes hand in hand with the antisymmetric stress contributions (33), which approves the assumption made further above. The ramifications will be more explicit when summarizing all components of the model.
with the symmetric part of a matrix A sym = (A + A )/2. The constitutive relations for the total stress tensor σ and its entropic part σ S are given by It can be shown [4] that both the total and the entropic stress tensor expressions are symmetric, since the internal energy U and entropy S may depend on the deformation gradient F only via right Cauchy-Green strain tensor, C = F · F. Furthermore, the additive form of the stress tensor expression (i.e., the elastomer matrix contribution plus the particle contribution) is similar to the Eindhoven glassy polymer (EGP) model [63,64]. While the matrix contribution, similar to the EGP spring element, allows to mimic even neo-Hookean material behavior, the glassy-network contribution corresponds to the phenomenological nonlinear Maxwell element in the EGP model. The irreversible p-current density is given by (31) and (30), where the 3 × 3-diffusion matrix depends on both ξ and ξ .
The above model makes the concurrent two-way coupling of the macroscopic scale and the particle-level evident. The macroscopic flow field v α affects the particle level p. This drives the particle arrangement away from the equilibrium state, which in turn results in nonequilibrium forces on the particle level. These forces are felt on the macroscopic scale by way of the constitutive relation for the stress tensor. In effect, one thus obtains the constitutive behavior of the nanocomposites, relating the stress to the deformation. For specific applications, the model is completed by specification of the functionals of internal energy U and entropy S, and diffusion tensor D, which is illustrated in the following section.
We close the general model development by reminding the reader that the model presented above is substantially different from common practice in polymer kinetic theory [33,35,36]. As discussed in Sect. 3.1, modeling filled elastomers in the representative pair setting requires the concept of an unloaded state, Q, which is unknown to polymer kinetic theory. As a consequence, the irreversible but nondissipative contributions to the dynamics in relation to the symmetrization of the stress tensor are absent in common polymer kinetic theory.

Energy and entropy
For illustration purposes, rather than out of necessity, we assume that both U and S can be split additively into contributions related to (i) the elastomer matrix and (ii) the particle configuration. For the potential that accounts for the glassy-bridge interactions discussed earlier in Sect. 3.4.2, we use a Hookean spring potential, with spring constant k. Since the modulus of a glassy polymer depends on temperature (see Ch. 13 in [43]), the same holds for the spring constant, and in turn for the potential in (37). Since plays the role of a Helmholtz free energy (per particle pair), its temperature dependence implies that it decomposes into energetic and entropic parts, namely respectively. This split of is reflected in the expressions for the functionals of internal energy and entropy (see also p. 135 in [33]), where the subscript "em" denotes the contributions of the pure elastomer matrix. The symbol k B stands for the Boltzmann constant, and the related contribution to the entropy is of statistical origin [44,45]. For details on this entropy contribution, the reader is referred to "Appendix C." As far as the contributions due to the elastomer matrix are concerned, u em and s em , we note that the modulus of an elastomer (in its rubbery state, i.e., above T g ) is influenced by temperature, most often approximated by a linear T -dependence (see Ch. 13 in [43]). Therefore, a split of the Helmholtz free energy density a em analogous to (39) can be used, for the energetic and entropic contributions, respectively, which yields For the example in this section, the Helmholtz free energy for compressible neo-Hookean systems is used [46], with shear modulus G em , bulk modulus κ em , I 1 = tr(F T · F) and J = det F. It is noted that by virtue of the expressions for U and S in (40) and the splits (39) and (41), the quantity (20) simplifies to = T . Furthermore, δ A/δp used in the flux (36) and the stress tensor (35a) becomes δ A/δp = + k B T (ln p + 1), while the stress tensor itself assumes the form with n defined in (9), the left Cauchy-Green strain tensor B αγ = F αν F γ ν , and given by (37). Apart from the isotropic contributions, the term proportional to the shear modulus G is typical for neo-Hookean materials. Most interesting is the last contribution in (44). It corresponds to known stress tensor expressions in statistical mechanics [47][48][49] and in polymer kinetic theory [33,35,36], a.k.a. Kramers-Kirkwood expression. However, the major difference consists in the potential . One should recall that for the filled elastomers discussed here, the potential depends not only on the current state R, but also on a unloaded state Q, e.g., according to (37).
As it was already mentioned in Sect. 3.1, the spring constant k needs to account for the collective effect of all the neighboring particles in order to mimic properly, in the representative particle-pair setting, the restoring mechanism as the glassy-spring network is deformed. If an effective particle pair is defined by a surfaceto-surface separation smaller than half the particle diameter d, the spring constant can be approximated by (see "Appendix D"), with G comp the shear modulus of the nanocomposite. G comp depends on the volume fraction φ, and might be determined experimentally; e.g., from [18], one gets G comp = {4, 10, 50}MPa for φ = {0.2, 0.3, 0.4}. The values for the number density of particle pairs n (see "Appendix D") and for the resulting spring constant k as a function of volume fraction φ are given in Table 1.

Dynamics
Since the material has no inherent anisotropy, it is reasonable to assume that the relaxation of the glassy bridges, j, occurs along the direction of the thermodynamic driving force d td . According to (31), this implies that D is a multiple of the unity tensor. Furthermore, we assume for illustrative purposes that the probability current j(ξ ) depends on the driving force d td (ξ ) only for ξ = ξ . One can thus write with a friction coefficient ζ > 0. Without loss of generality, the positive factors and p have been introduced for later convenience. The assumptions made for the energy and entropy (40), the glassy-spring potential (37), and the diffusion matrix (46) render the evolution equation (34d, 36) into the form This Fokker-Planck equation can be translated into an equivalent set of stochastic differential equations [36], which will be used for the numerical simulations in Sect. 5. The last term in (48b) represents the fluctuating Brownian contributions, which adds stochasticity to the Q-dynamics. The symbol dW t stands for increments in Wiener processes, with average and variance, d W tα = 0 and d W tα dW t β = δ αβ δ tt dt, respectively. In other words dW t represents white noise, being uncorrelated in time.

Relaxation of glassy bridges
In the dynamic model (48), the occurrence of a relaxation time τ for the Q-dynamics is evident, with This relaxation time contains all the information about the temperature and load dependence of the glassybridge dynamics. In the following, we will discuss a possible choice for the relaxation time τ , which for given spring constant k will imply a specific choice for the friction coefficient ζ .
To determine the dependence of the relaxation time τ on temperature and loading, we refer to [18][19][20]. In that study, the concept of a local glass transition temperature is invoked (see also Sect. 1). If the matrix material adheres well to the filler-particle surface, the elastomer chain segments have decreased mobility the closer they are located to the particle, which is captured by a (local) glass transition temperature that increases toward the filler-particle surface. It was proposed in [6,7] that the glass transition temperature T g (z) at distance z from a particle is given by T g (z) = T g,b (1 + β/z), which could fit their results well over a 100 K temperature range. Based on the relation for T g (z), the thickness of the glassy layer z th at temperature T can be estimated by T g (z th ) = T , i.e., z th = βT g,b / T − T g,b . The divergence in T g (z) is purely formal, since the physics which lead to this relation [50,51] has a natural cutoff at short distance of about 2 nm. The corresponding increase in the glass transition temperature, about 100 K at 2 nm in their system, corresponds to the increase of T g measured in [7], either in mechanical experiments or by using NMR. Further support of the glassy-layer picture has been given in [52][53][54], where glassy volume fractions have been measured for model samples using DSC, which were consistent with NMR and mechanical measurements.
If two particles are separated by a surface-to-surface separation of z, the most mobile chain segments are located in the middle, which will dominate the relative motion of the two particles. This relative motion of the two particles is influenced by the local glass transition temperature T g in the middle, that in general not only depends on the distance z but also on the local load loc [18,19], where T g,b is the bulk glass transition temperature of the pure elastomer, β and K are material parameters. K relates the yield stress σ y to the temperature below T g by σ y = K (T g − T ), and is a physically measurable quantity. The dependence of T g on loc mimics the plasticization of the glassy bridge upon loading, which resembles the well-known plasticization/yielding of bulk polymer glasses upon loading [55,56]. The change in T g due to the presence of fillers and local load is relevant for the α-relaxation time described by the Williams-Landel-Ferry (WLF) law, where C 1 and C 2 are the WLF parameters specific to the polymer, and τ g represents the relaxation time at T = T g .
In the context of the model presented in this paper, one can incorporate the approach of [18,19] by using the identifications τ = τ α (i.e., ζ = k τ α ), and with d the particle diameter. As for the quantity loc , it has been identified in [18,19] with the absolute value of the glassy-bridge force. In contrast, modeling of bulk polymer glasses suggests to identify loc with the stress [56]. However, in our setting of a representative pair, the stress as such in between the two particles is not well defined. Since resemblance to the concept of stress is nevertheless desirable, the following choice has been made in analogy to established von Mises J 2 -plasticity, e.g., [57,58], with the second invariant inspired by the stress tensor expression (44). The loc in (55) should account for all stress contributions, and therefore, t should include not only the glassy-bridge contribution in (37) but also an elastomer matrix contribution, denoted by em . In principle, the elastomer behavior is captured in our approach in a continuum sense rather than on the particle level [see, e.g., (40) and (44)]; however, for the sake of a meaningful definition of loc one may recast it on the particle level into the form em (ξ ) = 1 2 k em (R − Q 0 ) 2 , analogous to the procedure outlined in "Appendix D", with k em the spring constant representing the elastomer matrix. Since the elastomer matrix does not yield but is purely elastic, the corresponding strain energy depends on the difference between the current state R and Q 0 , the initial unloaded state.
It is noted that the description given above with a load-dependent relaxation time closely resembles the modeling of viscoplastic behavior in bulk polymers [56,59]. Particularly, yielding comes about by the strongly nonlinear dependence of the relaxation time on the applied load, and hence, the use of a classical yield criterion [57,58] is not required. However, in comparison with modeling bulk polymers in the post-yield regime, see, e.g., [63,64], the model in this paper is different. First, strain softening, being connected with physical aging, is not in the scope of the presented work. And second, the particles in the two-particle model should always remain neighbors, which limits the applicability of this model to the deformation regimes where strain hardening becomes relevant (about 30 % and higher [63]).

Simulation setup
Using all the system specifications of Sect. 4, the dynamics (48) is solved numerically (see [36]) to determine the stress response (44) upon imposed oscillatory deformation, under isothermal conditions, i.e., constant T .
The initial particle configuration for each trajectory corresponds to an unloaded configuration, R = Q, with R being isotropically distributed. The length R = |R| was sampled from a narrow distribution of the form w(R) with amplitude γ 0 and angular frequency ω.
The parameters used are listed in Table 1 and are the same as the ones considered in [18][19][20], which is representative of a filled elastomer consisting of silica particles dispersed in a poly(ethyl acrylate) matrix. All these parameters can be related to physical parameters, as explained in [18]. In particular, the value we chose for the parameter β = 0.8 nm, which sets the amplitude of the T g -shift in the vicinity of the fillers, has been considered in [18], and corresponds to that measured experimentally in [6,7]. The only parameter that is slightly adapted from the many-particle simulations in [18][19][20] to our two-particle model is the stress sensitivity K in (51). In [18], K is determined such that, when imposing an oscillatory strain with amplitude 0.1 and rate 0.1 s −1 , the stress-induced shift of the glass transition temperature is 100 K at the maximum strain, which would result in K = 1.2 · 10 −20 J K −1 . However, in order to get better agreement between our model and the results in [20] (see further below), it was advantageous to choose the value K = 2 · 10 −20 J K −1 . In the simulations, the stress-induced shift in the glass transition temperature (51) has been limited to ≤200 K, since above that value the system is already plasticized (see also [19]). Finally, for numerical reasons, the range of the relaxation time (52) has been limited to 0.01 s ≤ τ α ≤ 100 s [20].
To express the frequency relative to an internal rate of relaxation, it is below presented in terms of ω · τ init α , where τ init

Effect of deformation amplitude
The stress-strain curves for small-and large-amplitude oscillatory deformation are shown in Figs. 3 and 4, respectively. For small amplitude (Fig. 3), for the slow deformation rate (solid line), the stress-strain curve has a substantial viscous component, leading to hysteresis and thus dissipation, while for increasing deformation rate the hysteresis diminishes until almost pure (linear) elastic behavior is recovered (dash-dotted line). Such kind of behavior is typical for a linear viscoelasticity.
The above results suggest that the material responds elastically if the applied rate of deformation is fast with respect to the rate of structural re-arrangements. In contrast, the material responds viscous-like if the applied rate of deformation is slow with respect to the rate of structural re-arrangements. This argument is straightforward in the case that the internal timescale does not depend on the imposed load and leads to linear viscoelastic behavior. However, the discussion in the previous section has shown that the internal timescale in the case of filled elastomers, i.e., of the glassy bridges, depends strongly on the applied load. This will lead thus to a nonlinear material response. While elliptic curves in a Lissajous-Bowditch plot are a hallmark of linear viscoelasticity, it is observed upon increasing the deformation amplitude (e.g., γ 0 = 0.2) that the stress-strain curves are significantly nonelliptic in shape, see Fig. 4. However, even in the nonlinear case, increasing the frequency leads to a narrowing of the loop in the Lissajous-Bowditch plot, meaning that the behavior becomes more elastic.

Measure of the nonlinearity: Payne effect
Nonlinear effects have been observed in Fig. 4, where large-amplitude oscillatory shear (LAOS) deformation is imposed, resulting in a nonelliptical stress-strain curve. To quantify the nonlinearity, i.e., nonellipticity, two approaches are considered. For a given deformation γ (t) (56), the stress response obtained from the simulation is denoted by σ data (t). One can then attempt to represent this dataset by a the common linear viscoelastic relation σ fit (t) = γ 0 G (ω) sin(ωt) + G (ω) cos(ωt) . An optimal representation (see Fig. 5 for an example) can be defined as the specific choice for the parameters G (ω) and G (ω) that minimizes the (dimensionless) functional The quality of the resulting best fit can then be quantified by the corresponding value of N . Note that N = 0 for linear viscoelasticity.
A different way to quantify the nonlinearity in the system was developed in [60] and extended in [61,62]. In the latter, two different elastic moduli and two different viscosities have been defined based on the Lissajous-Bowditch plot [61,62], which in turn are used to define two dimensionless indexes of nonlinearity [61,62], For linear viscoelasticity, both G M = G L and η M = η L , and therefore, S S = S T = 0. The parameters S S and S T can hence be used as a dimensionless quantification of the significance of nonlinear effects. While the quantity N tells directly how the stress-strain curve differs from the elliptical shape (i.e., from linear viscoelasticity), S S > 0 corresponds to strain stiffening, and S S < 0 stands for strain softening [61]. Furthermore, S T > 0 represents intracycle shear thickening, and S T < 0 intracycle shear thinning [61]. In Figs. 6 and 7, the nonlinearity of the mechanical response of the model is presented in terms of the parameters N defined in (57), as well as S S and S T defined in (59). At fixed strain amplitude γ 0 , the nonlinearity parameter N decreases by two orders of magnitude upon increasing the deformation rate by two orders of magnitude. The quantities S S and S T also show the same overall trend, with decreasing nonellipticity with increasing deformation rate, however with significantly more scatter than for N . While the modulus-based quantity S S shows a strong dependence on the deformation rate, the viscosity-based S T does not permit such a conclusion.
A prominent effect of the nonlinear viscoelastic response of filled elastomers is the Payne effect, i.e., the dependence of G and G on the amplitude of the applied deformation. For a quantitative evaluation of this effect in the present model, the ratios G M /G 0 and G L /G 0 are employed, where G M , G L are the tangential and secant moduli defined in (58a), respectively, and G 0 is a shear modulus at small amplitude (linear regime), in this particular case for γ 0 = 0.005, leading to Fig. 8. For both definitions, there is a decrease with increasing deformation amplitude, where the decrease is stronger for the tangent-based G M /G 0 than for the secant-based G L /G 0 . However, the trend of both curves agrees qualitatively with the original work of Payne [11] and also corresponds closely to the results of the many-particle model [19,20]. The mechanism leading to the Payne effect in this model can be summarized as follows: An applied macroscopic deformation induces local interparticle stresses, which in turn reduce the effective inter-particle glass transition temperature (see expression (51)), thereby lowering the local relaxation time. This implies that the ability of the microstructure to relax (part of) the imposed load is increased, which is then observed as the Payne effect. The fact that the Payne effect is qualitatively reproduced underlines that the most essential aspects of the filled elastomer are appropriately represented in the model. 5.4 Comparison to many-particle simulations [20] As stated in Sect. 1, the goal of this paper is to design a two-scale model to describe the main effects in filled elastomers; however, it should be more compact than the many-particle approach in [18][19][20]. The above detailed discussion of system specification and parameter estimation makes clear that the model presented in this paper is closely related to the one in [18][19][20]. Correspondingly, the simulation results discussed above are analogous to the ones in [18][19][20]. Nevertheless, it is also clear that the reduction from a many-particle system to a representative particle pair, this paper, comes at a price. The model comparison is performed on the basis of the mechanical response under oscillatory deformation γ (t) in (56). The simulations are performed using the parameters in Table 1. The only exception is that for this comparison β = 0.4 nm is used, in agreement with the specific value used in Fig. 15 in [20]. The results of the comparison can be seen in Figs. 9 and 10. At moderate strains, the agreement with the many-particle approach is rather satisfactory, keeping in mind the drastic simplification of our model. The most apparent difference between the many-particle and two-particle models is the behavior in the outer strain regions. In the many-particle model, the stress-strain curves show strain hardening at large amplitudes of deformation. It has been argued that this is due to the particle heterogeneity (see [20] for details). Naturally, this is an effect that is not accounted for in our effective particle-pair approach, and therefore, this is a shortcoming that is inherent to the substantial reduction in model complexity.

Conclusion
In this paper, a two-scale model was developed for describing the mechanical behavior of hard particle-filled elastomers, in agreement within nonequilibrium thermodynamics principles, specifically the General Equation for the Non-Equilibrium Reversible-Irreversible Coupling (GENERIC). The two scales considered are on the one hand the macroscopic scale on which mechanical tests are performed, and on the other hand the scale of the filler particles, specifically the glassy-bridge network. So doing, the intricate mechanical behavior of these composites can be captured by a mutual coupling of both scales. This results in the set of concurrent evolution equations (for the momentum density, the temperature, the total macroscopic deformation gradient and the structural variable), supplemented by the constitutive relation for the stress tensor, which is summarized in Sect. 3.5. The presented model can be applied to hard particle-filled elastomer systems with not too low particle volume fraction, i.e., in the presence of a glassy-bridge network.
The proposed two-scale model has been solved numerically for macroscopically homogeneous oscillatory deformation, for both small and large amplitude. Of specific interest is the large-amplitude oscillatory test (LAOS). Under LAOS for strain amplitudes approximately equal or higher than 0.1, the glassy bridges between the particles yield. This results in a dramatic decrease in the effective elastic modulus, G , in the literature known as the Payne effect [11], and is related to the stress dependence of the glassy-bridge relaxation time. The nonlinearity of the mechanical response has been quantified and discussed in Sect. 5.3. The numerical results compare favorably with the literature.
Similar effects as modeled in this two-particle model have also been accounted for in the many-particle approach in [18][19][20]. While the general features of the stress-strain response of both models agree, there are also marked differences, as shown in Sect. 5.4. Namely, at high strains, it has been argued in [20] that the heterogeneity in particle arrangement leads to strain hardening. Such effects are not present in the twoparticle approach taken here. More detailed comparisons between the two approaches will be done in the future, also in relation to the Mullin's effect (see below). However, it is already clear now that both approaches have their benefits, depending on application. On the one hand, if one is interested in the material behavior in a macroscopically homogeneous sample, one may choose for using the many-particle approach, since it captures the many-particle aspects more appropriately. On the other hand, if one is interested in studying a macroscopically inhomogeneous situation, one would need to solve in each volume element the many-particle simulation to express the stress state in terms of the deformation (history). One can anticipate that this is computationally too costly, and alternatives must be sought. A possible solution is to use the two-particle approach developed here, which is computationally substantially more efficient than the many-particle model, and simplifies the calculation in each volume element of the stress in terms of the deformation. This might open the way toward two-scale modeling of inhomogeneous phenomena, such as wear and tear, of hard particle-filled elastomers.
In future work, the model presented here will be extended to describe also the Mullins effect [12,13], i.e., the recovery of the effective elastic modulus with increasing waiting time after application of strong deformation. It seems that the origin of this recovery behavior roots in the physical aging of the glassy bridges, akin to the physical aging of bulk glassy polymers. Therefore, the proper modeling of the Mullins effect will substantially benefit from modeling efforts to describe physical aging in bulk glassy polymers, e.g., as described in the Eindhoven glassy polymer (EGP) model [63,64]. On the one hand, the Poisson operator and friction matrix depend in general on both r and r (r is different from r), i.e., L(r, r ), M(r, r ). The evolution equation (1) thus implies thatẋ(r) is affected not only by the driving forces at the same space point r, but also by the driving forces at any other points r . However, for local field theories, there are only contributions for r = r , and therefore, no integration is needed, i.e., both the Poisson operator and the friction matrix vanish for any r = r. The situation is different for the variable ξ in (11) describing the internal configuration, as can be seen by the following argument. The momentum balance equation has the form ∂ t u(r) = −∇ r · [u(r) v(r)] + ∇ r · σ (r). Here the stress tensor σ depends implicitly on the variable ξ through the distribution function p(r, ξ ), which implies that an integration over the variable ξ occurs, in the sense of an averaging operation. One thus concludes that implies an integration over internal configurations, and therefore, the Poisson operator and the friction matrix depend in general on both ξ and ξ , L = L(r, ξ , ξ ) and M = M(r, ξ , ξ ). Therefore, the time evolution equation (1) can be written in the more explicite form where the indices i and j refer to state variables in the set x. We are thus dealing with a field theory that is local with respect to the macroscopic position r, but nonlocal with respect to the internal configuration ξ .
To simplify notation in this paper, we are going to use the following notation. Since the general evolution equation (A.1) is local in terms of r, the argument r is omitted, except in special cases. As for the dependence of a function f on the internal configuration ξ , we use the shorthand notation and f = f (ξ ), respectively, unless the explicit notation with arguments is required for clarity.
Following the arguments given in Sect. 3.3, and using the result obtained in [4,33], one can show that the Poisson operator takes the form Here, it should be pointed out that all ∇-operators occurring in (A.3) act also on functions multiplied to L from the right (e.g., δ E/δx), except when enclosed in parentheses like (∇ r f (r)), where ∇ r acts only on the function f (r) enclosed in the parentheses.
The two yet unspecified elements L (uT ) and L (T u) can be determined by using the degeneracy condition (3a) and the anti-symmetry condition (4a). From the degeneracy condition (3a), it follows where we have used δS/δu = 0. Note, that without violating the degeneracy condition, we can add a term of the form f (x)∇ r α 1 to (A.4), with any arbitrary function f (since ∇ r α 1 = 0). This implies that L The L (T u) -entry can be determined by way of the anti-symmetry of Poisson operator, leading to So far we have respected all GENERIC conditions when specifying the Poisson operator, except the Jacobi identity (5). To check the Jacobi identity, it is useful to note that a transformation of dynamical variables does not affect the Jacobi identity. Particularly, it can be shown that the Poisson operator is considerably simplified when using x = {u, s, F, p} instead of x = {u, T, F, p}. Using the x -formulation, the Jacobi identity can be checked and verified analytically, by a lengthy but straightforward calculation. For a computer code to symbolically calculate the Jacobi identity, the reader is referred to [66].

Dissipative matrixM
In this appendix, the elements of C in the decomposition (22) are determined. To that end, the reader is reminded (as in (A.1)) that the symbol in (22) also implies integrations, and therefore, the explicit form of (22) is given byM Since the relaxation dynamics, j, of the unloaded state Q toward the current state R originates from the viscoplastic deformation of glassy material and hence involves friction-related, dissipative effects, the temperature of the system will be affected by j. In view of the complete set of variables x = {u, T, F, p}, the general operator C in (26) thus assumes the form with adjoint operator C * . The symbol • denotes yet unknown elements, which are determined as follows. The p-evolution equation (28) can be re-written in the form 4 Use has been made of the relation (in analogy to Eq. (14) in [65]) which, in view of (26) and the ansatz (B.2), leads to with C * ( p) α the operator adjoint of C . It should be mentioned that, similarly to the discussion in relation to (A.4), a term of the form f (x)∇ r α (δU/δT ) −1 could be added to C * (T ) α without violating the degeneracy condition, with an arbitrary function f (x). However, it can be shown that such a contribution leads to heat conduction [28], in which we are not interested in this paper, as stated earlier. Finally, the explicit form of C and C * is given by In turn, the thermodynamic driving force (24) assumes the form where the abbreviation (19b) has been used. Once a relation for a symmetric and positive semi-definite D is specified, the expression for the thermodynamic flux (25) For the construction of all components, we start by assuming 5 that the anti-symmetry of the stress tensor is related to the rotational component in the applied flow field, i.e., to the antisymmetric contribution to the velocity gradient, ∇ r v. In order to eliminate the rotational effects of the flow field, we strive to replace the velocity gradient in (17d) by its symmetric part. This is achieved by adding to the dynamics of the structural variable p (17d) an additional contribution in the irreversible part, specifically accounting for the nonaffinity with The identification of the elements inM is done as follows. Since δS/δu = 0, one can use (B.9) to determinẽ M ( pT ) by virtue of ∂ t p| non-a =M ( pT ) (δS/δT ), from whichM ( pu) follows by using the degeneracy condition (3b). Consecutively,M (T p) andM (up) follow from the condition of anti-symmetry ofM. Again using the degeneracy condition (3b) leads toM (T u) andM (uT ) . It must be mentioned that at two points of this procedure a nonuniqueness of the elements ofM arises, leading to four possible solutions. However, when requiring that M (T u) andM (uT ) are such that the entireM is antisymmetric (which has not yet been enforced), only one of the four solutions remains, and is hence unique, namelỹ

Appendix C: Configurational entropy
In this Appendix, the expression for the entropy functional S in terms of the structural variable p is discussed. While the use of the p ln p-expression is commonly accepted for the case that p is the 1-or N -particle distribution function of an N -particle system, the situation is less obvious for the model discussed in this paper. In our case, p depends on two vectors, namely R and Q, where both of them are difference vectors between two particles. However, due to translation invariance, 6 in terms of degrees of freedom the pair (R, Q) can be seen as representing the absolute positions of two particles. The question is thus how the entropy can be represented in terms of the 2-particle distribution function, for which we closely follow the developments in [67] in the following. Consider a system wit N interacting particles with phase space coordinates = (x 1 , x 2 , . . . , x N ) ≡ (x N ). The probability density for N particles to be at point (at time t) is denoted by f (N ) ( ) and satisfies the normalization condition [67], Based on (C.1), one can introduce the reduced 1-and 2-particle distributions, respectively, as In [67], the entropy has been calculated neglecting higher than 2-particle correlations. The result can be written in the form (see [67] for details), with k B the Boltzmann constant, and "eq" denotes the equilibrium quantities. It was noticed by Wallace earlier that entropy expressions of this type are well convergent for dense fluids (see p. 2283 in [68]). Moreover, using computer simulations, it was shown in [69] that the 3-particle (and higher-order) terms are small in comparison with the 1-and 2-particle terms. Relations closely related to (C.3) are often used in the literature, e.g., [69][70][71][72].
Taking into account that based on (C.2), it is manifest that the 2-particle distribution contains also the 1-particle distribution. Therefore, choosing the 2-particle distribution function as the modeling quantity automatically allows to also account for the 1-particle contributions in (C.3). However, when seeking to describe the dynamics of such a f (2) -model that involves the maximization of entropy, the driving force f (1) eq (x 1 ) eq (x 1 , x 2 ) contains mean-field-type effects in the form of f (1) being an integral of f (2) . One can anticipate that the final evolution equation for f (2) will contain these mean-field contributions as well, and it can be shown that one deals effectively with a Fokker-Planck equation for f (2) with mean-field effects. Since the purpose of Sect. 4 is to only give an illustrative example of the general model, the above issues of the entropy expression are not examined any further in this paper, but rather left as a topic for further studies. It is pointed out that Fokker-Planck equations with mean-field effects need special care in their interpretation, and the reader is also referred to [73] for more details.

Appendix D: Spring constant and representative particle-pair number density
First, an estimate is given for the value of the number density of particle pairs, n, which has been introduced in Sect. 3.1, Eq. (9). We begin with the relationship between pair number density, n, and the particle number density m (both are considered per unit volume) with p the average number of neighboring particles per particle. In the following, both m and p, and thus n, will be related to the volume fraction φ. For the pair number density n, we proceed as follows. A particle 1 is considered to be a neighbor of another particle 2 if any fraction of particle 1 is located closer than a certain distance l/2 from the center-of-mass of particle 2. Considering a sphere of radius l/2 of volume V l , the (given) volume fraction can be expressed in the form φ = N l V p /V l , with V p the volume of a single particle and N l the number (including fractions) of particles contained in the sphere volume V l . Using V p /V l = (d/l) 3 , one obtains N l = φ (l/d) 3 , and thus, the number of neighboring particles is given by p = N l − 1, Combining the particle number density (D.2) with the average number of neighboring particles (D.3), the average number density of particle pairs becomes which is complete once a choice for the length l has been made. Knowing the number density of particle pairs, n, we are able to calculate the value of spring constant k, as follows. As it was discussed in Sect. 3.1, the spring constant depends on the connectivity and thus on n. Let us consider the stress tensor at equilibrium. For simplicity, only the glassy-bridge contribution to the stress tensor is considered, and since the force on each particle in mechanical equilibrium vanishes, one obtains Perturbing the equilibrium position R by a small amount δR, the first-order change in the stress tensor takes the form where we have used the split (8), and R α R ν eq = (1/3) |R| 2 eq δ αν . Relation (D.7) can be set equal to δ σ αβ = G comp αβ , with G comp strictly speaking the shear modulus of the glassy-bridge network. However, since the modulus of the "underlying" matrix material in the rubbery state is much smaller than that of the glassy material [43], one may approximate the modulus of the glassy-bridge network by that of the entire composite. One thus obtains the following relation for the spring constant, where in the second equality the relation (D.4) and the estimate |R| 2 eq d 2 /φ 2/3 for the square of the average neighbor distance has been used. For the case that two particles are called neighbors if their surface-to-surface distance is smaller than d/2 (meaning l/2 = d), one obtains